Added gitignore and fixed Ising to work with new paradigm

1 files changed, 151 insertions, 99 deletions

diff --git a/ising.cpp b/ising.cpp
index 1987a8f..38dce3b 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -1,14 +1,53 @@
 
+#include "smiley.hpp"
 #include "space_wolff.hpp"
+#include "torus_symmetries.hpp"
+#include <GL/glut.h>
+
+const unsigned D = 2;
+typedef Model<signed, D, TorusGroup<signed, D>, signed> model;
+
+class animation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
+private:
+  uint64_t t1;
+  uint64_t t2;
+  unsigned n;
+  unsigned tmp;
+
+public:
+  animation(double L, unsigned w, int argc, char* argv[]) {
+    t1 = 0;
+    t2 = 0;
+    n = 0;
+
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+    glutInitWindowSize(w, w);
+    glutCreateWindow("wolff");
+    glClearColor(0.0, 0.0, 0.0, 0.0);
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    gluOrtho2D(0, L, 0, L);
+  }
 
-std::function<double(spin<signed, 2, signed>)> B_sin(unsigned L, unsigned n, double H) {
-  return [n, H, L] (spin<signed, 2, signed> s) -> double {
-    return H * s.s * cos(2 * M_PI * n * s.x(0) / ((double)L));
-  };
-}
+  void post_cluster(const model& m) override {
+    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
+    glClear(GL_COLOR_BUFFER_BIT);
+    glClearColor(0.0, 1.0, 0.0, 0.0);
+    for (const Spin<signed, 2, signed>* s : m.s) {
+      if (s->s == 1) {
+        glColor3f(0.0, 0.0, 0.0);
+      } else if (s->s == -1) {
+        glColor3f(1.0, 1.0, 1.0);
+      }
+      Vector<signed, 2> xx = m.s0.inverse().act(*s).x;
+      glRecti(xx(0), xx(1), xx(0) + 1, xx(1) + 1);
+    }
+    glFlush();
+  }
+};
 
 int main(int argc, char* argv[]) {
-  const unsigned D = 2;
 
   unsigned L = 32;
   unsigned N = 1000;
@@ -23,120 +62,133 @@ int main(int argc, char* argv[]) {
 
   while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:n:")) != -1) {
     switch (opt) {
-      case 'N': 
-        N = (unsigned)atof(optarg);
-        break;
-      case 'L':
-        L = atoi(optarg);
-        break;
-      case 'T':
-        T = atof(optarg);
-        break;
-      case 'H':
-        H = atof(optarg);
-        break;
-      case 'e':
-        ε = atof(optarg);
-        break;
-      case 'm':
-        mod = atoi(optarg);
-        break;
-      case 'M':
-        multi = atoi(optarg);
-        break;
-      case 'n':
-        mag = atof(optarg);
-        break;
-      default:
-        exit(1);
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'L':
+      L = atoi(optarg);
+      break;
+    case 'T':
+      T = atof(optarg);
+      break;
+    case 'H':
+      H = atof(optarg);
+      break;
+    case 'e':
+      ε = atof(optarg);
+      break;
+    case 'm':
+      mod = atoi(optarg);
+      break;
+    case 'M':
+      multi = atoi(optarg);
+      break;
+    case 'n':
+      mag = atof(optarg);
+      break;
+    default:
+      exit(1);
     }
   }
 
   double pZ = 1.0 - exp(-1.0 / T);
 
-  std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z =
-    [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double {
-      bool old_one_one = false;
-      bool old_many_ones = false;
-      bool old_any_two = false;
-
-      vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x);
-      
-      for (unsigned i = 0; i < D; i++) {
-        if (old_diff(i) == 1 && !old_one_one) {
-          old_one_one = true;
-        } else if (old_diff(i) == 1 && old_one_one) {
-          old_many_ones = true;
-        } else if (old_diff(i) > 1) {
-          old_any_two = true;
-        }
+  std::function<double(const Spin<signed, D, signed>&, const Spin<signed, D, signed>&)> Z =
+      [L, pZ](const Spin<signed, D, signed>& s1, const Spin<signed, D, signed>& s2) -> double {
+    bool old_one_one = false;
+    bool old_many_ones = false;
+    bool old_any_two = false;
+
+    Vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x);
+
+    for (unsigned i = 0; i < D; i++) {
+      if (old_diff(i) == 1 && !old_one_one) {
+        old_one_one = true;
+      } else if (old_diff(i) == 1 && old_one_one) {
+        old_many_ones = true;
+      } else if (old_diff(i) > 1) {
+        old_any_two = true;
       }
+    }
 
-      bool were_on_someone = !old_one_one && !old_any_two;
-      bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;
-
-      if (were_on_someone) {
-        return -std::numeric_limits<double>::infinity();;
-      } else if (were_nearest_neighbors) {
-        return s1.s * s2.s;
-      } else {
-        return 0.0;
-      }
-    };
-
-  std::function<double(spin<signed, D, signed>)> B_face =
-    [L, H] (spin<signed, D, signed> s) -> double {
-      return H * s.s * smiley[s.x(0) * 16 / L][s.x(1) * 16 / L];
-    };
-
-  std::function<double(spin<signed, D, signed>)> B;
+    bool were_on_someone = !old_one_one && !old_any_two;
+    bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;
 
-  if (mod > 0) {
-    B = B_sin(L, mod, H);
-  } else {
-    B = B_face;
-  }
+    if (were_on_someone) {
+      return -std::numeric_limits<double>::infinity();
+      ;
+    } else if (were_nearest_neighbors) {
+      return s1.s * s2.s;
+    } else {
+      return 0.0;
+    }
+  };
 
-  std::function<std::set<unsigned>(model<signed, D, signed>&, unsigned, spin<signed, D, signed>)> neighbors =
-    [] (model<signed, D, signed>& m, unsigned i0, spin<signed, D, signed> s1) -> std::set<unsigned> {
-      std::set<unsigned> nn;
-      if (i0 < m.s.size()) {
-        std::set<unsigned> nn0 = m.dict.neighbors(m.s[i0].x, 1);
-        std::set<unsigned> nn1 = m.dict.neighbors(s1.x, 1);
-        nn.insert(nn0.begin(), nn0.end());
-        nn.insert(nn1.begin(), nn1.end());
-        nn.insert(m.s.size());
-      } else {
-        for (unsigned i = 0; i < m.s.size(); i++) {
-          nn.insert(i);
-        }
-      }
-      return nn;
-    };
+  std::function<double(const Spin<signed, D, signed>&)> B_face =
+      [L, H](const Spin<signed, D, signed>& s) -> double {
+    return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L];
+  };
 
