n-point correlations, wrote basic kernel function

fastlib/trunk/contrib/march/npt/build.py

@@ -0,0 +1,13 @@
+librule(
+    name = "naive_two_point",
+    sources = [],
+    headers = ["naive_two_point.h"],
+    deplibs = ["fastlib:fastlib"]
+)
+
+binrule(
+    name = "n_point_testing",
+    sources = ["n_point_testing.cc"],
+    headers = [],
+    deplibs = ["fastlib:fastlib", ":naive_two_point"]
+)

fastlib/trunk/contrib/march/npt/matcher.cc

@@ -0,0 +1,59 @@
+/*
+ *  matcher.cc
+ *  
+ *
+ *  Created by William March on 2/23/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "matcher.h"
+
+
+bool Matcher::CheckDistances_(double dist_sq, index_t ind1, index_t ind2) {
+  
+  double upper_bd = upper_bounds_sqr_.ref(ind1, ind2);
+  double lower_bd = lower_bounds_sqr_.ref(ind1, ind2);
+  
+  return (dist_sq <= upper_bd && dist_sq >= lower_bd);
+  
+} // CheckDistances_()
+
+
+
+
+// IMPORTANT: this assumes that the two points have already been checked for 
+// symmetry
+bool Matcher::TestPointPair(double dist_sq, index_t tuple_index_1, 
+                            index_t tuple_index_2, 
+                            ArrayList<bool> permutation_ok) {
+
+  DEBUG_ASSERT(permutation_ok.size() == factorial(tuple_size_));
+
+  bool any_matches = false;
+  
+  for (index_t i = 0; i < num_permutations_; i++) {
+    
+    if (!(permutation_ok[i])) {
+      continue;
+    } // does this permutation work?
+    
+    index_t template_index_1 = GetPermutationIndex_(i, tuple_index_1);
+    index_t template_index_2 = GetPermutationIndex_(i, tuple_index_2);
+    
+   if (CheckDistances_(dist_sq, template_index_1, 
+                       template_index_2)) {
+     any_matches = true;
+   } // this placement works
+   else {
+     permutation_ok[i] = false;
+   } // this placement doesn't work
+    
+  } // for i
+  
+  
+  
+  
+  return any_matches;
+  
+} // TestPointPair()

fastlib/trunk/contrib/march/npt/matcher.h

@@ -0,0 +1,131 @@
+/*
+ *  matcher.h
+ *  
+ *
+ *  Created by William March on 2/23/10.
+ *
+ *
+ *  Stores the upper and lower bounds on each pair in the n-tuple.
+ *  Can also evaluate whether a set of points or nodes violates these 
+ *  conditions.
+ */
+
+// TODO: where should I keep the list of acceptable permutations?  
+// I think it should just be a vector passed down in the recursion
+
+class Matcher {
+  
+private:
+  
+  class Permutations {
+    
+  private:
+    
+    GenMatrix<index_t> permutation_indices_;
+    
+    int num_perms_;
+    
+  public:
+    
+    void Init(index_t n) {
+      
+      num_perms_ = factorials(n);
+      permutation_indices_.Init(num_perms_, n);
+      
+      // TODO: fill in the permutations
+
+      
+      
+    } // Init()
+    
+    int num_perms() {
+      return num_perms_;
+    }
+    
+    index_t GetPermutation(index_t perm_index, index_t pt_index) {
+      
+      return permutation_indices_.get(perm_index, pt_index);
+      
+    } // GetPermutation
+    
+  }; // class Permutations
+  
+  // these are symmetric n \times n matrices
+  // the diagonals are not defined, since they are never accessed.
+  // for a tuple to work, we need L_{i,j} \leq d(x_i, x_j) \leq H_{i,j}
+  // for all pairs (x_i, x_j) in the tuple (under some permutation)
+  // these are the squared bounds to prevent square roots in the code
+  Matrix lower_bounds_sqr_;
+  Matrix upper_bounds_sqr_;
+  
+  Permutations perms_;
+  
+  int tuple_size_;
+  int num_permutations_;
+  
+  ///////////////////// functions ///////////////////////
+  
+  /**
+   *
+   */
+  index_t GetPermutationIndex_(index_t perm_index, index_t pt_index) {
+    
+    return perms_.GetPermutation(perm_index, pt_index);
+    
+  } // GetPermutation
+  
+  /**
+   *
+   */
+  bool CheckDistances_(double dist_sq, index_t ind1, index_t ind2);
+
+  
+public:
+  
+  /**
+   *
+   */
+  void Init(const Matrix& lower, const Matrix& upper, int tuple_size) {
+    
+    lower_bounds_.Copy(lower);
+    upper_bounds_.Copy(upper);
+    tuple_size_ = tuple_size;
+    
+    // need to set them to be the squares of the values read in 
+    for (index_t i = 0; i < tuple_size_; i++) {
+      for (index_t j = 0; j < i; j++) {
+        
+        double new_low = lower_bounds_.get(i, j) * lower_bounds_.get(i, j);
+        double new_hi = upper_bounds_.get(i, j) * upper_bounds_.get(i, j);
+        
+        lower_bounds_.set(i, j, new_low);
+        lower_bounds_.set(j, i, new_low);
+        
+        upper_bound_.set(i, j, new_hi);
+        upper_bound_.set(j, i, new_hi);
+        
+      } // for j
+    } // for i
+    
+    perms_.Init(tuple_size_);
+    num_permutations_ = perms_.num_perms();
+    
+  } // Init()
+  
+  
+  /**
+   *
+   */
+  bool TestPointPair(double dist_sq, index_t tuple_index_1, 
+                     index_t tuple_index_2, ArrayList<bool> permutation_ok);
+  
+  
+  // TODO: how will I know what the nodes are inside this function?
+  // I think the Auton code just passes the hrects
+  /**
+   *
+   */
+  bool TestNodePair();
+  
+}; // Matcher
+

