From 1c843f20a4d2fa5011b7a52fdb0f5df076b40cf3 Mon Sep 17 00:00:00 2001 From: Dongryeol Lee Date: Fri, 21 Sep 2007 23:41:15 +0000 Subject: [PATCH] Added the basic finite difference scheme, now debugging and unit testing remains --- fastlib/u/dongryel/multibody/main.cc | 18 +- fastlib/u/dongryel/multibody/multibody.h | 238 ++++++++++++++++++++--- 2 files changed, 210 insertions(+), 46 deletions(-) diff --git a/fastlib/u/dongryel/multibody/main.cc b/fastlib/u/dongryel/multibody/main.cc index 6a4ce8db97..3c205d28b0 100644 --- a/fastlib/u/dongryel/multibody/main.cc +++ b/fastlib/u/dongryel/multibody/main.cc @@ -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 mtmb; - mtmb.Init(data); - mtmb.Compute(); + mtmb.Init(0.1); + mtmb.Compute(0.1); fx_timer_stop(NULL, "multibody"); // NAIVE diff --git a/fastlib/u/dongryel/multibody/multibody.h b/fastlib/u/dongryel/multibody/multibody.h index 0b8574b1e7..900982b97f 100644 --- a/fastlib/u/dongryel/multibody/multibody.h +++ b/fastlib/u/dongryel/multibody/multibody.h @@ -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, Matrix, MultibodyStat> Tree; + + typedef TKernel Kernel; + typedef TKernelDerivative KernelDerivative; + MultitreeMultibody() {} ~MultitreeMultibody() { delete root_; } /** Main computation */ - void Compute() { + void Compute(double tau) { ArrayList 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 tmp_nodes_; @@ -132,14 +174,40 @@ private: /** Storage for holding onto temporary leaf indices */ ArrayList tmp_leaf_indices_; + /** Temporary space for storing indices selected for exhaustive computation + */ + ArrayList 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 nodes) { + int Prunable(ArrayList 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 nodes) { - + void MTMultibodyBase(ArrayList 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 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); + } } }