diff --git a/fastlib/u/garryb/nbr/affinity.cc b/fastlib/u/garryb/nbr/affinity.cc
index 4c896e8679..c32a36c1dc 100644
--- a/fastlib/u/garryb/nbr/affinity.cc
+++ b/fastlib/u/garryb/nbr/affinity.cc
@@ -1,8 +1,17 @@
+#include "spbounds.h"
+#include "gnp.h"
+#include "dfs.h"
+#include "nbr_utils.h"
+
+#include "fastlib/fastlib.h"
 
 struct AffinityCommon {
   /** The bounding type. Required by NBR. */
   typedef SpHrectBound<2> Bound;
 
+  /**
+   * Alpha corresponds to "maximum availability" with the != k condition.
+   */
   struct Alpha {
    public:
     double max1;
@@ -16,7 +25,7 @@ struct AffinityCommon {
     }
 
    public:
-    double get(index_t i) {
+    double get(index_t i) const {
       if (unlikely(i == max1_index)) {
         return max2;
       } else {
@@ -25,24 +34,19 @@ struct AffinityCommon {
     }
   };
 
-  struct AlphaInfo {
+  struct CombinedInfo {
+    /** Maximum availability of the point. */
     Alpha alpha;
-
-    OT_DEF(AlphaInfo) {
-      OT_MY_OBJECT(alpha);
-    }
-  };
-  
-  struct RhoInfo {
+    /** Sum of responsibilities. */
     double rho;
 
-    OT_DEF(RhoInfo) {
+    OT_DEF(CombinedInfo) {
+      OT_MY_OBJECT(alpha);
       OT_MY_OBJECT(rho);
     }
   };
-  
-  typedef SpVectorPoint<AlphaInfo> AlphaPoint;
-  typedef SpVectorPoint<RhoInfo> RhoPoint;
+
+  typedef SpVectorInfoPoint<CombinedInfo> CombinedPoint;
 
   struct Param {
    public:
@@ -52,130 +56,144 @@ struct AffinityCommon {
     index_t dim;
     /** Number of points */
     index_t n_points;
+    /** Self-pereference. */
+    double pref;
+    /** The damping factor. */
+    double lambda;
 
     OT_DEF(Param) {
       OT_MY_OBJECT(eps);
       OT_MY_OBJECT(dim);
       OT_MY_OBJECT(n_points);
+      OT_MY_OBJECT(pref);
+      OT_MY_OBJECT(lambda);
     }
 
    public:
     void Copy(const Param& other) {
+      pref = other.pref;
       eps = other.eps;
       dim = other.dim;
+      n_points = other.n_points;
+      lambda = other.lambda;
     }
 
     void Init(datanode *module) {
       dim = -1;
       eps = fx_param_double(module, "eps", 1.0e-2);
+      pref = fx_param_double_req(module, "pref");
+      lambda = fx_param_double(module, "lambda", 0.8);
     }
 
-    void AnalyzePoint(const AlphaPoint& q) {
-      if (dim == -1) {
-        dim = q.vec().length();
-      } else {
-        DEBUG_ASSERT_MSG(dim == q.length(), "Differing dimensionality");
-      }
+    void BootstrapMonochromatic(CombinedPoint *point, index_t count) {
+      dim = point->vec().length();
+      n_points = count;
+      // TODO: Realistic values
+      point->info().rho = 0;
+      point->info().alpha.max1 = 0;
+      point->info().alpha.max2 = 0;
+      point->info().alpha.max1_index = 0;
     }
-    void AnalyzePoint(const RhoPoint& r) {
-      if (dim == -1) {
-        dim = q.vec().length();
-      } else {
-        DEBUG_ASSERT_MSG(dim == q.length(), "Differing dimensionality");
-      }
+  };
+
+  struct CombinedStat {
+   public:
+    SpRange alpha;
+    SpRange rho;
+
+    OT_DEF(CombinedStat) {
+      OT_MY_OBJECT(alpha);
+      OT_MY_OBJECT(rho);
     }
 
    public:
-    double Similarity(const Vector& q, const Vector& r) const {
-      return 1.0 / sqrt(la::EuclideanDistanceSq(a, b));
+    void Init(const Param& param) {
+      Reset(param);
     }
-    double Similarity(
+    void Reset(const Param& param) {
+      alpha.InitEmptySet();
+      rho.InitEmptySet();
+    }
+    void Accumulate(const Param& param, const CombinedPoint& point) {
+      alpha |= SpRange(point.info().alpha.max2, point.info().alpha.max1);
+      rho |= point.info().rho;
+    }
+    void Accumulate(const Param& param,
+        const CombinedStat& stat, const Bound& bound, index_t n) {
+      alpha |= stat.alpha;
+      rho |= stat.rho;
+    }
+    void Postprocess(const Param& param, const Bound& bound, index_t n) {}
+  };
+
+  typedef SpNode<Bound, CombinedStat> CombinedNode;
+
+  struct Helpers {
+    static double Similarity(double distsq) {
+      return -distsq;
+    }
+    static double Similarity(const Vector& a, const Vector& b) {
+      //uint32 anum = (mem::PointerAbsoluteAddress(a.ptr()) * 315187727);
+      //uint32 bnum = (mem::PointerAbsoluteAddress(b.ptr()) * 210787727);
+      //uint32 val = ((anum ^ bnum) >> 16) & 0xfff;
+      //double noise = (1.0e-6 / 4096) * val;
+      return Similarity(la::DistanceSqEuclidean(a, b));
+    }
+    static double Similarity(
+        const Param& param,
         const Vector& q, index_t q_index,
-        const Vector& r, index_t r_index) const {
+        const Vector& r, index_t r_index) {
       if (unlikely(q_index == r_index)) {
-        return pref;
+        return param.pref;
+      } else {
+        return Similarity(q, r);
       }
-      return Similarity(q, r);
     }
-    double SimilarityHi(const AlphaNode& a, const AlphaNode& b) const {
-      double dist = sqrt(a->bound().MinDistanceSqToBound(b->bound()));
-      double hi = 1.0 / dist;
-      if (q->begin() < r->end() && r->begin() < q->end() && pref > upper_bound) {
-        hi = pref;
+    static double SimilarityHi(
+        const Param& param,
+        const CombinedNode& a, const CombinedNode& b) {
+      double distsq = a.bound().MinDistanceSqToBound(b.bound());
+      double hi = Similarity(distsq);
+      if (a.begin() < b.end() && b.begin() < a.end()
+          && param.pref > hi) {
+        hi = param.pref;
       }
       return hi;
     }
+    static double SimilarityLo(
+        const Param& param,
+        const CombinedNode& a, const CombinedNode& b) {
+      double distsq = a.bound().MaxDistanceSqToBound(b.bound());
+      double lo = Similarity(distsq);
+      if (a.begin() < b.end() && b.begin() < a.end()
+          && param.pref < lo) {
+        lo = param.pref;
+      }
+      return lo;
+    }
     
-    double ErrorShare(double abs_error_used, const RNode& r_node) {
-      return (eps - abs_error_used) * r_node.count() / n_points;
+    static double ErrorShare(const Param& param,
+        double abs_error_used, const CombinedNode& r_node) {
+      return (param.eps - abs_error_used) * r_node.count() / param.n_points;
     }
   };
-
-
-  struct AlphaStat {
-   public:
-    SpRange alpha;
-    
-    OT_DEF(AlphaStat) {
-      OT_MY_OBJECT(alpha);
-    }
-
-   public:
-    void Init(const Param& param) {
-      alpha.InitEmptySet();
-    }
-    void Accumulate(const Param& param, const AlphaPoint& point) {
-      alpha |= SpRange(point.info().max2, point.info().max1);
-    }
-    void Accumulate(const Param& param,
-        const AlphaStat& stat, const Bound& bound, index_t n) {
-      alpha |= stat.alpha;
-    }
-    void Postprocess(const Param& param, const Bound& bound, index_t n) {}
-  };
-
-  struct RhoStat {
-   public:
-    SpRange rho;
-    
-    OT_DEF(RhoStat) {
-      OT_MY_OBJECT(rho);
-    }
-    
-   public:
-    void Init(const Param& param) {
-      rhos.InitEmptySet();
-    }
-    void Accumulate(const Param& param, const RhoPoint& point) {
-      rhos |= SpRange(point.info().rho);
-    }
-    void Accumulate(const Param& param,
-        const RhoStat& stat, const Bound& bound, index_t n) {
-      rhos |= stat.rhos;
-    }
-    void Postprocess(const Param& param, const Bound& bound, index_t n) {}
-  };
-
-  typedef SpNode<Bound, AlphaStat> AlphaNode;
-  typedef SpNode<Bound, RhoStat> RhoNode;
 };
 
 class AffinityAlpha {
  public:
-  typedef AffinityCommon::AlphaPoint QPoint;
-  typedef AffinityCommon::RhoPoint RPoint;
+  typedef AffinityCommon::CombinedPoint QPoint;
+  typedef AffinityCommon::CombinedPoint RPoint;
 
+  typedef AffinityCommon::Alpha Alpha;
+  
   typedef AffinityCommon::Param Param;
 
-  typedef AffinityCommon::AlphaStat QStat;
-  typedef AffinityCommon::RhoStat RStat;
-
-  typedef AffinityCommon::AlphaNode QNode;
-  typedef AffinityCommon::RhoNode RNode;
+  typedef AffinityCommon::CombinedNode QNode;
+  typedef AffinityCommon::CombinedNode RNode;
 
   typedef BlankGlobalResult GlobalResult;
 
-  struct BlankPostponed QPostponed;
+  typedef BlankQPostponed QPostponed;
 
   struct Delta {
    public:
@@ -200,12 +218,12 @@ class AffinityAlpha {
 
    public:
     void Init(const Param& param) {
-      alpha.max1 = -DBL_MAX;
-      alpha.max2 = -DBL_MAX;
+      alpha.max1 = param.pref;
+      alpha.max2 = param.pref;
       alpha.max1_index = -1;
     }
     void Postprocess(const Param& param,
-        const QPoint& q, const RNode& r_root) {}
+        const QPoint& q, index_t q_index, const RNode& r_root) {}
     void ApplyPostponed(const Param& param,
         const QPostponed& postponed, const QPoint& q) {}
   };
@@ -215,9 +233,7 @@ class AffinityAlpha {
     SpRange alpha;
 
     OT_DEF(QMassResult) {
-      OT_MY_OBJECT(alpha_hi);
-      OT_MY_OBJECT(alpha_lo);
-      OT_MY_OBJECT(xyz);
+      OT_MY_OBJECT(alpha);
     }
 
    public:
@@ -230,13 +246,13 @@ class AffinityAlpha {
     void ApplyDelta(const Param& param, const Delta& delta) {
       alpha.MaxWith(delta.alpha);
     }
-    bool ApplyPostponed(const Param& param,
+    void ApplyPostponed(const Param& param,
         const QPostponed& postponed, const QNode& q_node) {}
     void StartReaccumulate(const Param& param, const QNode& q_node) {
       alpha.InitEmptySet();
     }
     void Accumulate(const Param& param, const QResult& result) {
-      alpha |= SpRange(result.max2, result.max1);
+      alpha |= SpRange(result.alpha.max2, result.alpha.max1);
     }
     void Accumulate(const Param& param,
         const QMassResult& result, index_t n_points) {
@@ -245,6 +261,31 @@ class AffinityAlpha {
     void FinishReaccumulate(const Param& param, const QNode& q_node) {}
   };
 
+
+  /*
+  
+  As two-variable functions:
+  
+  \rho(i, k) = \sum_{j != i, j != k} max(0, S(j,k) - \alpha(j,k))
+  
+  \alpha(i, k) = min(0, \max_{j != k} S(j,j) + \alpha(j,j) + \rho(i,j))
+  
+  As one-variable rho:
+  
+  \rho(k) = \sum_{j != k} max(0, S(j,k) - \alpha(j,k))
+  
+  \sum max(0, S(j,k) - \alpha(j,k)) - max(0, S(k,k) - alpha(k,k))
+  
+  \alpha(i, k) = min(0, \max_{j != k} S(j,j) + \alpha(j,j)
+       + \rho(i) - max(0, S(i, j) - \alpha(i, j)))
+
+  \alpha(i) = min(0, \max^2{j}
+        S(j,j) + \alpha(j,j) + \rho(j) - max(0, S(i, j) - \alpha(i, j)))
+   except when i = j in which case we don't need to do the second part
+
+        S(j,j) + \alpha(j,j) + \rho(j) - S(i, j) + min(S(i,j), \alpha(i, j))
+  */
+
   struct PairVisitor {
    public:
     Alpha alpha;
@@ -253,32 +294,36 @@ class AffinityAlpha {
     void Init(const Param& param) {}
 
     bool StartVisitingQueryPoint(const Param& param,
-        const Point& q,
+        const QPoint& q,
         const RNode& r_node, const QMassResult& unapplied_mass_results,
         QResult* q_result, GlobalResult* global_result) {
       alpha = q_result->alpha;
+      return true;
     }
     void VisitPair(const Param& param,
-        const QPoint& q, index_t q_index,
-        const RPoint& r, index_t r_index) {
-      double cur_alpha;
-      
-      if (unlikely(q_index == r_index)) {
-        cur_alpha = r.info().rho + q.info().alpha.get(r_index);
-      } else {
-        double sim = param.Similarity(q.vec(), r.vec());
-        cur_alpha = min(
+        const QPoint& q, index_t q_index, const RPoint& r, index_t r_index) {
+      double candidate_alpha;
+
+      if (likely(q_index != r_index)) {
+        double sim = AffinityCommon::Helpers::Similarity(q.vec(), r.vec());
+        candidate_alpha = min(
             min(sim, q.info().alpha.get(r_index)) + r.info().rho,
             sim);
+        //double availability = 
+        //    r.info().rho - math::ClampNonNegative(
+        //        sim - q.info().alpha.get(r_index));
+        //candidate_alpha = math::ClampNonPositive(availability) + sim;
+      } else {
+        candidate_alpha = r.info().rho + q.info().alpha.get(r_index);
       }
-      
-      if (unlikely(cur_alpha > alpha.max2)) {
-        if (unlikely(cur_alpha > alpha.max1)) {
+
+      if (unlikely(candidate_alpha > alpha.max2)) {
+        if (unlikely(candidate_alpha > alpha.max1)) {
           alpha.max2 = alpha.max1;
-          alpha.max1 = cur_alpha;
-          alpha.max_index = r_index;
+          alpha.max1 = candidate_alpha;
+          alpha.max1_index = r_index;
         } else {
-          alpha.max2 = cur_alpha;
+          alpha.max2 = candidate_alpha;
         }
       }
     }
@@ -296,14 +341,23 @@ class AffinityAlpha {
         const QNode& q_node, const RNode& r_node,
         Delta* delta,
         GlobalResult* global_result, QPostponed* q_postponed) {
-      double sim_hi = param.SimilarityHi(q_node, r_node);
-      double sim_lo = param.SimilarityLo(q_node, r_node);
+      //WALDO
+      double sim_lo = AffinityCommon::Helpers::SimilarityLo(
+          param, q_node, r_node);
 
       delta->alpha.lo = min(
           min(q_node.stat().alpha.lo, sim_lo) + r_node.stat().rho.lo,
           sim_lo);
-      delta->alpha.hi =
-          q_node.stat().alpha.hi + r_node.stat().rho.hi;
+      if (q_node.begin() < r_node.end() && r_node.begin() < q_node.end()) {
+        delta->alpha.hi =
+            q_node.stat().alpha.hi + r_node.stat().rho.hi;
+      } else {
+        double sim_hi = AffinityCommon::Helpers::SimilarityHi(
+            param, q_node, r_node);
+        delta->alpha.hi = min(
+            min(q_node.stat().alpha.hi, sim_hi) + r_node.stat().rho.hi,
+            sim_hi);
+      }
 
       return true;
     }
@@ -325,7 +379,7 @@ class AffinityAlpha {
     }
     static double Heuristic(const Param& param,
         const QNode& q_node, const RNode& r_node, const Delta& delta) {
-      return -delta.sim_hi;
+      return -delta.alpha.hi;
     }
   };
 };
@@ -334,16 +388,13 @@ class AffinityRho {
  public:
   typedef AffinityCommon::Alpha Alpha;
 
-  typedef AffinityCommon::RhoPoint QPoint;
-  typedef AffinityCommon::AlphaPoint RPoint;
+  typedef AffinityCommon::CombinedPoint QPoint;
+  typedef AffinityCommon::CombinedPoint RPoint;
 
   typedef AffinityCommon::Param Param;
 
-  typedef AffinityCommon::RhoStat QStat;
-  typedef AffinityCommon::AlphaStat RStat;
-
-  typedef AffinityCommon::AlphaNode QNode;
-  typedef AffinityCommon::RhoNode RNode;
+  typedef AffinityCommon::CombinedNode QNode;
+  typedef AffinityCommon::CombinedNode RNode;
 
   typedef BlankGlobalResult GlobalResult;
 
@@ -359,7 +410,7 @@ class AffinityRho {
     
    public:
     void Init(const Param& param) {
-      Reset();
+      Reset(param);
     }
 
     void Reset(const Param& param) {
@@ -402,7 +453,15 @@ class AffinityRho {
       abs_error_used = 0;
     }
     void Postprocess(const Param& param,
-        const QPoint& q, const RNode& r_root) {}
+        const QPoint& q, index_t q_index, const RNode& r_root) {
+      double responsibility =
+          param.pref - q.info().alpha.get(q_index);
+
+      // Make sure we count ourselves regardless of sign.
+      if (responsibility < 0) {
+        rho += responsibility;
+      }
+    }
     void ApplyPostponed(const Param& param,
         const QPostponed& postponed, const QPoint& q) {
       rho += postponed.d_rho;
@@ -432,7 +491,7 @@ class AffinityRho {
     void ApplyDelta(const Param& param, const Delta& delta) {
       rho += delta.d_rho;
     }
-    bool ApplyPostponed(const Param& param,
+    void ApplyPostponed(const Param& param,
         const QPostponed& postponed, const QNode& q_node) {
       rho += postponed.d_rho;
       abs_error_used += postponed.abs_error_used;
@@ -464,18 +523,19 @@ class AffinityRho {
         const RNode& r_node, const QMassResult& unapplied_mass_results,
         QResult* q_result, GlobalResult* global_result) {
       rho = q_result->rho;
+      return true;
     }
     void VisitPair(const Param& param,
         const QPoint& q, index_t q_index,
         const RPoint& r, index_t r_index) {
-      double sim = param.Similarity(q.vec(), r.vec())
+      double responsibility =
+          AffinityCommon::Helpers::Similarity(
+              param, q.vec(), q_index, r.vec(), r_index)
           - r.info().alpha.get(r_index);
 
-      if (sim < 0 && likely(q_index != r_index)) {
-        sim = 0;
+      if (responsibility > 0) {
+        rho += responsibility;
       }
-      
-      rho += sim;
     }
     void FinishVisitingQueryPoint(const Param& param,
         const QPoint& q,
@@ -491,28 +551,42 @@ class AffinityRho {
         const QNode& q_node, const RNode& r_node,
         Delta* delta,
         GlobalResult* global_result, QPostponed* q_postponed) {
-      double sim_hi = param.SimilarityHi(q_node, r_node);
-      double sim_lo = param.SimilarityLo(q_node, r_node);
+      double sim_hi = AffinityCommon::Helpers::SimilarityHi(
+          param, q_node, r_node);
+      double sim_lo = AffinityCommon::Helpers::SimilarityLo(
+          param, q_node, r_node);
 
-      delta->rho.lo = (sim_lo - r_node.stat().alpha.hi) * r_node.count();
-      delta->rho.hi = max(0, sim_hi - r_node.stat().alpha.lo) * r_node.count();
+      // fprintf(stderr, "(%d,%d) alpha.lo,hi = (%f, %f)\n",
+      //           r_node.begin(), r_node.end(),
+      //           r_node.stat().alpha.hi,
+      //           r_node.stat().alpha.lo
+      //           );
+      delta->d_rho.lo = max(0.0, sim_lo - r_node.stat().alpha.hi)
+          * r_node.count();
+      delta->d_rho.hi = max(0.0, sim_hi - r_node.stat().alpha.lo)
+          * r_node.count();
 
-      return true;
+      return delta->d_rho.hi != 0;
     }
     static bool ConsiderPairExtrinsic(const Param& param,
         const QNode& q_node, const RNode& r_node, const Delta& delta,
         const QMassResult& q_mass_result, const GlobalResult& global_result,
         QPostponed* q_postponed) {
-      double abs_error = delta->d_rho.width() / 2;
+      /*
+      double abs_error = delta.d_rho.width() / 2;
       double rel_error_hi = abs_error / q_mass_result.rho.lo;
       
-      if (rel_error_hi < param.ErrorShare(q_mass_result.abs_error_used, r_node)) {
+      if (rel_error_hi < AffinityCommon::Helpers::ErrorShare(
+          param, q_mass_result.abs_error_used, r_node)) {
         q_postponed->abs_error_used += abs_error;
-        q_postponed->d_rho += delta->d_rho.mid();
+        q_postponed->d_rho += delta.d_rho.mid();
         return false;
       } else {
         return true;
       }
+      */
+      
+      return true;
     }
     static bool ConsiderQueryTermination(const Param& param,
         const QNode& q_node,
@@ -523,35 +597,163 @@ class AffinityRho {
     static double Heuristic(const Param& param,
         const QNode& q_node, const RNode& r_node, const Delta& delta) {
       // favor whatever brings our lower bound up the fastest
-      return -delta.d_rho_lo;
+      return -delta.d_rho.lo;
     }
   };
 };
 
+void FindExemplars(index_t dimensionality, index_t n_points,
+    CacheArray<AffinityAlpha::QPoint> *data_points) {
+  ArrayList<Vector> exemplars;
+  CacheReadIterator<AffinityAlpha::QPoint> point(data_points, 0);
+
+  exemplars.Init();
+
+  for (index_t point_i = 0; point_i < n_points; point_i++, point.Next()) {
+    //ot::Print(point->info());
+    if (point->info().rho > 0) {
+      exemplars.AddBack()->Copy(point->vec());
+    }
+  }
+  
+  ot::Print(exemplars);
+
+  Matrix m;
+  m.Init(dimensionality, exemplars.size());
+
+  for (index_t i = 0; i < exemplars.size(); i++) {
+    Vector dest;
+    m.MakeColumnVector(i, &dest);
+    dest.CopyValues(exemplars[i]);
+  }
+
+  data::Save("exemplars.txt", m);
+}
 
 void AffinityMain(datanode *module, const char *gnp_name) {
-  typename GNP::Param param;
+  AffinityAlpha::Param param;
 
   param.Init(fx_submodule(module, gnp_name, gnp_name));
 
-  TempCacheArray<typename GNP::QData> q_points;
-  TempCacheArray<typename GNP::QNode> q_nodes;
-  TempCacheArray<typename GNP::RData> r_points;
-  TempCacheArray<typename GNP::RNode> r_nodes;
-  TempCacheArray<typename GNP::QResult> q_results;
+  TempCacheArray<AffinityAlpha::QPoint> data_points;
+  TempCacheArray<AffinityAlpha::QNode> data_nodes;
 
-  nbr_utils::LoadKdTree(fx_submodule(module, "q", "q"),
-      &param, &q_points, &q_nodes);
-  nbr_utils::LoadKdTree(fx_submodule(module, "r", "r"),
-      &param, &r_points, &r_nodes);
+  index_t n_block_points = fx_param_int(
+      module, "n_block_points", 1024);
+  index_t n_block_nodes = fx_param_int(
+      module, "n_block_nodes", 128);
 
-  typename GNP::QResult default_result;
-  default_result.Init(param);
-  q_results.Init(default_result, q_points.end_index(),
-      q_points.n_block_elems());
+  datanode *data_module = fx_submodule(module, "data", "data");
 
-  Solver solver;
-  solver.InitSolve(fx_submodule(module, "solver", "solver"), param, 0,
-      q_points.cache(), q_nodes.cache(),
-      r_points.cache(), r_nodes.cache(), q_results.cache());
+  fx_timer_start(data_module, "read");
+
+  Matrix data_matrix;
+  MUST_PASS(data::Load(fx_param_str_req(data_module, ""), &data_matrix));
+  index_t n_points = data_matrix.n_cols();
+  index_t dimensionality = data_matrix.n_rows();
+  AffinityAlpha::QPoint default_point;
+  default_point.vec().Init(data_matrix.n_rows());
+  param.BootstrapMonochromatic(&default_point, data_matrix.n_cols());
+  data_points.Init(default_point, data_matrix.n_cols(), n_block_points);
+  for (index_t i = 0; i < data_matrix.n_cols(); i++) {
+    CacheWrite<AffinityAlpha::QPoint> point(&data_points, i);
+    point->vec().CopyValues(data_matrix.GetColumnPtr(i));
+  }
+  data_matrix.Destruct();
+  data_matrix.Init(0, 0);
+
+  fx_timer_stop(data_module, "read");
+
+  AffinityAlpha::QNode data_example_node;
+  data_example_node.Init(dimensionality, param);
+  data_nodes.Init(data_example_node, 0, n_block_nodes);
+  KdTreeMidpointBuilder
+      <AffinityAlpha::QPoint, AffinityAlpha::QNode, AffinityAlpha::Param>
+      ::Build(data_module, param, &data_points, &data_nodes);
+
+  // All the above is not any different for affinity than for anything
+  // else.  Now, time for the iteration.
+  int n_iter = fx_param_int(module, "n_iter", 1000);
+  int stable_iterations = 0;
+  index_t n_changed = n_points;
+
+  for (int iter = 0; iter < n_iter && stable_iterations < 50; iter++) {
+    {
+      TempCacheArray<AffinityAlpha::QResult> q_results_alpha;
+      
+      AffinityAlpha::QResult default_result_alpha;
+      default_result_alpha.Init(param);
+      q_results_alpha.Init(default_result_alpha, data_points.end_index(),
+          data_points.n_block_elems());
+
+      nbr_utils::ThreadedDualTreeSolver
+               < AffinityAlpha, DualTreeDepthFirst<AffinityAlpha> >::Solve(
+          fx_submodule(module, "threads", "iter%d_alpha", iter), param,
+          &data_points, &data_nodes, &data_points, &data_nodes,
+          &q_results_alpha);
+
+      for (index_t i = 0; i < n_points; i++) {
+        CacheWrite<AffinityAlpha::QPoint> point(&data_points, i);
+        CacheRead<AffinityAlpha::QResult> alpha(&q_results_alpha, i);
+        point->info().alpha = alpha->alpha;
+      }
+      nbr_utils::StatFixer<
+          AffinityAlpha::Param, AffinityAlpha::QPoint, AffinityAlpha::QNode>
+          ::Fix(param, &data_points, &data_nodes);
+    }
+
+    {
+      TempCacheArray<AffinityRho::QResult> q_results_rho;
+
+      AffinityRho::QResult default_result_rho;
+      default_result_rho.Init(param);
+      q_results_rho.Init(default_result_rho, data_points.end_index(),
+          data_points.n_block_elems());
+
+      nbr_utils::ThreadedDualTreeSolver
+               < AffinityRho, DualTreeDepthFirst<AffinityRho> >::Solve(
+          fx_submodule(module, "threads", "iter%d_rho", iter), param,
+          &data_points, &data_nodes, &data_points, &data_nodes,
+          &q_results_rho);
+
+      double temperature = 0;//1.0 * n_changed / n_points;
+      double lambda = param.lambda * (1.0 - temperature) + 0.6 * temperature;
+      double nonlambda = 1.0 - lambda;
+      n_changed = 0;
+
+      for (index_t i = 0; i < n_points; i++) {
+        CacheWrite<AffinityAlpha::QPoint> point(&data_points, i);
+        CacheRead<AffinityRho::QResult> compute_rho(&q_results_rho, i);
+        double old_rho = point->info().rho;
+        
+        point->info().rho = old_rho*lambda + compute_rho->rho*nonlambda;
+        
+        if ((old_rho > 0) != (point->info().rho > 0)) {
+          n_changed++;
+        }
+      }
+      
+      if (n_changed == 0) {
+        stable_iterations++;
+      } else {
+        stable_iterations = 0;
+      }
+      
+      fprintf(stderr, "------------- Changed = %"LI"d\n", n_changed);
+      
+      nbr_utils::StatFixer<
+          AffinityAlpha::Param, AffinityAlpha::QPoint, AffinityAlpha::QNode>
+          ::Fix(param, &data_points, &data_nodes);
+    }
+  }
+
+  FindExemplars(dimensionality, n_points, &data_points);
+}
+
+int main(int argc, char *argv[]) {
+  fx_init(argc, argv);
+
+  AffinityMain(fx_root, "affinity");
+
+  fx_done();
 }
diff --git a/fastlib/u/garryb/nbr/allnn.cc b/fastlib/u/garryb/nbr/allnn.cc
index e5dade8133..e69d5cfe0a 100644
--- a/fastlib/u/garryb/nbr/allnn.cc
+++ b/fastlib/u/garryb/nbr/allnn.cc
@@ -84,7 +84,7 @@ class Allnn {
     }
 
     void Postprocess(const Param& param,
-        const QPoint& q_point,
+        const QPoint& q_point, index_t q_index,
         const RNode& r_root) {}
     void ApplyPostponed(const Param& param,
         const QPostponed& postponed,
diff --git a/fastlib/u/garryb/nbr/build.py b/fastlib/u/garryb/nbr/build.py
index d510fd6ef2..73c8332b44 100644
--- a/fastlib/u/garryb/nbr/build.py
+++ b/fastlib/u/garryb/nbr/build.py
@@ -32,6 +32,10 @@ binrule(name = "allnn",
    sources = ["allnn.cc"],
    deplibs = [":nbr"])
 
+binrule(name = "affinity",
+   sources = ["affinity.cc"],
+   deplibs = [":nbr"])
+
 
 
 binrule(name = "cache_test",
diff --git a/fastlib/u/garryb/nbr/cache.h b/fastlib/u/garryb/nbr/cache.h
index 8c106d0af7..49193575f6 100644
--- a/fastlib/u/garryb/nbr/cache.h
+++ b/fastlib/u/garryb/nbr/cache.h
@@ -352,13 +352,13 @@ class CacheArray {
   }
   
   void ReleaseBlock_(BlockDevice::blockid_t blockid) {
-    DEBUG_ONLY(--metadatas_[blockid].lock_count);
+    DEBUG_ONLY(--metadatas_[blockid - begin_block_fake_].lock_count);
   }
   
   void ReleaseElement_(index_t element_id) {
     DEBUG_ONLY(BoundsCheck_(element_id));
     DEBUG_ONLY(
-        ReleaseBlock_((element_id >> n_block_elems_log_) - begin_block_fake_));
+        ReleaseBlock_((element_id >> n_block_elems_log_)));
   }
 };
 
diff --git a/fastlib/u/garryb/nbr/dfs.h b/fastlib/u/garryb/nbr/dfs.h
index f0c2a49fb1..500f237a7b 100644
--- a/fastlib/u/garryb/nbr/dfs.h
+++ b/fastlib/u/garryb/nbr/dfs.h
@@ -185,7 +185,7 @@ void DualTreeDepthFirst<GNP>::PushDown_(
       const typename GNP::QPoint *q_point = q_points_.StartRead(q_i);
       
       q_result->ApplyPostponed(param_, q_node_mut->postponed, *q_point);
-      q_result->Postprocess(param_, *q_point, *r_root_);
+      q_result->Postprocess(param_, *q_point, q_i, *r_root_);
       q_results_.StopWrite(q_i);
       q_points_.StopRead(q_i);
     }
@@ -356,13 +356,12 @@ void DualTreeDepthFirst<GNP>::BaseCase_(
 
     q_result->ApplyPostponed(param_, q_node_mut->postponed, *q_point);
 
-    if (unlikely(
-        visitor.StartVisitingQueryPoint(param_, *q_point, *r_node,
-          exclusive_unvisited, q_result, &global_result_))) {
+    if (visitor.StartVisitingQueryPoint(param_, *q_point, *r_node,
+          exclusive_unvisited, q_result, &global_result_)) {
       CacheReadIterator<typename GNP::RPoint> r_iter(&r_points_, r_node->begin());
-      index_t r_end = r_node->end();
 
-      for (index_t r_i = r_node->begin(); r_i < r_end; ++r_i, r_iter.Next()) {
+      for (index_t r_i = r_node->begin(); r_i < r_node->end(); ++r_i,
+          r_iter.Next()) {
         const typename GNP::RPoint *r_point = r_iter;
         
         visitor.VisitPair(param_, *q_point, q_i, *r_point, r_i);
diff --git a/fastlib/u/garryb/nbr/gnp.h b/fastlib/u/garryb/nbr/gnp.h
index 2a5482004e..7d6606db48 100644
--- a/fastlib/u/garryb/nbr/gnp.h
+++ b/fastlib/u/garryb/nbr/gnp.h
@@ -47,15 +47,15 @@ struct BlankQResult {
  public:
   OT_DEF(BlankQResult) {}
  public:
-  template<typename Param, typename Point, typename QPointInfo, typename RNode>
+  template<typename Param, typename Point, typename RNode>
   void Init(const Param& param,
-      const Point& q_point, const QPointInfo& q_info,
+      const Point& q_point,
       const RNode& r_root) {}
-  template<typename Param, typename Point, typename QPointInfo, typename RNode>
+  template<typename Param, typename Point, typename RNode>
   void Postprocess(const Param& param,
-      const Point& q_point, const QPointInfo& q_info,
+      const Point& q_point,
       const RNode& r_root) {}
-  template<typename Param, typename QPostponed, typename QPointInfo, typename Point>
+  template<typename Param, typename QPostponed, typename Point>
   void ApplyPostponed(const Param& param,
       const QPostponed& postponed,
       const Point& q_point) {}
diff --git a/fastlib/u/garryb/nbr/nbr_utils.h b/fastlib/u/garryb/nbr/nbr_utils.h
index 287c59c02a..ef731eeb85 100644
--- a/fastlib/u/garryb/nbr/nbr_utils.h
+++ b/fastlib/u/garryb/nbr/nbr_utils.h
@@ -10,6 +10,68 @@
 
 namespace nbr_utils {
 
+template<typename Param, typename Point, typename Node>
+class StatFixer {
+ public:
+  static void Fix(const Param &param,
+      CacheArray<Point> *points, CacheArray<Node> *nodes) {
+    StatFixer fixer;
+    fixer.InitFix(&param, points, nodes);
+  }
+  
+ private:
+  const Param *param_;
+  CacheArray<Point> points_;
+  CacheArray<Node> nodes_;
+  
+ public:
+  void InitFix(const Param *param,
+      CacheArray<Point> *points, CacheArray<Node> *nodes);
+  
+ private:
+  void FixRecursively_(index_t node_index);
+};
+
+template<typename Param, typename Point, typename Node>
+void StatFixer<Param, Point, Node>::InitFix(
+    const Param *param, CacheArray<Point> *points, CacheArray<Node> *nodes) {
+  param_ = param;
+  points_.Init(points, BlockDevice::READ);
+  nodes_.Init(nodes, BlockDevice::MODIFY);
+  FixRecursively_(0);
+  nodes_.Flush();
+  points_.Flush();
+}
+
+template<typename Param, typename Point, typename Node>
+void StatFixer<Param, Point, Node>::FixRecursively_(index_t node_index) {
+  CacheWrite<Node> node(&nodes_, node_index);
+
+  node->stat().Reset(*param_);
+
+  if (!node->is_leaf()) {
+    for (index_t k = 0; k < 2; k++) {
+      index_t child_index = node->child(k);
+
+      FixRecursively_(child_index);
+
+      CacheRead<Node> child(&nodes_, child_index);
+      node->stat().Accumulate(*param_, child->stat(),
+          child->bound(), child->count());
+    }
+    node->stat().Postprocess(*param_, node->bound(),
+        node->count());
+  } else {
+    CacheReadIterator<Point> point(&points_, node->begin());
+    
+    for (index_t i = 0; i < node->count(); i++, point.Next()) {
+      node->stat().Accumulate(*param_, *point);
+    }
+  }
+
+  node->stat().Postprocess(*param_, node->bound(), node->count());
+}
+
 template<typename GNP, typename Solver>
 class ThreadedDualTreeSolver {
  private:
@@ -42,7 +104,7 @@ class ThreadedDualTreeSolver {
           name.InitSprintf("grain_%d", work[i]);
           base_->mutex_.Lock();
           struct datanode *submodule = fx_submodule(base_->module_,
-              name.c_str(), name.c_str());
+              name.c_str(), "solver");
           base_->mutex_.Unlock();
 
           solver.InitSolve(submodule, *base_->param_, q_root_index,
@@ -213,7 +275,7 @@ struct MpiDualTreeConfig {
 
   void Copy(const MpiDualTreeConfig& other) {
     *this = other;
-KK  }
+  }
 
   OT_DEF(MpiDualTreeConfig) {
     OT_MY_OBJECT(n_threads);
diff --git a/fastlib/u/garryb/nbr/tkde.cc b/fastlib/u/garryb/nbr/tkde.cc
index 08ff5d3f12..9c2648fff3 100644
--- a/fastlib/u/garryb/nbr/tkde.cc
+++ b/fastlib/u/garryb/nbr/tkde.cc
@@ -41,18 +41,14 @@ class Tkde {
  public:
   /** The bounding type. Required by NBR. */
   typedef SpHrectBound<2> Bound;
-  /** The type of point in use. Required by NBR. */
-  typedef Vector Point;
+
+  typedef SpVectorPoint QPoint;
+  typedef SpVectorPoint RPoint;
 
   /** The type of kernel in use.  NOT required by NBR. */
   typedef EpanKernel Kernel;
 
-  /** Per-query statistic.  Required by NBR. */
-  typedef BlankStat QStat;
   typedef BlankGlobalResult GlobalResult;
-  
-  typedef BlankPointInfo QPointInfo;
-  typedef BlankPointInfo RPointInfo;
 
   /**
    * All parameters required by the execution of the algorithm.
@@ -62,21 +58,25 @@ class Tkde {
   struct Param {
    public:
     /**
-     * The threshold in use.
-     * This is a range to allow for epsilon checking.
+     * The threshold in use, with upper and lower bounds to prevent
+     * roundoff error.
+     *
+     * This is also normalized for dimensionality and the number of reference
+     * points.
      */
     SpRange thresh;
     /** The kernel in use. */
     Kernel kernel;
     /** The dimensionality of the data sets. */
     index_t dim;
-    /** The original threshold */
-    double threshold_orig;
+    /** The original threshold, before normalization. */
+    double nominal_threshold;
 
     OT_DEF(Param) {
-      OT_MY_OBJECT(kernel);
       OT_MY_OBJECT(thresh);
+      OT_MY_OBJECT(kernel);
       OT_MY_OBJECT(dim);
+      OT_MY_OBJECT(nominal_threshold);
     }
 
    public:
@@ -90,20 +90,16 @@ class Tkde {
      * Initialize parameters from a data node (Req NBR).
      */
     void Init(datanode *module) {
-      dim = -1;
       kernel.Init(fx_param_double_req(module, "h"));
-      threshold_orig = fx_param_double_req(module, "threshold");
+      nominal_threshold = fx_param_double_req(module, "threshold");
     }
 
-    void AnalyzePoint(const Point& q_point) {
-      if (dim == -1) {
-        dim = q_point.length();
-        double t = threshold_orig * kernel.CalcNormConstant(dim);
-        thresh.lo = t * (1.0 - 1.0e-4);
-        thresh.hi = t * (1.0 + 1.0e-4);
-      } else {
-        DEBUG_ASSERT_MSG(dim == q_point.length(), "Differing dimensionality");
-      }
+    void BootstrapMonochromatic(QPoint* point, index_t count) {
+      dim = point->vec().length();
+      double normalized_threshold =
+          nominal_threshold * kernel.CalcNormConstant(dim) * count;
+      thresh.lo = normalized_threshold * (1.0 - 1.0e5);
+      thresh.hi = normalized_threshold * (1.0 + 1.0e5);
     }
 
    public:
@@ -114,11 +110,11 @@ class Tkde {
      * the actual query point (not NBR).
      */
     double ComputeKernelSum(
-        const Vector& q_point,
+        const Vector& q,
         index_t r_count, const Vector& r_mass, double r_sumsq) const {
       double quadratic_term =
-          + r_count * la::Dot(q_point, q_point)
-          - 2.0 * la::Dot(q_point, r_mass)
+          + r_count * la::Dot(q, q)
+          - 2.0 * la::Dot(q, r_mass)
           + r_sumsq;
       return r_count - quadratic_term * kernel.inv_bandwidth_sq();
     }
@@ -245,8 +241,8 @@ class Tkde {
     /**
      * Accumulate data from a single point (Req NBR).
      */
-    void Accumulate(const Param& param, const Vector& point) {
-      moment_info.Add(1, point, la::Dot(point, point));
+    void Accumulate(const Param& param, const QPoint& point) {
+      moment_info.Add(1, point.vec(), la::Dot(point.vec(), point.vec()));
     }
 
     /**
@@ -264,8 +260,18 @@ class Tkde {
     void Postprocess(const Param& param, const Bound& bound, index_t n) {
     }
   };
-  
-  
+
+  /**
+   * Query-tree statistic.
+   *
+   * Note that this statistic is not actually needed and a blank statistic
+   * is fine, but QStat must equal RStat in order for us to allow
+   * monochromatic execution.
+   *
+   * This limitation may be removed in a further version of NBR.
+   */
+  typedef RStat QStat;
+ 
   /**
    * Query node.
    */
@@ -350,7 +356,7 @@ class Tkde {
     }
 
     void Postprocess(const Param& param,
-        const Vector& q_point, const QPointInfo& q_info,
+        const QPoint& q,
         const RNode& r_root) {
       if (density > param.thresh.hi) {
         label &= LAB_HI;
@@ -362,12 +368,13 @@ class Tkde {
 
     void ApplyPostponed(const Param& param,
         const QPostponed& postponed,
-        const Vector& q_point) {
+        const QPoint& q) {
       label &= postponed.label; /* bitwise OR */
       DEBUG_ASSERT(label != LAB_NEITHER);
 
       if (!postponed.moment_info.is_empty()) {
-        density += postponed.moment_info.ComputeKernelSum(param, q_point);
+        density += postponed.moment_info.ComputeKernelSum(
+            param, q.vec());
       }
     }
   };
@@ -463,8 +470,7 @@ class Tkde {
     // - this function must assume that global_result is incomplete (which is
     // reasonable in allnn)
     bool StartVisitingQueryPoint(const Param& param,
-        const Vector& q_point,
-        const QPointInfo& q_info,
+        const QPoint& q,
         const RNode& r_node,
         const QMassResult& unapplied_mass_results,
         QResult* q_result,
@@ -473,17 +479,17 @@ class Tkde {
         return false;
       }
 
-      double distance_sq_lo = r_node.bound().MinDistanceSqToPoint(q_point);
+      double distance_sq_lo = r_node.bound().MinDistanceSqToPoint(q.vec());
 
       if (unlikely(distance_sq_lo > param.kernel.bandwidth_sq())) {
         return false;
       }
 
-      double distance_sq_hi = r_node.bound().MaxDistanceSqToPoint(q_point);
+      double distance_sq_hi = r_node.bound().MaxDistanceSqToPoint(q.vec());
 
       if (unlikely(distance_sq_hi < param.kernel.bandwidth_sq())) {
         q_result->density += r_node.stat().moment_info.ComputeKernelSum(
-            param, q_point);
+            param, q.vec());
         return false;
       }
 
@@ -493,15 +499,14 @@ class Tkde {
     }
 
     void VisitPair(const Param& param,
-        const Vector& q_point, const QPointInfo& q_info, index_t q_index,
-        const Vector& r_point, const RPointInfo& r_info, index_t r_index) {
-      double distance = la::DistanceSqEuclidean(q_point, r_point);
+        const QPoint& q, index_t q_index,
+        const RPoint& r, index_t r_index) {
+      double distance = la::DistanceSqEuclidean(q.vec(), r.vec());
       density += param.kernel.EvalUnnormOnSq(distance);
     }
 
     void FinishVisitingQueryPoint(const Param& param,
-        const Vector& q_point,
-        const QPointInfo& q_info,
+        const QPoint& q,
         const RNode& r_node,
         const QMassResult& unapplied_mass_results,
         QResult* q_result,
@@ -627,7 +632,7 @@ int main(int argc, char *argv[]) {
       fx_root, "tkde");
   MPI_Finalize();
 #else
-  nbr_utils::ThreadedDualTreeMain<Tkde, DualTreeDepthFirst<Tkde> >(
+  nbr_utils::MonochromaticDualTreeMain<Tkde, DualTreeDepthFirst<Tkde> >(
       fx_root, "tkde");
 #endif