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

Combined error bound computation and derivative computations into a single class

Browse Source
This commit is contained in:
Dongryeol Lee
2007-10-15 14:27:59 +00:00
parent 9f46e23044
commit a0e402b4cf
11 changed files with 706 additions and 447 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/dongryel/kde/kde.h
+8 -8
View File
@@ -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<typename TKernel>
class NaiveKde {
@@ -105,7 +105,7 @@ class NaiveKde {
};
template<typename TKernel, typename TKernelDerivative>
template<typename TKernel, typename TKernelAux>
class FastKde {
public:
@@ -156,12 +156,12 @@ class FastKde {
/**
* Far field expansion created by the reference points in this node.
*/
FarFieldExpansion<TKernel, TKernelDerivative> farfield_expansion_;
FarFieldExpansion<TKernel, TKernelAux> farfield_expansion_;
/**
* Local expansion stored in this node.
*/
LocalExpansion<TKernel, TKernelDerivative> local_expansion_;
LocalExpansion<TKernel, TKernelAux> local_expansion_;
/** Initialize the statistics */
void Init() {
@@ -213,7 +213,7 @@ class FastKde {
void PushDownTokens(KdeStat &left_stat, KdeStat &right_stat,
double *de,
LocalExpansion<TKernel, TKernelDerivative>
LocalExpansion<TKernel, TKernelAux>
*local_expansion, double *dt) {
if(de != NULL) {
@@ -434,9 +434,9 @@ class FastKde {
KdeStat &rstat = rnode->stat();
// expansion objects
FarFieldExpansion<TKernel, TKernelDerivative> &farfield_expansion
FarFieldExpansion<TKernel, TKernelAux> &farfield_expansion
= rstat.farfield_expansion_;
LocalExpansion<TKernel, TKernelDerivative> &local_expansion
LocalExpansion<TKernel, TKernelAux> &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);
fastlib/u/dongryel/kde/main.cc
+3 -2
View File
@@ -8,7 +8,7 @@ int main(int argc, char *argv[]) {
bool do_naive = fx_param_exists(NULL, "do_naive");
FastKde<GaussianKernel, GaussianKernelDerivative> fast_kde;
FastKde<GaussianKernel, GaussianKernelAux> 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;
}
fastlib/u/dongryel/multibody/main.cc
+2 -2
View File
@@ -22,7 +22,7 @@ int main(int argc, char *argv[])
if(!strcmp(kernel, "gaussianthreebody")) {
fx_timer_start(NULL, "multibody");
MultitreeMultibody<GaussianThreeBodyKernel,
GaussianKernel, GaussianKernelDerivative> 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<AxilrodTellerKernel,
GaussianKernel, GaussianKernelDerivative> mtmb;
GaussianKernel, GaussianKernelAux> mtmb;
mtmb.Init(bandwidth);
mtmb.Compute(tau);
fx_timer_stop(NULL, "multibody");
fastlib/u/dongryel/multibody/multibody_kernel.h
+1 -1
View File
@@ -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 {
fastlib/u/dongryel/proximity_project/pca_tree.h
+2 -4
View File
@@ -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);
}
fastlib/u/dongryel/series_expansion/build.py
+1 -1
View File
@@ -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
fastlib/u/dongryel/series_expansion/farfield_expansion.h
+80 -177
View File
@@ -10,16 +10,16 @@
#include <values.h>
#include "fastlib/fastlib.h"
#include "kernel_derivative.h"
#include "kernel_aux.h"
#include "series_expansion_aux.h"
template<typename TKernel, typename TKernelDerivative>
template<typename TKernel, typename TKernelAux>
class LocalExpansion;
/**
* Far field expansion class
*/
template<typename TKernel, typename TKernelDerivative>
template<typename TKernel, typename TKernelAux>
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<TKernel, TKernelDerivative> &fe2,
const FarFieldExpansion<TKernel, TKernelDerivative> &fe3,
const FarFieldExpansion<TKernel, TKernelAux> &fe2,
const FarFieldExpansion<TKernel, TKernelAux> &fe3,
int order2, int order3) const;
/**
@@ -125,8 +125,8 @@ class FarFieldExpansion {
* expansions
*/
double ConvolveField
(const FarFieldExpansion<TKernel, TKernelDerivative> &fe2,
const FarFieldExpansion<TKernel, TKernelDerivative> &fe3,
(const FarFieldExpansion<TKernel, TKernelAux> &fe2,
const FarFieldExpansion<TKernel, TKernelAux> &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<TKernel, TKernelDerivative> &se);
(LocalExpansion<TKernel, TKernelAux> &se);
};
template<typename TKernel, typename TKernelDerivative>
void FarFieldExpansion<TKernel, TKernelDerivative>::AccumulateCoeffs
template<typename TKernel, typename TKernelAux>
void FarFieldExpansion<TKernel, TKernelAux>::AccumulateCoeffs
(const Matrix& data, const Vector& weights, int begin, int end,
int order) {
@@ -192,7 +192,7 @@ void FarFieldExpansion<TKernel, TKernelDerivative>::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<TKernel, TKernelDerivative>::AccumulateCoeffs
}
}
template<typename TKernel, typename TKernelDerivative>
void FarFieldExpansion<TKernel, TKernelDerivative>::RefineCoeffs
template<typename TKernel, typename TKernelAux>
void FarFieldExpansion<TKernel, TKernelAux>::RefineCoeffs
(const Matrix& data, const Vector& weights, int begin, int end,
int order) {
@@ -276,7 +276,7 @@ void FarFieldExpansion<TKernel, TKernelDerivative>::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<TKernel, TKernelDerivative>::RefineCoeffs
}
}
template<typename TKernel, typename TKernelDerivative>
double FarFieldExpansion<TKernel, TKernelDerivative>::
template<typename TKernel, typename TKernelAux>
double FarFieldExpansion<TKernel, TKernelAux>::
EvaluateField(Matrix* data, int row_num, Vector* x_q) const {
// dimension
@@ -346,7 +346,7 @@ double FarFieldExpansion<TKernel, TKernelDerivative>::
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<TKernel, TKernelDerivative>::
}
// 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<int> 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<TKernel, TKernelDerivative>::
return multipole_sum;
}
template<typename TKernel, typename TKernelDerivative>
double FarFieldExpansion<TKernel, TKernelDerivative>::MixField
template<typename TKernel, typename TKernelAux>
double FarFieldExpansion<TKernel, TKernelAux>::MixField
(const Matrix &data, int node1_begin, int node1_end,
int node2_begin, int node2_end,
const FarFieldExpansion<TKernel, TKernelDerivative> &fe2,
const FarFieldExpansion<TKernel, TKernelDerivative> &fe3,
const FarFieldExpansion<TKernel, TKernelAux> &fe2,
const FarFieldExpansion<TKernel, TKernelAux> &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<int> *multiindex_mapping = sea_->get_multiindex_mapping();
const ArrayList<int> *lower_mapping_index = sea_->get_lower_mapping_index();
@@ -454,8 +454,8 @@ template<typename TKernel, typename TKernelDerivative>
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<typename TKernel, typename TKernelDerivative>
ArrayList <int> beta_mapping = multiindex_mapping[beta];
ArrayList <int> 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<typename TKernel, typename TKernelDerivative>
ArrayList <int> gamma_mapping = multiindex_mapping[gamma];
ArrayList <int> 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<typename TKernel, typename TKernelDerivative>
return sum;
}
template<typename TKernel, typename TKernelDerivative>
double FarFieldExpansion<TKernel, TKernelDerivative>::ConvolveField
(const FarFieldExpansion<TKernel, TKernelDerivative> &fe2,
const FarFieldExpansion<TKernel, TKernelDerivative> &fe3,
template<typename TKernel, typename TKernelAux>
double FarFieldExpansion<TKernel, TKernelAux>::ConvolveField
(const FarFieldExpansion<TKernel, TKernelAux> &fe2,
const FarFieldExpansion<TKernel, TKernelAux> &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<int> *multiindex_mapping = sea_->get_multiindex_mapping();
const ArrayList<int> *lower_mapping_index = sea_->get_lower_mapping_index();
@@ -642,9 +642,9 @@ double FarFieldExpansion<TKernel, TKernelDerivative>::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<TKernel, TKernelDerivative>::ConvolveField
ArrayList <int> alpha_mapping = multiindex_mapping[alpha];
ArrayList <int> 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<TKernel, TKernelDerivative>::ConvolveField
ArrayList <int> beta_mapping = multiindex_mapping[beta];
ArrayList <int> 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<TKernel, TKernelDerivative>::ConvolveField
ArrayList <int> gamma_mapping = multiindex_mapping[gamma];
ArrayList <int> 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<TKernel, TKernelDerivative>::ConvolveField
return sum;
}
template<typename TKernel, typename TKernelDerivative>
void FarFieldExpansion<TKernel, TKernelDerivative>::Init
template<typename TKernel, typename TKernelAux>
void FarFieldExpansion<TKernel, TKernelAux>::Init
(double bandwidth, const Vector& center, SeriesExpansionAux *sea) {
// copy kernel type, center, and bandwidth squared
@@ -773,13 +773,18 @@ template<typename TKernel, typename TKernelDerivative>
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<typename TKernel, typename TKernelDerivative>
void FarFieldExpansion<TKernel, TKernelDerivative>::Init
template<typename TKernel, typename TKernelAux>
void FarFieldExpansion<TKernel, TKernelAux>::Init
(double bandwidth, SeriesExpansionAux *sea) {
// copy kernel type, center, and bandwidth squared
@@ -789,144 +794,42 @@ template<typename TKernel, typename TKernelDerivative>
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<typename TKernel, typename TKernelDerivative>
int FarFieldExpansion<TKernel, TKernelDerivative>::OrderForEvaluating
template<typename TKernel, typename TKernelAux>
int FarFieldExpansion<TKernel, TKernelAux>::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<typename TKernel, typename TKernelDerivative>
int FarFieldExpansion<TKernel, TKernelDerivative>::
OrderForConvertingtoLocal(const DHrectBound<2> &far_field_region,
template<typename TKernel, typename TKernelAux>
int FarFieldExpansion<TKernel, TKernelAux>::
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<typename TKernel, typename TKernelDerivative>
void FarFieldExpansion<TKernel, TKernelDerivative>::PrintDebug
template<typename TKernel, typename TKernelAux>
void FarFieldExpansion<TKernel, TKernelAux>::PrintDebug
(const char *name, FILE *stream) const {
@@ -975,11 +878,11 @@ template<typename TKernel, typename TKernelDerivative>
fprintf(stream, "\n");
}
template<typename TKernel, typename TKernelDerivative>
void FarFieldExpansion<TKernel, TKernelDerivative>::TranslateFromFarField
template<typename TKernel, typename TKernelAux>
void FarFieldExpansion<TKernel, TKernelAux>::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<typename TKernel, typename TKernelDerivative>
} // end of j-loop
}
template<typename TKernel, typename TKernelDerivative>
void FarFieldExpansion<TKernel, TKernelDerivative>::TranslateToLocal
(LocalExpansion<TKernel, TKernelDerivative> &se) {
template<typename TKernel, typename TKernelAux>
void FarFieldExpansion<TKernel, TKernelAux>::TranslateToLocal
(LocalExpansion<TKernel, TKernelAux> &se) {
Vector pos_arrtmp, neg_arrtmp;
Matrix derivative_map;
@@ -1078,7 +981,7 @@ void FarFieldExpansion<TKernel, TKernelDerivative>::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<TKernel, TKernelDerivative>::TranslateToLocal
}
// compute required partial derivatives
kd_.ComputeDirectionalDerivatives(cent_diff, derivative_map);
ka_.ComputeDirectionalDerivatives(cent_diff, derivative_map);
ArrayList<int> beta_plus_alpha;
beta_plus_alpha.Init(dimension);
@@ -1119,7 +1022,7 @@ void FarFieldExpansion<TKernel, TKernelDerivative>::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;
fastlib/u/dongryel/series_expansion/kernel_aux.h
+545
View File
@@ -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<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, 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<int> 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
fastlib/u/dongryel/series_expansion/kernel_derivative.h
-144
View File
@@ -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<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;
}
};
/**
* 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<int> 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
fastlib/u/dongryel/series_expansion/local_expansion.h
+49 -93
View File
@@ -10,16 +10,16 @@
#include <values.h>
#include "fastlib/fastlib.h"
#include "kernel_derivative.h"
#include "kernel_aux.h"
#include "series_expansion_aux.h"
template<typename TKernel, typename TKernelDerivative>
template<typename TKernel, typename TKernelAux>
class FarFieldExpansion;
/**
* Local expansion class
*/
template<typename TKernel, typename TKernelDerivative>
template<typename TKernel, typename TKernelAux>
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<TKernel, TKernelDerivative> &se);
(const FarFieldExpansion<TKernel, TKernelAux> &se);
/**
* Translate to the given local expansion. The translated coefficients
* are added up to the passed-in local expansion coefficients.
*/
void TranslateToLocal
(LocalExpansion<TKernel, TKernelDerivative> &se);
(LocalExpansion<TKernel, TKernelAux> &se);
};
template<typename TKernel, typename TKernelDerivative>
void LocalExpansion<TKernel, TKernelDerivative>::AccumulateCoeffs
template<typename TKernel, typename TKernelAux>
void LocalExpansion<TKernel, TKernelAux>::AccumulateCoeffs
(const Matrix& data, const Vector& weights,
int begin, int end, int order) {
@@ -164,7 +166,7 @@ void LocalExpansion<TKernel, TKernelDerivative>::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<TKernel, TKernelDerivative>::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<int> 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<TKernel, TKernelDerivative>::AccumulateCoeffs
} // End of looping through each reference point.
}
template<typename TKernel, typename TKernelDerivative>
void LocalExpansion<TKernel, TKernelDerivative>::PrintDebug
template<typename TKernel, typename TKernelAux>
void LocalExpansion<TKernel, TKernelAux>::PrintDebug
(const char *name, FILE *stream) const {
@@ -231,8 +233,8 @@ template<typename TKernel, typename TKernelDerivative>
fprintf(stream, "\n");
}
template<typename TKernel, typename TKernelDerivative>
double LocalExpansion<TKernel, TKernelDerivative>::
template<typename TKernel, typename TKernelAux>
double LocalExpansion<TKernel, TKernelAux>::
EvaluateField(Matrix* data, int row_num, Vector* x_q) const {
index_t k, t, tail;
@@ -247,7 +249,7 @@ template<typename TKernel, typename TKernelDerivative>
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<typename TKernel, typename TKernelDerivative>
return sum;
}
template<typename TKernel, typename TKernelDerivative>
void LocalExpansion<TKernel, TKernelDerivative>::Init
template<typename TKernel, typename TKernelAux>
void LocalExpansion<TKernel, TKernelAux>::Init
(double bandwidth, const Vector& center, SeriesExpansionAux *sea) {
// copy kernel type, center, and bandwidth squared
@@ -304,13 +306,18 @@ template<typename TKernel, typename TKernelDerivative>
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<typename TKernel, typename TKernelDerivative>
void LocalExpansion<TKernel, TKernelDerivative>::Init
template<typename TKernel, typename TKernelAux>
void LocalExpansion<TKernel, TKernelAux>::Init
(double bandwidth, SeriesExpansionAux *sea) {
// copy kernel type, center, and bandwidth squared
@@ -319,79 +326,28 @@ template<typename TKernel, typename TKernelDerivative>
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<typename TKernel, typename TKernelDerivative>
int LocalExpansion<TKernel, TKernelDerivative>::OrderForEvaluating
template<typename TKernel, typename TKernelAux>
int LocalExpansion<TKernel, TKernelAux>::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<typename TKernel, typename TKernelDerivative>
void LocalExpansion<TKernel, TKernelDerivative>::TranslateFromFarField
(const FarFieldExpansion<TKernel, TKernelDerivative> &se) {
template<typename TKernel, typename TKernelAux>
void LocalExpansion<TKernel, TKernelAux>::TranslateFromFarField
(const FarFieldExpansion<TKernel, TKernelAux> &se) {
Vector pos_arrtmp, neg_arrtmp;
Matrix derivative_map;
@@ -402,7 +358,7 @@ template<typename TKernel, typename TKernelDerivative>
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<typename TKernel, typename TKernelDerivative>
}
// compute required partial derivatives
kd_.ComputeDirectionalDerivatives(cent_diff, derivative_map);
ka_.ComputeDirectionalDerivatives(cent_diff, derivative_map);
ArrayList<int> beta_plus_alpha;
beta_plus_alpha.Init(dimension);
@@ -443,7 +399,7 @@ template<typename TKernel, typename TKernelDerivative>
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<typename TKernel, typename TKernelDerivative>
}
}
template<typename TKernel, typename TKernelDerivative>
void LocalExpansion<TKernel, TKernelDerivative>::TranslateToLocal
(LocalExpansion<TKernel, TKernelDerivative> &se) {
template<typename TKernel, typename TKernelAux>
void LocalExpansion<TKernel, TKernelAux>::TranslateToLocal
(LocalExpansion<TKernel, TKernelAux> &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<typename TKernel, typename TKernelDerivative>
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;
fastlib/u/dongryel/series_expansion/main.cc
+15 -15
View File
@@ -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<EpanKernel, EpanKernelDerivative> se;
FarFieldExpansion<EpanKernel, EpanKernelAux> 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<GaussianKernel, GaussianKernelDerivative> se;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> 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<GaussianKernel, GaussianKernelDerivative> se;
LocalExpansion<GaussianKernel, GaussianKernelAux> 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<GaussianKernel, GaussianKernelDerivative> se;
FarFieldExpansion<GaussianKernel, GaussianKernelDerivative> se_translated;
FarFieldExpansion<GaussianKernel, GaussianKernelDerivative> se_cmp;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> se;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> se_translated;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> 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<GaussianKernel, GaussianKernelDerivative> se;
LocalExpansion<GaussianKernel, GaussianKernelDerivative> se_translated;
LocalExpansion<GaussianKernel, GaussianKernelAux> se;
LocalExpansion<GaussianKernel, GaussianKernelAux> 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<GaussianKernel, GaussianKernelDerivative> se;
FarFieldExpansion<GaussianKernel, GaussianKernelDerivative> se2;
FarFieldExpansion<GaussianKernel, GaussianKernelDerivative> se3;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> se;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> se2;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> 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<GaussianKernel, GaussianKernelDerivative> se;
FarFieldExpansion<GaussianKernel, GaussianKernelDerivative> se2;
FarFieldExpansion<GaussianKernel, GaussianKernelDerivative> se3;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> se;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> se2;
FarFieldExpansion<GaussianKernel, GaussianKernelAux> se3;
// initialize expansion objects with respective centers and the bandwidth
// squared of 0.5