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

Code upgrade to compute force vectors due to Axilrod-Teller potential in progress.

Browse Source
This commit is contained in:
Dongryeol Lee
2008-06-26 19:30:05 +00:00
parent 9829df2f2c
commit d9fdc968a8
8 changed files with 1117 additions and 910 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib2/contrib/dongryel/multibody/axilrod_teller_force.tex
+43 -5
View File
@@ -97,11 +97,49 @@ contribution on each $r_i \in R_I$ due to $R_J$ and $R_K$ is:
\begin{align*}
& \sum\limits_{r_j \in R_J} \sum\limits_{r_k \in R_K} F(r_i, \{ (r_j,
r_k) \}) \\
=& \sum\limits_{r_j \in R_J} \sum\limits_{r_k \in R_K} -(r_i[x] -
r_j[x]) V(r_i, r_j) - (r_i[x] - r_k[x]) V(r_i, r_k)\\
=& \sum\limits_{r_j \in R_J} \sum\limits_{r_k \in R_K} -(r_i[x] -
r_j[x]) V(r_i, r_j) - \sum\limits_{r_j \in R_J} \sum\limits_{r_k \in
R_K} (r_i[x] - r_k[x]) V(r_i, r_k)\\
=& \sum\limits_{r_j \in R_J} \sum\limits_{r_k \in R_K} -(r_i -
r_j) V(r_i, r_j) - (r_i - r_k) V(r_i, r_k)\\
=& \sum\limits_{r_j \in R_J} \sum\limits_{r_k \in R_K} -(r_i -
r_j) V(r_i, r_j) - \sum\limits_{r_j \in R_J} \sum\limits_{r_k \in
R_K} (r_i - r_k) V(r_i, r_k)\\
=&- r_i \left( \sum\limits_{r_j \in R_J} \sum\limits_{r_k \in R_K} V(r_i, r_j)
+ V(r_i, r_k) \right) +
\sum\limits_{r_k \in R_K} \left( \sum\limits_{r_j \in R_J} r_j \right)
V(r_i, r_j) +\\
& \sum\limits_{r_j \in R_J} \left( \sum\limits_{r_k \in R_K} r_k \right)
V(r_i, r_k)
\end{align*}
Similarly, the total force contribution on each $r_j \in R_J$ due to
$R_I$ and $R_K$ is:
\begin{align*}
& \sum\limits_{r_i \in R_I} \sum\limits_{r_k \in R_K} F(r_j, \{ (r_i,
r_k) \})\\
=& \sum\limits_{r_i \in R_I} \sum\limits_{r_k \in R_K} (r_i - r_j)
V(r_i, r_j) - (r_j - r_k) V(r_j, r_k)\\
=& \sum\limits_{r_i \in R_I} \sum\limits_{r_k \in R_K} (r_i - r_j)
V(r_i, r_j) - \sum\limits_{r_i \in R_I} \sum\limits_{r_k \in R_K} (r_j
- r_k) V(r_j, r_k)\\
=&- r_j \left( \sum\limits_{r_i \in R_I} \sum\limits_{r_k \in R_K}
V(r_i, r_j) + V(r_j, r_k) \right) + \sum\limits_{r_k \in R_K} \left(
\sum\limits_{r_i \in R_I} r_i \right) V(r_i, r_j) +\\
& \sum\limits_{r_i \in R_I} \left( \sum\limits_{r_k \in R_K} r_k \right)
V(r_j, r_k)
\end{align*}
Finally, the total force contribution on each $r_k \in R_K$ due to
$R_I$ and $R_J$ is:
\begin{align*}
& \sum\limits_{r_i \in R_I} \sum\limits_{r_j \in R_J} F(r_k, \{ (r_i,
r_j) \})\\
=& \sum\limits_{r_i \in R_I} \sum\limits_{r_j \in R_J} (r_i - r_k)
V(r_i, r_k) + (r_j - r_k) V(r_j, r_k)\\
=& \sum\limits_{r_i \in R_I} \sum\limits_{r_j \in R_J} (r_i - r_k)
V(r_i, r_k) + \sum\limits_{r_i \in R_I} \sum\limits_{r_j \in R_J} (r_j
- r_k) V(r_j, r_k)\\
=&- r_k \left( \sum\limits_{r_i \in R_I} \sum\limits_{r_j \in R_J}
V(r_i, r_k) + V(r_j, r_k) \right) + \sum\limits_{r_j \in R_J} \left(
\sum\limits_{r_i \in R_I} r_i \right) V(r_i, r_k) +\\
& \sum\limits_{r_i \in R_I} \left( \sum\limits_{r_j \in R_J} r_j \right)
V(r_j, r_k)
\end{align*}
\end{document}
fastlib2/contrib/dongryel/multibody/build.py
+2
View File
@@ -3,7 +3,9 @@ librule(
name = "multibody", # this line can be safely omitted
sources = [], # files that must be compiled
headers = ["multibody.h",
"multibody_impl.h",
"multibody_kernel.h",
"multibody_stat.h",
"naive_multibody.h"], # include files part of the 'lib'
deplibs = ["mlpack/series_expansion:series_expansion",
"fastlib:fastlib_int"] # dependency
fastlib2/contrib/dongryel/multibody/multibody.h
+165 -648
View File
@@ -9,738 +9,255 @@
#include "mlpack/series_expansion/series_expansion_aux.h"
#include "multibody_kernel.h"
template<typename TKernelAux>
class MultibodyStat {
public:
/**
* Far field expansion created by the reference points in this node.
*/
FarFieldExpansion<TKernelAux> farfield_expansion_;
/**
* Local expansion stored in this node.
*/
LocalExpansion<TKernelAux> local_expansion_;
// getters and setters
FarFieldExpansion<TKernelAux> &get_farfield_coeffs() {
return farfield_expansion_;
}
/** Initialize the statistics */
void Init() {
}
void Init(const TKernelAux &ka) {
farfield_expansion_.Init(ka);
local_expansion_.Init(ka);
}
void Init(const Matrix& dataset, index_t &start, index_t &count) {
Init();
}
void Init(const Matrix& dataset, index_t &start, index_t &count,
const MultibodyStat& left_stat,
const MultibodyStat& right_stat) {
Init();
}
void Init(const Vector& center, const TKernelAux &ka) {
farfield_expansion_.Init(center, ka);
local_expansion_.Init(center, ka);
}
MultibodyStat() { }
~MultibodyStat() {}
};
#define INSIDE_MULTIBODY_H
#include "multibody_stat.h"
template<typename TMultibodyKernel, typename TKernelAux>
template<typename TMultibodyKernel>
class MultitreeMultibody {
FORBID_ACCIDENTAL_COPIES(MultitreeMultibody);
public:
public:
typedef BinarySpaceTree<DHrectBound<2>, Matrix,
MultibodyStat<TKernelAux> > Tree;
typedef BinarySpaceTree<DHrectBound<2>, Matrix, MultibodyStat > Tree;
typedef TMultibodyKernel MultibodyKernel;
// constructor/destructor
MultitreeMultibody() {}
////////// Constructor/Destructor //////////
/** @brief The default constructor.
*/
MultitreeMultibody() {
}
/** @brief The default destructor.
*/
~MultitreeMultibody() {
delete root_;
}
// getters/setters
const Matrix &get_data() const { return data_; }
////////// User-level Functions //////////
// interesting functions...
/** Main computation */
void Compute(double tau) {
/** @brief The main computation procedure.
*/
void Compute(double relative_error) {
ArrayList<Tree *> root_nodes;
root_nodes.Init(mkernel_.order());
// store node pointers
// Set the root node pointers for starting the computation.
for(index_t i = 0; i < mkernel_.order(); i++) {
root_nodes[i] = root_;
}
total_num_tuples_ = math::BinomialCoefficient(data_.n_cols(),
total_num_tuples_ = math::BinomialCoefficient(data_.n_cols(),
mkernel_.order());
tau_ = tau;
relative_error_ = relative_error;
// run and do timing for multitree multibody
NumPrunes_ = 0;
NumNodesExpanded_ = 0;
// Initialize intermediate computation spaces to zero.
negative_force1_e_.SetZero();
negative_force1_u_.SetZero();
positive_force1_l_.SetZero();
positive_force1_e_.SetZero();
negative_force2_e_.SetZero();
negative_force2_u_.SetZero();
positive_force2_l_.SetZero();
positive_force2_e_.SetZero();
total_force_e_.SetZero();
// Run and do timing for multitree multibody
MTMultibody(root_nodes, total_num_tuples_);
printf("Negative potential %g\n", neg_potential_e_);
printf("Positive potential %g\n", pos_potential_e_);
printf("Total potential estimate: %g\n", pos_potential_e_ +
neg_potential_e_);
printf("Number of series approximations: %d\n", NumPrunes_);
printf("Number of tuples of nodes expanded: %d\n", NumNodesExpanded_);
PostProcess(root_);
}
void InitExpansionObjects(Tree *node) {
if(node != NULL) {
node->stat().Init(ka_);
node->bound().CalculateMidpoint
(node->stat().farfield_expansion_.get_center());
node->bound().CalculateMidpoint
(node->stat().local_expansion_.get_center());
}
if(!node->is_leaf()) {
InitExpansionObjects(node->left());
InitExpansionObjects(node->right());
}
}
/** Initialize the kernel object, and build the tree */
/** @brief Initialize the kernel object, and build the tree.
*/
void Init(double bandwidth) {
const char *fname = fx_param_str(NULL, "data", NULL);
int leaflen = fx_param_int(NULL, "leaflen", 20);
// read in the dataset and build a kd-tree
// Read in the dataset and build a kd-tree.
fx_timer_start(NULL, "tree_d");
Dataset dataset_;
dataset_.InitFromFile(fname);
data_.Own(&(dataset_.matrix()));
root_ = tree::MakeKdTreeMidpoint<Tree>(data_, leaflen, NULL);
weights_.Init(data_.n_cols());
// by default, each point has a uniform weight
weights_.SetAll(1);
// set the maximum order of approximation here!
ka_.Init(bandwidth, 4, data_.n_rows());
// initialize the multibody kernel and the series expansion objects
// for all nodes
// Initialize the multibody kernel.
mkernel_.Init(bandwidth);
InitExpansionObjects(root_);
fx_timer_stop(NULL, "tree_d");
// more temporary variables initialization
// More temporary variables initialization.
non_leaf_indices_.Init(mkernel_.order());
distmat_.Init(mkernel_.order(), mkernel_.order());
exhaustive_indices_.Init(mkernel_.order());
node_bounds_.Init(mkernel_.order());
// initialize the combination generator
combination_.Init(mkernel_.order());
for(index_t i = 0; i < mkernel_.order(); i++) {
combination_[i] = i;
// Initialize space for computation values.
negative_force1_e_.Init(data_.n_cols());
negative_force1_u_.Init(data_.n_cols());
positive_force1_l_.Init(data_.n_cols());
positive_force1_e_.Init(data_.n_cols());
negative_force2_e_.Init(data_.n_rows(), data_.n_cols());
negative_force2_u_.Init(data_.n_rows(), data_.n_cols());
positive_force2_l_.Init(data_.n_rows(), data_.n_cols());
positive_force2_e_.Init(data_.n_rows(), data_.n_cols());
total_force_e_.Init(data_.n_rows(), data_.n_cols());
}
/** @brief Outputs the force vectors to the file.
*/
void PrintDebug() {
FILE *stream = fopen("force_vectors.txt", "w+");
for(index_t q = 0; q < data_.n_cols(); q++) {
for(index_t d = 0; d < data_.n_rows(); d++) {
fprintf(stream, "%g ", total_force_e_.get(d, q));
}
fprintf(stream, "\n");
}
combination_rank_ = 0;
// potential bounds and tokens initialized to 0.
neg_potential_u_ = neg_potential_e_ = 0;
pos_potential_l_ = pos_potential_e_ = 0;
extra_token_ = 0;
}
private:
// member variables
////////// Private Member Variables //////////
/** hrect bounds passed to evaluation */
ArrayList<DHrectBound<2> *> node_bounds_;
/** @brief The total number of n-tuples.
*/
double total_num_tuples_;
/** @brief The multibody kernel function.
*/
MultibodyKernel mkernel_;
/** @brief The accuracy requirement: componentwise relative error
* bound.
*/
double relative_error_;
/** The current list of non-leaf indices */
ArrayList<int> non_leaf_indices_;
/** Temporary space for storing indices selected for exhaustive computation
/** @brief The temporary space for storing indices selected for
* exhaustive computation.
*/
ArrayList<int> exhaustive_indices_;
/** Temporary space for storing pairwise distances */
/** @brief The temporary space for storing pairwise distances.
*/
Matrix distmat_;
/** pointer to the root of the tree */
/** @brief The pointer to the root of the tree.
*/
Tree *root_;
/** dataset for the tree */
/** @brief The dataset for the tree.
*/
Matrix data_;
/** weight for each point */
Vector weights_;
/** series approximation auxiliary computations */
TKernelAux ka_;
/** multibody kernel function */
MultibodyKernel mkernel_;
/** the total number of n-tuples to consider */
double total_num_tuples_;
/** Extra amount of error that can be spent */
double extra_token_;
/** negative potential estimate */
double neg_potential_e_;
/** running lower bound on the negative potential */
double neg_potential_u_;
/** positive potential estimate */
double pos_potential_e_;
/** running lower bound on the positive potential */
double pos_potential_l_;
/** approximation relative error bound */
double tau_;
/** number of prunes made */
int NumPrunes_;
/** number of nodes expanded */
int NumNodesExpanded_;
/** index enumerating a combination from beginning to the end */
Vector combination_;
/** rank of the current combination */
int combination_rank_;
// functions
/** combination enumerator */
success_t generate_next_symmetric_index(Vector &index) {
int i, ok_so_far;
int n = index.length();
int top = n-1;
do {
index[top] += 1;
ok_so_far = 1;
if (index[top] >= data_.n_cols()) {
index[top] = -1;
top -= 1;
ok_so_far = 0;
if (top < 0) {
return SUCCESS_FAIL;
}
}
for (i = 0; i < top && ok_so_far; i++) {
if (index[top] <= index[i]) {
ok_so_far = 0;
}
}
if(ok_so_far) {
top += 1;
}
}
while (top < n);
return SUCCESS_PASS;
}
/** test whether node a is an ancestor node of node b */
int as_indexes_strictly_surround_bs(Tree *a, Tree *b) {
return (a->begin() < b->begin() && a->end() >= b->end()) ||
(a->begin() <= b->begin() && a->end() > b->end());
}
/**
* Compute the total number of n-tuples by recursively splitting up
* the i-th node
/** @brief The negative force due to the multibody potential on each
* particle. Each column is a force vector on each particle.
*/
double two_ttn(int b, ArrayList<Tree *> &nodes, int i) {
Vector negative_force1_e_;
double result = 0.0;
Tree *kni = nodes[i];
nodes[i] = kni->left();
result += ttn(b, nodes);
nodes[i] = kni->right();
result += ttn(b, nodes);
nodes[i] = kni;
return result;
}
/** @brief The upper bound on the negative force due to the
* multibody potential on each particle. Each column is a
* force vector on each particle.
*/
Vector negative_force1_u_;
/** Compute the total number of n-tuples */
double ttn(int b, ArrayList<Tree *> &nodes) {
Tree *bkn = nodes[b];
double result;
int n = nodes.size();
/** @brief The lower bound on the positive force due to the
* multibody potential on each particle. Each column is a
* force vector on each particle.
*/
Vector positive_force1_l_;
if(b == n - 1) {
result = (double) bkn->count();
}
else {
int j;
int conflict = 0;
int simple_product = 1;
result = (double) bkn->count();
for(j = b + 1 ; j < n && !conflict; j++) {
Tree *knj = nodes[j];
if (bkn->begin() >= knj->end() - 1) {
conflict = 1;
}
else if(nodes[j - 1]->end() - 1 > knj->begin()) {
simple_product = 0;
}
}
if(conflict) {
result = 0.0;
}
else if(simple_product) {
for(j = b + 1; j < n; j++) {
result *= nodes[j]->count();
}
}
else {
int jdiff = -1;
/** @brief The positive force due to the multibody potential on each
* particle. Each column is a force vector on each particle.
*/
Vector positive_force1_e_;
// undefined... will eventually point to the
// lowest j > b such that nodes[j] is different from
// bkn
for(j = b + 1; jdiff < 0 && j < n; j++) {
Tree *knj = nodes[j];
if(bkn->begin() != knj->begin() ||
bkn->end() - 1 != knj->end() - 1) {
jdiff = j;
}
}
/** @brief The negative force due to the multibody potential on each
* particle. Each column is a force vector on each particle.
*/
Matrix negative_force2_e_;
if(jdiff < 0) {
result = math::BinomialCoefficient(bkn->count(), n - b);
}
else {
Tree *dkn = nodes[jdiff];
/** @brief The upper bound on the negative force due to the
* multibody potential on each particle. Each column is a
* force vector on each particle.
*/
Matrix negative_force2_u_;
if(dkn->begin() >= bkn->end() - 1) {
result = math::BinomialCoefficient(bkn->count(), jdiff - b);
if(result > 0.0) {
result *= ttn(jdiff, nodes);
}
}
else if(as_indexes_strictly_surround_bs(bkn, dkn)) {
result = two_ttn(b, nodes, b);
}
else if(as_indexes_strictly_surround_bs(dkn, bkn)) {
result = two_ttn(b, nodes, jdiff);
}
}
}
}
return result;
}
/** @brief The lower bound on the positive force due to the
* multibody potential on each particle. Each column is a
* force vector on each particle.
*/
Matrix positive_force2_l_;
/** Heuristic for node splitting - find the node with most points */
int FindSplitNode(ArrayList<Tree *> &nodes) {
/** @brief The positive force due to the multibody potential on each
* particle. Each column is a force vector on each particle.
*/
Matrix positive_force2_e_;
int global_index = -1;
int global_min = 0;
/** @brief The total estimated force due to the multibody potential
* on each particle. Each column is a force vector on each
* particle.
*/
Matrix total_force_e_;
for(index_t i = 0; i < non_leaf_indices_.size(); i++) {
/////////// Helper Functions //////////
/*
int non_leaf_index = non_leaf_indices_[i];
double minimum_side_length = MAXDOUBLE;
/** @brief Adds the postponed information from a leaf node to the
* point's contribution.
*/
void AddPostponed(Tree *node, index_t destination);
// find out the minimum side length
for(index_t j = 0; j < data_.n_rows(); j++) {
DRange range = nodes[non_leaf_index]->bound().get(j);
double side_length = range.width();
if(side_length < minimum_side_length) {
minimum_side_length = side_length;
}
}
if(minimum_side_length > global_min) {
global_min = minimum_side_length;
global_index = non_leaf_index;
}
*/
int non_leaf_index = non_leaf_indices_[i];
if(nodes[non_leaf_index]->count() > global_min) {
global_min = nodes[non_leaf_index]->count();
global_index = non_leaf_index;
}
}
return global_index;
}
/** @brief Adds the postponed information from a node to another.
*/
void AddPostponed(Tree *source_node, Tree *destination_node);
/** @brief Tests whether node a is an ancestor node of node b.
*/
int as_indexes_strictly_surround_bs(Tree *a, Tree *b);
/** @brief Compute the total number of n-tuples by recursively
* splitting up the i-th node
*/
double two_ttn(int b, ArrayList<Tree *> &nodes, int i);
/** @brief Compute the total number of n-tuples.
*/
double ttn(int b, ArrayList<Tree *> &nodes);
/** @brief Heuristic for node splitting - find the node with most
* points.
*/
int FindSplitNode(ArrayList<Tree *> &nodes);
/** Pruning rule */
int Prunable(ArrayList<Tree *> &nodes, double num_tuples,
double *allowed_err) {
bool Prunable(ArrayList<Tree *> &nodes, double num_tuples,
double *allowed_err);
double pos_min_potential, pos_max_potential;
double neg_min_potential, neg_max_potential;
double pos_lower_change, neg_upper_change;
double pos_error, pos_estimate, neg_error, neg_estimate;
double error;
/** @brief The base exhaustive computations.
*/
void MTMultibodyBase(const ArrayList<Tree *> &nodes, int level);
/** @brief The post-processing function to push down all unclaimed
* approximations.
*/
void PostProcess(Tree *node);
// compute pairwise bounding box distances
for(index_t i = 0; i < mkernel_.order(); i++) {
node_bounds_[i] = &(nodes[i]->bound());
}
mkernel_.EvalNodes(node_bounds_, &neg_min_potential, &neg_max_potential,
&pos_min_potential, &pos_max_potential);
if(isnan(pos_max_potential) || isinf(pos_max_potential)) {
return 0;
}
pos_lower_change = num_tuples * pos_min_potential;
pos_error = 0.5 * num_tuples * (pos_max_potential - pos_min_potential);
pos_estimate = 0.5 * num_tuples * (pos_min_potential + pos_max_potential);
neg_upper_change = num_tuples * neg_max_potential;
neg_error = 0.5 * num_tuples * (neg_max_potential - neg_min_potential);
neg_estimate = 0.5 * num_tuples * (neg_min_potential + neg_max_potential);
// compute whether the error is below the threshold
*allowed_err = tau_ * (pos_potential_l_ + pos_lower_change -
(neg_potential_u_ + neg_upper_change)) *
((num_tuples + extra_token_) / total_num_tuples_);
error = max(pos_error, neg_error);
if(likely(error >= 0) && error <= (*allowed_err)) {
pos_potential_l_ += pos_lower_change;
pos_potential_e_ += pos_estimate;
neg_potential_u_ += neg_upper_change;
neg_potential_e_ += neg_estimate;
extra_token_ = num_tuples + extra_token_ - error * total_num_tuples_ /
(tau_ * (pos_potential_l_ -neg_potential_u_));
DEBUG_ASSERT(extra_token_ >= 0);
return 1;
}
return 0;
}
/** Pruning rule for series approxiamation approach */
int PrunableSeriesExpansion(ArrayList<Tree *> &nodes, double num_tuples,
double allowed_err) {
if(nodes[0] != nodes[1] && nodes[0] != nodes[2] && nodes[1] != nodes[2]) {
Matrix distmat;
double actual_error1 = 0;
double actual_error2 = 0;
double actual_error3 = 0;
distmat.Alias(mkernel_.EvalMinMaxDsqds(node_bounds_));
double max_ij = mkernel_.EvalUnnormOnSqOnePair(distmat.get(0, 1));
double max_ik = mkernel_.EvalUnnormOnSqOnePair(distmat.get(0, 2));
double max_jk = mkernel_.EvalUnnormOnSqOnePair(distmat.get(1, 2));
double min_ij = mkernel_.EvalUnnormOnSqOnePair(distmat.get(1, 0));
double min_ik = mkernel_.EvalUnnormOnSqOnePair(distmat.get(2, 0));
double min_jk = mkernel_.EvalUnnormOnSqOnePair(distmat.get(2, 1));
FarFieldExpansion<TKernelAux> &coeffs0 =
nodes[0]->stat().get_farfield_coeffs();
FarFieldExpansion<TKernelAux> &coeffs1 =
nodes[1]->stat().get_farfield_coeffs();
FarFieldExpansion<TKernelAux> &coeffs2 =
nodes[2]->stat().get_farfield_coeffs();
double total_relerr = allowed_err /
(num_tuples * max_ij * max_ik * max_jk);
double rel_err = max(pow(total_relerr + 1, 1.0 / 3.0) - 1, 0.0);
// compute the required number of terms
int order_ij = coeffs0.OrderForConvertingToLocal(nodes[0]->bound(),
nodes[1]->bound(),
distmat.get(0, 1),
distmat.get(1, 0),
min_ij * rel_err,
&actual_error1);
int order_ik = coeffs1.OrderForConvertingToLocal(nodes[0]->bound(),
nodes[2]->bound(),
distmat.get(0, 2),
distmat.get(2, 0),
min_ik * rel_err,
&actual_error2);
int order_jk = coeffs2.OrderForConvertingToLocal(nodes[1]->bound(),
nodes[2]->bound(),
distmat.get(1, 2),
distmat.get(2, 1),
min_jk * rel_err,
&actual_error3);
int max_order = coeffs0.get_max_order() / 2 - 1;
if(order_ij >= 0 && order_ik >= 0 && order_jk >= 0 &&
order_ij < max_order && order_ik < max_order &&
order_jk < max_order &&
ka_.sea_.get_total_num_coeffs(order_ij) *
ka_.sea_.get_total_num_coeffs(order_ik) *
ka_.sea_.get_total_num_coeffs(order_jk) <
nodes[0]->count() * nodes[1]->count() * nodes[2]->count()) {
coeffs0.RefineCoeffs(data_, weights_, nodes[0]->begin(),
nodes[0]->end(), order_ij);
coeffs1.RefineCoeffs(data_, weights_, nodes[1]->begin(),
nodes[1]->end(), order_ik);
coeffs2.RefineCoeffs(data_, weights_, nodes[2]->begin(),
nodes[2]->end(), order_jk);
pos_potential_l_ += num_tuples * min_ij * min_ik * min_jk;
pos_potential_e_ += coeffs0.ConvolveField(coeffs1, coeffs2, order_ij,
order_ik, order_jk);
// the maximum relative error incurred
double max_rel_err_incurred1 = actual_error1 / min_ij;
double max_rel_err_incurred2 = actual_error2 / min_ik;
double max_rel_err_incurred3 = actual_error3 / min_jk;
double max_rel_err_incurred =
max_rel_err_incurred1 + max_rel_err_incurred2 +
max_rel_err_incurred3 + max_rel_err_incurred1 *
max_rel_err_incurred2 + max_rel_err_incurred1 *
max_rel_err_incurred3 + max_rel_err_incurred2 *
max_rel_err_incurred3 + max_rel_err_incurred1 *
max_rel_err_incurred2 * max_rel_err_incurred3;
double error = max_rel_err_incurred * max_ij * max_ik * max_jk;
extra_token_ = num_tuples + extra_token_ - error * total_num_tuples_ /
(tau_ * pos_potential_l_);
DEBUG_ASSERT(extra_token_ >= 0);
return 1;
}
return 0;
}
return 0;
}
/** Pruning rule for series approxiamation second approach */
int PrunableSeriesExpansion2(ArrayList<Tree *> &nodes, double num_tuples,
double allowed_err) {
if(nodes[0] != nodes[1] && nodes[0] != nodes[2] && nodes[1] != nodes[2]) {
Matrix distmat;
double actual_error2 = 0;
double actual_error3 = 0;
distmat.Alias(mkernel_.EvalMinMaxDsqds(node_bounds_));
double max_ik = mkernel_.EvalUnnormOnSqOnePair(distmat.get(0, 2));
double max_jk = mkernel_.EvalUnnormOnSqOnePair(distmat.get(1, 2));
double min_ij = mkernel_.EvalUnnormOnSqOnePair(distmat.get(1, 0));
double min_ik = mkernel_.EvalUnnormOnSqOnePair(distmat.get(2, 0));
double min_jk = mkernel_.EvalUnnormOnSqOnePair(distmat.get(2, 1));
FarFieldExpansion<TKernelAux> &coeffs0 =
nodes[0]->stat().get_farfield_coeffs();
FarFieldExpansion<TKernelAux> &coeffs1 =
nodes[1]->stat().get_farfield_coeffs();
FarFieldExpansion<TKernelAux> &coeffs2 =
nodes[2]->stat().get_farfield_coeffs();
double total_relerr = allowed_err /
(num_tuples * max_ik * max_jk);
double rel_err = max(pow(total_relerr + 1, 1.0 / 2.0) - 1, 0.0);
// compute the required number of terms
int order_ik = -1;
int order_jk = -1;
if(min_ik * rel_err > 0 && min_jk * rel_err > 0) {
order_ik = coeffs1.OrderForConvertingToLocal(nodes[0]->bound(),
nodes[2]->bound(),
distmat.get(0, 2),
distmat.get(2, 0),
min_ik * rel_err,
&actual_error2);
order_jk = coeffs2.OrderForConvertingToLocal(nodes[1]->bound(),
nodes[2]->bound(),
distmat.get(1, 2),
distmat.get(2, 1),
min_jk * rel_err,
&actual_error3);
}
int max_order = coeffs1.get_max_order() / 2 - 1;
if(order_ik >= 0 && order_jk >= 0 && order_ik < max_order &&
order_jk < max_order &&
ka_.sea_.get_total_num_coeffs(order_ik) *
ka_.sea_.get_total_num_coeffs(order_jk) < 2 * nodes[2]->count()) {
coeffs2.RefineCoeffs(data_, weights_, nodes[2]->begin(),
nodes[2]->end(), order_jk);
pos_potential_l_ += num_tuples * min_ij * min_ik * min_jk;
pos_potential_e_ += coeffs0.MixField(data_, nodes[0]->begin(),
nodes[0]->end(),
nodes[1]->begin(),
nodes[1]->end(), coeffs1,
coeffs2, order_ik, order_jk);
// the maximum relative error incurred
double max_rel_err_incurred2 = actual_error2 / min_ik;
double max_rel_err_incurred3 = actual_error3 / min_jk;
double max_rel_err_incurred =
max_rel_err_incurred2 + max_rel_err_incurred3 +
max_rel_err_incurred2 * max_rel_err_incurred3;
double error = max_rel_err_incurred * max_ik * max_jk;
extra_token_ = num_tuples + extra_token_ - error * total_num_tuples_ /
(tau_ * pos_potential_l_);
DEBUG_ASSERT(extra_token_ >= 0);
return 1;
}
return 0;
}
return 0;
}
/** Base exhaustive case */
void MTMultibodyBase(ArrayList<Tree *> &nodes, int level) {
int start_index;
int num_nodes = nodes.size();
if(level < num_nodes) {
/* run over each point in this node */
if(level > 0) {
if(nodes[level - 1] == nodes[level]) {
start_index = exhaustive_indices_[level - 1] + 1;
}
else {
start_index = nodes[level]->begin();
}
}
else {
start_index = nodes[level]->begin();
}
for(index_t i = start_index; i < (nodes[level])->end(); i++) {
exhaustive_indices_[level] = i;
MTMultibodyBase(nodes, level + 1);
}
}
else {
double neg, pos;
// complete the table of distance computation
mkernel_.Eval(data_, exhaustive_indices_, &neg, &pos);
neg_potential_e_ += neg;
neg_potential_u_ += neg;
pos_potential_e_ += pos;
pos_potential_l_ += pos;
}
}
/** Main multitree recursion */
void MTMultibody(ArrayList<Tree *> &nodes, double num_tuples) {
double allowed_err = 0;
NumNodesExpanded_++;
if(Prunable(nodes, num_tuples, &allowed_err)) {
return;
}
else if(PrunableSeriesExpansion2(nodes, num_tuples, allowed_err)) {
NumPrunes_++;
return;
}
// figure out which ones are non-leaves
non_leaf_indices_.Resize(0);
for(index_t i = 0; i < 3; i++) {
if(!(nodes[i]->is_leaf())) {
non_leaf_indices_.PushBackCopy(i);
}
}
// all leaves, then base case
if(non_leaf_indices_.size() == 0) {
MTMultibodyBase(nodes, 0);
extra_token_ += num_tuples;
return;
}
// else, split an internal node and recurse
else {
int split_index;
double new_num_tuples;
// copy to new nodes
ArrayList<Tree *> new_nodes;
new_nodes.Init(3);
for(index_t i = 0; i < 3; i++) {
new_nodes[i] = nodes[i];
}
// apply splitting heuristic
split_index = FindSplitNode(nodes);
// recurse to the left
new_nodes[split_index] = nodes[split_index]->left();
new_num_tuples = ttn(0, new_nodes);
if(new_num_tuples > 0) {
MTMultibody(new_nodes, new_num_tuples);
}
// recurse to the right
new_nodes[split_index] = nodes[split_index]->right();
new_num_tuples = ttn(0, new_nodes);
if(new_num_tuples > 0) {
MTMultibody(new_nodes, new_num_tuples);
}
}
}
/** @brief The main multitree recursion.
*/
void MTMultibody(ArrayList<Tree *> &nodes, double num_tuples);
};
#include "multibody_impl.h"
#undef INSIDE_MULTIBODY_H
#endif
fastlib2/contrib/dongryel/multibody/multibody_impl.h
+368
View File
@@ -0,0 +1,368 @@
#ifndef INSIDE_MULTIBODY_H
#error "This is not a public header file!"
#endif
#ifndef MULTIBODY_IMPL_H
#define MULTIBODY_IMPL_H
template<typename TMultibodyKernel>
int MultitreeMultibody<TMultibodyKernel>::as_indexes_strictly_surround_bs
(Tree *a, Tree *b) {
return (a->begin() < b->begin() && a->end() >= b->end()) ||
(a->begin() <= b->begin() && a->end() > b->end());
}
template<typename TMultibodyKernel>
double MultitreeMultibody<TMultibodyKernel>::ttn(int b,
ArrayList<Tree *> &nodes) {
Tree *bkn = nodes[b];
double result;
int n = nodes.size();
if(b == n - 1) {
result = (double) bkn->count();
}
else {
int j;
int conflict = 0;
int simple_product = 1;
result = (double) bkn->count();
for(j = b + 1 ; j < n && !conflict; j++) {
Tree *knj = nodes[j];
if (bkn->begin() >= knj->end() - 1) {
conflict = 1;
}
else if(nodes[j - 1]->end() - 1 > knj->begin()) {
simple_product = 0;
}
}
if(conflict) {
result = 0.0;
}
else if(simple_product) {
for(j = b + 1; j < n; j++) {
result *= nodes[j]->count();
}
}
else {
int jdiff = -1;
// undefined... will eventually point to the
// lowest j > b such that nodes[j] is different from
// bkn
for(j = b + 1; jdiff < 0 && j < n; j++) {
Tree *knj = nodes[j];
if(bkn->begin() != knj->begin() ||
bkn->end() - 1 != knj->end() - 1) {
jdiff = j;
}
}
if(jdiff < 0) {
result = math::BinomialCoefficient(bkn->count(), n - b);
}
else {
Tree *dkn = nodes[jdiff];
if(dkn->begin() >= bkn->end() - 1) {
result = math::BinomialCoefficient(bkn->count(), jdiff - b);
if(result > 0.0) {
result *= ttn(jdiff, nodes);
}
}
else if(as_indexes_strictly_surround_bs(bkn, dkn)) {
result = two_ttn(b, nodes, b);
}
else if(as_indexes_strictly_surround_bs(dkn, bkn)) {
result = two_ttn(b, nodes, jdiff);
}
}
}
}
return result;
}
template<typename TMultibodyKernel>
double MultitreeMultibody<TMultibodyKernel>::two_ttn
(int b, ArrayList<Tree *> &nodes, int i) {
double result = 0.0;
Tree *kni = nodes[i];
nodes[i] = kni->left();
result += ttn(b, nodes);
nodes[i] = kni->right();
result += ttn(b, nodes);
nodes[i] = kni;
return result;
}
template<typename TMultibodyKernel>
int MultitreeMultibody<TMultibodyKernel>::FindSplitNode
(ArrayList<Tree *> &nodes) {
int global_index = -1;
int global_min = 0;
for(index_t i = 0; i < non_leaf_indices_.size(); i++) {
/*
int non_leaf_index = non_leaf_indices_[i];
double minimum_side_length = MAXDOUBLE;
// find out the minimum side length
for(index_t j = 0; j < data_.n_rows(); j++) {
DRange range = nodes[non_leaf_index]->bound().get(j);
double side_length = range.width();
if(side_length < minimum_side_length) {
minimum_side_length = side_length;
}
}
if(minimum_side_length > global_min) {
global_min = minimum_side_length;
global_index = non_leaf_index;
}
*/
int non_leaf_index = non_leaf_indices_[i];
if(nodes[non_leaf_index]->count() > global_min) {
global_min = nodes[non_leaf_index]->count();
global_index = non_leaf_index;
}
}
return global_index;
}
template<typename TMultibodyKernel>
bool MultitreeMultibody<TMultibodyKernel>::Prunable
(ArrayList<Tree *> &nodes, double num_tuples, double *allowed_err) {
return 0;
}
template<typename TMultibodyKernel>
void MultitreeMultibody<TMultibodyKernel>::AddPostponed
(Tree *source_node, Tree *destination_node) {
destination_node->stat().postponed_negative_gradient1_e +=
source_node->stat().postponed_negative_gradient1_e;
destination_node->stat().postponed_negative_gradient1_u +=
source_node->stat().postponed_negative_gradient1_u;
destination_node->stat().postponed_positive_gradient1_l +=
source_node->stat().postponed_positive_gradient1_l;
destination_node->stat().postponed_positive_gradient1_e +=
source_node->stat().postponed_positive_gradient1_e;
la::AddTo(data_.n_rows(),
source_node->stat().postponed_negative_gradient2_e.ptr(),
destination_node->stat().postponed_negative_gradient2_e.ptr());
la::AddTo(data_.n_rows(),
source_node->stat().postponed_negative_gradient2_u.ptr(),
destination_node->stat().postponed_negative_gradient2_u.ptr());
la::AddTo(data_.n_rows(),
source_node->stat().postponed_positive_gradient2_l.ptr(),
destination_node->stat().postponed_positive_gradient2_l.ptr());
la::AddTo(data_.n_rows(),
source_node->stat().postponed_positive_gradient2_e.ptr(),
destination_node->stat().postponed_positive_gradient2_e.ptr());
}
template<typename TMultibodyKernel>
void MultitreeMultibody<TMultibodyKernel>::AddPostponed(Tree *node,
index_t destination) {
negative_force1_e_[destination] +=
node->stat().postponed_negative_gradient1_e;
negative_force1_u_[destination] +=
node->stat().postponed_negative_gradient1_u;
positive_force1_l_[destination] +=
node->stat().postponed_positive_gradient1_l;
positive_force1_e_[destination] +=
node->stat().postponed_positive_gradient1_e;
la::AddTo(data_.n_rows(), node->stat().postponed_negative_gradient2_e.ptr(),
negative_force2_e_.GetColumnPtr(destination));
la::AddTo(data_.n_rows(), node->stat().postponed_negative_gradient2_u.ptr(),
negative_force2_u_.GetColumnPtr(destination));
la::AddTo(data_.n_rows(), node->stat().postponed_positive_gradient2_l.ptr(),
positive_force2_l_.GetColumnPtr(destination));
la::AddTo(data_.n_rows(), node->stat().postponed_positive_gradient2_e.ptr(),
positive_force2_e_.GetColumnPtr(destination));
}
template<typename TMultibodyKernel>
void MultitreeMultibody<TMultibodyKernel>::MTMultibodyBase
(const ArrayList<Tree *> &nodes, int level) {
int start_index;
int num_nodes = nodes.size();
// Recurse to get a $n$ tuple.
if(level < num_nodes) {
// Run over each point in this node.
if(level > 0) {
if(nodes[level - 1] == nodes[level]) {
start_index = exhaustive_indices_[level - 1] + 1;
}
else {
start_index = nodes[level]->begin();
}
}
else {
start_index = nodes[level]->begin();
}
for(index_t i = start_index; i < (nodes[level])->end(); i++) {
exhaustive_indices_[level] = i;
MTMultibodyBase(nodes, level + 1);
}
}
else {
// Incorporate postponed force contribution for the given triple
// of atoms.
for(index_t i = 0; i < nodes.size(); i++) {
AddPostponed(nodes[i], exhaustive_indices_[i]);
}
// Complete the contribution among three atoms.
mkernel_.Eval(data_, exhaustive_indices_,
negative_force1_e_, negative_force1_u_,
positive_force1_l_, positive_force1_e_,
negative_force2_e_, negative_force2_u_,
positive_force2_l_, positive_force2_e_);
}
// Clear all postponed force contribution after incorporating.
if(level == 0) {
for(index_t i = 0; i < nodes.size(); i++) {
nodes[i]->stat().SetZero();
}
}
}
template<typename TMultibodyKernel>
void MultitreeMultibody<TMultibodyKernel>::PostProcess(Tree *node) {
//
// For a leaf node,
if(node->is_leaf()) {
for(index_t q = node->begin(); q < node->end(); q++) {
// Add postponed contribution to each point's force vector.
AddPostponed(node, q);
// Now, reconstruct the force vector from the complete
// approximations.
double *query_total_force_e = total_force_e_.GetColumnPtr(q);
for(index_t d = 0; d < data_.n_rows(); d++) {
// First, add in the negative contributions then the positive
// contributions.
if(data_.get(d, q) < 0) {
query_total_force_e[d] += (-data_.get(d, q) * negative_force1_e_[q] +
negative_force2_e_.get(d, q));
query_total_force_e[d] += (-data_.get(d, q) * positive_force1_e_[q] +
positive_force2_e_.get(d, q));
}
else {
query_total_force_e[d] += (-data_.get(d, q) * positive_force1_e_[q] +
negative_force2_e_.get(d, q));
query_total_force_e[d] += (-data_.get(d, q) * negative_force1_e_[q] +
positive_force2_e_.get(d, q));
}
} // end of iterating over each dimension...
} // end of iterating over each query point...
// Clear postponed information.
node->stat().SetZero();
}
else {
// Push down postponed contributions to the nodes below and clear
// them.
AddPostponed(node, node->left());
AddPostponed(node, node->right());
node->stat().SetZero();
// Recurse.
PostProcess(node->left());
PostProcess(node->right());
}
}
template<typename TMultibodyKernel>
void MultitreeMultibody<TMultibodyKernel>::MTMultibody
(ArrayList<Tree *> &nodes, double num_tuples) {
double allowed_err = 0;
/*
if(Prunable(nodes, num_tuples, &allowed_err)) {
return;
}
*/
// Figure out which ones are non-leaves.
non_leaf_indices_.Resize(0);
for(index_t i = 0; i < 3; i++) {
if(!(nodes[i]->is_leaf())) {
non_leaf_indices_.PushBackCopy(i);
}
}
// All leaves, then base case.
if(non_leaf_indices_.size() == 0) {
MTMultibodyBase(nodes, 0);
return;
}
// Else, split an internal node and recurse.
else {
int split_index;
double new_num_tuples;
// Copy to new nodes list before recursing.
ArrayList<Tree *> new_nodes;
new_nodes.Init(mkernel_.order());
for(index_t i = 0; i < mkernel_.order(); i++) {
new_nodes[i] = nodes[i];
}
// Apply splitting heuristic.
split_index = FindSplitNode(nodes);
// Push down approximations downward for the node that is to be
// expanded.
AddPostponed(nodes[split_index], nodes[split_index]->left());
AddPostponed(nodes[split_index], nodes[split_index]->right());
nodes[split_index]->stat().SetZero();
// Recurse to the left.
new_nodes[split_index] = nodes[split_index]->left();
new_num_tuples = ttn(0, new_nodes);
// If the current node combination is valid, then recurse.
if(new_num_tuples > 0) {
MTMultibody(new_nodes, new_num_tuples);
}
// Recurse to the right.
new_nodes[split_index] = nodes[split_index]->right();
new_num_tuples = ttn(0, new_nodes);
// If the current node combination is valid, then recurse.
if(new_num_tuples > 0) {
MTMultibody(new_nodes, new_num_tuples);
}
}
}
#endif
fastlib2/contrib/dongryel/multibody/multibody_kernel.h
+302 -233
View File
@@ -4,275 +4,344 @@
#include "fastlib/fastlib.h"
#include "mlpack/series_expansion/kernel_aux.h"
class GaussianThreeBodyKernel {
class AxilrodTellerForceKernel {
private:
GaussianKernel kernel_;
Matrix distmat_;
public:
GaussianThreeBodyKernel() {}
~GaussianThreeBodyKernel() {}
// getters and setters
double bandwidth_sq() const { return kernel_.bandwidth_sq(); }
const Matrix &pairwise_dsqd() const { return distmat_; }
void Init(double bandwidth_in) {
kernel_.Init(bandwidth_in);
distmat_.Init(3, 3);
}
int order() {
return 3;
}
double EvalUnnormOnSqOnePair(double sqdist) const {
return kernel_.EvalUnnormOnSq(sqdist);
}
void EvalUnnormOnSq(const Matrix &sqdists, double *neg, double *pos) const {
*pos = 1;
*neg = 0;
for(index_t i = 0; i < sqdists.n_cols(); i++) {
for(index_t j = i + 1; j < sqdists.n_cols(); j++) {
(*pos) *= kernel_.EvalUnnormOnSq(sqdists.get(i, j));
}
}
}
void EvalMinMax(double *negmin, double *negmax,
double *posmin, double *posmax) const {
*negmin = *negmax = 0;
*posmin = 1.0;
*posmax = 1.0;
for(index_t i = 0; i < 3; i++) {
for(index_t j = i + 1; j < 3; j++) {
*posmin = (*posmin) * kernel_.EvalUnnormOnSq(distmat_.get(j, i));
*posmax = (*posmax) * kernel_.EvalUnnormOnSq(distmat_.get(i, j));
}
}
}
const Matrix &EvalMinMaxDsqds
(const ArrayList<DHrectBound<2> *> &node_bounds) {
int num_nodes = node_bounds.size();
for(index_t i = 0; i < num_nodes - 1; i++) {
DHrectBound<2> *node_i_bound = node_bounds[i];
for(index_t j = i + 1; j < num_nodes; j++) {
DHrectBound<2> *node_j_bound = node_bounds[j];
double dmin = node_i_bound->MinDistanceSq(*node_j_bound);
double dmax = node_i_bound->MaxDistanceSq(*node_j_bound);
distmat_.set(i, j, dmin);
distmat_.set(j, i, dmax);
}
}
return distmat_;
}
void Eval(const Matrix &data, const ArrayList<int> &indices, double *neg,
double *pos) {
for(index_t i = 0; i < indices.size(); i++) {
const double *i_col = data.GetColumnPtr(indices[i]);
for(index_t j = i + 1; j < indices.size(); j++) {
const double *j_col = data.GetColumnPtr(indices[j]);
distmat_.set(i, j, la::DistanceSqEuclidean(data.n_rows(), i_col,
j_col));
}
}
EvalUnnormOnSq(distmat_, neg, pos);
}
void EvalNodes(const ArrayList<DHrectBound<2> *> &node_bounds,
double *negmin, double *negmax, double *posmin,
double *posmax) {
int num_nodes = node_bounds.size();
for(index_t i = 0; i < num_nodes - 1; i++) {
DHrectBound<2> *node_i_bound = node_bounds[i];
for(index_t j = i + 1; j < num_nodes; j++) {
DHrectBound<2> *node_j_bound = node_bounds[j];
double dmin = node_i_bound->MinDistanceSq(*node_j_bound);
double dmax = node_i_bound->MaxDistanceSq(*node_j_bound);
distmat_.set(i, j, dmin);
distmat_.set(j, i, dmax);
}
}
EvalMinMax(negmin, negmax, posmin, posmax);
}
};
class AxilrodTellerKernel {
private:
Matrix distmat_;
////////// Private Member Constants //////////
/** @brief The "nu" constant in front of the potential.
*/
static const double AXILROD_TELLER_COEFF = 1e-18;
public:
////////// Private Member Variables //////////
AxilrodTellerKernel() {}
/** @brief The temporary matrix to store pairwise distances.
*/
Matrix distmat_;
~AxilrodTellerKernel() {}
// getters and setters
double bandwidth_sq() const { return 1; }
const Matrix &pairwise_dsqd() const { return distmat_; }
/** @brief The temporary ArrayList to store the mapped indices for
* computing the gradient.
*/
ArrayList<index_t> index_orders_;
void Init(double bandwidth_in) {
distmat_.Init(3, 3);
////////// Private Member Functions //////////
void force_(const Matrix &data, const ArrayList<index_t> &indices,
double &negative_gradient1, double &positive_gradient1,
double &negative_gradient2, double &positive_gradient2,
double &negative_gradient3, double &positive_gradient3,
Vector &negative_force1_e, Vector &negative_force1_u,
Vector &positive_force1_l, Vector &positive_force1_e,
Matrix &negative_force2_e, Matrix &negative_force2_u,
Matrix &positive_force2_l, Matrix &positive_force2_e) {
// Negative contribution to the first component.
negative_force1_e[indices[index_orders_[0]]] +=
negative_gradient1 + negative_gradient2;
negative_force1_u[indices[index_orders_[0]]] +=
negative_gradient1 + negative_gradient2;
// Positive contribution to the first component.
positive_force1_l[indices[index_orders_[0]]] +=
positive_gradient1 + positive_gradient2;
positive_force1_e[indices[index_orders_[0]]] +=
positive_gradient1 + positive_gradient2;
// Negative contribution to the second component.
la::AddExpert(data.n_rows(), negative_gradient1,
data.GetColumnPtr(indices[index_orders_[1]]),
negative_force2_e.GetColumnPtr(indices[index_orders_[0]]));
la::AddExpert(data.n_rows(), negative_gradient2,
data.GetColumnPtr(indices[index_orders_[2]]),
negative_force2_e.GetColumnPtr(indices[index_orders_[0]]));
la::AddExpert(data.n_rows(), negative_gradient1,
data.GetColumnPtr(indices[index_orders_[1]]),
negative_force2_u.GetColumnPtr(indices[index_orders_[0]]));
la::AddExpert(data.n_rows(), negative_gradient2,
data.GetColumnPtr(indices[index_orders_[2]]),
negative_force2_u.GetColumnPtr(indices[index_orders_[0]]));
// Positive contribution to the second component.
la::AddExpert(data.n_rows(), positive_gradient1,
data.GetColumnPtr(indices[index_orders_[1]]),
positive_force2_e.GetColumnPtr(indices[index_orders_[0]]));
la::AddExpert(data.n_rows(), positive_gradient2,
data.GetColumnPtr(indices[index_orders_[2]]),
positive_force2_e.GetColumnPtr(indices[index_orders_[0]]));
la::AddExpert(data.n_rows(), positive_gradient1,
data.GetColumnPtr(indices[index_orders_[1]]),
positive_force2_l.GetColumnPtr(indices[index_orders_[0]]));
la::AddExpert(data.n_rows(), positive_gradient2,
data.GetColumnPtr(indices[index_orders_[2]]),
positive_force2_l.GetColumnPtr(indices[index_orders_[0]]));
}
void gradient_(const ArrayList<index_t> &index_orders,
double &minimum_negative_gradient,
double *maximum_negative_gradient,
double &minimum_positive_gradient,
double *maximum_positive_gradient) {
double min_dsqd1 = distmat_.get(index_orders[0], index_orders[1]);
double min_dist1 = sqrt(min_dsqd1);
double min_dqrt1 = math::Sqr(min_dsqd1);
double min_dsix1 = min_dsqd1 * min_dqrt1;
double max_dsqd1 = distmat_.get(index_orders[1], index_orders[0]);
double max_dist1 = sqrt(max_dsqd1);
double max_dqrt1 = math::Sqr(max_dsqd1);
double max_dsix1 = max_dsqd1 * max_dqrt1;
double min_dsqd2 = distmat_.get(index_orders[0], index_orders[2]);
double min_dist2 = sqrt(min_dsqd2);
double min_dcub2 = min_dsqd2 * min_dist2;
double min_dqui2 = min_dsqd2 * min_dcub2;
double max_dsqd2 = distmat_.get(index_orders[2], index_orders[0]);
double max_dist2 = sqrt(max_dsqd2);
double max_dcub2 = max_dsqd2 * max_dist2;
double max_dqui2 = max_dsqd2 * max_dcub2;
double min_dsqd3 = distmat_.get(index_orders[1], index_orders[2]);
double min_dist3 = sqrt(min_dsqd3);
double min_dcub3 = min_dsqd3 * min_dist3;
double min_dqui3 = min_dsqd3 * min_dcub3;
double max_dsqd3 = distmat_.get(index_orders[2], index_orders[1]);
double max_dist3 = sqrt(max_dsqd3);
double max_dcub3 = max_dsqd3 * max_dist3;
double max_dqui3 = max_dsqd3 * max_dcub3;
double min_common_factor = 3.0 * AXILROD_TELLER_COEFF / (8.0 * max_dist1);
double max_common_factor = 3.0 * AXILROD_TELLER_COEFF / (8.0 * min_dist1);
minimum_negative_gradient = max_common_factor *
(-8.0 / (min_dqrt1 * min_dcub2 * min_dcub3)
- 1.0 / (min_dqui2 * min_dqui3)
- 1.0 / (min_dsqd1 * min_dcub2 * min_dqui3)
- 1.0 / (min_dsqd1 * min_dqui2 * min_dcub3)
- 3.0 / (min_dqrt1 * min_dist2 * min_dqui3)
- 3.0 / (min_dqrt1 * min_dqui2 * min_dist3)
- 5.0 / (min_dsix1 * min_dist2 * min_dcub3)
- 5.0 / (min_dsix1 * min_dcub2 * min_dist3));
if(maximum_negative_gradient) {
*maximum_negative_gradient = min_common_factor *
(-8.0 / (max_dqrt1 * max_dcub2 * max_dcub3)
- 1.0 / (max_dqui2 * max_dqui3)
- 1.0 / (max_dsqd1 * max_dcub2 * max_dqui3)
- 1.0 / (max_dsqd1 * max_dqui2 * max_dcub3)
- 3.0 / (max_dqrt1 * max_dist2 * max_dqui3)
- 3.0 / (max_dqrt1 * max_dqui2 * max_dist3)
- 5.0 / (max_dsix1 * max_dist2 * max_dcub3)
- 5.0 / (max_dsix1 * max_dcub2 * max_dist3));
}
minimum_positive_gradient = min_common_factor *
(5 * min_dist2 / (max_dsix1 * max_dqui3) +
5 * min_dist3 / (max_dsix1 * max_dqui2) +
6 / (max_dqrt1 * max_dcub2 * max_dcub3));
if(maximum_positive_gradient) {
*maximum_positive_gradient = max_common_factor *
(5 * max_dist2 / (min_dsix1 * min_dqui3) +
5 * max_dist3 / (min_dsix1 * min_dqui2) +
6 / (min_dqrt1 * min_dcub2 * min_dcub3));
}
}
public:
////////// Constructor/Destructor //////////
/** @brief The default constructor.
*/
AxilrodTellerForceKernel() {
}
/** @brief The default destructor.
*/
~AxilrodTellerForceKernel() {
}
////////// Getters/Setters //////////
/** @brief Gets the squared distance matrix.
*/
const Matrix &pairwise_squared_distances() const { return distmat_; }
/** @brief Gets the interaction order of the kernel.
*/
int order() {
return 3;
}
double EvalUnnormOnSqOnePair(double sqdist) const {
////////// User-level Functions //////////
// this is a place holder, needs to be corrected...
return 0;
/** @brief Initializes the kernel.
*/
void Init(double bandwidth_in) {
distmat_.Init(3, 3);
index_orders_.Init(3);
}
void EvalUnnormOnSq(const Matrix &sqdists, double *neg, double *pos) const {
/** @brief Computes the pairwise distance among FastLib tree nodes.
*/
template<typename TTree, typename TBound>
void EvalMinMaxSquaredDistances(const ArrayList<TTree *> &tree_nodes) {
*neg = -0.375 *
(sqdists.get(0, 1) * sqdists.get(0, 1) * sqdists.get(0, 1) +
sqdists.get(0, 2) * sqdists.get(0, 2) * sqdists.get(0, 2) +
sqdists.get(1, 2) * sqdists.get(1, 2) * sqdists.get(1, 2)) /
pow(sqdists.get(0, 1) * sqdists.get(0, 2) * distmat_.get(1, 2), 2.5);
*pos = (3 * sqdists.get(0, 1) * sqdists.get(0, 1) *
(sqdists.get(0, 2) + sqdists.get(1, 2)) +
3 * sqdists.get(0, 2) * sqdists.get(1, 2) *
(sqdists.get(0, 2) + sqdists.get(1, 2)) +
sqdists.get(0, 1) *
(3 * sqdists.get(0, 2) * sqdists.get(0, 2) +
2 * sqdists.get(0, 2) * sqdists.get(1, 2) +
3 * sqdists.get(1, 2) * sqdists.get(1, 2))) /
(8 * pow(sqdists.get(0, 1) * sqdists.get(0, 2) *
sqdists.get(1, 2), 2.5));
*neg = AXILROD_TELLER_COEFF * (*neg);
*pos = AXILROD_TELLER_COEFF * (*pos);
}
void EvalMinMax(double *negmin, double *negmax,
double *posmin, double *posmax) const {
*negmin = -0.375 *
(distmat_.get(1, 0) * distmat_.get(1, 0) * distmat_.get(1, 0) +
distmat_.get(2, 0) * distmat_.get(2, 0) * distmat_.get(2, 0) +
distmat_.get(2, 1) * distmat_.get(2, 1) * distmat_.get(2, 1)) /
pow(distmat_.get(0, 1) * distmat_.get(0, 2) * distmat_.get(1, 2), 2.5);
*posmin = (3 * distmat_.get(0, 1) * distmat_.get(0, 1) *
(distmat_.get(0, 2) + distmat_.get(1, 2)) +
3 * distmat_.get(0, 2) * distmat_.get(1, 2) *
(distmat_.get(0, 2) + distmat_.get(1, 2)) +
distmat_.get(0, 1) *
(3 * distmat_.get(0, 2) * distmat_.get(0, 2) +
2 * distmat_.get(0, 2) * distmat_.get(1, 2) +
3 * distmat_.get(1, 2) * distmat_.get(1, 2))) /
(8 * pow(distmat_.get(1, 0) * distmat_.get(2, 0) *
distmat_.get(2, 1), 2.5));
*negmax = -0.375 *
(distmat_.get(0, 1) * distmat_.get(0, 1) * distmat_.get(0, 1) +
distmat_.get(0, 2) * distmat_.get(0, 2) * distmat_.get(0, 2) +
distmat_.get(1, 2) * distmat_.get(1, 2) * distmat_.get(1, 2)) /
pow(distmat_.get(1, 0) * distmat_.get(2, 0) * distmat_.get(2, 1), 2.5);
*posmax = (3 * distmat_.get(1, 0) * distmat_.get(1, 0) *
(distmat_.get(2, 0) + distmat_.get(2, 1)) +
3 * distmat_.get(2, 0) * distmat_.get(2, 1) *
(distmat_.get(2, 0) + distmat_.get(2, 1)) +
distmat_.get(1, 0) *
(3 * distmat_.get(2, 0) * distmat_.get(2, 0) +
2 * distmat_.get(2, 0) * distmat_.get(2, 1) +
3 * distmat_.get(2, 1) * distmat_.get(2, 1))) /
(8 * pow(distmat_.get(0, 1) * distmat_.get(0, 2) *
distmat_.get(1, 2), 2.5));
*negmin = AXILROD_TELLER_COEFF * (*negmin);
*negmax = AXILROD_TELLER_COEFF * (*negmax);
*posmin = AXILROD_TELLER_COEFF * (*posmin);
*posmax = AXILROD_TELLER_COEFF * (*posmax);
}
const Matrix &EvalMinMaxDsqds
(const ArrayList<DHrectBound<2> *> &node_bounds) {
int num_nodes = node_bounds.size();
int num_nodes = tree_nodes.size();
for(index_t i = 0; i < num_nodes - 1; i++) {
DHrectBound<2> *node_i_bound = node_bounds[i];
const TBound &node_i_bound = tree_nodes[i]->bound();
for(index_t j = i + 1; j < num_nodes; j++) {
DHrectBound<2> *node_j_bound = node_bounds[j];
double dmin = node_i_bound->MinDistanceSq(*node_j_bound);
double dmax = node_i_bound->MaxDistanceSq(*node_j_bound);
const TBound &node_j_bound = tree_nodes[j]->bound();
double min_squared_distance = node_i_bound.MinDistanceSq(node_j_bound);
double max_squared_distance = node_i_bound.MaxDistanceSq(node_j_bound);
distmat_.set(i, j, dmin);
distmat_.set(j, i, dmax);
distmat_.set(i, j, min_squared_distance);
distmat_.set(j, i, max_squared_distance);
}
}
return distmat_;
}
void Eval(const Matrix &data, const ArrayList<int> &indices, double *neg,
double *pos) {
/** @brief Computes the pairwise distance among three points.
*/
void EvalMinMaxSquaredDistances(const Matrix &data,
const ArrayList<index_t> &indices) {
for(index_t i = 0; i < indices.size(); i++) {
const double *i_col = data.GetColumnPtr(indices[i]);
for(index_t j = i + 1; j < indices.size(); j++) {
const double *j_col = data.GetColumnPtr(indices[j]);
distmat_.set(i, j, la::DistanceSqEuclidean(data.n_rows(), i_col,
j_col));
int num_order = order();
for(index_t i = 0; i < num_order - 1; i++) {
const double *point_i = data.GetColumnPtr(indices[i]);
for(index_t j = i + 1; j < num_order; j++) {
const double *point_j = data.GetColumnPtr(indices[j]);
double squared_distance = la::DistanceSqEuclidean(data.n_rows(),
point_i, point_j);
distmat_.set(i, j, squared_distance);
distmat_.set(j, i, squared_distance);
}
}
EvalUnnormOnSq(distmat_, neg, pos);
}
void EvalNodes(const ArrayList<DHrectBound<2> *> &node_bounds,
double *negmin, double *negmax, double *posmin,
double *posmax) {
/** @brief Computes $\frac{\nu}{r_i - r_j} \frac{\partial
* u}{\partial (r_i - r_j)}$, $\frac{\nu}{r_i - r_k}
* \frac{\partial u}{\partial (r_i - r_k)}$ and
* $\frac{\nu}{r_j - r_k} \frac{\partial u}{\partial (r_j -
* r_k)}$.
*/
void EvalGradients(const Matrix &dsqd_matrix,
double &min_negative_gradient1,
double *max_negative_gradient1,
double &min_positive_gradient1,
double *max_positive_gradient1,
double &min_negative_gradient2,
double *max_negative_gradient2,
double &min_positive_gradient2,
double *max_positive_gradient2,
double &min_negative_gradient3,
double *max_negative_gradient3,
double &min_positive_gradient3,
double *max_positive_gradient3) {
index_orders_[0] = 0;
index_orders_[1] = 1;
index_orders_[2] = 2;
gradient_(index_orders_, min_negative_gradient1, max_negative_gradient1,
min_positive_gradient1, max_positive_gradient1);
index_orders_[0] = 0;
index_orders_[1] = 2;
index_orders_[2] = 1;
gradient_(index_orders_, min_negative_gradient2, max_negative_gradient2,
min_positive_gradient2, max_positive_gradient2);
int num_nodes = node_bounds.size();
index_orders_[0] = 1;
index_orders_[1] = 2;
index_orders_[2] = 0;
gradient_(index_orders_, min_negative_gradient3, max_negative_gradient3,
min_positive_gradient3, max_positive_gradient3);
}
for(index_t i = 0; i < num_nodes - 1; i++) {
DHrectBound<2> *node_i_bound = node_bounds[i];
for(index_t j = i + 1; j < num_nodes; j++) {
DHrectBound<2> *node_j_bound = node_bounds[j];
double dmin = node_i_bound->MinDistanceSq(*node_j_bound);
double dmax = node_i_bound->MaxDistanceSq(*node_j_bound);
distmat_.set(i, j, dmin);
distmat_.set(j, i, dmax);
}
}
void EvalContributions
(const Matrix &data, const ArrayList<index_t> &indices,
double &negative_gradient1, double &positive_gradient1,
double &negative_gradient2, double &positive_gradient2,
double &negative_gradient3, double &positive_gradient3,
Vector &negative_force1_e, Vector &negative_force1_u,
Vector &positive_force1_l, Vector &positive_force1_e,
Matrix &negative_force2_e, Matrix &negative_force2_u,
Matrix &positive_force2_l, Matrix &positive_force2_e) {
index_orders_[0] = 0;
index_orders_[1] = 1;
index_orders_[2] = 2;
force_(data, indices, negative_gradient1, positive_gradient1,
negative_gradient2, positive_gradient2,
negative_gradient3, positive_gradient3,
negative_force1_e, negative_force1_u,
positive_force1_l, positive_force1_e,
negative_force2_e, negative_force2_u,
positive_force2_l, positive_force2_e);
EvalMinMax(negmin, negmax, posmin, posmax);
index_orders_[0] = 1;
index_orders_[1] = 0;
index_orders_[2] = 2;
force_(data, indices, negative_gradient1, positive_gradient1,
negative_gradient2, positive_gradient2,
negative_gradient3, positive_gradient3,
negative_force1_e, negative_force1_u,
positive_force1_l, positive_force1_e,
negative_force2_e, negative_force2_u,
positive_force2_l, positive_force2_e);
index_orders_[0] = 2;
index_orders_[1] = 1;
index_orders_[2] = 0;
force_(data, indices, negative_gradient2, positive_gradient2,
negative_gradient3, positive_gradient3,
negative_gradient1, positive_gradient1,
negative_force1_e, negative_force1_u,
positive_force1_l, positive_force1_e,
negative_force2_e, negative_force2_u,
positive_force2_l, positive_force2_e);
}
/** @brief Computes the first/second components of the
* negative/positive force components.
*/
void Eval(const Matrix &data, const ArrayList<index_t> &indices,
Vector &negative_force1_e, Vector &negative_force1_u,
Vector &positive_force1_l, Vector &positive_force1_e,
Matrix &negative_force2_e, Matrix &negative_force2_u,
Matrix &positive_force2_l, Matrix &positive_force2_e) {
double negative_gradient1, positive_gradient1;
double negative_gradient2, positive_gradient2;
double negative_gradient3, positive_gradient3;
// Evaluate the pairwise distances among all points.
EvalMinMaxSquaredDistances(data, indices);
// Evaluate the required components of the force vector.
EvalGradients(distmat_, negative_gradient1, NULL, positive_gradient1, NULL,
negative_gradient2, NULL, positive_gradient2, NULL,
negative_gradient3, NULL, positive_gradient3, NULL);
// Contributions to all three particles in the list.
EvalContributions(data, indices,
negative_gradient1, positive_gradient1,
negative_gradient2, positive_gradient2,
negative_gradient3, positive_gradient3,
negative_force1_e, negative_force1_u,
positive_force1_l, positive_force1_e,
negative_force2_e, negative_force2_u,
positive_force2_l, positive_force2_e);
}
};
fastlib2/contrib/dongryel/multibody/multibody_main.cc
+8 -24
View File
@@ -20,48 +20,32 @@ int main(int argc, char *argv[])
do_naive = fx_param_exists(NULL, "do_naive");
bandwidth = fx_param_double(NULL, "bandwidth", 0.1);
tau = fx_param_double(NULL, "tau", 0.1);
kernel = fx_param_str(NULL, "kernel", "gaussianthreebody");
kernel = fx_param_str(NULL, "kernel", "axilrodteller");
// Multibody computation
printf("Starting multitree multibody...\n");
if(!strcmp(kernel, "gaussianthreebody")) {
if(!strcmp(kernel, "axilrodteller")) {
fx_timer_start(NULL, "multibody");
MultitreeMultibody<GaussianThreeBodyKernel, GaussianKernelAux> mtmb;
MultitreeMultibody<AxilrodTellerForceKernel> mtmb;
mtmb.Init(bandwidth);
mtmb.Compute(tau);
fx_timer_stop(NULL, "multibody");
printf("multitree multibody completed...\n");
mtmb.PrintDebug();
// NAIVE
/*
if (do_naive) {
printf("Starting naive multibody...\n");
fx_timer_start(NULL, "naive_multibody");
NaiveMultibody<GaussianThreeBodyKernel> nmb;
nmb.Init(mtmb.get_data(), bandwidth);
nmb.Compute();
fx_timer_stop(NULL, "naive_multibody");
printf("finished naive multibody...\n");
}
}
else if(!strcmp(kernel, "axilrodteller")) {
fx_timer_start(NULL, "multibody");
MultitreeMultibody<AxilrodTellerKernel, GaussianKernelAux> mtmb;
mtmb.Init(bandwidth);
mtmb.Compute(tau);
fx_timer_stop(NULL, "multibody");
printf("multitree multibody completed...\n");
// NAIVE
if (do_naive) {
printf("Starting naive multibody...\n");
fx_timer_start(NULL, "naive_multibody");
NaiveMultibody<AxilrodTellerKernel> nmb;
NaiveMultibody<AxilrodTellerForceKernel> nmb;
nmb.Init(mtmb.get_data(), bandwidth);
nmb.Compute();
fx_timer_stop(NULL, "naive_multibody");
printf("finished naive multibody...\n");
}
*/
}
fx_done(NULL);
fastlib2/contrib/dongryel/multibody/multibody_stat.h
+151
View File
@@ -0,0 +1,151 @@
#ifndef INSIDE_MULTIBODY_H
#error "This is not a public header file!"
#endif
#ifndef MULTIBODY_STAT_H
#define MULTIBODY_STAT_H
class MultibodyStat {
public:
////////// Member Variables //////////
/** @brief The total coordinate sum.
*/
Vector coordinate_sum_;
/** @brief The maximum negative gradient (first component).
*/
double negative_gradient1_u;
/** @brief The minimum positive gradient (first component).
*/
double positive_gradient1_l;
/** @brief The maximum negative gradient (second component).
*/
Vector negative_gradient2_u;
/** @brief The minimum positive gradient (second component).
*/
Vector positive_gradient2_l;
/** @brief The postponed estimate of the first component of the
* negative gradient.
*/
double postponed_negative_gradient1_e;
/** @brief The postponed lower bound change to the first component
* of the negative gradient.
*/
double postponed_negative_gradient1_u;
/** @brief The postponed lower bound change to the first component
* of the positive gradient.
*/
double postponed_positive_gradient1_l;
/** @brief The postponed estimate of the first component of the
* positive gradient.
*/
double postponed_positive_gradient1_e;
/** @brief The postponed estimate of the second component of the
* negative gradient.
*/
Vector postponed_negative_gradient2_e;
/** @brief The postponed lower bound change to the second component
* of the negative gradient.
*/
Vector postponed_negative_gradient2_u;
/** @brief The postponed lower bound change to the second component
* of the positive gradient.
*/
Vector postponed_positive_gradient2_l;
/** @brief The postponed estimate of the second component of the
* positive gradient.
*/
Vector postponed_positive_gradient2_e;
/** @brief Resets the statistics to zero.
*/
void SetZero() {
negative_gradient1_u = 0;
positive_gradient1_l = 0;
negative_gradient2_u.SetZero();
positive_gradient2_l.SetZero();
postponed_negative_gradient1_e = 0;
postponed_negative_gradient1_u = 0;
postponed_positive_gradient1_l = 0;
postponed_positive_gradient1_e = 0;
postponed_negative_gradient2_e.SetZero();
postponed_negative_gradient2_u.SetZero();
postponed_positive_gradient2_l.SetZero();
postponed_positive_gradient2_e.SetZero();
}
/** @brief Initialize the statistics.
*/
void Init() {
coordinate_sum_.Init(3);
negative_gradient1_u = 0;
positive_gradient1_l = 0;
negative_gradient2_u.Init(3);
positive_gradient2_l.Init(3);
postponed_negative_gradient1_e = 0;
postponed_negative_gradient1_u = 0;
postponed_positive_gradient1_l = 0;
postponed_positive_gradient1_e = 0;
postponed_negative_gradient2_e.Init(3);
postponed_negative_gradient2_u.Init(3);
postponed_positive_gradient2_l.Init(3);
postponed_positive_gradient2_e.Init(3);
}
/** @brief The initialization for leaf stats.
*/
void Init(const Matrix& dataset, index_t &start, index_t &count) {
Init();
coordinate_sum_.SetZero();
for(index_t i = start; i < start + count; i++) {
la::AddTo(dataset.n_rows(), dataset.GetColumnPtr(i),
coordinate_sum_.ptr());
}
SetZero();
}
/** @brief The initialization for non-leaf stats based on the stats
* owned by the child nodes.
*/
void Init(const Matrix& dataset, index_t &start, index_t &count,
const MultibodyStat& left_stat, const MultibodyStat& right_stat) {
Init();
la::AddOverwrite(left_stat.coordinate_sum_, right_stat.coordinate_sum_,
&coordinate_sum_);
SetZero();
}
////////// Constructor/Destructor //////////
/** @brief The default constructor.
*/
MultibodyStat() {
}
/** @brief The default destructor.
*/
~MultibodyStat() {
}
};
#endif
fastlib2/contrib/dongryel/multibody/naive_multibody.h
+78
View File
@@ -0,0 +1,78 @@
#ifndef NAIVE_MULTIBODY_H
#define NAIVE_MULTIBODY_H
template<typename TMultibodyKernel>
class NaiveMultibody {
FORBID_ACCIDENTAL_COPIES(NaiveMultibody);
private:
/** Temporary space for storing indices selected for exhaustive computation
*/
ArrayList<int> exhaustive_indices_;
/** dataset for the tree */
Matrix data_;
/** multibody kernel function */
TMultibodyKernel mkernel_;
/** potential estimate */
double neg_potential_e_;
double pos_potential_e_;
/** exhaustive computer */
void NMultibody(int level) {
int num_nodes = mkernel_.order();
int start_index = 0;
double neg, pos;
if(level < num_nodes) {
if(level == 0) {
start_index = 0;
}
else {
start_index = exhaustive_indices_[level - 1] + 1;
}
for(index_t i = start_index; i < data_.n_cols() -
(num_nodes - level - 1); i++) {
exhaustive_indices_[level] = i;
NMultibody(level + 1);
}
}
else {
mkernel_.Eval(data_, exhaustive_indices_, &neg, &pos);
neg_potential_e_ += neg;
pos_potential_e_ += pos;
}
}
public:
NaiveMultibody() {}
~NaiveMultibody() {}
void Init(const Matrix &data, double bandwidth) {
data_.Alias(data);
exhaustive_indices_.Init(3);
mkernel_.Init(bandwidth);
neg_potential_e_ = pos_potential_e_ = 0;
}
void Compute() {
NMultibody(0);
printf("Negative potential sum %g\n", neg_potential_e_);
printf("Positive potential sum %g\n", pos_potential_e_);
printf("Got potential sum %g\n", neg_potential_e_ + pos_potential_e_);
}
};
#endif