1 files changed, 54 insertions, 2 deletions

diff --git a/spheres.hpp b/spheres.hpp
index 161ef9c..dd2a189 100644
--- a/spheres.hpp
+++ b/spheres.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <assert.h>
+#include <random>
 #include <set>
 #include <vector>
 #include <queue>
@@ -224,7 +225,6 @@ public:
   }
 };
 
-
 template <int D, Particle<D> T> class Model {
 private:
   bool metropolisAccept(double β, double ΔE, Rng& r, double a = 1.0) const {
@@ -356,7 +356,7 @@ public:
         }
 
         if (pMax != NULL) {
-          if (1 - 1e2 * exp(-β * ΔEmax) > r.uniform(0.0, 1.0)) {
+          if (1 - exp(-β * ΔEmax) > r.uniform(0.0, 1.0)) {
             m.insert(pMax);
           }
         }
@@ -404,6 +404,58 @@ public:
     return clusterFlip(β, g, m, r);
   }
 
+  unsigned randomCluster(double β, Rng& r) {
+    Vector<2> t = {r.uniform(0.0, 1.0), r.uniform(0.0, 1.0)};
+
+    Matrix<2> mat;
+
+    mat(0, 0) = -1;
+    mat(1, 1) = -1;
+    mat(0, 1) = 0;
+    mat(1, 0) = 0;
+
+    Euclidean<2> g(t - mat * t, mat);
+
+    Move<D, T> m(g, {&r.pick(particles)});
+
+    return clusterFlip(β, g, m, r);
+  }
+
+  unsigned clusterNudge(double β, T& s, Rng& r) {
+    Vector<2> t = s.x + (Vector<2>){r.variate<double, std::normal_distribution>(0, s.r / 100),
+        r.variate<double, std::normal_distribution>(0, s.r / 100)}
+      ;
+
+    Matrix<2> mat;
+
+    mat(0, 0) = -1;
+    mat(1, 1) = -1;
+    mat(0, 1) = 0;
+    mat(1, 0) = 0;
+
+    Euclidean<2> g(t - mat * t, mat);
+
+    Move<D, T> m(g, {&s});
+
+    return clusterFlip(β, g, m, r);
+  }
+
+  unsigned clusterNeighborSwap(double β, T& s, Rng& r) {
+    T* s2 = &s;
+    while (s2 == &s) {
+      s2 = r.pick(nl.neighborsOf(s.x));
+    }
+    return clusterSwap(β, s, *s2, r);
+  }
+
+  unsigned randomClusterNeighborSwap(double β, Rng& r) {
+    return clusterNeighborSwap(β, r.pick(particles), r);
+  }
+
+  unsigned randomClusterNudge(double β, Rng& r) {
+    return clusterNudge(β, r.pick(particles), r);
+  }
+
   unsigned randomClusterSwap(double β, Rng& r) {
     return clusterSwap(β, r.pick(particles), r.pick(particles), r);
   }