fastlib/trunk/contrib/march/npt/n_point_testing.cc

@@ -0,0 +1,51 @@
+/*
+ *  n_point_testing.cc
+ *  
+ *
+ *  Created by William March on 2/16/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "naive_two_point.h"
+#include "fastlib/fastlib.h"
+
+const fx_entry_doc n_point_testing_main_entries[] = {
+{"data", FX_REQUIRED, FX_STR, NULL,
+"A file containing the data.\n"},
+};
+
+const fx_submodule_doc n_point_testing_main_submodules[] = {
+{"naive_two_point_module", &naive_two_point_doc,
+"Naive two-point module.\n"},
+FX_SUBMODULE_DOC_DONE
+};
+
+const fx_module_doc n_point_testing_main_doc = {
+n_point_testing_main_entries, n_point_testing_main_submodules,
+"Runs and compares different methods for computing n point correlations.\n"
+};
+
+int main(int argc, char* argv[]) {
+  
+  fx_init(argc, argv, &n_point_testing_main_doc);
+  
+  const char* data_name = fx_param_str_req(NULL, "data");
+  
+  Matrix data;
+  data::Load(data_name, &data);
+  
+  //////////// Naive two point //////////////////
+  
+  fx_module* naive_two_point_mod = fx_submodule(NULL, "naive_two_point_module");
+  
+  NaiveTwoPoint naive_two_point_alg;
+  naive_two_point_alg.Init(data, naive_two_point_mod);
+
+  naive_two_point_alg.Compute();
+  
+  fx_done(NULL);
+  
+  return 0;
+  
+} 

fastlib/trunk/contrib/march/npt/naive_two_point.h

@@ -0,0 +1,105 @@
+/*
+ *  two_point.h
+ *  
+ *
+ *  Created by William March on 2/16/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "fastlib/fastlib.h"
+
+const fx_entry_doc naive_two_point_entries[] = {
+{"total_runtime", FX_TIMER, FX_CUSTOM, NULL, 
+"Total time required to compute the 2-point correlation.\n"},
+{"radius", FX_REQUIRED, FX_DOUBLE, NULL,
+  "The radius to compute the correlation for.\n"},
+{"num_pairs", FX_RESULT, FX_INT, NULL,
+  "The number of pairs found.\n"},
+FX_ENTRY_DOC_DONE
+};
+
+const fx_submodule_doc naive_two_point_submodules[] = {
+FX_SUBMODULE_DOC_DONE
+};
+
+const fx_module_doc naive_two_point_doc = {
+naive_two_point_entries, naive_two_point_submodules,
+"Algorithm module for naive serial two-point correlation.\n"
+};
+
+
+class NaiveTwoPoint {
+
+private:
+  
+  Matrix data_points_;
+
+  // TODO: replace with upper and lower bounds?
+  // TODO: replace with a range of values and compute the correlation for all
+  // of them?
+  double radius_;
+  
+  int num_pairs_;
+  index_t num_points_;
+  
+  fx_module* mod_;
+  
+public:
+  
+  void Init(const Matrix& data, fx_module* mod) {
+    
+    data_points_.Copy(data);
+    mod_ = mod;
+    
+    num_points_ = data_points_.n_cols();
+    
+    radius_ = fx_param_double_req(mod_, "radius");
+    
+    if (radius_ <= 0.0) {
+      FATAL("Negative radii not allowed.\n");
+    }
+    
+    num_pairs_ = 0;
+    
+  } // Init()
+  
+  
+  void Compute() {
+    
+    fx_timer_start(mod_, "total_runtime");
+
+    for (index_t i = 0; i < num_points_ - 1; i++) {
+      
+      Vector i_vec;
+      data_points_.MakeColumnVector(i, &i_vec);
+      
+      for (index_t j = i + 1; j < num_points_; j++) {
+        
+        Vector j_vec;
+        data_points_.MakeColumnVector(j, &j_vec);
+        
+        double dist = sqrt(la::DistanceSqEuclidean(i_vec, j_vec));
+        
+        //printf("distance = %g\n", dist);
+        
+        if (dist < radius_) {
+          
+          num_pairs_++;
+          
+        } // is dist small
+        
+      } // for j
+      
+    } // for i
+    
+    fx_timer_stop(mod_, "total_runtime");
+    
+    fx_result_int(mod_, "num_pairs", num_pairs_);
+    
+    printf("Number of pairs: %d\n\n", num_pairs_);
+    
+  } // Compute()
+  
+  
+}; // NaiveTwoPoint

fastlib/trunk/contrib/march/npt/test.csv

@@ -0,0 +1,4 @@
+0, 0, 0
+1.0, 0, 0
+2.1, 0, 0
+3.2, 0, 0