From 5909ee739b56b905970b0a988a0f3fe6a2908b82 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 26 Nov 2019 13:52:48 -0500
Subject: overhaul now working, adding measurements

---
 space_wolff.hpp | 701 +++++++++++++++++++++++++++-----------------------------
 spheres.cpp     |  50 ++--
 2 files changed, 359 insertions(+), 392 deletions(-)

diff --git a/space_wolff.hpp b/space_wolff.hpp
index c87697e..cc4cd4f 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -1,45 +1,42 @@
 
-#include <vector>
-#include <list>
-#include <set>
-#include <iostream>
+#include "randutils/randutils.hpp"
+#include <cmath>
+#include <eigen3/Eigen/Dense>
 #include <fstream>
 #include <functional>
-#include <random>
+#include <iostream>
+#include <list>
 #include <queue>
-#include <cmath>
-#include <eigen3/Eigen/Dense>
-#include "randutils/randutils.hpp"
+#include <random>
+#include <set>
+#include <vector>
 
-const std::array<std::array<unsigned, 16>, 16> smiley = 
-  {{
-    {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}},
-    {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
-    {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}},
-    {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}},
-    {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}},
-    {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}},
-    {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}},
-    {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
-    {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
-    {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}},
-    {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}},
-    {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}},
-    {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}},
-    {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}},
-    {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
-    {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}
-  }};
-
-template <class T, int D>
-using vector = Eigen::Matrix<T, D, 1>;
-
-template <class T, int D>
-using matrix = Eigen::Matrix<T, D, D>;
-
-template <class T, int D>
-vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) {
-  vector<T, D> v;
+const std::array<std::array<unsigned, 16>, 16> smiley = {
+    {{{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}},
+     {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
+     {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}},
+     {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}},
+     {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}},
+     {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}},
+     {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}},
+     {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
+     {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
+     {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}},
+     {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}},
+     {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}},
+     {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}},
+     {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}},
+     {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
+     {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}}};
+
+template <class T, int D> using Vector = Eigen::Matrix<T, D, 1>;
+
+template <class T, int D> using Matrix = Eigen::Matrix<T, D, D>;
+
+/** brief diff - periodic subtraction of vectors
+ */
+template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) {
+  Vector<T, D> v;
 
   for (unsigned i = 0; i < D; i++) {
     v(i) = std::abs(v1(i) - v2(i));
@@ -51,425 +48,411 @@ vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) {
   return v;
 }
 
-template <class T, int D, class state>
-class spin {
-  public:
-    vector<T, D> x;
-    state s;
+template <class T, int D, class S> class Spin {
+public:
+  Vector<T, D> x;
+  S s;
 };
 
-template <class T, int D>
-class euclidean {
-  private:
-    T L;
-
-  public:
-    vector<T, D> t;
-    matrix<T, D> r;
-    euclidean(T L) : L(L) {
-      for (unsigned i = 0; i < D; i++) {
-        t(i) = 0;
-        r(i, i) = 1;
-        for (unsigned j = 1; j < D; j++) {
-          r(i, (i + j) % D) = 0;
-        }
+template <class T, int D> class Euclidean {
+private:
+  T L;
+
+public:
+  Vector<T, D> t;
+  Matrix<T, D> r;
+
+  /** brief Euclidean - default constructor, constructs the identity
+   */
+  Euclidean(T L) : L(L) {
+    for (unsigned i = 0; i < D; i++) {
+      t(i) = 0;
+      r(i, i) = 1;
+      for (unsigned j = 1; j < D; j++) {
+        r(i, (i + j) % D) = 0;
       }
     }
+  }
 
-    euclidean(T L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
-      t = t0;
-      r = r0;
-    }
+  Euclidean(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
+    t = t0;
+    r = r0;
+  }
 
-    template <class state>
-    spin<T, D, state> act(const spin<T, D, state>& s) const {
-      spin<T, D, state> s_new;
+  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+    Spin<T, D, S> s_new;
 
-      s_new.x = t + r * s.x;
-      s_new.s = s.s;
+    s_new.x = t + r * s.x;
+    s_new.s = s.s;
 
-      for (unsigned i = 0; i < D; i++) {
-        s_new.x(i) = fmod(L + s_new.x(i), L);
-      }
-
-      return s_new;
+    for (unsigned i = 0; i < D; i++) {
+      s_new.x(i) = fmod(L + s_new.x(i), L);
     }
 
-    euclidean act(const euclidean& x) const {
-      vector<T, D> tnew = r * x.t + t;
-      matrix<T, D> rnew = r * x.r;
-
-      for (unsigned i = 0; i < D; i++) {
-        tnew(i) = fmod(L + tnew(i), L);
-      }
+    return s_new;
+  }
 
-      euclidean pnew(this->L, tnew, rnew);
+  Euclidean act(const Euclidean& x) const {
+    Vector<T, D> tnew = r * x.t + t;
+    Matrix<T, D> rnew = r * x.r;
 
-      return pnew;
+    for (unsigned i = 0; i < D; i++) {
+      tnew(i) = fmod(L + tnew(i), L);
     }
 
-    euclidean inverse() const {
-      vector<T, D> tnew = - r.transpose() * t;
-      matrix<T, D> rnew = r.transpose();
+    Euclidean pnew(this->L, tnew, rnew);
+
+    return pnew;
+  }
 
-      euclidean pnew(this->L, tnew, rnew);
+  Euclidean inverse() const {
+    Vector<T, D> tnew = -r.transpose() * t;
+    Matrix<T, D> rnew = r.transpose();
 
-      return pnew;
-    }
-};
+    Euclidean pnew(this->L, tnew, rnew);
 
-template <int D>
-class dictionary {
-  private:
-    unsigned L;
-    std::vector<std::set<unsigned>> d;
+    return pnew;
+  }
+};
 
-  public:
-    dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {};
+template <class T, int D, class S> class Octree {
+private:
+  T L;
+  unsigned N;
+  std::vector<std::set<Spin<T, D, S>*>> data;
 
-    template <class T>
-    unsigned dictionary_index(vector<T, D> x) const {
-      unsigned pos_ind = 0;
+public:
+  Octree(T L, unsigned depth) : L(L), N(pow(2, depth)), data(pow(N, D)){};
 
-      for (unsigned i = 0; i < D; i++) {
-        pos_ind += pow(L, i) * (unsigned)x(i);
-      };
+  unsigned ind(Vector<T, D> x) const {
+    unsigned pos_ind = 0;
 
-      return pos_ind;
+    for (unsigned i = 0; i < D; i++) {
+      pos_ind += pow(N, i) * (unsigned)std::floor(N * x(i) / L);
     }
 
-    template <class T>
-    void record(vector<T, D> x, unsigned ind) {
-      d[dictionary_index<T>(x)].insert(ind);
-    };
+    assert(pos_ind < data.size());
 
-    template <class T>
-    void remove(vector<T, D> x, unsigned ind) {
-      d[dictionary_index<T>(x)].erase(ind);
-    };
+    return pos_ind;
+  }
 
-    template <class T>
-    std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const {
-      return nearest_neighbors_of(dictionary_index<T>(x), depth, {});
-    };
+  void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); };
 
-    std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, std::list<unsigned> ignores) const {
-      std::set<unsigned> ns = d[ind];
+  void remove( Spin<T, D, S>* s) { data[ind(s->x)].erase(s); };
 
-      if (depth > 0) {
-        for (unsigned i = 0; i < D; i++) {
-          if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k){return i==k;})) {
-            std::list<unsigned> ignores_new = ignores;
-            ignores_new.push_back(i);
-            for (signed j : {-1, 1}) {
-              unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));
+  std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x, unsigned depth) const {
+    std::set<Spin<T, D, S>*> ns;
+    std::set<unsigned> seen;
+    nearest_neighbors_of(ind(x), depth, ns, seen);
+    return ns;
+  };
 
-              std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new);
+  void nearest_neighbors_of(unsigned ind, unsigned depth, std::set<Spin<T, D, S>*>& ns,
+                            std::set<unsigned>& seen) const {
+    ns.insert(data[ind].begin(), data[ind].end());
+    seen.insert(ind);
 
-              for (unsigned guy : nns) {
-                ns.insert(guy);
-              }
-            }
+    if (depth > 0) {
+      for (unsigned i = 0; i < D; i++) {
+        for (signed j : {-1, 1}) {
+          unsigned nind = pow(N, i + 1) * (ind / ((unsigned)pow(N, i + 1))) +
+                          fmod(pow(N, i + 1) + ind + j * pow(N, i), pow(N, i + 1));
+
+          if (!seen.contains(nind)) {
+            seen.insert(nind);
+            nearest_neighbors_of(nind, depth - 1, ns, seen);
           }
         }
       }
-
-      return ns;
-    };
+    }
+  }
 };
 
 class quantity {
-  private:
-    double total;
-    double total2;
-    std::vector<double> C;
-    unsigned wait;
-
-  public:
-    unsigned n;
-    std::list<double> hist;
-
-    quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
-      n = 0;
-      total = 0;
-      total2 = 0;
-    }
+private:
+  double total;
+  double total2;
+  std::vector<double> C;
+  unsigned wait;
+
+public:
+  unsigned n;
+  std::list<double> hist;
+
+  quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
+    n = 0;
+    total = 0;
+    total2 = 0;
+  }
 
