Fixed the periodic sphere code to work with the new backend and implemented saving cluster flip information

3 files changed, 148 insertions, 105 deletions

diff --git a/euclidean.hpp b/euclidean.hpp
index 49ec2d9..b1ce2c6 100644
--- a/euclidean.hpp
+++ b/euclidean.hpp
@@ -95,6 +95,16 @@ public:
     return s_new;
   }
 
+  Vector<T, D> act(const Vector<T, D>& s) const {
+    Vector<T, D> s_new = t + r * s;
+
+    for (unsigned i = 0; i < D; i++) {
+      s_new(i) = fmod(L + s_new(i), L);
+    }
+
+    return s_new;
+  }
+
   TorusGroup act(const TorusGroup& x) const {
     Vector<T, D> tnew = r * x.t + t;
     Matrix<T, D> rnew = r * x.r;
diff --git a/spheres.cpp b/spheres.cpp
index b7a9d94..4235d66 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -15,9 +15,13 @@ private:
   uint64_t t2;
   unsigned n;
   unsigned tmp;
+  unsigned wait;
+  double L;
+  TorusGroup<double, D> s0_tmp;
+  std::ofstream& file;
 
 public:
-  animation(double L, unsigned w, int argc, char* argv[]) {
+  animation(double L, unsigned w, int argc, char* argv[], unsigned wait, std::ofstream& file) : s0_tmp(0), wait(wait), L(50 * L), file(file) {
     t1 = 0;
     t2 = 0;
     n = 0;
@@ -25,37 +29,80 @@ public:
     glutInit(&argc, argv);
     glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
     glutInitWindowSize(w, w);
-    glutCreateWindow("wolff");
+    glutCreateWindow("wolffWindow");
     glClearColor(0.0, 0.0, 0.0, 0.0);
     glMatrixMode(GL_PROJECTION);
     glLoadIdentity();
-    gluOrtho2D(-1, L + 1, -1, L + 1);
+    gluOrtho2D(-1, L+1, -1, L+1);
   }
 
-  void pre_cluster(const model&, unsigned, const TorusGroup<double, D>&) override { tmp = 0; }
-
-  void plain_site_transformed(const model&, const Spin<double, D, double>*,
-                              const Spin<double, D, double>&) override {
-    tmp++;
-  }
-
-  void post_cluster(const model& m) override {
+  void pre_cluster(const model& m, unsigned, const Transformation<double, D, TorusGroup<double, D>, double>* t) override {
+    s0_tmp = m.s0.inverse();
+    tmp = 0;
     glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
     glClear(GL_COLOR_BUFFER_BIT);
-    for (const Spin<double, 2, double>& s : m.s) {
+    glColor3d(1.0, 0.0, 0.0);
+    glBegin(GL_LINES);
+      Vector<double, 2> r = t->r.t / 2;
+      double θ = atan2(r(1), r(0));
+      Vector<double, 2> v1 = s0_tmp.act({r(0) + L * sin(θ), r(1) - L * cos(θ)});
+      Vector<double, 2> v2 = s0_tmp.act({r(0) - L * sin(θ), r(1) + L * cos(θ)});
+      glVertex2d(v1(0), v1(1));
+      glVertex2d(v2(0), v2(1));
+      if (n > wait)
+        file << v1.transpose() << " " << v2.transpose() << std::endl;
+    glEnd();
+    for (const Spin<double, 2, double>* s : m.s) {
       glBegin(GL_POLYGON);
       unsigned n_points = 50;
       glColor3f(0.0f, 0.0f, 0.0f);
+      if (n > wait)
+        file << s << " " << s->s << " " << s0_tmp.act(*s).x.transpose() << " ";
       for (unsigned i = 0; i < n_points; i++) {
-        glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points),
-                   m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points));
+        glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
+                   s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
       }
       glEnd();
     }
+    if (n > wait)
+      file << std::endl;
+    glLineWidth(3);
     glFlush();
+  }
+
+  void plain_site_transformed(const model& m, const Transformation<double, D, TorusGroup<double, D>, double>& t) override {
+    if (n > wait) {
+      if (t.current().size() > 1) {
+        glColor3f(0.0f, 0.0f, 1.0f);
+      } else {
+        glColor3f(0.0f, 1.0f, 0.0f);
+      }
+      for (const Spin<double, 2, double>* s : t.current()) {
+        if (s != NULL) {
+          file << s << " " << s0_tmp.act(t.r.act(*s)).x.transpose() << " ";
+          glBegin(GL_POLYGON);
+          unsigned n_points = 50;
+          for (unsigned i = 0; i < n_points; i++) {
+            glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
+                       s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
+          }
+          glEnd();
+        } else {
+          file << s << " " << s0_tmp.act({0, 0}).transpose() << " ";
+        }
+      }
+    }
+    tmp++;
+  }
 
+  void post_cluster(const model& m) override {
     t1 += tmp;
     t2 += tmp * tmp;
+    if (n > wait) {
+      file << std::endl;
+      glFlush();
+      sleep(2);
+    }
     n++;
   }
 
@@ -68,40 +115,29 @@ public:
   double var() { return (t2 - t1 * t1 / (double)n) / (double)n; }
 };
 
