refactoring and generalization in preparation for pressure ensemble

1 files changed, 45 insertions, 4 deletions

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");