checkin of thor kde (initial)

fastlib/u/dongryel/kde/build.py

@@ -35,9 +35,9 @@ librule(
     deplibs = ["fastlib:fastlib_int"]        # dependency
     )
 
-librule(
-    name = "thor_kde",
-    sources = [],
+binrule(
+    name = "thor_kde_bin",
+    sources = ["thor_kde_main.cc"],
     headers = ["thor_kde.h"],
     deplibs = ["fastlib:fastlib_int"]
     )
@@ -49,7 +49,6 @@ binrule(
     deplibs = [":fft_kde",
                ":fgt_kde",
                ":kde",
-               ":thor_kde",
                "u/dongryel/series_expansion:series_expansion",
                "fastlib:fastlib_int"]
     )

fastlib/u/dongryel/kde/main.cc

@@ -2,7 +2,6 @@
 #include "fft_kde.h"
 #include "fgt_kde.h"
 #include "kde.h"
-#include "thor_kde.h"
 
 int main(int argc, char *argv[]) {
 

fastlib/u/dongryel/kde/thor_kde_main.cc

@@ -0,0 +1,10 @@
+#include "fastlib/fastlib_int.h"
+#include "u/dongryel/series_expansion/kernel_aux.h"
+#include "thor_kde.h"
+
+int main(int argc, char *argv[]) {
+  fx_init(argc, argv);
+  ThorKde<GaussianKernel, GaussianKernelAux> thor_kde_;
+  fx_done();
+  return 0;
+}

fastlib/u/dongryel/series_expansion/farfield_expansion.h

@@ -90,6 +90,12 @@ class FarFieldExpansion {
 
   // interesting functions...
   
+  /**
+   * Accumulates the contribution of a single reference point as a 
+   * far-field moment.
+   */
+  void Accumulate(const Vector &v, double weight, int order);
+
   /**
    * Accumulates the far field moment represented by the given reference
    * data into the coefficients
@@ -184,6 +190,79 @@ class FarFieldExpansion {
 
 };
 
+template<typename TKernel, typename TKernelAux>
+void FarFieldExpansion<TKernel, TKernelAux>::Accumulate(const Vector &v,
+							double weight,
+							int order) {
+
+  int dim = v.length();
+  int total_num_coeffs = sea_->get_total_num_coeffs(order);
+  Vector tmp;
+  int r, i, j, k, t, tail;
+  Vector heads;
+  Vector x_r;
+  double bandwidth_factor = ka_.BandwidthFactor(kernel_.bandwidth_sq());
+
+  // initialize temporary variables
+  tmp.Init(total_num_coeffs);
+  heads.Init(dim + 1);
+  x_r.Init(dim);
+  Vector pos_coeffs;
+  Vector neg_coeffs;
+  pos_coeffs.Init(total_num_coeffs);
+  pos_coeffs.SetZero();
+  neg_coeffs.Init(total_num_coeffs);
+  neg_coeffs.SetZero();
+
+  // set to new order if greater
+  if(order_ < order) {
+    order_ = order;
+  }
+  Vector C_k;
+    
+  // Calculate the coordinate difference between the ref point and the 
+  // centroid.
+  for(i = 0; i < dim; i++) {
+    x_r[i] = (v[i] - center_[i]) / bandwidth_factor;
+  }
+  
+  // initialize heads
+  heads.SetZero();
+  heads[dim] = MAXINT;
+  
+  tmp[0] = 1.0;
+  
+  for(k = 1, t = 1, tail = 1; k <= order; k++, tail = t) {
+    for(i = 0; i < dim; i++) {
+      int head = (int) heads[i];
+      heads[i] = t;
+      
+      for(j = head; j < tail; j++, t++) {
+	tmp[t] = tmp[j] * x_r[i];
+      }
+    }
+  }
+  
+  // Tally up the result in A_k.
+  for(i = 0; i < total_num_coeffs; i++) {
+    double prod = weight * tmp[i];
+    
+    if(prod > 0) {
+      pos_coeffs[i] += prod;
+    }
+    else {
+      neg_coeffs[i] += prod;
+    }
+  }
+
+  // get multiindex factors
+  C_k.Alias(sea_->get_inv_multiindex_factorials());
+  
+  for(r = 0; r < total_num_coeffs; r++) {
+    coeffs_[r] += (pos_coeffs[r] + neg_coeffs[r]) * C_k[r];
+  }
+}
+
 template<typename TKernel, typename TKernelAux>
 void FarFieldExpansion<TKernel, TKernelAux>::AccumulateCoeffs
 (const Matrix& data, const Vector& weights, int begin, int end,