From 452600249b16eb1239931b49f3d0fa795a611f5a Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 14 Feb 2020 07:44:50 -0500
Subject: Broken: attempt to implement dynamic matching of pair-flipped
 spheres.

---
 space_wolff.hpp    |  2 +-
 transformation.hpp | 56 +++++++++++++++++++++++++++++++++++++++++++-----------
 2 files changed, 46 insertions(+), 12 deletions(-)

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 9d529a3..141396a 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -62,7 +62,7 @@ public:
         for (Spin<U, D, S>* s : t->toConsider()) {
           double ΔE = t->ΔE(s);
           if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) {
-            queue.push(t->createNew(s));
+            queue.push(t->createNew(T, visited, s, rng));
           }
         }
 
diff --git a/transformation.hpp b/transformation.hpp
index 0731ec9..fc1833e 100644
--- a/transformation.hpp
+++ b/transformation.hpp
@@ -12,7 +12,7 @@ public:
   virtual std::set<Spin<U, D, S>*> current() const { return {NULL}; }
   virtual std::set<Spin<U, D, S>*> toConsider() const { return {}; }
   virtual double ΔE(const Spin<U, D, S>*) const { return 0; }
-  virtual Transformation* createNew(Spin<U, D, S>*) const { return new Transformation(); }
+  virtual Transformation* createNew(double T, const std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(); }
   virtual void apply(){};
 };
 
@@ -56,7 +56,7 @@ public:
     }
   }
 
-  Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override;
+  Transformation<U, D, R, S>* createNew(double, const std::set<Spin<U, D, S>*>&, Spin<U, D, S>* s, Rng&) const override;
 
   void apply() override {
     M.dict.remove(sOld);
@@ -112,8 +112,8 @@ public:
     }
   }
 
-  Transformation<U, D, R, S>* createNew(Spin<double, D, S>* s) const;
-  Transformation<U, D, R, S>* createNew(Spin<signed, D, S>* s) const;
+  Transformation<U, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const;
+  Transformation<U, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<signed, D, S>* s, Rng&) const;
 
   void apply() override {
     M.dict.remove(s1Old);
@@ -150,11 +150,11 @@ public:
     return M.B(s0s_new) - M.B(s0s_old);
   }
 
-  Transformation<double, D, R, S>* createNew(Spin<double, D, S>* s) const {
+  Transformation<double, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const {
     return new SpinFlip<U, D, R, S>(M, r, s);
   }
 
-  Transformation<signed, D, R, S>* createNew(Spin<signed, D, S>* s) const {
+  Transformation<signed, D, R, S>* createNew(double, const std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
     Vector<signed, 2> v = r.act(*s).x;
     std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
     if (on_site.empty()) {
@@ -168,25 +168,59 @@ public:
 };
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(Spin<U, D, S>* s) const {
+Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {
-    return new SpinFlip(M, r, s);
+    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, r, s);
+    Spin<U, D, S>* sMax = NULL;
+    double ΔEMax = 0.0;
+    for (Spin<U, D, S>* ss : t->toConsider()) {
+      if (ss != NULL && ss != s && !v.contains(ss)) {
+        double ΔE = t->ΔE(ss);
+        if (ΔE > ΔEMax) {
+          sMax = ss;
+          ΔEMax = ΔE;
+        }
+      }
+    }
+    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1)) {
+      return t;
+    } else {
+      delete t;
+      return new PairFlip<U, D, R, S>(M, r, s, sMax);
+    }
   }
 }
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<double, D, S>* s) const {
+Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {
-    return new SpinFlip<U, D, R, S>(M, r, s);
+    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, r, s);
+    Spin<U, D, S>* sMax = NULL;
+    double ΔEMax = 0.0;
+    for (Spin<U, D, S>* ss : t->toConsider()) {
+      if (ss != NULL && ss != s && !v.contains(ss)) {
+        double ΔE = t->ΔE(ss);
+        if (ΔE > ΔEMax) {
+          sMax = ss;
+          ΔEMax = ΔE;
+        }
+      }
+    }
+    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1)) {
+      return t;
+    } else {
+      delete t;
+      return new PairFlip<U, D, R, S>(M, r, s, sMax);
+    }
   }
 }
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<signed, D, S>* s) const {
+Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double, const std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {
-- 
cgit v1.2.3-54-g00ecf


From 9437ed889390e04a1e693a5dfa5290cd90ae150e Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 14 Feb 2020 17:04:17 -0500
Subject: Dynamic pairing method appears to be working with spheres, needs to
 be cleaned up and ising code needs to be updated with new transformation
 interfaces

---
 space_wolff.hpp      |  8 ++++++--
 spheres_infinite.cpp | 27 +++++++++++++++++++++++++--
 transformation.hpp   | 46 ++++++++++++++++++----------------------------
 3 files changed, 49 insertions(+), 32 deletions(-)

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 141396a..63bf118 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -50,19 +50,23 @@ public:
     queue.push(t0);
 
     std::set<Spin<U, D, S>*> visited;
+    std::set<Spin<U, D, S>*> fused;
 
     while (!queue.empty()) {
       Transformation<U, D, R, S>* t = queue.front();
       queue.pop();
 
       std::set<Spin<U, D, S>*> c = t->current();
-      if (!std::any_of(c.begin(), c.end(), [&visited](Spin<U, D, S>* s) { return visited.contains(s); })) {
+      if (!std::any_of(c.begin(), c.end(), [&visited](Spin<U, D, S>* s) { return visited.contains(s); }) && !(c.size() == 1 && std::any_of(c.begin(), c.end(), [&fused](Spin<U, D, S>* s) { return fused.contains(s); }))) {
         visited.insert(c.begin(), c.end());
+        fused.insert(c.begin(), c.end());
 
         for (Spin<U, D, S>* s : t->toConsider()) {
+          if (!fused.contains(s)) {
           double ΔE = t->ΔE(s);
           if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) {
-            queue.push(t->createNew(T, visited, s, rng));
+            queue.push(t->createNew(T, fused, s, rng));
+          }
           }
         }
 
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 5287a84..b39f164 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -117,6 +117,28 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
   };
 }
 
+Gen<double, D, Euclidean<double, D>, double> rGen(double ε) {
+  return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
+    Vector<double, D> t;
+    Matrix<double, D> m;
+
+    double θ = rng.uniform((double)0.0, 2 * M_PI);
+    m(0, 0) = -cos(2 * θ);
+    m(1, 1) = cos(2 * θ);
+    m(0, 1) = -2 * cos(θ) * sin(θ);
+    m(1, 0) = -2 * cos(θ) * sin(θ);
+
+    Spin<double, D, double>* s = rng.pick(M.s);
+    t = M.s0.t;
+    for (unsigned j = 0; j < D; j++) {
+      t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
+    }
+
+    Euclidean<double, D> g(t - m * t, m);
+    return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s);
+  };
+}
+
 int main(int argc, char* argv[]) {
   const unsigned D = 2;
 
@@ -179,8 +201,9 @@ int main(int argc, char* argv[]) {
     return H * s.x.norm();
   };
 
-  auto g1 = eGen(1);
+  auto g1 = eGen(2);
   auto g2 = mGen(0.5);
+  auto g3 = rGen(5);
   animation A(L, 750, argc, argv);
   model sphere(1.0, Z, B);
 
@@ -197,7 +220,7 @@ int main(int argc, char* argv[]) {
     sphere.dict.insert(ss);
   }
 
-  sphere.wolff(T, {g1, g2}, A, N);
+  sphere.wolff(T, {g1, g2, g3}, A, N);
   std::ofstream outfile;
   outfile.open("test.dat");
 
diff --git a/transformation.hpp b/transformation.hpp
index fc1833e..ad3b49f 100644
--- a/transformation.hpp
+++ b/transformation.hpp
@@ -12,7 +12,7 @@ public:
   virtual std::set<Spin<U, D, S>*> current() const { return {NULL}; }
   virtual std::set<Spin<U, D, S>*> toConsider() const { return {}; }
   virtual double ΔE(const Spin<U, D, S>*) const { return 0; }
-  virtual Transformation* createNew(double T, const std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(); }
+  virtual Transformation* createNew(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(); }
   virtual void apply(){};
 };
 
