Separated out the implementations of series expansion into separate impl files
This commit is contained in:
Dongryeol Lee
@@ -5,9 +5,12 @@ librule(
|
||||
headers = ["farfield_expansion.h",
|
||||
"farfield_expansion_impl.h",
|
||||
"mult_farfield_expansion.h",
|
||||
"mult_farfield_expansion_impl.h",
|
||||
"kernel_aux.h",
|
||||
"local_expansion.h",
|
||||
"local_expansion_impl.h"
|
||||
"mult_local_expansion.h",
|
||||
"mult_local_expansion_impl.h",
|
||||
"mult_series_expansion_aux.h",
|
||||
"series_expansion_aux.h",
|
||||
"mult_series_expansion_aux.h"],
|
||||
|
||||
@@ -135,435 +135,8 @@ class LocalExpansion {
|
||||
|
||||
};
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
|
||||
const Vector& weights,
|
||||
int begin, int end,
|
||||
int order) {
|
||||
|
||||
if(order > order_) {
|
||||
order_ = order;
|
||||
}
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
|
||||
// get inverse factorials (precomputed)
|
||||
Vector neg_inv_multiindex_factorials;
|
||||
neg_inv_multiindex_factorials.Alias
|
||||
(sea_->get_neg_inv_multiindex_factorials());
|
||||
|
||||
// declare deritave mapping
|
||||
Matrix derivative_map;
|
||||
derivative_map.Init(dim, order + 1);
|
||||
|
||||
// some temporary variables
|
||||
Vector arrtmp;
|
||||
arrtmp.Init(total_num_coeffs);
|
||||
Vector x_r_minus_x_Q;
|
||||
x_r_minus_x_Q.Init(dim);
|
||||
|
||||
// sqrt two times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// for each data point,
|
||||
for(index_t r = begin; r < end; r++) {
|
||||
|
||||
// calculate x_r - x_Q
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
x_r_minus_x_Q[d] = (center_[d] - data.get(d, r)) /
|
||||
bandwidth_factor;
|
||||
}
|
||||
|
||||
// precompute necessary partial derivatives based on coordinate difference
|
||||
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++) {
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(j);
|
||||
arrtmp[j] = ka_->ComputePartialDerivative(derivative_map, mapping);
|
||||
}
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
coeffs_[j] += neg_inv_multiindex_factorials[j] * weights[r] *
|
||||
arrtmp[j];
|
||||
}
|
||||
} // End of looping through each reference point.
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::PrintDebug(const char *name,
|
||||
FILE *stream) const {
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
|
||||
fprintf(stream, "Local expansion\n");
|
||||
fprintf(stream, "Center: ");
|
||||
|
||||
for (index_t i = 0; i < center_.length(); i++) {
|
||||
fprintf(stream, "%g ", center_[i]);
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
|
||||
fprintf(stream, "f(");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "x_q%d", d);
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, ",");
|
||||
}
|
||||
fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
|
||||
|
||||
for (index_t i = 0; i < total_num_coeffs; i++) {
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(i);
|
||||
fprintf(stream, "%g", coeffs_[i]);
|
||||
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "(x_q%d - (%g))^%d ", d, center_[d], mapping[d]);
|
||||
}
|
||||
|
||||
if(i < total_num_coeffs - 1) {
|
||||
fprintf(stream, " + ");
|
||||
}
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double LocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
|
||||
int row_num) const {
|
||||
|
||||
// if there are no local expansion here, then return 0
|
||||
if(order_ < 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
index_t k, t, tail;
|
||||
|
||||
// total number of coefficient
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
// number of dimensions
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// evaluated sum to be returned
|
||||
double sum = 0;
|
||||
|
||||
// sqrt two bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// temporary variable
|
||||
Vector x_Q_to_x_q;
|
||||
x_Q_to_x_q.Init(dim);
|
||||
Vector tmp;
|
||||
tmp.Init(total_num_coeffs);
|
||||
ArrayList<int> heads;
|
||||
heads.Init(dim + 1);
|
||||
|
||||
// compute (x_q - x_Q) / (sqrt(2h^2))
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_Q_to_x_q[i] = (data.get(i, row_num) - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < dim; i++)
|
||||
heads[i] = 0;
|
||||
heads[dim] = MAXINT;
|
||||
|
||||
tmp[0] = 1.0;
|
||||
|
||||
for(k = 1, t = 1, tail = 1; k <= order_; k++, tail = t) {
|
||||
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
int head = heads[i];
|
||||
heads[i] = t;
|
||||
|
||||
for(index_t j = head; j < tail; j++, t++) {
|
||||
tmp[t] = tmp[j] * x_Q_to_x_q[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
sum += coeffs_[i] * tmp[i];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double LocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
|
||||
|
||||
// if there are no local expansion here, then return 0
|
||||
if(order_ < 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
index_t k, t, tail;
|
||||
|
||||
// total number of coefficient
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
// number of dimensions
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// evaluated sum to be returned
|
||||
double sum = 0;
|
||||
|
||||
// sqrt two bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// temporary variable
|
||||
Vector x_Q_to_x_q;
|
||||
x_Q_to_x_q.Init(dim);
|
||||
Vector tmp;
|
||||
tmp.Init(total_num_coeffs);
|
||||
ArrayList<int> heads;
|
||||
heads.Init(dim + 1);
|
||||
|
||||
// compute (x_q - x_Q) / (sqrt(2h^2))
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_Q_to_x_q[i] = (x_q[i] - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < dim; i++)
|
||||
heads[i] = 0;
|
||||
heads[dim] = MAXINT;
|
||||
|
||||
tmp[0] = 1.0;
|
||||
|
||||
for(k = 1, t = 1, tail = 1; k <= order_; k++, tail = t) {
|
||||
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
int head = heads[i];
|
||||
heads[i] = t;
|
||||
|
||||
for(index_t j = head; j < tail; j++, t++) {
|
||||
tmp[t] = tmp[j] * x_Q_to_x_q[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
sum += coeffs_[i] * tmp[i];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::Init(const Vector& center,
|
||||
const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
center_.Copy(center);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
ka_ = &ka;
|
||||
|
||||
// initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
center_.Init(sea_->get_dimension());
|
||||
ka_ = &ka;
|
||||
|
||||
// initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int LocalExpansion<TKernelAux>::OrderForEvaluating
|
||||
(const TBound &far_field_region,
|
||||
const TBound &local_field_region, double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions, double max_error, double *actual_error) const {
|
||||
|
||||
return ka_->OrderForEvaluatingLocal(far_field_region, local_field_region,
|
||||
min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error,
|
||||
actual_error);
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::TranslateFromFarField
|
||||
(const FarFieldExpansion<TKernelAux> &se) {
|
||||
|
||||
Vector pos_arrtmp, neg_arrtmp;
|
||||
Matrix derivative_map;
|
||||
Vector far_center;
|
||||
Vector cent_diff;
|
||||
Vector far_coeffs;
|
||||
int dimension = sea_->get_dimension();
|
||||
int far_order = se.get_order();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
|
||||
int limit;
|
||||
double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
|
||||
|
||||
// get center and coefficients for far field expansion
|
||||
far_center.Alias(*(se.get_center()));
|
||||
far_coeffs.Alias(se.get_coeffs());
|
||||
cent_diff.Init(dimension);
|
||||
|
||||
// if the order of the far field expansion is greater than the
|
||||
// local one we are adding onto, then increase the order.
|
||||
if(far_order > order_) {
|
||||
order_ = far_order;
|
||||
}
|
||||
|
||||
// compute Gaussian derivative
|
||||
limit = 2 * order_ + 1;
|
||||
derivative_map.Init(dimension, limit);
|
||||
pos_arrtmp.Init(total_num_coeffs);
|
||||
neg_arrtmp.Init(total_num_coeffs);
|
||||
|
||||
// compute center difference divided by bw_times_sqrt_two;
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
cent_diff[j] = (center_[j] - far_center[j]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute required partial derivatives
|
||||
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
|
||||
ArrayList<int> beta_plus_alpha;
|
||||
beta_plus_alpha.Init(dimension);
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
const ArrayList<int> &beta_mapping = sea_->get_multiindex(j);
|
||||
pos_arrtmp[j] = neg_arrtmp[j] = 0;
|
||||
|
||||
for(index_t k = 0; k < total_num_coeffs; k++) {
|
||||
|
||||
const ArrayList<int> &alpha_mapping = sea_->get_multiindex(k);
|
||||
for(index_t d = 0; d < dimension; d++) {
|
||||
beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
|
||||
}
|
||||
double derivative_factor =
|
||||
ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
|
||||
|
||||
double prod = far_coeffs[k] * derivative_factor;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_arrtmp[j] += prod;
|
||||
}
|
||||
else {
|
||||
neg_arrtmp[j] += prod;
|
||||
}
|
||||
} // end of k-loop
|
||||
} // end of j-loop
|
||||
|
||||
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
coeffs_[j] += (pos_arrtmp[j] + neg_arrtmp[j]) * C_k_neg[j];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::TranslateToLocal(LocalExpansion &se) {
|
||||
|
||||
// if there are no local coefficients to translate, return
|
||||
if(order_ < 0) {
|
||||
return;
|
||||
}
|
||||
|
||||
// get the center and the order and the total number of coefficients of
|
||||
// the expansion we are translating from. Also get coefficients we
|
||||
// are translating
|
||||
Vector new_center;
|
||||
new_center.Alias(*(se.get_center()));
|
||||
int prev_order = se.get_order();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
const ArrayList < int > *upper_mapping_index =
|
||||
sea_->get_upper_mapping_index();
|
||||
Vector new_coeffs;
|
||||
new_coeffs.Alias(se.get_coeffs());
|
||||
|
||||
// dimension
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// temporary variable
|
||||
ArrayList<int> tmp_storage;
|
||||
tmp_storage.Init(dim);
|
||||
|
||||
// sqrt two times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// center difference between the old center and the new one
|
||||
Vector center_diff;
|
||||
center_diff.Init(dim);
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// set to the new order if the order of the expansion we are translating
|
||||
// from is higher
|
||||
if(prev_order < order_) {
|
||||
se.set_order(order_);
|
||||
}
|
||||
|
||||
// inverse multiindex factorials
|
||||
Vector C_k;
|
||||
C_k.Alias(sea_->get_inv_multiindex_factorials());
|
||||
|
||||
// do the actual translation
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
const ArrayList<int> &alpha_mapping = sea_->get_multiindex(j);
|
||||
const ArrayList<int> &upper_mappings_for_alpha = upper_mapping_index[j];
|
||||
double pos_coeffs = 0;
|
||||
double neg_coeffs = 0;
|
||||
|
||||
for(index_t k = 0; k < upper_mappings_for_alpha.size(); k++) {
|
||||
|
||||
if(upper_mappings_for_alpha[k] >= total_num_coeffs) {
|
||||
break;
|
||||
}
|
||||
|
||||
const ArrayList<int> &beta_mapping =
|
||||
sea_->get_multiindex(upper_mappings_for_alpha[k]);
|
||||
int flag = 0;
|
||||
double diff1 = 1.0;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
tmp_storage[l] = beta_mapping[l] - alpha_mapping[l];
|
||||
|
||||
if(tmp_storage[l] < 0) {
|
||||
flag = 1;
|
||||
break;
|
||||
}
|
||||
} // end of looping over dimension
|
||||
|
||||
if(flag)
|
||||
continue;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
diff1 *= pow(center_diff[l], tmp_storage[l]);
|
||||
}
|
||||
|
||||
double prod = coeffs_[upper_mappings_for_alpha[k]] * diff1 *
|
||||
sea_->get_n_multichoose_k_by_pos(upper_mappings_for_alpha[k], j);
|
||||
|
||||
if(prod > 0) {
|
||||
pos_coeffs += prod;
|
||||
}
|
||||
else {
|
||||
neg_coeffs += prod;
|
||||
}
|
||||
|
||||
} // end of k loop
|
||||
|
||||
new_coeffs[j] += pos_coeffs + neg_coeffs;
|
||||
} // end of j loop
|
||||
}
|
||||
#define INSIDE_LOCAL_EXPANSION_H
|
||||
#include "local_expansion_impl.h"
|
||||
#undef INSIDE_LOCAL_EXPANSION_H
|
||||
|
||||
#endif
|
||||
|
||||
@@ -0,0 +1,441 @@
|
||||
#ifndef INSIDE_LOCAL_EXPANSION_H
|
||||
#error "This is not a public header file!"
|
||||
#endif
|
||||
|
||||
#ifndef LOCAL_EXPANSION_IMPL_H
|
||||
#define LOCAL_EXPANSION_IMPL_H
|
||||
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
|
||||
const Vector& weights,
|
||||
int begin, int end,
|
||||
int order) {
|
||||
|
||||
if(order > order_) {
|
||||
order_ = order;
|
||||
}
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
|
||||
// get inverse factorials (precomputed)
|
||||
Vector neg_inv_multiindex_factorials;
|
||||
neg_inv_multiindex_factorials.Alias
|
||||
(sea_->get_neg_inv_multiindex_factorials());
|
||||
|
||||
// declare deritave mapping
|
||||
Matrix derivative_map;
|
||||
derivative_map.Init(dim, order + 1);
|
||||
|
||||
// some temporary variables
|
||||
Vector arrtmp;
|
||||
arrtmp.Init(total_num_coeffs);
|
||||
Vector x_r_minus_x_Q;
|
||||
x_r_minus_x_Q.Init(dim);
|
||||
|
||||
// sqrt two times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// for each data point,
|
||||
for(index_t r = begin; r < end; r++) {
|
||||
|
||||
// calculate x_r - x_Q
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
x_r_minus_x_Q[d] = (center_[d] - data.get(d, r)) /
|
||||
bandwidth_factor;
|
||||
}
|
||||
|
||||
// precompute necessary partial derivatives based on coordinate difference
|
||||
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++) {
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(j);
|
||||
arrtmp[j] = ka_->ComputePartialDerivative(derivative_map, mapping);
|
||||
}
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
coeffs_[j] += neg_inv_multiindex_factorials[j] * weights[r] *
|
||||
arrtmp[j];
|
||||
}
|
||||
} // End of looping through each reference point.
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::PrintDebug(const char *name,
|
||||
FILE *stream) const {
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
|
||||
fprintf(stream, "Local expansion\n");
|
||||
fprintf(stream, "Center: ");
|
||||
|
||||
for (index_t i = 0; i < center_.length(); i++) {
|
||||
fprintf(stream, "%g ", center_[i]);
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
|
||||
fprintf(stream, "f(");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "x_q%d", d);
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, ",");
|
||||
}
|
||||
fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
|
||||
|
||||
for (index_t i = 0; i < total_num_coeffs; i++) {
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(i);
|
||||
fprintf(stream, "%g", coeffs_[i]);
|
||||
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "(x_q%d - (%g))^%d ", d, center_[d], mapping[d]);
|
||||
}
|
||||
|
||||
if(i < total_num_coeffs - 1) {
|
||||
fprintf(stream, " + ");
|
||||
}
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double LocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
|
||||
int row_num) const {
|
||||
|
||||
// if there are no local expansion here, then return 0
|
||||
if(order_ < 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
index_t k, t, tail;
|
||||
|
||||
// total number of coefficient
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
// number of dimensions
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// evaluated sum to be returned
|
||||
double sum = 0;
|
||||
|
||||
// sqrt two bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// temporary variable
|
||||
Vector x_Q_to_x_q;
|
||||
x_Q_to_x_q.Init(dim);
|
||||
Vector tmp;
|
||||
tmp.Init(total_num_coeffs);
|
||||
ArrayList<int> heads;
|
||||
heads.Init(dim + 1);
|
||||
|
||||
// compute (x_q - x_Q) / (sqrt(2h^2))
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_Q_to_x_q[i] = (data.get(i, row_num) - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < dim; i++)
|
||||
heads[i] = 0;
|
||||
heads[dim] = MAXINT;
|
||||
|
||||
tmp[0] = 1.0;
|
||||
|
||||
for(k = 1, t = 1, tail = 1; k <= order_; k++, tail = t) {
|
||||
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
int head = heads[i];
|
||||
heads[i] = t;
|
||||
|
||||
for(index_t j = head; j < tail; j++, t++) {
|
||||
tmp[t] = tmp[j] * x_Q_to_x_q[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
sum += coeffs_[i] * tmp[i];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double LocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
|
||||
|
||||
// if there are no local expansion here, then return 0
|
||||
if(order_ < 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
index_t k, t, tail;
|
||||
|
||||
// total number of coefficient
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
// number of dimensions
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// evaluated sum to be returned
|
||||
double sum = 0;
|
||||
|
||||
// sqrt two bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// temporary variable
|
||||
Vector x_Q_to_x_q;
|
||||
x_Q_to_x_q.Init(dim);
|
||||
Vector tmp;
|
||||
tmp.Init(total_num_coeffs);
|
||||
ArrayList<int> heads;
|
||||
heads.Init(dim + 1);
|
||||
|
||||
// compute (x_q - x_Q) / (sqrt(2h^2))
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_Q_to_x_q[i] = (x_q[i] - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < dim; i++)
|
||||
heads[i] = 0;
|
||||
heads[dim] = MAXINT;
|
||||
|
||||
tmp[0] = 1.0;
|
||||
|
||||
for(k = 1, t = 1, tail = 1; k <= order_; k++, tail = t) {
|
||||
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
int head = heads[i];
|
||||
heads[i] = t;
|
||||
|
||||
for(index_t j = head; j < tail; j++, t++) {
|
||||
tmp[t] = tmp[j] * x_Q_to_x_q[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
sum += coeffs_[i] * tmp[i];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::Init(const Vector& center,
|
||||
const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
center_.Copy(center);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
ka_ = &ka;
|
||||
|
||||
// initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
center_.Init(sea_->get_dimension());
|
||||
ka_ = &ka;
|
||||
|
||||
// initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int LocalExpansion<TKernelAux>::OrderForEvaluating
|
||||
(const TBound &far_field_region,
|
||||
const TBound &local_field_region, double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions, double max_error, double *actual_error) const {
|
||||
|
||||
return ka_->OrderForEvaluatingLocal(far_field_region, local_field_region,
|
||||
min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error,
|
||||
actual_error);
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::TranslateFromFarField
|
||||
(const FarFieldExpansion<TKernelAux> &se) {
|
||||
|
||||
Vector pos_arrtmp, neg_arrtmp;
|
||||
Matrix derivative_map;
|
||||
Vector far_center;
|
||||
Vector cent_diff;
|
||||
Vector far_coeffs;
|
||||
int dimension = sea_->get_dimension();
|
||||
int far_order = se.get_order();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
|
||||
int limit;
|
||||
double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
|
||||
|
||||
// get center and coefficients for far field expansion
|
||||
far_center.Alias(*(se.get_center()));
|
||||
far_coeffs.Alias(se.get_coeffs());
|
||||
cent_diff.Init(dimension);
|
||||
|
||||
// if the order of the far field expansion is greater than the
|
||||
// local one we are adding onto, then increase the order.
|
||||
if(far_order > order_) {
|
||||
order_ = far_order;
|
||||
}
|
||||
|
||||
// compute Gaussian derivative
|
||||
limit = 2 * order_ + 1;
|
||||
derivative_map.Init(dimension, limit);
|
||||
pos_arrtmp.Init(total_num_coeffs);
|
||||
neg_arrtmp.Init(total_num_coeffs);
|
||||
|
||||
// compute center difference divided by bw_times_sqrt_two;
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
cent_diff[j] = (center_[j] - far_center[j]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute required partial derivatives
|
||||
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
|
||||
ArrayList<int> beta_plus_alpha;
|
||||
beta_plus_alpha.Init(dimension);
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
const ArrayList<int> &beta_mapping = sea_->get_multiindex(j);
|
||||
pos_arrtmp[j] = neg_arrtmp[j] = 0;
|
||||
|
||||
for(index_t k = 0; k < total_num_coeffs; k++) {
|
||||
|
||||
const ArrayList<int> &alpha_mapping = sea_->get_multiindex(k);
|
||||
for(index_t d = 0; d < dimension; d++) {
|
||||
beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
|
||||
}
|
||||
double derivative_factor =
|
||||
ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
|
||||
|
||||
double prod = far_coeffs[k] * derivative_factor;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_arrtmp[j] += prod;
|
||||
}
|
||||
else {
|
||||
neg_arrtmp[j] += prod;
|
||||
}
|
||||
} // end of k-loop
|
||||
} // end of j-loop
|
||||
|
||||
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
coeffs_[j] += (pos_arrtmp[j] + neg_arrtmp[j]) * C_k_neg[j];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void LocalExpansion<TKernelAux>::TranslateToLocal(LocalExpansion &se) {
|
||||
|
||||
// if there are no local coefficients to translate, return
|
||||
if(order_ < 0) {
|
||||
return;
|
||||
}
|
||||
|
||||
// get the center and the order and the total number of coefficients of
|
||||
// the expansion we are translating from. Also get coefficients we
|
||||
// are translating
|
||||
Vector new_center;
|
||||
new_center.Alias(*(se.get_center()));
|
||||
int prev_order = se.get_order();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
const ArrayList < int > *upper_mapping_index =
|
||||
sea_->get_upper_mapping_index();
|
||||
Vector new_coeffs;
|
||||
new_coeffs.Alias(se.get_coeffs());
|
||||
|
||||
// dimension
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// temporary variable
|
||||
ArrayList<int> tmp_storage;
|
||||
tmp_storage.Init(dim);
|
||||
|
||||
// sqrt two times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// center difference between the old center and the new one
|
||||
Vector center_diff;
|
||||
center_diff.Init(dim);
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// set to the new order if the order of the expansion we are translating
|
||||
// from is higher
|
||||
if(prev_order < order_) {
|
||||
se.set_order(order_);
|
||||
}
|
||||
|
||||
// inverse multiindex factorials
|
||||
Vector C_k;
|
||||
C_k.Alias(sea_->get_inv_multiindex_factorials());
|
||||
|
||||
// do the actual translation
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
const ArrayList<int> &alpha_mapping = sea_->get_multiindex(j);
|
||||
const ArrayList<int> &upper_mappings_for_alpha = upper_mapping_index[j];
|
||||
double pos_coeffs = 0;
|
||||
double neg_coeffs = 0;
|
||||
|
||||
for(index_t k = 0; k < upper_mappings_for_alpha.size(); k++) {
|
||||
|
||||
if(upper_mappings_for_alpha[k] >= total_num_coeffs) {
|
||||
break;
|
||||
}
|
||||
|
||||
const ArrayList<int> &beta_mapping =
|
||||
sea_->get_multiindex(upper_mappings_for_alpha[k]);
|
||||
int flag = 0;
|
||||
double diff1 = 1.0;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
tmp_storage[l] = beta_mapping[l] - alpha_mapping[l];
|
||||
|
||||
if(tmp_storage[l] < 0) {
|
||||
flag = 1;
|
||||
break;
|
||||
}
|
||||
} // end of looping over dimension
|
||||
|
||||
if(flag)
|
||||
continue;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
diff1 *= pow(center_diff[l], tmp_storage[l]);
|
||||
}
|
||||
|
||||
double prod = coeffs_[upper_mappings_for_alpha[k]] * diff1 *
|
||||
sea_->get_n_multichoose_k_by_pos(upper_mappings_for_alpha[k], j);
|
||||
|
||||
if(prod > 0) {
|
||||
pos_coeffs += prod;
|
||||
}
|
||||
else {
|
||||
neg_coeffs += prod;
|
||||
}
|
||||
|
||||
} // end of k loop
|
||||
|
||||
new_coeffs[j] += pos_coeffs + neg_coeffs;
|
||||
} // end of j loop
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
@@ -207,571 +207,8 @@ class MultFarFieldExpansion {
|
||||
|
||||
};
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
|
||||
const Vector& weights,
|
||||
int begin, int end,
|
||||
int order) {
|
||||
|
||||
int dim = data.n_rows();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
int max_total_num_coeffs = sea_->get_max_total_num_coeffs();
|
||||
Vector x_r, tmp;
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// initialize temporary variables
|
||||
x_r.Init(dim);
|
||||
tmp.Init(max_total_num_coeffs);
|
||||
Vector pos_coeffs;
|
||||
Vector neg_coeffs;
|
||||
pos_coeffs.Init(max_total_num_coeffs);
|
||||
pos_coeffs.SetZero();
|
||||
neg_coeffs.Init(max_total_num_coeffs);
|
||||
neg_coeffs.SetZero();
|
||||
|
||||
// set to new order if greater
|
||||
if(order_ < order) {
|
||||
order_ = order;
|
||||
}
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
// Repeat for each reference point in this reference node.
|
||||
for(index_t r = begin; r < end; r++) {
|
||||
|
||||
// Calculate the coordinate difference between the ref point and the
|
||||
// centroid.
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_r[i] = (data.get(i, r) - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
tmp.SetZero();
|
||||
tmp[0] = 1.0;
|
||||
|
||||
for(index_t i = 1; i < total_num_coeffs; i++) {
|
||||
|
||||
int index = traversal_order[i];
|
||||
const ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
|
||||
|
||||
// from the direct descendant, recursively compute the multipole moments
|
||||
int direct_ancestor_mapping_pos =
|
||||
lower_mappings[lower_mappings.size() - 2];
|
||||
|
||||
int position = 0;
|
||||
const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
|
||||
const ArrayList<int> &direct_ancestor_mapping =
|
||||
sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
if(mapping[i] != direct_ancestor_mapping[i]) {
|
||||
position = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
tmp[index] = tmp[direct_ancestor_mapping_pos] * x_r[position];
|
||||
}
|
||||
|
||||
// Tally up the result in A_k.
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
|
||||
int index = traversal_order[i];
|
||||
double prod = weights[r] * tmp[index];
|
||||
|
||||
if(prod > 0) {
|
||||
pos_coeffs[index] += prod;
|
||||
}
|
||||
else {
|
||||
neg_coeffs[index] += prod;
|
||||
}
|
||||
}
|
||||
|
||||
} // End of looping through each reference point
|
||||
|
||||
for(index_t r = 0; r < total_num_coeffs; r++) {
|
||||
int index = traversal_order[r];
|
||||
coeffs_[index] += (pos_coeffs[index] + neg_coeffs[index]) *
|
||||
sea_->inv_multiindex_factorials_[index];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::
|
||||
ComputeStratifiedLengthSquareDistribution() {
|
||||
|
||||
squared_coeffs_[0] = coeffs_[0] * coeffs_[0];
|
||||
for(index_t i = 1; i < coeffs_.length(); i++) {
|
||||
squared_coeffs_[i] = squared_coeffs_[i - 1] + coeffs_[i] * coeffs_[i];
|
||||
}
|
||||
squared_coeffs_.PrintDebug();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
|
||||
const Vector& weights,
|
||||
int begin, int end,
|
||||
int order) {
|
||||
|
||||
// if we already have the order of approximation, then return.
|
||||
if(order_ >= order) {
|
||||
return;
|
||||
}
|
||||
|
||||
// otherwise, recompute from scratch... this could be improved potentially
|
||||
// but I believe it will not squeeze out more performance (as in O(D^p)
|
||||
// expansions).
|
||||
else {
|
||||
order_ = order;
|
||||
|
||||
coeffs_.SetZero();
|
||||
AccumulateCoeffs(data, weights, begin, end, order);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Matrix& data,
|
||||
int row_num,
|
||||
int order) const {
|
||||
|
||||
// dimension
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// total number of coefficients
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
|
||||
// square root times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// the evaluated sum
|
||||
double pos_multipole_sum = 0;
|
||||
double neg_multipole_sum = 0;
|
||||
double multipole_sum = 0;
|
||||
|
||||
// computed derivative map
|
||||
Matrix derivative_map;
|
||||
derivative_map.Init(dim, order_ + 1);
|
||||
|
||||
// temporary variable
|
||||
Vector arrtmp;
|
||||
arrtmp.Init(total_num_coeffs);
|
||||
|
||||
// (x_q - x_R) scaled by bandwidth
|
||||
Vector x_q_minus_x_R;
|
||||
x_q_minus_x_R.Init(dim);
|
||||
|
||||
// compute (x_q - x_R) / (sqrt(2h^2))
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
x_q_minus_x_R[d] = (data.get(d, row_num) - center_[d]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute deriative maps based on coordinate difference.
|
||||
ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map);
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
// 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++) {
|
||||
|
||||
int index = traversal_order[j];
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(index);
|
||||
double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
|
||||
double prod = coeffs_[index] * arrtmp;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_multipole_sum += prod;
|
||||
}
|
||||
else {
|
||||
neg_multipole_sum += prod;
|
||||
}
|
||||
}
|
||||
|
||||
multipole_sum = pos_multipole_sum + neg_multipole_sum;
|
||||
return multipole_sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Vector& x_q,
|
||||
int order) const {
|
||||
|
||||
// dimension
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// total number of coefficients
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
|
||||
// square root times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// the evaluated sum
|
||||
double pos_multipole_sum = 0;
|
||||
double neg_multipole_sum = 0;
|
||||
double multipole_sum = 0;
|
||||
|
||||
// computed derivative map
|
||||
Matrix derivative_map;
|
||||
derivative_map.Init(dim, order_ + 1);
|
||||
|
||||
// temporary variable
|
||||
Vector arrtmp;
|
||||
arrtmp.Init(total_num_coeffs);
|
||||
|
||||
// (x_q - x_R) scaled by bandwidth
|
||||
Vector x_q_minus_x_R;
|
||||
x_q_minus_x_R.Init(dim);
|
||||
|
||||
// compute (x_q - x_R) / (sqrt(2h^2))
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
x_q_minus_x_R[d] = (x_q[d] - center_[d]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute deriative maps based on coordinate difference.
|
||||
ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map);
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
ArrayList<int> traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
// 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++) {
|
||||
|
||||
int index = traversal_order[j];
|
||||
ArrayList<int> mapping = sea_->get_multiindex(index);
|
||||
double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
|
||||
double prod = coeffs_[index] * arrtmp;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_multipole_sum += prod;
|
||||
}
|
||||
else {
|
||||
neg_multipole_sum += prod;
|
||||
}
|
||||
}
|
||||
|
||||
multipole_sum = pos_multipole_sum + neg_multipole_sum;
|
||||
return multipole_sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultFarFieldExpansion<TKernelAux>::EvaluateFieldByMonteCarlo
|
||||
(const Matrix& data, int row_num, int order, int num_samples) const {
|
||||
|
||||
// I need to implement this...
|
||||
return 0;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::Init(const Vector& center,
|
||||
const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
center_.Copy(center);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
ka_ = &ka;
|
||||
|
||||
// Initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
|
||||
// Initialize the list of squared coefficients for sampling.
|
||||
squared_coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
squared_coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::Init(const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
center_.Init(sea_->get_dimension());
|
||||
center_.SetZero();
|
||||
ka_ = &ka;
|
||||
|
||||
// Initialize coefficient array.
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
|
||||
// Initialize the list of squared coefficients for sampling.
|
||||
squared_coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
squared_coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int MultFarFieldExpansion<TKernelAux>::OrderForEvaluating
|
||||
(const TBound &far_field_region,
|
||||
const TBound &local_field_region, double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions, double max_error, double *actual_error) const {
|
||||
|
||||
return ka_->OrderForEvaluatingFarField(far_field_region,
|
||||
local_field_region,
|
||||
min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error,
|
||||
actual_error);
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int MultFarFieldExpansion<TKernelAux>::OrderForEvaluatingByMonteCarlo
|
||||
(const TBound &far_field_region,
|
||||
const TBound &local_field_region, double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions, double max_error, double *actual_error,
|
||||
int *num_samples) const {
|
||||
|
||||
int order = ka_->OrderForEvaluatingFarField(far_field_region,
|
||||
local_field_region,
|
||||
min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error,
|
||||
actual_error);
|
||||
*num_samples = std::max((int) sqrt(coeffs_.length()), order + 1);
|
||||
|
||||
return order;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int MultFarFieldExpansion<TKernelAux>::
|
||||
OrderForConvertingToLocal(const TBound &far_field_region,
|
||||
const TBound &local_field_region,
|
||||
double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions,
|
||||
double max_error,
|
||||
double *actual_error) const {
|
||||
|
||||
return ka_->OrderForConvertingFromFarFieldToLocal
|
||||
(far_field_region, local_field_region, min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error, actual_error);
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::PrintDebug
|
||||
(const char *name, FILE *stream) const {
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
|
||||
fprintf(stream, "Far field expansion\n");
|
||||
fprintf(stream, "Center: ");
|
||||
|
||||
for (index_t i = 0; i < center_.length(); i++) {
|
||||
fprintf(stream, "%g ", center_[i]);
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
|
||||
fprintf(stream, "f(");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "x_q%d", d);
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, ",");
|
||||
}
|
||||
fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
|
||||
|
||||
for (index_t i = 0; i < total_num_coeffs; i++) {
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(i);
|
||||
fprintf(stream, "%g ", coeffs_[i]);
|
||||
|
||||
fprintf(stream, "(-1)^(");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "%d", mapping[d]);
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, " + ");
|
||||
}
|
||||
fprintf(stream, ") D^((");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "%d", mapping[d]);
|
||||
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, ",");
|
||||
}
|
||||
fprintf(stream, ")) f(x_q - x_R)");
|
||||
if(i < total_num_coeffs - 1) {
|
||||
fprintf(stream, " + ");
|
||||
}
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::TranslateFromFarField
|
||||
(const MultFarFieldExpansion &se) {
|
||||
|
||||
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);
|
||||
Vector prev_coeffs;
|
||||
Vector prev_center;
|
||||
const ArrayList < int > *multiindex_mapping = sea_->get_multiindex_mapping();
|
||||
const ArrayList < int > *lower_mapping_index =
|
||||
sea_->get_lower_mapping_index();
|
||||
|
||||
ArrayList <int> tmp_storage;
|
||||
Vector center_diff;
|
||||
Vector inv_multiindex_factorials;
|
||||
|
||||
center_diff.Init(dim);
|
||||
|
||||
// retrieve coefficients to be translated and helper mappings
|
||||
prev_coeffs.Alias(se.get_coeffs());
|
||||
prev_center.Alias(*(se.get_center()));
|
||||
tmp_storage.Init(sea_->get_dimension());
|
||||
inv_multiindex_factorials.Alias(sea_->get_inv_multiindex_factorials());
|
||||
|
||||
// no coefficients can be translated
|
||||
if(order == -1) {
|
||||
return;
|
||||
}
|
||||
else {
|
||||
order_ = order;
|
||||
}
|
||||
|
||||
// compute center difference
|
||||
for(index_t j = 0; j < dim; j++) {
|
||||
center_diff[j] = prev_center[j] - center_[j];
|
||||
}
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order];
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
int index = traversal_order[j];
|
||||
const ArrayList <int> &gamma_mapping = multiindex_mapping[index];
|
||||
const ArrayList <int> &lower_mappings_for_gamma =
|
||||
lower_mapping_index[index];
|
||||
double pos_coeff = 0;
|
||||
double neg_coeff = 0;
|
||||
|
||||
for(index_t k = 0; k < lower_mappings_for_gamma.size(); k++) {
|
||||
|
||||
const ArrayList <int> &inner_mapping =
|
||||
multiindex_mapping[lower_mappings_for_gamma[k]];
|
||||
|
||||
int flag = 0;
|
||||
double diff1;
|
||||
|
||||
// compute gamma minus alpha
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
tmp_storage[l] = gamma_mapping[l] - inner_mapping[l];
|
||||
|
||||
if(tmp_storage[l] < 0) {
|
||||
flag = 1;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(flag) {
|
||||
continue;
|
||||
}
|
||||
|
||||
diff1 = 1.0;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
diff1 *= pow(center_diff[l] / bandwidth_factor, tmp_storage[l]);
|
||||
}
|
||||
|
||||
double prod = prev_coeffs[lower_mappings_for_gamma[k]] * diff1 *
|
||||
inv_multiindex_factorials
|
||||
[sea_->ComputeMultiindexPosition(tmp_storage)];
|
||||
|
||||
if(prod > 0) {
|
||||
pos_coeff += prod;
|
||||
}
|
||||
else {
|
||||
neg_coeff += prod;
|
||||
}
|
||||
|
||||
} // end of k-loop
|
||||
|
||||
coeffs_[j] += pos_coeff + neg_coeff;
|
||||
|
||||
} // end of j-loop
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::TranslateToLocal
|
||||
(MultLocalExpansion<TKernelAux> &se, int truncation_order) {
|
||||
|
||||
Vector pos_arrtmp, neg_arrtmp;
|
||||
Matrix derivative_map;
|
||||
Vector local_center;
|
||||
Vector cent_diff;
|
||||
Vector local_coeffs;
|
||||
int local_order = se.get_order();
|
||||
int dimension = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(truncation_order);
|
||||
int limit;
|
||||
double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
|
||||
|
||||
// get center and coefficients for local expansion
|
||||
local_center.Alias(*(se.get_center()));
|
||||
local_coeffs.Alias(se.get_coeffs());
|
||||
cent_diff.Init(dimension);
|
||||
|
||||
// if the order of the far field expansion is greater than the
|
||||
// local one we are adding onto, then increase the order.
|
||||
if(local_order < truncation_order) {
|
||||
se.set_order(truncation_order);
|
||||
}
|
||||
|
||||
// compute Gaussian derivative
|
||||
limit = 2 * truncation_order + 1;
|
||||
derivative_map.Init(dimension, limit);
|
||||
pos_arrtmp.Init(sea_->get_max_total_num_coeffs());
|
||||
neg_arrtmp.Init(sea_->get_max_total_num_coeffs());
|
||||
|
||||
// compute center difference divided by bw_times_sqrt_two;
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
cent_diff[j] = (local_center[j] - center_[j]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute required partial derivatives
|
||||
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
|
||||
ArrayList<int> beta_plus_alpha;
|
||||
beta_plus_alpha.Init(dimension);
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order =
|
||||
sea_->traversal_mapping_[truncation_order];
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
int index = traversal_order[j];
|
||||
const ArrayList<int> &beta_mapping = sea_->get_multiindex(index);
|
||||
pos_arrtmp[index] = neg_arrtmp[index] = 0;
|
||||
|
||||
for(index_t k = 0; k < total_num_coeffs; k++) {
|
||||
|
||||
int index_k = traversal_order[k];
|
||||
|
||||
const ArrayList<int> &alpha_mapping = sea_->get_multiindex(index_k);
|
||||
for(index_t d = 0; d < dimension; d++) {
|
||||
beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
|
||||
}
|
||||
double derivative_factor =
|
||||
ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
|
||||
|
||||
double prod = coeffs_[index_k] * derivative_factor;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_arrtmp[index] += prod;
|
||||
}
|
||||
else {
|
||||
neg_arrtmp[index] += prod;
|
||||
}
|
||||
} // end of k-loop
|
||||
} // end of j-loop
|
||||
|
||||
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
int index = traversal_order[j];
|
||||
local_coeffs[index] += (pos_arrtmp[index] + neg_arrtmp[index]) *
|
||||
C_k_neg[index];
|
||||
}
|
||||
}
|
||||
#define INSIDE_MULT_FARFIELD_EXPANSION_H
|
||||
#include "mult_farfield_expansion_impl.h"
|
||||
#undef INSIDE_MULT_FARFIELD_EXPANSION_H
|
||||
|
||||
#endif
|
||||
|
||||
@@ -0,0 +1,575 @@
|
||||
#ifndef INSIDE_MULT_FARFIELD_EXPANSION_H
|
||||
#error "This file is not a public header file!"
|
||||
#endif
|
||||
|
||||
#ifndef MULT_FARFIELD_EXPANSION_IMPL_H
|
||||
#define MULT_FARFIELD_EXPANSION_IMPL_H
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
|
||||
const Vector& weights,
|
||||
int begin, int end,
|
||||
int order) {
|
||||
|
||||
int dim = data.n_rows();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
int max_total_num_coeffs = sea_->get_max_total_num_coeffs();
|
||||
Vector x_r, tmp;
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// initialize temporary variables
|
||||
x_r.Init(dim);
|
||||
tmp.Init(max_total_num_coeffs);
|
||||
Vector pos_coeffs;
|
||||
Vector neg_coeffs;
|
||||
pos_coeffs.Init(max_total_num_coeffs);
|
||||
pos_coeffs.SetZero();
|
||||
neg_coeffs.Init(max_total_num_coeffs);
|
||||
neg_coeffs.SetZero();
|
||||
|
||||
// set to new order if greater
|
||||
if(order_ < order) {
|
||||
order_ = order;
|
||||
}
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
// Repeat for each reference point in this reference node.
|
||||
for(index_t r = begin; r < end; r++) {
|
||||
|
||||
// Calculate the coordinate difference between the ref point and the
|
||||
// centroid.
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_r[i] = (data.get(i, r) - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
tmp.SetZero();
|
||||
tmp[0] = 1.0;
|
||||
|
||||
for(index_t i = 1; i < total_num_coeffs; i++) {
|
||||
|
||||
int index = traversal_order[i];
|
||||
const ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
|
||||
|
||||
// from the direct descendant, recursively compute the multipole moments
|
||||
int direct_ancestor_mapping_pos =
|
||||
lower_mappings[lower_mappings.size() - 2];
|
||||
|
||||
int position = 0;
|
||||
const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
|
||||
const ArrayList<int> &direct_ancestor_mapping =
|
||||
sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
if(mapping[i] != direct_ancestor_mapping[i]) {
|
||||
position = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
tmp[index] = tmp[direct_ancestor_mapping_pos] * x_r[position];
|
||||
}
|
||||
|
||||
// Tally up the result in A_k.
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
|
||||
int index = traversal_order[i];
|
||||
double prod = weights[r] * tmp[index];
|
||||
|
||||
if(prod > 0) {
|
||||
pos_coeffs[index] += prod;
|
||||
}
|
||||
else {
|
||||
neg_coeffs[index] += prod;
|
||||
}
|
||||
}
|
||||
|
||||
} // End of looping through each reference point
|
||||
|
||||
for(index_t r = 0; r < total_num_coeffs; r++) {
|
||||
int index = traversal_order[r];
|
||||
coeffs_[index] += (pos_coeffs[index] + neg_coeffs[index]) *
|
||||
sea_->inv_multiindex_factorials_[index];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::
|
||||
ComputeStratifiedLengthSquareDistribution() {
|
||||
|
||||
squared_coeffs_[0] = coeffs_[0] * coeffs_[0];
|
||||
for(index_t i = 1; i < coeffs_.length(); i++) {
|
||||
squared_coeffs_[i] = squared_coeffs_[i - 1] + coeffs_[i] * coeffs_[i];
|
||||
}
|
||||
squared_coeffs_.PrintDebug();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
|
||||
const Vector& weights,
|
||||
int begin, int end,
|
||||
int order) {
|
||||
|
||||
// if we already have the order of approximation, then return.
|
||||
if(order_ >= order) {
|
||||
return;
|
||||
}
|
||||
|
||||
// otherwise, recompute from scratch... this could be improved potentially
|
||||
// but I believe it will not squeeze out more performance (as in O(D^p)
|
||||
// expansions).
|
||||
else {
|
||||
order_ = order;
|
||||
|
||||
coeffs_.SetZero();
|
||||
AccumulateCoeffs(data, weights, begin, end, order);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Matrix& data,
|
||||
int row_num,
|
||||
int order) const {
|
||||
|
||||
// dimension
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// total number of coefficients
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
|
||||
// square root times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// the evaluated sum
|
||||
double pos_multipole_sum = 0;
|
||||
double neg_multipole_sum = 0;
|
||||
double multipole_sum = 0;
|
||||
|
||||
// computed derivative map
|
||||
Matrix derivative_map;
|
||||
derivative_map.Init(dim, order_ + 1);
|
||||
|
||||
// temporary variable
|
||||
Vector arrtmp;
|
||||
arrtmp.Init(total_num_coeffs);
|
||||
|
||||
// (x_q - x_R) scaled by bandwidth
|
||||
Vector x_q_minus_x_R;
|
||||
x_q_minus_x_R.Init(dim);
|
||||
|
||||
// compute (x_q - x_R) / (sqrt(2h^2))
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
x_q_minus_x_R[d] = (data.get(d, row_num) - center_[d]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute deriative maps based on coordinate difference.
|
||||
ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map);
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
// 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++) {
|
||||
|
||||
int index = traversal_order[j];
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(index);
|
||||
double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
|
||||
double prod = coeffs_[index] * arrtmp;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_multipole_sum += prod;
|
||||
}
|
||||
else {
|
||||
neg_multipole_sum += prod;
|
||||
}
|
||||
}
|
||||
|
||||
multipole_sum = pos_multipole_sum + neg_multipole_sum;
|
||||
return multipole_sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Vector& x_q,
|
||||
int order) const {
|
||||
|
||||
// dimension
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// total number of coefficients
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
|
||||
// square root times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// the evaluated sum
|
||||
double pos_multipole_sum = 0;
|
||||
double neg_multipole_sum = 0;
|
||||
double multipole_sum = 0;
|
||||
|
||||
// computed derivative map
|
||||
Matrix derivative_map;
|
||||
derivative_map.Init(dim, order_ + 1);
|
||||
|
||||
// temporary variable
|
||||
Vector arrtmp;
|
||||
arrtmp.Init(total_num_coeffs);
|
||||
|
||||
// (x_q - x_R) scaled by bandwidth
|
||||
Vector x_q_minus_x_R;
|
||||
x_q_minus_x_R.Init(dim);
|
||||
|
||||
// compute (x_q - x_R) / (sqrt(2h^2))
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
x_q_minus_x_R[d] = (x_q[d] - center_[d]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute deriative maps based on coordinate difference.
|
||||
ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map);
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
ArrayList<int> traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
// 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++) {
|
||||
|
||||
int index = traversal_order[j];
|
||||
ArrayList<int> mapping = sea_->get_multiindex(index);
|
||||
double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
|
||||
double prod = coeffs_[index] * arrtmp;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_multipole_sum += prod;
|
||||
}
|
||||
else {
|
||||
neg_multipole_sum += prod;
|
||||
}
|
||||
}
|
||||
|
||||
multipole_sum = pos_multipole_sum + neg_multipole_sum;
|
||||
return multipole_sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultFarFieldExpansion<TKernelAux>::EvaluateFieldByMonteCarlo
|
||||
(const Matrix& data, int row_num, int order, int num_samples) const {
|
||||
|
||||
// I need to implement this...
|
||||
return 0;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::Init(const Vector& center,
|
||||
const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
center_.Copy(center);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
ka_ = &ka;
|
||||
|
||||
// Initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
|
||||
// Initialize the list of squared coefficients for sampling.
|
||||
squared_coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
squared_coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::Init(const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
center_.Init(sea_->get_dimension());
|
||||
center_.SetZero();
|
||||
ka_ = &ka;
|
||||
|
||||
// Initialize coefficient array.
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
|
||||
// Initialize the list of squared coefficients for sampling.
|
||||
squared_coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
squared_coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int MultFarFieldExpansion<TKernelAux>::OrderForEvaluating
|
||||
(const TBound &far_field_region,
|
||||
const TBound &local_field_region, double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions, double max_error, double *actual_error) const {
|
||||
|
||||
return ka_->OrderForEvaluatingFarField(far_field_region,
|
||||
local_field_region,
|
||||
min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error,
|
||||
actual_error);
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int MultFarFieldExpansion<TKernelAux>::OrderForEvaluatingByMonteCarlo
|
||||
(const TBound &far_field_region,
|
||||
const TBound &local_field_region, double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions, double max_error, double *actual_error,
|
||||
int *num_samples) const {
|
||||
|
||||
int order = ka_->OrderForEvaluatingFarField(far_field_region,
|
||||
local_field_region,
|
||||
min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error,
|
||||
actual_error);
|
||||
*num_samples = std::max((int) sqrt(coeffs_.length()), order + 1);
|
||||
|
||||
return order;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int MultFarFieldExpansion<TKernelAux>::
|
||||
OrderForConvertingToLocal(const TBound &far_field_region,
|
||||
const TBound &local_field_region,
|
||||
double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions,
|
||||
double max_error,
|
||||
double *actual_error) const {
|
||||
|
||||
return ka_->OrderForConvertingFromFarFieldToLocal
|
||||
(far_field_region, local_field_region, min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error, actual_error);
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::PrintDebug
|
||||
(const char *name, FILE *stream) const {
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
|
||||
fprintf(stream, "Far field expansion\n");
|
||||
fprintf(stream, "Center: ");
|
||||
|
||||
for (index_t i = 0; i < center_.length(); i++) {
|
||||
fprintf(stream, "%g ", center_[i]);
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
|
||||
fprintf(stream, "f(");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "x_q%d", d);
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, ",");
|
||||
}
|
||||
fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
|
||||
|
||||
for (index_t i = 0; i < total_num_coeffs; i++) {
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(i);
|
||||
fprintf(stream, "%g ", coeffs_[i]);
|
||||
|
||||
fprintf(stream, "(-1)^(");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "%d", mapping[d]);
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, " + ");
|
||||
}
|
||||
fprintf(stream, ") D^((");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "%d", mapping[d]);
|
||||
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, ",");
|
||||
}
|
||||
fprintf(stream, ")) f(x_q - x_R)");
|
||||
if(i < total_num_coeffs - 1) {
|
||||
fprintf(stream, " + ");
|
||||
}
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::TranslateFromFarField
|
||||
(const MultFarFieldExpansion &se) {
|
||||
|
||||
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);
|
||||
Vector prev_coeffs;
|
||||
Vector prev_center;
|
||||
const ArrayList < int > *multiindex_mapping = sea_->get_multiindex_mapping();
|
||||
const ArrayList < int > *lower_mapping_index =
|
||||
sea_->get_lower_mapping_index();
|
||||
|
||||
ArrayList <int> tmp_storage;
|
||||
Vector center_diff;
|
||||
Vector inv_multiindex_factorials;
|
||||
|
||||
center_diff.Init(dim);
|
||||
|
||||
// retrieve coefficients to be translated and helper mappings
|
||||
prev_coeffs.Alias(se.get_coeffs());
|
||||
prev_center.Alias(*(se.get_center()));
|
||||
tmp_storage.Init(sea_->get_dimension());
|
||||
inv_multiindex_factorials.Alias(sea_->get_inv_multiindex_factorials());
|
||||
|
||||
// no coefficients can be translated
|
||||
if(order == -1) {
|
||||
return;
|
||||
}
|
||||
else {
|
||||
order_ = order;
|
||||
}
|
||||
|
||||
// compute center difference
|
||||
for(index_t j = 0; j < dim; j++) {
|
||||
center_diff[j] = prev_center[j] - center_[j];
|
||||
}
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order];
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
int index = traversal_order[j];
|
||||
const ArrayList <int> &gamma_mapping = multiindex_mapping[index];
|
||||
const ArrayList <int> &lower_mappings_for_gamma =
|
||||
lower_mapping_index[index];
|
||||
double pos_coeff = 0;
|
||||
double neg_coeff = 0;
|
||||
|
||||
for(index_t k = 0; k < lower_mappings_for_gamma.size(); k++) {
|
||||
|
||||
const ArrayList <int> &inner_mapping =
|
||||
multiindex_mapping[lower_mappings_for_gamma[k]];
|
||||
|
||||
int flag = 0;
|
||||
double diff1;
|
||||
|
||||
// compute gamma minus alpha
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
tmp_storage[l] = gamma_mapping[l] - inner_mapping[l];
|
||||
|
||||
if(tmp_storage[l] < 0) {
|
||||
flag = 1;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(flag) {
|
||||
continue;
|
||||
}
|
||||
|
||||
diff1 = 1.0;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
diff1 *= pow(center_diff[l] / bandwidth_factor, tmp_storage[l]);
|
||||
}
|
||||
|
||||
double prod = prev_coeffs[lower_mappings_for_gamma[k]] * diff1 *
|
||||
inv_multiindex_factorials
|
||||
[sea_->ComputeMultiindexPosition(tmp_storage)];
|
||||
|
||||
if(prod > 0) {
|
||||
pos_coeff += prod;
|
||||
}
|
||||
else {
|
||||
neg_coeff += prod;
|
||||
}
|
||||
|
||||
} // end of k-loop
|
||||
|
||||
coeffs_[j] += pos_coeff + neg_coeff;
|
||||
|
||||
} // end of j-loop
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultFarFieldExpansion<TKernelAux>::TranslateToLocal
|
||||
(MultLocalExpansion<TKernelAux> &se, int truncation_order) {
|
||||
|
||||
Vector pos_arrtmp, neg_arrtmp;
|
||||
Matrix derivative_map;
|
||||
Vector local_center;
|
||||
Vector cent_diff;
|
||||
Vector local_coeffs;
|
||||
int local_order = se.get_order();
|
||||
int dimension = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(truncation_order);
|
||||
int limit;
|
||||
double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
|
||||
|
||||
// get center and coefficients for local expansion
|
||||
local_center.Alias(*(se.get_center()));
|
||||
local_coeffs.Alias(se.get_coeffs());
|
||||
cent_diff.Init(dimension);
|
||||
|
||||
// if the order of the far field expansion is greater than the
|
||||
// local one we are adding onto, then increase the order.
|
||||
if(local_order < truncation_order) {
|
||||
se.set_order(truncation_order);
|
||||
}
|
||||
|
||||
// compute Gaussian derivative
|
||||
limit = 2 * truncation_order + 1;
|
||||
derivative_map.Init(dimension, limit);
|
||||
pos_arrtmp.Init(sea_->get_max_total_num_coeffs());
|
||||
neg_arrtmp.Init(sea_->get_max_total_num_coeffs());
|
||||
|
||||
// compute center difference divided by bw_times_sqrt_two;
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
cent_diff[j] = (local_center[j] - center_[j]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute required partial derivatives
|
||||
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
|
||||
ArrayList<int> beta_plus_alpha;
|
||||
beta_plus_alpha.Init(dimension);
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order =
|
||||
sea_->traversal_mapping_[truncation_order];
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
int index = traversal_order[j];
|
||||
const ArrayList<int> &beta_mapping = sea_->get_multiindex(index);
|
||||
pos_arrtmp[index] = neg_arrtmp[index] = 0;
|
||||
|
||||
for(index_t k = 0; k < total_num_coeffs; k++) {
|
||||
|
||||
int index_k = traversal_order[k];
|
||||
|
||||
const ArrayList<int> &alpha_mapping = sea_->get_multiindex(index_k);
|
||||
for(index_t d = 0; d < dimension; d++) {
|
||||
beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
|
||||
}
|
||||
double derivative_factor =
|
||||
ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
|
||||
|
||||
double prod = coeffs_[index_k] * derivative_factor;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_arrtmp[index] += prod;
|
||||
}
|
||||
else {
|
||||
neg_arrtmp[index] += prod;
|
||||
}
|
||||
} // end of k-loop
|
||||
} // end of j-loop
|
||||
|
||||
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
int index = traversal_order[j];
|
||||
local_coeffs[index] += (pos_arrtmp[index] + neg_arrtmp[index]) *
|
||||
C_k_neg[index];
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
@@ -135,471 +135,8 @@ class MultLocalExpansion {
|
||||
|
||||
};
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
|
||||
const Vector& weights,
|
||||
int begin, int end,
|
||||
int order) {
|
||||
|
||||
if(order > order_) {
|
||||
order_ = order;
|
||||
}
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
|
||||
// get inverse factorials (precomputed)
|
||||
Vector neg_inv_multiindex_factorials;
|
||||
neg_inv_multiindex_factorials.Alias
|
||||
(sea_->get_neg_inv_multiindex_factorials());
|
||||
|
||||
// declare deritave mapping
|
||||
Matrix derivative_map;
|
||||
derivative_map.Init(dim, order + 1);
|
||||
|
||||
// some temporary variables
|
||||
Vector x_r_minus_x_Q;
|
||||
x_r_minus_x_Q.Init(dim);
|
||||
|
||||
// sqrt two times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order];
|
||||
|
||||
// for each data point,
|
||||
for(index_t r = begin; r < end; r++) {
|
||||
|
||||
// calculate x_r - x_Q
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
x_r_minus_x_Q[d] = (center_[d] - data.get(d, r)) /
|
||||
bandwidth_factor;
|
||||
}
|
||||
|
||||
// precompute necessary partial derivatives based on coordinate difference
|
||||
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++) {
|
||||
int index = traversal_order[j];
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(index);
|
||||
double partial_derivative =
|
||||
ka_->ComputePartialDerivative(derivative_map, mapping);
|
||||
coeffs_[index] += neg_inv_multiindex_factorials[index] * weights[r] *
|
||||
partial_derivative;
|
||||
}
|
||||
} // End of looping through each reference point.
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::PrintDebug(const char *name,
|
||||
FILE *stream) const {
|
||||
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
|
||||
fprintf(stream, "Local expansion\n");
|
||||
fprintf(stream, "Center: ");
|
||||
|
||||
for (index_t i = 0; i < center_.length(); i++) {
|
||||
fprintf(stream, "%g ", center_[i]);
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
|
||||
fprintf(stream, "f(");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "x_q%d", d);
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, ",");
|
||||
}
|
||||
fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
|
||||
|
||||
for (index_t i = 0; i < total_num_coeffs; i++) {
|
||||
ArrayList<int> mapping = sea_->get_multiindex(i);
|
||||
fprintf(stream, "%g", coeffs_[i]);
|
||||
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "(x_q%d - (%g))^%d ", d, center_[d], mapping[d]);
|
||||
}
|
||||
|
||||
if(i < total_num_coeffs - 1) {
|
||||
fprintf(stream, " + ");
|
||||
}
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultLocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
|
||||
int row_num) const {
|
||||
|
||||
// if there are no local coefficients, then return 0
|
||||
if(order_ < 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
// total number of coefficient
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
// number of dimensions
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// evaluated sum to be returned
|
||||
double sum = 0;
|
||||
|
||||
// sqrt two bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// temporary variable
|
||||
Vector x_Q_to_x_q;
|
||||
x_Q_to_x_q.Init(dim);
|
||||
Vector tmp;
|
||||
tmp.Init(sea_->get_max_total_num_coeffs());
|
||||
ArrayList<int> heads;
|
||||
heads.Init(dim + 1);
|
||||
|
||||
// compute (x_q - x_Q) / (sqrt(2h^2))
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_Q_to_x_q[i] = (data.get(i, row_num) - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < dim; i++)
|
||||
heads[i] = 0;
|
||||
heads[dim] = MAXINT;
|
||||
|
||||
tmp[0] = 1.0;
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
for(index_t i = 1; i < total_num_coeffs; i++) {
|
||||
|
||||
int index = traversal_order[i];
|
||||
const ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
|
||||
|
||||
// from the direct descendant, recursively compute the multipole moments
|
||||
int direct_ancestor_mapping_pos =
|
||||
lower_mappings[lower_mappings.size() - 2];
|
||||
int position = 0;
|
||||
const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
|
||||
const ArrayList<int> &direct_ancestor_mapping =
|
||||
sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
if(mapping[i] != direct_ancestor_mapping[i]) {
|
||||
position = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
tmp[index] = tmp[direct_ancestor_mapping_pos] * x_Q_to_x_q[position];
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
int index = traversal_order[i];
|
||||
sum += coeffs_[index] * tmp[index];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultLocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
|
||||
|
||||
// if there are no local coefficients, then return 0
|
||||
if(order_ < 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
// total number of coefficient
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
// number of dimensions
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// evaluated sum to be returned
|
||||
double sum = 0;
|
||||
|
||||
// sqrt two bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_.bandwidth_sq());
|
||||
|
||||
// temporary variable
|
||||
Vector x_Q_to_x_q;
|
||||
x_Q_to_x_q.Init(dim);
|
||||
Vector tmp;
|
||||
tmp.Init(sea_->get_max_total_num_coeffs());
|
||||
ArrayList<int> heads;
|
||||
heads.Init(dim + 1);
|
||||
|
||||
// compute (x_q - x_Q) / (sqrt(2h^2))
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_Q_to_x_q[i] = (x_q[i] - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < dim; i++)
|
||||
heads[i] = 0;
|
||||
heads[dim] = MAXINT;
|
||||
|
||||
tmp[0] = 1.0;
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
for(index_t i = 1; i < total_num_coeffs; i++) {
|
||||
|
||||
int index = traversal_order[i];
|
||||
ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
|
||||
|
||||
// from the direct descendant, recursively compute the multipole moments
|
||||
int direct_ancestor_mapping_pos =
|
||||
lower_mappings[lower_mappings.size() - 2];
|
||||
int position = 0;
|
||||
const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
|
||||
const ArrayList<int> &direct_ancestor_mapping =
|
||||
sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
if(mapping[i] != direct_ancestor_mapping[i]) {
|
||||
position = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
tmp[index] = tmp[direct_ancestor_mapping_pos] * x_Q_to_x_q[position];
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
int index = traversal_order[i];
|
||||
sum += coeffs_[index] * tmp[index];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::Init(const Vector& center,
|
||||
const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
center_.Copy(center);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
ka_ = &ka;
|
||||
|
||||
// initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
sea_ = &(ka.sea_);
|
||||
center_.Init(sea_->get_dimension());
|
||||
order_ = -1;
|
||||
ka_ = &ka;
|
||||
|
||||
// initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int MultLocalExpansion<TKernelAux>::OrderForEvaluating
|
||||
(const TBound &far_field_region,
|
||||
const TBound &local_field_region, double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions, double max_error, double *actual_error) const {
|
||||
|
||||
return ka_->OrderForEvaluatingLocal(far_field_region, local_field_region,
|
||||
min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error,
|
||||
actual_error);
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::TranslateFromFarField
|
||||
(const MultFarFieldExpansion<TKernelAux> &se) {
|
||||
|
||||
Vector pos_arrtmp, neg_arrtmp;
|
||||
Matrix derivative_map;
|
||||
Vector far_center;
|
||||
Vector cent_diff;
|
||||
Vector far_coeffs;
|
||||
int dimension = sea_->get_dimension();
|
||||
int far_order = se.get_order();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
|
||||
int limit;
|
||||
double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
|
||||
|
||||
// get center and coefficients for far field expansion
|
||||
far_center.Alias(*(se.get_center()));
|
||||
far_coeffs.Alias(se.get_coeffs());
|
||||
cent_diff.Init(dimension);
|
||||
|
||||
// if the order of the far field expansion is greater than the
|
||||
// local one we are adding onto, then increase the order.
|
||||
if(far_order > order_) {
|
||||
order_ = far_order;
|
||||
}
|
||||
|
||||
// compute Gaussian derivative
|
||||
limit = 2 * order_ + 1;
|
||||
derivative_map.Init(dimension, limit);
|
||||
pos_arrtmp.Init(total_num_coeffs);
|
||||
neg_arrtmp.Init(total_num_coeffs);
|
||||
|
||||
// compute center difference divided by bw_times_sqrt_two;
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
cent_diff[j] = (center_[j] - far_center[j]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute required partial derivatives
|
||||
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
|
||||
ArrayList<int> beta_plus_alpha;
|
||||
beta_plus_alpha.Init(dimension);
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
ArrayList<int> &traversal_order = sea_->traversal_mapping_[far_order];
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
int index_j = traversal_order[j];
|
||||
ArrayList<int> beta_mapping = sea_->get_multiindex(index_j);
|
||||
pos_arrtmp[index_j] = neg_arrtmp[index_j] = 0;
|
||||
|
||||
for(index_t k = 0; k < total_num_coeffs; k++) {
|
||||
|
||||
int index_k = traversal_order[k];
|
||||
ArrayList<int> alpha_mapping = sea_->get_multiindex(index_k);
|
||||
for(index_t d = 0; d < dimension; d++) {
|
||||
beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
|
||||
}
|
||||
double derivative_factor =
|
||||
ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
|
||||
|
||||
double prod = far_coeffs[index_k] * derivative_factor;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_arrtmp[index_j] += prod;
|
||||
}
|
||||
else {
|
||||
neg_arrtmp[index_j] += prod;
|
||||
}
|
||||
} // end of k-loop
|
||||
} // end of j-loop
|
||||
|
||||
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
int index_j = traversal_order[j];
|
||||
coeffs_[index_j] += (pos_arrtmp[index_j] + neg_arrtmp[index_j]) *
|
||||
C_k_neg[index_j];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::TranslateToLocal(MultLocalExpansion &se) {
|
||||
|
||||
// if no local coefficients have formed, then nothing to translate
|
||||
if(order_ < 0) {
|
||||
return;
|
||||
}
|
||||
|
||||
// get the center and the order and the total number of coefficients of
|
||||
// the expansion we are translating from. Also get coefficients we
|
||||
// are translating
|
||||
Vector new_center;
|
||||
new_center.Alias(*(se.get_center()));
|
||||
int prev_order = se.get_order();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
const ArrayList < int > *upper_mapping_index =
|
||||
sea_->get_upper_mapping_index();
|
||||
Vector new_coeffs;
|
||||
new_coeffs.Alias(se.get_coeffs());
|
||||
|
||||
// dimension
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// temporary variable
|
||||
ArrayList<int> tmp_storage;
|
||||
tmp_storage.Init(dim);
|
||||
|
||||
// sqrt two times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// center difference between the old center and the new one
|
||||
Vector center_diff;
|
||||
center_diff.Init(dim);
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// set to the new order if the order of the expansion we are translating
|
||||
// from is higher
|
||||
if(prev_order < order_) {
|
||||
se.set_order(order_);
|
||||
}
|
||||
|
||||
// inverse multiindex factorials
|
||||
Vector C_k;
|
||||
C_k.Alias(sea_->get_inv_multiindex_factorials());
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
// do the actual translation
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
int index_j = traversal_order[j];
|
||||
const ArrayList<int> &alpha_mapping = sea_->get_multiindex(index_j);
|
||||
const ArrayList <int> &upper_mappings_for_alpha =
|
||||
upper_mapping_index[index_j];
|
||||
double pos_coeffs = 0;
|
||||
double neg_coeffs = 0;
|
||||
|
||||
for(index_t k = 0; k < upper_mappings_for_alpha.size(); k++) {
|
||||
|
||||
if(upper_mappings_for_alpha[k] >= total_num_coeffs) {
|
||||
break;
|
||||
}
|
||||
|
||||
const ArrayList<int> &beta_mapping =
|
||||
sea_->get_multiindex(upper_mappings_for_alpha[k]);
|
||||
int flag = 0;
|
||||
double diff1 = 1.0;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
tmp_storage[l] = beta_mapping[l] - alpha_mapping[l];
|
||||
|
||||
if(tmp_storage[l] < 0) {
|
||||
flag = 1;
|
||||
break;
|
||||
}
|
||||
} // end of looping over dimension
|
||||
|
||||
if(flag)
|
||||
continue;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
diff1 *= pow(center_diff[l], tmp_storage[l]);
|
||||
}
|
||||
|
||||
double prod = coeffs_[upper_mappings_for_alpha[k]] * diff1 *
|
||||
sea_->get_n_multichoose_k_by_pos
|
||||
(upper_mappings_for_alpha[k], index_j);
|
||||
|
||||
if(prod > 0) {
|
||||
pos_coeffs += prod;
|
||||
}
|
||||
else {
|
||||
neg_coeffs += prod;
|
||||
}
|
||||
|
||||
} // end of k loop
|
||||
|
||||
new_coeffs[index_j] += pos_coeffs + neg_coeffs;
|
||||
} // end of j loop
|
||||
}
|
||||
#define INSIDE_MULT_LOCAL_EXPANSION_H
|
||||
#include "mult_local_expansion_impl.h"
|
||||
#undef INSIDE_MULT_LOCAL_EXPANSION_H
|
||||
|
||||
#endif
|
||||
|
||||
@@ -0,0 +1,475 @@
|
||||
#ifndef INSIDE_MULT_LOCAL_EXPANSION_H
|
||||
#error "This is not a public header file!"
|
||||
#endif
|
||||
|
||||
#ifndef MULT_LOCAL_EXPANSION_IMPL_H
|
||||
#define MULT_LOCAL_EXPANSION_IMPL_H
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
|
||||
const Vector& weights,
|
||||
int begin, int end,
|
||||
int order) {
|
||||
|
||||
if(order > order_) {
|
||||
order_ = order;
|
||||
}
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order);
|
||||
|
||||
// get inverse factorials (precomputed)
|
||||
Vector neg_inv_multiindex_factorials;
|
||||
neg_inv_multiindex_factorials.Alias
|
||||
(sea_->get_neg_inv_multiindex_factorials());
|
||||
|
||||
// declare deritave mapping
|
||||
Matrix derivative_map;
|
||||
derivative_map.Init(dim, order + 1);
|
||||
|
||||
// some temporary variables
|
||||
Vector x_r_minus_x_Q;
|
||||
x_r_minus_x_Q.Init(dim);
|
||||
|
||||
// sqrt two times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order];
|
||||
|
||||
// for each data point,
|
||||
for(index_t r = begin; r < end; r++) {
|
||||
|
||||
// calculate x_r - x_Q
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
x_r_minus_x_Q[d] = (center_[d] - data.get(d, r)) /
|
||||
bandwidth_factor;
|
||||
}
|
||||
|
||||
// precompute necessary partial derivatives based on coordinate difference
|
||||
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++) {
|
||||
int index = traversal_order[j];
|
||||
const ArrayList<int> &mapping = sea_->get_multiindex(index);
|
||||
double partial_derivative =
|
||||
ka_->ComputePartialDerivative(derivative_map, mapping);
|
||||
coeffs_[index] += neg_inv_multiindex_factorials[index] * weights[r] *
|
||||
partial_derivative;
|
||||
}
|
||||
} // End of looping through each reference point.
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::PrintDebug(const char *name,
|
||||
FILE *stream) const {
|
||||
|
||||
|
||||
int dim = sea_->get_dimension();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
fprintf(stream, "----- SERIESEXPANSION %s ------\n", name);
|
||||
fprintf(stream, "Local expansion\n");
|
||||
fprintf(stream, "Center: ");
|
||||
|
||||
for (index_t i = 0; i < center_.length(); i++) {
|
||||
fprintf(stream, "%g ", center_[i]);
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
|
||||
fprintf(stream, "f(");
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "x_q%d", d);
|
||||
if(d < dim - 1)
|
||||
fprintf(stream, ",");
|
||||
}
|
||||
fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
|
||||
|
||||
for (index_t i = 0; i < total_num_coeffs; i++) {
|
||||
ArrayList<int> mapping = sea_->get_multiindex(i);
|
||||
fprintf(stream, "%g", coeffs_[i]);
|
||||
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
fprintf(stream, "(x_q%d - (%g))^%d ", d, center_[d], mapping[d]);
|
||||
}
|
||||
|
||||
if(i < total_num_coeffs - 1) {
|
||||
fprintf(stream, " + ");
|
||||
}
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultLocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
|
||||
int row_num) const {
|
||||
|
||||
// if there are no local coefficients, then return 0
|
||||
if(order_ < 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
// total number of coefficient
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
// number of dimensions
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// evaluated sum to be returned
|
||||
double sum = 0;
|
||||
|
||||
// sqrt two bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// temporary variable
|
||||
Vector x_Q_to_x_q;
|
||||
x_Q_to_x_q.Init(dim);
|
||||
Vector tmp;
|
||||
tmp.Init(sea_->get_max_total_num_coeffs());
|
||||
ArrayList<int> heads;
|
||||
heads.Init(dim + 1);
|
||||
|
||||
// compute (x_q - x_Q) / (sqrt(2h^2))
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_Q_to_x_q[i] = (data.get(i, row_num) - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < dim; i++)
|
||||
heads[i] = 0;
|
||||
heads[dim] = MAXINT;
|
||||
|
||||
tmp[0] = 1.0;
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
for(index_t i = 1; i < total_num_coeffs; i++) {
|
||||
|
||||
int index = traversal_order[i];
|
||||
const ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
|
||||
|
||||
// from the direct descendant, recursively compute the multipole moments
|
||||
int direct_ancestor_mapping_pos =
|
||||
lower_mappings[lower_mappings.size() - 2];
|
||||
int position = 0;
|
||||
const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
|
||||
const ArrayList<int> &direct_ancestor_mapping =
|
||||
sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
if(mapping[i] != direct_ancestor_mapping[i]) {
|
||||
position = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
tmp[index] = tmp[direct_ancestor_mapping_pos] * x_Q_to_x_q[position];
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
int index = traversal_order[i];
|
||||
sum += coeffs_[index] * tmp[index];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
double MultLocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
|
||||
|
||||
// if there are no local coefficients, then return 0
|
||||
if(order_ < 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
// total number of coefficient
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
|
||||
// number of dimensions
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// evaluated sum to be returned
|
||||
double sum = 0;
|
||||
|
||||
// sqrt two bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_.bandwidth_sq());
|
||||
|
||||
// temporary variable
|
||||
Vector x_Q_to_x_q;
|
||||
x_Q_to_x_q.Init(dim);
|
||||
Vector tmp;
|
||||
tmp.Init(sea_->get_max_total_num_coeffs());
|
||||
ArrayList<int> heads;
|
||||
heads.Init(dim + 1);
|
||||
|
||||
// compute (x_q - x_Q) / (sqrt(2h^2))
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
x_Q_to_x_q[i] = (x_q[i] - center_[i]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < dim; i++)
|
||||
heads[i] = 0;
|
||||
heads[dim] = MAXINT;
|
||||
|
||||
tmp[0] = 1.0;
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
for(index_t i = 1; i < total_num_coeffs; i++) {
|
||||
|
||||
int index = traversal_order[i];
|
||||
ArrayList<int> &lower_mappings = sea_->lower_mapping_index_[index];
|
||||
|
||||
// from the direct descendant, recursively compute the multipole moments
|
||||
int direct_ancestor_mapping_pos =
|
||||
lower_mappings[lower_mappings.size() - 2];
|
||||
int position = 0;
|
||||
const ArrayList<int> &mapping = sea_->multiindex_mapping_[index];
|
||||
const ArrayList<int> &direct_ancestor_mapping =
|
||||
sea_->multiindex_mapping_[direct_ancestor_mapping_pos];
|
||||
for(index_t i = 0; i < dim; i++) {
|
||||
if(mapping[i] != direct_ancestor_mapping[i]) {
|
||||
position = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
tmp[index] = tmp[direct_ancestor_mapping_pos] * x_Q_to_x_q[position];
|
||||
}
|
||||
|
||||
for(index_t i = 0; i < total_num_coeffs; i++) {
|
||||
int index = traversal_order[i];
|
||||
sum += coeffs_[index] * tmp[index];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::Init(const Vector& center,
|
||||
const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
center_.Copy(center);
|
||||
order_ = -1;
|
||||
sea_ = &(ka.sea_);
|
||||
ka_ = &ka;
|
||||
|
||||
// initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
|
||||
|
||||
// copy kernel type, center, and bandwidth squared
|
||||
kernel_ = &(ka.kernel_);
|
||||
sea_ = &(ka.sea_);
|
||||
center_.Init(sea_->get_dimension());
|
||||
order_ = -1;
|
||||
ka_ = &ka;
|
||||
|
||||
// initialize coefficient array
|
||||
coeffs_.Init(sea_->get_max_total_num_coeffs());
|
||||
coeffs_.SetZero();
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
template<typename TBound>
|
||||
int MultLocalExpansion<TKernelAux>::OrderForEvaluating
|
||||
(const TBound &far_field_region,
|
||||
const TBound &local_field_region, double min_dist_sqd_regions,
|
||||
double max_dist_sqd_regions, double max_error, double *actual_error) const {
|
||||
|
||||
return ka_->OrderForEvaluatingLocal(far_field_region, local_field_region,
|
||||
min_dist_sqd_regions,
|
||||
max_dist_sqd_regions, max_error,
|
||||
actual_error);
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::TranslateFromFarField
|
||||
(const MultFarFieldExpansion<TKernelAux> &se) {
|
||||
|
||||
Vector pos_arrtmp, neg_arrtmp;
|
||||
Matrix derivative_map;
|
||||
Vector far_center;
|
||||
Vector cent_diff;
|
||||
Vector far_coeffs;
|
||||
int dimension = sea_->get_dimension();
|
||||
int far_order = se.get_order();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
|
||||
int limit;
|
||||
double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
|
||||
|
||||
// get center and coefficients for far field expansion
|
||||
far_center.Alias(*(se.get_center()));
|
||||
far_coeffs.Alias(se.get_coeffs());
|
||||
cent_diff.Init(dimension);
|
||||
|
||||
// if the order of the far field expansion is greater than the
|
||||
// local one we are adding onto, then increase the order.
|
||||
if(far_order > order_) {
|
||||
order_ = far_order;
|
||||
}
|
||||
|
||||
// compute Gaussian derivative
|
||||
limit = 2 * order_ + 1;
|
||||
derivative_map.Init(dimension, limit);
|
||||
pos_arrtmp.Init(total_num_coeffs);
|
||||
neg_arrtmp.Init(total_num_coeffs);
|
||||
|
||||
// compute center difference divided by bw_times_sqrt_two;
|
||||
for(index_t j = 0; j < dimension; j++) {
|
||||
cent_diff[j] = (center_[j] - far_center[j]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// compute required partial derivatives
|
||||
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map);
|
||||
ArrayList<int> beta_plus_alpha;
|
||||
beta_plus_alpha.Init(dimension);
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
ArrayList<int> &traversal_order = sea_->traversal_mapping_[far_order];
|
||||
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
int index_j = traversal_order[j];
|
||||
ArrayList<int> beta_mapping = sea_->get_multiindex(index_j);
|
||||
pos_arrtmp[index_j] = neg_arrtmp[index_j] = 0;
|
||||
|
||||
for(index_t k = 0; k < total_num_coeffs; k++) {
|
||||
|
||||
int index_k = traversal_order[k];
|
||||
ArrayList<int> alpha_mapping = sea_->get_multiindex(index_k);
|
||||
for(index_t d = 0; d < dimension; d++) {
|
||||
beta_plus_alpha[d] = beta_mapping[d] + alpha_mapping[d];
|
||||
}
|
||||
double derivative_factor =
|
||||
ka_->ComputePartialDerivative(derivative_map, beta_plus_alpha);
|
||||
|
||||
double prod = far_coeffs[index_k] * derivative_factor;
|
||||
|
||||
if(prod > 0) {
|
||||
pos_arrtmp[index_j] += prod;
|
||||
}
|
||||
else {
|
||||
neg_arrtmp[index_j] += prod;
|
||||
}
|
||||
} // end of k-loop
|
||||
} // end of j-loop
|
||||
|
||||
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
int index_j = traversal_order[j];
|
||||
coeffs_[index_j] += (pos_arrtmp[index_j] + neg_arrtmp[index_j]) *
|
||||
C_k_neg[index_j];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TKernelAux>
|
||||
void MultLocalExpansion<TKernelAux>::TranslateToLocal(MultLocalExpansion &se) {
|
||||
|
||||
// if no local coefficients have formed, then nothing to translate
|
||||
if(order_ < 0) {
|
||||
return;
|
||||
}
|
||||
|
||||
// get the center and the order and the total number of coefficients of
|
||||
// the expansion we are translating from. Also get coefficients we
|
||||
// are translating
|
||||
Vector new_center;
|
||||
new_center.Alias(*(se.get_center()));
|
||||
int prev_order = se.get_order();
|
||||
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
|
||||
const ArrayList < int > *upper_mapping_index =
|
||||
sea_->get_upper_mapping_index();
|
||||
Vector new_coeffs;
|
||||
new_coeffs.Alias(se.get_coeffs());
|
||||
|
||||
// dimension
|
||||
int dim = sea_->get_dimension();
|
||||
|
||||
// temporary variable
|
||||
ArrayList<int> tmp_storage;
|
||||
tmp_storage.Init(dim);
|
||||
|
||||
// sqrt two times bandwidth
|
||||
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
|
||||
|
||||
// center difference between the old center and the new one
|
||||
Vector center_diff;
|
||||
center_diff.Init(dim);
|
||||
for(index_t d = 0; d < dim; d++) {
|
||||
center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
|
||||
}
|
||||
|
||||
// set to the new order if the order of the expansion we are translating
|
||||
// from is higher
|
||||
if(prev_order < order_) {
|
||||
se.set_order(order_);
|
||||
}
|
||||
|
||||
// inverse multiindex factorials
|
||||
Vector C_k;
|
||||
C_k.Alias(sea_->get_inv_multiindex_factorials());
|
||||
|
||||
// get the order of traversal for the given order of approximation
|
||||
const ArrayList<int> &traversal_order = sea_->traversal_mapping_[order_];
|
||||
|
||||
// do the actual translation
|
||||
for(index_t j = 0; j < total_num_coeffs; j++) {
|
||||
|
||||
int index_j = traversal_order[j];
|
||||
const ArrayList<int> &alpha_mapping = sea_->get_multiindex(index_j);
|
||||
const ArrayList <int> &upper_mappings_for_alpha =
|
||||
upper_mapping_index[index_j];
|
||||
double pos_coeffs = 0;
|
||||
double neg_coeffs = 0;
|
||||
|
||||
for(index_t k = 0; k < upper_mappings_for_alpha.size(); k++) {
|
||||
|
||||
if(upper_mappings_for_alpha[k] >= total_num_coeffs) {
|
||||
break;
|
||||
}
|
||||
|
||||
const ArrayList<int> &beta_mapping =
|
||||
sea_->get_multiindex(upper_mappings_for_alpha[k]);
|
||||
int flag = 0;
|
||||
double diff1 = 1.0;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
tmp_storage[l] = beta_mapping[l] - alpha_mapping[l];
|
||||
|
||||
if(tmp_storage[l] < 0) {
|
||||
flag = 1;
|
||||
break;
|
||||
}
|
||||
} // end of looping over dimension
|
||||
|
||||
if(flag)
|
||||
continue;
|
||||
|
||||
for(index_t l = 0; l < dim; l++) {
|
||||
diff1 *= pow(center_diff[l], tmp_storage[l]);
|
||||
}
|
||||
|
||||
double prod = coeffs_[upper_mappings_for_alpha[k]] * diff1 *
|
||||
sea_->get_n_multichoose_k_by_pos
|
||||
(upper_mappings_for_alpha[k], index_j);
|
||||
|
||||
if(prod > 0) {
|
||||
pos_coeffs += prod;
|
||||
}
|
||||
else {
|
||||
neg_coeffs += prod;
|
||||
}
|
||||
|
||||
} // end of k loop
|
||||
|
||||
new_coeffs[index_j] += pos_coeffs + neg_coeffs;
|
||||
} // end of j loop
|
||||
}
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user