-std::function<TorusGroup<double, D>(const model&, randutils::mt19937_rng&)>
-eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs,
-     double ε, double L) {
-  return
-      [&mats, &vecs, L, ε](const model& M, randutils::mt19937_rng& rng) -> TorusGroup<double, 2> {
-        Vector<double, D> t;
-        Matrix<double, D> m;
-        unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1));
-        if (flip < mats.size()) {
-          unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size());
-          t = M.s[f_ind].x;
-          for (unsigned j = 0; j < D; j++) {
-            t(j) = fmod(10 * L + t(j) + rng.variate<double, std::normal_distribution>(0.0, ε), L);
-          }
-          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;
-              }
-            }
-          }
+Gen<double, D, TorusGroup<double, D>, double> eGen(double L) {
+  std::vector<Vector<double, 2>> torusVectors = torus_vecs<double, 2>(L);
+  std::vector<Matrix<double, 2>> torusMatrices = torus_mats<double, 2>();
+  return [L, torusVectors,
+          torusMatrices](Model<double, D, TorusGroup<double, D>, double>& M,
+                         Rng& r) -> Transformation<double, D, TorusGroup<double, D>, double>* {
+    Matrix<double, 2> m;
+    Vector<double, 2> t;
 
-          t = vecs[flip - mats.size()];
-        }
+    m = r.pick(torusMatrices);
+    t(0) = r.uniform<double>(0, L);
+    t(1) = r.uniform<double>(0, L);
+    t = t - m * t;
+
+    TorusGroup<double, 2> g = TorusGroup<double, 2>({(double)L, t, m});
 
-        TorusGroup<double, D> g(M.L, t, m);
-        return g;
-      };
+    Spin<double, 2, double>* ss = r.pick(M.s);
+
+    return new SpinFlip<double, 2, TorusGroup<double, 2>, double>(M, g, ss);
+  };
 }
 
+
 int main(int argc, char* argv[]) {
   const unsigned D = 2;
 
@@ -135,8 +171,8 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  double k = 1e2;
-  double a = 0.05;
+  double k = 1e8;
+  double a = 0.0;
 
   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 {
@@ -158,47 +194,26 @@ int main(int argc, char* argv[]) {
     return H * s.x(1);
   };
 
-  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, L, L);
-  animation A(L, 750, argc, argv);
+  auto g = eGen(L);
+  std::ofstream ofile("test.dat");
+  animation A(L, 750, argc, argv, 1000, ofile);
   model sphere(L, Z, B);
 
   randutils::mt19937_rng rng;
 
-  sphere.s.reserve(n);
+  sphere.s.resize(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, 0.5});
-    sphere.dict.insert(&sphere.s.back());
-  }
-
-  sphere.wolff(T, g, A, N);
-  std::ofstream outfile;
-  outfile.open("test.dat");
-
-  for (signed i = -10; i <= 10; i++) {
-    A.clear();
-    double ε = pow(2, -4 + i / 2.0);
-    auto gn = eGen(mats, vecs, ε, L);
-    sphere.wolff(T, gn, A, N);
-    outfile << ε << " " << A.var() / sphere.s.size() << std::endl;
-    std::cout << ε << " " << A.var() / sphere.s.size() << std::endl;
-  }
-
-  outfile.close();
-
-  std::ofstream snapfile;
-  snapfile.open("sphere_snap.dat");
-
-  for (Spin<double, D, double> s : sphere.s) {
-    Spin<double, D, double> rs = sphere.s0.inverse().act(s);
-    snapfile << rs.x.transpose() << "\n";
+  for (unsigned i = 0; i < sphere.s.size(); i++) {
+    Spin<double, 2, double>* ss = new Spin<double, 2, double>();
+    ss->x = {(i / nx) * L / nx, (i % nx) * L / nx};
+    ss->s = rng.pick({0.45, 0.45});
+    sphere.s[i] = ss;
+    sphere.dict.insert(ss);
   }
 
-  snapfile.close();
+  sphere.wolff(T, {g}, A, N);
+  ofile.close();
 
   return 0;
 }
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 3f8747b..56abe77 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -36,28 +36,6 @@ public:
     gluOrtho2D(-L, L, -L, L);
   }
 