@@ -43,6 +43,8 @@ public:
     neighbors.insert(new_neighbors.begin(), new_neighbors.end());
     neighbors.insert(NULL);
 
+    neighbors.erase(sOld);
+
     return neighbors;
   }
 
@@ -56,7 +58,7 @@ public:
     }
   }
 
-  Transformation<U, D, R, S>* createNew(double, const std::set<Spin<U, D, S>*>&, Spin<U, D, S>* s, Rng&) const override;
+  Transformation<U, D, R, S>* createNew(double, std::set<Spin<U, D, S>*>&, Spin<U, D, S>* s, Rng&) const override;
 
   void apply() override {
     M.dict.remove(sOld);
@@ -97,6 +99,9 @@ public:
     neighbors.insert(new_neighbors_2.begin(), new_neighbors_2.end());
     neighbors.insert(NULL);
 
+    neighbors.erase(s1Old);
+    neighbors.erase(s2Old);
+
     return neighbors;
   }
 
@@ -112,8 +117,8 @@ public:
     }
   }
 
-  Transformation<U, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const;
-  Transformation<U, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<signed, D, S>* s, Rng&) const;
+  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const;
+  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<signed, D, S>* s, Rng&) const;
 
   void apply() override {
     M.dict.remove(s1Old);
@@ -150,11 +155,11 @@ public:
     return M.B(s0s_new) - M.B(s0s_old);
   }
 
-  Transformation<double, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const {
+  Transformation<double, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const {
     return new SpinFlip<U, D, R, S>(M, r, s);
   }
 
-  Transformation<signed, D, R, S>* createNew(double, const std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
+  Transformation<signed, D, R, S>* createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
     Vector<signed, 2> v = r.act(*s).x;
     std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
     if (on_site.empty()) {
@@ -168,33 +173,16 @@ public:
 };
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const {
+Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {
-    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, r, s);
-    Spin<U, D, S>* sMax = NULL;
-    double ΔEMax = 0.0;
-    for (Spin<U, D, S>* ss : t->toConsider()) {
-      if (ss != NULL && ss != s && !v.contains(ss)) {
-        double ΔE = t->ΔE(ss);
-        if (ΔE > ΔEMax) {
-          sMax = ss;
-          ΔEMax = ΔE;
-        }
-      }
-    }
-    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1)) {
-      return t;
-    } else {
-      delete t;
-      return new PairFlip<U, D, R, S>(M, r, s, sMax);
-    }
+    return new SpinFlip<U, D, R, S>(M, r, s);
   }
 }
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const {
+Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {
@@ -210,17 +198,19 @@ Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, const std:
         }
       }
     }
-    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1)) {
+    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || v.contains(sMax)) {
       return t;
     } else {
       delete t;
+      v.insert(s);
+      v.insert(sMax);
       return new PairFlip<U, D, R, S>(M, r, s, sMax);
     }
   }
 }
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double, const std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
+Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {
-- 
cgit v1.2.3-54-g00ecf


From 2ebf2f181edac37bfb932dbb353101e37e97223a Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sun, 16 Feb 2020 21:02:25 -0500
Subject: More work on new transformations, and animation of cluster formation

---
 space_wolff.hpp      |  8 +++---
 spheres_infinite.cpp | 69 ++++++++++++++++++++++++++++++++++++++--------------
 transformation.hpp   | 58 ++++++++++++++++++++++++++++---------------
 3 files changed, 94 insertions(+), 41 deletions(-)

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 63bf118..dc06011 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -18,11 +18,10 @@ template <class U, int D, class R, class S> class Model;
 
 template <class U, int D, class R, class S> class measurement {
 public:
-  virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const R&){};
+  virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const Transformation<double, D, Euclidean<double, D>, double>*){};
   virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*,
                                   const Spin<U, D, S>*, const Spin<U, D, S>&, double){};
