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

Going the wrapper approach in dense point to use armadillo operations.

Browse Source
This commit is contained in:
Dongryeol Lee
2010-12-22 04:22:55 +00:00
parent 78fd08868a
commit 766f545a67
10 changed files with 107 additions and 282 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/trunk/contrib/dongryel/thesis_research/core/math/CMakeLists.txt
-1
View File
@@ -6,7 +6,6 @@ set(SOURCES
global.cc
math_lib.h
math_lib_impl.h
linear_algebra.h
range.h
)
fastlib/trunk/contrib/dongryel/thesis_research/core/math/linear_algebra.h
-225
View File
@@ -1,225 +0,0 @@
/** @file linear_algebra.h
*
* @author Dongryeol Lee (dongryel@cc.gatech.edu)
*/
#ifndef CORE_MATH_LINEAR_ALGEBRA_H
#define CORE_MATH_LINEAR_ALGEBRA_H
#include <armadillo>
namespace core {
namespace table {
class DensePoint;
class DenseMatrix;
};
};
namespace core {
namespace math {
template<typename MatrixType>
static void MatrixTripleProduct(
const MatrixType &left, const MatrixType &mid,
core::table::DenseMatrix *product) {
product->Init(left.n_rows(), left.n_cols());
// Use armadillo matrices to compute the triple product. This makes
// a light copy, so there is very little performance lost.
arma::mat left_alias(left.ptr(), left.n_rows(), left.n_cols());
arma::mat mid_alias(mid.ptr(), mid.n_rows(), mid.n_cols());
arma::mat product_alias(
product->ptr(), product->n_rows(), product->n_cols(), false);
product_alias = left_alias * mid_alias * arma::trans(left_alias);
}
template<typename MatrixType>
static void MatrixTripleProduct(
const MatrixType &left, const MatrixType &mid, const MatrixType &right,
core::table::DenseMatrix *product) {
product->Init(left.n_rows(), right.n_cols());
// Use armadillo matrices to compute the triple product. This makes
// a light copy, so there is very little performance lost.
arma::mat left_alias(left.ptr(), left.n_rows(), left.n_cols());
arma::mat mid_alias(mid.ptr(), mid.n_rows(), mid.n_cols());
arma::mat right_alias(right.ptr(), right.n_rows(), right.n_cols());
arma::mat product_alias(
product->ptr(), product->n_rows(), product->n_cols(), false);
product_alias = left_alias * mid_alias * right_alias;
}
template<typename VectorType>
static double Dot(const VectorType &a, const VectorType &b) {
arma::vec a_mat(a.ptr(), a.length());
arma::vec b_mat(b.ptr(), b.length());
return arma::dot(a_mat, b_mat);
}
template<typename VectorType>
static double LengthEuclidean(const VectorType &a) {
return sqrt(core::math::Dot(a, a));
}
template<typename VectorType>
static void CopyValues(
const VectorType &vec_in, VectorType *vec_out) {
for(unsigned int i = 0; i < vec_in.n_elem; i++) {
(*vec_out)[i] = vec_in[i];
}
}
template<typename VectorType>
static void SubFrom(
const VectorType &vec_in, VectorType *vec_out) {
for(unsigned int i = 0; i < vec_in.n_elem; i++) {
(*vec_out)[i] -= vec_in[i];
}
}
template<typename VectorType>
static void SubOverwrite(
const VectorType &sub, const VectorType &sub_from, VectorType *vec_out) {
for(int i = 0; i < sub.length(); i++) {
(*vec_out)[i] = sub_from[i] - sub[i];
}
}
template<typename VectorType>
static void SubInit(
const VectorType &sub, const VectorType &sub_from, VectorType *vec_out) {
vec_out->Init(sub.length());
SubOverwrite(sub, sub_from, vec_out);
}
template<typename T>
class AddExpertTrait {
public:
static void Compute(double scale, const T &vec_scaled, T *vec_add_to);
};
template<>
class AddExpertTrait<core::table::DenseMatrix> {
public:
static void Compute(
double scale, const core::table::DenseMatrix &mat_scaled,
core::table::DenseMatrix *mat_add_to) {
arma::mat mat_scaled_alias(
mat_scaled.ptr(), mat_scaled.n_rows(), mat_scaled.n_cols());
arma::mat mat_add_to_alias(
mat_add_to->ptr(), mat_add_to->n_rows(), mat_add_to->n_cols(), false);
mat_add_to_alias = mat_add_to_alias + scale * mat_scaled_alias;
}
};
template<>
class AddExpertTrait<core::table::DensePoint> {
public:
static void Compute(
double scale, const core::table::DensePoint &vec_scaled,
core::table::DensePoint *vec_add_to) {
arma::vec vec_scaled_alias(
vec_scaled.ptr(), vec_scaled.length());
arma::vec vec_add_to_alias(
vec_add_to->ptr(), vec_add_to->length(), false);
vec_add_to_alias = vec_add_to_alias + scale * vec_scaled_alias;
}
};
template<typename T>
static void AddExpert(
double scale, const T &vec_scaled, T *vec_add_to) {
core::math::AddExpertTrait<T>::Compute(scale, vec_scaled, vec_add_to);
}
template<typename VectorType>
static void AddTo(
const VectorType &vec_in, VectorType *vec_out) {
arma::vec vec_in_alias(vec_in.ptr(), vec_in.length());
arma::vec vec_out_alias(vec_out->ptr(), vec_in.length(), false);
vec_out_alias = vec_out_alias + vec_in_alias;
}
/** @brief Computes $c = c + \alpha * a b^T$.
*/
template<typename VectorType, typename MatrixType>
static void MulExpert(
double alpha, const VectorType &a, const VectorType &b, MatrixType *c) {
arma::vec a_alias(a.ptr(), a.length());
arma::vec b_alias(b.ptr(), b.length());
arma::mat c_alias(c->ptr(), c->n_rows(), c->n_cols(), false);
c_alias = c_alias + alpha * a_alias * arma::trans(b_alias);
}
template<typename MatrixType, typename VectorType>
static void MulInit(
const MatrixType &a, const VectorType &b, VectorType *c) {
arma::mat a_alias(a.ptr(), a.n_rows(), a.n_cols());
arma::vec b_alias(b.ptr(), b.length());
c->Init(a.n_rows());
arma::vec c_alias(c->ptr(), c.length(), false);
c_alias = a_alias * b_alias;
}
template<typename VectorType>
class ScaleTrait {
public:
static void Compute(double scale, VectorType *vec);
};
template<>
class ScaleTrait<core::table::DensePoint> {
public:
static void Compute(double scale, core::table::DensePoint *vec) {
for(int i = 0; i < vec->length(); i++) {
(*vec)[i] *= scale;
}
}
};
template<>
class ScaleTrait<core::table::DenseMatrix> {
public:
static void Compute(double scale, core::table::DenseMatrix *vec) {
for(int j = 0; j < vec->n_cols(); j++) {
for(int i = 0; i < vec->n_rows(); i++) {
vec->set(i, j, vec->get(i, j) * scale);
}
}
}
};
template<typename VectorType>
static void Scale(double scale, VectorType *vec) {
ScaleTrait<VectorType>::Compute(scale, vec);
}
template<typename VectorType>
static void ScaleOverwrite(
double scale, const VectorType &vec_in, VectorType *vec_out) {
arma::vec vec_in_alias(vec_in.ptr(), vec_in.length());
arma::vec vec_out_alias(vec_out->ptr(), vec_in.length(), false);
for(unsigned int i = 0; i < vec_in_alias.n_elem; i++) {
vec_out_alias[i] = vec_in_alias[i] * scale;
}
}
template<typename VectorType>
static void ScaleInit(
double scale, const VectorType &vec_in, VectorType *vec_out) {
vec_out->set_size(vec_in.n_elem);
ScaleOverwrite(scale, vec_in, vec_out);
}
};
};
#endif
fastlib/trunk/contrib/dongryel/thesis_research/core/monte_carlo/mean_variance_pair_matrix.h
+6 -5
View File
@@ -6,6 +6,7 @@
#ifndef CORE_MONTE_CARLO_MEAN_VARIANCE_PAIR_MATRIX_H
#define CORE_MONTE_CARLO_MEAN_VARIANCE_PAIR_MATRIX_H
#include <armadillo>
#include "core/monte_carlo/mean_variance_pair.h"
namespace core {
@@ -55,8 +56,8 @@ class MeanVariancePairVector {
new core::monte_carlo::MeanVariancePair[n_elements_];
}
void sample_means(core::table::DensePoint *point_out) const {
point_out->Init(n_elements_);
void sample_means(arma::vec *point_out) const {
point_out->set_size(n_elements_);
for(int i = 0; i < n_elements_; i++) {
(*point_out)[i] = ptr_[i].sample_mean();
}
@@ -110,11 +111,11 @@ class MeanVariancePairMatrix {
}
}
void sample_means(core::table::DenseMatrix *point_out) const {
point_out->Init(n_rows_, n_cols_);
void sample_means(arma::mat *point_out) const {
point_out->set_size(n_rows_, n_cols_);
for(int j = 0; j < n_cols_; j++) {
for(int i = 0; i < n_rows_; i++) {
point_out->set(i, j, this->get(i, j).sample_mean());
point_out->at(i, j) = this->get(i, j).sample_mean();
}
}
}
fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h
+7 -1
View File
@@ -14,11 +14,16 @@ namespace core {
namespace optimization {
template<typename FunctionType>
class TrustRegion {
public:
enum TrustRegionSearchMethod {CAUCHY, DOGLEG, STEIHAUG};
private:
double max_radius_;
FunctionType *function_;
TrustRegionSearchMethod search_method_;
private:
void Evaluate_(const arma::vec &iterate);
@@ -53,7 +58,8 @@ class TrustRegion {
max_radius_ = max_radius_in;
}
void Init(FunctionType &function_in);
void Init(
FunctionType &function_in, TrustRegionSearchMethod search_method_in);
void Optimize(int num_iterations, arma::vec *iterate);
};
fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h
+11 -1
View File
@@ -16,6 +16,12 @@ template<typename FunctionType>
TrustRegion<FunctionType>::TrustRegion() {
max_radius_ = 10.0;
function_ = NULL;
search_method_ = CAUCHY;
}
template<typename FunctionType>
double TrustRegion<FunctionType>::ReductionRatio_(const arma::vec &iterate, const arma::vec &step,) {
}
template<typename FunctionType>
@@ -266,12 +272,16 @@ template<typename FunctionType>
void TrustRegion<FunctionType>::ObtainStepDirection_(
arma::vec *step_direction, double *step_direction_norm) {
switch(search_method_) {
}
}
template<typename FunctionType>
void TrustRegion<FunctionType>::Init(FunctionType &function_in) {
void TrustRegion<FunctionType>::Init(
FunctionType &function_in, TrustRegionSearchMethod search_method_in) {
function_ = &function_in;
search_method_ = search_method_in;
}
template<typename FunctionType>
fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h
+11
View File
@@ -8,6 +8,7 @@
#ifndef CORE_TABLE_DENSE_POINT_H
#define CORE_TABLE_DENSE_POINT_H
#include <armadillo>
#include <vector>
#include <boost/interprocess/offset_ptr.hpp>
#include <boost/serialization/serialization.hpp>
@@ -175,6 +176,16 @@ class DensePoint {
printf("\n");
}
};
template<typename DensePointType>
static void DensePointToArmaVec(
const DensePointType &point_in, arma::vec *vec_out) {
const_cast<arma::u32 &>(vec_out->n_rows) = point_in.length();
const_cast<arma::u32 &>(vec_out->n_cols) = 1;
const_cast<arma::u32 &>(vec_out->n_elem) = point_in.length();
const_cast<bool &>(vec_out->use_aux_mem) = true;
const_cast<double *&>(vec_out->mem) = const_cast<double *>(point_in.ptr());
}
};
};
fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h
+9 -5
View File
@@ -152,7 +152,9 @@ class DistributedTable: public boost::noncopyable {
void ReplenishNodes_(std::vector<TreeType *> &top_leaf_nodes) {
core::table::DensePoint tmp_point;
arma::vec tmp_point_alias;
tmp_point.Init(top_leaf_nodes[0]->bound().center().length());
core::table::DensePointToArmaVec(tmp_point, &tmp_point_alias);
int num_additional = table_outbox_group_comm_size_ -
top_leaf_nodes.size();
int num_samples = std::max(1, core::math::RandInt(top_leaf_nodes.size()));
@@ -160,15 +162,17 @@ class DistributedTable: public boost::noncopyable {
// Randomly add new dummy nodes with randomly chosen centroids
// averaged.
for(int j = 0; j < num_additional; j++) {
tmp_point.SetZero();
tmp_point_alias.zeros();
for(int i = 0; i < num_samples; i++) {
tmp_point.SetZero();
core::math::AddTo(
arma::vec random_node_center;
core::table::DensePointToArmaVec(
top_leaf_nodes[
core::math::RandInt(top_leaf_nodes.size())]->bound().center(),
&tmp_point);
&random_node_center);
tmp_point_alias += random_node_center;
}
core::math::Scale(1.0 / static_cast<double>(num_samples), &tmp_point);
tmp_point_alias = (1.0 / static_cast<double>(num_samples)) *
tmp_point_alias;
top_leaf_nodes.push_back(new TreeType());
top_leaf_nodes[ top_leaf_nodes.size() - 1 ]->bound().center().Copy(
tmp_point);
fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h
+14 -7
View File
@@ -6,9 +6,9 @@
#ifndef CORE_TREE_DISTRIBUTED_LOCAL_KMEANS_H
#define CORE_TREE_DISTRIBUTED_LOCAL_KMEANS_H
#include <armadillo>
#include <boost/mpi.hpp>
#include <boost/serialization/string.hpp>
#include "core/math/linear_algebra.h"
#include "core/table/dense_point.h"
namespace core {
@@ -25,12 +25,15 @@ class DistributedLocalKMeans {
core::table::DensePoint centroid_;
arma::vec centroid_alias_;
int num_points_;
public:
void Init(int length_in) {
centroid_.Init(length_in);
core::table::DensePointToArmaVec(centroid_, &centroid_alias_);
}
const core::table::DensePoint &centroid() const {
@@ -53,22 +56,26 @@ class DistributedLocalKMeans {
double factor =
static_cast<double>(num_points_) /
static_cast<double>(num_points_ + centroid_in.num_points());
core::math::Scale(factor, &centroid_);
core::math::AddExpert(
1.0 - factor, centroid_in.centroid(), &centroid_);
arma::vec centroid_in_alias;
core::table::DensePointToArmaVec(
centroid_in.centroid(), &centroid_in_alias);
centroid_alias_ = factor * centroid_alias_ +
(1.0 - factor) * centroid_in_alias;
num_points_ = num_points_ + centroid_in.num_points();
}
void Add(const core::table::DensePoint &point_in) {
double factor = static_cast<double>(num_points_) /
static_cast<double>(num_points_ + 1);
core::math::Scale(factor, &centroid_);
core::math::AddExpert(1.0 - factor, point_in, &centroid_);
arma::vec point_in_alias;
core::table::DensePointToArmaVec(point_in, &point_in_alias);
centroid_alias_ = factor * centroid_alias_;
centroid_alias_ += (1.0 - factor) * point_in_alias;
num_points_++;
}
void Reset() {
centroid_.SetZero();
centroid_alias_.zeros();
num_points_ = 0;
}
fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h
+17 -9
View File
@@ -9,7 +9,6 @@
#include <vector>
#include "ball_bound.h"
#include "general_spacetree.h"
#include "core/math/linear_algebra.h"
#include "core/metric_kernels/abstract_metric.h"
#include "core/table/dense_matrix.h"
#include "core/table/memory_mapped_file.h"
@@ -69,11 +68,15 @@ class GenMetricTree {
int end = begin + count;
core::table::DensePoint col_point;
arma::vec bound_ref;
core::table::DensePointToArmaVec(bounds->center(), &bound_ref);
for(int i = begin; i < end; i++) {
matrix.MakeColumnVector(i, &col_point);
core::math::AddTo(col_point, &(bounds->center()));
arma::vec col_point_ref;
core::table::DensePointToArmaVec(col_point, &col_point_ref);
bound_ref += col_point_ref;
}
core::math::Scale(1.0 / static_cast<double>(count), &(bounds->center()));
bound_ref = (1.0 / static_cast<double>(count)) * bound_ref;
double furthest_distance;
FurthestColumnIndex_(
@@ -88,12 +91,17 @@ class GenMetricTree {
TreeType *node, TreeType *left, TreeType *right) {
// Compute the weighted sum of the two pivots
node->bound().center().CopyValues(left->bound().center());
core::math::Scale(left->count(), &(node->bound().center()));
core::math::AddExpert(
right->count(), right->bound().center(), & (node->bound().center()));
core::math::Scale(
1.0 / static_cast<double>(node->count()), & (node->bound().center()));
arma::vec bound_ref;
core::table::DensePointToArmaVec(node->bound().center(), &bound_ref);
arma::vec left_bound_ref;
core::table::DensePointToArmaVec(left->bound().center(), &left_bound_ref);
arma::vec right_bound_ref;
core::table::DensePointToArmaVec(
right->bound().center(), &right_bound_ref);
bound_ref = left->count() * left_bound_ref +
right->count() * right_bound_ref;
bound_ref =
(1.0 / static_cast<double>(node->count())) * bound_ref;
double left_max_dist, right_max_dist;
FurthestColumnIndex_(
fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h
+32 -28
View File
@@ -6,10 +6,10 @@
#ifndef MLPACK_MIXED_LOGIT_DCM_DCM_TABLE_H
#define MLPACK_MIXED_LOGIT_DCM_DCM_TABLE_H
#include <armadillo>
#include <algorithm>
#include <vector>
#include "core/table/table.h"
#include "core/math/linear_algebra.h"
#include "core/monte_carlo/mean_variance_pair.h"
#include "core/monte_carlo/mean_variance_pair_matrix.h"
#include "mlpack/mixed_logit_dcm/mixed_logit_dcm_distribution.h"
@@ -92,8 +92,8 @@ class DCMTable {
* form.
*/
void ComputeChoiceProbabilities_(
int person_index, const core::table::DensePoint &parameter_vector,
core::table::DensePoint *choice_probabilities) {
int person_index, const arma::vec &parameter_vector,
arma::vec *choice_probabilities) {
int num_discrete_choices = this->num_discrete_choices(person_index);
choice_probabilities->Init(num_discrete_choices);
@@ -109,8 +109,12 @@ class DCMTable {
core::table::DensePoint attribute_for_discrete_choice;
this->get_attribute_vector(
person_index, discrete_choice_index, &attribute_for_discrete_choice);
double dot_product = core::math::Dot(
parameter_vector, attribute_for_discrete_choice);
arma::vec attribute_for_discrete_choice_alias;
core::table::DensePointToArmaVec(
attribute_for_discrete_choice, attribute_for_discrete_choice_alias);
double dot_product =
arma::dot(
parameter_vector, attribute_for_discrete_choice_alias);
double unnormalized_probability = exp(dot_product);
normalizing_sum += unnormalized_probability;
(*choice_probabilities)[discrete_choice_index] =
@@ -152,11 +156,11 @@ class DCMTable {
* computation of Equation 8.14 in the paper.
*/
void SimulatedLoglikelihoodHessian(
core::table::DenseMatrix *likelihood_hessian) const {
arma::mat *likelihood_hessian) const {
likelihood_hessian->Init(
likelihood_hessian->set_size(
distribution_->num_parameters(), distribution_->num_parameters());
likelihood_hessian->SetZero();
likelihood_hessian->zeros();
// For each active person,
for(int i = 0; i < num_active_people_; i++) {
@@ -178,23 +182,23 @@ class DCMTable {
simulated_loglikelihood_hessians_[person_index].first;
const core::monte_carlo::MeanVariancePairVector &hessian_second_part =
simulated_loglikelihood_hessians_[person_index].second;
core::table::DenseMatrix hessian_first;
core::table::DensePoint hessian_second;
arma::mat hessian_first;
arma::vec hessian_second;
hessian_first_part.sample_means(&hessian_first);
hessian_second_part.sample_means(&hessian_second);
// Construct the contribution on the fly.
core::math::AddExpert(
inverse_simulated_choice_probability,
hessian_first, likelihood_hessian);
core::math::MulExpert(
- core::math::Sqr(inverse_simulated_choice_probability),
hessian_second, hessian_second, likelihood_hessian);
(*likelihood_hessian) += inverse_simulated_choice_probability *
hessian_first;
(*likelihood_hessian) +=
(- core::math::Sqr(inverse_simulated_choice_probability)) *
hessian_second * arma::trans(hessian_second);
}
// Divide by the number of people.
core::math::Scale(
1.0 / static_cast<double>(num_active_people_), likelihood_hessian);
(*likelihood_hessian) =
(1.0 / static_cast<double>(num_active_people_)) *
(*likelihood_hessian);
}
/** @brief Return the gradient of the current simulated log
@@ -202,10 +206,10 @@ class DCMTable {
* in the paper.
*/
void SimulatedLoglikelihoodGradient(
core::table::DensePoint *likelihood_gradient) const {
arma::vec *likelihood_gradient) const {
likelihood_gradient->Init(distribution_->num_parameters());
likelihood_gradient->SetZero();
likelihood_gradient->set_size(distribution_->num_parameters());
likelihood_gradient->zeros();
// For each active person,
for(int i = 0; i < num_active_people_; i++) {
@@ -230,14 +234,14 @@ class DCMTable {
// Add the inverse probability weighted gradient vector for
// the current person to the total tally.
core::math::AddExpert(
inverse_simulated_choice_probability,
gradient_vector, likelihood_gradient);
(*likelihood_gradient) +=
inverse_simulated_choice_probability * gradient_vector;
}
// Divide by the number of people.
core::math::Scale(
1.0 / static_cast<double>(num_active_people_), likelihood_gradient);
(*likelihood_gradient) =
(1.0 / static_cast<double>(num_active_people_)) *
(*likelihood_gradient);
}
/** @brief Return the current simulated log likelihood score.
@@ -378,8 +382,8 @@ class DCMTable {
// Simulated log-likelihood gradient update by the simulated
// choice probabilty scaled gradient product.
core::math::Scale(
choice_probabilities[discrete_choice_index], &beta_gradient_product);
beta_gradient_product = choice_probabilities[discrete_choice_index] *
beta_gradient_product;
simulated_loglikelihood_gradients_[person_index].push_back(
beta_gradient_product);