-  void plain_site_transformed(const model& m, const Transformation<double, D, Euclidean<double, D>, double>& t) override {
-    if (n > wait) {
-    if (t.current().size() > 1) {
-      glColor3f(0.0f, 0.0f, 1.0f);
-    } else {
-      glColor3f(0.0f, 1.0f, 0.0f);
-    }
-    for (const Spin<double, 2, double>* s : t.current()) {
-      if (s != NULL) {
-      glBegin(GL_POLYGON);
-      unsigned n_points = 50;
-      for (unsigned i = 0; i < n_points; i++) {
-        glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
-                   s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
-      }
-      glEnd();
-      }
-    }
-    }
-    tmp++;
-  }
-
   void pre_cluster(const model& m, unsigned, const Transformation<double, D, Euclidean<double, D>, double>* t) override {
     s0_tmp = m.s0.inverse();
     tmp = 0;
@@ -71,26 +49,58 @@ public:
       Vector<double, 2> v2 = s0_tmp.act({r(0) - L * sin(θ), r(1) + L * cos(θ)});
       glVertex2d(v1(0), v1(1));
       glVertex2d(v2(0), v2(1));
+      if (n > wait)
+        file << v1.transpose() << " " << v2.transpose() << std::endl;
     glEnd();
     for (const Spin<double, 2, double>* s : m.s) {
       glBegin(GL_POLYGON);
       unsigned n_points = 50;
       glColor3f(0.0f, 0.0f, 0.0f);
+      if (n > wait)
+        file << s << " " << s->s << " " << s0_tmp.act(*s).x.transpose() << " ";
       for (unsigned i = 0; i < n_points; i++) {
         glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
                    s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
       }
       glEnd();
     }
+    if (n > wait)
+      file << std::endl;
     glLineWidth(3);
     glFlush();
   }
 
+  void plain_site_transformed(const model& m, const Transformation<double, D, Euclidean<double, D>, double>& t) override {
+    if (n > wait) {
+      if (t.current().size() > 1) {
+        glColor3f(0.0f, 0.0f, 1.0f);
+      } else {
+        glColor3f(0.0f, 1.0f, 0.0f);
+      }
+      for (const Spin<double, 2, double>* s : t.current()) {
+        if (s != NULL) {
+          file << s << " " << s0_tmp.act(t.r.act(*s)).x.transpose() << " ";
+          glBegin(GL_POLYGON);
+          unsigned n_points = 50;
+          for (unsigned i = 0; i < n_points; i++) {
+            glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
+                       s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
+          }
+          glEnd();
+        } else {
+          file << s << " " << s0_tmp.act({0, 0}).transpose() << " ";
+        }
+      }
+    }
+    tmp++;
+  }
+
   void post_cluster(const model& m) override {
-    glFlush();
     t1 += tmp;
     t2 += tmp * tmp;
     if (n > wait) {
+      file << std::endl;
+      glFlush();
       sleep(2);
     }
     n++;
@@ -280,10 +290,18 @@ int main(int argc, char* argv[]) {
     return H * s.x.norm();
   };
 
+  std::function<double(Spin<double, D, double>)> B_hard = [L, H](Spin<double, D, double> s) -> double {
+    if (fabs(s.x(0)) < 3 * L / 4 && fabs(s.x(1)) < 3 * L / 4) {
+      return 0;
+    } else {
+      return std::numeric_limits<double>::infinity();
+    }
+  };
+
   auto g1 = eGen(1);
   auto g2 = mGen(0.1);
   auto g3 = tGen(0.1);
-  auto g4 = rGen(1);
+  auto g4 = rGen(0);
 
   auto tag = std::chrono::high_resolution_clock::now();
 
@@ -303,8 +321,8 @@ int main(int argc, char* argv[]) {
   unsigned nx = floor(sqrt(n));
   for (unsigned i = 0; i < sphere.s.size(); i++) {
     Spin<double, 2, double>* ss = new Spin<double, 2, double>();
-    ss->x = {(i / nx) * L / nx, (i % nx) * L / nx};
-    ss->s = rng.pick({0.5, 0.5});
+    ss->x = {(i / nx) * L / nx - L / 2, (i % nx) * L / nx - L / 2};
+    ss->s = rng.pick({0.45, 0.45});
     sphere.s[i] = ss;
     sphere.dict.insert(ss);
   }