-  virtual void plain_site_transformed(const Model<U, D, R, S>&, const Spin<U, D, S>*,
-                                      const Spin<U, D, S>&){};
+  virtual void plain_site_transformed(const Model<U, D, R, S>&, const Transformation<U, D, R, S>&){};
 
   virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&,
                                   const Spin<U, D, S>&, double){};
@@ -70,6 +69,7 @@ public:
           }
         }
 
+        A.plain_site_transformed(*this, *t);
         t->apply();
       }
       delete t;
@@ -84,6 +84,8 @@ public:
       Gen<U, D, R, S>& g = rng.pick(gs);
       Transformation<U, D, R, S> *t = g(*this, rng);
 
+      A.pre_cluster(*this, i, t);
+
       this->step(T, t, A);
 
       A.post_cluster(*this);
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index b39f164..77a3582 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -14,9 +14,10 @@ private:
   uint64_t t2;
   unsigned n;
   unsigned tmp;
+  Euclidean<double, D> s0_tmp;
 
 public:
-  animation(double L, unsigned w, int argc, char* argv[]) {
+  animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0) {
     t1 = 0;
     t2 = 0;
     n = 0;
@@ -31,30 +32,58 @@ public:
     gluOrtho2D(-L, L, -L, L);
   }
 
-  void pre_cluster(const model&, unsigned, const Euclidean<double, D>&) override { tmp = 0; }
-
-  void plain_site_transformed(const model&, const Spin<double, D, double>*,
-                              const Spin<double, D, double>&) override {
+  void plain_site_transformed(const model& m, const Transformation<double, D, Euclidean<double, D>, double>& t) override {
+    if (n > 5000) {
+    if (t.current().size() > 1) {
+      glColor3f(0.0f, 0.0f, 1.0f);
+    } else {
+      glColor3f(0.0f, 1.0f, 0.0f);
+    }
+    for (const Spin<double, 2, double>* s : t.current()) {
+      if (s != NULL) {
+      glBegin(GL_POLYGON);
+      unsigned n_points = 50;
+      for (unsigned i = 0; i < n_points; i++) {
+        glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
+                   s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
+      }
+      glEnd();
+      }
+    }
+    }
     tmp++;
   }
 
-  void post_cluster(const model& m) override {
+  void pre_cluster(const model& m, unsigned, const Transformation<double, D, Euclidean<double, D>, double>* t) override {
+    s0_tmp = m.s0.inverse();
+    tmp = 0;
     glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
     glClear(GL_COLOR_BUFFER_BIT);
+    /*
+    glBegin(GL_LINE);
+    Vector<double, 2> r_center = s0_tmp.act({t->r.t, 0});
+    glEnd();
+    */
     for (const Spin<double, 2, double>* s : m.s) {
       glBegin(GL_POLYGON);
       unsigned n_points = 50;
       glColor3f(0.0f, 0.0f, 0.0f);
       for (unsigned i = 0; i < n_points; i++) {
-        glVertex2d(m.s0.inverse().act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
-                   m.s0.inverse().act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
+        glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
+                   s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
       }
       glEnd();
     }
     glFlush();
+  }
 
+  void post_cluster(const model& m) override {
+    glFlush();
     t1 += tmp;
     t2 += tmp * tmp;
+    if (n > 5000) {
+      sleep(2);
+    }
     n++;
   }
 
@@ -95,17 +124,21 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
 
     Spin<double, D, double>* s1 = rng.pick(M.s);
     Spin<double, D, double>* s2 = rng.pick(M.s);
-
-    while (s1 == s2) {
-      s2 = rng.pick(M.s);
-    }
-
     Vector<double, D> t1 = s1->x;
     Vector<double, D> t2 = s2->x;
     Vector<double, D> t12 = t1 - t2;
+
+    while (s1 == s2 || 1 / t12.norm() < rng.uniform<double>(0, 1)) {
+    s1 = rng.pick(M.s);
+    s2 = rng.pick(M.s);
+    t1 = s1->x;
+    t2 = s2->x;
+    t12 = t1 - t2;
+    }
+
     Vector<double, D> t = (t1 + t2) / 2;
 
-    double θ = atan2(t12[1], t12[0]) + rng.variate<double, std::normal_distribution>(0.0, ε);
+    double θ = atan2(t12[1], t12[0]) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
 
     m(0, 0) = -cos(2 * θ);
     m(1, 1) = cos(2 * θ);
@@ -201,9 +234,9 @@ int main(int argc, char* argv[]) {
     return H * s.x.norm();
   };
 
-  auto g1 = eGen(2);
-  auto g2 = mGen(0.5);
-  auto g3 = rGen(5);
+  auto g1 = eGen(0.5);
+  auto g2 = mGen(0.05);
+  auto g3 = rGen(1);
   animation A(L, 750, argc, argv);
   model sphere(1.0, Z, B);
 
@@ -220,7 +253,7 @@ int main(int argc, char* argv[]) {
     sphere.dict.insert(ss);
   }
 
-  sphere.wolff(T, {g1, g2, g3}, A, N);
+  sphere.wolff(T, {g1, g2}, A, N);
   std::ofstream outfile;
   outfile.open("test.dat");
 
diff --git a/transformation.hpp b/transformation.hpp
index ad3b49f..717a033 100644
--- a/transformation.hpp
+++ b/transformation.hpp
@@ -9,10 +9,12 @@ template <class U, int D, class R, class S> class Model;
 
 template <class U, int D, class R, class S> class Transformation {
 public:
+  const R r;
+  Transformation(const R& r) : r(r) {}
   virtual std::set<Spin<U, D, S>*> current() const { return {NULL}; }
   virtual std::set<Spin<U, D, S>*> toConsider() const { return {}; }
   virtual double ΔE(const Spin<U, D, S>*) const { return 0; }
-  virtual Transformation* createNew(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(); }
+  virtual Transformation* createNew(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(r); }
   virtual void apply(){};
 };
 
@@ -23,11 +25,10 @@ template <class U, int D, class R, class S> class SpinFlip : public Transformati
 private:
   Model<U, D, R, S>& M;
   Spin<U, D, S>* sOld;
-  const R r;
   Spin<U, D, S> sNew;
 
 public:
-  SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s) : M(M), r(r) {
+  SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s) : M(M), Transformation<U, D, R, S>(r) {
     sOld = s;
     sNew = r.act(*s);
   }
@@ -72,12 +73,11 @@ private:
   Model<U, D, R, S>& M;
   Spin<U, D, S>* s1Old;
   Spin<U, D, S>* s2Old;
-  const R r;
   Spin<U, D, S> s1New;
   Spin<U, D, S> s2New;
 
 public:
-  PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2) : M(M), r(r) {
+  PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2) : M(M), Transformation<U, D, R, S>(r) {
     s1Old = s1;
     s2Old = s2;
     s1New = r.act(*s1);
@@ -133,11 +133,10 @@ public:
 template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> {
 private:
   Model<U, D, R, S>& M;
-  const R r;
   R s0New;
 
 public:
-  FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), r(r), s0New(r.act(M.s0)) {}
+  FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), Transformation<U, D, R, S>(r), s0New(r.act(M.s0)) {}
 
   std::set<Spin<U, D, S>*> toConsider() const override {
     std::set<Spin<U, D, S>*> neighbors;
@@ -156,16 +155,16 @@ public:
   }
 
   Transformation<double, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const {
-    return new SpinFlip<U, D, R, S>(M, r, s);
+    return new SpinFlip<U, D, R, S>(M, this->r, s);
   }
 
   Transformation<signed, D, R, S>* createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
-    Vector<signed, 2> v = r.act(*s).x;
+    Vector<signed, 2> v = this->r.act(*s).x;
     std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
     if (on_site.empty()) {
-      return new SpinFlip<U, D, R, S>(M, r, s);
+      return new SpinFlip<U, D, R, S>(M, this->r, s);
     } else {
-      return new PairFlip<U, D, R, S>(M, r, s, *on_site.begin());
+      return new PairFlip<U, D, R, S>(M, this->r, s, *on_site.begin());
     }
   }
 
