refactoring and generalization in preparation for pressure ensemble

2 files changed, 144 insertions, 120 deletions

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 02c544d..1335f07 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -244,125 +244,136 @@ public:
   unsigned num_added() const { return n - wait; }
 };
 
-template <class U, int D, class S> class Model;
+template <class U, int D, class R, class S> class Model;
 
-template <class U, int D, class S> class measurement {
+template <class U, int D, class R, 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 Spin<U, D, S>*, const Spin<U, D, S>*,
+  virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const R&){};
+  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, S>&, const Spin<U, D, S>*,
+  virtual void plain_site_transformed(const Model<U, D, R, S>&, const Spin<U, D, S>*,
                                       const Spin<U, D, S>&){};
 
-  virtual void ghost_bond_visited(const Model<U, D, S>&, const Spin<U, D, S>&, const Spin<U, D, 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, S>&, const Euclidean<U, D>&){};
+  virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){};
 
-  virtual void post_cluster(const Model<U, D, S>&){};
+  virtual void post_cluster(const Model<U, D, R, S>&){};
 };
 
-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;
-  measurement<U, D, S>& A;
-
-  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<int D>
+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(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z,
-        std::function<double(const Spin<U, D, S>&)> B, unsigned dDepth, unsigned nDepth,
-        measurement<U, D, S>& A)
-      : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), dDepth(dDepth), nDepth(nDepth), A(A) {
-    std::array<double, D> ini_sequence;
-    ini_sequence.fill(1);
-    std::list<std::array<double, D>> sequences;
-    sequences.push_back(ini_sequence);
+template <class U, int D>
+std::vector<Vector<U, D>> torus_vecs(U L) {
+  std::vector<Vector<U, D>> vecs;
+  std::array<double, D> ini_sequence;
+  ini_sequence.fill(1);
+  std::list<std::array<double, D>> sequences;
+  sequences.push_back(ini_sequence);
+  sequences.pop_front(); // don't want the identity matrix!
+ 
+  for (std::array<double, D> sequence : sequences) {
+    Vector<U, D> v;
+    for (unsigned i = 0; i < D; i++) {
+      if (sequence[i] == 1) {
+        v(i) = 0;
+       } else {
+        v(i) = L / 2;
+       }
+     }
+
+     vecs.push_back(v);
+   }
 
-    one_sequences(sequences, D);
+  return vecs;
+}
 
-    sequences.pop_front(); // don't want the identity matrix!
+template <class U, int D>
+std::vector<Matrix<U, D>> torus_mats() {
+  std::vector<Matrix<U, D>> mats;
 
-    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;
-          }
-        }
-      }
+  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;
+  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 {
-          v(i) = L / 2;
+          m(i, j) = 0;
         }
       }
-
-      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;
-                }
+    mats.push_back(m);
+  }
+
+  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) = 0;
+                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]);
-      }
+  return mats;
+}
 
-      E -= B(s0.inverse().act(s[i]));
-    }
-  }
+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;
+  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::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g;
+  randutils::mt19937_rng rng;
+  measurement<U, D, R, S>& A;
+
+
+  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, std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g, unsigned dDepth, unsigned nDepth,
+        measurement<U, D, R, S>& A)
+      : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), g(g), dDepth(dDepth), nDepth(nDepth), rng(), A(A) {}
 
-  void step(double T, unsigned ind, Euclidean<U, D> r, std::mt19937& rng) {
+  void step(double T, unsigned ind, Euclidean<U, D> r) {
     unsigned cluster_size = 0;
     std::uniform_real_distribution<double> dist(0.0, 1.0);
 
@@ -427,37 +438,11 @@ public:
     }
   }
 
-  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);
+  void wolff(double T, unsigned N) {
     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);
-
+      R g = g(*this, rng);
       unsigned ind = ind_dist(rng);
 
       A.pre_cluster(*this, ind, g);
@@ -465,8 +450,6 @@ public:
       this->step(T, ind, g, rng);
 
       A.post_cluster(*this);
-
-      this->update_energy();
     }
   }
 };
diff --git a/spheres.cpp b/spheres.cpp
index 25c5a67..9e67b87 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -2,7 +2,10 @@
 #include "space_wolff.hpp"
 #include <GL/glut.h>
 
-class animation : public measurement<double, 2, double> {
+const unsigned D = 2;
+typedef Model<double, D, Euclidean<double, D>, double> model;
+
+class animation : public measurement<double, 2, Euclidean<double, D>, double> {
   public:
     animation(double L, unsigned w, int argc, char *argv[]) {
       glutInit(&argc, argv);
@@ -15,7 +18,7 @@ class animation : public measurement<double, 2, double> {
       gluOrtho2D(-1, L + 1, -1 , L + 1);
     }
 
-    void post_cluster(const Model<double, 2, double>& m) override {
+    void post_cluster(const model& m) override {
       glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
       glClear(GL_COLOR_BUFFER_BIT);
       for (const Spin<double, 2, double>& s : m.s) {
@@ -32,6 +35,41 @@ class animation : public measurement<double, 2, double> {
     }
 };
 
+std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs) {
+  return [&mats, &vecs] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+    std::uniform_real_distribution<double> t_dist(0, M.L);
+    std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
+    std::uniform_int_distribution<unsigned> ind_dist(0, M.s.size() - 1);
+    std::uniform_int_distribution<unsigned> coin(0, mats.size() + vecs.size() - 1);
+
+    Vector<double, D> t;
+    Matrix<double, D> m;
+    unsigned flip = coin(rng);
+    if (flip < mats.size()) {
+      for (unsigned j = 0; j < D; j++) {
+        t(j) = (double)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 = vecs[flip - mats.size()];
+    }
+
+    Euclidean<double, D> g(M.L, t, m);
+    return g;
+  };
+
+}
+
 int main(int argc, char* argv[]) {
   const unsigned D = 2;
 
@@ -85,8 +123,11 @@ int main(int argc, char* argv[]) {
       return H * sin(2 * M_PI * 3 * s.x(0) / L);
     };
 
+  std::vector<Matrix<double, D>> mats = torus_mats<double, D>();
+  std::vector<Vector<double, D>> vecs = torus_vecs<double, D>(L);
+  auto g = eGen(mats, vecs);
   animation A(L, 750, argc, argv);
-  Model<double, D, double> sphere(L, Z, B, std::floor(log2(L)), 2, A);
+  model sphere(L, Z, B, g, std::floor(log2(L)), 2, A);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -102,7 +143,7 @@ int main(int argc, char* argv[]) {
   }
 
 
-  sphere.wolff(T, N, rng);
+  sphere.wolff(T, N);
 
   std::ofstream snapfile;
   snapfile.open("sphere_snap.dat");