From 5d129247586bb646007fda68d3a84d3d87e07eee Mon Sep 17 00:00:00 2001 From: Bill March Date: Mon, 7 Jun 2010 19:56:02 +0000 Subject: [PATCH] wrote multi-matcher with inefficient base case --- fastlib/trunk/contrib/march/npt/build.py | 14 + .../trunk/contrib/march/npt/multi_matcher.cc | 52 ++- .../trunk/contrib/march/npt/multi_matcher.h | 49 ++- .../trunk/contrib/march/npt/n_point_multi.cc | 348 ++++++++++++++---- .../trunk/contrib/march/npt/n_point_multi.h | 117 ++++++ .../contrib/march/npt/n_point_multi_main.cc | 40 ++ .../contrib/march/npt/n_point_perm_free.cc | 2 +- .../trunk/contrib/march/npt/results_tensor.cc | 102 ++++- .../trunk/contrib/march/npt/results_tensor.h | 153 +++++--- 9 files changed, 737 insertions(+), 140 deletions(-) create mode 100644 fastlib/trunk/contrib/march/npt/n_point_multi_main.cc diff --git a/fastlib/trunk/contrib/march/npt/build.py b/fastlib/trunk/contrib/march/npt/build.py index 8d1f4df429..d53b96529a 100644 --- a/fastlib/trunk/contrib/march/npt/build.py +++ b/fastlib/trunk/contrib/march/npt/build.py @@ -12,6 +12,20 @@ librule( deplibs = ["fastlib:fastlib", ":n_point_impl"] ) +librule( + name = "n_point_multi", + sources = ["n_point_multi.cc", "results_tensor.cc", "multi_matcher.cc", "n_point_nodes.cc"], + headers = ["n_point_multi.h", "results_tensor.h", "multi_matcher.h", "n_point_nodes.h"], + deplibs = ["fastlib:fastlib", ":n_point_impl"] +) + +binrule( + name = "n_point_multi_main", + sources = ["n_point_multi_main.cc"], + headers = [], + deplibs = ["fastlib:fastlib", ":n_point_multi"] +) + binrule( name = "n_point_testing", sources = ["n_point_testing.cc"], diff --git a/fastlib/trunk/contrib/march/npt/multi_matcher.cc b/fastlib/trunk/contrib/march/npt/multi_matcher.cc index 26a23316b9..3c96914d92 100644 --- a/fastlib/trunk/contrib/march/npt/multi_matcher.cc +++ b/fastlib/trunk/contrib/march/npt/multi_matcher.cc @@ -9,19 +9,53 @@ #include "multi_matcher.h" -// fills in the ranges of indices in the bandwidth -void MultiMatcher::FindBandwidths_(double min_dist_sq, double max_dist_sq, - double* max_subsume, double* min_exclude) { +bool MultiMatcher::TestPointPair(double dist_sq, index_t tuple_index_1, + index_t tuple_index_2, + ArrayList& permutation_ok, + ArrayList >& permutation_ranges) { + + bool this_point_works = false; + + DEBUG_ASSERT(tuple_index_1 < tuple_index_2); + + for (index_t perm_ind = 0; perm_ind < num_permutations_; perm_ind++) { + + // this permutation is already bad + if (! permutation_ok[perm_ind]) { + continue; + } + + if (dist_sq >= distances_[num_bins_ - 1]) { + + // this permutation is bad + permutation_ok[perm_ind] = false; + continue; + + } + + // figure out what the largest index that works is here + //permutation_ranges[perm_ind].set(tuple_index_1, tuple_index_2, dist_sq); + // TODO: do I need to set the other side of the diagonal? + + // TODO: double check this + double* ind_ptr = std::upper_bound(distances_.begin(), distances_.end(), + dist_sq); + int ind = (int)(ind_ptr - distances_.begin()); + permutation_ranges[perm_ind].set(tuple_index_1, tuple_index_2, ind); + + this_point_works = true; + + } // loop over permutations + + return this_point_works; -} // FindBandwidths_() +} // TestPointPair -// this needs to return the correct range of inconclusive bandwidths for the -// pair of nodes (i.e. everything inside is subsumed, everything outside is -// excluded) -// PROBLEM: the range for one permutation and that for another need not overlap + +/* void MultiMatcher::TestNodes_(const DHrectBound<2>& box1, const DHrectBound<2>& box2, index_t tuple_index_1, index_t tuple_index_2, @@ -40,5 +74,5 @@ void MultiMatcher::TestNodes_(const DHrectBound<2>& box1, range_out.Init(max_subsume, min_exclude); } // TestNodes_ - +*/ diff --git a/fastlib/trunk/contrib/march/npt/multi_matcher.h b/fastlib/trunk/contrib/march/npt/multi_matcher.h index 464f72b69b..09eaecfe43 100644 --- a/fastlib/trunk/contrib/march/npt/multi_matcher.h +++ b/fastlib/trunk/contrib/march/npt/multi_matcher.h @@ -10,8 +10,8 @@ #ifndef MULTI_MATCHER_H #define MULTI_MATCHER_H -#include "n_point_results.h" #include "fastlib/fastlib.h" +#include "n_point_impl.h" class MultiMatcher { @@ -20,16 +20,57 @@ private: ////////// variables ///////////// - ArrayList ranges_; + ArrayList distances_; + int num_bins_; Permutations perms_; + int num_permutations_; + + int tuple_size_; public: + index_t GetPermutationIndex_(index_t perm_index, index_t pt_index) { + + // these needed to be swapped to match matcher code + return perms_.GetPermutation(perm_index, pt_index); + + } // GetPermutation + int num_permutations() { + return num_permutations_; + } - void TestNodes_(const DHrectBound<2>& box1, const DHrectBound<2>& box2, - ResultsTensor& results); + int num_bins() { + return num_bins_; + } + + double max_dist() { + return distances_[num_bins_ - 1]; + } + + ArrayList& distances() { + return distances_; + } + + bool TestPointPair(double dist_sq, index_t tuple_index_1, + index_t tuple_index_2, + ArrayList& permutation_ok, + ArrayList >& permutation_ranges); + + + void Init(ArrayList& dists, int n) { + + tuple_size_ = n; + + // TODO: are these squared? + distances_.InitCopy(dists); + num_bins_ = distances_.size(); + + perms_.Init(tuple_size_); + num_permutations_ = perms_.num_perms(); + + } // Init() }; // class MultiMatcher diff --git a/fastlib/trunk/contrib/march/npt/n_point_multi.cc b/fastlib/trunk/contrib/march/npt/n_point_multi.cc index 20f4b922bf..76484d1833 100644 --- a/fastlib/trunk/contrib/march/npt/n_point_multi.cc +++ b/fastlib/trunk/contrib/march/npt/n_point_multi.cc @@ -10,7 +10,7 @@ #include "n_point_multi.h" -bool SymmetryCorrect_(ArrayList& nodes) { +bool NPointMulti::SymmetryCorrect_(ArrayList& nodes) { for (index_t i = 0; i < tuple_size_; i++) { @@ -31,7 +31,46 @@ bool SymmetryCorrect_(ArrayList& nodes) { } // SymmetryCorrect_ -index_t CheckBaseCase_(ArrayList& nodes) { +// fills inds with the indices in the range list that need to be recomputed +void NPointMulti::FindInvalidIndices_() { + + invalid_indices_.Init(tuple_size_); + + for (index_t split_ind = 0; split_ind < tuple_size_; split_ind++) { + + invalid_indices_[split_ind].Init(); + + // inserted this easy fix, not sure if the rest is right yet + if (tuple_size_ == 2) { + invalid_indices_[split_ind].PushBackCopy(0); + } + else { + + index_t bad_ind = split_ind - 1; + index_t bad_ind2 = 0; + + for (index_t i = 0; i < split_ind; i++) { + + invalid_indices_[split_ind].PushBackCopy(bad_ind); + bad_ind += tuple_size_ - 1 - (i+1); + bad_ind2 += tuple_size_ - i - 1; + + } // horizontal + + for (index_t i = split_ind+1; i < tuple_size_; i++) { + + invalid_indices_[split_ind].PushBackCopy(bad_ind2); + bad_ind2++; + + } + + } // n > 2 + + } // for split_ind + +} // FindInvalidIndices_() + +index_t NPointMulti::CheckBaseCase_(ArrayList& nodes) { index_t split_ind = -1; int split_size = 0; @@ -55,7 +94,7 @@ index_t CheckBaseCase_(ArrayList& nodes) { // needs to find the right index among (n choose 2) quantities, where i is // less than j -index_t FindInd_(index_t i, index_t j) { +index_t NPointMulti::FindInd_(index_t i, index_t j) { DEBUG_ASSERT(i < j); @@ -63,88 +102,253 @@ index_t FindInd_(index_t i, index_t j) { } // FindInd_() -void BaseCase_(ArrayList& nodes, ResultsTensor& tuple_status, - ResultsTensor& results) { +// returns true if the indices violate the symmetry requirement +bool NPointMulti::PointsViolateSymmetry_(index_t ind1, index_t ind2) { + DEBUG_ASSERT(ind1 >= 0); + DEBUG_ASSERT(ind2 >= 0); + return (ind2 <= ind1); +} // PointsViolateSymmetry_() + + + +void NPointMulti::BaseCaseHelper_(ArrayList >& point_sets, + ArrayList& permutation_ok, + ArrayList& points_in_tuple, + int k, + ArrayList >& permutation_ranges) { - // iterate over status tensor, skip things that are excluded + ArrayList permutation_ok_copy; + permutation_ok_copy.InitCopy(permutation_ok); - for (int i = 0; i < tensor_size; i++) { + ArrayList > permutation_ranges_copy; + permutation_ranges_copy.InitCopy(permutation_ranges); + + ArrayList k_rows; + k_rows.InitAlias(point_sets[k]); + + bool bad_symmetry = false; + + // loop over possible points for the kth member of the tuple + //for (index_t i = 0; !bad_symmetry && i < k_rows.size(); i++) { + // IMPORTANT: can't exit here for bad symmetry, it can get better as + // i increases + for (index_t i = 0; i < k_rows.size(); i++) { - if (! exclude) { - - matcher_.TestPointSets(); + index_t point_index_i = k_rows[i]; + + bool this_point_works = true; + + Vector point_i; + data_points_.MakeColumnVector(point_index_i, &point_i); + + // TODO: is this too inefficient? + permutation_ok_copy.Clear(); + permutation_ok_copy.AppendCopy(permutation_ok); + + permutation_ranges_copy.Clear(); + permutation_ranges_copy.AppendCopy(permutation_ranges); + + + // TODO: figure out a way to handle the bad symmetry more elegantly + // I should be able to avoid a bit more work + for (index_t j = 0; this_point_works && j < k; j++) { - } + index_t point_index_j = points_in_tuple[j]; + + // j should come before i since j comes first + bad_symmetry = PointsViolateSymmetry_(point_index_j, point_index_i); + //printf("point_j: %d, point_i: %d, bad_symmetry: %d\n", point_index_j, + // point_index_i, bad_symmetry); + + // don't compute the distances if we don't have to + if (!bad_symmetry) { + Vector point_j; + data_points_.MakeColumnVector(point_index_j, &point_j); + + double point_dist_sq = la::DistanceSqEuclidean(point_i, point_j); + + //printf("Testing point pair (%d, %d)\n", j, k); + // This needs to fill in the permutation_ok_copy for each matcher + this_point_works = matcher_.TestPointPair(point_dist_sq, j, k, + permutation_ok_copy, + permutation_ranges_copy); + //printf("this_point_works: %d\n", this_point_works); + + } // compute the distances and check the matcher + + } // for j - } + /* + printf("Considering point %d in position %d. bad_symmetry: %d, works: %d\n", + point_index_i, k, bad_symmetry, this_point_works); + */ + + // now, if the point passed, we put it in place and recurse + if (this_point_works && !bad_symmetry) { + + points_in_tuple[k] = point_index_i; + + // base case of the recursion + if (k == tuple_size_ - 1) { + + results_.ClearFilledResults(); + + for (index_t perm_index = 0; perm_index < matcher_.num_permutations(); + perm_index++) { + + // this one won't fit anywhere anyway + if (! permutation_ok_copy[perm_index]) { + continue; + } + + results_.IncrementRange(permutation_ranges[perm_index]); + + + } // iterate over permutations + + } // base case + else { + + BaseCaseHelper_(point_sets, permutation_ok_copy, + points_in_tuple, k+1, permutation_ranges_copy); + + } // recurse + + //DEBUG_ONLY(points_in_tuple[k] = -1); + + } // did the point work + + } // for i + + +} // BaseCaseHelper_() + + + +// Collect the indices of the valid matchers, check against each one? +// How to re-use info? +void NPointMulti::BaseCase_(NodeTuple& nodes, + ArrayList >& valid_ranges) { + + // Create the lists of points + ArrayList > point_sets; + point_sets.Init(tuple_size_); + + for (index_t i = 0; i < tuple_size_; i++) { + point_sets[i].Init(nodes.node_list(i)->count()); + + for (index_t j = 0; j < nodes.node_list(i)->count(); j++) { + point_sets[i][j] = j + nodes.node_list(i)->begin(); + } // for j + + } // for i + + + ArrayList > permutation_ranges; + permutation_ranges.Init(matcher_.num_permutations()); + + for (int i = 0; i < permutation_ranges.size(); i++) { + + permutation_ranges[i].Init(tuple_size_, tuple_size_); + // TODO: do I need to initialize it to some safe value? + + } // fill in permutation matrices + + + ArrayList permutation_ok; + permutation_ok.InitRepeat(true, matcher_.num_permutations()); + + ArrayList points_in_tuple; + points_in_tuple.InitRepeat(-1, tuple_size_); + + // TODO: figure out which matchers we need to worry about here + + BaseCaseHelper_(point_sets, permutation_ok, points_in_tuple, 0, + permutation_ranges); } // BaseCase_() -void DepthFirstRecursion_(ArrayList& nodes, - StatusTensor& status, - ResultsTensor& results) { +// valid_ranges are the ranges of indices in the distances_ array in the matcher +// it has length (n choose 2), the lower ends should be strictly non-decreasing +void NPointMulti::DepthFirstRecursion_(NodeTuple& nodes, + ArrayList >& valid_ranges) { + + bool can_prune = false; + + // update valid_ranges + // valid_ranges holds the range of distances that WON'T prune + // i.e. the only matchers that can't be pruned are ones that have a non-empty + // overlap with valid_ranges[i] for all i + for (index_t i = 0; i < valid_ranges.size(); i++) { - // check symmetry - if (!SymmetryCorrect_(nodes)) { + // IMPORTANT: first is lo, second is hi + + // TODO: how to account for upper and lower bounds in the matcher? + valid_ranges[i].second = min(valid_ranges[i].second, nodes.upper_bound(i)); + valid_ranges[i].first = max(valid_ranges[i].first, nodes.lower_bound(i)); + + if (valid_ranges[i].first >= valid_ranges[i].second) { + can_prune = true; + break; + } // check if the range is empty + + // TODO: make sure that it's not too small or large for any matcher + if (valid_ranges[i].first > matcher_.max_dist()) { + can_prune = true; + break; + } // too large + + // add lower bounds here later + + } // update ranges + + // check prune - i.e. check if it's still possible to contribute to anything + if (can_prune) { + num_total_prunes_++; return; - } - index_t split_ind = CheckBaseCase_(nodes); - - if (split_ind < 0) { - BaseCase_(nodes, status, results); - } + } // check prune + else if (nodes.all_leaves()) { + BaseCase_(nodes, valid_ranges); + } // base case else { + + NodeTuple left_node; + NodeTuple* left_node_ptr = &left_node; + + NodeTuple right_node; + NodeTuple* right_node_ptr = &right_node; + + // just pass in the invalid indices here + nodes.PerformSplit(left_node_ptr, right_node_ptr, invalid_indices_); + + // check if the list of bandwidths is still sorted here + + + if (left_node_ptr) { + + //printf("Left node\n"); + //left_node.Print(); + ArrayList > left_ranges; + left_ranges.InitCopy(valid_ranges); + DEBUG_ASSERT(left_node.node_list(0)); + DepthFirstRecursion_(left_node, left_ranges); + + } + if (right_node_ptr) { + + ArrayList > right_ranges; + right_ranges.InitCopy(valid_ranges); + //printf("Right node\n"); + //right_node.Print(); + DEBUG_ASSERT(right_node.node_list(0)); + DepthFirstRecursion_(right_node, right_ranges); + + } + + } // recurse - for (int perm_ind = 0; perm_ind < matcher_.num_permutations(); perm_ind++) { - - - ArrayList ranges; - ranges.Init(n_point_impl::NChooseR(tuple_size_, 2)); - - for (index_t i = 0; i < tuple_size_; i++) { - - NPointNode* node_i = nodes[i]; - - for (index_t j = i+1; j < tuple_size_; j++) { - - NPointNode* node_j = nodes[j]; - - DRange& range_ij = ranges[FindInd_(i, j)]; - matcher_.TestNodes_(node_i->bound(), node_j->bound(), i, j, perm_ind, - range_ij); - - } // for j - - } // for i - - // now, we have the ranges for this permutation, fill in the results status - - status.FillResults(ranges); - - } // for permutations - - // TODO: need to be able to tell if something subsumes, and perform it here - - // now, split and recurse - - NPointNode* split_node = nodes[split_ind]; - - ResultsTensor right_status; - right_status.Copy(status); - - nodes[split_ind] = split_node->left(); - - DepthFirstRecursion_(nodes, status, results); - - nodes[split_ind] = split_node->right(); - DepthFirstRecursion_(nodes, right_status, results); - - nodes[split_ind] = split_node; - - } // not a base case - -} // DepthFirstRecursion +} // DepthFirstRecursion_() diff --git a/fastlib/trunk/contrib/march/npt/n_point_multi.h b/fastlib/trunk/contrib/march/npt/n_point_multi.h index 255d291aa1..3338543ae1 100644 --- a/fastlib/trunk/contrib/march/npt/n_point_multi.h +++ b/fastlib/trunk/contrib/march/npt/n_point_multi.h @@ -14,6 +14,7 @@ #include "fastlib/fastlib.h" #include "n_point_impl.h" #include "results_tensor.h" +#include "n_point_nodes.h" class NPointMulti { @@ -24,18 +25,134 @@ private: MultiMatcher matcher_; + ResultsTensor results_; + fx_module* mod_; + Matrix data_points_; + int leaf_size_; + NPointNode* tree_; + + int tuple_size_; + + int num_total_prunes_; + + ArrayList > invalid_indices_; ///////////////// functions //////////////////// + bool SymmetryCorrect_(ArrayList& nodes); + + index_t CheckBaseCase_(ArrayList& nodes); + + index_t FindInd_(index_t i, index_t j); + + void BaseCaseHelper_(ArrayList >& point_sets, + ArrayList& permutations_ok, + ArrayList& points_in_tuple, + int k, + ArrayList >& permutation_ranges); + + + void BaseCase_(NodeTuple& nodes, + ArrayList >& valid_ranges); + + void DepthFirstRecursion_(NodeTuple& nodes, + ArrayList >& valid_ranges); + bool PointsViolateSymmetry_(index_t ind1, index_t ind2); + + void FindInvalidIndices_(); + + +public: + void Init(const Matrix& data, double band_min, double band_max, + int num_bands, int n, fx_module* mod) { + + mod_ = mod; + data_points_.Copy(data); + + tuple_size_ = n; + + leaf_size_ = fx_param_int(mod_, "leaf_size", 1); + + + // initialize results tensor + results_.Init(tuple_size_, num_bands); + + + // initialize matcher + ArrayList dists_sq; + dists_sq.Init(num_bands); + + double this_dist = band_min; + double dist_step = (band_max - band_min) / (double)num_bands; + + // TODO: double check this + for (index_t i = 0; i < num_bands; i++) { + + dists_sq[i] = this_dist * this_dist; + this_dist += dist_step; + + } + matcher_.Init(dists_sq, tuple_size_); + + ArrayList old_from_new; + tree_ = tree::MakeKdTreeMidpoint (data_points_, leaf_size_, + &old_from_new, NULL); + + num_total_prunes_ = 0; + + FindInvalidIndices_(); + + } // Init() + + void Compute() { + + fx_timer_start(mod_, "n_point_time"); + + NodeTuple nodes; + ArrayList node_list; + node_list.Init(tuple_size_); + + + for (index_t i = 0; i < tuple_size_; i++) { + + node_list[i] = tree_; + + } // for i + + nodes.Init(node_list); + + ArrayList > valid_ranges; + valid_ranges.Init(n_point_impl::NChooseR(tuple_size_, 2)); + + for (index_t i = 0; i < valid_ranges.size(); i++) { + valid_ranges[i].first = 0.0; + valid_ranges[i].second = DBL_MAX; + } + + DepthFirstRecursion_(nodes, valid_ranges); + + + fx_timer_stop(mod_, "n_point_time"); + + const char* filename = fx_param_str(mod_, "output_file", "output.txt"); + + FILE* fp; + fp = fopen(filename, "w"); + + results_.Output(matcher_.distances(), fp); + + fclose(fp); + + } // Compute }; // NPointMulti diff --git a/fastlib/trunk/contrib/march/npt/n_point_multi_main.cc b/fastlib/trunk/contrib/march/npt/n_point_multi_main.cc new file mode 100644 index 0000000000..bb2cd9abfd --- /dev/null +++ b/fastlib/trunk/contrib/march/npt/n_point_multi_main.cc @@ -0,0 +1,40 @@ +/* + * n_point_multi_main.cc + * + * + * Created by William March on 6/7/10. + * Copyright 2010 __MyCompanyName__. All rights reserved. + * + */ + + +#include "fastlib/fastlib.h" +#include "n_point_multi.h" + +int main(int argc, char* argv[]) { + + fx_init(argc, argv, NULL); + + Matrix data; + const char* data_file = fx_param_str_req(NULL, "data"); + data::Load(data_file, &data); + + double min_band, max_band; + min_band = fx_param_double_req(NULL, "min_band"); + max_band = fx_param_double_req(NULL, "max_band"); + int num_bands = fx_param_int_req(NULL, "num_bands"); + + int n = fx_param_int_req(NULL, "n"); + + fx_module* mod = fx_submodule(NULL, "n_point_multi"); + + NPointMulti alg; + alg.Init(data, min_band, max_band, num_bands, n, mod); + + alg.Compute(); + + fx_done(NULL); + + return 0; + +} // main() \ No newline at end of file diff --git a/fastlib/trunk/contrib/march/npt/n_point_perm_free.cc b/fastlib/trunk/contrib/march/npt/n_point_perm_free.cc index 961eb5ced1..b4d7f3b2b8 100644 --- a/fastlib/trunk/contrib/march/npt/n_point_perm_free.cc +++ b/fastlib/trunk/contrib/march/npt/n_point_perm_free.cc @@ -152,7 +152,7 @@ int NPointPermFree::BaseCaseHelper_(ArrayList >& point_sets, DEBUG_ONLY(points_in_tuple[k] = -1); - points_in_tuple[k] = -1; + //points_in_tuple[k] = -1; } // did the point work diff --git a/fastlib/trunk/contrib/march/npt/results_tensor.cc b/fastlib/trunk/contrib/march/npt/results_tensor.cc index b1f2ad6241..e791856787 100644 --- a/fastlib/trunk/contrib/march/npt/results_tensor.cc +++ b/fastlib/trunk/contrib/march/npt/results_tensor.cc @@ -7,17 +7,15 @@ * */ -#include "n_point_results.h" +#include "results_tensor.h" -// TODO: test me! +// This is the strictly upper triangular version +/* index_t ResultsTensor::FindIndex_(const ArrayList& indices) { - ArrayList sort_ind; - sort_ind.InitCopy(indices); - - // TODO: double check that this is correct - std::sort(sort_ind.begin(), sort_ind.end()); - + + // This should be unnecessary + // assuming the smallest index is first index_t return_ind = 0; for (index_t i = 1; i <= tensor_rank_; i++) { @@ -40,14 +38,94 @@ index_t ResultsTensor::FindIndex_(const ArrayList& indices) { } // FindIndex_ +*/ -// There is a range for each of the n! permutations -void ResultsTensor::FillRanges(ArrayList& ranges) { +// Full tensor version +index_t ResultsTensor::FindIndex_(const ArrayList& indices) { + + index_t result = indices[0]; + index_t power = tensor_rank_; + + for (index_t i = 1; i < indices.size(); i++) { + + result += indices[i] * power; + power = power * tensor_rank_; + + } // for i + + return result; + +} // FindIndex_() + +// TODO: test me! +bool ResultsTensor::IncrementIndex_(ArrayList& new_ind, + const ArrayList& orig_ind, + index_t k) { + + if (k >= tensor_rank_) { + return true; + } + + new_ind[k]++; + if (new_ind[k] >= lengths_) { + new_ind[k] = orig_ind[k]; + return IncrementIndex_(new_ind, orig_ind, k+1); + } + else { + return false; + } + +} // IncrementIndex + + +void ResultsTensor::IncrementRange(const GenMatrix& lower_inds) { + + // map row major into the index array + // TODO: would column major be more efficient? -} // FillRanges - + ArrayList this_result; + this_result.Init(lengths_); + + index_t row_ind = 0; + index_t col_ind = 1; + + for (index_t i = 0; i < this_result.size(); i++) { + + this_result[i] = lower_inds.get(row_ind, col_ind); + DEBUG_ASSERT(this_result[i] < lengths_ && this_result[i] >= 0); + col_ind++; + + if (col_ind >= lengths_) { + row_ind++; + col_ind = row_ind+1; + } // are we at the end of the row? + + } // for i + + ArrayList this_result_orig; + this_result_orig.InitCopy(this_result); + + bool done = false; + + while (!done) { + + // fill in the entry + + index_t this_ind = FindIndex_(this_result); + + if (!filled_results_[this_ind]) { + filled_results_[this_ind] = true; + results_[this_ind]++; + } + + // increment the array + done = IncrementIndex_(this_result, this_result_orig, 0); + + } // while + +} // IncrementRange() diff --git a/fastlib/trunk/contrib/march/npt/results_tensor.h b/fastlib/trunk/contrib/march/npt/results_tensor.h index 96f9e7d839..a4af40ec29 100644 --- a/fastlib/trunk/contrib/march/npt/results_tensor.h +++ b/fastlib/trunk/contrib/march/npt/results_tensor.h @@ -22,41 +22,43 @@ class ResultsTensor { private: int tensor_rank_; + int tuple_size_; - int num_bandwidths_; + int lengths_; - ArrayList bandwidths_; - - // TODO: how is this organized? ArrayList results_; + int num_results_; + + ArrayList filled_results_; //////////// functions ////////////////////// index_t FindIndex_(const ArrayList& indices); + bool IncrementIndex_(ArrayList& new_ind, + const ArrayList& orig_ind, + index_t k); + public: - void Init(int n, double min_band, double max_band, int num_bands) { - - DEBUG_ASSERT(max_band > min_band); - DEBUG_ASSERT(num_bands > 0); - + void Init(int n, int length) { + tuple_size_ = n; + tensor_rank_ = n_point_impl::NChooseR(n, 2); - num_bandwidths_ = num_bands; + lengths_ = length; - bandwidths_.Init(num_bandwidths_); + num_results_ = 1; + for (index_t i = 0; i < tensor_rank_; i++) { + num_results_ *= lengths_; + } - double bandwidth_step = (max_band - min_band) / (double)num_bandwidths_; - - for (index_t i = 0; i < num_bandwidths_; i++) { - - bandwidths_[i] = min_band + (double)i * bandwidth_step; - - } // fill in bandwidths - - results_.Init(n_point_impl::NChooseR(num_bandwidths_ + tuple_size_ + 1, - tuple_size_)); + // The strictly upper triangular version + //results_.InitRepeat(initial_result, + // n_point_impl::NChooseR(lengths_ + tuple_size_ + 1, + // tuple_size_)); + results_.InitRepeat(0, num_results_); + filled_results_.InitRepeat(false, num_results_); } // Init() @@ -75,39 +77,106 @@ public: results_[ind] = val; } // set() + + void SetAll(int val) { + + results_.Clear(); + results_.InitRepeat(val, num_results_); + + } - void AddTo(const ArrayList& indices, int val) { - - index_t ind = FindIndex_(indices); - - results_[ind] += val; - - } // AddTo() + /* + ArrayList& results() const { + return results_; + } + + ArrayList& filled_results() const { + return filled_results_; + } + */ + + int tensor_rank() { + return tensor_rank_; + } + + int lengths() { + return lengths_; + } + + void ClearFilledResults() { + for (index_t i = 0; i < num_results_; i++) { + filled_results_[i] = false; + } + } + + void IncrementRange(const GenMatrix& lower_inds); +/* void SetRange(const ArrayList& lower_ind, const ArrayList& upper_ind, int val); void AddToRange(const ArrayList& lower_ind, const ArrayList& upper_ind, int val); - - - void Print() { +*/ -// ArrayList indices; -// indices.InitRepeat(0, tuple_size_); + /* + void Copy(ResultsTensor& other) { + + results_.InitCopy(other.results()); + filled_results_.InitCopy(other.filled_results()); + + tensor_rank_ = other.tensor_rank(); + lengths_ = other.lengths(); + + } // Copy() + */ - // TODO: Make this format better - for (index_t i = 0; i < results_.size(); i++) { - - printf("%d\n", results_[i]); + void Output(ArrayList& distances_, FILE* fp) { + + ArrayList indices; + indices.InitRepeat(0, tensor_rank_); + + ArrayList indices_copy; + indices_copy.InitCopy(indices); + + bool done = false; + + while(!done) { - } // initialize + Matrix this_matcher; + this_matcher.Init(tuple_size_, tuple_size_); + + index_t row_ind = 0; + index_t col_ind = 1; + + for (index_t i = 0; i < indices.size(); i++) { + + this_matcher.set(row_ind, col_ind, distances_[indices_copy[i]]); + + row_ind++; + if (row_ind >= lengths_) { + row_ind = 0; + col_ind = row_ind + 1; + } + + } // fill in the matcher's matrix + + index_t ind = FindIndex_(indices_copy); + + int this_result = results_[ind]; + + // now do the printing + + this_matcher.PrintDebug("Matcher", fp); + fprintf(fp, "==Result: %d==\n\n", this_result); + + + done = IncrementIndex_(indices_copy, indices, 0); + + } // while - - - } // Print() - + } // Output() }; // NPointResults