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

Added the basic finite difference scheme, now debugging and unit testing remains

Browse Source
This commit is contained in:
Dongryeol Lee
2007-09-21 23:41:15 +00:00
parent 763bf0230e
commit 1c843f20a4
2 changed files with 210 additions and 46 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/dongryel/multibody/main.cc
+2 -16
View File
@@ -2,34 +2,20 @@
int main(int argc, char *argv[])
{
const char *datafile_name;
int leaflen;
Dataset dataset;
Matrix data;
bool do_naive;
fx_init(argc, argv);
// PARSE INPUTS
datafile_name = fx_param_str(NULL, "data", "small.arff");
leaflen = fx_param_int(NULL, "leaflen", 20);
do_naive = fx_param_exists(NULL, "do_naive");
// READING DATA
fx_timer_start(NULL, "read_d");
// read the dataset and get the matrix
if (!PASSED(dataset.InitFromFile(datafile_name))) {
fprintf(stderr, "main: Couldn't open file '%s'.\n", datafile_name);
return 1;
}
data.Alias(dataset.matrix());
fx_timer_stop(NULL, "read_d");
// Multibody computation
fx_timer_start(NULL,"multibody");
MultitreeMultibody<GaussianKernel, GaussianKernelDerivative> mtmb;
mtmb.Init(data);
mtmb.Compute();
mtmb.Init(0.1);
mtmb.Compute(0.1);
fx_timer_stop(NULL, "multibody");
// NAIVE
fastlib/u/dongryel/multibody/multibody.h
+208 -30
View File
@@ -18,11 +18,6 @@ public:
class MultibodyStat {
public:
/** lower index of the depth first order */
int lo_index;
/** high index of the depth first order */
int hi_index;
/** Summed up potential for query points in this node */
double potential_;
@@ -45,7 +40,13 @@ public:
/** Initialize the statistics */
void Init() {
}
void Init(double bandwidth, SeriesExpansionAux *sea) {
potential_ = 0;
extra_token_ = 0;
farfield_expansion_.Init(bandwidth, sea);
local_expansion_.Init(bandwidth, sea);
}
void Init(const Matrix& dataset, index_t &start, index_t &count) {
@@ -66,18 +67,26 @@ public:
farfield_expansion_.Init(bandwidth, center, sea);
local_expansion_.Init(bandwidth, center, sea);
}
MultibodyStat() { }
~MultibodyStat() {}
};
typedef BinarySpaceTree<DHrectBound<2>, Matrix, MultibodyStat> Tree;
typedef TKernel Kernel;
typedef TKernelDerivative KernelDerivative;
MultitreeMultibody() {}
~MultitreeMultibody() { delete root_; }
/** Main computation */
void Compute() {
void Compute(double tau) {
ArrayList<Tree *> root_nodes;
double total_n_tuples;
// Warning, I should fix this so that it generalizes to any number of
// tuples... May involves coming up with a solution that involves
@@ -88,7 +97,9 @@ public:
tmp_non_leaf_indices_.Init(0);
leaf_indices_.Init(0);
tmp_leaf_indices_.Init(0);
distmat_.Init(3, 3);
exhaustive_indices_.Init(3);
// store node pointers
for(index_t i = 0; i < 3; i++) {
root_nodes[i] = root_;
@@ -104,19 +115,50 @@ public:
}
}
total_n_tuples = ttn(0, root_nodes);
MTMultibody(root_nodes, total_n_tuples);
total_num_tuples_ = ttn(0, root_nodes);
tau_ = tau;
MTMultibody(root_nodes, total_num_tuples_);
printf("Total potential estimate: %g\n", potential_e_);
}
/** Initialize the tree */
void Init(Matrix& data) {
void InitExpansionObjects(Tree *node) {
if(node != NULL) {
Vector far_center;
Vector local_center;
far_center.Alias(node->stat().farfield_expansion_.get_center());
local_center.Alias(node->stat().local_expansion_.get_center());
node->bound().CalculateMidpoint(&far_center);
node->bound().CalculateMidpoint(&local_center);
node->stat().Init(sqrt(kernel_.bandwidth_sq()), &sea_);
}
if(!node->is_leaf()) {
InitExpansionObjects(node->left());
InitExpansionObjects(node->right());
}
}
/** Initialize the kernel object, and build the tree */
void Init(double bandwidth) {
fx_timer_start(NULL, "tree_d");
tree::LoadKdTree(NULL, &data, &root_, NULL);
fx_timer_stop(NULL, "tree_d");
tree::LoadKdTree(NULL, &data_, &root_, NULL);
sea_.Init(10, data_.n_rows());
kernel_.Init(bandwidth);
InitExpansionObjects(root_);
fx_timer_stop(NULL, "timer_d");
}
private:
/** Temporary storage space for holding onto the node pointers */
ArrayList<Tree *> tmp_nodes_;
@@ -132,14 +174,40 @@ private:
/** Storage for holding onto temporary leaf indices */
ArrayList<int> tmp_leaf_indices_;
/** Temporary space for storing indices selected for exhaustive computation
*/
ArrayList<int> exhaustive_indices_;
/** Temporary space for storing pairwise distances */
Matrix distmat_;
/** pointer to the root of the tree */
Tree *root_;
/** dataset for the tree */
Matrix data_;
/** series approximation auxiliary computations */
SeriesExpansionAux sea_;
/** kernel function */
Kernel kernel_;
/** the total number of n-tuples to consider */
double total_num_tuples_;
/** potential estimate */
double potential_e_;
/** Running lower bound on the potential */
double potential_l_;
/** approximation relative error bound */
double tau_;
bool as_indexes_strictly_surround_bs(Tree *a, Tree *b) {
return (a->stat().lo_index < b->stat().lo_index &&
a->stat().hi_index >= b->stat().hi_index) ||
(a->stat().lo_index <= b->stat().lo_index &&
a->stat().hi_index > b->stat().hi_index);
return (a->begin() < b->begin() && a->end() >= b->end()) ||
(a->begin() <= b->begin() && a->end() > b->end());
}
/**
@@ -177,10 +245,10 @@ private:
for(j = b+1 ; j < n && !conflict; j++) {
Tree *knj = nodes[j];
if (bkn->stat().lo_index >= knj->stat().hi_index) {
if (bkn->begin() >= knj->end() - 1) {
conflict = 1;
}
else if(nodes[j-1]->stat().hi_index > knj->stat().lo_index) {
else if(nodes[j-1]->end() - 1 > knj->begin()) {
simple_product = 0;
}
}
@@ -201,8 +269,8 @@ private:
for ( j = b+1 ; jdiff < 0 && j < n ; j++ ) {
Tree *knj = nodes[j];
if(bkn->stat().lo_index != knj->stat().lo_index ||
bkn->stat().hi_index != knj->stat().hi_index) {
if(bkn->begin() != knj->begin() ||
bkn->end() - 1 != knj->end() - 1) {
jdiff = j;
}
}
@@ -213,7 +281,7 @@ private:
else {
Tree *dkn = nodes[jdiff];
if(dkn->stat().lo_index >= bkn->stat().hi_index) {
if(dkn->begin() >= bkn->end() - 1) {
result = math::BinomialCoefficient(bkn->count(), jdiff - b);
if(result > 0.0) {
result *= ttn(jdiff, nodes);
@@ -250,30 +318,129 @@ private:
}
/** Pruning rule */
int Prunable(ArrayList<Tree *> nodes) {
int Prunable(ArrayList<Tree *> nodes, double num_tuples) {
int i, j;
double min_potential, max_potential;
double dsqd_ij_min, dsqd_ij_max, dsqd_ik_min, dsqd_ik_max, dsqd_jk_min,
dsqd_jk_max;
int num_nodes = nodes.size();
double lower_change;
double error, estimate;
double dmin = 0, dmax = 0;
// compute pairwise bounding box distances
for(i = 0; i < num_nodes - 1; i++) {
Tree *node_i = nodes[i];
for(j = i + 1; j < num_nodes; j++) {
Tree *node_j = nodes[j];
dmin = node_i->bound().MinDistanceSq(node_j->bound());
dmax = node_i->bound().MaxDistanceSq(node_j->bound());
distmat_.set(i, j, dmin);
distmat_.set(j, i, dmax);
}
}
dsqd_ij_min = distmat_.get(0, 1);
dsqd_ij_max = distmat_.get(1, 0);
dsqd_ik_min = distmat_.get(0, 2);
dsqd_ik_max = distmat_.get(2, 0);
dsqd_jk_min = distmat_.get(1, 2);
dsqd_jk_max = distmat_.get(2, 1);
min_potential = kernel_.EvalUnnormOnSq(dsqd_ij_max) *
kernel_.EvalUnnormOnSq(dsqd_ik_max) *
kernel_.EvalUnnormOnSq(dsqd_jk_max);
max_potential = kernel_.EvalUnnormOnSq(dsqd_ij_min) *
kernel_.EvalUnnormOnSq(dsqd_ik_min) *
kernel_.EvalUnnormOnSq(dsqd_jk_min);
lower_change = num_tuples * min_potential;
error = num_tuples * 0.5 * (max_potential - min_potential);
estimate = 0.5 * num_tuples * (min_potential + max_potential);
// compute whether the error is below the threshold
if(max_potential - min_potential <=
2 * tau_ * (potential_l_ + lower_change) / total_num_tuples_) {
potential_l_ += lower_change;
potential_e_ += estimate;
return 1;
}
return 0;
}
/** Base exhaustive case */
void MTMultibodyBase(ArrayList<Tree *> nodes) {
void MTMultibodyBase(ArrayList<Tree *> nodes, int level) {
int start_index;
double result;
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])->count(); i++) {
exhaustive_indices_[level] = i;
MTMultibodyBase(nodes, level + 1);
}
}
else {
/* complete the table of distance computation */
for(index_t i = 0; i < num_nodes; i++) {
const double *i_col = data_.GetColumnPtr(exhaustive_indices_[i]);
for(index_t j = i + 1; j < num_nodes; j++) {
const double *j_col = data_.GetColumnPtr(exhaustive_indices_[j]);
distmat_.set(i, j, la::DistanceSqEuclidean(data_.n_rows(), i_col,
j_col));
}
}
result = kernel_.EvalUnnormOnSq(distmat_.get(0, 1)) *
kernel_.EvalUnnormOnSq(distmat_.get(0, 2)) *
kernel_.EvalUnnormOnSq(distmat_.get(1, 2));
potential_e_ += result;
potential_l_ += result;
}
}
/** Main multitree recursion */
void MTMultibody(ArrayList<Tree *> nodes, double num_tuples) {
if(Prunable(nodes)) {
if(Prunable(nodes, num_tuples)) {
return;
}
// all leaves, then base case
else if(non_leaf_indices_.size() == 0) {
MTMultibodyBase(nodes);
MTMultibodyBase(nodes, 0);
}
// else, split an internal node and recurse
else {
int split_index;
double new_num_tuples;
tmp_non_leaf_indices_.Resize(0);
tmp_leaf_indices_.Resize(0);
@@ -292,7 +459,18 @@ private:
split_index = FindSplitNode(nodes);
nodes[split_index] = tmp_nodes_[split_index]->left();
new_num_tuples = ttn(0, nodes);
if(new_num_tuples > 0) {
MTMultibody(nodes, new_num_tuples);
}
nodes[split_index] = tmp_nodes_[split_index]->right();
new_num_tuples = ttn(0, nodes);
if(new_num_tuples > 0) {
MTMultibody(nodes, new_num_tuples);
}
}
}