-  model<signed, D, signed> ising(L, L, Z, B, neighbors);
+  std::function<double(Spin<signed, D, signed>)> B = B_face;
 
-  randutils::auto_seed_128 seeds;
-  std::mt19937 rng{seeds};
+  Model<signed, D, TorusGroup<signed, D>, signed> ising(L, Z, B);
 
-  std::uniform_int_distribution<unsigned> coin(0, 1);
+  Rng rng;
 
   unsigned n = 0;
   unsigned up = 0;
   unsigned down = 0;
   for (unsigned i = 0; i < L; i++) {
     for (unsigned j = 0; j < L; j++) {
-      if ((coin(rng) && up < pow(L, 2) * mag) || down >= pow(L, 2) * mag) {
-        ising.s.push_back({{i, j}, 1});
+      Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
+      ss->x = {i, j};
+      if (rng.uniform<double>(0, 1) < mag) {
+        ss->s = 1;
         up++;
       } else {
-        ising.s.push_back({{i, j}, -1});
+        ss->s = -1;
         down++;
       }
-      ising.dict.record({i, j}, n);
+      ising.s.push_back(ss);
+      ising.dict.insert(ss);
       n++;
     }
   }
+
+  std::vector<Vector<signed, 2>> torusVectors = torus_vecs<signed, 2>(L);
+  std::vector<Matrix<signed, 2>> torusMatrices = torus_mats<signed, 2>();
+
+  Gen<signed, D, TorusGroup<signed, D>, signed> g =
+      [L, &torusVectors,
+       &torusMatrices](Model<signed, D, TorusGroup<signed, D>, signed>& M,
+                       Rng& r) -> Transformation<signed, D, TorusGroup<signed, D>, signed>* {
+    Matrix<signed, 2> m;
+    Vector<signed, 2> t;
+
+    m = r.pick(torusMatrices);
+    t(0) = r.uniform<signed>(0, L);
+    t(1) = r.uniform<signed>(0, L);
+    t = t - m * t;
+
+    TorusGroup<signed, 2> g = TorusGroup<signed, 2>({(signed)L, t, m});
+
+    Spin<signed, 2, signed>* ss = r.pick(M.s);
+    Spin<signed, 2, signed> s2 = g.act(*ss);
+
+    while (ss->x(0) == s2.x(0) && ss->x(1) == s2.x(1)) {
+      ss = r.pick(M.s);
+      s2 = g.act(*ss);
+    }
+
+    std::set<Spin<signed, 2, signed>*> s2s = M.dict.at(s2.x);
+
+    return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
+  };
+
+  animation A(L, 750, argc, argv);
+
+  ising.wolff(T, {g}, A, N);
+
   /*
   for (unsigned i = 0; i < L; i++) {
     for (unsigned j = 0; j < L; j++) {
@@ -151,8 +203,7 @@ int main(int argc, char* argv[]) {
   }
   */
 
-  ising.update_energy();
-
+  /*
   while (true) {
     ising.wolff(T, N, rng);
     std::array<double, 2> τ = ising.Eq.τ();
@@ -168,8 +219,8 @@ int main(int argc, char* argv[]) {
   std::array<double, 2> act = ising.Eq.τ();
   std::vector<double> ρ = ising.Eq.ρ();
 
-  outfile << L << " " << T << " " << mod << " " << H << " " << ising.Eq.num_added() << " " << ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1];
-  for (double ρi : ρ) {
+  outfile << L << " " << T << " " << mod << " " << H << " " << ising.Eq.num_added() << " " <<
+  ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1]; for (double ρi : ρ) {
     outfile << " " << ρi;
   }
   outfile << "\n";
@@ -195,6 +246,7 @@ int main(int argc, char* argv[]) {
   }
 
   snapfile.close();
+  */
 
   return 0;
 }