@@ -175,18 +174,37 @@ public:
 template <class U, int D, class R, class S>
 Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const {
   if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, r);
+    return new FieldFlip<U, D, R, S>(M, this->r);
   } else {
-    return new SpinFlip<U, D, R, S>(M, r, s);
+    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, this->r, s);
+    Spin<U, D, S>* sMax = NULL;
+    double ΔEMax = 0.0;
+    for (Spin<U, D, S>* ss : t->toConsider()) {
+      if (ss != NULL && ss != s && !v.contains(ss)) {
+        double ΔE = t->ΔE(ss);
+        if (ΔE > ΔEMax) {
+          sMax = ss;
+          ΔEMax = ΔE;
+        }
+      }
+    }
+    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || v.contains(sMax)) {
+      return t;
+    } else {
+      delete t;
+      v.insert(s);
+      v.insert(sMax);
+      return new PairFlip<U, D, R, S>(M, this->r, s, sMax);
+    }
   }
 }
 
 template <class U, int D, class R, class S>
 Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const {
   if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, r);
+    return new FieldFlip<U, D, R, S>(M, this->r);
   } else {
-    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, r, s);
+    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, this->r, s);
     Spin<U, D, S>* sMax = NULL;
     double ΔEMax = 0.0;
     for (Spin<U, D, S>* ss : t->toConsider()) {
@@ -204,7 +222,7 @@ Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, std::set<S
       delete t;
       v.insert(s);
       v.insert(sMax);
-      return new PairFlip<U, D, R, S>(M, r, s, sMax);
+      return new PairFlip<U, D, R, S>(M, this->r, s, sMax);
     }
   }
 }
@@ -212,14 +230,14 @@ Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, std::set<S
 template <class U, int D, class R, class S>
 Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
   if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, r);
+    return new FieldFlip<U, D, R, S>(M, this->r);
   } else {
-    Vector<signed, 2> v = r.act(*s).x;
+    Vector<signed, 2> v = this->r.act(*s).x;
     std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
     if (on_site.empty()) {
-      return new SpinFlip<U, D, R, S>(M, r, s);
+      return new SpinFlip<U, D, R, S>(M, this->r, s);
     } else {
-      return new PairFlip<U, D, R, S>(M, r, s, *on_site.begin());
+      return new PairFlip<U, D, R, S>(M, this->r, s, *on_site.begin());
     }
   }
 }
-- 
cgit v1.2.3-54-g00ecf


From 3a3f2dd3b2c47d9d89ed29b7039e39626f2dcf72 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 18 Feb 2020 18:59:05 -0500
Subject: Cleaned up Transformation system and sped things up generally

---
 measurement.hpp      |  20 ++++
 space_wolff.hpp      |  49 +++-------
 spheres_infinite.cpp |  38 ++++----
 transformation.hpp   | 252 +++++++++++++++++++++++----------------------------
 4 files changed, 169 insertions(+), 190 deletions(-)
 create mode 100644 measurement.hpp

