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

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 {