diff --git a/fastlib/u/garryb/superpar/new/allnn.cc b/fastlib/u/garryb/superpar/new/allnn.cc
index 40317d0a80..6bb1fa6c54 100644
--- a/fastlib/u/garryb/superpar/new/allnn.cc
+++ b/fastlib/u/garryb/superpar/new/allnn.cc
@@ -171,7 +171,7 @@ class Allnn {
 
     void FinishReaccumulate(const Param& param,
         const QNode& q_node) {
-      if (distance_sq_hi < 0) abort();
+      DEBUG_ASSERT_MSG(distance_sq_hi >= 0, "%f", distance_sq_hi);
       /* no post-processing steps necessary */
     }
 
diff --git a/fastlib/u/garryb/superpar/new/gnp.h b/fastlib/u/garryb/superpar/new/gnp.h
index ceb6b969a6..2cb229703b 100644
--- a/fastlib/u/garryb/superpar/new/gnp.h
+++ b/fastlib/u/garryb/superpar/new/gnp.h
@@ -89,10 +89,10 @@ void DualTreeDepthFirst<GNP>::Init(datanode *datanode) {
   do_naive_ = fx_param_bool(datanode, "do_naive", 0);
   
   Matrix q_matrix;
-  data::Load(fx_param_str_req(datanode, "q"), &q_matrix);
+  ASSERT_PASS(data::Load(fx_param_str_req(datanode, "q"), &q_matrix));
   
   Matrix r_matrix;
-  data::Load(fx_param_str_req(datanode, "r"), &r_matrix);
+  ASSERT_PASS(data::Load(fx_param_str_req(datanode, "r"), &r_matrix));
   
   param_.Init(fx_submodule(datanode, "param", "param"),
       q_matrix, r_matrix);
@@ -154,7 +154,9 @@ void DualTreeDepthFirst<GNP>::Begin() {
   }
   fx_timer_stop(datanode_, "execute");
   
-  //ot::Print(q_results_);
+  if (fx_param_bool(datanode_, "print", 0)) {
+    ot::Print(q_results_);
+  }
 }
 
 template<typename GNP>
@@ -166,7 +168,9 @@ void DualTreeDepthFirst<GNP>::PushDown_(index_t q_node_i) {
     for (index_t q_i = q_node->begin(); q_i < q_node->end(); q_i++) {
       typename GNP::QResult *q_result = &q_results_[q_i];
       typename GNP::Point *q_point = &q_tree_.points()[q_i];
+      typename GNP::QPointInfo *q_info = &q_tree_.point_info()[q_i];
       q_result->ApplyPostponed(param_, q_node_mut->postponed, *q_point);
+      q_result->Postprocess(param_, *q_point, *q_info, r_tree_.nodes()[0]);
     }
   } else {
     for (index_t k = 0; k < 2; k++) {
@@ -229,6 +233,7 @@ void DualTreeDepthFirst<GNP>::Pair_(index_t q_node_i, index_t r_node_i,
     }
 
     // Phase 2: Explore children, and reincorporate their results.
+    q_node_mut->postponed.Reset(param_);
     q_node_mut->mass_result.StartReaccumulate(param_, *q_node);
 
     for (index_t k = 0; k < 2; k++) {
@@ -248,7 +253,6 @@ void DualTreeDepthFirst<GNP>::Pair_(index_t q_node_i, index_t r_node_i,
     }
 
     q_node_mut->mass_result.FinishReaccumulate(param_, *q_node);
-    q_node_mut->postponed.Reset(param_);
   } else {
     index_t r_child1_i = r_node->child(0);
     index_t r_child2_i = r_node->child(1);
@@ -306,6 +310,8 @@ void DualTreeDepthFirst<GNP>::BaseCase_(
   visitor.Init(param_);
 
   q_node_mut->mass_result.StartReaccumulate(param_, *q_node);
+  DEBUG_ASSERT_MSG(q_node->count() != 0, "%d, count = %d, begin = %d",
+      q_node - qnode_(0), q_node->begin(), q_node->count());
 
   for (index_t q_i = q_node->begin(); q_i < q_node->end(); ++q_i) {
     typename GNP::Point *q_point = &q_tree_.points()[q_i];
@@ -331,6 +337,7 @@ void DualTreeDepthFirst<GNP>::BaseCase_(
     q_node_mut->mass_result.Accumulate(param_, *q_result);
   }
   q_node_mut->mass_result.FinishReaccumulate(param_, *q_node);
+  q_node_mut->postponed.Reset(param_);
 }
 
 
diff --git a/fastlib/u/garryb/superpar/new/kdtree.h b/fastlib/u/garryb/superpar/new/kdtree.h
index 94a98ce716..82c50abcc2 100644
--- a/fastlib/u/garryb/superpar/new/kdtree.h
+++ b/fastlib/u/garryb/superpar/new/kdtree.h
@@ -173,7 +173,7 @@ void KdTreeMidpointBuilder<TPointInfo, TNode, TParam>::KdTreeMidpointBuilder::Bu
           nodes_[node_i].begin(), nodes_[node_i].count(),
           &nodes_[left_i].bound(),
           &nodes_[right_i].bound());
-
+      
       DEBUG_MSG(3.0,"split (%d,[%d],%d) split_dim %d on %f (between %f, %f)",
           nodes_[node_i].begin(), split_col,
           nodes_[node_i].begin() + nodes_[node_i].count(), split_dim, split_val,
@@ -183,7 +183,7 @@ void KdTreeMidpointBuilder<TPointInfo, TNode, TParam>::KdTreeMidpointBuilder::Bu
       nodes_[left_i].Init(nodes_[node_i].begin(),
           split_col - nodes_[node_i].begin());
       nodes_[right_i].Init(split_col,
-          nodes_[node_i].begin() + nodes_[node_i].count() - split_col);
+          nodes_[node_i].end() - split_col);
 
       // This should never happen if max_width > 0
       DEBUG_ASSERT(nodes_[left_i].count() != 0 && nodes_[right_i].count() != 0);
diff --git a/fastlib/u/garryb/superpar/new/spbounds.h b/fastlib/u/garryb/superpar/new/spbounds.h
index f1a29cf8e2..0cf1556b21 100644
--- a/fastlib/u/garryb/superpar/new/spbounds.h
+++ b/fastlib/u/garryb/superpar/new/spbounds.h
@@ -268,7 +268,7 @@ class SpHrectBound {
       
       double v = (v1 + fabs(v1)) + (v2 + fabs(v2));
 
-      sumsq += math::PowAbs<t_pow, 1>(v * v);
+      sumsq += math::PowAbs<t_pow, 1>(v);
     }
 
     return math::Pow<2, t_pow>(sumsq) / 4;
diff --git a/fastlib/u/garryb/superpar/new/tkde.cc b/fastlib/u/garryb/superpar/new/tkde.cc
index accb634d09..c9c700a67c 100644
--- a/fastlib/u/garryb/superpar/new/tkde.cc
+++ b/fastlib/u/garryb/superpar/new/tkde.cc
@@ -47,8 +47,10 @@ class Tkde {
       kernel.Init(fx_param_double_req(datanode, "h"));
       double t = fx_param_double_req(datanode, "threshold");
       t = t * kernel.CalcNormConstant(dim);
-      thresh.lo = t * (1.0 - 1.0e-7);
-      thresh.hi = t * (1.0 + 1.0e-7);
+      fx_format_result(datanode, "norm_constant", "%f",
+          kernel.CalcNormConstant(dim));
+      thresh.lo = t * (1.0 - 1.0e-4);
+      thresh.hi = t * (1.0 + 1.0e-4);
       // WALDO: Fix me
     }
 
@@ -81,11 +83,13 @@ class Tkde {
     double ComputeKernelSum(
         double distance_squared,
         index_t r_count, const Vector& r_center, double r_sumsq) const {
+      //q*q - 2qr + rsumsq
+      //q*q - 2qr + r*r - r*r
       double quadratic_term =
           (distance_squared - la::Dot(r_center, r_center)) * r_count
           + r_sumsq;
 
-      return r_count - quadratic_term * kernel.inv_bandwidth_sq();
+      return -quadratic_term * kernel.inv_bandwidth_sq() + r_count;
     }
   };
 
@@ -273,10 +277,12 @@ class Tkde {
    public:
     void Init(const Param& param) {
       moment_info.Init(param);
+      label = LAB_UNKNOWN;
     }
 
     void Reset(const Param& param) {
       moment_info.Reset();
+      label = LAB_UNKNOWN;
     }
 
     void ApplyPostponed(const Param& param, const QPostponed& other) {
@@ -302,10 +308,6 @@ class Tkde {
     void Init(const Param& param) {
       d_density.Init(0, 0);
     }
-
-    void ApplyDelta(const Param& param, const Delta& other) {
-      d_density += other.d_density;
-    }
   };
 
   // rho
@@ -330,13 +332,19 @@ class Tkde {
     void Postprocess(const Param& param,
         const Vector& q_point, const QPointInfo& q_info,
         const RNode& r_root) {
-      /* nothing special to do */
+      if (density > param.thresh.hi) {
+        label |= LAB_HI;
+      } else if (density < param.thresh.lo) {
+        label |= LAB_LO;
+      }
+      DEBUG_ASSERT(label != LAB_CONFLICT);
     }
 
     void ApplyPostponed(const Param& param,
         const QPostponed& postponed,
         const Vector& q_point) {
       label |= postponed.label; /* bitwise OR */
+      DEBUG_ASSERT_MSG(label >= 0 && label < 3, "%d", label);
 
       if (!postponed.moment_info.is_empty()) {
         density += postponed.moment_info.ComputeKernelSum(param, q_point);
@@ -372,7 +380,7 @@ class Tkde {
     void Init(const Param& param) {
       /* horizontal init */
       density.Init(0, 0);
-      label = 0;
+      label = LAB_UNKNOWN;
     }
 
     void StartReaccumulate(const Param& param, const QNode& q_node) {
@@ -386,12 +394,14 @@ class Tkde {
       // but in some cases may require a copy/undo stage
       density |= result.density;
       label &= result.label;
+      DEBUG_ASSERT(result.label != LAB_CONFLICT);
     }
 
     void Accumulate(const Param& param,
         const QMassResult& result, index_t n_points) {
       density |= result.density;
-      density &= result.label;
+      label &= result.label;
+      DEBUG_ASSERT(result.label != LAB_CONFLICT);
     }
 
     void FinishReaccumulate(const Param& param,
@@ -403,8 +413,9 @@ class Tkde {
     void ApplyMassResult(const Param& param,
         const QMassResult& mass_result) {
       density += mass_result.density;
+      DEBUG_ASSERT_MSG((label | mass_result.label) != LAB_CONFLICT,
+          "%d and %d", label, mass_result.label);
       label |= mass_result.label;
-      DEBUG_ASSERT(label != LAB_CONFLICT);
     }
 
     void ApplyDelta(const Param& param,
@@ -417,6 +428,7 @@ class Tkde {
       bool change_made;
 
       if (unlikely(postponed.label)) {
+        DEBUG_ASSERT((label | postponed.label) != LAB_CONFLICT);
         label = postponed.label;
         change_made = true;
       } else  if (unlikely(!postponed.moment_info.is_empty())) {
@@ -520,24 +532,38 @@ class Tkde {
       double distance_sq_lo =
           q_node.bound().MinDistanceSqToBound(r_node.bound());
       bool need_expansion;
+      
+      //printf("%f %f %f\n",
+      //    q_node.bound().MinDistanceSqToBound(r_node.bound()),
+      //    q_node.bound().MidDistanceSqToBound(r_node.bound()),
+      //    q_node.bound().MaxDistanceSqToBound(r_node.bound())
+      //    );
 
-      if (distance_sq_lo >= param.kernel.bandwidth_sq()) {
+      if (distance_sq_lo > param.kernel.bandwidth_sq()) {
         need_expansion = false;
       } else {
         double distance_sq_hi =
             q_node.bound().MaxDistanceSqToBound(r_node.bound());
 
-        if (distance_sq_hi <= param.kernel.bandwidth_sq()) {
+        if (distance_sq_hi < param.kernel.bandwidth_sq()) {
           q_postponed->moment_info.Add(r_node.stat().moment_info);
           need_expansion = false;
         } else {
-          delta->d_density = r_node.stat().moment_info.ComputeKernelSumRange(
-              param, q_node.bound());
+#ifdef BIGPRUNE
+#endif
+          //delta->d_density = r_node.stat().moment_info.ComputeKernelSumRange(
+          //    param, q_node.bound());
           // we computed the lower bound of the quadratic.  if it is positive
           // it means we have a better-than-nothing bound; if it is not, then
           // we can resort to saying the min contribution is zero.
-          delta->d_density.lo = max(delta->d_density.lo, 0.0);
+          //    - problem: the upper bound is no good.
+          //max(delta->d_density.lo, 0.0);
+          delta->d_density.lo = 0;
+          delta->d_density.hi = r_node.count() *
+              param.kernel.EvalUnnormOnSq(distance_sq_lo);
           need_expansion = true;
+#ifdef BIGPRUNE
+#endif
         }
       }