diff --git a/measurement.hpp b/measurement.hpp
new file mode 100644
index 0000000..02fde9c
--- /dev/null
+++ b/measurement.hpp
@@ -0,0 +1,20 @@
+#pragma once
+
+#include "spin.hpp"
+#include "transformation.hpp"
+
+template <class U, int D, class R, class S> class Model;
+
+template <class U, int D, class R, class S> class measurement {
+public:
+  virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const Transformation<U, D, R, S>*){};
+  virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*,
+                                  const Spin<U, D, S>*, const Spin<U, D, S>&, double){};
+  virtual void plain_site_transformed(const Model<U, D, R, S>&, const Transformation<U, D, R, S>&){};
+
+  virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&,
+                                  const Spin<U, D, S>&, double){};
+  virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){};
+
+  virtual void post_cluster(const Model<U, D, R, S>&){};
+};
diff --git a/space_wolff.hpp b/space_wolff.hpp
index dc06011..3314ef9 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -8,31 +8,15 @@
 #include <vector>
 
 #include "euclidean.hpp"
+#include "measurement.hpp"
 #include "octree.hpp"
 #include "quantity.hpp"
-#include "spin.hpp"
 #include "random.hpp"
+#include "spin.hpp"
 #include "transformation.hpp"
 
-template <class U, int D, class R, class S> class Model;
-
-template <class U, int D, class R, class S> class measurement {
-public:
-  virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const Transformation<double, D, Euclidean<double, D>, double>*){};
-  virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*,
-                                  const Spin<U, D, S>*, const Spin<U, D, S>&, double){};
-  virtual void plain_site_transformed(const Model<U, D, R, S>&, const Transformation<U, D, R, S>&){};
-
-  virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&,
-                                  const Spin<U, D, S>&, double){};
-  virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){};
-
-  virtual void post_cluster(const Model<U, D, R, S>&){};
-};
-
 template <class U, int D, class R, class S> class Model {
 public:
-  U L;
   R s0;
   std::vector<Spin<U, D, S>*> s;
   Octree<U, D, S> dict;
@@ -42,47 +26,40 @@ public:
 
   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)
