Changed the RNG and implemented Transformations as the method for parameterizing movement.

4 files changed, 189 insertions, 85 deletions

diff --git a/.gitmodules b/.gitmodules
index 0a97926..15de333 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,3 +1,6 @@
 [submodule "randutils"]
 	path = randutils
 	url = https://gist.github.com/imneme/540829265469e673d045/
+[submodule "pcg-cpp"]
+	path = pcg-cpp
+	url = https://github.com/imneme/pcg-cpp
diff --git a/pcg-cpp b/pcg-cpp
new file mode 160000
+Subproject 5b5cac8d61339e810c5dbb4692d868a1d7ca1b2
diff --git a/space_wolff.hpp b/space_wolff.hpp
index 4c935db..8ea7dfb 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -7,6 +7,7 @@
 #include <queue>
 #include <vector>
 
+#include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
 
 #include "euclidean.hpp"
@@ -14,6 +15,8 @@
 #include "quantity.hpp"
 #include "spin.hpp"
 
+using Rng = randutils::random_generator<pcg32>;
+
 template <class U, int D, class R, class S> class Model;
 
 template <class U, int D, class R, class S> class measurement {
@@ -31,42 +34,171 @@ public:
   virtual void post_cluster(const Model<U, D, R, S>&){};
 };
 
-template <class U, int D, class R, class S> using Gen = std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)>;
-
-template <class U, int D, class R, class S> class Model;
 
 template <class U, int D, class R, class S> class Transformation {
   public:
-    virtual void ready() {};
-    virtual double p(const Transformation&) const { return 0; }
-    virtual void apply() {};
-    virtual std::set<Transformation*> to_consider() const {};
+    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 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 SpinFlip : public Transformation<U, D, R, S> {
   private:
     Model<U, D, R, S>& M;
-    Spin<U, D, S>& s;
-    const R& r;
-    Spin<U, D, S> s_new;
+    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), sOld(s) {
+      sNew = r.act(s);
+    }
+
+    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);
+
+      neighbors.insert(current_neighbors.begin(), current_neighbors.end());
+      neighbors.insert(new_neighbors.begin(), new_neighbors.end());
+      neighbors.insert(NULL);
+
+      return neighbors;
+    }
+
+    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);
+      } else {
+        return M.Z(sOld, *s) - M.Z(sNew, *s);
+      }
+    }
+
+    Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override;
+
+    void apply() override {
+      M.dict.remove(&sOld);
+      sOld = sNew;
+      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;
+    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), s1Old(s1), s2Old(s2) {
+      s1New = r.act(s1);
+      s2New = r.act(s2);
+    }
+
+    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);
+
+      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);
+
+      return neighbors;
+    }
+
+    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);
+      } else {
+        return M.Z(s1Old, *s) + M.Z(s2Old, *s) - M.Z(s1New, *s) - M.Z(s2New, *s);
+      }
+    }
+
+    Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override;
+
+    void apply() override {
+      M.dict.remove(&s1Old);
+      M.dict.remove(&s2Old);
+      s1Old = s1New;
+      s2Old = s2New;
+      M.dict.insert(&s1Old);
+      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;
+    const R r;
+    R s0New;
   public:
-    SpinFlip(Model<U, D, R, S>& M, Spin<U, D, S>& s, const R& r) : M(M), s(s), r(r) {}
+    FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), r(r), s0New(r.act(M.s0)) {}
+
+    std::set<Spin<U, D, S>*> toConsider() const override {
+      std::set<Spin<U, D, S>*> neighbors;
 
-    void ready() override {
-      s_new = r.act(s);
-    };
+      for (Spin<U, D, S>& s : M.s) {
+        neighbors.insert(&s);
+      }
 
-    double p(const SpinFlip& sf, double T) const override {
-      return 1.0 - exp(-(M.Z(s, sf.s) - M.Z(s_new, sf.s)) / T);
+      return neighbors;
     }
 
-    double p(const SpinFlip& sf, double T) const override {
-      return 1.0 - exp(-(M.Z(s, sf.s) - M.Z(s_new, sf.s)) / T);
+    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<U, D, R, S>* createNew(Spin<U, D, S>* s) const override {
+      return new SpinFlip<U, D, R, S>(M, r, *s);
+    }
 
+    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(Spin<U, D, S>* s) const {
+  if (s == NULL) {
+    return new FieldFlip<U, D, R, S>(M, r);
+  } else {
+    return new SpinFlip(M, r, *s);
+  }
+}
+
+template<class U, int D, class R, class S> Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<U, D, S>* s) const {
+  if (s == NULL) {
+    return new FieldFlip<U, D, R, S>(M, r);
+  } else {
+    return new SpinFlip<U, D, R, S>(M, r, *s);
+  }
+}
+
 template <class U, int D, class R, class S> class Model {
 public:
   U L;
@@ -75,70 +207,36 @@ public:
   Octree<U, D, S> dict;
   std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z;
   std::function<double(const Spin<U, D, S>&)> B;
-  randutils::mt19937_rng rng;
+  Rng rng;
 
   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)) {}
 
