From 3c93fbf30f3e8a35a2d6c3191a2e18a25005a85c Mon Sep 17 00:00:00 2001 From: Dongryeol Lee Date: Tue, 16 Oct 2007 21:59:12 +0000 Subject: [PATCH] generalized plug-in series-expansion framework completed --- fastlib/u/dongryel/kde/kde.h | 34 +- fastlib/u/dongryel/kde/main.cc | 49 +- .../u/dongryel/series_expansion/kernel_aux.h | 459 +++++++++++------- fastlib/u/dongryel/series_expansion/main.cc | 2 +- .../mult_farfield_expansion.h | 16 +- .../series_expansion/mult_local_expansion.h | 10 +- 6 files changed, 346 insertions(+), 224 deletions(-) diff --git a/fastlib/u/dongryel/kde/kde.h b/fastlib/u/dongryel/kde/kde.h index 27c7a0f06b..bfd4a03b54 100644 --- a/fastlib/u/dongryel/kde/kde.h +++ b/fastlib/u/dongryel/kde/kde.h @@ -4,6 +4,8 @@ #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/mult_farfield_expansion.h" +#include "u/dongryel/series_expansion/mult_local_expansion.h" #include "u/dongryel/series_expansion/kernel_aux.h" template @@ -109,8 +111,11 @@ template class FastKde { public: - + + // forward declaration of KdeStat class class KdeStat; + + // our tree type using the KdeStat typedef BinarySpaceTree, Matrix, KdeStat > Tree; class KdeStat { @@ -156,12 +161,12 @@ class FastKde { /** * Far field expansion created by the reference points in this node. */ - FarFieldExpansion farfield_expansion_; + typename TKernelAux::TFarFieldExpansion farfield_expansion_; /** * Local expansion stored in this node. */ - LocalExpansion local_expansion_; + typename TKernelAux::TLocalExpansion local_expansion_; /** Initialize the statistics */ void Init() { @@ -175,7 +180,9 @@ class FastKde { mass_t_ = 0; } - void Init(double bandwidth, SeriesExpansionAux *sea) { + void Init(double bandwidth, + typename TKernelAux::TSeriesExpansionAux *sea) { + farfield_expansion_.Init(bandwidth, sea); local_expansion_.Init(bandwidth, sea); } @@ -191,7 +198,7 @@ class FastKde { } void Init(double bandwidth, const Vector& center, - SeriesExpansionAux *sea) { + typename TKernelAux::TSeriesExpansionAux *sea) { farfield_expansion_.Init(bandwidth, center, sea); local_expansion_.Init(bandwidth, center, sea); @@ -211,10 +218,9 @@ class FastKde { right_stat.mass_t_ -= min_mass_t; } - void PushDownTokens(KdeStat &left_stat, KdeStat &right_stat, - double *de, - LocalExpansion - *local_expansion, double *dt) { + void PushDownTokens + (KdeStat &left_stat, KdeStat &right_stat, double *de, + typename TKernelAux::TLocalExpansion *local_expansion, double *dt) { if(de != NULL) { double de_ref = *de; @@ -244,7 +250,7 @@ class FastKde { private: /** series expansion auxililary object */ - SeriesExpansionAux sea_; + typename TKernelAux::TSeriesExpansionAux sea_; /** query dataset */ Matrix qset_; @@ -436,10 +442,10 @@ class FastKde { KdeStat &rstat = rnode->stat(); // expansion objects - FarFieldExpansion &farfield_expansion - = rstat.farfield_expansion_; - LocalExpansion &local_expansion - = qstat.local_expansion_; + typename TKernelAux::TFarFieldExpansion &farfield_expansion = + rstat.farfield_expansion_; + typename TKernelAux::TLocalExpansion &local_expansion = + qstat.local_expansion_; // number of reference points int num_references = rnode->count(); diff --git a/fastlib/u/dongryel/kde/main.cc b/fastlib/u/dongryel/kde/main.cc index 2e4f6d5045..7afc416d1d 100644 --- a/fastlib/u/dongryel/kde/main.cc +++ b/fastlib/u/dongryel/kde/main.cc @@ -10,21 +10,48 @@ int main(int argc, char *argv[]) { const char *kernel_name = fx_param_str_req(NULL, "kernel"); if(!strcmp(kernel_name, "gaussian")) { - FastKde fast_kde; - fast_kde.Init(); - fast_kde.Compute(fx_param_double(NULL, "tau", 0.1)); - - if(fx_param_exists(NULL, "fast_kde_output")) { - fast_kde.PrintDebug(); - } - + Vector fast_kde_results; - fast_kde_results.Alias(fast_kde.get_density_estimates()); + Matrix query_dataset; + Matrix reference_dataset; + + // for O(p^D) expansion + if(fx_param_exists(NULL, "multiplicative_expansion")) { + + printf("O(p^D) expansion KDE\n"); + FastKde fast_kde; + fast_kde.Init(); + fast_kde.Compute(fx_param_double(NULL, "tau", 0.1)); + + if(fx_param_exists(NULL, "fast_kde_output")) { + fast_kde.PrintDebug(); + } + + fast_kde_results.Copy(fast_kde.get_density_estimates()); + query_dataset.Copy(fast_kde.get_query_dataset()); + reference_dataset.Copy(fast_kde.get_reference_dataset()); + } + + // otherwise do O(D^p) expansion + else { + + printf("O(D^p) expansion KDE\n"); + FastKde fast_kde; + fast_kde.Init(); + fast_kde.Compute(fx_param_double(NULL, "tau", 0.1)); + + if(fx_param_exists(NULL, "fast_kde_output")) { + fast_kde.PrintDebug(); + } + + fast_kde_results.Copy(fast_kde.get_density_estimates()); + query_dataset.Copy(fast_kde.get_query_dataset()); + reference_dataset.Copy(fast_kde.get_reference_dataset()); + } if(do_naive) { NaiveKde naive_kde; - naive_kde.Init(fast_kde.get_query_dataset(), - fast_kde.get_reference_dataset()); + naive_kde.Init(query_dataset, reference_dataset); naive_kde.Compute(); if(fx_param_exists(NULL, "naive_kde_output")) { diff --git a/fastlib/u/dongryel/series_expansion/kernel_aux.h b/fastlib/u/dongryel/series_expansion/kernel_aux.h index 85bf374239..005804414a 100644 --- a/fastlib/u/dongryel/series_expansion/kernel_aux.h +++ b/fastlib/u/dongryel/series_expansion/kernel_aux.h @@ -10,8 +10,267 @@ #include "fastlib/fastlib.h" -#include "mult_series_expansion_aux.h" -#include "series_expansion_aux.h" +#include "farfield_expansion.h" +#include "local_expansion.h" +#include "mult_farfield_expansion.h" +#include "mult_local_expansion.h" + +/** + * Auxiliary computer class for multiplicative p^D expansion for Gaussian + * kernel + */ +/** + * Auxiliary computer class for Gaussian kernel + */ +class GaussianKernelMultAux { + FORBID_COPY(GaussianKernelMultAux); + + public: + + typedef MultSeriesExpansionAux TSeriesExpansionAux; + + typedef MultFarFieldExpansion + TFarFieldExpansion; + + typedef MultLocalExpansion + TLocalExpansion; + + /** pointer to the Gaussian kernel */ + GaussianKernel *kernel_; + + /** pointer to the series expansion auxiliary object */ + MultSeriesExpansionAux *sea_; + + GaussianKernelMultAux() {} + + ~GaussianKernelMultAux() {} + + 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, + const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, + double max_dist_sqd_regions, double max_error, + double *actual_error) const { + + double max_far_field_length = 0; + + for(index_t d = 0; d < sea_->get_dimension(); d++) { + DRange far_range = far_field_region.get(d); + max_far_field_length = max(max_far_field_length, far_range.width()); + } + + double two_times_bandwidth = sqrt(kernel_->bandwidth_sq()) * 2; + double r = max_far_field_length / two_times_bandwidth; + + int dim = sea_->get_dimension(); + double r_raised_to_p_alpha = 1.0; + double ret, ret2; + int p_alpha = 0; + double factorialvalue = 1.0; + double first_factor, second_factor; + double one_minus_r; + + // In this case, it is "impossible" to prune for the Gaussian kernel. + if(r >= 1.0) { + return -1; + } + one_minus_r = 1.0 - r; + ret = 1.0 / pow(one_minus_r, dim); + + do { + factorialvalue *= (p_alpha + 1); + + if(factorialvalue < 0.0 || p_alpha > sea_->get_max_order() - 1) { + return -1; + } + + r_raised_to_p_alpha *= r; + first_factor = 1.0 - r_raised_to_p_alpha; + second_factor = r_raised_to_p_alpha / sqrt(factorialvalue); + + ret2 = ret * (pow((first_factor + second_factor), dim) - + pow(first_factor, dim)); + + if(ret2 <= max_error) { + break; + } + + p_alpha++; + + } while(1); + + *actual_error = ret2; + return p_alpha; + } + + int OrderForConvertingFromFarFieldToLocal + (const DHrectBound<2> &far_field_region, + const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, + double max_dist_sqd_regions, double max_error, + double *actual_error) const { + + double max_far_field_length = 0; + double max_local_field_length = 0; + + for(index_t d = 0; d < sea_->get_dimension(); d++) { + DRange far_range = far_field_region.get(d); + DRange local_range = local_field_region.get(d); + max_far_field_length = max(max_far_field_length, far_range.width()); + max_local_field_length = max(max_local_field_length, + local_range.width()); + } + + double two_times_bandwidth = sqrt(kernel_->bandwidth_sq()) * 2; + double r = max_far_field_length / two_times_bandwidth; + double r2 = max_local_field_length / two_times_bandwidth; + + int dim = sea_->get_dimension(); + double r_raised_to_p_alpha = 1.0; + double ret, ret2; + int p_alpha = 0; + double factorialvalue = 1.0; + double first_factor, second_factor; + double one_minus_two_r, two_r; + + // In this case, it is "impossible" to prune for the Gaussian kernel. + if(r >= 0.5 || r2 >= 0.5) + return -1; + + r = max(r, r2); + two_r = 2.0 * r; + one_minus_two_r = 1.0 - two_r; + ret = 1.0 / pow(one_minus_two_r * one_minus_two_r, dim); + + do { + factorialvalue *= (p_alpha + 1); + + if(factorialvalue < 0.0 || p_alpha > sea_->get_max_order() - 1) { + return -1; + } + + r_raised_to_p_alpha *= two_r; + first_factor = 1.0 - r_raised_to_p_alpha; + first_factor *= first_factor; + second_factor = r_raised_to_p_alpha * (2.0 - r_raised_to_p_alpha) + / sqrt(factorialvalue); + + ret2 = ret * (pow((first_factor + second_factor), dim) - + pow(first_factor, dim)); + + if(ret2 <= max_error) { + break; + } + + p_alpha++; + + } while(1); + + *actual_error = ret2; + return p_alpha; + } + + int OrderForEvaluatingLocal + (const DHrectBound<2> &far_field_region, + const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, + double max_dist_sqd_regions, double max_error, + double *actual_error) const { + + double max_local_field_length = 0; + + for(index_t d = 0; d < sea_->get_dimension(); d++) { + DRange local_range = local_field_region.get(d); + max_local_field_length = max(max_local_field_length, + local_range.width()); + } + + double two_times_bandwidth = sqrt(kernel_->bandwidth_sq()) * 2; + double r = max_local_field_length / two_times_bandwidth; + + int dim = sea_->get_dimension(); + double r_raised_to_p_alpha = 1.0; + double ret, ret2; + int p_alpha = 0; + double factorialvalue = 1.0; + double first_factor, second_factor; + double one_minus_r; + + // In this case, it is "impossible" to prune for the Gaussian kernel. + if(r >= 1.0) { + return -1; + } + one_minus_r = 1.0 - r; + ret = 1.0 / pow(one_minus_r, dim); + + do { + factorialvalue *= (p_alpha + 1); + + if(factorialvalue < 0.0 || p_alpha > sea_->get_max_order() - 1) { + return -1; + } + + r_raised_to_p_alpha *= r; + first_factor = 1.0 - r_raised_to_p_alpha; + second_factor = r_raised_to_p_alpha / sqrt(factorialvalue); + + ret2 = ret * (pow((first_factor + second_factor), dim) - + pow(first_factor, dim)); + + if(ret2 <= max_error) { + break; + } + + p_alpha++; + + } while(1); + + *actual_error = ret2; + return p_alpha; + } +}; /** * Auxiliary computer class for Gaussian kernel @@ -21,14 +280,19 @@ class GaussianKernelAux { public: + typedef SeriesExpansionAux TSeriesExpansionAux; + + typedef FarFieldExpansion + TFarFieldExpansion; + + typedef LocalExpansion TLocalExpansion; + /** pointer to the Gaussian kernel */ GaussianKernel *kernel_; /** pointer to the series expansion auxiliary object */ SeriesExpansionAux *sea_; - MultSeriesExpansionAux *msea_; - GaussianKernelAux() {} ~GaussianKernelAux() {} @@ -277,187 +541,6 @@ class GaussianKernelAux { *actual_error = ret; return p_alpha; } - - int OrderForEvaluatingMultFarField - (const DHrectBound<2> &far_field_region, - const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, - double max_dist_sqd_regions, double max_error, - double *actual_error) const { - - double max_far_field_length = 0; - - for(index_t d = 0; d < msea_->get_dimension(); d++) { - DRange far_range = far_field_region.get(d); - max_far_field_length = max(max_far_field_length, far_range.width()); - } - - double two_times_bandwidth = sqrt(kernel_->bandwidth_sq()) * 2; - double r = max_far_field_length / two_times_bandwidth; - - int dim = msea_->get_dimension(); - double r_raised_to_p_alpha = 1.0; - double ret, ret2; - int p_alpha = 0; - double factorialvalue = 1.0; - double first_factor, second_factor; - double one_minus_r; - - // In this case, it is "impossible" to prune for the Gaussian kernel. - if(r >= 1.0) { - return -1; - } - one_minus_r = 1.0 - r; - ret = 1.0 / pow(one_minus_r, dim); - - do { - factorialvalue *= (p_alpha + 1); - - if(factorialvalue < 0.0 || p_alpha > msea_->get_max_order() - 1) { - return -1; - } - - r_raised_to_p_alpha *= r; - first_factor = 1.0 - r_raised_to_p_alpha; - second_factor = r_raised_to_p_alpha / sqrt(factorialvalue); - - ret2 = ret * (pow((first_factor + second_factor), dim) - - pow(first_factor, dim)); - - if(ret2 <= max_error) { - break; - } - - p_alpha++; - - } while(1); - - *actual_error = ret2; - return p_alpha; - } - - int OrderForConvertingFromMultFarFieldToMultLocal - (const DHrectBound<2> &far_field_region, - const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, - double max_dist_sqd_regions, double max_error, - double *actual_error) const { - - double max_far_field_length = 0; - double max_local_field_length = 0; - - for(index_t d = 0; d < msea_->get_dimension(); d++) { - DRange far_range = far_field_region.get(d); - DRange local_range = local_field_region.get(d); - max_far_field_length = max(max_far_field_length, far_range.width()); - max_local_field_length = max(max_local_field_length, - local_range.width()); - } - - double two_times_bandwidth = sqrt(kernel_->bandwidth_sq()) * 2; - double r = max_far_field_length / two_times_bandwidth; - double r2 = max_local_field_length / two_times_bandwidth; - - int dim = msea_->get_dimension(); - double r_raised_to_p_alpha = 1.0; - double ret, ret2; - int p_alpha = 0; - double factorialvalue = 1.0; - double first_factor, second_factor; - double one_minus_two_r, two_r; - - // In this case, it is "impossible" to prune for the Gaussian kernel. - if(r >= 0.5 || r2 >= 0.5) - return -1; - - r = max(r, r2); - two_r = 2.0 * r; - one_minus_two_r = 1.0 - two_r; - ret = 1.0 / pow(one_minus_two_r * one_minus_two_r, dim); - - do { - factorialvalue *= (p_alpha + 1); - - if(factorialvalue < 0.0 || p_alpha > msea_->get_max_order() - 1) { - return -1; - } - - r_raised_to_p_alpha *= two_r; - first_factor = 1.0 - r_raised_to_p_alpha; - first_factor *= first_factor; - second_factor = r_raised_to_p_alpha * (2.0 - r_raised_to_p_alpha) - / sqrt(factorialvalue); - - ret2 = ret * (pow((first_factor + second_factor), dim) - - pow(first_factor, dim)); - - if(ret2 <= max_error) { - break; - } - - p_alpha++; - - } while(1); - - *actual_error = ret2; - return p_alpha; - } - - int OrderForEvaluatingMultLocal - (const DHrectBound<2> &far_field_region, - const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, - double max_dist_sqd_regions, double max_error, - double *actual_error) const { - - double max_local_field_length = 0; - - for(index_t d = 0; d < msea_->get_dimension(); d++) { - DRange local_range = local_field_region.get(d); - max_local_field_length = max(max_local_field_length, - local_range.width()); - } - - double two_times_bandwidth = sqrt(kernel_->bandwidth_sq()) * 2; - double r = max_local_field_length / two_times_bandwidth; - - int dim = msea_->get_dimension(); - double r_raised_to_p_alpha = 1.0; - double ret, ret2; - int p_alpha = 0; - double factorialvalue = 1.0; - double first_factor, second_factor; - double one_minus_r; - - // In this case, it is "impossible" to prune for the Gaussian kernel. - if(r >= 1.0) { - return -1; - } - one_minus_r = 1.0 - r; - ret = 1.0 / pow(one_minus_r, dim); - - do { - factorialvalue *= (p_alpha + 1); - - if(factorialvalue < 0.0 || p_alpha > msea_->get_max_order() - 1) { - return -1; - } - - r_raised_to_p_alpha *= r; - first_factor = 1.0 - r_raised_to_p_alpha; - second_factor = r_raised_to_p_alpha / sqrt(factorialvalue); - - ret2 = ret * (pow((first_factor + second_factor), dim) - - pow(first_factor, dim)); - - if(ret2 <= max_error) { - break; - } - - p_alpha++; - - } while(1); - - *actual_error = ret2; - return p_alpha; - } }; /** @@ -468,6 +551,12 @@ class EpanKernelAux { public: + typedef SeriesExpansionAux TSeriesExpansionAux; + + typedef FarFieldExpansion TFarFieldExpansion; + + typedef LocalExpansion TLocalExpansion; + EpanKernel *kernel_; SeriesExpansionAux *sea_; diff --git a/fastlib/u/dongryel/series_expansion/main.cc b/fastlib/u/dongryel/series_expansion/main.cc index 68b541f521..bee118255e 100644 --- a/fastlib/u/dongryel/series_expansion/main.cc +++ b/fastlib/u/dongryel/series_expansion/main.cc @@ -562,7 +562,7 @@ int TestMultEvaluateFarField(const Matrix &data, const Vector &weights, evaluate_here[0] = evaluate_here[1] = 3; // declare expansion objects at (0,0) and other centers - MultFarFieldExpansion se; + MultFarFieldExpansion se; // initialize expansion objects with respective centers and the bandwidth // squared of 0.5 diff --git a/fastlib/u/dongryel/series_expansion/mult_farfield_expansion.h b/fastlib/u/dongryel/series_expansion/mult_farfield_expansion.h index 83789be09e..b780871739 100644 --- a/fastlib/u/dongryel/series_expansion/mult_farfield_expansion.h +++ b/fastlib/u/dongryel/series_expansion/mult_farfield_expansion.h @@ -371,7 +371,7 @@ template // pass in the pointer to the kernel and the series expansion auxiliary // object ka_.kernel_ = &kernel_; - ka_.msea_ = sea_; + ka_.sea_ = sea_; // initialize coefficient array coeffs_.Init(sea_->get_max_total_num_coeffs()); @@ -392,7 +392,7 @@ template // pass in the pointer to the kernel and the series expansion auxiliary // object ka_.kernel_ = &kernel_; - ka_.msea_ = sea_; + ka_.sea_ = sea_; // initialize coefficient array coeffs_.Init(sea_->get_max_total_num_coeffs()); @@ -405,11 +405,11 @@ template const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, double max_dist_sqd_regions, double max_error, double *actual_error) const { - return ka_.OrderForEvaluatingMultFarField(far_field_region, - local_field_region, - min_dist_sqd_regions, - max_dist_sqd_regions, max_error, - actual_error); + return ka_.OrderForEvaluatingFarField(far_field_region, + local_field_region, + min_dist_sqd_regions, + max_dist_sqd_regions, max_error, + actual_error); } template @@ -421,7 +421,7 @@ template double max_error, double *actual_error) const { - return ka_.OrderForConvertingFromMultFarFieldToMultLocal + return ka_.OrderForConvertingFromFarFieldToLocal (far_field_region, local_field_region, min_dist_sqd_regions, max_dist_sqd_regions, max_error, actual_error); } diff --git a/fastlib/u/dongryel/series_expansion/mult_local_expansion.h b/fastlib/u/dongryel/series_expansion/mult_local_expansion.h index 0fcbb866f0..763e567f24 100644 --- a/fastlib/u/dongryel/series_expansion/mult_local_expansion.h +++ b/fastlib/u/dongryel/series_expansion/mult_local_expansion.h @@ -346,7 +346,7 @@ template // pass in the pointer to the kernel and the series expansion auxiliary // object ka_.kernel_ = &kernel_; - ka_.msea_ = sea_; + ka_.sea_ = sea_; // initialize coefficient array coeffs_.Init(sea_->get_max_total_num_coeffs()); @@ -359,10 +359,10 @@ template const DHrectBound<2> &local_field_region, double min_dist_sqd_regions, double max_dist_sqd_regions, double max_error, double *actual_error) const { - return ka_.OrderForEvaluatingMultLocal(far_field_region, local_field_region, - min_dist_sqd_regions, - max_dist_sqd_regions, max_error, - actual_error); + return ka_.OrderForEvaluatingLocal(far_field_region, local_field_region, + min_dist_sqd_regions, + max_dist_sqd_regions, max_error, + actual_error); } template