-      : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {}
+      : s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {}
 
   void step(double T, Transformation<U, D, R, S>* t0, measurement<U, D, R, S>& A) {
     std::queue<Transformation<U, D, R, S>*> queue;
     queue.push(t0);
 
-    std::set<Spin<U, D, S>*> visited;
-    std::set<Spin<U, D, S>*> fused;
+    std::set<Spin<U, D, S>*> cluster = t0->current();
 
     while (!queue.empty()) {
       Transformation<U, D, R, S>* t = queue.front();
       queue.pop();
 
-      std::set<Spin<U, D, S>*> c = t->current();
-      if (!std::any_of(c.begin(), c.end(), [&visited](Spin<U, D, S>* s) { return visited.contains(s); }) && !(c.size() == 1 && std::any_of(c.begin(), c.end(), [&fused](Spin<U, D, S>* s) { return fused.contains(s); }))) {
-        visited.insert(c.begin(), c.end());
-        fused.insert(c.begin(), c.end());
-
-        for (Spin<U, D, S>* s : t->toConsider()) {
-          if (!fused.contains(s)) {
+      for (Spin<U, D, S>* s : t->toConsider()) {
+        if (!cluster.contains(s)) {
           double ΔE = t->ΔE(s);
           if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) {
-            queue.push(t->createNew(T, fused, s, rng));
-          }
+            cluster.insert(s);
+            queue.push(t->createNew(T, cluster, s, rng));
           }
         }
-
-        A.plain_site_transformed(*this, *t);
-        t->apply();
       }
+
+      A.plain_site_transformed(*this, *t);
+      t->apply();
       delete t;
     }
   }
 
   void resize(double T, double P) {}
 
-  void wolff(double T, std::vector<Gen<U, D, R, S>> gs,
-             measurement<U, D, R, S>& A, unsigned N) {
+  void wolff(double T, std::vector<Gen<U, D, R, S>> gs, measurement<U, D, R, S>& A, unsigned N) {
     for (unsigned i = 0; i < N; i++) {
       Gen<U, D, R, S>& g = rng.pick(gs);
-      Transformation<U, D, R, S> *t = g(*this, rng);
+      Transformation<U, D, R, S>* t = g(*this, rng);
 
       A.pre_cluster(*this, i, t);
 
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 77a3582..e41ea17 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -14,10 +14,11 @@ private:
   uint64_t t2;
   unsigned n;
   unsigned tmp;
+  unsigned wait;
   Euclidean<double, D> s0_tmp;
 
 public:
-  animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0) {
+  animation(double L, unsigned w, int argc, char* argv[]) : s0_tmp(0), wait(1000) {
     t1 = 0;
     t2 = 0;
     n = 0;
@@ -25,7 +26,7 @@ public:
     glutInit(&argc, argv);
     glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
     glutInitWindowSize(w, w);
-    glutCreateWindow("wolff");
+    glutCreateWindow("wolffWindow");
     glClearColor(0.0, 0.0, 0.0, 0.0);
     glMatrixMode(GL_PROJECTION);
     glLoadIdentity();
@@ -33,7 +34,7 @@ public:
   }
 
   void plain_site_transformed(const model& m, const Transformation<double, D, Euclidean<double, D>, double>& t) override {
-    if (n > 5000) {
+    if (n > wait) {
     if (t.current().size() > 1) {
       glColor3f(0.0f, 0.0f, 1.0f);
     } else {
@@ -81,7 +82,7 @@ public:
     glFlush();
     t1 += tmp;
     t2 += tmp * tmp;
-    if (n > 5000) {
+    if (n > wait) {
       sleep(2);
     }
     n++;
@@ -124,28 +125,31 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
 
     Spin<double, D, double>* s1 = rng.pick(M.s);
     Spin<double, D, double>* s2 = rng.pick(M.s);
-    Vector<double, D> t1 = s1->x;
-    Vector<double, D> t2 = s2->x;
-    Vector<double, D> t12 = t1 - t2;
 
-    while (s1 == s2 || 1 / t12.norm() < rng.uniform<double>(0, 1)) {
-    s1 = rng.pick(M.s);
-    s2 = rng.pick(M.s);
-    t1 = s1->x;
-    t2 = s2->x;
-    t12 = t1 - t2;
+    while (s1 == s2) {
+      s2 = rng.pick(M.s);
     }
 
+    Vector<double, D> t1 = s1->x;
+    Vector<double, D> t2 = s2->x;
+    Vector<double, D> t12 = t1 - t2;
     Vector<double, D> t = (t1 + t2) / 2;
 
-    double θ = atan2(t12[1], t12[0]) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
+    double θ = atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
 
     m(0, 0) = -cos(2 * θ);
     m(1, 1) = cos(2 * θ);
     m(0, 1) = -2 * cos(θ) * sin(θ);
     m(1, 0) = -2 * cos(θ) * sin(θ);
 
-    Euclidean<double, D> g(t - m * t, m);
+    Vector<double, D> t3 = t - m * t;
+
+    if (t3(0) != t3(0)) {
+      std::cout << t3 << "\n" << t << "\n" << m << "\n" << t12 << "\n" ;
+      getchar();
+    }
+
+    Euclidean<double, D> g(t3, m);
     return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, s1, s2);
   };
 }
@@ -234,8 +238,8 @@ int main(int argc, char* argv[]) {
     return H * s.x.norm();
   };
 
-  auto g1 = eGen(0.5);
-  auto g2 = mGen(0.05);
+  auto g1 = eGen(1);
+  auto g2 = mGen(0.1);
   auto g3 = rGen(1);
   animation A(L, 750, argc, argv);
   model sphere(1.0, Z, B);
diff --git a/transformation.hpp b/transformation.hpp
index 717a033..b37b79e 100644
--- a/transformation.hpp
+++ b/transformation.hpp
@@ -9,26 +9,63 @@ template <class U, int D, class R, class S> class Model;
 
 template <class U, int D, class R, class S> class Transformation {
 public:
+  Model<U, D, R, S>& M;
   const R r;
-  Transformation(const R& r) : r(r) {}
+
+  Transformation(Model<U, D, R, S>& M, const R& r) : r(r), M(M) {}
   virtual std::set<Spin<U, D, S>*> current() const { return {NULL}; }
   virtual std::set<Spin<U, D, S>*> toConsider() const { return {}; }
   virtual double ΔE(const Spin<U, D, S>*) const { return 0; }
-  virtual Transformation* createNew(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(r); }
+  virtual Transformation* createNew(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*,
+                                    Rng& rng) const {
+    return new Transformation(M, r);
+  }
   virtual void apply(){};
 };
 
 template <class U, int D, class R, class S>
 using Gen = std::function<Transformation<U, D, R, S>*(Model<U, D, R, S>&, Rng&)>;
 
+template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> {
+private:
+  R s0New;
+
+public:
+  FieldFlip(Model<U, D, R, S>& M, const R& r)
+      : Transformation<U, D, R, S>(M, r), s0New(r.act(M.s0)) {}
+
+  std::set<Spin<U, D, S>*> toConsider() const override {
+    std::set<Spin<U, D, S>*> neighbors;
+
+    for (Spin<U, D, S>* s : this->M.s) {
+      neighbors.insert(s);
+    }
+
+    return neighbors;
+  }
+
+  double ΔE(const Spin<U, D, S>* s) const override {
+    Spin<U, D, S> s0s_old = this->M.s0.inverse().act(*s);
+    Spin<U, D, S> s0s_new = s0New.inverse().act(*s);
+    return this->M.B(s0s_new) - this->M.B(s0s_old);
+  }
+
+  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster,
+                                             Spin<U, D, S>* s, Rng& rng) const override {
+    return pairGenerator<D, R, S>(this, T, cluster, s, rng);
+  }
+
+  void apply() override { this->M.s0 = s0New; }
+};
+
 template <class U, int D, class R, class S> class SpinFlip : public Transformation<U, D, R, S> {
 private:
-  Model<U, D, R, S>& M;
   Spin<U, D, S>* sOld;
   Spin<U, D, S> sNew;
 
 public:
-  SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s) : M(M), Transformation<U, D, R, S>(r) {
+  SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s)
+      : Transformation<U, D, R, S>(M, r) {
     sOld = s;
     sNew = r.act(*s);
   }
@@ -36,13 +73,12 @@ public:
   std::set<Spin<U, D, S>*> current() const override { return {sOld}; }
 
   std::set<Spin<U, D, S>*> toConsider() const override {
-    std::set<Spin<U, D, S>*> neighbors;
-    std::set<Spin<U, D, S>*> current_neighbors = M.dict.neighbors(sOld->x);
-    std::set<Spin<U, D, S>*> new_neighbors = M.dict.neighbors(sNew.x);
+    std::set<Spin<U, D, S>*> neighbors = {NULL};
+    std::set<Spin<U, D, S>*> oldNeighbors = this->M.dict.neighbors(sOld->x);
+    std::set<Spin<U, D, S>*> newNeighbors = this->M.dict.neighbors(sNew.x);
 
-    neighbors.insert(current_neighbors.begin(), current_neighbors.end());
-    neighbors.insert(new_neighbors.begin(), new_neighbors.end());
-    neighbors.insert(NULL);
+    neighbors.insert(oldNeighbors.begin(), oldNeighbors.end());
+    neighbors.insert(newNeighbors.begin(), newNeighbors.end());
 
     neighbors.erase(sOld);
 
@@ -51,33 +87,40 @@ public:
 
   double ΔE(const Spin<U, D, S>* s) const override {
     if (s == NULL) {
-      Spin<U, D, S> s0s_old = M.s0.inverse().act(*sOld);
-      Spin<U, D, S> s0s_new = M.s0.inverse().act(sNew);
-      return M.B(s0s_new) - M.B(s0s_old);
+      Spin<U, D, S> s0s_old = this->M.s0.inverse().act(*sOld);
+      Spin<U, D, S> s0s_new = this->M.s0.inverse().act(sNew);
+      return this->M.B(s0s_new) - this->M.B(s0s_old);
     } else {
-      return M.Z(*sOld, *s) - M.Z(sNew, *s);
+      return this->M.Z(*sOld, *s) - this->M.Z(sNew, *s);
     }
   }
 
-  Transformation<U, D, R, S>* createNew(double, std::set<Spin<U, D, S>*>&, Spin<U, D, S>* s, Rng&) const override;
+  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster,
+                                             Spin<U, D, S>* s, Rng& rng) const override {
+    if (s == NULL) {
+      return new FieldFlip<U, D, R, S>(this->M, this->r);
+    } else {
+      return pairGenerator<D, R, S>(this, T, cluster, s, rng);
+    }
+  }
 
   void apply() override {
-    M.dict.remove(sOld);
+    this->M.dict.remove(sOld);
     *sOld = sNew;
-    M.dict.insert(sOld);
+    this->M.dict.insert(sOld);
   }
 };
 
 template <class U, int D, class R, class S> class PairFlip : public Transformation<U, D, R, S> {
 private:
-  Model<U, D, R, S>& M;
   Spin<U, D, S>* s1Old;
   Spin<U, D, S>* s2Old;
   Spin<U, D, S> s1New;
   Spin<U, D, S> s2New;
 
 public:
-  PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2) : M(M), Transformation<U, D, R, S>(r) {
+  PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2)
+      : Transformation<U, D, R, S>(M, r) {
     s1Old = s1;
     s2Old = s2;
     s1New = r.act(*s1);
@@ -87,17 +130,17 @@ public:
   std::set<Spin<U, D, S>*> current() const override { return {s1Old, s2Old}; }
 
   std::set<Spin<U, D, S>*> toConsider() const override {
-    std::set<Spin<U, D, S>*> neighbors;
-    std::set<Spin<U, D, S>*> current_neighbors_1 = M.dict.neighbors(s1Old->x);
-    std::set<Spin<U, D, S>*> current_neighbors_2 = M.dict.neighbors(s2Old->x);
-    std::set<Spin<U, D, S>*> new_neighbors_1 = M.dict.neighbors(s1New.x);
-    std::set<Spin<U, D, S>*> new_neighbors_2 = M.dict.neighbors(s2New.x);
+    std::set<Spin<U, D, S>*> neighbors = {NULL};
+
+    std::set<Spin<U, D, S>*> current_neighbors_1 = this->M.dict.neighbors(s1Old->x);
+    std::set<Spin<U, D, S>*> current_neighbors_2 = this->M.dict.neighbors(s2Old->x);
+    std::set<Spin<U, D, S>*> new_neighbors_1 = this->M.dict.neighbors(s1New.x);
+    std::set<Spin<U, D, S>*> new_neighbors_2 = this->M.dict.neighbors(s2New.x);
 
     neighbors.insert(current_neighbors_1.begin(), current_neighbors_1.end());
     neighbors.insert(current_neighbors_2.begin(), current_neighbors_2.end());
     neighbors.insert(new_neighbors_1.begin(), new_neighbors_1.end());
     neighbors.insert(new_neighbors_2.begin(), new_neighbors_2.end());
-    neighbors.insert(NULL);
 
     neighbors.erase(s1Old);
     neighbors.erase(s2Old);
@@ -107,137 +150,72 @@ public:
 
   double ΔE(const Spin<U, D, S>* s) const override {
     if (s == NULL) {
-      Spin<U, D, S> s0s1_old = M.s0.inverse().act(*s1Old);
-      Spin<U, D, S> s0s1_new = M.s0.inverse().act(s1New);
-      Spin<U, D, S> s0s2_old = M.s0.inverse().act(*s2Old);
-      Spin<U, D, S> s0s2_new = M.s0.inverse().act(s2New);
-      return M.B(s0s1_new) + M.B(s0s2_new) - M.B(s0s1_old) - M.B(s0s2_old);
+      Spin<U, D, S> s0s1_old = this->M.s0.inverse().act(*s1Old);
+      Spin<U, D, S> s0s1_new = this->M.s0.inverse().act(s1New);
+      Spin<U, D, S> s0s2_old = this->M.s0.inverse().act(*s2Old);
+      Spin<U, D, S> s0s2_new = this->M.s0.inverse().act(s2New);
+      return this->M.B(s0s1_new) + this->M.B(s0s2_new) - this->M.B(s0s1_old) - this->M.B(s0s2_old);
     } else {
-      return M.Z(*s1Old, *s) + M.Z(*s2Old, *s) - M.Z(s1New, *s) - M.Z(s2New, *s);
+      return this->M.Z(*s1Old, *s) + this->M.Z(*s2Old, *s) - this->M.Z(s1New, *s) - this->M.Z(s2New, *s);
     }
   }
 
-  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const;
-  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<signed, D, S>* s, Rng&) const;
+  Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& cluster,
+                                             Spin<U, D, S>* s, Rng& rng) const override {
+    if (s == NULL) {
+      return new FieldFlip<U, D, R, S>(this->M, this->r);
+    } else {
+      return pairGenerator<D, R, S>(this, T, cluster, s, rng);
+    }
+  }
 
   void apply() override {
-    M.dict.remove(s1Old);
-    M.dict.remove(s2Old);
+    this->M.dict.remove(s1Old);
+    this->M.dict.remove(s2Old);
     *s1Old = s1New;
     *s2Old = s2New;
-    M.dict.insert(s1Old);
-    M.dict.insert(s2Old);
+    this->M.dict.insert(s1Old);
+    this->M.dict.insert(s2Old);
   }
 };
 
-template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> {
-private:
-  Model<U, D, R, S>& M;
-  R s0New;
-
-public:
-  FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), Transformation<U, D, R, S>(r), s0New(r.act(M.s0)) {}
-
-  std::set<Spin<U, D, S>*> toConsider() const override {
-    std::set<Spin<U, D, S>*> neighbors;
-
-    for (Spin<U, D, S>* s : M.s) {
-      neighbors.insert(s);
-    }
-
-    return neighbors;
-  }
-
-  double ΔE(const Spin<U, D, S>* s) const override {
-    Spin<U, D, S> s0s_old = M.s0.inverse().act(*s);
-    Spin<U, D, S> s0s_new = s0New.inverse().act(*s);
-    return M.B(s0s_new) - M.B(s0s_old);
-  }
-
-  Transformation<double, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const {
-    return new SpinFlip<U, D, R, S>(M, this->r, s);
-  }
-
-  Transformation<signed, D, R, S>* createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
-    Vector<signed, 2> v = this->r.act(*s).x;
-    std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
-    if (on_site.empty()) {
-      return new SpinFlip<U, D, R, S>(M, this->r, s);
-    } else {
-      return new PairFlip<U, D, R, S>(M, this->r, s, *on_site.begin());
-    }
-  }
-
-  void apply() override { M.s0 = s0New; }
-};
 
-template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const {
-  if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, this->r);
-  } else {
-    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, this->r, s);
-    Spin<U, D, S>* sMax = NULL;
-    double ΔEMax = 0.0;
-    for (Spin<U, D, S>* ss : t->toConsider()) {
-      if (ss != NULL && ss != s && !v.contains(ss)) {
-        double ΔE = t->ΔE(ss);
-        if (ΔE > ΔEMax) {
-          sMax = ss;
-          ΔEMax = ΔE;
-        }
+template <int D, class R, class S>
+Transformation<double, D, R, S>*
+pairGenerator(const Transformation<double, D, R, S>* tOld, double T,
+                        std::set<Spin<double, D, S>*>& cluster, Spin<double, D, S>* s, Rng& rng) {
+  SpinFlip<double, D, R, S>* t = new SpinFlip<double, D, R, S>(tOld->M, tOld->r, s);
+  Spin<double, D, S>* sMax = NULL;
+  double ΔEMax = 0.0;
+  for (Spin<double, D, S>* ss : t->toConsider()) {
+    if (ss != NULL && ss != s && !cluster.contains(ss)) {
+      double ΔE = t->ΔE(ss);
+      if (ΔE > ΔEMax) {
+        sMax = ss;
+        ΔEMax = ΔE;
       }
     }
-    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || v.contains(sMax)) {
-      return t;
-    } else {
-      delete t;
-      v.insert(s);
-      v.insert(sMax);
-      return new PairFlip<U, D, R, S>(M, this->r, s, sMax);
-    }
   }
-}
-
-template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const {
-  if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, this->r);
+  if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || cluster.contains(sMax)) {
+    return t;
   } else {
-    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, this->r, s);
-    Spin<U, D, S>* sMax = NULL;
-    double ΔEMax = 0.0;
-    for (Spin<U, D, S>* ss : t->toConsider()) {
-      if (ss != NULL && ss != s && !v.contains(ss)) {
-        double ΔE = t->ΔE(ss);
-        if (ΔE > ΔEMax) {
-          sMax = ss;
-          ΔEMax = ΔE;
-        }
-      }
-    }
-    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || v.contains(sMax)) {
-      return t;
-    } else {
-      delete t;
-      v.insert(s);
-      v.insert(sMax);
-      return new PairFlip<U, D, R, S>(M, this->r, s, sMax);
-    }
+    delete t;
+    cluster.insert(sMax);
+    return new PairFlip<double, D, R, S>(tOld->M, tOld->r, s, sMax);
   }
 }
 
-template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
-  if (s == NULL) {
-    return new FieldFlip<U, D, R, S>(M, this->r);
+template <int D, class R, class S>
+Transformation<signed, D, R, S>*
+pairGenerator(const Transformation<signed, D, R, S>* tOld, double T,
+                        std::set<Spin<signed, D, S>*>& cluster, Spin<signed, D, S>* s, Rng& rng) {
+  Vector<signed, 2> v = tOld->r.act(*s).x;
+  std::set<Spin<signed, D, S>*> on_site = tOld->M.dict.at(v);
+  if (on_site.empty()) {
+    return new SpinFlip<signed, D, R, S>(tOld->M, tOld->r, s);
   } else {
-    Vector<signed, 2> v = this->r.act(*s).x;
-    std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
-    if (on_site.empty()) {
-      return new SpinFlip<U, D, R, S>(M, this->r, s);
-    } else {
-      return new PairFlip<U, D, R, S>(M, this->r, s, *on_site.begin());
-    }
+    cluster.insert(*on_site.begin());
+    return new PairFlip<signed, D, R, S>(tOld->M, tOld->r, s, *on_site.begin());
   }
 }
+
-- 
cgit v1.2.3-54-g00ecf