wrote multi-matcher with inefficient base case
This commit is contained in:
Bill March
@@ -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"],
|
||||
|
||||
@@ -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<bool>& permutation_ok,
|
||||
ArrayList<GenMatrix<index_t> >& 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_
|
||||
|
||||
*/
|
||||
|
||||
|
||||
@@ -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<double> ranges_;
|
||||
ArrayList<double> 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<double>& distances() {
|
||||
return distances_;
|
||||
}
|
||||
|
||||
bool TestPointPair(double dist_sq, index_t tuple_index_1,
|
||||
index_t tuple_index_2,
|
||||
ArrayList<bool>& permutation_ok,
|
||||
ArrayList<GenMatrix<index_t> >& permutation_ranges);
|
||||
|
||||
|
||||
void Init(ArrayList<double>& 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
|
||||
|
||||
|
||||
@@ -10,7 +10,7 @@
|
||||
#include "n_point_multi.h"
|
||||
|
||||
|
||||
bool SymmetryCorrect_(ArrayList<NPointNode*>& nodes) {
|
||||
bool NPointMulti::SymmetryCorrect_(ArrayList<NPointNode*>& nodes) {
|
||||
|
||||
for (index_t i = 0; i < tuple_size_; i++) {
|
||||
|
||||
@@ -31,7 +31,46 @@ bool SymmetryCorrect_(ArrayList<NPointNode*>& nodes) {
|
||||
|
||||
} // SymmetryCorrect_
|
||||
|
||||
index_t CheckBaseCase_(ArrayList<NPointNode*>& 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<NPointNode*>& nodes) {
|
||||
|
||||
index_t split_ind = -1;
|
||||
int split_size = 0;
|
||||
@@ -55,7 +94,7 @@ index_t CheckBaseCase_(ArrayList<NPointNode*>& 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<NPointNode*>& 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<ArrayList<index_t> >& point_sets,
|
||||
ArrayList<bool>& permutation_ok,
|
||||
ArrayList<index_t>& points_in_tuple,
|
||||
int k,
|
||||
ArrayList<GenMatrix<index_t> >& permutation_ranges) {
|
||||
|
||||
// iterate over status tensor, skip things that are excluded
|
||||
ArrayList<bool> permutation_ok_copy;
|
||||
permutation_ok_copy.InitCopy(permutation_ok);
|
||||
|
||||
for (int i = 0; i < tensor_size; i++) {
|
||||
ArrayList<GenMatrix<index_t> > permutation_ranges_copy;
|
||||
permutation_ranges_copy.InitCopy(permutation_ranges);
|
||||
|
||||
ArrayList<index_t> 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<std::pair<double, double> >& valid_ranges) {
|
||||
|
||||
// Create the lists of points
|
||||
ArrayList<ArrayList<index_t> > 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<GenMatrix<index_t> > 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<bool> permutation_ok;
|
||||
permutation_ok.InitRepeat(true, matcher_.num_permutations());
|
||||
|
||||
ArrayList<index_t> 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<NPointNode*>& 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<std::pair<double, double> >& 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<std::pair<double, double> > left_ranges;
|
||||
left_ranges.InitCopy(valid_ranges);
|
||||
DEBUG_ASSERT(left_node.node_list(0));
|
||||
DepthFirstRecursion_(left_node, left_ranges);
|
||||
|
||||
}
|
||||
if (right_node_ptr) {
|
||||
|
||||
ArrayList<std::pair<double, double> > 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<DRange> 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_()
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -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<ArrayList<index_t> > invalid_indices_;
|
||||
|
||||
|
||||
///////////////// functions ////////////////////
|
||||
|
||||
|
||||
bool SymmetryCorrect_(ArrayList<NPointNode*>& nodes);
|
||||
|
||||
index_t CheckBaseCase_(ArrayList<NPointNode*>& nodes);
|
||||
|
||||
index_t FindInd_(index_t i, index_t j);
|
||||
|
||||
void BaseCaseHelper_(ArrayList<ArrayList<index_t> >& point_sets,
|
||||
ArrayList<bool>& permutations_ok,
|
||||
ArrayList<index_t>& points_in_tuple,
|
||||
int k,
|
||||
ArrayList<GenMatrix<index_t> >& permutation_ranges);
|
||||
|
||||
|
||||
void BaseCase_(NodeTuple& nodes,
|
||||
ArrayList<std::pair<double, double> >& valid_ranges);
|
||||
|
||||
void DepthFirstRecursion_(NodeTuple& nodes,
|
||||
ArrayList<std::pair<double, double> >& 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<double> 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<index_t> old_from_new;
|
||||
tree_ = tree::MakeKdTreeMidpoint<NPointNode> (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<NPointNode*> 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<std::pair<double, double> > 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
|
||||
|
||||
@@ -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()
|
||||
@@ -152,7 +152,7 @@ int NPointPermFree::BaseCaseHelper_(ArrayList<ArrayList<index_t> >& point_sets,
|
||||
|
||||
DEBUG_ONLY(points_in_tuple[k] = -1);
|
||||
|
||||
points_in_tuple[k] = -1;
|
||||
//points_in_tuple[k] = -1;
|
||||
|
||||
} // did the point work
|
||||
|
||||
|
||||
@@ -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<index_t>& indices) {
|
||||
|
||||
ArrayList<index_t> 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<index_t>& indices) {
|
||||
|
||||
} // FindIndex_
|
||||
|
||||
*/
|
||||
|
||||
// There is a range for each of the n! permutations
|
||||
void ResultsTensor::FillRanges(ArrayList<DRange>& ranges) {
|
||||
// Full tensor version
|
||||
index_t ResultsTensor::FindIndex_(const ArrayList<index_t>& 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<index_t>& new_ind,
|
||||
const ArrayList<index_t>& 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<index_t>& lower_inds) {
|
||||
|
||||
// map row major into the index array
|
||||
// TODO: would column major be more efficient?
|
||||
|
||||
|
||||
|
||||
} // FillRanges
|
||||
|
||||
ArrayList<index_t> 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<index_t> 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()
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -22,41 +22,43 @@ class ResultsTensor {
|
||||
private:
|
||||
|
||||
int tensor_rank_;
|
||||
int tuple_size_;
|
||||
|
||||
int num_bandwidths_;
|
||||
int lengths_;
|
||||
|
||||
ArrayList<double> bandwidths_;
|
||||
|
||||
// TODO: how is this organized?
|
||||
ArrayList<int> results_;
|
||||
int num_results_;
|
||||
|
||||
ArrayList<bool> filled_results_;
|
||||
|
||||
//////////// functions //////////////////////
|
||||
|
||||
index_t FindIndex_(const ArrayList<index_t>& indices);
|
||||
|
||||
bool IncrementIndex_(ArrayList<index_t>& new_ind,
|
||||
const ArrayList<index_t>& 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<index_t>& indices, int val) {
|
||||
|
||||
index_t ind = FindIndex_(indices);
|
||||
|
||||
results_[ind] += val;
|
||||
|
||||
} // AddTo()
|
||||
/*
|
||||
ArrayList<int>& results() const {
|
||||
return results_;
|
||||
}
|
||||
|
||||
ArrayList<bool>& 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<index_t>& lower_inds);
|
||||
|
||||
|
||||
/*
|
||||
void SetRange(const ArrayList<index_t>& lower_ind,
|
||||
const ArrayList<index_t>& upper_ind, int val);
|
||||
|
||||
void AddToRange(const ArrayList<index_t>& lower_ind,
|
||||
const ArrayList<index_t>& upper_ind, int val);
|
||||
|
||||
|
||||
void Print() {
|
||||
*/
|
||||
|
||||
// ArrayList<index_t> 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<double>& distances_, FILE* fp) {
|
||||
|
||||
ArrayList<index_t> indices;
|
||||
indices.InitRepeat(0, tensor_rank_);
|
||||
|
||||
ArrayList<index_t> 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
|
||||
|
||||
Reference in New Issue
Block a user