wrote and tested single matcher

fastlib/branches/fastlib-stl/contrib/march/n_point/CMakeLists.txt

@@ -7,6 +7,8 @@ set(SOURCES
 	n_point_main.cc
 	permutations.h
 	permutations.cc
+	single_matcher.h
+	single_matcher.cc
 )
 
 # add directory name to sources

fastlib/branches/fastlib-stl/contrib/march/n_point/n_point_main.cc

@@ -8,7 +8,7 @@
  */
 
 // take this out later, just debugging with it
-#include "permutations.h"
+#include "single_matcher.h"
 
 using namespace npt;
 
@@ -22,8 +22,29 @@ int main(int argc, char* argv[]) {
   
   // output results
   
-  Permutations perms(3);
-  perms.Print();
+  
+  // testing stuff
+  
+  //Permutations perms(3);
+  //perms.Print();
+  
+  arma::mat lower_bds;
+  arma::mat upper_bds;
+  
+  lower_bds << 0.0 << 0.0 << 0.0 << arma::endr
+  << 0.0 << 0.0 << 0.0 << arma::endr
+  << 0.0 << 0.0 << 0.0 << arma::endr;
+  
+  upper_bds.load("test_upper_bds.csv");
+  
+  SingleMatcher matcher(3, lower_bds, upper_bds);
+  
+  std::vector<bool> perm_ok1(6, true);
+  std::cout << "Testing points: " << matcher.TestPointPair(0.5, 0, 1, perm_ok1)
+  << "\n";
+
+  //std::vector<bool> perm_ok2(6, true);
+  //std::cout << "Testing boxes: " << TestHrectPair(5.0, );
   
   return 0;
   

fastlib/branches/fastlib-stl/contrib/march/n_point/permutations.cc

@@ -9,6 +9,7 @@
 
 #include "permutations.h"
 
+// TODO: try a reference?
 void npt::Permutations::GeneratePermutations_(int k, int* perm_index,
                                               arma::Col<index_t>& trial_perm) {
 
@@ -17,14 +18,15 @@ void npt::Permutations::GeneratePermutations_(int k, int* perm_index,
     return;
   } 
 
-  // fill in all of trial perm
+  // Iterate over all points (i.e. everything that might be in the permutation)
   for (index_t i = 0; i < tuple_size_; i++) {
     
     bool perm_ok = true;
     
-    // check to see if i has already been used in this permutation
+    // Iterate over everything already in trial_perm
     for (index_t j = 0; perm_ok && j < k; j++) {
       
+      // Did we already use j in this one?  Then don't use it again.
       if (trial_perm[j] == i) {
         perm_ok = false;
       }
@@ -33,28 +35,27 @@ void npt::Permutations::GeneratePermutations_(int k, int* perm_index,
     
     // go to the next i if this one didn't work, otherwise proceed
     if (perm_ok) {
+      
+      // add i to the trial permutation
       trial_perm(k) = i;
     
-    
       // if the whole permutation is filled, put it in the matrix
       if (k == tuple_size_ - 1) {
         
-        // not sure if this works yet, needs to be tested
-        //permutation_indices_.insert_cols(permutation_indices_.n_cols, 
-        //                                 trial_perm);
         permutation_indices_.col(*perm_index) = trial_perm;
         (*perm_index)++;
         
       } // is the permutation filled?
       else {
-        
+        // move on to the next spot in the permutation
         GeneratePermutations_(k+1, perm_index, trial_perm);
         
-      } 
+      } // permutation not filled
         
     } // if perm_ok
     
   } // for i
   
   
-} // GeneratePermutations_
+} // GeneratePermutations_
+

fastlib/branches/fastlib-stl/contrib/march/n_point/permutations.h

@@ -33,6 +33,9 @@ namespace npt {
     // tuple_size_!
     index_t num_perms_;
     
+    
+    //////////////// functions //////////////////
+    
     // helper function at startup, actually forms all the permutations
     void GeneratePermutations_(int k, int* perm_index,
                                arma::Col<index_t>& trial_perm);
@@ -42,6 +45,9 @@ namespace npt {
     
   public:
     
+    // Dummy empty constructor
+    Permutations() {}
+    
     // The constructor needs to fill in permuation_indices_
     Permutations(index_t n) {
 
@@ -53,10 +59,11 @@ namespace npt {
         num_perms_ = num_perms_ * i;
       } // for i
       
+      // make the trial permutation
       arma::Col<index_t> trial_perm(tuple_size_);
       trial_perm.fill(-1);
       
-      
+      // allocate the matrix
       permutation_indices_.set_size(tuple_size_, num_perms_);
       permutation_indices_.fill(-1);
       //permutation_indices_(tuple_size_, 0);
@@ -71,6 +78,10 @@ namespace npt {
       
     } // constructor
     
+    int num_permutations() const {
+      return num_perms_;
+    }
+    
     // just accesses elements of permutation_indices_
     index_t GetPermutation(index_t perm_index, index_t point_index) const {
       

fastlib/branches/fastlib-stl/contrib/march/n_point/single_bandwidth_alg.h

@@ -11,10 +11,40 @@
  *  This is basically the same as the old auton code.
  */
 
+#ifndef SINGLE_BANDWIDTH_ALG_H
+#define SINGLE_BANDWIDTH_ALG_H
+
+#include "single_matcher.h"
+
 namespace npt {
 
   class SingleBandwidthAlg {
     
+  private:
+    
+    // the data and weights
+    arma::mat data_points_;
+    arma::colvec data_weights_;
+
+    // input params
+    index_t num_points_;
+    index_t tuple_size_;
+    
+    // the matcher
+    SingleMatcher matcher_;
+    
+    // the answer: num_tuples_ is the raw count, weighted is the sum of products
+    // of weights of all matching tuples
+    int num_tuples_;
+    double weighted_num_tuples_;
+    
+    // the number of times we pruned a tuple
+    int num_prunes_;
+    
+    NPointNode* tree_;
+    
+  public:
+    
     
     
     
@@ -23,4 +53,7 @@ namespace npt {
   
   
   
-} // namespace
+} // namespace
+
+#endif
+

fastlib/branches/fastlib-stl/contrib/march/n_point/single_matcher.cc

@@ -9,3 +9,95 @@
 
 #include "single_matcher.h"
 
+bool npt::SingleMatcher::CheckDistances_(double dist_sq, index_t ind1, 
+                                         index_t ind2) {
+  
+  return (dist_sq <= upper_bounds_sqr_(ind1, ind2) && 
+          dist_sq >= lower_bounds_sqr_(ind1, ind2));
+  
+} //CheckDistances_
+
+
+// note that this assumes that the points have been checked for symmetry
+bool npt::SingleMatcher::TestPointPair(double dist_sq, index_t tuple_ind_1, 
+                                       index_t tuple_ind_2,
+                                       std::vector<bool>& permutation_ok) {
+
+  bool any_matches = false;
+  
+  // iterate over all the permutations
+  for (index_t i = 0; i < num_permutations_; i++) {
+    
+    // did we already invalidate this one?
+    if (!(permutation_ok[i])) {
+      continue;
+    }
+    
+    index_t template_index_1 = GetPermIndex_(i, tuple_ind_1);
+    index_t template_index_2 = GetPermIndex_(i, tuple_ind_2);
+    
+    // Do the distances work?
+    if (CheckDistances_(dist_sq, template_index_1, template_index_2)) {
+      any_matches = true;
+    }
+    else {
+      permutation_ok[i] = false;
+    }
+    
+    // IMPORTANT: we can't exit here if any_matches is true
+    // This is because the ok permutation might get invalidated later, but we
+    // could still end up believing that unchecked ones are ok for this pair
+    
+  } // for i
+  
+  return any_matches;
+  
+} // TestPointPair
+
+// note that for now, there is no subsuming
+// this function will need to change if I want to add it
+bool npt::SingleMatcher::TestHrectPair(const DHrectBound<2>& box1, 
+                                       const DHrectBound<2>& box2,
+                                       index_t tuple_ind_1, index_t tuple_ind_2,
+                                       std::vector<bool>& permutation_ok) {
+  
+  bool any_matches = false;
+  
+  double max_dist_sq = box1.MaxDistanceSq(box2);
+  double min_dist_sq = box1.MinDistanceSq(box2);
+  
+  // iterate over all the permutations
+  for (index_t i = 0; i < num_permutations_; i++) {
+    
+    // did we already invalidate this one?
+    if (!(permutation_ok[i])) {
+      continue;
+    }
+    
+    index_t template_index_1 = GetPermIndex_(i, tuple_ind_1);
+    index_t template_index_2 = GetPermIndex_(i, tuple_ind_2);
+    
+    double upper_bound_sqr = upper_bounds_sqr_(template_index_1, 
+                                               template_index_2);
+    double lower_bound_sqr = lower_bounds_sqr_(template_index_1, 
+                                               template_index_2);
+    
+    // are they too far or too close?
+    if (max_dist_sq < lower_bound_sqr || min_dist_sq > upper_bound_sqr) {
+     
+      // this permutation doesn't work
+      permutation_ok[i] = false;
+      
+    }
+    else {
+     
+      // this permutation might work
+      any_matches = true;
+      
+    } // end if
+    
+  } // for i
+  
+  return any_matches;
+  
+} // TestHrectPair()

fastlib/branches/fastlib-stl/contrib/march/n_point/single_matcher.h

@@ -6,17 +6,78 @@
  *  Copyright 2011 __MyCompanyName__. All rights reserved.
  *
  *
- *  A single matcher.  
+ *  A single matcher.  This will use the same matrix formulation as before.
  *
  */
 
+#ifndef SINGLE_MATCHER_H
+#define SINGLE_MATCHER_H
+
+#include "permutations.h"
+
 namespace npt {
 
   class SingleMatcher {
     
+  private:
     
+    // stores the permutations, make sure to always reuse it instead of making
+    // more
+    Permutations perms_;
+    
+    // The upper and lower bounds for the matcher, stored in an upper triangular
+    // matrix.  They're squared to avoid dealing with square roots of distancess
+    arma::mat lower_bounds_sqr_;
+    arma::mat upper_bounds_sqr_;
+    
+    // n
+    index_t tuple_size_;
+    
+    // n!
+    index_t num_permutations_;
+    
+    /////////////////// functions ////////////////////////
+    
+    /**
+     * Just accesses the Permutations class
+     */
+    index_t GetPermIndex_(index_t perm_index, index_t pt_index) {
+      return perms_.GetPermutation(perm_index, pt_index);
+    } // GetPermIndex_
+    
+    
+    /**
+     * Helper function for checking points or bounds.
+     */
+    bool CheckDistances_(double dist_sq, index_t ind1, index_t ind2);
+    
+    
+  public:
+    
+    // Dummy empty constructor
+    //SingleMatcher() {}
+    
+    // constructor
+    SingleMatcher(index_t n, arma::mat& lower_bds, arma::mat& upper_bds) :
+        tuple_size_(n), perms_(n) {
+      
+      lower_bounds_sqr_ = arma::square(lower_bds);
+      upper_bounds_sqr_ = arma::square(upper_bds);
+      
+      num_permutations_ = perms_.num_permutations();
+      
+    } // constructor
+    
+    bool TestPointPair(double dist_sq, index_t tuple_ind_1, index_t tuple_ind_2,
+                       std::vector<bool>& permutation_ok);
+
+    bool TestHrectPair(const DHrectBound<2>& box1, const DHrectBound<2>& box2,
+                       index_t tuple_ind_1, index_t tuple_ind_2,
+                       std::vector<bool>& permutation_ok);
     
     
   }; // class
 
-} // namespace
+} // namespace
+
+#endif