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

Update series_expansion methods to use armadillo and std:: components.

Browse Source
This commit is contained in:
Ryan Curtin
2011-01-31 19:16:13 +00:00
parent 6bd0293fa2
commit 4e98f49684
13 changed files with 626 additions and 803 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion.h
+18 -26
View File
@@ -12,7 +12,8 @@
#ifndef FARFIELD_EXPANSION
#define FARFIELD_EXPANSION
#include "fastlib/fastlib.h"
#include <fastlib/fastlib.h>
#include "kernel_aux.h"
#include "series_expansion_aux.h"
@@ -38,10 +39,10 @@ class FarFieldExpansion {
////////// Private Member Variables //////////
/** @brief The center of the expansion. */
Vector center_;
arma::vec center_;
/** @brief The coefficients. */
Vector coeffs_;
arma::vec coeffs_;
/** @brief The order of the expansion. */
int order_;
@@ -61,12 +62,6 @@ class FarFieldExpansion {
*/
const typename TKernelAux::TSeriesExpansionAux *sea_;
OT_DEF(FarFieldExpansion) {
OT_MY_OBJECT(center_);
OT_MY_OBJECT(coeffs_);
OT_MY_OBJECT(order_);
}
public:
////////// Getters/Setters //////////
@@ -82,16 +77,16 @@ class FarFieldExpansion {
*
* @return The center of expansion for the current far-field expansion.
*/
Vector* get_center() { return &center_; }
const Vector* get_center() const { return &center_; }
arma::vec& get_center() { return center_; }
const arma::vec& get_center() const { return center_; }
/** @brief Gets the set of far-field coefficients.
*
* @return The const reference to the vector containing the
* far-field coefficients.
*/
const Vector& get_coeffs() const { return coeffs_; }
arma::vec& get_coeffs() { return coeffs_; }
const arma::vec& get_coeffs() const { return coeffs_; }
/** @brief Gets the approximation order.
*
@@ -122,11 +117,8 @@ class FarFieldExpansion {
* will be copied to the center of the given far-field
* expansion object.
*/
void set_center(const Vector &center) {
for(index_t i = 0; i < center.length(); i++) {
center_[i] = center[i];
}
void set_center(const arma::vec& center) {
center_ = center;
}
////////// User-level Functions //////////
@@ -151,7 +143,7 @@ class FarFieldExpansion {
* @param order The order up to which the far-field moments should be
* accumulated up to.
*/
void Accumulate(const Vector &reference_point, double weight, int order);
void Accumulate(const arma::vec& reference_point, double weight, int order);
/** @brief Accumulates the far field moment represented by the given
* reference data into the coefficients.
@@ -178,24 +170,24 @@ class FarFieldExpansion {
* @param order The order up to which the far-field moments should be
* accumulated up to.
*/
void AccumulateCoeffs(const Matrix& data, const Vector& weights,
void AccumulateCoeffs(const arma::mat& data, const arma::vec& weights,
int begin, int end, int order);
/** @brief Refine the far field moment that has been computed before
* up to a new order.
*/
void RefineCoeffs(const Matrix& data, const Vector& weights,
void RefineCoeffs(const arma::mat& data, const arma::vec& weights,
int begin, int end, int order);
/** @brief Evaluates the far-field coefficients at the given point.
*/
double EvaluateField(const Matrix& data, int row_num, int order) const;
double EvaluateField(const arma::mat& data, int row_num, int order) const;
double EvaluateField(const double *x_q, int order) const;
/** @brief Evaluates the two-way convolution mixed with exhaustive
* computations with two other far field expansions.
*/
double MixField(const Matrix &data, int node1_begin, int node1_end,
double MixField(const arma::mat& data, int node1_begin, int node1_end,
int node2_begin, int node2_end, const FarFieldExpansion &fe2,
const FarFieldExpansion &fe3, int order2, int order3) const;
@@ -214,14 +206,14 @@ class FarFieldExpansion {
/** @brief Initializes the current far field expansion object with
* the given center.
*/
void Init(const Vector& center, const TKernelAux &ka);
void Init(const arma::vec& center, const TKernelAux &ka);
void Init(const TKernelAux &ka);
template<typename TBound>
int OrderForConvolving(const TBound &far_field_region,
const Vector &far_field_region_centroid,
const arma::vec &far_field_region_centroid,
const TBound &local_field_region,
const Vector &local_field_region_centroid,
const arma::vec &local_field_region_centroid,
double min_dist_sqd_regions,
double max_dist_dsqd_regons,
double max_error, double *actual_error) const;
fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion_impl.h
+176 -217
View File
@@ -6,33 +6,31 @@
#define FARFIELD_EXPANSION_IMPL_H
template<typename TKernelAux>
void FarFieldExpansion<TKernelAux>::Accumulate(const Vector &v, double weight,
void FarFieldExpansion<TKernelAux>::Accumulate(const arma::vec& v, double weight,
int order) {
int dim = v.length();
int dim = v.n_elem;
int total_num_coeffs = sea_->get_total_num_coeffs(order);
Vector tmp;
arma::vec tmp;
int r, i, j, k, t, tail;
GenVector<short int> heads;
Vector x_r;
arma::Col<short int> heads;
arma::vec x_r;
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
// initialize temporary variables
tmp.Init(total_num_coeffs);
heads.Init(dim + 1);
x_r.Init(dim);
Vector pos_coeffs;
Vector neg_coeffs;
pos_coeffs.Init(total_num_coeffs);
pos_coeffs.SetZero();
neg_coeffs.Init(total_num_coeffs);
neg_coeffs.SetZero();
tmp.set_size(total_num_coeffs);
heads.set_size(dim + 1);
x_r.set_size(dim);
arma::vec pos_coeffs;
arma::vec neg_coeffs;
pos_coeffs.zeros(total_num_coeffs);
neg_coeffs.zeros(total_num_coeffs);
// set to new order if greater
if(order_ < order) {
order_ = order;
}
Vector C_k;
arma::vec& C_k;
// Calculate the coordinate difference between the ref point and the
// centroid.
@@ -41,7 +39,7 @@ void FarFieldExpansion<TKernelAux>::Accumulate(const Vector &v, double weight,
}
// initialize heads
heads.SetZero();
heads.zeros();
heads[dim] = SHRT_MAX;
tmp[0] = 1.0;
@@ -70,7 +68,7 @@ void FarFieldExpansion<TKernelAux>::Accumulate(const Vector &v, double weight,
}
// get multiindex factors
C_k.Alias(sea_->get_inv_multiindex_factorials());
C_k = sea_->get_inv_multiindex_factorials();
for(r = 0; r < total_num_coeffs; r++) {
coeffs_[r] += (pos_coeffs[r] + neg_coeffs[r]) * C_k[r];
@@ -78,47 +76,44 @@ void FarFieldExpansion<TKernelAux>::Accumulate(const Vector &v, double weight,
}
template<typename TKernelAux>
void FarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
const Vector& weights,
void FarFieldExpansion<TKernelAux>::AccumulateCoeffs(const arma::mat& data,
const arma::vec& weights,
int begin, int end,
int order) {
int dim = data.n_rows();
int dim = data.n_rows;
int total_num_coeffs = sea_->get_total_num_coeffs(order);
Vector tmp;
arma::vec tmp;
int r, i, j, k, t, tail;
GenVector<short int> heads;
Vector x_r;
arma::Col<short int> heads;
arma::vec x_r;
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
// initialize temporary variables
tmp.Init(total_num_coeffs);
heads.Init(dim + 1);
x_r.Init(dim);
Vector pos_coeffs;
Vector neg_coeffs;
pos_coeffs.Init(total_num_coeffs);
pos_coeffs.SetZero();
neg_coeffs.Init(total_num_coeffs);
neg_coeffs.SetZero();
tmp.set_size(total_num_coeffs);
heads.set_size(dim + 1);
x_r.set_size(dim);
arma::vec pos_coeffs;
arma::vec neg_coeffs;
pos_coeffs.zeros(total_num_coeffs);
neg_coeffs.zeros(total_num_coeffs);
// set to new order if greater
if(order_ < order) {
order_ = order;
}
Vector C_k;
// Repeat for each reference point in this reference node.
for(r = begin; r < end; r++) {
// Calculate the coordinate difference between the ref point and the
// centroid.
for(i = 0; i < dim; i++) {
x_r[i] = (data.get(i, r) - center_[i]) / bandwidth_factor;
x_r[i] = (data(i, r) - center_[i]) / bandwidth_factor;
}
// initialize heads
heads.SetZero();
heads.zeros();
heads[dim] = SHRT_MAX;
tmp[0] = 1.0;
@@ -149,7 +144,7 @@ void FarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
} // End of looping through each reference point
// get multiindex factors
C_k.Alias(sea_->get_inv_multiindex_factorials());
const arma::vec& C_k = sea_->get_inv_multiindex_factorials();
for(r = 0; r < total_num_coeffs; r++) {
coeffs_[r] += (pos_coeffs[r] + neg_coeffs[r]) * C_k[r];
@@ -157,8 +152,8 @@ void FarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
}
template<typename TKernelAux>
void FarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
const Vector& weights,
void FarFieldExpansion<TKernelAux>::RefineCoeffs(const arma::mat& data,
const arma::vec& weights,
int begin, int end,
int order) {
@@ -168,22 +163,19 @@ void FarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
return;
}
int dim = data.n_rows();
int dim = data.n_rows;
int old_total_num_coeffs = sea_->get_total_num_coeffs(order_);
int total_num_coeffs = sea_->get_total_num_coeffs(order);
double tmp;
int r, i, j;
Vector x_r;
arma::vec x_r(dim);
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
// initialize temporary variables
x_r.Init(dim);
Vector pos_coeffs;
Vector neg_coeffs;
pos_coeffs.Init(total_num_coeffs);
pos_coeffs.SetZero();
neg_coeffs.Init(total_num_coeffs);
neg_coeffs.SetZero();
arma::vec pos_coeffs;
arma::vec neg_coeffs;
pos_coeffs.zeros(total_num_coeffs);
neg_coeffs.zeros(total_num_coeffs);
// if we already have the order of approximation, then return.
if(order_ >= order) {
@@ -193,7 +185,7 @@ void FarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
order_ = order;
}
Vector C_k;
const arma::vec& C_k = sea_->get_inv_multiindex_factorials();
// Repeat for each reference point in this reference node.
for(r = begin; r < end; r++) {
@@ -201,12 +193,12 @@ void FarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
// Calculate the coordinate difference between the ref point and the
// centroid.
for(i = 0; i < dim; i++) {
x_r[i] = (data.get(i, r) - center_[i]) / bandwidth_factor;
x_r[i] = (data(i, r) - center_[i]) / bandwidth_factor;
}
// compute in bruteforce way
for(i = old_total_num_coeffs; i < total_num_coeffs; i++) {
const ArrayList<short int> &mapping = sea_->get_multiindex(i);
const std::vector<short int> &mapping = sea_->get_multiindex(i);
tmp = 1;
for(j = 0; j < dim; j++) {
@@ -225,19 +217,16 @@ void FarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
} // End of looping through each reference point
// get multiindex factors
C_k.Alias(sea_->get_inv_multiindex_factorials());
for(r = old_total_num_coeffs; r < total_num_coeffs; r++) {
coeffs_[r] = (pos_coeffs[r] + neg_coeffs[r]) * C_k[r];
}
}
template<typename TKernelAux>
double FarFieldExpansion<TKernelAux>::EvaluateField(const Matrix& data,
double FarFieldExpansion<TKernelAux>::EvaluateField(const arma::mat& data,
int row_num,
int order) const {
return EvaluateField(data.GetColumnPtr(row_num), order);
return EvaluateField(data.unsafe_col(row_num).memptr(), order);
}
template<typename TKernelAux>
@@ -259,16 +248,14 @@ double FarFieldExpansion<TKernelAux>::EvaluateField(const double *x_q,
double multipole_sum = 0;
// computed derivative map
Matrix derivative_map;
ka_->AllocateDerivativeMap(dim, order, &derivative_map);
arma::mat derivative_map;
ka_->AllocateDerivativeMap(dim, order, derivative_map);
// temporary variable
Vector arrtmp;
arrtmp.Init(total_num_coeffs);
arma::vec arrtmp(total_num_coeffs);
// (x_q - x_R) scaled by bandwidth
Vector x_q_minus_x_R;
x_q_minus_x_R.Init(dim);
arma::vec x_q_minus_x_R(dim);
// compute (x_q - x_R) / (sqrt(2h^2))
for(index_t d = 0; d < dim; d++) {
@@ -276,11 +263,11 @@ double FarFieldExpansion<TKernelAux>::EvaluateField(const double *x_q,
}
// compute deriative maps based on coordinate difference.
ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, &derivative_map, order);
ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map, 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++) {
const ArrayList<short int> &mapping = sea_->get_multiindex(j);
const std::vector<short int> &mapping = sea_->get_multiindex(j);
double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
double prod = coeffs_[j] * arrtmp;
@@ -297,7 +284,7 @@ double FarFieldExpansion<TKernelAux>::EvaluateField(const double *x_q,
}
template<typename TKernelAux>
double FarFieldExpansion<TKernelAux>::MixField(const Matrix &data,
double FarFieldExpansion<TKernelAux>::MixField(const arma::mat &data,
int node1_begin, int node1_end,
int node2_begin, int node2_end,
const FarFieldExpansion &fe2,
@@ -307,18 +294,17 @@ double FarFieldExpansion<TKernelAux>::MixField(const Matrix &data,
// bandwidth factor and multiindex mapping stuffs
double result;
double bandwidth_factor = ka_->BandwidthFactor(bandwidth_sq());
const ArrayList<short int> *multiindex_mapping =
const std::vector<short int>& multiindex_mapping =
sea_->get_multiindex_mapping();
const ArrayList<short int> *lower_mapping_index =
const std::vector<short int>& lower_mapping_index =
sea_->get_lower_mapping_index();
// get the total number of coefficients and coefficients
int total_num_coeffs2 = sea_->get_total_num_coeffs(order2);
int total_num_coeffs3 = sea_->get_total_num_coeffs(order3);
int dim = sea_->get_dimension();
Vector coeffs2, coeffs3;
coeffs2.Alias(fe2.get_coeffs());
coeffs3.Alias(fe3.get_coeffs());
arma::vec& coeffs2 = fe2.get_coeffs();
arma::vec& coeffs3 = fe3.get_coeffs();
// actual accumulated sum
double neg_sum = 0;
@@ -330,65 +316,58 @@ double FarFieldExpansion<TKernelAux>::MixField(const Matrix &data,
double xi_xI, xj_xJ, diff;
// temporary array
ArrayList<short int> beta_gamma_nu_eta_mapping;
ArrayList<short int> beta_nu_mapping;
ArrayList<short int> gamma_eta_mapping;
beta_nu_mapping.Init(dim);
gamma_eta_mapping.Init(dim);
beta_gamma_nu_eta_mapping.Init(dim);
std::vector<short int> beta_gamma_nu_eta_mapping;
std::vector<short int> beta_nu_mapping;
std::vector<short int> gamma_eta_mapping;
beta_nu_mapping.reserve(dim);
gamma_eta_mapping.reserve(dim);
beta_gamma_nu_eta_mapping.reserve(dim);
// partial derivatives table
Matrix derivative_map_beta;
ka_->AllocateDerivativeMap(dim, order2, &derivative_map_beta);
Matrix derivative_map_gamma;
ka_->AllocateDerivativeMap(dim, order3, &derivative_map_gamma);
arma::mat derivative_map_beta;
ka_->AllocateDerivativeMap(dim, order2, derivative_map_beta);
arma::mat derivative_map_gamma;
ka_->AllocateDerivativeMap(dim, order3, derivative_map_gamma);
// compute center differences and complete the table of partial derivatives
Vector xI_xK, xJ_xK;
xI_xK.Init(dim);
xJ_xK.Init(dim);
Vector xJ_center, xK_center;
xJ_center.Alias(*(fe2.get_center()));
xK_center.Alias(*(fe3.get_center()));
arma::vec xI_xK(dim), xJ_xK(dim);
arma::vec& xJ_center = fe2.get_center();
arma::vec& xK_center = fe3.get_center();
for(index_t d = 0; d < dim; d++) {
xI_xK[d] = (center_[d] - xK_center[d]) / bandwidth_factor;
xJ_xK[d] = (xJ_center[d] - xK_center[d]) / bandwidth_factor;
}
ka_->ComputeDirectionalDerivatives(xI_xK, &derivative_map_beta, order2);
ka_->ComputeDirectionalDerivatives(xJ_xK, &derivative_map_gamma, order3);
ka_->ComputeDirectionalDerivatives(xI_xK, derivative_map_beta, order2);
ka_->ComputeDirectionalDerivatives(xJ_xK, derivative_map_gamma, order3);
// inverse factorials
Vector inv_multiindex_factorials;
inv_multiindex_factorials.Alias(sea_->get_inv_multiindex_factorials());
arma::vec& inv_multiindex_factorials = sea_->get_inv_multiindex_factorials();
// precompute pairwise kernel values between node i and node j
Matrix exhaustive_ij;
exhaustive_ij.Init(node1_end - node1_begin, node2_end - node2_begin);
arma::mat exhaustive_ij(node1_end - node1_begin, node2_end - node2_begin);
for(index_t i = node1_begin; i < node1_end; i++) {
const double *i_col = data.GetColumnPtr(i);
arma::vec i_col = data.unsafe_col(i);
for(index_t j = node2_begin; j < node2_end; j++) {
const double *j_col = data.GetColumnPtr(j);
arma::vec j_col = data.unsafe_col(j);
exhaustive_ij.set
(i - node1_begin, j - node2_begin,
kernel_->EvalUnnormOnSq(la::DistanceSqEuclidean(data.n_rows(),
i_col, j_col)));
exhaustive_ij(i - node1_begin, j - node2_begin) =
kernel_->EvalUnnormOnSq(la::DistanceSqEuclidean(i_col, j_col));
}
}
// main loop
for(index_t beta = 0; beta < total_num_coeffs2; beta++) {
const ArrayList <short int> &beta_mapping = multiindex_mapping[beta];
const ArrayList <short int> &lower_mappings_for_beta =
const std::vector<short int> &beta_mapping = multiindex_mapping[beta];
const std::vector<short int> &lower_mappings_for_beta =
lower_mapping_index[beta];
double beta_derivative = ka_->ComputePartialDerivative
(derivative_map_beta, beta_mapping);
for(index_t nu = 0; nu < lower_mappings_for_beta.size(); nu++) {
const ArrayList<short int> &nu_mapping =
const std::vector<short int> &nu_mapping =
multiindex_mapping[lower_mappings_for_beta[nu]];
// beta - nu
@@ -398,19 +377,19 @@ double FarFieldExpansion<TKernelAux>::MixField(const Matrix &data,
for(index_t gamma = 0; gamma < total_num_coeffs3; gamma++) {
const ArrayList <short int> &gamma_mapping = multiindex_mapping[gamma];
const ArrayList <short int> &lower_mappings_for_gamma =
const std::vector<short int> &gamma_mapping = multiindex_mapping[gamma];
const std::vector<short int> &lower_mappings_for_gamma =
lower_mapping_index[gamma];
double gamma_derivative = ka_->ComputePartialDerivative
(derivative_map_gamma, gamma_mapping);
for(index_t eta = 0; eta < lower_mappings_for_gamma.size();
eta++){
eta++) {
// add up alpha, mu, eta and beta, gamma, nu, eta
int sign = 0;
const ArrayList<short int> &eta_mapping =
const std::vector<short int> &eta_mapping =
multiindex_mapping[lower_mappings_for_gamma[eta]];
for(index_t d = 0; d < dim; d++) {
@@ -439,7 +418,7 @@ double FarFieldExpansion<TKernelAux>::MixField(const Matrix &data,
inv_multiindex_factorials
[sea_->ComputeMultiindexPosition(nu_mapping)];
for(index_t d = 0; d < dim; d++) {
diff = (data.get(d, i) - center_[d]) / bandwidth_factor;
diff = (data(d, i) - center_[d]) / bandwidth_factor;
xi_xI *= pow(diff, nu_mapping[d]);
}
@@ -448,7 +427,7 @@ double FarFieldExpansion<TKernelAux>::MixField(const Matrix &data,
xj_xJ = inv_multiindex_factorials
[sea_->ComputeMultiindexPosition(eta_mapping)];
for(index_t d = 0; d < dim; d++) {
diff = (data.get(d, j) - xJ_center[d]) / bandwidth_factor;
diff = (data(d, j) - xJ_center[d]) / bandwidth_factor;
xj_xJ *= pow(diff, eta_mapping[d]);
}
@@ -457,8 +436,7 @@ double FarFieldExpansion<TKernelAux>::MixField(const Matrix &data,
(sea_->ComputeMultiindexPosition(beta_gamma_nu_eta_mapping),
sea_->ComputeMultiindexPosition(beta_nu_mapping)) *
beta_derivative * gamma_derivative * xi_xI * xj_xJ *
moment_k * exhaustive_ij.get
(i - node1_begin, j - node2_begin);
moment_k * exhaustive_ij(i - node1_begin, j - node2_begin);
if(result > 0) {
pos_sum += result;
@@ -485,16 +463,15 @@ double FarFieldExpansion<TKernelAux>::ConvolveField
// The bandwidth factor and the multiindex mapping stuffs.
double bandwidth_factor = ka_->BandwidthFactor(bandwidth_sq());
const ArrayList<short int> *multiindex_mapping =
const std::vector<short int>& multiindex_mapping =
sea_->get_multiindex_mapping();
const ArrayList<short int> *lower_mapping_index =
const std::vector<short int>& lower_mapping_index =
sea_->get_lower_mapping_index();
// Get the total number of coefficients and the coefficient themselves.
int total_num_coeffs = sea_->get_total_num_coeffs(order);
int dim = sea_->get_dimension();
Vector coeffs2;
coeffs2.Alias(fe.get_coeffs());
arma::vec& coeffs2 = fe.get_coeffs();
// Actual accumulated sum.
double neg_sum = 0;
@@ -502,41 +479,38 @@ double FarFieldExpansion<TKernelAux>::ConvolveField
double sum = 0;
// The partial derivatives table.
Matrix derivative_map_alpha;
ka_->AllocateDerivativeMap(dim, order, &derivative_map_alpha);
arma::mat derivative_map_alpha;
ka_->AllocateDerivativeMap(dim, order, derivative_map_alpha);
// Compute the center difference and its table of partial
// derivatives.
Vector xI_xJ;
xI_xJ.Init(dim);
Vector xJ_center;
xJ_center.Alias(*(fe.get_center()));
arma::vec xI_xJ(dim);
arma::vec& xJ_center = fe.get_center();
for(index_t d = 0; d < dim; d++) {
xI_xJ[d] = (center_[d] - xJ_center[d]) / bandwidth_factor;
}
ka_->ComputeDirectionalDerivatives(xI_xJ, &derivative_map_alpha, order);
ka_->ComputeDirectionalDerivatives(xI_xJ, derivative_map_alpha, order);
// The inverse factorials.
Vector inv_multiindex_factorials;
inv_multiindex_factorials.Alias(sea_->get_inv_multiindex_factorials());
arma::vec& inv_multiindex_factorials = sea_->get_inv_multiindex_factorials();
// The temporary space for computing the difference of two mappings.
ArrayList<short int> alpha_minus_beta_mapping;
alpha_minus_beta_mapping.Init(dim);
std::vector<short int> alpha_minus_beta_mapping;
alpha_minus_beta_mapping.reserve(dim);
// The main loop.
for(index_t alpha = 0; alpha < total_num_coeffs; alpha++) {
const ArrayList <short int> &alpha_mapping = multiindex_mapping[alpha];
const ArrayList <short int> &lower_mappings_for_alpha =
const std::vector<short int> &alpha_mapping = multiindex_mapping[alpha];
const std::vector<short int> &lower_mappings_for_alpha =
lower_mapping_index[alpha];
double alpha_derivative = ka_->ComputePartialDerivative
(derivative_map_alpha, alpha_mapping);
for(index_t beta = 0; beta < lower_mappings_for_alpha.size(); beta++) {
const ArrayList <short int> &beta_mapping =
const std::vector<short int> &beta_mapping =
multiindex_mapping[lower_mappings_for_alpha[beta]];
double n_choose_k_factor = sea_->get_n_multichoose_k_by_pos
@@ -583,9 +557,9 @@ double FarFieldExpansion<TKernelAux>::ConvolveField
// bandwidth factor and multiindex mapping stuffs
double result;
double bandwidth_factor = ka_->BandwidthFactor(bandwidth_sq());
const ArrayList<short int> *multiindex_mapping =
const std::vector<short int>& multiindex_mapping =
sea_->get_multiindex_mapping();
const ArrayList<short int> *lower_mapping_index =
const std::vector<short int>& lower_mapping_index =
sea_->get_lower_mapping_index();
// get the total number of coefficients and coefficients
@@ -593,9 +567,8 @@ double FarFieldExpansion<TKernelAux>::ConvolveField
int total_num_coeffs2 = sea_->get_total_num_coeffs(order2);
int total_num_coeffs3 = sea_->get_total_num_coeffs(order3);
int dim = sea_->get_dimension();
Vector coeffs2, coeffs3;
coeffs2.Alias(fe2.get_coeffs());
coeffs3.Alias(fe3.get_coeffs());
arma::vec& coeffs2 = fe2.get_coeffs();
arma::vec& coeffs3 = fe3.get_coeffs();
// actual accumulated sum
double neg_sum = 0;
@@ -606,61 +579,56 @@ double FarFieldExpansion<TKernelAux>::ConvolveField
double moment_i, moment_j, moment_k;
// temporary array
ArrayList<short int> mu_nu_mapping;
ArrayList<short int> alpha_mu_eta_mapping;
ArrayList<short int> beta_gamma_nu_eta_mapping;
ArrayList<short int> alpha_mu_mapping;
ArrayList<short int> beta_nu_mapping;
ArrayList<short int> gamma_eta_mapping;
alpha_mu_mapping.Init(dim);
beta_nu_mapping.Init(dim);
gamma_eta_mapping.Init(dim);
mu_nu_mapping.Init(dim);
alpha_mu_eta_mapping.Init(dim);
beta_gamma_nu_eta_mapping.Init(dim);
std::vector<short int> mu_nu_mapping;
std::vector<short int> alpha_mu_eta_mapping;
std::vector<short int> beta_gamma_nu_eta_mapping;
std::vector<short int> alpha_mu_mapping;
std::vector<short int> beta_nu_mapping;
std::vector<short int> gamma_eta_mapping;
alpha_mu_mapping.reserve(dim);
beta_nu_mapping.reserve(dim);
gamma_eta_mapping.reserve(dim);
mu_nu_mapping.reserve(dim);
alpha_mu_eta_mapping.reserve(dim);
beta_gamma_nu_eta_mapping.reserve(dim);
// partial derivatives table
Matrix derivative_map_alpha;
ka_->AllocateDerivativeMap(dim, order1, &derivative_map_alpha);
Matrix derivative_map_beta;
ka_->AllocateDerivativeMap(dim, order2, &derivative_map_beta);
Matrix derivative_map_gamma;
ka_->AllocateDerivativeMap(dim, order3, &derivative_map_gamma);
arma::mat derivative_map_alpha;
ka_->AllocateDerivativeMap(dim, order1, derivative_map_alpha);
arma::mat derivative_map_beta;
ka_->AllocateDerivativeMap(dim, order2, derivative_map_beta);
arma::mat derivative_map_gamma;
ka_->AllocateDerivativeMap(dim, order3, derivative_map_gamma);
// compute center differences and complete the table of partial derivatives
Vector xI_xJ, xI_xK, xJ_xK;
xI_xJ.Init(dim);
xI_xK.Init(dim);
xJ_xK.Init(dim);
Vector xJ_center, xK_center;
xJ_center.Alias(*(fe2.get_center()));
xK_center.Alias(*(fe3.get_center()));
arma::vec xI_xJ(dim), xI_xK(dim), xJ_xK(dim);
arma::vec& xJ_center = fe2.get_center();
arma::vec& xK_center = fe3.get_center();
for(index_t d = 0; d < dim; d++) {
xI_xJ[d] = (center_[d] - xJ_center[d]) / bandwidth_factor;
xI_xK[d] = (center_[d] - xK_center[d]) / bandwidth_factor;
xJ_xK[d] = (xJ_center[d] - xK_center[d]) / bandwidth_factor;
}
ka_->ComputeDirectionalDerivatives(xI_xJ, &derivative_map_alpha, order1);
ka_->ComputeDirectionalDerivatives(xI_xK, &derivative_map_beta, order2);
ka_->ComputeDirectionalDerivatives(xJ_xK, &derivative_map_gamma, order3);
ka_->ComputeDirectionalDerivatives(xI_xJ, derivative_map_alpha, order1);
ka_->ComputeDirectionalDerivatives(xI_xK, derivative_map_beta, order2);
ka_->ComputeDirectionalDerivatives(xJ_xK, derivative_map_gamma, order3);
// inverse factorials
Vector inv_multiindex_factorials;
inv_multiindex_factorials.Alias(sea_->get_inv_multiindex_factorials());
arma::vec& inv_multiindex_factorials = sea_->get_inv_multiindex_factorials();
// main loop
for(index_t alpha = 0; alpha < total_num_coeffs1; alpha++) {
const ArrayList <short int> &alpha_mapping = multiindex_mapping[alpha];
const ArrayList <short int> &lower_mappings_for_alpha =
const std::vector<short int>& alpha_mapping = multiindex_mapping[alpha];
const std::vector<short int>& lower_mappings_for_alpha =
lower_mapping_index[alpha];
double alpha_derivative = ka_->ComputePartialDerivative
(derivative_map_alpha, alpha_mapping);
for(index_t mu = 0; mu < lower_mappings_for_alpha.size(); mu++) {
const ArrayList <short int> &mu_mapping =
const std::vector<short int>& mu_mapping =
multiindex_mapping[lower_mappings_for_alpha[mu]];
// alpha - mu
@@ -670,15 +638,15 @@ double FarFieldExpansion<TKernelAux>::ConvolveField
for(index_t beta = 0; beta < total_num_coeffs2; beta++) {
const ArrayList <short int> &beta_mapping = multiindex_mapping[beta];
const ArrayList <short int> &lower_mappings_for_beta =
const std::vector<short int>& beta_mapping = multiindex_mapping[beta];
const std::vector<short int>& lower_mappings_for_beta =
lower_mapping_index[beta];
double beta_derivative = ka_->ComputePartialDerivative
(derivative_map_beta, beta_mapping);
for(index_t nu = 0; nu < lower_mappings_for_beta.size(); nu++) {
const ArrayList<short int> &nu_mapping =
const std::vector<short int> &nu_mapping =
multiindex_mapping[lower_mappings_for_beta[nu]];
// mu + nu and beta - nu
@@ -689,20 +657,20 @@ double FarFieldExpansion<TKernelAux>::ConvolveField
for(index_t gamma = 0; gamma < total_num_coeffs3; gamma++) {
const ArrayList <short int> &gamma_mapping =
const std::vector<short int> &gamma_mapping =
multiindex_mapping[gamma];
const ArrayList <short int> &lower_mappings_for_gamma =
const std::vector<short int> &lower_mappings_for_gamma =
lower_mapping_index[gamma];
double gamma_derivative = ka_->ComputePartialDerivative
(derivative_map_gamma, gamma_mapping);
for(index_t eta = 0; eta < lower_mappings_for_gamma.size();
eta++){
eta++) {
// add up alpha, mu, eta and beta, gamma, nu, eta
int sign = 0;
const ArrayList<short int> &eta_mapping =
const std::vector<short int>& eta_mapping =
multiindex_mapping[lower_mappings_for_gamma[eta]];
for(index_t d = 0; d < dim; d++) {
@@ -766,19 +734,18 @@ double FarFieldExpansion<TKernelAux>::ConvolveField
}
template<typename TKernelAux>
void FarFieldExpansion<TKernelAux>::Init(const Vector& center,
void FarFieldExpansion<TKernelAux>::Init(const arma::vec& center,
const TKernelAux &ka) {
// copy kernel type, center, and bandwidth squared
kernel_ = &(ka.kernel_);
center_.Copy(center);
center_ = center;
order_ = -1;
sea_ = &(ka.sea_);
ka_ = &ka;
// initialize coefficient array
coeffs_.Init(sea_->get_max_total_num_coeffs());
coeffs_.SetZero();
coeffs_.zeros(sea_->get_max_total_num_coeffs());
}
template<typename TKernelAux>
@@ -788,20 +755,18 @@ void FarFieldExpansion<TKernelAux>::Init(const TKernelAux &ka) {
kernel_ = &(ka.kernel_);
order_ = -1;
sea_ = &(ka.sea_);
center_.Init(sea_->get_dimension());
center_.SetZero();
center_.zeros(sea_->get_dimension());
ka_ = &ka;
// initialize coefficient array
coeffs_.Init(sea_->get_max_total_num_coeffs());
coeffs_.SetZero();
coeffs_.zeros(sea_->get_max_total_num_coeffs());
}
template<typename TKernelAux>
template<typename TBound>
int FarFieldExpansion<TKernelAux>::OrderForConvolving
(const TBound &far_field_region, const Vector &far_field_region_centroid,
const TBound &local_field_region, const Vector &local_field_region_centroid,
(const TBound &far_field_region, const arma::vec& far_field_region_centroid,
const TBound &local_field_region, const arma::vec& local_field_region_centroid,
double min_dist_sqd_regions, double max_dist_sqd_regions, double max_error,
double *actual_error) const {
@@ -853,7 +818,7 @@ void FarFieldExpansion<TKernelAux>::PrintDebug(const char *name,
fprintf(stream, "Far field expansion\n");
fprintf(stream, "Center: ");
for (index_t i = 0; i < center_.length(); i++) {
for (index_t i = 0; i < center_.n_elem; i++) {
fprintf(stream, "%g ", center_[i]);
}
fprintf(stream, "\n");
@@ -867,7 +832,7 @@ void FarFieldExpansion<TKernelAux>::PrintDebug(const char *name,
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<short int> &mapping = sea_->get_multiindex(i);
const std::vector<short int> &mapping = sea_->get_multiindex(i);
fprintf(stream, "%g ", coeffs_[i]);
fprintf(stream, "(-1)^(");
@@ -899,24 +864,20 @@ void FarFieldExpansion<TKernelAux>::TranslateFromFarField
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 <short int> *multiindex_mapping =
const arma::vec& prev_coeffs = se.get_coeffs();
const arma::vec& prev_center = se.get_center();
const std::vector<std::vector<short int> >& multiindex_mapping =
sea_->get_multiindex_mapping();
const ArrayList <short int> *lower_mapping_index =
const std::vector<std::vector<short int> >& lower_mapping_index =
sea_->get_lower_mapping_index();
ArrayList <short int> tmp_storage;
Vector center_diff;
Vector inv_multiindex_factorials;
center_diff.Init(dim);
std::vector<short int> tmp_storage;
arma::vec center_diff(dim);
const arma::vec& inv_multiindex_factorials =
sea_->get_inv_multiindex_factorials();
// 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());
tmp_storage.reserve(sea_->get_dimension());
// no coefficients can be translated
if(order == -1)
@@ -931,15 +892,15 @@ void FarFieldExpansion<TKernelAux>::TranslateFromFarField
for(index_t j = 0; j < total_num_coeffs; j++) {
const ArrayList <short int> &gamma_mapping = multiindex_mapping[j];
const ArrayList <short int> &lower_mappings_for_gamma =
const std::vector<short int>& gamma_mapping = multiindex_mapping[j];
const std::vector<short int>& lower_mappings_for_gamma =
lower_mapping_index[j];
double pos_coeff = 0;
double neg_coeff = 0;
for(index_t k = 0; k < lower_mappings_for_gamma.size(); k++) {
const ArrayList <short int> &inner_mapping =
const std::vector<short int>& inner_mapping =
multiindex_mapping[lower_mappings_for_gamma[k]];
int flag = 0;
@@ -987,22 +948,20 @@ template<typename TKernelAux>
void FarFieldExpansion<TKernelAux>::TranslateToLocal
(LocalExpansion<TKernelAux> &se, int truncation_order) {
Vector pos_arrtmp, neg_arrtmp;
Matrix derivative_map;
arma::vec pos_arrtmp, neg_arrtmp;
arma::mat derivative_map;
ka_->AllocateDerivativeMap(sea_->get_dimension(), 2 * truncation_order,
&derivative_map);
Vector local_center;
Vector cent_diff;
Vector local_coeffs;
derivative_map);
const arma::vec& local_center = se.get_center();
arma::vec cent_diff;
arma::vec& local_coeffs = se.get_coeffs();
int local_order = se.get_order();
int dimension = sea_->get_dimension();
int total_num_coeffs = sea_->get_total_num_coeffs(truncation_order);
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);
cent_diff.set_size(dimension);
// if the order of the far field expansion is greater than the
// local one we are adding onto, then increase the order.
@@ -1011,8 +970,8 @@ void FarFieldExpansion<TKernelAux>::TranslateToLocal
}
// Compute derivatives.
pos_arrtmp.Init(total_num_coeffs);
neg_arrtmp.Init(total_num_coeffs);
pos_arrtmp.set_size(total_num_coeffs);
neg_arrtmp.set_size(total_num_coeffs);
// Compute center difference divided by the bandwidth factor.
for(index_t j = 0; j < dimension; j++) {
@@ -1020,19 +979,19 @@ void FarFieldExpansion<TKernelAux>::TranslateToLocal
}
// Compute required partial derivatives.
ka_->ComputeDirectionalDerivatives(cent_diff, &derivative_map,
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map,
2 * truncation_order);
ArrayList<short int> beta_plus_alpha;
beta_plus_alpha.Init(dimension);
std::vector<short int> beta_plus_alpha;
beta_plus_alpha.reserve(dimension);
for(index_t j = 0; j < total_num_coeffs; j++) {
const ArrayList<short int> &beta_mapping = sea_->get_multiindex(j);
const std::vector<short 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<short int> &alpha_mapping = sea_->get_multiindex(k);
const std::vector<short 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];
}
@@ -1050,7 +1009,7 @@ void FarFieldExpansion<TKernelAux>::TranslateToLocal
} // end of k-loop
} // end of j-loop
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
arma::vec C_k_neg = sea_->get_neg_inv_multiindex_factorials();
for(index_t j = 0; j < total_num_coeffs; j++) {
local_coeffs[j] += (pos_arrtmp[j] + neg_arrtmp[j]) * C_k_neg[j];
}
fastlib/branches/fastlib-stl/mlpack/series_expansion/inverse_pow_dist_kernel_aux.h
+53 -59
View File
@@ -11,8 +11,8 @@ class InversePowDistGradientKernelAux {
private:
void SubFrom_(index_t dimension, int decrement,
const ArrayList<short int> &subtract_from,
ArrayList<short int> &result) const {
const std::vector<short int> &subtract_from,
std::vector<short int> &result) const {
for(index_t d = 0; d < subtract_from.size(); d++) {
if(d == dimension) {
@@ -30,8 +30,7 @@ class InversePowDistGradientKernelAux {
typedef SeriesExpansionAux TSeriesExpansionAux;
typedef FarFieldExpansion<InversePowDistGradientKernelAux> \
TFarFieldExpansion;
typedef FarFieldExpansion<InversePowDistGradientKernelAux> TFarFieldExpansion;
typedef LocalExpansion<InversePowDistGradientKernelAux> TLocalExpansion;
@@ -43,11 +42,6 @@ class InversePowDistGradientKernelAux {
*/
TSeriesExpansionAux sea_;
OT_DEF_BASIC(InversePowDistGradientKernelAux) {
OT_MY_OBJECT(kernel_);
OT_MY_OBJECT(sea_);
}
public:
void Init(double bandwidth, int max_order, int dim) {
@@ -56,40 +50,40 @@ class InversePowDistGradientKernelAux {
}
void AllocateDerivativeMap(int dim, int order,
Matrix *derivative_map) const {
derivative_map->Init(sea_.get_total_num_coeffs(order), 1);
arma::mat& derivative_map) const {
derivative_map.set_size(sea_.get_total_num_coeffs(order), 1);
}
void ComputeDirectionalDerivatives(const Vector &x,
Matrix *derivative_map, int order) const {
void ComputeDirectionalDerivatives(const arma::vec &x,
arma::mat& derivative_map, int order) const {
derivative_map->SetZero();
derivative_map.zeros();
// Squared L2 norm of the vector.
double squared_l2_norm = la::Dot(x, x);
double squared_l2_norm = dot(x, x);
// Temporary variable to look for arithmetic operations on
// multiindex.
ArrayList<short int> tmp_multiindex;
tmp_multiindex.Init(sea_.get_dimension());
std::vector<short int> tmp_multiindex;
tmp_multiindex.reserve(sea_.get_dimension());
for(index_t i = 0; i < derivative_map->n_rows(); i++) {
for(index_t i = 0; i < derivative_map.n_rows; i++) {
// Contribution to the current multiindex position.
double contribution = 0;
// Retrieve the multiindex mapping.
const ArrayList<short int> &multiindex = sea_.get_multiindex(i);
const std::vector<short int>& multiindex = sea_.get_multiindex(i);
// $D_{x}^{0} \phi_{\nu, d}(x)$ should be computed normally.
if(i == 0) {
derivative_map->set(0, 0, kernel_.EvalUnnorm(x.ptr()));
derivative_map(0, 0) = kernel_.EvalUnnorm(x.memptr());
continue;
}
// Compute the contribution of $D_{x}^{n - e_d} \phi_{\nu,
// d}(x)$ component for each $d$.
for(index_t d = 0; d < x.length(); d++) {
for(index_t d = 0; d < x.n_elem; d++) {
// Subtract 1 from the given dimension.
SubFrom_(d, 1, multiindex, tmp_multiindex);
@@ -102,7 +96,7 @@ class InversePowDistGradientKernelAux {
factor += (kernel_.lambda_ - 2);
}
contribution += factor *
derivative_map->get(n_minus_e_d_position, 0);
derivative_map(n_minus_e_d_position, 0);
}
// Subtract 2 from the given dimension.
@@ -118,39 +112,39 @@ class InversePowDistGradientKernelAux {
}
contribution += factor *
derivative_map->get(n_minus_two_e_d_position, 0);
derivative_map(n_minus_two_e_d_position, 0);
}
} // end of iterating over each dimension.
// Set the final contribution for this multiindex.
derivative_map->set(i, 0, -contribution / squared_l2_norm);
derivative_map(i, 0) = -contribution / squared_l2_norm;
} // end of iterating over all required multiindex positions...
// Iterate again, and invert the sum if the sum of the indices of
// the current mapping is odd.
for(index_t i = 1; i < derivative_map->n_rows(); i++) {
for(index_t i = 1; i < derivative_map.n_rows; i++) {
// Retrieve the multiindex mapping.
const ArrayList<short int> &multiindex = sea_.get_multiindex(i);
const std::vector<short int>& multiindex = sea_.get_multiindex(i);
// The sum of the indices.
index_t sum_of_indices = 0;
for(index_t d = 0; d < x.length(); d++) {
for(index_t d = 0; d < x.n_elem; d++) {
sum_of_indices += multiindex[d];
}
if(sum_of_indices % 2 == 1) {
derivative_map->set(i, 0, -derivative_map->get(i, 0));
derivative_map(i, 0) *= -1;
}
}
}
double ComputePartialDerivative(const Matrix &derivative_map,
const ArrayList<short int> &mapping) const {
double ComputePartialDerivative(const arma::mat& derivative_map,
const std::vector<short int>& mapping) const {
return derivative_map.get(sea_.ComputeMultiindexPosition(mapping), 0);
return derivative_map(sea_.ComputeMultiindexPosition(mapping), 0);
}
};
@@ -162,8 +156,8 @@ class InversePowDistKernelAux {
private:
void SubFrom_(index_t dimension, int decrement,
const ArrayList<short int> &subtract_from,
ArrayList<short int> &result) const {
const std::vector<short int>& subtract_from,
std::vector<short int>& result) const {
for(index_t d = 0; d < subtract_from.size(); d++) {
if(d == dimension) {
@@ -201,44 +195,44 @@ class InversePowDistKernelAux {
}
void AllocateDerivativeMap(int dim, int order,
Matrix *derivative_map) const {
derivative_map->Init(sea_.get_total_num_coeffs(order), 1);
arma::mat& derivative_map) const {
derivative_map.set_size(sea_.get_total_num_coeffs(order), 1);
}
void ComputeDirectionalDerivatives(const Vector &x,
Matrix *derivative_map, int order) const {
void ComputeDirectionalDerivatives(const arma::vec& x,
arma::mat& derivative_map, int order) const {
derivative_map->SetZero();
derivative_map.zeros();
// Squared L2 norm of the vector.
double squared_l2_norm = la::Dot(x, x);
double squared_l2_norm = dot(x, x);
// Temporary variable to look for arithmetic operations on
// multiindex.
ArrayList<short int> tmp_multiindex;
tmp_multiindex.Init(sea_.get_dimension());
std::vector<short int> tmp_multiindex;
tmp_multiindex.reserve(sea_.get_dimension());
// Get the inverse multiindex factorial factors.
const Vector &inv_multiindex_factorials =
const arma::vec& inv_multiindex_factorials =
sea_.get_inv_multiindex_factorials();
for(index_t i = 0; i < derivative_map->n_rows(); i++) {
for(index_t i = 0; i < derivative_map.n_rows; i++) {
// Contribution to the current multiindex position.
double contribution = 0;
// Retrieve the multiindex mapping.
const ArrayList<short int> &multiindex = sea_.get_multiindex(i);
const std::vector<short int>& multiindex = sea_.get_multiindex(i);
// $D_{x}^{0} \phi_{\nu, d}(x)$ should be computed normally.
if(i == 0) {
derivative_map->set(0, 0, kernel_.EvalUnnorm(x.ptr()));
derivative_map(0, 0) = kernel_.EvalUnnorm(x.memptr());
continue;
}
// The sum of the indices.
index_t sum_of_indices = 0;
for(index_t d = 0; d < x.length(); d++) {
for(index_t d = 0; d < x.n_elem; d++) {
sum_of_indices += multiindex[d];
}
@@ -250,7 +244,7 @@ class InversePowDistKernelAux {
// Compute the contribution of $D_{x}^{n - e_d} \phi_{\nu,
// d}(x)$ component for each $d$.
for(index_t d = 0; d < x.length(); d++) {
for(index_t d = 0; d < x.n_elem; d++) {
// Subtract 1 from the given dimension.
SubFrom_(d, 1, multiindex, tmp_multiindex);
@@ -258,7 +252,7 @@ class InversePowDistKernelAux {
sea_.ComputeMultiindexPosition(tmp_multiindex);
if(n_minus_e_d_position >= 0) {
contribution += first_factor * x[d] *
derivative_map->get(n_minus_e_d_position, 0) *
derivative_map(n_minus_e_d_position, 0) *
inv_multiindex_factorials[n_minus_e_d_position];
}
@@ -268,7 +262,7 @@ class InversePowDistKernelAux {
sea_.ComputeMultiindexPosition(tmp_multiindex);
if(n_minus_two_e_d_position >= 0) {
contribution += second_factor *
derivative_map->get(n_minus_two_e_d_position, 0) *
derivative_map(n_minus_two_e_d_position, 0) *
inv_multiindex_factorials[n_minus_two_e_d_position];
}
@@ -276,38 +270,38 @@ class InversePowDistKernelAux {
// Set the final contribution for this multiindex.
if(squared_l2_norm == 0) {
derivative_map->set(i, 0, 0);
derivative_map(i, 0) = 0;
}
else {
derivative_map->set(i, 0, -contribution / squared_l2_norm /
sum_of_indices / inv_multiindex_factorials[i]);
derivative_map(i, 0) = -contribution / squared_l2_norm /
sum_of_indices / inv_multiindex_factorials[i];
}
} // end of iterating over all required multiindex positions...
// Iterate again, and invert the sum if the sum of the indices of
// the current mapping is odd.
for(index_t i = 1; i < derivative_map->n_rows(); i++) {
for(index_t i = 1; i < derivative_map.n_rows; i++) {
// Retrieve the multiindex mapping.
const ArrayList<short int> &multiindex = sea_.get_multiindex(i);
const std::vector<short int>& multiindex = sea_.get_multiindex(i);
// The sum of the indices.
index_t sum_of_indices = 0;
for(index_t d = 0; d < x.length(); d++) {
for(index_t d = 0; d < x.n_elem; d++) {
sum_of_indices += multiindex[d];
}
if(sum_of_indices % 2 == 1) {
derivative_map->set(i, 0, -derivative_map->get(i, 0));
derivative_map(i, 0) *= -1;
}
}
}
double ComputePartialDerivative(const Matrix &derivative_map,
const ArrayList<short int> &mapping) const {
double ComputePartialDerivative(const arma::mat& derivative_map,
const std::vector<short int>& mapping) const {
return derivative_map.get(sea_.ComputeMultiindexPosition(mapping), 0);
return derivative_map(sea_.ComputeMultiindexPosition(mapping), 0);
}
};
fastlib/branches/fastlib-stl/mlpack/series_expansion/kernel_aux.h
+37 -55
View File
@@ -43,11 +43,6 @@ class GaussianKernelMultAux {
/** pointer to the series expansion auxiliary object */
TSeriesExpansionAux sea_;
OT_DEF_BASIC(GaussianKernelMultAux) {
OT_MY_OBJECT(kernel_);
OT_MY_OBJECT(sea_);
}
public:
void Init(double bandwidth, int max_order, int dim) {
@@ -60,14 +55,14 @@ class GaussianKernelMultAux {
}
void AllocateDerivativeMap(int dim, int order,
Matrix *derivative_map) const {
derivative_map->Init(dim, order + 1);
arma::mat& derivative_map) const {
derivative_map.set_size(dim, order + 1);
}
void ComputeDirectionalDerivatives(const Vector &x,
Matrix *derivative_map, int order) const {
void ComputeDirectionalDerivatives(const arma::vec &x,
arma::mat& derivative_map, int order) const {
int dim = x.length();
int dim = x.n_elem;
// precompute necessary Hermite polynomials based on coordinate difference
for(index_t d = 0; d < dim; d++) {
@@ -76,30 +71,30 @@ class GaussianKernelMultAux {
double d2 = 2 * coord_div_band;
double facj = exp(-coord_div_band * coord_div_band);
derivative_map->set(d, 0, facj);
derivative_map(d, 0) = facj;
if(order > 0) {
derivative_map->set(d, 1, d2 * facj);
derivative_map(d, 1) = d2 * facj;
if(order > 1) {
for(index_t k = 1; k < order; k++) {
int k2 = k * 2;
derivative_map->set(d, k + 1, d2 * derivative_map->get(d, k) -
k2 * derivative_map->get(d, k - 1));
derivative_map(d, k + 1) = d2 * derivative_map(d, k) -
k2 * derivative_map(d, k - 1);
}
}
}
} // end of looping over each dimension
}
double ComputePartialDerivative(const Matrix &derivative_map,
const ArrayList<short int> &mapping) const {
double ComputePartialDerivative(const arma::mat& derivative_map,
const std::vector<short int>& mapping) const {
double partial_derivative = 1.0;
for(index_t d = 0; d < mapping.size(); d++) {
partial_derivative *= derivative_map.get(d, mapping[d]);
partial_derivative *= derivative_map(d, mapping[d]);
}
return partial_derivative;
}
@@ -299,11 +294,6 @@ class GaussianKernelAux {
/** pointer to the series expansion auxiliary object */
TSeriesExpansionAux sea_;
OT_DEF_BASIC(GaussianKernelAux) {
OT_MY_OBJECT(kernel_);
OT_MY_OBJECT(sea_);
}
public:
void Init(double bandwidth, int max_order, int dim) {
@@ -316,14 +306,14 @@ class GaussianKernelAux {
}
void AllocateDerivativeMap(int dim, int order,
Matrix *derivative_map) const {
derivative_map->Init(dim, order + 1);
arma::mat& derivative_map) const {
derivative_map.set_size(dim, order + 1);
}
void ComputeDirectionalDerivatives(const Vector &x,
Matrix *derivative_map, int order) const {
void ComputeDirectionalDerivatives(const arma::vec& x,
arma::mat& derivative_map, int order) const {
int dim = x.length();
int dim = x.n_elem;
// precompute necessary Hermite polynomials based on coordinate difference
for(index_t d = 0; d < dim; d++) {
@@ -332,48 +322,47 @@ class GaussianKernelAux {
double d2 = 2 * coord_div_band;
double facj = exp(-coord_div_band * coord_div_band);
derivative_map->set(d, 0, facj);
derivative_map(d, 0) = facj;
if(order > 0) {
derivative_map->set(d, 1, d2 * facj);
derivative_map(d, 1) = d2 * facj;
if(order > 1) {
for(index_t k = 1; k < order; k++) {
int k2 = k * 2;
derivative_map->set(d, k + 1, d2 * derivative_map->get(d, k) -
k2 * derivative_map->get(d, k - 1));
derivative_map(d, k + 1) = d2 * derivative_map(d, k) -
k2 * derivative_map(d, k - 1);
}
}
}
} // end of looping over each dimension
}
double ComputePartialDerivative(const Matrix &derivative_map,
const ArrayList<short int> &mapping) const {
double ComputePartialDerivative(const arma::mat& derivative_map,
const std::vector<short int>& mapping) const {
double partial_derivative = 1.0;
for(index_t d = 0; d < mapping.size(); d++) {
partial_derivative *= derivative_map.get(d, mapping[d]);
partial_derivative *= derivative_map(d, mapping[d]);
}
return partial_derivative;
}
template<typename TBound>
int OrderForConvolvingFarField(const TBound &far_field_region,
const Vector &far_field_region_centroid,
const arma::vec &far_field_region_centroid,
const TBound &local_field_region,
const Vector &local_field_region_centroid,
const arma::vec &local_field_region_centroid,
double min_dist_sqd_regions,
double max_dist_sqd_regions,
double max_error,
double *actual_error) const {
double squared_distance_between_two_centroids =
la::DistanceSqEuclidean(far_field_region_centroid.length(),
far_field_region_centroid.ptr(),
local_field_region_centroid.ptr());
la::DistanceSqEuclidean(far_field_region_centroid,
local_field_region_centroid);
double frontfactor =
exp(-squared_distance_between_two_centroids /
(4 * kernel_.bandwidth_sq()));
@@ -586,12 +575,6 @@ class EpanKernelAux {
InversePowDistKernelAux squared_component_;
OT_DEF_BASIC(EpanKernelAux) {
OT_MY_OBJECT(kernel_);
OT_MY_OBJECT(sea_);
OT_MY_OBJECT(squared_component_);
}
public:
void Init(double bandwidth, int max_order, int dim) {
@@ -607,26 +590,25 @@ class EpanKernelAux {
}
void AllocateDerivativeMap(int dim, int order,
Matrix *derivative_map) const {
derivative_map->Init(sea_.get_total_num_coeffs(order), 1);
arma::mat& derivative_map) const {
derivative_map.set_size(sea_.get_total_num_coeffs(order), 1);
}
void ComputeDirectionalDerivatives(const Vector &x,
Matrix *derivative_map, int order) const {
void ComputeDirectionalDerivatives(const arma::vec& x,
arma::mat& derivative_map, int order) const {
// Compute the derivatives for $||x||^2$ and negate it. Then, add
// $(1, 0, 0, ... 0)$ to it.
squared_component_.ComputeDirectionalDerivatives(x, derivative_map, order);
la::Scale(derivative_map->n_rows(), -1, derivative_map->GetColumnPtr(0));
(derivative_map->GetColumnPtr(0))[0] += 1.0;
derivative_map.unsafe_col(0) *= -1;
derivative_map(0, 0) += 1.0;
}
double ComputePartialDerivative(const Matrix &derivative_map,
const ArrayList<short int> &mapping) const {
double ComputePartialDerivative(const arma::mat& derivative_map,
const std::vector<short int>& mapping) const {
return derivative_map.get(sea_.ComputeMultiindexPosition(mapping), 0);
return derivative_map(sea_.ComputeMultiindexPosition(mapping), 0);
}
template<typename TBound>
fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion.h
+11 -17
View File
@@ -12,7 +12,7 @@
#ifndef LOCAL_EXPANSION
#define LOCAL_EXPANSION
#include "fastlib/fastlib.h"
#include <fastlib/fastlib.h>
#include "kernel_aux.h"
#include "series_expansion_aux.h"
@@ -28,10 +28,10 @@ class LocalExpansion {
private:
/** The center of the expansion */
Vector center_;
arma::vec center_;
/** The coefficients */
Vector coeffs_;
arma::vec coeffs_;
/** order */
int order_;
@@ -45,12 +45,6 @@ class LocalExpansion {
/** pointer to the precomputed constants inside kernel auxiliary object */
const typename TKernelAux::TSeriesExpansionAux *sea_;
OT_DEF(LocalExpansion) {
OT_MY_OBJECT(center_);
OT_MY_OBJECT(coeffs_);
OT_MY_OBJECT(order_);
}
public:
// getters and setters
@@ -59,12 +53,12 @@ class LocalExpansion {
double bandwidth_sq() const { return kernel_->bandwidth_sq(); }
/** Get the center of expansion */
Vector* get_center() { return &center_; }
const Vector* get_center() const { return &center_; }
arma::vec& get_center() { return center_; }
const arma::vec& get_center() const { return center_; }
/** Get the coefficients */
const Vector& get_coeffs() const { return coeffs_; }
arma::vec& get_coeffs() { return coeffs_; }
const arma::vec& get_coeffs() const { return coeffs_; }
/** Get the approximation order */
int get_order() const { return order_; }
@@ -81,26 +75,26 @@ class LocalExpansion {
* Accumulates the local moment represented by the given reference
* data into the coefficients
*/
void AccumulateCoeffs(const Matrix& data, const Vector& weights,
void AccumulateCoeffs(const arma::mat& data, const arma::vec& weights,
int begin, int end, int order);
/**
* This does not apply for local coefficients.
*/
void RefineCoeffs(const Matrix& data, const Vector& weights,
void RefineCoeffs(const arma::mat& data, const arma::vec& weights,
int begin, int end, int order) { }
/**
* Evaluates the local coefficients at the given point
*/
double EvaluateField(const Matrix& data, int row_num) const;
double EvaluateField(const arma::mat& data, int row_num) const;
double EvaluateField(const double *x_q) const;
/**
* Initializes the current local expansion object with the given
* center.
*/
void Init(const Vector& center, const TKernelAux &ka);
void Init(const arma::vec& center, const TKernelAux &ka);
void Init(const TKernelAux &ka);
/**
fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion_impl.h
+55 -74
View File
@@ -7,32 +7,27 @@
template<typename TKernelAux>
void LocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
const Vector& weights,
void LocalExpansion<TKernelAux>::AccumulateCoeffs(const arma::mat& data,
const arma::vec& weights,
int begin, int end,
int order) {
if(order > 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());
// TODO: this is supposed to be an alias (just get it compiling for now)
arma::vec neg_inv_multiindex_factorials = sea_->get_neg_inv_multiindex_factorials();
// declare deritave mapping
Matrix derivative_map;
ka_->AllocateDerivativeMap(dim, order, &derivative_map);
arma::mat derivative_map;
ka_->AllocateDerivativeMap(dim, order, derivative_map);
// some temporary variables
Vector arrtmp;
arrtmp.Init(total_num_coeffs);
Vector x_r_minus_x_Q;
x_r_minus_x_Q.Init(dim);
arma::vec arrtmp(total_num_coeffs), x_r_minus_x_Q(dim);
// The bandwidth factor to be divided along each dimension.
double bandwidth_factor = ka_->BandwidthFactor(kernel_->bandwidth_sq());
@@ -42,22 +37,22 @@ void LocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
// 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)) /
x_r_minus_x_Q[d] = (center_[d] - data(d, r)) /
bandwidth_factor;
}
// precompute necessary partial derivatives based on coordinate difference
ka_->ComputeDirectionalDerivatives(x_r_minus_x_Q, &derivative_map, order);
ka_->ComputeDirectionalDerivatives(x_r_minus_x_Q, derivative_map, order);
// compute h_{beta}((x_r - x_Q) / sqrt(2h^2))
for(index_t j = 0; j < total_num_coeffs; j++) {
const ArrayList<short int> &mapping = sea_->get_multiindex(j);
const std::vector<short 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];
arrtmp[j];
}
} // End of looping through each reference point.
}
@@ -73,7 +68,7 @@ void LocalExpansion<TKernelAux>::PrintDebug(const char *name,
fprintf(stream, "Local expansion\n");
fprintf(stream, "Center: ");
for (index_t i = 0; i < center_.length(); i++) {
for (index_t i = 0; i < center_.n_elem; i++) {
fprintf(stream, "%g ", center_[i]);
}
fprintf(stream, "\n");
@@ -87,7 +82,7 @@ void LocalExpansion<TKernelAux>::PrintDebug(const char *name,
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<short int> &mapping = sea_->get_multiindex(i);
const std::vector<short int>& mapping = sea_->get_multiindex(i);
fprintf(stream, "%g", coeffs_[i]);
for(index_t d = 0; d < dim; d++) {
@@ -102,18 +97,17 @@ void LocalExpansion<TKernelAux>::PrintDebug(const char *name,
}
template<typename TKernelAux>
double LocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
double LocalExpansion<TKernelAux>::EvaluateField(const arma::mat& data,
int row_num) const {
return EvaluateField(data.GetColumnPtr(row_num));
return EvaluateField(data.colptr(row_num));
}
template<typename TKernelAux>
double LocalExpansion<TKernelAux>::EvaluateField(const double *x_q) const {
// if there are no local expansion here, then return 0
if(order_ < 0) {
if(order_ < 0)
return 0;
}
index_t k, t, tail;
@@ -130,12 +124,9 @@ double LocalExpansion<TKernelAux>::EvaluateField(const double *x_q) const {
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<short int> heads;
heads.Init(dim + 1);
arma::vec x_Q_to_x_q(dim), tmp(total_num_coeffs);
std::vector<short int> heads;
heads.reserve(dim + 1);
// compute (x_q - x_Q) / (sqrt(2h^2))
for(index_t i = 0; i < dim; i++) {
@@ -169,19 +160,18 @@ double LocalExpansion<TKernelAux>::EvaluateField(const double *x_q) const {
}
template<typename TKernelAux>
void LocalExpansion<TKernelAux>::Init(const Vector& center,
void LocalExpansion<TKernelAux>::Init(const arma::vec& center,
const TKernelAux &ka) {
// copy kernel type, center, and bandwidth squared
kernel_ = &(ka.kernel_);
center_.Copy(center);
center_ = center;
order_ = -1;
sea_ = &(ka.sea_);
ka_ = &ka;
// initialize coefficient array
coeffs_.Init(sea_->get_max_total_num_coeffs());
coeffs_.SetZero();
coeffs_.zeros(sea_->get_max_total_num_coeffs());
}
template<typename TKernelAux>
@@ -191,12 +181,11 @@ void LocalExpansion<TKernelAux>::Init(const TKernelAux &ka) {
kernel_ = &(ka.kernel_);
order_ = -1;
sea_ = &(ka.sea_);
center_.Init(sea_->get_dimension());
center_.set_size(sea_->get_dimension());
ka_ = &ka;
// initialize coefficient array
coeffs_.Init(sea_->get_max_total_num_coeffs());
coeffs_.SetZero();
coeffs_.zeros(sea_->get_max_total_num_coeffs());
}
template<typename TKernelAux>
@@ -216,13 +205,12 @@ 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;
arma::vec pos_arrtmp, neg_arrtmp;
arma::mat derivative_map;
arma::vec cent_diff;
int dimension = sea_->get_dimension();
ka_->AllocateDerivativeMap(dimension, 2 * order_, &derivative_map);
ka_->AllocateDerivativeMap(dimension, 2 * order_, derivative_map);
int far_order = se.get_order();
int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
@@ -230,9 +218,9 @@ void LocalExpansion<TKernelAux>::TranslateFromFarField
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);
arma::vec& far_center = se.get_center();
arma::vec& far_coeffs = se.get_coeffs();
cent_diff.set_size(dimension);
// if the order of the far field expansion is greater than the
// local one we are adding onto, then increase the order.
@@ -241,8 +229,8 @@ void LocalExpansion<TKernelAux>::TranslateFromFarField
}
// compute Gaussian derivative
pos_arrtmp.Init(total_num_coeffs);
neg_arrtmp.Init(total_num_coeffs);
pos_arrtmp.set_size(total_num_coeffs);
neg_arrtmp.set_size(total_num_coeffs);
// compute center difference divided by bw_times_sqrt_two;
for(index_t j = 0; j < dimension; j++) {
@@ -250,18 +238,18 @@ void LocalExpansion<TKernelAux>::TranslateFromFarField
}
// compute required partial derivatives
ka_->ComputeDirectionalDerivatives(cent_diff, &derivative_map, 2 * order_);
ArrayList<short int> beta_plus_alpha;
beta_plus_alpha.Init(dimension);
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map, 2 * order_);
std::vector<short int> beta_plus_alpha;
beta_plus_alpha.reserve(dimension);
for(index_t j = 0; j < total_num_coeffs; j++) {
const ArrayList<short int> &beta_mapping = sea_->get_multiindex(j);
const std::vector<short 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<short int> &alpha_mapping = sea_->get_multiindex(k);
const std::vector<short 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];
}
@@ -270,16 +258,14 @@ void LocalExpansion<TKernelAux>::TranslateFromFarField
double prod = far_coeffs[k] * derivative_factor;
if(prod > 0) {
if(prod > 0)
pos_arrtmp[j] += prod;
}
else {
else
neg_arrtmp[j] += prod;
}
} // end of k-loop
} // end of j-loop
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
arma::vec& 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];
}
@@ -296,28 +282,25 @@ void LocalExpansion<TKernelAux>::TranslateToLocal(LocalExpansion &se) {
// get the center and the order and the total number of coefficients of
// the expansion we are translating from. Also get coefficients we
// are translating
Vector new_center;
new_center.Alias(*(se.get_center()));
const arma::vec& new_center = se.get_center();
int prev_order = se.get_order();
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
const ArrayList<short int> *upper_mapping_index =
const std::vector<std::vector<short int> >& upper_mapping_index =
sea_->get_upper_mapping_index();
Vector new_coeffs;
new_coeffs.Alias(se.get_coeffs());
arma::vec& new_coeffs = se.get_coeffs();
// dimension
int dim = sea_->get_dimension();
// temporary variable
ArrayList<short int> tmp_storage;
tmp_storage.Init(dim);
std::vector<short int> tmp_storage;
tmp_storage.reserve(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);
arma::vec center_diff(dim);
for(index_t d = 0; d < dim; d++) {
center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
}
@@ -329,25 +312,24 @@ void LocalExpansion<TKernelAux>::TranslateToLocal(LocalExpansion &se) {
}
// inverse multiindex factorials
Vector C_k;
C_k.Alias(sea_->get_inv_multiindex_factorials());
const arma::vec& C_k = sea_->get_inv_multiindex_factorials();
// do the actual translation
for(index_t j = 0; j < total_num_coeffs; j++) {
const ArrayList<short int> &alpha_mapping = sea_->get_multiindex(j);
const ArrayList<short int> &upper_mappings_for_alpha =
upper_mapping_index[j];
const std::vector<short int>& alpha_mapping = sea_->get_multiindex(j);
const std::vector<short 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<short int> &beta_mapping =
const std::vector<short int>& beta_mapping =
sea_->get_multiindex(upper_mappings_for_alpha[k]);
int flag = 0;
double diff1 = 1.0;
@@ -385,4 +367,3 @@ void LocalExpansion<TKernelAux>::TranslateToLocal(LocalExpansion &se) {
}
#endif
fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion.h
+26 -32
View File
@@ -12,8 +12,7 @@
#ifndef MULT_FARFIELD_EXPANSION
#define MULT_FARFIELD_EXPANSION
#include "fastlib/fastlib.h"
#include <fastlib/fastlib.h>
#include "kernel_aux.h"
#include "mult_series_expansion_aux.h"
@@ -29,10 +28,10 @@ class MultFarFieldExpansion {
private:
/** @brief The center of the expansion. */
Vector center_;
arma::vec center_;
/** @brief The coefficients. */
Vector coeffs_;
arma::vec coeffs_;
/** @brief The order of approximation. */
int order_;
@@ -46,12 +45,6 @@ class MultFarFieldExpansion {
/** pointer to the precomputed constants inside kernel auxiliary object */
const typename TKernelAux::TSeriesExpansionAux *sea_;
OT_DEF(MultFarFieldExpansion) {
OT_MY_OBJECT(center_);
OT_MY_OBJECT(coeffs_);
OT_MY_OBJECT(order_);
}
public:
// getters and setters
@@ -60,12 +53,13 @@ class MultFarFieldExpansion {
double bandwidth_sq() const { return kernel_->bandwidth_sq(); }
/** Get the center of expansion */
Vector *get_center() { return &center_; }
arma::vec& get_center() { return center_; }
const Vector *get_center() const { return &center_; }
const arma::vec& get_center() const { return center_; }
/** Get the coefficients */
const Vector& get_coeffs() const { return coeffs_; }
arma::vec& get_coeffs() { return coeffs_; }
const arma::vec& get_coeffs() const { return coeffs_; }
/** Get the approximation order */
int get_order() const { return order_; }
@@ -84,9 +78,9 @@ class MultFarFieldExpansion {
* Set the center of the expansion - assumes that the center has been
* initialized before...
*/
void set_center(const Vector &center) {
void set_center(const arma::vec& center) {
for(index_t i = 0; i < center.length(); i++) {
for(index_t i = 0; i < center.n_elem; i++) {
center_[i] = center[i];
}
}
@@ -97,56 +91,56 @@ class MultFarFieldExpansion {
* Accumulates the far field moment represented by the given reference
* data into the coefficients
*/
void AccumulateCoeffs(const Matrix& data, const Vector& weights,
void AccumulateCoeffs(const arma::mat& data, const arma::vec& weights,
int begin, int end, int order);
/**
* Refine the far field moment that has been computed before up to
* a new order.
*/
void RefineCoeffs(const Matrix& data, const Vector& weights,
void RefineCoeffs(const arma::mat& data, const arma::vec& weights,
int begin, int end, int order);
/**
* Evaluates the far-field coefficients at the given point
*/
double EvaluateField(const Matrix& data, int row_num, int order) const;
double EvaluateField(const arma::mat& data, int row_num, int order) const;
double EvaluateField(const double *x_q, int order) const;
/**
* Evaluates the two-way convolution mixed with exhaustive computations
* with two other far field expansions
*/
double MixField(const Matrix &data, int node1_begin, int node1_end,
double MixField(const arma::mat& data, int node1_begin, int node1_end,
int node2_begin, int node2_end,
const MultFarFieldExpansion &fe2,
const MultFarFieldExpansion &fe3,
const MultFarFieldExpansion& fe2,
const MultFarFieldExpansion& fe3,
int order2, int order3) const;
double ConvolveField(const MultFarFieldExpansion &fe, int order) const;
double ConvolveField(const MultFarFieldExpansion& fe, int order) const;
/**
* Evaluates the three-way convolution with two other far field
* expansions
*/
double ConvolveField(const MultFarFieldExpansion &fe2,
const MultFarFieldExpansion &fe3,
double ConvolveField(const MultFarFieldExpansion& fe2,
const MultFarFieldExpansion& fe3,
int order1, int order2, int order3) const;
/**
* Initializes the current far field expansion object with the given
* center.
*/
void Init(const Vector& center, const TKernelAux &ka);
void Init(const TKernelAux &ka);
void Init(const arma::vec& center, const TKernelAux& ka);
void Init(const TKernelAux& ka);
/** @brief Computes the required order for evaluating the far field
* expansion for any query point within the specified region
* for a given bound.
*/
template<typename TBound>
int OrderForEvaluating(const TBound &far_field_region,
const TBound &local_field_region,
int 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;
@@ -161,8 +155,8 @@ class MultFarFieldExpansion {
* -1 if approximation up to the maximum order is not possible
*/
template<typename TBound>
int OrderForConvertingToLocal(const TBound &far_field_region,
const TBound &local_field_region,
int OrderForConvertingToLocal(const TBound& far_field_region,
const TBound& local_field_region,
double min_dist_sqd_regions,
double max_dist_sqd_regions,
double required_bound,
@@ -177,13 +171,13 @@ class MultFarFieldExpansion {
* here. The translated coefficients are added up to the
* ones here.
*/
void TranslateFromFarField(const MultFarFieldExpansion &se);
void TranslateFromFarField(const MultFarFieldExpansion& se);
/**
* Translate to the given local expansion. The translated coefficients
* are added up to the passed-in local expansion coefficients.
*/
void TranslateToLocal(MultLocalExpansion<TKernelAux> &se,
void TranslateToLocal(MultLocalExpansion<TKernelAux>& se,
int truncation_order);
};
fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion_impl.h
+80 -100
View File
@@ -6,26 +6,24 @@
#define MULT_FARFIELD_EXPANSION_IMPL_H
template<typename TKernelAux>
void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
const Vector& weights,
void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const arma::mat& data,
const arma::vec& weights,
int begin, int end,
int order) {
int dim = data.n_rows();
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;
arma::vec 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();
x_r.set_size(dim);
tmp.set_size(max_total_num_coeffs);
arma::vec pos_coeffs;
arma::vec neg_coeffs;
pos_coeffs.zeros(max_total_num_coeffs);
neg_coeffs.zeros(max_total_num_coeffs);
// set to new order if greater
if(order_ < order) {
@@ -33,7 +31,7 @@ void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
}
// get the order of traversal for the given order of approximation
const ArrayList<short int> &traversal_order =
const std::vector<short int>& traversal_order =
sea_->traversal_mapping_[order_];
// Repeat for each reference point in this reference node.
@@ -42,16 +40,16 @@ void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
// 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;
x_r[i] = (data(i, r) - center_[i]) / bandwidth_factor;
}
tmp.SetZero();
tmp.zeros();
tmp[0] = 1.0;
for(index_t i = 1; i < total_num_coeffs; i++) {
int index = traversal_order[i];
const ArrayList<short int> &lower_mappings =
const std::vector<short int>& lower_mappings =
sea_->lower_mapping_index_[index];
// from the direct descendant, recursively compute the multipole moments
@@ -59,8 +57,8 @@ void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
lower_mappings[lower_mappings.size() - 2];
int position = 0;
const ArrayList<short int> &mapping = sea_->multiindex_mapping_[index];
const ArrayList<short int> &direct_ancestor_mapping =
const std::vector<short int>& mapping = sea_->multiindex_mapping_[index];
const std::vector<short 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]) {
@@ -97,20 +95,19 @@ void MultFarFieldExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
template<typename TKernelAux>
double MultFarFieldExpansion<TKernelAux>::ConvolveField
(const MultFarFieldExpansion &fe, int order) const {
(const MultFarFieldExpansion& fe, int order) const {
// The bandwidth factor and the multiindex mapping stuffs.
double bandwidth_factor = ka_->BandwidthFactor(bandwidth_sq());
const ArrayList<short int> *multiindex_mapping =
const std::vector<short int>& multiindex_mapping =
sea_->get_multiindex_mapping();
const ArrayList<short int> *lower_mapping_index =
const std::vector<short int>& lower_mapping_index =
sea_->get_lower_mapping_index();
// Get the total number of coefficients and the coefficient themselves.
int total_num_coeffs = sea_->get_total_num_coeffs(order);
int dim = sea_->get_dimension();
Vector coeffs2;
coeffs2.Alias(fe.get_coeffs());
arma::vec& coeffs2 = fe.get_coeffs();
// Actual accumulated sum.
double neg_sum = 0;
@@ -118,41 +115,38 @@ double MultFarFieldExpansion<TKernelAux>::ConvolveField
double sum = 0;
// The partial derivatives table.
Matrix derivative_map_alpha;
ka_->AllocateDerivativeMap(dim, order, &derivative_map_alpha);
arma::mat derivative_map_alpha;
ka_->AllocateDerivativeMap(dim, order, derivative_map_alpha);
// Compute the center difference and its table of partial
// derivatives.
Vector xI_xJ;
xI_xJ.Init(dim);
Vector xJ_center;
xJ_center.Alias(*(fe.get_center()));
arma::vec xI_xJ(dim);
arma::vec& xJ_center = fe.get_center();
for(index_t d = 0; d < dim; d++) {
xI_xJ[d] = (center_[d] - xJ_center[d]) / bandwidth_factor;
}
ka_->ComputeDirectionalDerivatives(xI_xJ, &derivative_map_alpha, order);
ka_->ComputeDirectionalDerivatives(xI_xJ, derivative_map_alpha, order);
// The inverse factorials.
Vector inv_multiindex_factorials;
inv_multiindex_factorials.Alias(sea_->get_inv_multiindex_factorials());
arma::vec& inv_multiindex_factorials = sea_->get_inv_multiindex_factorials();
// The temporary space for computing the difference of two mappings.
ArrayList<short int> alpha_minus_beta_mapping;
alpha_minus_beta_mapping.Init(dim);
std::vector<short int> alpha_minus_beta_mapping;
alpha_minus_beta_mapping.reserve(dim);
// The main loop.
for(index_t alpha = 0; alpha < total_num_coeffs; alpha++) {
const ArrayList<short int> &alpha_mapping = multiindex_mapping[alpha];
const ArrayList<short int> &lower_mappings_for_alpha =
const std::vector<short int>& alpha_mapping = multiindex_mapping[alpha];
const std::vector<short int>& lower_mappings_for_alpha =
lower_mapping_index[alpha];
double alpha_derivative = ka_->ComputePartialDerivative
(derivative_map_alpha, alpha_mapping);
for(index_t beta = 0; beta < lower_mappings_for_alpha.size(); beta++) {
const ArrayList<short int> &beta_mapping =
const std::vector<short int>& beta_mapping =
multiindex_mapping[lower_mappings_for_alpha[beta]];
double n_choose_k_factor = sea_->get_n_multichoose_k_by_pos
@@ -192,8 +186,8 @@ double MultFarFieldExpansion<TKernelAux>::ConvolveField
}
template<typename TKernelAux>
void MultFarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
const Vector& weights,
void MultFarFieldExpansion<TKernelAux>::RefineCoeffs(const arma::mat& data,
const arma::vec& weights,
int begin, int end,
int order) {
@@ -208,16 +202,17 @@ void MultFarFieldExpansion<TKernelAux>::RefineCoeffs(const Matrix& data,
else {
order_ = order;
coeffs_.SetZero();
coeffs_.zeros();
AccumulateCoeffs(data, weights, begin, end, order);
}
}
template<typename TKernelAux>
double MultFarFieldExpansion<TKernelAux>::EvaluateField(const Matrix& data,
double MultFarFieldExpansion<TKernelAux>::EvaluateField(const arma::mat& data,
int row_num,
int order) const {
return EvaluateField(data.GetColumnPtr(row_num), order);
// TODO: use of "row" is misleading and wrong
return EvaluateField(data.unsafe_col(row_num).memptr(), order);
}
template<typename TKernelAux>
@@ -239,16 +234,14 @@ double MultFarFieldExpansion<TKernelAux>::EvaluateField(const double *x_q,
double multipole_sum = 0;
// computed derivative map
Matrix derivative_map;
ka_->AllocateDerivativeMap(dim, order_, &derivative_map);
arma::mat derivative_map;
ka_->AllocateDerivativeMap(dim, order_, derivative_map);
// temporary variable
Vector arrtmp;
arrtmp.Init(total_num_coeffs);
arma::vec arrtmp(total_num_coeffs);
// (x_q - x_R) scaled by bandwidth
Vector x_q_minus_x_R;
x_q_minus_x_R.Init(dim);
arma::vec x_q_minus_x_R(dim);
// compute (x_q - x_R) / (sqrt(2h^2))
for(index_t d = 0; d < dim; d++) {
@@ -256,17 +249,17 @@ double MultFarFieldExpansion<TKernelAux>::EvaluateField(const double *x_q,
}
// compute deriative maps based on coordinate difference.
ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, &derivative_map, order_);
ka_->ComputeDirectionalDerivatives(x_q_minus_x_R, derivative_map, order_);
// get the order of traversal for the given order of approximation
const ArrayList<short int> &traversal_order =
const std::vector<short 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<short int> &mapping = sea_->get_multiindex(index);
const std::vector<short int>& mapping = sea_->get_multiindex(index);
double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping);
double prod = coeffs_[index] * arrtmp;
@@ -283,19 +276,18 @@ double MultFarFieldExpansion<TKernelAux>::EvaluateField(const double *x_q,
}
template<typename TKernelAux>
void MultFarFieldExpansion<TKernelAux>::Init(const Vector& center,
void MultFarFieldExpansion<TKernelAux>::Init(const arma::vec& center,
const TKernelAux &ka) {
// copy kernel type, center, and bandwidth squared
kernel_ = &(ka.kernel_);
center_.Copy(center);
center_ = center;
order_ = -1;
sea_ = &(ka.sea_);
ka_ = &ka;
// Initialize coefficient array
coeffs_.Init(sea_->get_max_total_num_coeffs());
coeffs_.SetZero();
coeffs_.zeros(sea_->get_max_total_num_coeffs());
}
template<typename TKernelAux>
@@ -305,13 +297,11 @@ template<typename TKernelAux>
kernel_ = &(ka.kernel_);
order_ = -1;
sea_ = &(ka.sea_);
center_.Init(sea_->get_dimension());
center_.SetZero();
center_.zeros(sea_->get_dimension());
ka_ = &ka;
// Initialize coefficient array.
coeffs_.Init(sea_->get_max_total_num_coeffs());
coeffs_.SetZero();
coeffs_.zeros(sea_->get_max_total_num_coeffs());
}
@@ -355,7 +345,7 @@ void MultFarFieldExpansion<TKernelAux>::PrintDebug
fprintf(stream, "Far field expansion\n");
fprintf(stream, "Center: ");
for (index_t i = 0; i < center_.length(); i++) {
for (index_t i = 0; i < center_.n_elem; i++) {
fprintf(stream, "%g ", center_[i]);
}
fprintf(stream, "\n");
@@ -369,7 +359,7 @@ void MultFarFieldExpansion<TKernelAux>::PrintDebug
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<short int> &mapping = sea_->get_multiindex(i);
const std::vector<short int>& mapping = sea_->get_multiindex(i);
fprintf(stream, "%g ", coeffs_[i]);
fprintf(stream, "(-1)^(");
@@ -401,24 +391,19 @@ void MultFarFieldExpansion<TKernelAux>::TranslateFromFarField
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<short int> *multiindex_mapping =
const std::vector<std::vector<short int> >& multiindex_mapping =
sea_->get_multiindex_mapping();
const ArrayList<short int> *lower_mapping_index =
const std::vector<std::vector<short int> >& lower_mapping_index =
sea_->get_lower_mapping_index();
ArrayList<short int> tmp_storage;
Vector center_diff;
Vector inv_multiindex_factorials;
center_diff.Init(dim);
std::vector<short int> tmp_storage;
arma::vec center_diff(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());
const arma::vec& prev_coeffs = se.get_coeffs();
const arma::vec& prev_center = se.get_center();
tmp_storage.reserve(sea_->get_dimension());
const arma::vec& inv_multiindex_factorials = sea_->get_inv_multiindex_factorials();
// no coefficients can be translated
if(order == -1) {
@@ -434,22 +419,21 @@ void MultFarFieldExpansion<TKernelAux>::TranslateFromFarField
}
// get the order of traversal for the given order of approximation
const ArrayList<short int> &traversal_order =
const std::vector<short int>& traversal_order =
sea_->traversal_mapping_[order];
for(index_t j = 0; j < total_num_coeffs; j++) {
int index = traversal_order[j];
const ArrayList<short int> &gamma_mapping = multiindex_mapping[index];
const ArrayList<short int> &lower_mappings_for_gamma =
const std::vector<short int>& gamma_mapping = multiindex_mapping[index];
const std::vector<short 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<short int> &inner_mapping =
multiindex_mapping[lower_mappings_for_gamma[k]];
const std::vector<short int>& inner_mapping = multiindex_mapping[lower_mappings_for_gamma[k]];
int flag = 0;
double diff1;
@@ -496,22 +480,18 @@ 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();
arma::vec pos_arrtmp, neg_arrtmp;
arma::mat derivative_map;
// get center and coefficients for local expansion
const arma::vec& local_center = se.get_center();
int dimension = sea_->get_dimension();
arma::vec cent_diff(dimension);
arma::vec& local_coeffs = se.get_coeffs();
int local_order = se.get_order();
int total_num_coeffs = sea_->get_total_num_coeffs(truncation_order);
double bandwidth_factor = ka_->BandwidthFactor(se.bandwidth_sq());
ka_->AllocateDerivativeMap(dimension, 2 * truncation_order, &derivative_map);
// get center and coefficients for local expansion
local_center.Alias(*(se.get_center()));
local_coeffs.Alias(se.get_coeffs());
cent_diff.Init(dimension);
ka_->AllocateDerivativeMap(dimension, 2 * truncation_order, derivative_map);
// if the order of the far field expansion is greater than the
// local one we are adding onto, then increase the order.
@@ -520,8 +500,8 @@ void MultFarFieldExpansion<TKernelAux>::TranslateToLocal
}
// compute Gaussian derivative
pos_arrtmp.Init(sea_->get_max_total_num_coeffs());
neg_arrtmp.Init(sea_->get_max_total_num_coeffs());
pos_arrtmp.set_size(sea_->get_max_total_num_coeffs());
neg_arrtmp.set_size(sea_->get_max_total_num_coeffs());
// compute center difference divided by bw_times_sqrt_two;
for(index_t j = 0; j < dimension; j++) {
@@ -529,26 +509,26 @@ void MultFarFieldExpansion<TKernelAux>::TranslateToLocal
}
// compute required partial derivatives
ka_->ComputeDirectionalDerivatives(cent_diff, &derivative_map,
ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map,
2 * truncation_order);
ArrayList<short int> beta_plus_alpha;
beta_plus_alpha.Init(dimension);
std::vector<short int> beta_plus_alpha;
beta_plus_alpha.reserve(dimension);
// get the order of traversal for the given order of approximation
const ArrayList<short int> &traversal_order =
const std::vector<short 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<short int> &beta_mapping = sea_->get_multiindex(index);
const std::vector<short 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<short int> &alpha_mapping =
const std::vector<short 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];
@@ -567,7 +547,7 @@ void MultFarFieldExpansion<TKernelAux>::TranslateToLocal
} // end of k-loop
} // end of j-loop
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
arma::vec 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]) *
fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion.h
+13 -18
View File
@@ -12,7 +12,8 @@
#ifndef MULT_LOCAL_EXPANSION
#define MULT_LOCAL_EXPANSION
#include "fastlib/fastlib.h"
#include <fastlib/fastlib.h>
#include "kernel_aux.h"
#include "mult_series_expansion_aux.h"
@@ -28,10 +29,10 @@ class MultLocalExpansion {
private:
/** The center of the expansion */
Vector center_;
arma::vec center_;
/** The coefficients */
Vector coeffs_;
arma::vec coeffs_;
/** order */
int order_;
@@ -45,12 +46,6 @@ class MultLocalExpansion {
/** pointer to the precomputed constants inside kernel auxiliary object */
const typename TKernelAux::TSeriesExpansionAux *sea_;
OT_DEF(MultLocalExpansion) {
OT_MY_OBJECT(center_);
OT_MY_OBJECT(coeffs_);
OT_MY_OBJECT(order_);
}
public:
// getters and setters
@@ -59,12 +54,12 @@ class MultLocalExpansion {
double bandwidth_sq() const { return kernel_->bandwidth_sq(); }
/** Get the center of expansion */
Vector* get_center() { return &center_; }
const Vector* get_center() const { return &center_; }
arma::vec& get_center() { return center_; }
const arma::vec& get_center() const { return center_; }
/** Get the coefficients */
const Vector& get_coeffs() const { return coeffs_; }
arma::vec& get_coeffs() { return coeffs_; }
const arma::vec& get_coeffs() const { return coeffs_; }
/** Get the approximation order */
int get_order() const { return order_; }
@@ -81,26 +76,26 @@ class MultLocalExpansion {
* Accumulates the local moment represented by the given reference
* data into the coefficients
*/
void AccumulateCoeffs(const Matrix& data, const Vector& weights,
void AccumulateCoeffs(const arma::mat& data, const arma::vec& weights,
int begin, int end, int order);
/**
* This does not apply for local coefficients.
*/
void RefineCoeffs(const Matrix& data, const Vector& weights,
void RefineCoeffs(const arma::mat& data, const arma::vec& weights,
int begin, int end, int order) { }
/**
* Evaluates the local coefficients at the given point
*/
double EvaluateField(const Matrix& data, int row_num) const;
double EvaluateField(const Vector& x_q) const;
double EvaluateField(const arma::mat& data, int row_num) const;
double EvaluateField(const arma::vec& x_q) const;
/**
* Initializes the current local expansion object with the given
* center.
*/
void Init(const Vector& center, const TKernelAux &sea);
void Init(const arma::vec& center, const TKernelAux &sea);
void Init(const TKernelAux &sea);
/**
fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion_impl.h
+61 -78
View File
@@ -6,8 +6,8 @@
#define MULT_LOCAL_EXPANSION_IMPL_H
template<typename TKernelAux>
void MultLocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
const Vector& weights,
void MultLocalExpansion<TKernelAux>::AccumulateCoeffs(const arma::mat& data,
const arma::vec& weights,
int begin, int end,
int order) {
@@ -19,23 +19,20 @@ void MultLocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
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());
const arma::vec& neg_inv_multiindex_factorials = sea_->get_neg_inv_multiindex_factorials();
// declare deritave mapping
Matrix derivative_map;
ka_->AllocateDerivativeMap(dim, order, &derivative_map);
arma::mat derivative_map;
ka_->AllocateDerivativeMap(dim, order, derivative_map);
// some temporary variables
Vector x_r_minus_x_Q;
x_r_minus_x_Q.Init(dim);
arma::vec x_r_minus_x_Q(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<short int> &traversal_order =
const std::vector<short int> &traversal_order =
sea_->traversal_mapping_[order];
// for each data point,
@@ -43,17 +40,17 @@ void MultLocalExpansion<TKernelAux>::AccumulateCoeffs(const Matrix& data,
// 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)) /
x_r_minus_x_Q[d] = (center_[d] - data(d, r)) /
bandwidth_factor;
}
// precompute necessary partial derivatives based on coordinate difference
ka_->ComputeDirectionalDerivatives(x_r_minus_x_Q, &derivative_map, order);
ka_->ComputeDirectionalDerivatives(x_r_minus_x_Q, derivative_map, order);
// 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<short int> &mapping = sea_->get_multiindex(index);
const std::vector<short int>& mapping = sea_->get_multiindex(index);
double partial_derivative =
ka_->ComputePartialDerivative(derivative_map, mapping);
coeffs_[index] += neg_inv_multiindex_factorials[index] * weights[r] *
@@ -74,7 +71,7 @@ void MultLocalExpansion<TKernelAux>::PrintDebug(const char *name,
fprintf(stream, "Local expansion\n");
fprintf(stream, "Center: ");
for (index_t i = 0; i < center_.length(); i++) {
for (index_t i = 0; i < center_.n_elem; i++) {
fprintf(stream, "%g ", center_[i]);
}
fprintf(stream, "\n");
@@ -88,7 +85,7 @@ void MultLocalExpansion<TKernelAux>::PrintDebug(const char *name,
fprintf(stream, ") = \\sum\\limits_{x_r \\in R} K(||x_q - x_r||) = ");
for (index_t i = 0; i < total_num_coeffs; i++) {
ArrayList<short int> mapping = sea_->get_multiindex(i);
std::vector<short int> mapping = sea_->get_multiindex(i);
fprintf(stream, "%g", coeffs_[i]);
for(index_t d = 0; d < dim; d++) {
@@ -103,7 +100,7 @@ void MultLocalExpansion<TKernelAux>::PrintDebug(const char *name,
}
template<typename TKernelAux>
double MultLocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
double MultLocalExpansion<TKernelAux>::EvaluateField(const arma::mat& data,
int row_num) const {
// if there are no local coefficients, then return 0
@@ -124,16 +121,13 @@ double MultLocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
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<short int> heads;
heads.Init(dim + 1);
arma::vec x_Q_to_x_q(dim);
arma::vec tmp(sea_->get_max_total_num_coeffs());
std::vector<short int> heads(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;
x_Q_to_x_q[i] = (data(i, row_num) - center_[i]) / bandwidth_factor;
}
for(index_t i = 0; i < dim; i++)
@@ -143,21 +137,21 @@ double MultLocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
tmp[0] = 1.0;
// get the order of traversal for the given order of approximation
const ArrayList<short int> &traversal_order =
const std::vector<short int>& traversal_order =
sea_->traversal_mapping_[order_];
for(index_t i = 1; i < total_num_coeffs; i++) {
int index = traversal_order[i];
const ArrayList<short int> &lower_mappings =
const std::vector<short 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<short int> &mapping = sea_->multiindex_mapping_[index];
const ArrayList<short int> &direct_ancestor_mapping =
const std::vector<short int>& mapping = sea_->multiindex_mapping_[index];
const std::vector<short 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]) {
@@ -177,7 +171,7 @@ double MultLocalExpansion<TKernelAux>::EvaluateField(const Matrix& data,
}
template<typename TKernelAux>
double MultLocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
double MultLocalExpansion<TKernelAux>::EvaluateField(const arma::vec& x_q) const {
// if there are no local coefficients, then return 0
if(order_ < 0) {
@@ -197,12 +191,9 @@ double MultLocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
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<short int> heads;
heads.Init(dim + 1);
arma::vec x_Q_to_x_q(dim);
arma::vec tmp(sea_->get_max_total_num_coeffs());
std::vector<short int> heads(dim + 1);
// compute (x_q - x_Q) / (sqrt(2h^2))
for(index_t i = 0; i < dim; i++) {
@@ -216,19 +207,19 @@ double MultLocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
tmp[0] = 1.0;
// get the order of traversal for the given order of approximation
ArrayList<short int> &traversal_order = sea_->traversal_mapping_[order_];
std::vector<short int>& traversal_order = sea_->traversal_mapping_[order_];
for(index_t i = 1; i < total_num_coeffs; i++) {
int index = traversal_order[i];
ArrayList<short int> &lower_mappings = sea_->lower_mapping_index_[index];
std::vector<short 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<short int> &mapping = sea_->multiindex_mapping_[index];
const ArrayList<short int> &direct_ancestor_mapping =
const std::vector<short int> &mapping = sea_->multiindex_mapping_[index];
const std::vector<short 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]) {
@@ -248,19 +239,18 @@ double MultLocalExpansion<TKernelAux>::EvaluateField(const Vector& x_q) const {
}
template<typename TKernelAux>
void MultLocalExpansion<TKernelAux>::Init(const Vector& center,
void MultLocalExpansion<TKernelAux>::Init(const arma::vec& center,
const TKernelAux &ka) {
// copy kernel type, center, and bandwidth squared
kernel_ = &(ka.kernel_);
center_.Copy(center);
center_ = center;
order_ = -1;
sea_ = &(ka.sea_);
ka_ = &ka;
// initialize coefficient array
coeffs_.Init(sea_->get_max_total_num_coeffs());
coeffs_.SetZero();
coeffs_.zeros(sea_->get_max_total_num_coeffs());
}
template<typename TKernelAux>
@@ -269,13 +259,10 @@ 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());
// initialize coefficient array
center_.zeros(sea_->get_dimension());
order_ = -1;
ka_ = &ka;
// initialize coefficient array
coeffs_.Init(sea_->get_max_total_num_coeffs());
coeffs_.SetZero();
}
template<typename TKernelAux>
@@ -295,11 +282,11 @@ 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;
arma::vec pos_arrtmp, neg_arrtmp;
arma::mat derivative_map;
arma::vec& far_center;
arma::vec cent_diff;
arma::vec& far_coeffs;
int dimension = sea_->get_dimension();
int far_order = se.get_order();
int total_num_coeffs = sea_->get_total_num_coeffs(far_order);
@@ -309,9 +296,9 @@ void MultLocalExpansion<TKernelAux>::TranslateFromFarField
ka_->AllocateDerivativeMap(dimension, 2 * order_, &derivative_map);
// get center and coefficients for far field expansion
far_center.Alias(*(se.get_center()));
far_coeffs.Alias(se.get_coeffs());
cent_diff.Init(dimension);
far_center = se.get_center();
far_coeffs = se.get_coeffs();
cent_diff.set_size(dimension);
// if the order of the far field expansion is greater than the
// local one we are adding onto, then increase the order.
@@ -320,8 +307,8 @@ void MultLocalExpansion<TKernelAux>::TranslateFromFarField
}
// compute Gaussian derivative
pos_arrtmp.Init(total_num_coeffs);
neg_arrtmp.Init(total_num_coeffs);
pos_arrtmp.set_size(total_num_coeffs);
neg_arrtmp.set_size(total_num_coeffs);
// compute center difference divided by bw_times_sqrt_two;
for(index_t j = 0; j < dimension; j++) {
@@ -330,22 +317,22 @@ void MultLocalExpansion<TKernelAux>::TranslateFromFarField
// compute required partial derivatives
ka_->ComputeDirectionalDerivatives(cent_diff, &derivative_map, 2 * order_);
ArrayList<short int> beta_plus_alpha;
beta_plus_alpha.Init(dimension);
std::vector<short int> beta_plus_alpha;
beta_plus_alpha.reserve(dimension);
// get the order of traversal for the given order of approximation
ArrayList<short int> &traversal_order = sea_->traversal_mapping_[far_order];
std::vector<short 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<short int> beta_mapping = sea_->get_multiindex(index_j);
std::vector<short 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<short int> alpha_mapping = sea_->get_multiindex(index_k);
std::vector<short 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];
}
@@ -363,7 +350,7 @@ void MultLocalExpansion<TKernelAux>::TranslateFromFarField
} // end of k-loop
} // end of j-loop
Vector C_k_neg = sea_->get_neg_inv_multiindex_factorials();
arma::vec 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]) *
@@ -382,28 +369,25 @@ void MultLocalExpansion<TKernelAux>::TranslateToLocal(MultLocalExpansion &se) {
// 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()));
arma::vec& new_center = se.get_center();
int prev_order = se.get_order();
int total_num_coeffs = sea_->get_total_num_coeffs(order_);
const ArrayList<short int> *upper_mapping_index =
const std::vector<std::vector<short int> >& upper_mapping_index =
sea_->get_upper_mapping_index();
Vector new_coeffs;
new_coeffs.Alias(se.get_coeffs());
arma::vec& new_coeffs = se.get_coeffs();
// dimension
int dim = sea_->get_dimension();
// temporary variable
ArrayList<short int> tmp_storage;
tmp_storage.Init(dim);
std::vector<short int> tmp_storage;
tmp_storage.reserve(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);
arma::vec center_diff(dim);
for(index_t d = 0; d < dim; d++) {
center_diff[d] = (new_center[d] - center_[d]) / bandwidth_factor;
}
@@ -415,19 +399,18 @@ void MultLocalExpansion<TKernelAux>::TranslateToLocal(MultLocalExpansion &se) {
}
// inverse multiindex factorials
Vector C_k;
C_k.Alias(sea_->get_inv_multiindex_factorials());
const arma::vec& C_k = sea_->get_inv_multiindex_factorials();
// get the order of traversal for the given order of approximation
const ArrayList<short int> &traversal_order =
const std::vector<short 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<short int> &alpha_mapping = sea_->get_multiindex(index_j);
const ArrayList<short int> &upper_mappings_for_alpha =
const std::vector<short int> &alpha_mapping = sea_->get_multiindex(index_j);
const std::vector<short int> &upper_mappings_for_alpha =
upper_mapping_index[index_j];
double pos_coeffs = 0;
double neg_coeffs = 0;
@@ -438,7 +421,7 @@ void MultLocalExpansion<TKernelAux>::TranslateToLocal(MultLocalExpansion &se) {
break;
}
const ArrayList<short int> &beta_mapping =
const std::vector<short int> &beta_mapping =
sea_->get_multiindex(upper_mappings_for_alpha[k]);
int flag = 0;
double diff1 = 1.0;
fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_series_expansion_aux.h
+62 -78
View File
@@ -4,7 +4,8 @@
#ifndef MULT_SERIES_EXPANSION_AUX_H
#define MULT_SERIES_EXPANSION_AUX_H
#include "fastlib/fastlib.h"
#include <fastlib/fastlib.h>
#include <armadillo>
/**
* Series expansion class for multiplicative kernel functions
@@ -18,63 +19,48 @@ class MultSeriesExpansionAux {
int max_order_;
Vector factorials_;
arma::vec factorials_;
ArrayList<int> list_total_num_coeffs_;
std::vector<int> list_total_num_coeffs_;
Vector inv_multiindex_factorials_;
arma::vec inv_multiindex_factorials_;
Vector neg_inv_multiindex_factorials_;
arma::vec neg_inv_multiindex_factorials_;
Matrix multiindex_combination_;
arma::mat multiindex_combination_;
ArrayList< ArrayList<short int> > multiindex_mapping_;
std::vector< std::vector<short int> > multiindex_mapping_;
/**
* for each i-th multiindex m_i, store the positions of the j-th
* multiindex mapping such that m_i - m_j >= 0 (the difference in
* all coordinates is nonnegative).
*/
ArrayList< ArrayList<short int> > lower_mapping_index_;
std::vector< std::vector<short int> > lower_mapping_index_;
/**
* for each i-th multiindex m_i, store the positions of the j-th
* multiindex mapping such that m_i - m_j <= 0 (the difference in
* all coordinates is nonpositive).
*/
ArrayList< ArrayList<short int> > upper_mapping_index_;
std::vector< std::vector<short int> > upper_mapping_index_;
/** row index is for n, column index is for k */
Matrix n_choose_k_;
arma::mat n_choose_k_;
/**
* For each i-th order, store the positions of the coefficient
* array to traverse.
*/
ArrayList< ArrayList<short int> > traversal_mapping_;
OT_DEF_BASIC(MultSeriesExpansionAux) {
OT_MY_OBJECT(dim_);
OT_MY_OBJECT(max_order_);
OT_MY_OBJECT(factorials_);
OT_MY_OBJECT(list_total_num_coeffs_);
OT_MY_OBJECT(inv_multiindex_factorials_);
OT_MY_OBJECT(neg_inv_multiindex_factorials_);
OT_MY_OBJECT(multiindex_combination_);
OT_MY_OBJECT(multiindex_mapping_);
OT_MY_OBJECT(lower_mapping_index_);
OT_MY_OBJECT(upper_mapping_index_);
OT_MY_OBJECT(n_choose_k_);
OT_MY_OBJECT(traversal_mapping_);
}
std::vector< std::vector<short int> > traversal_mapping_;
public:
void ComputeFactorials() {
factorials_.Init(2 * max_order_ + 1);
factorials_.set_size(2 * max_order_ + 1);
factorials_[0] = 1;
for(index_t t = 1; t < factorials_.length(); t++) {
for(index_t t = 1; t < factorials_.n_elem; t++) {
factorials_[t] = t * factorials_[t - 1];
}
}
@@ -83,15 +69,15 @@ class MultSeriesExpansionAux {
// initialize the index
int limit = 2 * max_order_;
traversal_mapping_.Init(limit + 1);
traversal_mapping_.reserve(limit + 1);
for(index_t i = 0; i <= max_order_; i++) {
traversal_mapping_[i].Init();
traversal_mapping_[i].clear();
for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) {
const ArrayList<short int> &mapping = multiindex_mapping_[j];
const std::vector<short int>& mapping = multiindex_mapping_[j];
int flag = 0;
for(index_t d = 0; d < dim_; d++) {
@@ -102,7 +88,7 @@ class MultSeriesExpansionAux {
}
if(flag == 0) {
(traversal_mapping_[i]).PushBackCopy(j);
(traversal_mapping_[i]).push_back(j);
}
} // end of j-loop
} // end of i-loop
@@ -110,19 +96,19 @@ class MultSeriesExpansionAux {
void ComputeLowerMappingIndex() {
ArrayList<short int> diff;
diff.Init(dim_);
std::vector<short int> diff;
diff.reserve(dim_);
// initialize the index
int limit = 2 * max_order_;
lower_mapping_index_.Init(list_total_num_coeffs_[limit]);
lower_mapping_index_.reserve(list_total_num_coeffs_[limit]);
for(index_t i = 0; i < list_total_num_coeffs_[limit]; i++) {
const ArrayList<short int> &outer_mapping = multiindex_mapping_[i];
lower_mapping_index_[i].Init();
const std::vector<short int>& outer_mapping = multiindex_mapping_[i];
lower_mapping_index_[i].clear();
for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) {
const ArrayList<short int> &inner_mapping = multiindex_mapping_[j];
const std::vector<short int>& inner_mapping = multiindex_mapping_[j];
int flag = 0;
for(index_t d = 0; d < dim_; d++) {
@@ -135,7 +121,7 @@ class MultSeriesExpansionAux {
}
if(flag == 0) {
(lower_mapping_index_[i]).PushBackCopy(j);
(lower_mapping_index_[i]).push_back(j);
}
} // end of j-loop
} // end of i-loop
@@ -144,28 +130,27 @@ class MultSeriesExpansionAux {
void ComputeMultiindexCombination() {
int limit = 2 * max_order_;
multiindex_combination_.Init(list_total_num_coeffs_[limit],
list_total_num_coeffs_[limit]);
multiindex_combination_.set_size(list_total_num_coeffs_[limit],
list_total_num_coeffs_[limit]);
for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) {
// beta mapping
const ArrayList<short int> &beta_mapping = multiindex_mapping_[j];
const std::vector<short int>& beta_mapping = multiindex_mapping_[j];
for(index_t k = 0; k < list_total_num_coeffs_[limit]; k++) {
// alpha mapping
const ArrayList<short int> &alpha_mapping = multiindex_mapping_[k];
const std::vector<short int>& alpha_mapping = multiindex_mapping_[k];
// initialize the factor to 1
multiindex_combination_.set(j, k, 1);
multiindex_combination_(j, k) = 1;
for(index_t i = 0; i < dim_; i++) {
multiindex_combination_.set
(j, k, multiindex_combination_.get(j, k) *
n_choose_k_.get(beta_mapping[i], alpha_mapping[i]));
multiindex_combination_(j, k) *=
n_choose_k_(beta_mapping[i], alpha_mapping[i]);
if(multiindex_combination_.get(j, k) == 0)
if(multiindex_combination_(j, k) == 0)
break;
}
}
@@ -174,19 +159,19 @@ class MultSeriesExpansionAux {
void ComputeUpperMappingIndex() {
ArrayList<short int> diff;
diff.Init(dim_);
std::vector<short int> diff;
diff.reserve(dim_);
// initialize the index
int limit = 2 * max_order_;
upper_mapping_index_.Init(list_total_num_coeffs_[limit]);
upper_mapping_index_.reserve(list_total_num_coeffs_[limit]);
for(index_t i = 0; i < list_total_num_coeffs_[limit]; i++) {
const ArrayList<short int> &outer_mapping = multiindex_mapping_[i];
upper_mapping_index_[i].Init();
const std::vector<short int>& outer_mapping = multiindex_mapping_[i];
upper_mapping_index_[i].clear();
for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) {
const ArrayList<short int> &inner_mapping = multiindex_mapping_[j];
const std::vector<short int>& inner_mapping = multiindex_mapping_[j];
int flag = 0;
for(index_t d = 0; d < dim_; d++) {
@@ -199,7 +184,7 @@ class MultSeriesExpansionAux {
}
if(flag == 0) {
(upper_mapping_index_[i]).PushBackCopy(j);
(upper_mapping_index_[i]).push_back(j);
}
} // end of j-loop
} // end of i-loop
@@ -218,40 +203,40 @@ class MultSeriesExpansionAux {
return list_total_num_coeffs_[max_order_];
}
const Vector& get_inv_multiindex_factorials() const {
const arma::vec& get_inv_multiindex_factorials() const {
return inv_multiindex_factorials_;
}
const ArrayList<short int> * get_lower_mapping_index() const {
return lower_mapping_index_.begin();
const std::vector<std::vector<short int> >& get_lower_mapping_index() const {
return lower_mapping_index_;
}
int get_max_order() const {
return max_order_;
}
const ArrayList<short int> & get_multiindex(int pos) const {
const std::vector<short int>& get_multiindex(int pos) const {
return multiindex_mapping_[pos];
}
const ArrayList<short int> * get_multiindex_mapping() const {
return multiindex_mapping_.begin();
const std::vector<std::vector<short int> >& get_multiindex_mapping() const {
return multiindex_mapping_;
}
const Vector& get_neg_inv_multiindex_factorials() const {
const arma::vec& get_neg_inv_multiindex_factorials() const {
return neg_inv_multiindex_factorials_;
}
double get_n_choose_k(int n, int k) const {
return n_choose_k_.get(n, (int) math::ClampNonNegative(k));
return n_choose_k_(n, (int) math::ClampNonNegative(k));
}
double get_n_multichoose_k_by_pos(int n, int k) const {
return multiindex_combination_.get(n, k);
return multiindex_combination_(n, k);
}
const ArrayList<short int> * get_upper_mapping_index() const {
return upper_mapping_index_.begin();
const std::vector<std::vector<short int> >& get_upper_mapping_index() const {
return upper_mapping_index_;
}
// interesting functions
@@ -259,7 +244,7 @@ class MultSeriesExpansionAux {
/**
* Computes the position of the given multiindex
*/
int ComputeMultiindexPosition(const ArrayList<short int> &multiindex) const {
int ComputeMultiindexPosition(const std::vector<short int>& multiindex) const {
int index = 0;
// using Horner's rule
@@ -304,7 +289,7 @@ class MultSeriesExpansionAux {
// compute the list of total number of coefficients for p-th order
// expansion
int limit = 2 * max_order_;
list_total_num_coeffs_.Init(limit + 1);
list_total_num_coeffs_.reserve(limit + 1);
list_total_num_coeffs_[0] = 1;
for(index_t p = 1; p <= limit; p++) {
list_total_num_coeffs_[p] = (int) pow(p + 1, dim);
@@ -316,15 +301,14 @@ class MultSeriesExpansionAux {
// allocate space for inverse factorial and
// negative inverse factorials and multiindex mapping and n_choose_k
// and multiindex_combination precomputed factors
inv_multiindex_factorials_.Init(list_total_num_coeffs_[limit]);
neg_inv_multiindex_factorials_.Init(list_total_num_coeffs_[limit]);
multiindex_mapping_.Init(list_total_num_coeffs_[limit]);
(multiindex_mapping_[0]).Init(dim_);
inv_multiindex_factorials_.set_size(list_total_num_coeffs_[limit]);
neg_inv_multiindex_factorials_.set_size(list_total_num_coeffs_[limit]);
multiindex_mapping_.reserve(list_total_num_coeffs_[limit]);
(multiindex_mapping_[0]).reserve(dim_);
for(index_t j = 0; j < dim; j++) {
(multiindex_mapping_[0])[j] = 0;
}
n_choose_k_.Init(dim * (limit + 1), dim * (limit + 1));
n_choose_k_.SetZero();
n_choose_k_.zeros(dim * (limit + 1), dim * (limit + 1));
// compute inverse factorial and negative inverse factorials and
// multiindex mappings...
@@ -353,7 +337,7 @@ class MultSeriesExpansionAux {
div++;
// copy multiindex from old to the new position
multiindex_mapping_[i].InitCopy(multiindex_mapping_[i - step]);
multiindex_mapping_[i] = multiindex_mapping_[i - step];
(multiindex_mapping_[i])[k] = (multiindex_mapping_[i])[k] + 1;
}
}
@@ -361,9 +345,9 @@ class MultSeriesExpansionAux {
}
// compute n choose k's
for(index_t j = 0; j < n_choose_k_.n_rows(); j++) {
for(index_t k = 0; k < n_choose_k_.n_cols(); k++) {
n_choose_k_.set(j, k, math::BinomialCoefficient(j, k));
for(index_t j = 0; j < n_choose_k_.n_rows; j++) {
for(index_t k = 0; k < n_choose_k_.n_cols; k++) {
n_choose_k_(j, k) = math::BinomialCoefficient(j, k);
}
}
fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.cc
+14 -15
View File
@@ -1,6 +1,6 @@
#include "series_expansion_aux.h"
const Vector& SeriesExpansionAux::get_inv_multiindex_factorials() const {
const arma::vec& SeriesExpansionAux::get_inv_multiindex_factorials() const {
return inv_multiindex_factorials_;
}
@@ -8,9 +8,9 @@ int SeriesExpansionAux::get_max_total_num_coeffs() const {
return list_total_num_coeffs_[max_order_];
}
const std::vector < short int > &SeriesExpansionAux::get_lower_mapping_index()
const std::vector < std::vector<short int> >& SeriesExpansionAux::get_lower_mapping_index()
const {
return lower_mapping_index_.front();
return lower_mapping_index_;
}
int SeriesExpansionAux::get_max_order() const {
@@ -22,21 +22,21 @@ const std::vector < short int > &SeriesExpansionAux::get_multiindex(int pos)
return multiindex_mapping_[pos];
}
const std::vector < short int > &SeriesExpansionAux::get_multiindex_mapping()
const std::vector < std::vector<short int> >& SeriesExpansionAux::get_multiindex_mapping()
const {
return multiindex_mapping_.front();
return multiindex_mapping_;
}
const Vector& SeriesExpansionAux::get_neg_inv_multiindex_factorials() const {
const arma::vec& SeriesExpansionAux::get_neg_inv_multiindex_factorials() const {
return neg_inv_multiindex_factorials_;
}
double SeriesExpansionAux::get_n_choose_k(int n, int k) const {
return n_choose_k_.get(n, (int) math::ClampNonNegative(k));
return n_choose_k_(n, (int) math::ClampNonNegative(k));
}
double SeriesExpansionAux::get_n_multichoose_k_by_pos(int n, int k) const {
return multiindex_combination_.get(n, k);
return multiindex_combination_(n, k);
}
int SeriesExpansionAux::get_total_num_coeffs(int order) const {
@@ -44,10 +44,10 @@ int SeriesExpansionAux::get_total_num_coeffs(int order) const {
return list_total_num_coeffs_[order];
}
const std::vector < short int > &SeriesExpansionAux::get_upper_mapping_index()
const std::vector < std::vector< short int > >& SeriesExpansionAux::get_upper_mapping_index()
const {
return upper_mapping_index_.front();
return upper_mapping_index_;
}
int SeriesExpansionAux::ComputeMultiindexPosition
@@ -112,15 +112,14 @@ void SeriesExpansionAux::Init(int max_order, int dim) {
// allocate space for inverse factorial and
// negative inverse factorials and multiindex mapping and n_choose_k
// and multiindex_combination precomputed factors
inv_multiindex_factorials_.Init(list_total_num_coeffs_[limit - 1]);
neg_inv_multiindex_factorials_.Init(list_total_num_coeffs_[limit - 1]);
inv_multiindex_factorials_.set_size(list_total_num_coeffs_[limit - 1]);
neg_inv_multiindex_factorials_.set_size(list_total_num_coeffs_[limit - 1]);
multiindex_mapping_.reserve(list_total_num_coeffs_[limit - 1]);
(multiindex_mapping_[0]).reserve(dim_);
for(j = 0; j < dim; j++) {
(multiindex_mapping_[0])[j] = 0;
}
n_choose_k_.Init((limit - 1) + dim + 1, (limit - 1) + dim + 1);
n_choose_k_.SetZero();
n_choose_k_.zeros((limit - 1) + dim + 1, (limit - 1) + dim + 1);
// initialization of temporary variables for computation...
heads.reserve(dim + 1);
@@ -156,7 +155,7 @@ void SeriesExpansionAux::Init(int max_order, int dim) {
// compute n choose k's
for(j = 0; j <= 2 * max_order + dim; j++) {
for(k = 0; k <= 2 * max_order + dim; k++) {
n_choose_k_.set(j, k, math::BinomialCoefficient(j, k));
n_choose_k_(j, k) = math::BinomialCoefficient(j, k);
}
}
fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.h
+20 -34
View File
@@ -4,7 +4,8 @@
#ifndef SERIES_EXPANSION_AUX
#define SERIES_EXPANSION_AUX
#include "fastlib/fastlib.h"
#include <fastlib/fastlib.h>
#include <armadillo>
/**
* Series expansion class.
@@ -17,15 +18,15 @@ class SeriesExpansionAux {
int max_order_;
Vector factorials_;
arma::vec factorials_;
std::vector<int> list_total_num_coeffs_;
Vector inv_multiindex_factorials_;
arma::vec inv_multiindex_factorials_;
Vector neg_inv_multiindex_factorials_;
arma::vec neg_inv_multiindex_factorials_;
Matrix multiindex_combination_;
arma::mat multiindex_combination_;
std::vector< std::vector<short int> > multiindex_mapping_;
@@ -44,29 +45,15 @@ class SeriesExpansionAux {
std::vector< std::vector<short int> > upper_mapping_index_;
/** row index is for n, column index is for k */
Matrix n_choose_k_;
OT_DEF_BASIC(SeriesExpansionAux) {
OT_MY_OBJECT(dim_);
OT_MY_OBJECT(max_order_);
OT_MY_OBJECT(factorials_);
OT_MY_OBJECT(list_total_num_coeffs_);
OT_MY_OBJECT(inv_multiindex_factorials_);
OT_MY_OBJECT(neg_inv_multiindex_factorials_);
OT_MY_OBJECT(multiindex_combination_);
OT_MY_OBJECT(multiindex_mapping_);
OT_MY_OBJECT(lower_mapping_index_);
OT_MY_OBJECT(upper_mapping_index_);
OT_MY_OBJECT(n_choose_k_);
}
arma::mat n_choose_k_;
public:
void ComputeFactorials() {
factorials_.Init(2 * max_order_ + 1);
factorials_.set_size(2 * max_order_ + 1);
factorials_[0] = 1;
for(index_t t = 1; t < factorials_.length(); t++) {
for(index_t t = 1; t < factorials_.n_elem; t++) {
factorials_[t] = t * factorials_[t - 1];
}
}
@@ -107,8 +94,8 @@ class SeriesExpansionAux {
void ComputeMultiindexCombination() {
int limit = 2 * max_order_;
multiindex_combination_.Init(list_total_num_coeffs_[limit],
list_total_num_coeffs_[limit]);
multiindex_combination_.set_size(list_total_num_coeffs_[limit],
list_total_num_coeffs_[limit]);
for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) {
@@ -121,14 +108,13 @@ class SeriesExpansionAux {
const std::vector<short int> &alpha_mapping = multiindex_mapping_[k];
// initialize the factor to 1
multiindex_combination_.set(j, k, 1);
multiindex_combination_(j, k) = 1;
for(index_t i = 0; i < dim_; i++) {
multiindex_combination_.set
(j, k, multiindex_combination_.get(j, k) *
n_choose_k_.get(beta_mapping[i], alpha_mapping[i]));
multiindex_combination_(j, k) *=
n_choose_k_(beta_mapping[i], alpha_mapping[i]);
if(multiindex_combination_.get(j, k) == 0)
if(multiindex_combination_(j, k) == 0)
break;
}
}
@@ -176,23 +162,23 @@ class SeriesExpansionAux {
int get_max_total_num_coeffs() const;
const Vector& get_inv_multiindex_factorials() const;
const arma::vec& get_inv_multiindex_factorials() const;
const std::vector< short int > & get_lower_mapping_index() const;
const std::vector< std::vector<short int> > & get_lower_mapping_index() const;
int get_max_order() const;
const std::vector< short int > & get_multiindex(int pos) const;
const std::vector< short int > & get_multiindex_mapping() const;
const std::vector< std::vector<short int> > & get_multiindex_mapping() const;
const Vector& get_neg_inv_multiindex_factorials() const;
const arma::vec& get_neg_inv_multiindex_factorials() const;
double get_n_choose_k(int n, int k) const;
double get_n_multichoose_k_by_pos(int n, int k) const;
const std::vector< short int > & get_upper_mapping_index() const;
const std::vector< std::vector<short int> >& get_upper_mapping_index() const;
// interesting functions