From d2a4268ff2eeea179439ff760fe232c11da2e249 Mon Sep 17 00:00:00 2001 From: Dongryeol Lee Date: Wed, 26 Sep 2007 01:32:46 +0000 Subject: [PATCH] Test case added for Three way convolution of Gaussian sums --- .../series_expansion/farfield_expansion.h | 22 ++++--- fastlib/u/dongryel/series_expansion/main.cc | 64 ++++++++++++++++--- 2 files changed, 69 insertions(+), 17 deletions(-) diff --git a/fastlib/u/dongryel/series_expansion/farfield_expansion.h b/fastlib/u/dongryel/series_expansion/farfield_expansion.h index 54e9e6fa80..3f41f617de 100644 --- a/fastlib/u/dongryel/series_expansion/farfield_expansion.h +++ b/fastlib/u/dongryel/series_expansion/farfield_expansion.h @@ -397,15 +397,20 @@ double FarFieldExpansion::ConvolveField const FarFieldExpansion &fe3, int order1, int order2, int order3) const { + // bandwidth factor and multiindex mapping stuffs + double bandwidth_factor = kd_.BandwidthFactor(bandwidth_sq()); const ArrayList *multiindex_mapping = sea_->get_multiindex_mapping(); const ArrayList *lower_mapping_index = sea_->get_lower_mapping_index(); - // get the total number of coefficients + // get the total number of coefficients and coefficients int total_num_coeffs1 = sea_->get_total_num_coeffs(order1); int total_num_coeffs2 = sea_->get_total_num_coeffs(order2); int total_num_coeffs3 = sea_->get_total_num_coeffs(order3); int dim = sea_->get_dimension(); - + Vector coeffs2, coeffs3; + coeffs2.Alias(fe2.get_coeffs()); + coeffs3.Alias(fe3.get_coeffs()); + // actual accumulated sum double neg_sum = 0; double pos_sum = 0; @@ -446,9 +451,9 @@ double FarFieldExpansion::ConvolveField xK_center.Alias(fe3.get_center()); for(index_t d = 0; d < dim; d++) { - xI_xJ[d] = center_[d] - xJ_center[d]; - xI_xK[d] = center_[d] - xK_center[d]; - xJ_xK[d] = xJ_center[d] - xK_center[d]; + xI_xJ[d] = (center_[d] - xJ_center[d]) / bandwidth_factor; + xI_xK[d] = (center_[d] - xK_center[d]) / bandwidth_factor; + xJ_xK[d] = (xJ_center[d] - xK_center[d]) / bandwidth_factor; } kd_.ComputeDirectionalDerivatives(xI_xJ, derivative_map_alpha); kd_.ComputeDirectionalDerivatives(xI_xK, derivative_map_beta); @@ -520,9 +525,8 @@ double FarFieldExpansion::ConvolveField sign += 2 * (alpha_mapping[d] + beta_mapping[d] + gamma_mapping[d]) - mu_mapping[d] - nu_mapping[d] - eta_mapping[d]; - sign = sign % 2; } - if(sign == 1) { + if(sign % 2 == 1) { sign = -1; } else { @@ -533,9 +537,9 @@ double FarFieldExpansion::ConvolveField moment_i = coeffs_[sea_->ComputeMultiindexPosition(mu_nu_mapping)]; moment_j = - coeffs_[sea_->ComputeMultiindexPosition(alpha_mu_eta_mapping)]; + coeffs2[sea_->ComputeMultiindexPosition(alpha_mu_eta_mapping)]; moment_k = - coeffs_[sea_->ComputeMultiindexPosition + coeffs3[sea_->ComputeMultiindexPosition (beta_gamma_nu_eta_mapping)]; pos_sum += sign * diff --git a/fastlib/u/dongryel/series_expansion/main.cc b/fastlib/u/dongryel/series_expansion/main.cc index 6747405bae..db82543c6b 100644 --- a/fastlib/u/dongryel/series_expansion/main.cc +++ b/fastlib/u/dongryel/series_expansion/main.cc @@ -351,35 +351,83 @@ int TestConvolveFarField(const Matrix &data, const Vector &weights, printf("\n----- TestConvolveFarField -----\n"); - // bandwidth of sqrt(0.5) Gaussian kernel - double bandwidth = sqrt(0.5); + // bandwidth of 5 Gaussian kernel + double bandwidth = 5; GaussianKernel kernel; - kernel.Init(sqrt(0.5)); + kernel.Init(bandwidth); // declare auxiliary object and initialize SeriesExpansionAux sea; sea.Init(20, data.n_rows()); - // declare center at the origin + // declare center at the origin, (10, -10) and (-10, -10) Vector center; center.Init(2); center.SetZero(); + Vector center2; + center2.Init(2); + center2[0] = 10; center2[1] = -10; + Vector center3; + center3.Init(2); + center3[0] = center3[1] = -10; + + // create fake data + Matrix data2, data3; + data2.Copy(data); + data3.Copy(data); + for(index_t c = 0; c < data.n_cols(); c++) { + data2.set(0, c, data2.get(0, c) + center2[0]); + data2.set(1, c, data2.get(1, c) + center2[1]); + data3.set(0, c, data3.get(0, c) + center3[0]); + data3.set(1, c, data3.get(1, c) + center3[1]); + } + + data.PrintDebug(); + data2.PrintDebug(); + data3.PrintDebug(); // declare expansion objects at (0,0) and other centers FarFieldExpansion se; + FarFieldExpansion se2; + FarFieldExpansion se3; // initialize expansion objects with respective centers and the bandwidth // squared of 0.5 se.Init(bandwidth, center, &sea); + se2.Init(bandwidth, center2, &sea); + se3.Init(bandwidth, center3, &sea); - // compute up to 4-th order multivariate polynomial. + // compute up to 20-th order multivariate polynomial. se.AccumulateCoeffs(data, weights, rows, 20); + se2.AccumulateCoeffs(data2, weights, rows, 20); + se3.AccumulateCoeffs(data3, weights, rows, 20); - // print out the objects - se.PrintDebug(); // expansion at (0, 0) + printf("Convolution: %g\n", se.ConvolveField(se2, se3, 6, 6, 6)); - printf("Convolution: %g\n", se.ConvolveField(se, se, 5, 5, 5)); + // compare with naive + double naive_result = 0; + for(index_t i = 0; i < data.n_cols(); i++) { + const double *i_col = data.GetColumnPtr(i); + for(index_t j = 0; j < data2.n_cols(); j++) { + const double *j_col = data2.GetColumnPtr(j); + for(index_t k = 0; k < data3.n_cols(); k++) { + const double *k_col = data3.GetColumnPtr(k); + // compute pairwise distances + double dsqd_ij = 0; + double dsqd_ik = 0; + double dsqd_jk = 0; + for(index_t d = 0; d < sea.get_dimension(); d++) { + dsqd_ij += (i_col[d] - j_col[d]) * (i_col[d] - j_col[d]); + dsqd_ik += (i_col[d] - k_col[d]) * (i_col[d] - k_col[d]); + dsqd_jk += (j_col[d] - k_col[d]) * (j_col[d] - k_col[d]); + } + naive_result += kernel.EvalUnnormOnSq(dsqd_ij + dsqd_ik + + dsqd_jk); + } + } + } + printf("Naive algorithm: %g\n", naive_result); return 1; }