-  void step(double T, unsigned ind, R& r, measurement<U, D, R, S>& A) {
-    std::queue<Spin<U, D, S>*> queue;
-    queue.push(&(s[ind]));
+  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;
 
     while (!queue.empty()) {
-      Spin<U, D, S>* si = queue.front();
+      Transformation<U, D, R, S>* t = queue.front();
       queue.pop();
 
-      if (!visited.contains(si)) {
-        visited.insert(si);
-
-        if (si == NULL) {
-          R 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);
-            A.ghost_bond_visited(*this, s0s_old, s0s_new, p);
-            if (rng.uniform(0.0, 1.0) < p) {
-              queue.push(&ss);
-            }
-          }
-          A.ghost_site_transformed(*this, s0_new);
-          s0 = s0_new;
-        } 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);
-          std::set<Spin<U, D, S>*> new_neighbors = dict.neighbors(si_new.x);
-
-          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;
-              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);
-                A.ghost_bond_visited(*this, s0s_old, s0s_new, p);
-              } else {
-                p = 1.0 - exp(-(Z(*si, *sj) - Z(si_new, *sj)) / T);
-                A.plain_bond_visited(*this, si, sj, si_new, p);
-              }
-              if (rng.uniform(0.0, 1.0) < p) {
-                queue.push(sj);
-              }
-            }
+      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); })) {
+        visited.insert(c.begin(), c.end());
+
+        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));
           }
-          A.plain_site_transformed(*this, si, si_new);
-          dict.remove(si);
-          *si = si_new;
-          dict.insert(si);
         }
+
+        t->apply();
       }
+      delete t;
     }
   }
 
@@ -148,12 +246,9 @@ public:
              measurement<U, D, R, S>& A, unsigned N) {
     for (unsigned i = 0; i < N; i++) {
       Gen<U, D, R, S>& g = rng.pick(gs);
-      R r = g(*this, rng);
-      unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1));
-
-      A.pre_cluster(*this, ind, r);
+      Transformation<U, D, R, S> *t = g(*this, rng);
 
-      this->step(T, ind, r, A);
+      this->step(T, t, A);
 
       A.post_cluster(*this);
     }
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index a9cd7d0..3505d6d 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -68,7 +68,7 @@ public:
 };
 
 Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
-  return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+  return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
     Vector<double, D> t;
     Matrix<double, D> m;
 
@@ -85,12 +85,12 @@ Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
     }
 
     Euclidean<double, D> g(t - m * t, m);
-    return g;
+    return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, M.s[f_ind]);
   };
 }
 
 Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
-  return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+  return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
     Matrix<double, D> m;
 
     unsigned f_ind1 = rng.uniform((unsigned)0, (unsigned)M.s.size());
@@ -110,7 +110,7 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
     m(1, 0) = -2 * cos(θ) * sin(θ);
 
     Euclidean<double, D> g(t - m * t, m);
-    return g;
+    return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, M.s[f_ind1], M.s[f_ind2]);
   };
 }
 
@@ -123,9 +123,12 @@ int main(int argc, char* argv[]) {
   double H = 1.0;
   unsigned n = 25;
 
+  double k = 1e2;
+  double a = 0.1;
+
   int opt;
 
-  while ((opt = getopt(argc, argv, "n:N:L:T:H:")) != -1) {
+  while ((opt = getopt(argc, argv, "n:N:L:T:H:a:k:")) != -1) {
     switch (opt) {
     case 'n':
       n = (unsigned)atof(optarg);
@@ -142,14 +145,17 @@ int main(int argc, char* argv[]) {
     case 'H':
       H = atof(optarg);
       break;
+    case 'a':
+      a = atof(optarg);
+      break;
+    case 'k':
+      k = atof(optarg);
+      break;
     default:
       exit(1);
     }
   }
 
-  double k = 1000;
-  double a = 0.05;
-
   std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
       [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
     Vector<double, D> d = s1.x - s2.x;
@@ -170,19 +176,19 @@ int main(int argc, char* argv[]) {
     return H * s.x.norm();
   };
 
-  auto g1 = eGen(0.25);
-  auto g2 = mGen(0.1);
+  auto g1 = eGen(0.5);
+  auto g2 = mGen(0.2);
   animation A(L, 750, argc, argv);
   model sphere(1, Z, B);
 
-  randutils::mt19937_rng rng;
+  Rng rng;
 
   sphere.s.reserve(n);
 
   unsigned nx = floor(sqrt(n));
   for (unsigned i = 0; i < n; i++) {
     Vector<double, D> pos = {(i / nx) * L / nx, (i % nx) * L / nx};
-    sphere.s.push_back({pos, rng.uniform<double>(0.45, 0.45)});
+    sphere.s.push_back({pos, rng.uniform<double>(0.25, 0.45)});
     sphere.dict.insert(&sphere.s.back());
   }