apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

generalized plug-in series-expansion framework completed

Browse Source
This commit is contained in:
Dongryeol Lee
2007-10-16 21:59:12 +00:00
parent 16838b6e4d
commit 3c93fbf30f
6 changed files with 346 additions and 224 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/dongryel/kde/kde.h
+20 -14
View File
@@ -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<typename TKernel>
@@ -109,8 +111,11 @@ template<typename TKernel, typename TKernelAux>
class FastKde {
public:
// forward declaration of KdeStat class
class KdeStat;
// our tree type using the KdeStat
typedef BinarySpaceTree<DHrectBound<2>, Matrix, KdeStat > Tree;
class KdeStat {
@@ -156,12 +161,12 @@ class FastKde {
/**
* Far field expansion created by the reference points in this node.
*/
FarFieldExpansion<TKernel, TKernelAux> farfield_expansion_;
typename TKernelAux::TFarFieldExpansion farfield_expansion_;
/**
* Local expansion stored in this node.
*/
LocalExpansion<TKernel, TKernelAux> 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<TKernel, TKernelAux>
*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<TKernel, TKernelAux> &farfield_expansion
= rstat.farfield_expansion_;
LocalExpansion<TKernel, TKernelAux> &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();
fastlib/u/dongryel/kde/main.cc
+38 -11
View File
@@ -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<GaussianKernel, GaussianKernelAux> 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<GaussianKernel, GaussianKernelMultAux> 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<GaussianKernel, GaussianKernelAux> 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<GaussianKernel> 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")) {
fastlib/u/dongryel/series_expansion/kernel_aux.h
+274 -185
View File
@@ -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<GaussianKernel, GaussianKernelMultAux>
TFarFieldExpansion;
typedef MultLocalExpansion<GaussianKernel, GaussianKernelMultAux>
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<int> 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<GaussianKernel, GaussianKernelAux>
TFarFieldExpansion;
typedef LocalExpansion<GaussianKernel, GaussianKernelAux> 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<EpanKernel, EpanKernelAux> TFarFieldExpansion;
typedef LocalExpansion<EpanKernel, EpanKernelAux> TLocalExpansion;
EpanKernel *kernel_;
SeriesExpansionAux *sea_;
fastlib/u/dongryel/series_expansion/main.cc
+1 -1
View File
@@ -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<GaussianKernel, GaussianKernelAux> se;
MultFarFieldExpansion<GaussianKernel, GaussianKernelMultAux> se;
// initialize expansion objects with respective centers and the bandwidth
// squared of 0.5
fastlib/u/dongryel/series_expansion/mult_farfield_expansion.h
+8 -8
View File
@@ -371,7 +371,7 @@ template<typename TKernel, typename TKernelAux>
// 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<typename TKernel, typename TKernelAux>
// 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<typename TKernel, typename TKernelAux>
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<typename TKernel, typename TKernelAux>
@@ -421,7 +421,7 @@ template<typename TKernel, typename TKernelAux>
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);
}
fastlib/u/dongryel/series_expansion/mult_local_expansion.h
+5 -5
View File
@@ -346,7 +346,7 @@ template<typename TKernel, typename TKernelAux>
// 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<typename TKernel, typename TKernelAux>
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<typename TKernel, typename TKernelAux>