diff --git a/fastlib/u/dongryel/kde/kde.h b/fastlib/u/dongryel/kde/kde.h index 9dff59123b..6abe3751e8 100644 --- a/fastlib/u/dongryel/kde/kde.h +++ b/fastlib/u/dongryel/kde/kde.h @@ -4,7 +4,7 @@ #include "fastlib/fastlib_int.h" #include "u/dongryel/series_expansion/farfield_expansion.h" #include "u/dongryel/series_expansion/local_expansion.h" -#include "u/dongryel/series_expansion/kernel_derivative.h" +#include "u/dongryel/series_expansion/kernel_aux.h" template class NaiveKde { @@ -105,7 +105,7 @@ class NaiveKde { }; -template +template class FastKde { public: @@ -156,12 +156,12 @@ class FastKde { /** * Far field expansion created by the reference points in this node. */ - FarFieldExpansion farfield_expansion_; + FarFieldExpansion farfield_expansion_; /** * Local expansion stored in this node. */ - LocalExpansion local_expansion_; + LocalExpansion local_expansion_; /** Initialize the statistics */ void Init() { @@ -213,7 +213,7 @@ class FastKde { void PushDownTokens(KdeStat &left_stat, KdeStat &right_stat, double *de, - LocalExpansion + LocalExpansion *local_expansion, double *dt) { if(de != NULL) { @@ -434,9 +434,9 @@ class FastKde { KdeStat &rstat = rnode->stat(); // expansion objects - FarFieldExpansion &farfield_expansion + FarFieldExpansion &farfield_expansion = rstat.farfield_expansion_; - LocalExpansion &local_expansion + LocalExpansion &local_expansion = qstat.local_expansion_; // number of reference points @@ -460,7 +460,7 @@ class FastKde { // get the order of approximations order_farfield_to_local_ = - farfield_expansion.OrderForConvertingtoLocal(rnode->bound(), + farfield_expansion.OrderForConvertingToLocal(rnode->bound(), qnode->bound(), dsqd_range.lo, allowed_err, &actual_err); diff --git a/fastlib/u/dongryel/kde/main.cc b/fastlib/u/dongryel/kde/main.cc index 58ea1d64bf..5ca3f2f8b8 100644 --- a/fastlib/u/dongryel/kde/main.cc +++ b/fastlib/u/dongryel/kde/main.cc @@ -8,7 +8,7 @@ int main(int argc, char *argv[]) { bool do_naive = fx_param_exists(NULL, "do_naive"); - FastKde fast_kde; + FastKde fast_kde; fast_kde.Init(); fast_kde.Compute(fx_param_double(NULL, "tau", 0.1)); @@ -31,10 +31,11 @@ int main(int argc, char *argv[]) { naive_kde.ComputeMaximumRelativeError(fast_kde_results); } + /* FFTKde fft_kde; fft_kde.Init(fast_kde.get_query_dataset(), fast_kde.get_reference_dataset()); - + */ fx_done(); return 0; } diff --git a/fastlib/u/dongryel/multibody/main.cc b/fastlib/u/dongryel/multibody/main.cc index 70e8313fa0..06c47c35ce 100644 --- a/fastlib/u/dongryel/multibody/main.cc +++ b/fastlib/u/dongryel/multibody/main.cc @@ -22,7 +22,7 @@ int main(int argc, char *argv[]) if(!strcmp(kernel, "gaussianthreebody")) { fx_timer_start(NULL, "multibody"); MultitreeMultibody mtmb; + GaussianKernel, GaussianKernelAux> mtmb; mtmb.Init(bandwidth); mtmb.Compute(tau); fx_timer_stop(NULL, "multibody"); @@ -42,7 +42,7 @@ int main(int argc, char *argv[]) else if(!strcmp(kernel, "axilrodteller")) { fx_timer_start(NULL, "multibody"); MultitreeMultibody mtmb; + GaussianKernel, GaussianKernelAux> mtmb; mtmb.Init(bandwidth); mtmb.Compute(tau); fx_timer_stop(NULL, "multibody"); diff --git a/fastlib/u/dongryel/multibody/multibody_kernel.h b/fastlib/u/dongryel/multibody/multibody_kernel.h index 6a6376f05f..b0231c3c0d 100644 --- a/fastlib/u/dongryel/multibody/multibody_kernel.h +++ b/fastlib/u/dongryel/multibody/multibody_kernel.h @@ -2,7 +2,7 @@ #define MULTIBODY_KERNEL_H #include "fastlib/fastlib_int.h" -#include "u/dongryel/series_expansion/kernel_derivative.h" +#include "u/dongryel/series_expansion/kernel_aux.h" class GaussianThreeBodyKernel { diff --git a/fastlib/u/dongryel/proximity_project/pca_tree.h b/fastlib/u/dongryel/proximity_project/pca_tree.h index a67be2666d..58234be2e1 100644 --- a/fastlib/u/dongryel/proximity_project/pca_tree.h +++ b/fastlib/u/dongryel/proximity_project/pca_tree.h @@ -135,7 +135,7 @@ class PCAStat { // extract the relevant part of the dataset and mean-center it ExtractSubMatrix(dataset, start, count, orig_mean_centered_); ComputeColumnMeanVector(orig_mean_centered_, means_); - SubtractVectorFromMatrix(orig_mean_centered_, means_, orig_mean_centered_); + //SubtractVectorFromMatrix(orig_mean_centered_, means_, orig_mean_centered_); // compute PCA on the extracted submatrix Matrix U, VT; @@ -144,7 +144,7 @@ class PCAStat { // reduce the dimension in half Matrix U_trunc; - int new_dimension = U.n_cols(); + int new_dimension = U.n_cols() / 2; U_trunc.Init(new_dimension, U.n_rows()); for(index_t i = 0; i < new_dimension; i++) { Vector s; @@ -280,8 +280,6 @@ class PCAStat { right_pca_transformed.PrintDebug(); pca_transformed_.PrintDebug(); - printf("Checking!\n"); - Init(dataset, start, count); exit(0); } diff --git a/fastlib/u/dongryel/series_expansion/build.py b/fastlib/u/dongryel/series_expansion/build.py index d71ea13547..382c163f48 100644 --- a/fastlib/u/dongryel/series_expansion/build.py +++ b/fastlib/u/dongryel/series_expansion/build.py @@ -3,7 +3,7 @@ librule( name = "series_expansion", # this line can be safely omitted sources = ["series_expansion_aux.cc"], # files that must be compiled headers = ["farfield_expansion.h", - "kernel_derivative.h", + "kernel_aux.h", "local_expansion.h", "series_expansion_aux.h"], # include files part of the 'lib' deplibs = ["fastlib:fastlib_int"] # depends on fastlib core diff --git a/fastlib/u/dongryel/series_expansion/farfield_expansion.h b/fastlib/u/dongryel/series_expansion/farfield_expansion.h index 4d8799c65d..e504d23620 100644 --- a/fastlib/u/dongryel/series_expansion/farfield_expansion.h +++ b/fastlib/u/dongryel/series_expansion/farfield_expansion.h @@ -10,16 +10,16 @@ #include #include "fastlib/fastlib.h" -#include "kernel_derivative.h" +#include "kernel_aux.h" #include "series_expansion_aux.h" -template +template class LocalExpansion; /** * Far field expansion class */ -template +template class FarFieldExpansion { FORBID_COPY(FarFieldExpansion); @@ -27,7 +27,7 @@ class FarFieldExpansion { typedef TKernel Kernel; - typedef TKernelDerivative KernelDerivative; + typedef TKernelAux KernelAux; private: @@ -46,8 +46,8 @@ class FarFieldExpansion { /** precomputed quantities */ SeriesExpansionAux *sea_; - /** Derivative computer based on the kernel passed in */ - KernelDerivative kd_; + /** auxilirary methods for the kernel (derivative, truncation error bound) */ + KernelAux ka_; public: @@ -116,8 +116,8 @@ class FarFieldExpansion { */ double MixField(const Matrix &data, int node1_begin, int node1_end, int node2_begin, int node2_end, - const FarFieldExpansion &fe2, - const FarFieldExpansion &fe3, + const FarFieldExpansion &fe2, + const FarFieldExpansion &fe3, int order2, int order3) const; /** @@ -125,8 +125,8 @@ class FarFieldExpansion { * expansions */ double ConvolveField - (const FarFieldExpansion &fe2, - const FarFieldExpansion &fe3, + (const FarFieldExpansion &fe2, + const FarFieldExpansion &fe3, int order1, int order2, int order3) const; /** @@ -155,7 +155,7 @@ class FarFieldExpansion { * @return the minimum approximation order required for the error, * -1 if approximation up to the maximum order is not possible */ - int OrderForConvertingtoLocal(const DHrectBound<2> &far_field_region, + int OrderForConvertingToLocal(const DHrectBound<2> &far_field_region, const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, double required_bound, @@ -177,12 +177,12 @@ class FarFieldExpansion { * are added up to the passed-in local expansion coefficients. */ void TranslateToLocal - (LocalExpansion &se); + (LocalExpansion &se); }; -template -void FarFieldExpansion::AccumulateCoeffs +template +void FarFieldExpansion::AccumulateCoeffs (const Matrix& data, const Vector& weights, int begin, int end, int order) { @@ -192,7 +192,7 @@ void FarFieldExpansion::AccumulateCoeffs int r, i, j, k, t, tail; Vector heads; Vector x_r; - double bandwidth_factor = kd_.BandwidthFactor(kernel_.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(kernel_.bandwidth_sq()); // initialize temporary variables tmp.Init(total_num_coeffs); @@ -259,8 +259,8 @@ void FarFieldExpansion::AccumulateCoeffs } } -template -void FarFieldExpansion::RefineCoeffs +template +void FarFieldExpansion::RefineCoeffs (const Matrix& data, const Vector& weights, int begin, int end, int order) { @@ -276,7 +276,7 @@ void FarFieldExpansion::RefineCoeffs double tmp; int r, i, j; Vector x_r; - double bandwidth_factor = kd_.BandwidthFactor(kernel_.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(kernel_.bandwidth_sq()); // initialize temporary variables x_r.Init(dim); @@ -335,8 +335,8 @@ void FarFieldExpansion::RefineCoeffs } } -template -double FarFieldExpansion:: +template +double FarFieldExpansion:: EvaluateField(Matrix* data, int row_num, Vector* x_q) const { // dimension @@ -346,7 +346,7 @@ double FarFieldExpansion:: int total_num_coeffs = sea_->get_total_num_coeffs(order_); // square root times bandwidth - double bandwidth_factor = kd_.BandwidthFactor(kernel_.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(kernel_.bandwidth_sq()); // the evaluated sum double pos_multipole_sum = 0; @@ -377,12 +377,12 @@ double FarFieldExpansion:: } // compute deriative maps based on coordinate difference. - kd_.ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map); + ka_.ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map); // compute h_{\alpha}((x_q - x_R)/sqrt(2h^2)) ((x_r - x_R)/h)^{\alpha} for(index_t j = 0; j < total_num_coeffs; j++) { ArrayList mapping = sea_->get_multiindex(j); - double arrtmp = kd_.ComputePartialDerivative(derivative_map, mapping); + double arrtmp = ka_.ComputePartialDerivative(derivative_map, mapping); double prod = coeffs_[j] * arrtmp; if(prod > 0) { @@ -397,17 +397,17 @@ double FarFieldExpansion:: return multipole_sum; } -template - double FarFieldExpansion::MixField +template + double FarFieldExpansion::MixField (const Matrix &data, int node1_begin, int node1_end, int node2_begin, int node2_end, - const FarFieldExpansion &fe2, - const FarFieldExpansion &fe3, + const FarFieldExpansion &fe2, + const FarFieldExpansion &fe3, int order2, int order3) const { // bandwidth factor and multiindex mapping stuffs double result; - double bandwidth_factor = kd_.BandwidthFactor(bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(bandwidth_sq()); const ArrayList *multiindex_mapping = sea_->get_multiindex_mapping(); const ArrayList *lower_mapping_index = sea_->get_lower_mapping_index(); @@ -454,8 +454,8 @@ template xI_xK[d] = (center_[d] - xK_center[d]) / bandwidth_factor; xJ_xK[d] = (xJ_center[d] - xK_center[d]) / bandwidth_factor; } - kd_.ComputeDirectionalDerivatives(xI_xK, derivative_map_beta); - kd_.ComputeDirectionalDerivatives(xJ_xK, derivative_map_gamma); + ka_.ComputeDirectionalDerivatives(xI_xK, derivative_map_beta); + ka_.ComputeDirectionalDerivatives(xJ_xK, derivative_map_gamma); // inverse factorials Vector inv_multiindex_factorials; @@ -481,7 +481,7 @@ template ArrayList beta_mapping = multiindex_mapping[beta]; ArrayList lower_mappings_for_beta = lower_mapping_index[beta]; - double beta_derivative = kd_.ComputePartialDerivative + double beta_derivative = ka_.ComputePartialDerivative (derivative_map_beta, beta_mapping); for(index_t nu = 0; nu < lower_mappings_for_beta.size(); nu++) { @@ -499,7 +499,7 @@ template ArrayList gamma_mapping = multiindex_mapping[gamma]; ArrayList lower_mappings_for_gamma = lower_mapping_index[gamma]; - double gamma_derivative = kd_.ComputePartialDerivative + double gamma_derivative = ka_.ComputePartialDerivative (derivative_map_gamma, gamma_mapping); for(index_t eta = 0; eta < lower_mappings_for_gamma.size(); @@ -577,15 +577,15 @@ template return sum; } -template -double FarFieldExpansion::ConvolveField - (const FarFieldExpansion &fe2, - const FarFieldExpansion &fe3, +template +double FarFieldExpansion::ConvolveField + (const FarFieldExpansion &fe2, + const FarFieldExpansion &fe3, int order1, int order2, int order3) const { // bandwidth factor and multiindex mapping stuffs double result; - double bandwidth_factor = kd_.BandwidthFactor(bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(bandwidth_sq()); const ArrayList *multiindex_mapping = sea_->get_multiindex_mapping(); const ArrayList *lower_mapping_index = sea_->get_lower_mapping_index(); @@ -642,9 +642,9 @@ double FarFieldExpansion::ConvolveField 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); - kd_.ComputeDirectionalDerivatives(xJ_xK, derivative_map_gamma); + ka_.ComputeDirectionalDerivatives(xI_xJ, derivative_map_alpha); + ka_.ComputeDirectionalDerivatives(xI_xK, derivative_map_beta); + ka_.ComputeDirectionalDerivatives(xJ_xK, derivative_map_gamma); // inverse factorials Vector inv_multiindex_factorials; @@ -655,7 +655,7 @@ double FarFieldExpansion::ConvolveField ArrayList alpha_mapping = multiindex_mapping[alpha]; ArrayList lower_mappings_for_alpha = lower_mapping_index[alpha]; - double alpha_derivative = kd_.ComputePartialDerivative + double alpha_derivative = ka_.ComputePartialDerivative (derivative_map_alpha, alpha_mapping); for(index_t mu = 0; mu < lower_mappings_for_alpha.size(); mu++) { @@ -672,7 +672,7 @@ double FarFieldExpansion::ConvolveField ArrayList beta_mapping = multiindex_mapping[beta]; ArrayList lower_mappings_for_beta = lower_mapping_index[beta]; - double beta_derivative = kd_.ComputePartialDerivative + double beta_derivative = ka_.ComputePartialDerivative (derivative_map_beta, beta_mapping); for(index_t nu = 0; nu < lower_mappings_for_beta.size(); nu++) { @@ -691,7 +691,7 @@ double FarFieldExpansion::ConvolveField ArrayList gamma_mapping = multiindex_mapping[gamma]; ArrayList lower_mappings_for_gamma = lower_mapping_index[gamma]; - double gamma_derivative = kd_.ComputePartialDerivative + double gamma_derivative = ka_.ComputePartialDerivative (derivative_map_gamma, gamma_mapping); for(index_t eta = 0; eta < lower_mappings_for_gamma.size(); @@ -763,8 +763,8 @@ double FarFieldExpansion::ConvolveField return sum; } -template - void FarFieldExpansion::Init +template + void FarFieldExpansion::Init (double bandwidth, const Vector& center, SeriesExpansionAux *sea) { // copy kernel type, center, and bandwidth squared @@ -773,13 +773,18 @@ template order_ = -1; sea_ = sea; + // pass in the pointer to the kernel and the series expansion auxiliary + // object + ka_.kernel_ = &kernel_; + ka_.sea_ = sea_; + // initialize coefficient array coeffs_.Init(sea_->get_max_total_num_coeffs()); coeffs_.SetZero(); } -template - void FarFieldExpansion::Init +template + void FarFieldExpansion::Init (double bandwidth, SeriesExpansionAux *sea) { // copy kernel type, center, and bandwidth squared @@ -789,144 +794,42 @@ template order_ = -1; sea_ = sea; + // pass in the pointer to the kernel and the series expansion auxiliary + // object + ka_.kernel_ = &kernel_; + ka_.sea_ = sea_; + // initialize coefficient array coeffs_.Init(sea_->get_max_total_num_coeffs()); coeffs_.SetZero(); } -template - int FarFieldExpansion::OrderForEvaluating +template + int FarFieldExpansion::OrderForEvaluating (const DHrectBound<2> &far_field_region, double min_dist_sqd_regions, double max_error, double *actual_error) const { - double frontfactor = - exp(-min_dist_sqd_regions / (4 * kernel_.bandwidth_sq())); - double widest_width = 0; - int dim = far_field_region.dim(); - int max_order = sea_->get_max_order(); - - // find out the widest dimension and its length - for(index_t d = 0; d < dim; d++) { - DRange range = far_field_region.get(d); - widest_width = max(widest_width, range.width()); - } - - double two_bandwidth = 2 * sqrt(kernel_.bandwidth_sq()); - double r = widest_width / two_bandwidth; - - // This is not really necessary for O(D^p) expansion, but it is for - // speeding up the convergence of the Taylor expansion. - if(r >= 1.0) - return -1; - - double r_raised_to_p_alpha = 1.0; - double ret; - int p_alpha = 0; - double floor_fact, ceil_fact; - int remainder; - - do { - - if(p_alpha > max_order - 1) - return -1; - - r_raised_to_p_alpha *= r; - - floor_fact = - sea_->factorial((int)floor(((double) p_alpha) / ((double) dim))); - ceil_fact = - sea_->factorial((int)ceil(((double) p_alpha) / ((double) dim))); - - if(floor_fact < 0.0 || ceil_fact < 0.0) - return -1; - - remainder = p_alpha % dim; - - ret = frontfactor * - (sea_->get_total_num_coeffs(p_alpha + 1) - - sea_->get_total_num_coeffs(p_alpha)) * r_raised_to_p_alpha / - sqrt(pow(floor_fact, dim - remainder) * pow(ceil_fact, remainder)); - - if(ret > max_error) { - p_alpha++; - } - else { - break; - } - } while(1); - - *actual_error = ret; - return p_alpha; + return ka_.OrderForEvaluatingFarField(far_field_region, + min_dist_sqd_regions, max_error, + actual_error); } -template - int FarFieldExpansion:: - OrderForConvertingtoLocal(const DHrectBound<2> &far_field_region, +template + int FarFieldExpansion:: + OrderForConvertingToLocal(const DHrectBound<2> &far_field_region, const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, double max_error, double *actual_error) const { - double max_ref_length = 0; - double max_query_length = 0; - int dim = sea_->get_dimension(); - - for(index_t i = 0; i < dim; i++) { - DRange far_field_range = far_field_region.get(i); - DRange local_range = local_field_region.get(i); - max_ref_length = max(max_ref_length, far_field_range.width()); - max_query_length = max(max_query_length, local_range.width()); - } - - double two_times_bandwidth = sqrt(kernel_.bandwidth_sq()) * 2; - double r_R = max_ref_length / two_times_bandwidth; - double r_Q = max_query_length / two_times_bandwidth; - double sqrt_two_r_R = sqrt(2.0) * r_R; - double sqrt_two_r_Q = sqrt(2.0) * r_Q; - - if(sqrt_two_r_R >= 1.0 || sqrt_two_r_Q >= 1.0) { - return -1; - } - - int p_alpha = -1; - double sqrt_two_r_R_raised_to_p = 1.0; - double r_Q_raised_to_p = 1.0; - int remainder; - double ret2; - double frontfactor = - exp(-min_dist_sqd_regions / (4.0 * kernel_.bandwidth_sq())); - double floor_fact, ceil_fact; - - do { - p_alpha++; - - r_Q_raised_to_p *= r_Q; - sqrt_two_r_R_raised_to_p *= sqrt_two_r_R; - floor_fact = - sea_->factorial((int) floor((double) p_alpha / (double) dim)); - ceil_fact = - sea_->factorial((int) ceil((double)p_alpha / (double)dim)); - - if(floor_fact < 0 || ceil_fact < 0 || p_alpha > sea_->get_max_order() - 1) - return -1; - - remainder = p_alpha % dim; - - ret2 = (sea_->get_total_num_coeffs(p_alpha + 1) - - sea_->get_total_num_coeffs(p_alpha)) - / sqrt(pow(floor_fact, dim - remainder) * - pow(ceil_fact, remainder)); - ret2 *= (r_Q_raised_to_p + sqrt_two_r_R_raised_to_p * - sea_->get_total_num_coeffs(p_alpha)) * frontfactor; - - } while(ret2 >= max_error); - - *actual_error = ret2; - return p_alpha; + return ka_.OrderForConvertingFromFarFieldToLocal(far_field_region, + local_field_region, + min_dist_sqd_regions, + max_error, actual_error); } -template - void FarFieldExpansion::PrintDebug +template + void FarFieldExpansion::PrintDebug (const char *name, FILE *stream) const { @@ -975,11 +878,11 @@ template fprintf(stream, "\n"); } -template - void FarFieldExpansion::TranslateFromFarField +template + void FarFieldExpansion::TranslateFromFarField (const FarFieldExpansion &se) { - double bandwidth_factor = kd_.BandwidthFactor(se.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(se.bandwidth_sq()); int dim = sea_->get_dimension(); int order = se.get_order(); int total_num_coeffs = sea_->get_total_num_coeffs(order); @@ -1065,9 +968,9 @@ template } // end of j-loop } -template -void FarFieldExpansion::TranslateToLocal - (LocalExpansion &se) { +template +void FarFieldExpansion::TranslateToLocal + (LocalExpansion &se) { Vector pos_arrtmp, neg_arrtmp; Matrix derivative_map; @@ -1078,7 +981,7 @@ void FarFieldExpansion::TranslateToLocal int dimension = sea_->get_dimension(); int total_num_coeffs = sea_->get_total_num_coeffs(order_); int limit; - double bandwidth_factor = kd_.BandwidthFactor(se.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(se.bandwidth_sq()); // get center and coefficients for local expansion local_center.Alias(se.get_center()); @@ -1103,7 +1006,7 @@ void FarFieldExpansion::TranslateToLocal } // compute required partial derivatives - kd_.ComputeDirectionalDerivatives(cent_diff, derivative_map); + ka_.ComputeDirectionalDerivatives(cent_diff, derivative_map); ArrayList beta_plus_alpha; beta_plus_alpha.Init(dimension); @@ -1119,7 +1022,7 @@ void FarFieldExpansion::TranslateToLocal beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d]; } double derivative_factor = - kd_.ComputePartialDerivative(derivative_map, beta_plus_alpha); + ka_.ComputePartialDerivative(derivative_map, beta_plus_alpha); double prod = coeffs_[k] * derivative_factor; diff --git a/fastlib/u/dongryel/series_expansion/kernel_aux.h b/fastlib/u/dongryel/series_expansion/kernel_aux.h new file mode 100644 index 0000000000..2b0f7bbdaf --- /dev/null +++ b/fastlib/u/dongryel/series_expansion/kernel_aux.h @@ -0,0 +1,545 @@ +/** + * @file kernel_aux.h + * + * The header file for the class for computing auxiliary stuffs for the kernel + * functions (derivative, truncation error bound) + */ + +#ifndef KERNEL_AUX +#define KERNEL_AUX + +#include "fastlib/fastlib.h" + +#include "series_expansion_aux.h" + +/** + * Auxiliary computer class for Gaussian kernel + */ +class GaussianKernelAux { + FORBID_COPY(GaussianKernelAux); + + public: + + /** pointer to the Gaussian kernel */ + GaussianKernel *kernel_; + + /** pointer to the series expansion auxiliary object */ + SeriesExpansionAux *sea_; + + GaussianKernelAux() {} + + ~GaussianKernelAux() {} + + double BandwidthFactor(double bandwidth_sq) const { + return sqrt(2 * bandwidth_sq); + } + + void ComputeDirectionalDerivatives(const Vector &x, + Matrix &derivative_map) const { + + int dim = derivative_map.n_rows(); + int order = derivative_map.n_cols() - 1; + + // precompute necessary Hermite polynomials based on coordinate difference + for(index_t d = 0; d < dim; d++) { + + double coord_div_band = x[d]; + double d2 = 2 * coord_div_band; + double facj = exp(-coord_div_band * coord_div_band); + + derivative_map.set(d, 0, facj); + + if(order > 0) { + + derivative_map.set(d, 1, d2 * facj); + + if(order > 1) { + for(index_t k = 1; k < order; k++) { + int k2 = k * 2; + derivative_map.set(d, k + 1, d2 * derivative_map.get(d, k) - + k2 * derivative_map.get(d, k - 1)); + } + } + } + } // end of looping over each dimension + } + + double ComputePartialDerivative(const Matrix &derivative_map, + ArrayList mapping) const { + + double partial_derivative = 1.0; + + for(index_t d = 0; d < mapping.size(); d++) { + partial_derivative *= derivative_map.get(d, mapping[d]); + } + return partial_derivative; + } + + int OrderForEvaluatingFarField + (const DHrectBound<2> &far_field_region, double min_dist_sqd_regions, + double max_error, double *actual_error) const { + + double frontfactor = + exp(-min_dist_sqd_regions / (4 * kernel_->bandwidth_sq())); + double widest_width = 0; + int dim = far_field_region.dim(); + int max_order = sea_->get_max_order(); + + // find out the widest dimension and its length + for(index_t d = 0; d < dim; d++) { + DRange range = far_field_region.get(d); + widest_width = max(widest_width, range.width()); + } + + double two_bandwidth = 2 * sqrt(kernel_->bandwidth_sq()); + double r = widest_width / two_bandwidth; + + // This is not really necessary for O(D^p) expansion, but it is for + // speeding up the convergence of the Taylor expansion. + if(r >= 1.0) { + return -1; + } + + double r_raised_to_p_alpha = 1.0; + double ret; + int p_alpha = 0; + double floor_fact, ceil_fact; + int remainder; + + do { + + if(p_alpha > max_order - 1) + return -1; + + r_raised_to_p_alpha *= r; + + floor_fact = + sea_->factorial((int)floor(((double) p_alpha) / ((double) dim))); + ceil_fact = + sea_->factorial((int)ceil(((double) p_alpha) / ((double) dim))); + + if(floor_fact < 0.0 || ceil_fact < 0.0) { + return -1; + } + + remainder = p_alpha % dim; + + ret = frontfactor * + (sea_->get_total_num_coeffs(p_alpha + 1) - + sea_->get_total_num_coeffs(p_alpha)) * r_raised_to_p_alpha / + sqrt(pow(floor_fact, dim - remainder) * pow(ceil_fact, remainder)); + + if(ret > max_error) { + p_alpha++; + } + else { + break; + } + } while(1); + + *actual_error = ret; + return p_alpha; + } + + int OrderForConvertingFromFarFieldToLocal + (const DHrectBound<2> &far_field_region, + const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, + double max_error, double *actual_error) const { + + double max_ref_length = 0; + double max_query_length = 0; + int dim = sea_->get_dimension(); + + for(index_t i = 0; i < dim; i++) { + DRange far_field_range = far_field_region.get(i); + DRange local_range = local_field_region.get(i); + max_ref_length = max(max_ref_length, far_field_range.width()); + max_query_length = max(max_query_length, local_range.width()); + } + + double two_times_bandwidth = sqrt(kernel_->bandwidth_sq()) * 2; + double r_R = max_ref_length / two_times_bandwidth; + double r_Q = max_query_length / two_times_bandwidth; + double sqrt_two_r_R = sqrt(2.0) * r_R; + double sqrt_two_r_Q = sqrt(2.0) * r_Q; + + if(sqrt_two_r_R >= 1.0 || sqrt_two_r_Q >= 1.0) { + return -1; + } + + int p_alpha = -1; + double sqrt_two_r_R_raised_to_p = 1.0; + double r_Q_raised_to_p = 1.0; + int remainder; + double ret2; + double frontfactor = + exp(-min_dist_sqd_regions / (4.0 * kernel_->bandwidth_sq())); + double floor_fact, ceil_fact; + + do { + p_alpha++; + + r_Q_raised_to_p *= r_Q; + sqrt_two_r_R_raised_to_p *= sqrt_two_r_R; + floor_fact = + sea_->factorial((int) floor((double) p_alpha / (double) dim)); + ceil_fact = + sea_->factorial((int) ceil((double)p_alpha / (double)dim)); + + if(floor_fact < 0 || ceil_fact < 0 || + p_alpha > sea_->get_max_order() - 1) { + return -1; + } + + remainder = p_alpha % dim; + + ret2 = (sea_->get_total_num_coeffs(p_alpha + 1) - + sea_->get_total_num_coeffs(p_alpha)) + / sqrt(pow(floor_fact, dim - remainder) * + pow(ceil_fact, remainder)); + ret2 *= (r_Q_raised_to_p + sqrt_two_r_R_raised_to_p * + sea_->get_total_num_coeffs(p_alpha)) * frontfactor; + + } while(ret2 >= max_error); + + *actual_error = ret2; + return p_alpha; + } + + int OrderForEvaluatingLocal + (const DHrectBound<2> &local_region, double min_dist_sqd_regions, + double max_error, double *actual_error) const { + + double frontfactor = + exp(-min_dist_sqd_regions / (4 * kernel_->bandwidth_sq())); + double widest_width = 0; + int dim = local_region.dim(); + int max_order = sea_->get_max_order(); + + // find out the widest dimension and its length + for(index_t d = 0; d < dim; d++) { + DRange range = local_region.get(d); + widest_width = max(widest_width, range.width()); + } + + double two_bandwidth = 2 * sqrt(kernel_->bandwidth_sq()); + double r = widest_width / two_bandwidth; + + // This is not really necessary for O(D^p) expansion, but it is for + // speeding up the convergence of the Taylor expansion. + if(r >= 1.0) + return -1; + + double r_raised_to_p_alpha = 1.0; + double ret; + int p_alpha = 0; + double floor_fact, ceil_fact; + int remainder; + + do { + + if(p_alpha > max_order - 1) + return -1; + + r_raised_to_p_alpha *= r; + + floor_fact = + sea_->factorial((int)floor(((double) p_alpha) / ((double) dim))); + ceil_fact = + sea_->factorial((int)ceil(((double) p_alpha) / ((double) dim))); + + if(floor_fact < 0.0 || ceil_fact < 0.0) + return -1; + + remainder = p_alpha % dim; + + ret = frontfactor * + (sea_->get_total_num_coeffs(p_alpha + 1) - + sea_->get_total_num_coeffs(p_alpha)) * r_raised_to_p_alpha / + sqrt(pow(floor_fact, dim - remainder) * pow(ceil_fact, remainder)); + + if(ret > max_error) { + p_alpha++; + } + else { + break; + } + } while(1); + + *actual_error = ret; + return p_alpha; + } +}; + +/** + * Auxilairy computer class for Epanechnikov kernel + */ +class EpanKernelAux { + FORBID_COPY(EpanKernelAux); + + public: + + EpanKernel *kernel_; + + SeriesExpansionAux *sea_; + + EpanKernelAux() {} + + ~EpanKernelAux() {} + + double BandwidthFactor(double bandwidth_sq) const { + return sqrt(bandwidth_sq); + } + + void ComputeDirectionalDerivatives(const Vector &x, + Matrix &derivative_map) const { + + int dim = derivative_map.n_rows(); + int order = derivative_map.n_cols() - 1; + + // precompute necessary Hermite polynomials based on coordinate difference + for(index_t d = 0; d < dim; d++) { + + double coord_div_band = x[d]; + + derivative_map.set(d, 0, coord_div_band * coord_div_band); + + if(order > 0) { + derivative_map.set(d, 1, 2 * coord_div_band); + + if(order > 1) { + derivative_map.set(d, 2, -2); + + for(index_t k = 3; k <= order; k++) { + derivative_map.set(d, k, 0); + } + } + } + + } // end of looping over each dimension + } + + double ComputePartialDerivative(const Matrix &derivative_map, + ArrayList mapping) const { + + int nonzero_count = 0; + int nonzero_index = 0; + + for(index_t d = 0; d < mapping.size(); d++) { + if(mapping[d] > 0) { + nonzero_count++; + nonzero_index = d; + } + + if(nonzero_count > 1) { + return 0; + } + } + if(nonzero_count == 0) { + double prod = 0; + for(index_t d = 0; d < mapping.size(); d++) { + prod += derivative_map.get(d, 0); + } + return 1.0 - prod; + } + + return derivative_map.get(nonzero_index, mapping[nonzero_index]); + } + + int OrderForEvaluatingFarField + (const DHrectBound<2> &far_field_region, double min_dist_sqd_regions, + double max_error, double *actual_error) const { + + double frontfactor = + exp(-min_dist_sqd_regions / (4 * kernel_->bandwidth_sq())); + double widest_width = 0; + int dim = far_field_region.dim(); + int max_order = sea_->get_max_order(); + + // find out the widest dimension and its length + for(index_t d = 0; d < dim; d++) { + DRange range = far_field_region.get(d); + widest_width = max(widest_width, range.width()); + } + + double two_bandwidth = 2 * sqrt(kernel_->bandwidth_sq()); + double r = widest_width / two_bandwidth; + + // This is not really necessary for O(D^p) expansion, but it is for + // speeding up the convergence of the Taylor expansion. + if(r >= 1.0) { + return -1; + } + + double r_raised_to_p_alpha = 1.0; + double ret; + int p_alpha = 0; + double floor_fact, ceil_fact; + int remainder; + + do { + + if(p_alpha > max_order - 1) + return -1; + + r_raised_to_p_alpha *= r; + + floor_fact = + sea_->factorial((int)floor(((double) p_alpha) / ((double) dim))); + ceil_fact = + sea_->factorial((int)ceil(((double) p_alpha) / ((double) dim))); + + if(floor_fact < 0.0 || ceil_fact < 0.0) { + return -1; + } + + remainder = p_alpha % dim; + + ret = frontfactor * + (sea_->get_total_num_coeffs(p_alpha + 1) - + sea_->get_total_num_coeffs(p_alpha)) * r_raised_to_p_alpha / + sqrt(pow(floor_fact, dim - remainder) * pow(ceil_fact, remainder)); + + if(ret > max_error) { + p_alpha++; + } + else { + break; + } + } while(1); + + *actual_error = ret; + return p_alpha; + } + + int OrderForConvertingFromFarFieldToLocal + (const DHrectBound<2> &far_field_region, + const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, + double max_error, double *actual_error) const { + + double max_ref_length = 0; + double max_query_length = 0; + int dim = sea_->get_dimension(); + + for(index_t i = 0; i < dim; i++) { + DRange far_field_range = far_field_region.get(i); + DRange local_range = local_field_region.get(i); + max_ref_length = max(max_ref_length, far_field_range.width()); + max_query_length = max(max_query_length, local_range.width()); + } + + double two_times_bandwidth = sqrt(kernel_->bandwidth_sq()) * 2; + double r_R = max_ref_length / two_times_bandwidth; + double r_Q = max_query_length / two_times_bandwidth; + double sqrt_two_r_R = sqrt(2.0) * r_R; + double sqrt_two_r_Q = sqrt(2.0) * r_Q; + + if(sqrt_two_r_R >= 1.0 || sqrt_two_r_Q >= 1.0) { + return -1; + } + + int p_alpha = -1; + double sqrt_two_r_R_raised_to_p = 1.0; + double r_Q_raised_to_p = 1.0; + int remainder; + double ret2; + double frontfactor = + exp(-min_dist_sqd_regions / (4.0 * kernel_->bandwidth_sq())); + double floor_fact, ceil_fact; + + do { + p_alpha++; + + r_Q_raised_to_p *= r_Q; + sqrt_two_r_R_raised_to_p *= sqrt_two_r_R; + floor_fact = + sea_->factorial((int) floor((double) p_alpha / (double) dim)); + ceil_fact = + sea_->factorial((int) ceil((double)p_alpha / (double)dim)); + + if(floor_fact < 0 || ceil_fact < 0 || + p_alpha > sea_->get_max_order() - 1) { + return -1; + } + + remainder = p_alpha % dim; + + ret2 = (sea_->get_total_num_coeffs(p_alpha + 1) - + sea_->get_total_num_coeffs(p_alpha)) + / sqrt(pow(floor_fact, dim - remainder) * + pow(ceil_fact, remainder)); + ret2 *= (r_Q_raised_to_p + sqrt_two_r_R_raised_to_p * + sea_->get_total_num_coeffs(p_alpha)) * frontfactor; + + } while(ret2 >= max_error); + + *actual_error = ret2; + return p_alpha; + } + + int OrderForEvaluatingLocal + (const DHrectBound<2> &local_region, double min_dist_sqd_regions, + double max_error, double *actual_error) const { + + double frontfactor = + exp(-min_dist_sqd_regions / (4 * kernel_->bandwidth_sq())); + double widest_width = 0; + int dim = local_region.dim(); + int max_order = sea_->get_max_order(); + + // find out the widest dimension and its length + for(index_t d = 0; d < dim; d++) { + DRange range = local_region.get(d); + widest_width = max(widest_width, range.width()); + } + + double two_bandwidth = 2 * sqrt(kernel_->bandwidth_sq()); + double r = widest_width / two_bandwidth; + + // This is not really necessary for O(D^p) expansion, but it is for + // speeding up the convergence of the Taylor expansion. + if(r >= 1.0) + return -1; + + double r_raised_to_p_alpha = 1.0; + double ret; + int p_alpha = 0; + double floor_fact, ceil_fact; + int remainder; + + do { + + if(p_alpha > max_order - 1) + return -1; + + r_raised_to_p_alpha *= r; + + floor_fact = + sea_->factorial((int)floor(((double) p_alpha) / ((double) dim))); + ceil_fact = + sea_->factorial((int)ceil(((double) p_alpha) / ((double) dim))); + + if(floor_fact < 0.0 || ceil_fact < 0.0) + return -1; + + remainder = p_alpha % dim; + + ret = frontfactor * + (sea_->get_total_num_coeffs(p_alpha + 1) - + sea_->get_total_num_coeffs(p_alpha)) * r_raised_to_p_alpha / + sqrt(pow(floor_fact, dim - remainder) * pow(ceil_fact, remainder)); + + if(ret > max_error) { + p_alpha++; + } + else { + break; + } + } while(1); + + *actual_error = ret; + return p_alpha; + } +}; + +#endif diff --git a/fastlib/u/dongryel/series_expansion/kernel_derivative.h b/fastlib/u/dongryel/series_expansion/kernel_derivative.h deleted file mode 100644 index 2136daeee2..0000000000 --- a/fastlib/u/dongryel/series_expansion/kernel_derivative.h +++ /dev/null @@ -1,144 +0,0 @@ -/** - * @file kernel_derivative.h - * - * The header file for the class for computing derivatives of kernel functions - */ - -#ifndef KERNEL_DERIVATIVE -#define KERNEL_DERIVATIVE - -#include "fastlib/fastlib.h" - -/** - * Derivative computer class for Gaussian kernel - */ -class GaussianKernelDerivative { - FORBID_COPY(GaussianKernelDerivative); - - public: - - GaussianKernelDerivative() {} - - ~GaussianKernelDerivative() {} - - double BandwidthFactor(double bandwidth_sq) const { - return sqrt(2 * bandwidth_sq); - } - - void ComputeDirectionalDerivatives(const Vector &x, - Matrix &derivative_map) const { - - int dim = derivative_map.n_rows(); - int order = derivative_map.n_cols() - 1; - - // precompute necessary Hermite polynomials based on coordinate difference - for(index_t d = 0; d < dim; d++) { - - double coord_div_band = x[d]; - double d2 = 2 * coord_div_band; - double facj = exp(-coord_div_band * coord_div_band); - - derivative_map.set(d, 0, facj); - - if(order > 0) { - - derivative_map.set(d, 1, d2 * facj); - - if(order > 1) { - for(index_t k = 1; k < order; k++) { - int k2 = k * 2; - derivative_map.set(d, k + 1, d2 * derivative_map.get(d, k) - - k2 * derivative_map.get(d, k - 1)); - } - } - } - } // end of looping over each dimension - } - - double ComputePartialDerivative(const Matrix &derivative_map, - ArrayList mapping) const { - - double partial_derivative = 1.0; - - for(index_t d = 0; d < mapping.size(); d++) { - partial_derivative *= derivative_map.get(d, mapping[d]); - } - return partial_derivative; - } - -}; - -/** - * Derivative computer class for Epanechnikov kernel - */ -class EpanKernelDerivative { - FORBID_COPY(EpanKernelDerivative); - - public: - - EpanKernelDerivative() {} - - ~EpanKernelDerivative() {} - - double BandwidthFactor(double bandwidth_sq) const { - return sqrt(bandwidth_sq); - } - - void ComputeDirectionalDerivatives(const Vector &x, - Matrix &derivative_map) const { - - int dim = derivative_map.n_rows(); - int order = derivative_map.n_cols() - 1; - - // precompute necessary Hermite polynomials based on coordinate difference - for(index_t d = 0; d < dim; d++) { - - double coord_div_band = x[d]; - - derivative_map.set(d, 0, coord_div_band * coord_div_band); - - if(order > 0) { - derivative_map.set(d, 1, 2 * coord_div_band); - - if(order > 1) { - derivative_map.set(d, 2, -2); - - for(index_t k = 3; k <= order; k++) { - derivative_map.set(d, k, 0); - } - } - } - - } // end of looping over each dimension - } - - double ComputePartialDerivative(const Matrix &derivative_map, - ArrayList mapping) const { - - int nonzero_count = 0; - int nonzero_index = 0; - - for(index_t d = 0; d < mapping.size(); d++) { - if(mapping[d] > 0) { - nonzero_count++; - nonzero_index = d; - } - - if(nonzero_count > 1) { - return 0; - } - } - if(nonzero_count == 0) { - double prod = 0; - for(index_t d = 0; d < mapping.size(); d++) { - prod += derivative_map.get(d, 0); - } - return 1.0 - prod; - } - - return derivative_map.get(nonzero_index, mapping[nonzero_index]); - } - -}; - -#endif diff --git a/fastlib/u/dongryel/series_expansion/local_expansion.h b/fastlib/u/dongryel/series_expansion/local_expansion.h index 33835b7870..fd4c488021 100644 --- a/fastlib/u/dongryel/series_expansion/local_expansion.h +++ b/fastlib/u/dongryel/series_expansion/local_expansion.h @@ -10,16 +10,16 @@ #include #include "fastlib/fastlib.h" -#include "kernel_derivative.h" +#include "kernel_aux.h" #include "series_expansion_aux.h" -template +template class FarFieldExpansion; /** * Local expansion class */ -template +template class LocalExpansion { FORBID_COPY(LocalExpansion); @@ -27,7 +27,7 @@ class LocalExpansion { typedef TKernel Kernel; - typedef TKernelDerivative KernelDerivative; + typedef TKernelAux KernelAux; private: @@ -46,9 +46,11 @@ class LocalExpansion { /** precomputed quantities */ SeriesExpansionAux *sea_; - /** Derivative computer based on the kernel passed in */ - KernelDerivative kd_; + /** auxiliary methods for the kernel (derivative, truncation error bound) */ + KernelAux ka_; + /** error bound computer for the kernel */ + public: LocalExpansion() {} @@ -125,19 +127,19 @@ class LocalExpansion { * The translated coefficients are added up to the ones here. */ void TranslateFromFarField - (const FarFieldExpansion &se); + (const FarFieldExpansion &se); /** * Translate to the given local expansion. The translated coefficients * are added up to the passed-in local expansion coefficients. */ void TranslateToLocal - (LocalExpansion &se); + (LocalExpansion &se); }; -template -void LocalExpansion::AccumulateCoeffs +template +void LocalExpansion::AccumulateCoeffs (const Matrix& data, const Vector& weights, int begin, int end, int order) { @@ -164,7 +166,7 @@ void LocalExpansion::AccumulateCoeffs x_r_minus_x_Q.Init(dim); // sqrt two times bandwidth - double bandwidth_factor = kd_.BandwidthFactor(kernel_.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(kernel_.bandwidth_sq()); // for each data point, for(index_t r = begin; r < end; r++) { @@ -176,12 +178,12 @@ void LocalExpansion::AccumulateCoeffs } // precompute necessary partial derivatives based on coordinate difference - kd_.ComputeDirectionalDerivatives(x_r_minus_x_Q, derivative_map); + ka_.ComputeDirectionalDerivatives(x_r_minus_x_Q, derivative_map); // compute h_{beta}((x_r - x_Q) / sqrt(2h^2)) for(index_t j = 0; j < total_num_coeffs; j++) { ArrayList mapping = sea_->get_multiindex(j); - arrtmp[j] = kd_.ComputePartialDerivative(derivative_map, mapping); + arrtmp[j] = ka_.ComputePartialDerivative(derivative_map, mapping); } for(index_t j = 0; j < total_num_coeffs; j++) { @@ -191,8 +193,8 @@ void LocalExpansion::AccumulateCoeffs } // End of looping through each reference point. } -template - void LocalExpansion::PrintDebug +template + void LocalExpansion::PrintDebug (const char *name, FILE *stream) const { @@ -231,8 +233,8 @@ template fprintf(stream, "\n"); } -template - double LocalExpansion:: +template + double LocalExpansion:: EvaluateField(Matrix* data, int row_num, Vector* x_q) const { index_t k, t, tail; @@ -247,7 +249,7 @@ template double sum = 0; // sqrt two bandwidth - double bandwidth_factor = kd_.BandwidthFactor(kernel_.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(kernel_.bandwidth_sq()); // temporary variable Vector x_Q_to_x_q; @@ -294,8 +296,8 @@ template return sum; } -template - void LocalExpansion::Init +template + void LocalExpansion::Init (double bandwidth, const Vector& center, SeriesExpansionAux *sea) { // copy kernel type, center, and bandwidth squared @@ -304,13 +306,18 @@ template order_ = 0; sea_ = sea; + // pass in the pointer to the kernel and the series expansion auxiliary + // object + ka_.kernel_ = &kernel_; + ka_.sea_ = sea_; + // initialize coefficient array coeffs_.Init(sea_->get_max_total_num_coeffs()); coeffs_.SetZero(); } -template - void LocalExpansion::Init +template + void LocalExpansion::Init (double bandwidth, SeriesExpansionAux *sea) { // copy kernel type, center, and bandwidth squared @@ -319,79 +326,28 @@ template order_ = 0; sea_ = sea; + // pass in the pointer to the kernel and the series expansion auxiliary + // object + ka_.kernel_ = &kernel_; + ka_.sea_ = sea_; + // initialize coefficient array coeffs_.Init(sea_->get_max_total_num_coeffs()); coeffs_.SetZero(); } -template - int LocalExpansion::OrderForEvaluating +template + int LocalExpansion::OrderForEvaluating (const DHrectBound<2> &local_region, double min_dist_sqd_regions, double max_error, double *actual_error) const { - - double frontfactor = - exp(-min_dist_sqd_regions / (4 * kernel_.bandwidth_sq())); - double widest_width = 0; - int dim = local_region.dim(); - int max_order = sea_->get_max_order(); - - // find out the widest dimension and its length - for(index_t d = 0; d < dim; d++) { - DRange range = local_region.get(d); - widest_width = max(widest_width, range.width()); - } - - double two_bandwidth = 2 * sqrt(kernel_.bandwidth_sq()); - double r = widest_width / two_bandwidth; - - // This is not really necessary for O(D^p) expansion, but it is for - // speeding up the convergence of the Taylor expansion. - if(r >= 1.0) - return -1; - - double r_raised_to_p_alpha = 1.0; - double ret; - int p_alpha = 0; - double floor_fact, ceil_fact; - int remainder; - - do { - - if(p_alpha > max_order - 1) - return -1; - - r_raised_to_p_alpha *= r; - - floor_fact = - sea_->factorial((int)floor(((double) p_alpha) / ((double) dim))); - ceil_fact = - sea_->factorial((int)ceil(((double) p_alpha) / ((double) dim))); - - if(floor_fact < 0.0 || ceil_fact < 0.0) - return -1; - - remainder = p_alpha % dim; - - ret = frontfactor * - (sea_->get_total_num_coeffs(p_alpha + 1) - - sea_->get_total_num_coeffs(p_alpha)) * r_raised_to_p_alpha / - sqrt(pow(floor_fact, dim - remainder) * pow(ceil_fact, remainder)); - - if(ret > max_error) { - p_alpha++; - } - else { - break; - } - } while(1); - - *actual_error = ret; - return p_alpha; + + return ka_.OrderForEvaluatingLocal(local_region, min_dist_sqd_regions, + max_error, actual_error); } -template - void LocalExpansion::TranslateFromFarField - (const FarFieldExpansion &se) { +template + void LocalExpansion::TranslateFromFarField + (const FarFieldExpansion &se) { Vector pos_arrtmp, neg_arrtmp; Matrix derivative_map; @@ -402,7 +358,7 @@ template int far_order = se.get_order(); int total_num_coeffs = sea_->get_total_num_coeffs(far_order); int limit; - double bandwidth_factor = kd_.BandwidthFactor(se.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(se.bandwidth_sq()); // get center and coefficients for far field expansion far_center.Alias(se.get_center()); @@ -427,7 +383,7 @@ template } // compute required partial derivatives - kd_.ComputeDirectionalDerivatives(cent_diff, derivative_map); + ka_.ComputeDirectionalDerivatives(cent_diff, derivative_map); ArrayList beta_plus_alpha; beta_plus_alpha.Init(dimension); @@ -443,7 +399,7 @@ template beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d]; } double derivative_factor = - kd_.ComputePartialDerivative(derivative_map, beta_plus_alpha); + ka_.ComputePartialDerivative(derivative_map, beta_plus_alpha); double prod = far_coeffs[k] * derivative_factor; @@ -462,9 +418,9 @@ template } } -template - void LocalExpansion::TranslateToLocal - (LocalExpansion &se) { +template + void LocalExpansion::TranslateToLocal + (LocalExpansion &se) { // get the center and the order and the total number of coefficients of // the expansion we are translating from. Also get coefficients we @@ -486,7 +442,7 @@ template tmp_storage.Init(dim); // sqrt two times bandwidth - double bandwidth_factor = kd_.BandwidthFactor(kernel_.bandwidth_sq()); + double bandwidth_factor = ka_.BandwidthFactor(kernel_.bandwidth_sq()); // center difference between the old center and the new one Vector center_diff; diff --git a/fastlib/u/dongryel/series_expansion/main.cc b/fastlib/u/dongryel/series_expansion/main.cc index 6d5c97effc..ed7e787f5e 100644 --- a/fastlib/u/dongryel/series_expansion/main.cc +++ b/fastlib/u/dongryel/series_expansion/main.cc @@ -6,7 +6,7 @@ */ #include "fastlib/fastlib.h" -#include "kernel_derivative.h" +#include "kernel_aux.h" #include "farfield_expansion.h" #include "local_expansion.h" #include "series_expansion_aux.h" @@ -36,7 +36,7 @@ int TestEpanKernelEvaluateFarField(const Matrix &data, const Vector &weights, evaluate_here[0] = evaluate_here[1] = 0.1; // declare expansion object - FarFieldExpansion se; + FarFieldExpansion se; // initialize expansion objects with respective center and the bandwidth se.Init(bandwidth, center, &sea); @@ -124,7 +124,7 @@ int TestEvaluateFarField(const Matrix &data, const Vector &weights, evaluate_here[0] = evaluate_here[1] = 3; // declare expansion objects at (0,0) and other centers - FarFieldExpansion se; + FarFieldExpansion se; // initialize expansion objects with respective centers and the bandwidth // squared of 0.5 @@ -180,7 +180,7 @@ int TestEvaluateLocalField(const Matrix &data, const Vector &weights, evaluate_here[0] = evaluate_here[1] = 3.5; // declare expansion objects at (0,0) and other centers - LocalExpansion se; + LocalExpansion se; // initialize expansion objects with respective centers and the bandwidth // squared of 1 @@ -247,9 +247,9 @@ int TestTransFarToFar(const Matrix &data, const Vector &weights, new_center[1] = -2; // declare expansion objects at (0,0) and other centers - FarFieldExpansion se; - FarFieldExpansion se_translated; - FarFieldExpansion se_cmp; + FarFieldExpansion se; + FarFieldExpansion se_translated; + FarFieldExpansion se_cmp; // initialize expansion objects with respective centers and the bandwidth // squared of 0.1 @@ -309,8 +309,8 @@ int TestTransLocalToLocal(const Matrix &data, const Vector &weights, new_center[0] = new_center[1] = 3.5; // declare expansion objects at (0,0) and other centers - LocalExpansion se; - LocalExpansion se_translated; + LocalExpansion se; + LocalExpansion se_translated; // initialize expansion objects with respective centers and the bandwidth // squared of 0.1 @@ -389,9 +389,9 @@ int TestMixFarField(const Matrix &data, const Vector &weights, data_comb.PrintDebug(); // declare expansion objects at (0,0) and other centers - FarFieldExpansion se; - FarFieldExpansion se2; - FarFieldExpansion se3; + FarFieldExpansion se; + FarFieldExpansion se2; + FarFieldExpansion se3; // initialize expansion objects with respective centers and the bandwidth // squared of 0.5 @@ -479,9 +479,9 @@ int TestConvolveFarField(const Matrix &data, const Vector &weights, data3.PrintDebug(); // declare expansion objects at (0,0) and other centers - FarFieldExpansion se; - FarFieldExpansion se2; - FarFieldExpansion se3; + FarFieldExpansion se; + FarFieldExpansion se2; + FarFieldExpansion se3; // initialize expansion objects with respective centers and the bandwidth // squared of 0.5