diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion.h index a7317d3dc8..acb292a53a 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion.h @@ -12,7 +12,8 @@ #ifndef FARFIELD_EXPANSION #define FARFIELD_EXPANSION -#include "fastlib/fastlib.h" +#include + #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 ¢er_; } - - const Vector* get_center() const { return ¢er_; } + 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 ¢er) { - - 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 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; diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion_impl.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion_impl.h index 6e87ee8a1d..217e112e65 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion_impl.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/farfield_expansion_impl.h @@ -6,33 +6,31 @@ #define FARFIELD_EXPANSION_IMPL_H template -void FarFieldExpansion::Accumulate(const Vector &v, double weight, +void FarFieldExpansion::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 heads; - Vector x_r; + arma::Col 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::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::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::Accumulate(const Vector &v, double weight, } template -void FarFieldExpansion::AccumulateCoeffs(const Matrix& data, - const Vector& weights, +void FarFieldExpansion::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 heads; - Vector x_r; + arma::Col 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::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::AccumulateCoeffs(const Matrix& data, } template -void FarFieldExpansion::RefineCoeffs(const Matrix& data, - const Vector& weights, +void FarFieldExpansion::RefineCoeffs(const arma::mat& data, + const arma::vec& weights, int begin, int end, int order) { @@ -168,22 +163,19 @@ void FarFieldExpansion::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::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::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 &mapping = sea_->get_multiindex(i); + const std::vector &mapping = sea_->get_multiindex(i); tmp = 1; for(j = 0; j < dim; j++) { @@ -225,19 +217,16 @@ void FarFieldExpansion::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 -double FarFieldExpansion::EvaluateField(const Matrix& data, +double FarFieldExpansion::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 @@ -259,16 +248,14 @@ double FarFieldExpansion::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::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 &mapping = sea_->get_multiindex(j); + const std::vector &mapping = sea_->get_multiindex(j); double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping); double prod = coeffs_[j] * arrtmp; @@ -297,7 +284,7 @@ double FarFieldExpansion::EvaluateField(const double *x_q, } template -double FarFieldExpansion::MixField(const Matrix &data, +double FarFieldExpansion::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::MixField(const Matrix &data, // bandwidth factor and multiindex mapping stuffs double result; double bandwidth_factor = ka_->BandwidthFactor(bandwidth_sq()); - const ArrayList *multiindex_mapping = + const std::vector& multiindex_mapping = sea_->get_multiindex_mapping(); - const ArrayList *lower_mapping_index = + const std::vector& 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::MixField(const Matrix &data, double xi_xI, xj_xJ, diff; // temporary array - ArrayList beta_gamma_nu_eta_mapping; - ArrayList beta_nu_mapping; - ArrayList gamma_eta_mapping; - beta_nu_mapping.Init(dim); - gamma_eta_mapping.Init(dim); - beta_gamma_nu_eta_mapping.Init(dim); + std::vector beta_gamma_nu_eta_mapping; + std::vector beta_nu_mapping; + std::vector 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 &beta_mapping = multiindex_mapping[beta]; - const ArrayList &lower_mappings_for_beta = + const std::vector &beta_mapping = multiindex_mapping[beta]; + const std::vector &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 &nu_mapping = + const std::vector &nu_mapping = multiindex_mapping[lower_mappings_for_beta[nu]]; // beta - nu @@ -398,19 +377,19 @@ double FarFieldExpansion::MixField(const Matrix &data, for(index_t gamma = 0; gamma < total_num_coeffs3; gamma++) { - const ArrayList &gamma_mapping = multiindex_mapping[gamma]; - const ArrayList &lower_mappings_for_gamma = + const std::vector &gamma_mapping = multiindex_mapping[gamma]; + const std::vector &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 &eta_mapping = + const std::vector &eta_mapping = multiindex_mapping[lower_mappings_for_gamma[eta]]; for(index_t d = 0; d < dim; d++) { @@ -439,7 +418,7 @@ double FarFieldExpansion::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::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::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::ConvolveField // The bandwidth factor and the multiindex mapping stuffs. double bandwidth_factor = ka_->BandwidthFactor(bandwidth_sq()); - const ArrayList *multiindex_mapping = + const std::vector& multiindex_mapping = sea_->get_multiindex_mapping(); - const ArrayList *lower_mapping_index = + const std::vector& 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::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 alpha_minus_beta_mapping; - alpha_minus_beta_mapping.Init(dim); + std::vector 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 &alpha_mapping = multiindex_mapping[alpha]; - const ArrayList &lower_mappings_for_alpha = + const std::vector &alpha_mapping = multiindex_mapping[alpha]; + const std::vector &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 &beta_mapping = + const std::vector &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::ConvolveField // bandwidth factor and multiindex mapping stuffs double result; double bandwidth_factor = ka_->BandwidthFactor(bandwidth_sq()); - const ArrayList *multiindex_mapping = + const std::vector& multiindex_mapping = sea_->get_multiindex_mapping(); - const ArrayList *lower_mapping_index = + const std::vector& lower_mapping_index = sea_->get_lower_mapping_index(); // get the total number of coefficients and coefficients @@ -593,9 +567,8 @@ double FarFieldExpansion::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::ConvolveField double moment_i, moment_j, moment_k; // temporary array - ArrayList mu_nu_mapping; - ArrayList alpha_mu_eta_mapping; - ArrayList beta_gamma_nu_eta_mapping; - ArrayList alpha_mu_mapping; - ArrayList beta_nu_mapping; - ArrayList 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 mu_nu_mapping; + std::vector alpha_mu_eta_mapping; + std::vector beta_gamma_nu_eta_mapping; + std::vector alpha_mu_mapping; + std::vector beta_nu_mapping; + std::vector 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 &alpha_mapping = multiindex_mapping[alpha]; - const ArrayList &lower_mappings_for_alpha = + const std::vector& alpha_mapping = multiindex_mapping[alpha]; + const std::vector& 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 &mu_mapping = + const std::vector& mu_mapping = multiindex_mapping[lower_mappings_for_alpha[mu]]; // alpha - mu @@ -670,15 +638,15 @@ double FarFieldExpansion::ConvolveField for(index_t beta = 0; beta < total_num_coeffs2; beta++) { - const ArrayList &beta_mapping = multiindex_mapping[beta]; - const ArrayList &lower_mappings_for_beta = + const std::vector& beta_mapping = multiindex_mapping[beta]; + const std::vector& 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 &nu_mapping = + const std::vector &nu_mapping = multiindex_mapping[lower_mappings_for_beta[nu]]; // mu + nu and beta - nu @@ -689,20 +657,20 @@ double FarFieldExpansion::ConvolveField for(index_t gamma = 0; gamma < total_num_coeffs3; gamma++) { - const ArrayList &gamma_mapping = + const std::vector &gamma_mapping = multiindex_mapping[gamma]; - const ArrayList &lower_mappings_for_gamma = + const std::vector &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 &eta_mapping = + const std::vector& eta_mapping = multiindex_mapping[lower_mappings_for_gamma[eta]]; for(index_t d = 0; d < dim; d++) { @@ -766,19 +734,18 @@ double FarFieldExpansion::ConvolveField } template -void FarFieldExpansion::Init(const Vector& center, +void FarFieldExpansion::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 @@ -788,20 +755,18 @@ void FarFieldExpansion::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 template int FarFieldExpansion::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::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::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 &mapping = sea_->get_multiindex(i); + const std::vector &mapping = sea_->get_multiindex(i); fprintf(stream, "%g ", coeffs_[i]); fprintf(stream, "(-1)^("); @@ -899,24 +864,20 @@ void FarFieldExpansion::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 *multiindex_mapping = + const arma::vec& prev_coeffs = se.get_coeffs(); + const arma::vec& prev_center = se.get_center(); + const std::vector >& multiindex_mapping = sea_->get_multiindex_mapping(); - const ArrayList *lower_mapping_index = + const std::vector >& lower_mapping_index = sea_->get_lower_mapping_index(); - ArrayList tmp_storage; - Vector center_diff; - Vector inv_multiindex_factorials; - - center_diff.Init(dim); + std::vector 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::TranslateFromFarField for(index_t j = 0; j < total_num_coeffs; j++) { - const ArrayList &gamma_mapping = multiindex_mapping[j]; - const ArrayList &lower_mappings_for_gamma = + const std::vector& gamma_mapping = multiindex_mapping[j]; + const std::vector& 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 &inner_mapping = + const std::vector& inner_mapping = multiindex_mapping[lower_mappings_for_gamma[k]]; int flag = 0; @@ -987,22 +948,20 @@ template void FarFieldExpansion::TranslateToLocal (LocalExpansion &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::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::TranslateToLocal } // Compute required partial derivatives. - ka_->ComputeDirectionalDerivatives(cent_diff, &derivative_map, + ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map, 2 * truncation_order); - ArrayList beta_plus_alpha; - beta_plus_alpha.Init(dimension); + std::vector beta_plus_alpha; + beta_plus_alpha.reserve(dimension); for(index_t j = 0; j < total_num_coeffs; j++) { - const ArrayList &beta_mapping = sea_->get_multiindex(j); + const std::vector& 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 &alpha_mapping = sea_->get_multiindex(k); + const std::vector& 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::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]; } diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/inverse_pow_dist_kernel_aux.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/inverse_pow_dist_kernel_aux.h index 11959f9a3d..0d86a9cece 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/inverse_pow_dist_kernel_aux.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/inverse_pow_dist_kernel_aux.h @@ -11,8 +11,8 @@ class InversePowDistGradientKernelAux { private: void SubFrom_(index_t dimension, int decrement, - const ArrayList &subtract_from, - ArrayList &result) const { + const std::vector &subtract_from, + std::vector &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 \ - TFarFieldExpansion; + typedef FarFieldExpansion TFarFieldExpansion; typedef LocalExpansion 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 tmp_multiindex; - tmp_multiindex.Init(sea_.get_dimension()); + std::vector 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 &multiindex = sea_.get_multiindex(i); + const std::vector& 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 &multiindex = sea_.get_multiindex(i); + const std::vector& 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 &mapping) const { + double ComputePartialDerivative(const arma::mat& derivative_map, + const std::vector& 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 &subtract_from, - ArrayList &result) const { + const std::vector& subtract_from, + std::vector& 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 tmp_multiindex; - tmp_multiindex.Init(sea_.get_dimension()); + std::vector 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 &multiindex = sea_.get_multiindex(i); + const std::vector& 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 &multiindex = sea_.get_multiindex(i); + const std::vector& 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 &mapping) const { + double ComputePartialDerivative(const arma::mat& derivative_map, + const std::vector& mapping) const { - return derivative_map.get(sea_.ComputeMultiindexPosition(mapping), 0); + return derivative_map(sea_.ComputeMultiindexPosition(mapping), 0); } }; diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/kernel_aux.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/kernel_aux.h index 1118174c74..979114173d 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/kernel_aux.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/kernel_aux.h @@ -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 &mapping) const { + double ComputePartialDerivative(const arma::mat& derivative_map, + const std::vector& 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 &mapping) const { + double ComputePartialDerivative(const arma::mat& derivative_map, + const std::vector& 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 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 &mapping) const { + double ComputePartialDerivative(const arma::mat& derivative_map, + const std::vector& mapping) const { - return derivative_map.get(sea_.ComputeMultiindexPosition(mapping), 0); + return derivative_map(sea_.ComputeMultiindexPosition(mapping), 0); } template diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion.h index 87a9554a0a..4fe227b48f 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion.h @@ -12,7 +12,7 @@ #ifndef LOCAL_EXPANSION #define LOCAL_EXPANSION -#include "fastlib/fastlib.h" +#include #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 ¢er_; } - - const Vector* get_center() const { return ¢er_; } + 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); /** diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion_impl.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion_impl.h index c044d55afb..8c64eb6e23 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion_impl.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/local_expansion_impl.h @@ -7,32 +7,27 @@ template -void LocalExpansion::AccumulateCoeffs(const Matrix& data, - const Vector& weights, +void LocalExpansion::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::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 &mapping = sea_->get_multiindex(j); + const std::vector& 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::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::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 &mapping = sea_->get_multiindex(i); + const std::vector& mapping = sea_->get_multiindex(i); fprintf(stream, "%g", coeffs_[i]); for(index_t d = 0; d < dim; d++) { @@ -102,18 +97,17 @@ void LocalExpansion::PrintDebug(const char *name, } template -double LocalExpansion::EvaluateField(const Matrix& data, +double LocalExpansion::EvaluateField(const arma::mat& data, int row_num) const { - return EvaluateField(data.GetColumnPtr(row_num)); + return EvaluateField(data.colptr(row_num)); } template double LocalExpansion::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::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 heads; - heads.Init(dim + 1); + arma::vec x_Q_to_x_q(dim), tmp(total_num_coeffs); + std::vector 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::EvaluateField(const double *x_q) const { } template -void LocalExpansion::Init(const Vector& center, +void LocalExpansion::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 @@ -191,12 +181,11 @@ void LocalExpansion::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 @@ -216,13 +205,12 @@ template void LocalExpansion::TranslateFromFarField (const FarFieldExpansion &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::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::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::TranslateFromFarField } // compute required partial derivatives - ka_->ComputeDirectionalDerivatives(cent_diff, &derivative_map, 2 * order_); - ArrayList beta_plus_alpha; - beta_plus_alpha.Init(dimension); + ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map, 2 * order_); + std::vector beta_plus_alpha; + beta_plus_alpha.reserve(dimension); for(index_t j = 0; j < total_num_coeffs; j++) { - const ArrayList &beta_mapping = sea_->get_multiindex(j); + const std::vector& 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 &alpha_mapping = sea_->get_multiindex(k); + const std::vector& 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::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::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 *upper_mapping_index = + const std::vector >& 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 tmp_storage; - tmp_storage.Init(dim); + std::vector 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::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 &alpha_mapping = sea_->get_multiindex(j); - const ArrayList &upper_mappings_for_alpha = - upper_mapping_index[j]; + const std::vector& alpha_mapping = sea_->get_multiindex(j); + const std::vector& 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 &beta_mapping = + const std::vector& beta_mapping = sea_->get_multiindex(upper_mappings_for_alpha[k]); int flag = 0; double diff1 = 1.0; @@ -385,4 +367,3 @@ void LocalExpansion::TranslateToLocal(LocalExpansion &se) { } #endif - diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion.h index d3054e9449..084dd90302 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion.h @@ -12,8 +12,7 @@ #ifndef MULT_FARFIELD_EXPANSION #define MULT_FARFIELD_EXPANSION - -#include "fastlib/fastlib.h" +#include #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 ¢er_; } + arma::vec& get_center() { return center_; } - const Vector *get_center() const { return ¢er_; } + 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 ¢er) { + 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 - 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 - 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 &se, + void TranslateToLocal(MultLocalExpansion& se, int truncation_order); }; diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion_impl.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion_impl.h index 325aa2d6d3..27b9b7139d 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion_impl.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_farfield_expansion_impl.h @@ -6,26 +6,24 @@ #define MULT_FARFIELD_EXPANSION_IMPL_H template -void MultFarFieldExpansion::AccumulateCoeffs(const Matrix& data, - const Vector& weights, +void MultFarFieldExpansion::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::AccumulateCoeffs(const Matrix& data, } // get the order of traversal for the given order of approximation - const ArrayList &traversal_order = + const std::vector& traversal_order = sea_->traversal_mapping_[order_]; // Repeat for each reference point in this reference node. @@ -42,16 +40,16 @@ void MultFarFieldExpansion::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 &lower_mappings = + const std::vector& lower_mappings = sea_->lower_mapping_index_[index]; // from the direct descendant, recursively compute the multipole moments @@ -59,8 +57,8 @@ void MultFarFieldExpansion::AccumulateCoeffs(const Matrix& data, lower_mappings[lower_mappings.size() - 2]; int position = 0; - const ArrayList &mapping = sea_->multiindex_mapping_[index]; - const ArrayList &direct_ancestor_mapping = + const std::vector& mapping = sea_->multiindex_mapping_[index]; + const std::vector& 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::AccumulateCoeffs(const Matrix& data, template double MultFarFieldExpansion::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 *multiindex_mapping = + const std::vector& multiindex_mapping = sea_->get_multiindex_mapping(); - const ArrayList *lower_mapping_index = + const std::vector& 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::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 alpha_minus_beta_mapping; - alpha_minus_beta_mapping.Init(dim); + std::vector 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 &alpha_mapping = multiindex_mapping[alpha]; - const ArrayList &lower_mappings_for_alpha = + const std::vector& alpha_mapping = multiindex_mapping[alpha]; + const std::vector& 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 &beta_mapping = + const std::vector& 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::ConvolveField } template -void MultFarFieldExpansion::RefineCoeffs(const Matrix& data, - const Vector& weights, +void MultFarFieldExpansion::RefineCoeffs(const arma::mat& data, + const arma::vec& weights, int begin, int end, int order) { @@ -208,16 +202,17 @@ void MultFarFieldExpansion::RefineCoeffs(const Matrix& data, else { order_ = order; - coeffs_.SetZero(); + coeffs_.zeros(); AccumulateCoeffs(data, weights, begin, end, order); } } template -double MultFarFieldExpansion::EvaluateField(const Matrix& data, +double MultFarFieldExpansion::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 @@ -239,16 +234,14 @@ double MultFarFieldExpansion::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::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 &traversal_order = + const std::vector& 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 &mapping = sea_->get_multiindex(index); + const std::vector& mapping = sea_->get_multiindex(index); double arrtmp = ka_->ComputePartialDerivative(derivative_map, mapping); double prod = coeffs_[index] * arrtmp; @@ -283,19 +276,18 @@ double MultFarFieldExpansion::EvaluateField(const double *x_q, } template -void MultFarFieldExpansion::Init(const Vector& center, +void MultFarFieldExpansion::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 @@ -305,13 +297,11 @@ template 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::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::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 &mapping = sea_->get_multiindex(i); + const std::vector& mapping = sea_->get_multiindex(i); fprintf(stream, "%g ", coeffs_[i]); fprintf(stream, "(-1)^("); @@ -401,24 +391,19 @@ void MultFarFieldExpansion::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 *multiindex_mapping = + const std::vector >& multiindex_mapping = sea_->get_multiindex_mapping(); - const ArrayList *lower_mapping_index = + const std::vector >& lower_mapping_index = sea_->get_lower_mapping_index(); - ArrayList tmp_storage; - Vector center_diff; - Vector inv_multiindex_factorials; - - center_diff.Init(dim); + std::vector 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::TranslateFromFarField } // get the order of traversal for the given order of approximation - const ArrayList &traversal_order = + const std::vector& traversal_order = sea_->traversal_mapping_[order]; for(index_t j = 0; j < total_num_coeffs; j++) { int index = traversal_order[j]; - const ArrayList &gamma_mapping = multiindex_mapping[index]; - const ArrayList &lower_mappings_for_gamma = + const std::vector& gamma_mapping = multiindex_mapping[index]; + const std::vector& 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 &inner_mapping = - multiindex_mapping[lower_mappings_for_gamma[k]]; + const std::vector& inner_mapping = multiindex_mapping[lower_mappings_for_gamma[k]]; int flag = 0; double diff1; @@ -496,22 +480,18 @@ template void MultFarFieldExpansion::TranslateToLocal (MultLocalExpansion &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::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::TranslateToLocal } // compute required partial derivatives - ka_->ComputeDirectionalDerivatives(cent_diff, &derivative_map, + ka_->ComputeDirectionalDerivatives(cent_diff, derivative_map, 2 * truncation_order); - ArrayList beta_plus_alpha; - beta_plus_alpha.Init(dimension); + std::vector beta_plus_alpha; + beta_plus_alpha.reserve(dimension); // get the order of traversal for the given order of approximation - const ArrayList &traversal_order = + const std::vector &traversal_order = sea_->traversal_mapping_[truncation_order]; for(index_t j = 0; j < total_num_coeffs; j++) { int index = traversal_order[j]; - const ArrayList &beta_mapping = sea_->get_multiindex(index); + const std::vector &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 &alpha_mapping = + const std::vector &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::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]) * diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion.h index 133c05e64b..c2554dc81b 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion.h @@ -12,7 +12,8 @@ #ifndef MULT_LOCAL_EXPANSION #define MULT_LOCAL_EXPANSION -#include "fastlib/fastlib.h" +#include + #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 ¢er_; } - - const Vector* get_center() const { return ¢er_; } + 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); /** diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion_impl.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion_impl.h index 481f82f36c..5415abb1cd 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion_impl.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_local_expansion_impl.h @@ -6,8 +6,8 @@ #define MULT_LOCAL_EXPANSION_IMPL_H template -void MultLocalExpansion::AccumulateCoeffs(const Matrix& data, - const Vector& weights, +void MultLocalExpansion::AccumulateCoeffs(const arma::mat& data, + const arma::vec& weights, int begin, int end, int order) { @@ -19,23 +19,20 @@ void MultLocalExpansion::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 &traversal_order = + const std::vector &traversal_order = sea_->traversal_mapping_[order]; // for each data point, @@ -43,17 +40,17 @@ void MultLocalExpansion::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 &mapping = sea_->get_multiindex(index); + const std::vector& 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::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::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 mapping = sea_->get_multiindex(i); + std::vector mapping = sea_->get_multiindex(i); fprintf(stream, "%g", coeffs_[i]); for(index_t d = 0; d < dim; d++) { @@ -103,7 +100,7 @@ void MultLocalExpansion::PrintDebug(const char *name, } template -double MultLocalExpansion::EvaluateField(const Matrix& data, +double MultLocalExpansion::EvaluateField(const arma::mat& data, int row_num) const { // if there are no local coefficients, then return 0 @@ -124,16 +121,13 @@ double MultLocalExpansion::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 heads; - heads.Init(dim + 1); + arma::vec x_Q_to_x_q(dim); + arma::vec tmp(sea_->get_max_total_num_coeffs()); + std::vector 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::EvaluateField(const Matrix& data, tmp[0] = 1.0; // get the order of traversal for the given order of approximation - const ArrayList &traversal_order = + const std::vector& traversal_order = sea_->traversal_mapping_[order_]; for(index_t i = 1; i < total_num_coeffs; i++) { int index = traversal_order[i]; - const ArrayList &lower_mappings = + const std::vector &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 &mapping = sea_->multiindex_mapping_[index]; - const ArrayList &direct_ancestor_mapping = + const std::vector& mapping = sea_->multiindex_mapping_[index]; + const std::vector& 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::EvaluateField(const Matrix& data, } template -double MultLocalExpansion::EvaluateField(const Vector& x_q) const { +double MultLocalExpansion::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::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 heads; - heads.Init(dim + 1); + arma::vec x_Q_to_x_q(dim); + arma::vec tmp(sea_->get_max_total_num_coeffs()); + std::vector 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::EvaluateField(const Vector& x_q) const { tmp[0] = 1.0; // get the order of traversal for the given order of approximation - ArrayList &traversal_order = sea_->traversal_mapping_[order_]; + std::vector& traversal_order = sea_->traversal_mapping_[order_]; for(index_t i = 1; i < total_num_coeffs; i++) { int index = traversal_order[i]; - ArrayList &lower_mappings = sea_->lower_mapping_index_[index]; + std::vector &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 &mapping = sea_->multiindex_mapping_[index]; - const ArrayList &direct_ancestor_mapping = + const std::vector &mapping = sea_->multiindex_mapping_[index]; + const std::vector &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::EvaluateField(const Vector& x_q) const { } template -void MultLocalExpansion::Init(const Vector& center, +void MultLocalExpansion::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 @@ -269,13 +259,10 @@ void MultLocalExpansion::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 @@ -295,11 +282,11 @@ template void MultLocalExpansion::TranslateFromFarField (const MultFarFieldExpansion &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::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::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::TranslateFromFarField // compute required partial derivatives ka_->ComputeDirectionalDerivatives(cent_diff, &derivative_map, 2 * order_); - ArrayList beta_plus_alpha; - beta_plus_alpha.Init(dimension); + std::vector beta_plus_alpha; + beta_plus_alpha.reserve(dimension); // get the order of traversal for the given order of approximation - ArrayList &traversal_order = sea_->traversal_mapping_[far_order]; + std::vector &traversal_order = sea_->traversal_mapping_[far_order]; for(index_t j = 0; j < total_num_coeffs; j++) { int index_j = traversal_order[j]; - ArrayList beta_mapping = sea_->get_multiindex(index_j); + std::vector 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 alpha_mapping = sea_->get_multiindex(index_k); + std::vector 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::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::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 *upper_mapping_index = + const std::vector >& 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 tmp_storage; - tmp_storage.Init(dim); + std::vector 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::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 &traversal_order = + const std::vector &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 &alpha_mapping = sea_->get_multiindex(index_j); - const ArrayList &upper_mappings_for_alpha = + const std::vector &alpha_mapping = sea_->get_multiindex(index_j); + const std::vector &upper_mappings_for_alpha = upper_mapping_index[index_j]; double pos_coeffs = 0; double neg_coeffs = 0; @@ -438,7 +421,7 @@ void MultLocalExpansion::TranslateToLocal(MultLocalExpansion &se) { break; } - const ArrayList &beta_mapping = + const std::vector &beta_mapping = sea_->get_multiindex(upper_mappings_for_alpha[k]); int flag = 0; double diff1 = 1.0; diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_series_expansion_aux.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_series_expansion_aux.h index 806ed281af..72eb069341 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_series_expansion_aux.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/mult_series_expansion_aux.h @@ -4,7 +4,8 @@ #ifndef MULT_SERIES_EXPANSION_AUX_H #define MULT_SERIES_EXPANSION_AUX_H -#include "fastlib/fastlib.h" +#include +#include /** * Series expansion class for multiplicative kernel functions @@ -18,63 +19,48 @@ class MultSeriesExpansionAux { int max_order_; - Vector factorials_; + arma::vec factorials_; - ArrayList list_total_num_coeffs_; + std::vector 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 > multiindex_mapping_; + std::vector< std::vector > 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 > lower_mapping_index_; + std::vector< std::vector > 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 > upper_mapping_index_; + std::vector< std::vector > 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 > 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 > 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 &mapping = multiindex_mapping_[j]; + const std::vector& 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 diff; - diff.Init(dim_); + std::vector 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 &outer_mapping = multiindex_mapping_[i]; - lower_mapping_index_[i].Init(); + const std::vector& outer_mapping = multiindex_mapping_[i]; + lower_mapping_index_[i].clear(); for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) { - const ArrayList &inner_mapping = multiindex_mapping_[j]; + const std::vector& 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 &beta_mapping = multiindex_mapping_[j]; + const std::vector& beta_mapping = multiindex_mapping_[j]; for(index_t k = 0; k < list_total_num_coeffs_[limit]; k++) { // alpha mapping - const ArrayList &alpha_mapping = multiindex_mapping_[k]; + const std::vector& 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 diff; - diff.Init(dim_); + std::vector 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 &outer_mapping = multiindex_mapping_[i]; - upper_mapping_index_[i].Init(); + const std::vector& outer_mapping = multiindex_mapping_[i]; + upper_mapping_index_[i].clear(); for(index_t j = 0; j < list_total_num_coeffs_[limit]; j++) { - const ArrayList &inner_mapping = multiindex_mapping_[j]; + const std::vector& 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 * get_lower_mapping_index() const { - return lower_mapping_index_.begin(); + const std::vector >& get_lower_mapping_index() const { + return lower_mapping_index_; } int get_max_order() const { return max_order_; } - const ArrayList & get_multiindex(int pos) const { + const std::vector& get_multiindex(int pos) const { return multiindex_mapping_[pos]; } - const ArrayList * get_multiindex_mapping() const { - return multiindex_mapping_.begin(); + const std::vector >& 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 * get_upper_mapping_index() const { - return upper_mapping_index_.begin(); + const std::vector >& 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 &multiindex) const { + int ComputeMultiindexPosition(const std::vector& 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); } } diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.cc b/fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.cc index bb5a78f266..b02ef8b3cb 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.cc +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.cc @@ -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 >& 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 >& 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); } } diff --git a/fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.h b/fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.h index b14e23cc6a..b1c5a8b1d6 100644 --- a/fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.h +++ b/fastlib/branches/fastlib-stl/mlpack/series_expansion/series_expansion_aux.h @@ -4,7 +4,8 @@ #ifndef SERIES_EXPANSION_AUX #define SERIES_EXPANSION_AUX -#include "fastlib/fastlib.h" +#include +#include /** * Series expansion class. @@ -17,15 +18,15 @@ class SeriesExpansionAux { int max_order_; - Vector factorials_; + arma::vec factorials_; std::vector 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 > multiindex_mapping_; @@ -44,29 +45,15 @@ class SeriesExpansionAux { std::vector< std::vector > 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 &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 > & 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 > & 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 >& get_upper_mapping_index() const; // interesting functions