-    void add(double x) {
-      if (n > wait) {
-        hist.push_front(x);
-        if (hist.size() > C.size()) {
-          hist.pop_back();
-          unsigned t = 0;
-          for (double a : hist) {
-            C[t] += a * x;
-            t++;
-          }
-          total += x;
-          total2 += pow(x, 2);
+  void add(double x) {
+    if (n > wait) {
+      hist.push_front(x);
+      if (hist.size() > C.size()) {
+        hist.pop_back();
+        unsigned t = 0;
+        for (double a : hist) {
+          C[t] += a * x;
+          t++;
         }
+        total += x;
+        total2 += pow(x, 2);
       }
-      n++;
-    }
-
-    double avg() const {
-      return total / (n - wait);
     }
+    n++;
+  }
 
-    double avg2() const {
-      return total2 / (n - wait);
-    }
+  double avg() const { return total / (n - wait); }
 
-    std::vector<double> ρ() const {
-      double C0 = C.front() / (n - wait);
-      double avg2 = pow(total / (n - wait), 2);
+  double avg2() const { return total2 / (n - wait); }
 
-      std::vector<double> ρtmp;
+  std::vector<double> ρ() const {
+    double C0 = C.front() / (n - wait);
+    double avg2 = pow(total / (n - wait), 2);
 
-      for (double Ct : C) {
-        ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
-      }
+    std::vector<double> ρtmp;
 
-      return ρtmp;
+    for (double Ct : C) {
+      ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
     }
 
-    std::array<double, 2> τ() const {
-      double τtmp = 0.5;
-      unsigned M = 1;
-      double c = 8.0;
+    return ρtmp;
+  }
 
-      std::vector<double> ρ_tmp = this->ρ();
+  std::array<double, 2> τ() const {
+    double τtmp = 0.5;
+    unsigned M = 1;
+    double c = 8.0;
 
-      while (c * τtmp > M && M < C.size()) {
-        τtmp += ρ_tmp[M];
-        M++;
-      }
+    std::vector<double> ρ_tmp = this->ρ();
 
-      return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+    while (c * τtmp > M && M < C.size()) {
+      τtmp += ρ_tmp[M];
+      M++;
     }
 
-    double σ() const {
-      return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
-    }
+    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+  }
 
-    double serr() const {
-      return sqrt(this->σ());
-    }
+  double σ() const {
+    return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+  }
 
-    unsigned num_added() const {
-      return n - wait;
-    }
+  double serr() const { return sqrt(this->σ()); }
+
+  unsigned num_added() const { return n - wait; }
 };
 
-template <class U, int D, class state>
-class model {
-  public:
-    U L;
-    euclidean<U, D> s0;
-    std::vector<spin<U, D, state>> s;
-    dictionary<D> dict;
-    std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors;
-    std::function<double(spin<U, D, state>, spin<U, D, state>, spin<U, D, state>)> Z;
-    std::function<double(spin<U, D, state>)> B;
-    std::vector<matrix<U, D>> mats;
-    std::vector<vector<U, D>> steps;
-    long double E;
-    quantity Eq;
-    quantity Cq;
-
-    void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
-      if (level > 0) {
-        unsigned new_level = level - 1;
-        unsigned old_length = sequences.size();
-        for (std::array<double, D>& sequence : sequences) {
-          std::array<double, D> new_sequence = sequence;
-          new_sequence[new_level] = -1;
-          sequences.push_front(new_sequence);
-        }
-        one_sequences(sequences, new_level);
+template <class U, int D, class S> class Model;
+
+/*
+template <class U, int D, class S>
+  class measurement {
+    public:
+      virtual void pre_cluster(const Model<U, D, S>&, unsigned, const Euclidean<U, D>&) {};
+      virtual void plain_bond_visited(const Model<U, D, S>&, const X_t&, double) {};
+      virtual void plain_site_transformed(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const X_t&) {};
+
+      virtual void ghost_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const X_t&, const X_t&, double) {};
+      virtual void ghost_site_transformed(const system<R_t, X_t, G_t>&, const R_t&) {};
+
+      virtual void post_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&) {};
+  };
+  */
+
+template <class U, int D, class S> class Model {
+public:
+  U L;
+  Euclidean<U, D> s0;
+  std::vector<Spin<U, D, S>> s;
+  Octree<U, D, S> dict;
+  unsigned dDepth;
+  unsigned nDepth;
+  std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z;
+  std::function<double(const Spin<U, D, S>&)> B;
+  std::vector<Matrix<U, D>> mats;
+  std::vector<Vector<U, D>> steps;
+  double E;
+
+  void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+    if (level > 0) {
+      unsigned new_level = level - 1;
+      unsigned old_length = sequences.size();
+      for (std::array<double, D>& sequence : sequences) {
+        std::array<double, D> new_sequence = sequence;
+        new_sequence[new_level] = -1;
+        sequences.push_front(new_sequence);
       }
+      one_sequences(sequences, new_level);
     }
+  }
 
-    model(U L, std::function<double(spin<U, D, state>, spin<U, D, state>, spin<U, D, state>)> Z,
-        std::function<double(spin<U, D, state>)> B,
-        std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) : 
-      L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) {
-        std::array<double, D> ini_sequence;
-        ini_sequence.fill(1);
-        std::list<std::array<double, D>> sequences;
-        sequences.push_back(ini_sequence);
-
-        one_sequences(sequences, D);
-
-        sequences.pop_front(); // don't want the identity matrix!
-
-        for (std::array<double, D> sequence : sequences) {
-          matrix<U, D> m;
-          for (unsigned i = 0; i < D; i++) {
-            for (unsigned j = 0; j < D; j++) {
-              if (i == j) {
-                m(i, j) = sequence[i];
-              } else {
-                m(i, j) = 0;
-              }
-            }
-          }
+  Model(U L, std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z,
+        std::function<double(const Spin<U, D, S>&)> B, unsigned dDepth, unsigned nDepth)
+      : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), dDepth(dDepth), nDepth(nDepth) {
+    std::array<double, D> ini_sequence;
+    ini_sequence.fill(1);
+    std::list<std::array<double, D>> sequences;
+    sequences.push_back(ini_sequence);
 
-          mats.push_back(m);
+    one_sequences(sequences, D);
 
-          vector<U, D> v;
-          for (unsigned i = 0; i < D; i++) {
-            if (sequence[i] == 1) {
-              v(i) = 0;
-            } else {
-              v(i) = L / 2;
-            }
+    sequences.pop_front(); // don't want the identity matrix!
+
+    for (std::array<double, D> sequence : sequences) {
+      Matrix<U, D> m;
+      for (unsigned i = 0; i < D; i++) {
+        for (unsigned j = 0; j < D; j++) {
+          if (i == j) {
+            m(i, j) = sequence[i];
+          } else {
+            m(i, j) = 0;
           }
+        }
+      }
+
+      mats.push_back(m);
 
-          steps.push_back(v);
+      Vector<U, D> v;
+      for (unsigned i = 0; i < D; i++) {
+        if (sequence[i] == 1) {
+          v(i) = 0;
+        } else {
+          v(i) = L / 2;
         }
+      }
+
+      steps.push_back(v);
+    }
 
-        for (unsigned i = 0; i < D; i++) {
-          for (unsigned j = 0; j < D; j++) {
-            if (i != j) {
-              matrix<U, D> m;
-              for (unsigned k = 0; k < D; k++) {
-                for (unsigned l = 0; l < D; l++) {
-                  if ((k == i && l == j) || (k == j && l == i)) {
-                    if (i < j) {
-                      m(k, l) = 1;
-                    } else {
-                      m(k, l) = -1;
-                    }
-                  } else {
-                    m(k, l) = 0;
-                  }
+    for (unsigned i = 0; i < D; i++) {
+      for (unsigned j = 0; j < D; j++) {
+        if (i != j) {
+          Matrix<U, D> m;
+          for (unsigned k = 0; k < D; k++) {
+            for (unsigned l = 0; l < D; l++) {
+              if ((k == i && l == j) || (k == j && l == i)) {
+                if (i < j) {
+                  m(k, l) = 1;
+                } else {
+                  m(k, l) = -1;
                 }
+              } else {
+                m(k, l) = 0;
               }
-              mats.push_back(m);
             }
           }
+          mats.push_back(m);
         }
       }
-
-    void update_energy() {
-      E = 0;
-      for (unsigned i = 0; i < s.size(); i++) {
-        for (unsigned j = 0; j < i; j++) {
-        //  E -= Z(s[i], s[j]);
-        }
-
-        E -= B(s0.inverse().act(s[i]));
-      }
     }
+  }
 
-    void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) {
-      unsigned cluster_size = 0;
-      std::uniform_real_distribution<double> dist(0.0, 1.0);
+  void update_energy() {
+    E = 0;
+    for (unsigned i = 0; i < s.size(); i++) {
+      for (unsigned j = 0; j < i; j++) {
+          E -= Z(s[i], s[j]);
+      }
 
-      std::queue<unsigned> queue;
-      queue.push(ind);
+      E -= B(s0.inverse().act(s[i]));
+    }
+  }
 
-      std::vector<bool> visited(s.size() + 1, false);
+  void step(double T, unsigned ind, Euclidean<U, D> r, std::mt19937& rng) {
+    unsigned cluster_size = 0;
+    std::uniform_real_distribution<double> dist(0.0, 1.0);
 
-      while (!queue.empty()) {
-        unsigned i = queue.front();
-        queue.pop();
+    std::queue<Spin<U, D, S>*> queue;
+    queue.push(&(s[ind]));
 
-        if (!visited[i]) {
-          visited[i] = true;
+    std::set<Spin<U, D, S>*> visited;
 
-          bool we_are_ghost = i == s.size();
+    while (!queue.empty()) {
+      Spin<U, D, S>* si = queue.front();
+      queue.pop();
 
-          spin<U, D, state> si_new;
-          euclidean<U, D> s0_new(L);
+      if (!visited.contains(si)) {
+        visited.insert(si);
 
-          if (we_are_ghost) {
-            s0_new = r.act(s0);
-          } else {
-            si_new = r.act(s[i]);
+        if (si == NULL) {
+          Euclidean<U, D> s0_new = r.act(s0);
+          for (Spin<U, D, S>& ss : s) {
+            Spin<U, D, S> s0s_old = s0.inverse().act(ss);
+            Spin<U, D, S> s0s_new = s0_new.inverse().act(ss);
+            double p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
+            if (dist(rng) < p) {
+              queue.push(&ss);
+            }
+            s0 = s0_new;
           }
-
-          for (unsigned j : neighbors(*this, i, si_new)) {
-            if (j != i) {
+        } else {
+          Spin<U, D, S> si_new = r.act(*si);
+          std::set<Spin<U, D, S>*> all_neighbors;
+          std::set<Spin<U, D, S>*> current_neighbors = dict.neighbors(si->x, nDepth);
+          std::set<Spin<U, D, S>*> new_neighbors = dict.neighbors(si_new.x, nDepth);
+
+          all_neighbors.insert(current_neighbors.begin(), current_neighbors.end());
+          all_neighbors.insert(new_neighbors.begin(), new_neighbors.end());
+          all_neighbors.insert(NULL);
+          for (Spin<U, D, S>* sj : all_neighbors) {
+            if (sj != si) {
               double p;
-              bool neighbor_is_ghost = j == s.size();
-
-              if (we_are_ghost || neighbor_is_ghost) {
-                spin<U, D, state> s0s_old, s0s_new;
-                unsigned non_ghost;
-
-                if (neighbor_is_ghost) {
-                  non_ghost = i;
-                  s0s_old = s0.inverse().act(s[i]);
-                  s0s_new = s0.inverse().act(si_new);
-                } else {
-                  non_ghost = j;
-                  s0s_old = s0.inverse().act(s[j]);
-                  s0s_new = s0_new.inverse().act(s[j]);
-                }
-
+              if (sj == NULL) {
+                Spin<U, D, S> s0s_old = s0.inverse().act(*si);
+                Spin<U, D, S> s0s_new = s0.inverse().act(si_new);
                 p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
               } else {
-                p = Z(s[i], s[j], si_new);
+                p = 1.0 - exp(-(Z(si_new, *sj) - Z(*si, *sj)) / T);
               }
-
               if (dist(rng) < p) {
-                queue.push(j);
+                queue.push(sj);
               }
             }
           }
-
-          if (we_are_ghost) {
-            s0 = s0_new;
-          } else {
-            dict.remove(s[i].x, i);
-            s[i] = si_new;
-            dict.record(s[i].x, i);
-            cluster_size++;
-          }
+          dict.remove(si);
+          *si = si_new;
+          dict.insert(si);
+          cluster_size++;
         }
       }
-
-      Cq.add(cluster_size);
     }
+  }
 
-    void wolff(double T, unsigned N, std::mt19937& rng) {
-      std::uniform_real_distribution<double> t_dist(0, L);
-      std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
-      std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
-      std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1);
-
-      for (unsigned i = 0; i < N; i++) {
-        vector<U, D> t;
-        matrix<U, D> m;
-        unsigned flip = coin(rng);
-        if (flip < mats.size()) {
-          for (unsigned j = 0; j < D; j++) {
-            t(j) = (U)t_dist(rng);
-          }
-          m = mats[flip];
-        } else {
-          for (unsigned j = 0; j < D; j++) {
-            for (unsigned k = 0; k < D; k++) {
-              if (j == k) {
-                m(j, k) = 1;
-              } else {
-                m(j, k) = 0;
-              }
+  void wolff(double T, unsigned N, std::mt19937& rng) {
+    std::uniform_real_distribution<double> t_dist(0, L);
+    std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
+    std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
+    std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1);
+
+    for (unsigned i = 0; i < N; i++) {
+      Vector<U, D> t;
+      Matrix<U, D> m;
+      unsigned flip = coin(rng);
+      if (flip < mats.size()) {
+        for (unsigned j = 0; j < D; j++) {
+          t(j) = (U)t_dist(rng);
+        }
+        m = mats[flip];
+      } else {
+        for (unsigned j = 0; j < D; j++) {
+          for (unsigned k = 0; k < D; k++) {
+            if (j == k) {
+              m(j, k) = 1;
+            } else {
+              m(j, k) = 0;
             }
           }
-
-          t = steps[flip - mats.size()];
         }
 
-        euclidean<U, D> g(L, t, m);
+        t = steps[flip - mats.size()];
+      }
 
-        this->step(T, ind_dist(rng), g, rng);
+      Euclidean<U, D> g(L, t, m);
 
-        this->update_energy();
-        Eq.add(E);
-      }
+      this->step(T, ind_dist(rng), g, rng);
+
+      this->update_energy();
     }
+  }
 };
-
diff --git a/spheres.cpp b/spheres.cpp
index b0f4fb4..633a26a 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -34,65 +34,49 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  std::function<double(spin<double, D, double>, spin<double, D, double>, spin<double, D, double>)> Z =
-    [L] (spin<double, D, double> s1, spin<double, D, double> s2, spin<double, D, double> s1_new) -> double {
-      vector<double, D> diff_old = diff(L, s1.x, s2.x);
-      vector<double, D> diff_new = diff(L, s1_new.x, s2.x);
+  std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
+    [L] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
+      Vector<double, D> d = diff(L, s1.x, s2.x);
 
       double rad_sum = pow(s1.s + s2.s, 2);
 
-      bool old_overlap = diff_old.transpose() * diff_old < rad_sum;
-      bool new_overlap = diff_new.transpose() * diff_new < rad_sum;
+      bool overlap = d.transpose() * d < rad_sum;
 
-      if (new_overlap) {
-        return 1.0;
+      if (overlap) {
+        return -1e8;
       } else {
-        return 0.0;
+        return 0;
       }
     };
 
-  std::function<double(spin<double, D, double>)> B =
-    [L, H] (spin<double, D, double> s) -> double {
+  std::function<double(Spin<double, D, double>)> B =
+    [L, H] (Spin<double, D, double> s) -> double {
       return H * sin(2 * M_PI * 3 * s.x(0) / L);
     };
 
-  std::function<std::set<unsigned>(model<double, D, double>&, unsigned, spin<double, D, double>)> neighbors =
-    [] (model<double, D, double>& m, unsigned i0, spin<double, D, double> s1) -> std::set<unsigned> {
-      std::set<unsigned> nn;
-      if (i0 < m.s.size()) {
-        std::set<unsigned> n_old = m.dict.neighbors(m.s[i0].x, 1);
-        std::set<unsigned> n_new = m.dict.neighbors(s1.x, 1);;
-        nn.insert(n_old.begin(), n_old.end());
-        nn.insert(n_new.begin(), n_new.end());
-        nn.insert(m.s.size());
-      } else {
-        for (unsigned i = 0; i < m.s.size(); i++) {
-          nn.insert(i);
-        }
-      }
-      return nn;
-    };
-
-  model<double, D, double> sphere(L, Z, B, neighbors);
+  Model<double, D, double> sphere(L, Z, B, std::floor(log2(L)), 2);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
   std::uniform_real_distribution<double> dist(0.0, L);
 
+  sphere.s.reserve(n);
+
   for (unsigned i = 0; i < n; i++) {
-    vector<double, D> pos = {dist(rng), dist(rng)};
+    Vector<double, D> pos = {dist(rng), dist(rng)};
     sphere.s.push_back({pos, 0.5});
-    sphere.dict.record<double>(pos, i);
+    sphere.dict.insert(&sphere.s.back());
   }
 
+
   sphere.wolff(T, N, rng);
 
   std::ofstream snapfile;
   snapfile.open("sphere_snap.dat");
 
-  for (spin<double, D, double> s : sphere.s) {
-    spin<double, D, double> rs = sphere.s0.inverse().act(s);
+  for (Spin<double, D, double> s : sphere.s) {
+    Spin<double, D, double> rs = sphere.s0.inverse().act(s);
     snapfile << rs.x.transpose() << "\n";
   }
 
-- 
cgit v1.2.3-54-g00ecf