Split up Ising-related functions and added an executable for measuring Ising energy autocorrelation time.

6 files changed, 348 insertions, 307 deletions

diff --git a/ising.cpp b/ising.cpp
deleted file mode 100644
index 838ec9e..0000000
--- a/ising.cpp
+++ /dev/null
@@ -1,280 +0,0 @@
-
-#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;
-  bool color;
-
-public:
-  animation(double L, unsigned w, bool tcolor, int argc, char* argv[]) {
-    color = tcolor;
-    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);
-  }
-
-  void post_cluster(const model& m) override {
-    glClearColor(1.0, 1.0, 1.0, 1.0);
-    glClear(GL_COLOR_BUFFER_BIT);
-    for (const Spin<signed, 2, signed>* s : m.s) {
-      if (s->s == 1) {
-        if (color)
-          glColor3f(1.0, 0.0, 0.0);
-        else
-          glColor3f(1.0, 1.0, 1.0);
-      } else if (s->s == -1) {
-        if (color)
-          glColor3f(0.0, 0.0, 1.0);
-        else 
-          glColor3f(0.0, 0.0, 0.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[]) {
-
-  unsigned L = 32;
-  unsigned N = 1000;
-  unsigned mod = 0;
-  unsigned multi = 1e4;
-  double mag = 0.5;
-  double pop = 1.0;
-  double T = 2.0 / log(1.0 + sqrt(2.0));
-  double H = 1.0;
-  double ε = 0.1;
-  bool color = false;
-
-  int opt;
-
-  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:r:p:c")) != -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 'r':
-      mag = atof(optarg);
-      break;
-    case 'p':
-      pop = atof(optarg);
-      break;
-    case 'c':
-      color = true;
-      break;
-    default:
-      exit(1);
-    }
-  }
-
-  double pZ = 1.0 - exp(-1.0 / T);
-
-  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;
-
-  if (mod == 0) {
-    B = [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];
-    };
-  } else {
-    B = [mod, L, H](const Spin<signed, D, signed>& s) -> double {
-      return H * s.s * sin(s.x(0) * mod * 2 * M_PI / L);
-    };
-  }
-
-
-  Model<signed, D, TorusGroup<signed, D>, signed> ising(L, Z, B);
-
-  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 (rng.uniform<double>(0, 1) < pop) {
-      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 {
-        ss->s = -1;
-        down++;
-      }
-      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);
-
-    if (s2s.empty()) {
-      return new SpinFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss);
-    } else {
-      return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
-    }
-  };
-
-  animation A(L, 750, color, argc, argv);
-
-  ising.wolff(T, {g}, A, N);
-
-  /*
-  for (unsigned i = 0; i < L; i++) {
-    for (unsigned j = 0; j < L; j++) {
-      if (i < L / 2) {
-        ising.s.push_back({{i, j}, 1});
-      } else {
-        ising.s.push_back({{i, j}, -1});
-      }
-      ising.dict.record<signed>({i, j}, n);
-      n++;
-    }
-  }
-  */
-
-  /*
-  while (true) {
-    ising.wolff(T, N, rng);
-    std::array<double, 2> τ = ising.Eq.τ();
-    std::cout << ising.Eq.num_added() << " " << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
-    if (τ[1] / τ[0] < ε && τ[0] * multi < ising.Eq.num_added()) {
-      break;
-    }
-  }
-
-  std::ofstream outfile;
-  outfile.open("out.dat", std::ios::app);
-
-  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 << " " << ρi;
-  }
-  outfile << "\n";
-
-  std::ofstream snapfile;
-  snapfile.open("snap.dat");
-
-  std::vector<std::vector<signed>> snap(L);
-  for (std::vector<signed>& line : snap) {
-    line.resize(L);
-  }
-
-  for (spin<signed, D, signed> s : ising.s) {
-    spin<signed, D, signed> snew = ising.s0.inverse().act(s);
-    snap[snew.x(0)][snew.x(1)] = snew.s;
-  }
-
-  for (std::vector<signed> row : snap) {
-    for (signed s : row) {
-      snapfile << s << " ";
-    }
-    snapfile << "\n";
-  }
-
-  snapfile.close();
-  */
-
-  return 0;
-}
-
diff --git a/ising.hpp b/ising.hpp
new file mode 100644
index 0000000..6f140be
--- /dev/null
+++ b/ising.hpp
@@ -0,0 +1,106 @@
+
+#include "smiley.hpp"
+#include "space_wolff.hpp"
+#include "torus_symmetries.hpp"
+
+const unsigned D = 2;
+typedef Model<signed, D, TorusGroup<signed, D>, signed> isingModel;
+typedef Spin<signed, D, signed> isingSpin;
+
+std::function<double(const isingSpin&, const isingSpin&)> isingZ(unsigned L) {
+  return [L](const isingSpin& s1, const isingSpin& 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(isingSpin)> isingBFace(unsigned L, double H) {
+  return [L, H](const isingSpin& s) -> double {
+    return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L];
+  };
+}
+
+std::function<double(isingSpin)> isingBMod(unsigned L, unsigned mod, double H) {
+  return [mod, L, H](const isingSpin& s) -> double {
+    return H * s.s * sin(s.x(0) * mod * 2 * M_PI / L);
+  };
+}
+
+Gen<signed, D, TorusGroup<signed, D>, signed> isingGen(unsigned L) {
+  std::vector<Vector<signed, 2>> torusVectors = torus_vecs<signed, 2>(L);
+  std::vector<Matrix<signed, 2>> torusMatrices = torus_mats<signed, 2>();
+  return [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);
+
+    if (s2s.empty()) {
+      return new SpinFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss);
+    } else {
+      return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
+    }
+  };
+}
+
+void isingPopulate(isingModel& m, unsigned L, double pop, double mag) {
+  Rng rng;
+
+  for (unsigned i = 0; i < L; i++) {
+    for (unsigned j = 0; j < L; j++) {
+      if (rng.uniform<double>(0, 1) < pop) {
+      Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
+      ss->x = {i, j};
+      if (rng.uniform<double>(0, 1) < mag) {
+        ss->s = 1;
+      } else {
+        ss->s = -1;
+      }
+      m.s.push_back(ss);
+      m.dict.insert(ss);
+      }
+    }
+  }
+}
diff --git a/ising_animate.cpp b/ising_animate.cpp
new file mode 100644
index 0000000..e6586d9
--- /dev/null
+++ b/ising_animate.cpp
@@ -0,0 +1,105 @@
+#include "ising.hpp"
+#include <GL/glut.h>
+
+class Animation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
+private:
+  bool color;
+
+public:
+  Animation(double L, unsigned w, bool tcolor, int argc, char* argv[]) {
+    color = tcolor;
+
+    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);
+  }
+
+  void post_cluster(const isingModel& m) override {
+    glClearColor(1.0, 1.0, 1.0, 1.0);
+    glClear(GL_COLOR_BUFFER_BIT);
+    for (const Spin<signed, 2, signed>* s : m.s) {
+      if (s->s == 1) {
+        if (color)
+          glColor3f(1.0, 0.0, 0.0);
+        else
+          glColor3f(1.0, 1.0, 1.0);
+      } else if (s->s == -1) {
+        if (color)
+          glColor3f(0.0, 0.0, 1.0);
+        else
+          glColor3f(0.0, 0.0, 0.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[]) {
+  unsigned L = 32;
+  unsigned N = 1000;
+  unsigned mod = 0;
+  double mag = 0.5;
+  double pop = 1.0;
+  double T = 2.0 / log(1.0 + sqrt(2.0));
+  double H = 1.0;
+  bool color = false;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:L:T:H:m:r:p:c")) != -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 'm':
+      mod = atoi(optarg);
+      break;
+    case 'r':
+      mag = atof(optarg);
+      break;
+    case 'p':
+      pop = atof(optarg);
+      break;
+    case 'c':
+      color = true;
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  std::function<double(Spin<signed, D, signed>)> B;
+
+  if (mod == 0) {
+    B = isingBFace(L, H);
+  } else {
+    B = isingBMod(L, mod, H);
+  }
+
+  isingModel ising(L, isingZ(L), B);
+  isingPopulate(ising, L, pop, mag);
+
+  auto g = isingGen(L);
+
+  Animation A(L, 750, color, argc, argv);
+
+  ising.wolff(T, {g}, A, N);
+
+  return 0;
+}
diff --git a/ising_autocorrelation.cpp b/ising_autocorrelation.cpp
new file mode 100644
index 0000000..053be01
--- /dev/null
+++ b/ising_autocorrelation.cpp
@@ -0,0 +1,115 @@
+#include <iostream>
+#include <fstream>
+
+#include "ising.hpp"
+#include "quantity.hpp"
+
+class Autocorrelation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
+public:
+  Quantity E;
+  Autocorrelation(unsigned lag) : E(lag) {}
+
+  void post_cluster(const isingModel& m) override {
+    double E_tmp = 0;
+
+    TorusGroup<signed, D> s0Inv = m.s0.inverse();
+
+    for (const isingSpin* si : m.s) {
+      E_tmp -= m.B(s0Inv.act(*si));
+
+      for (const isingSpin* sj : m.dict.neighbors(si->x)) {
+        if (si != sj) {
+          E_tmp -= m.Z(*si, *sj) / 2;
+        }
+      }
+    }
+
+    E.add(E_tmp);
+  }
+};
+
+int main(int argc, char* argv[]) {
+  unsigned L = 32;
+  unsigned N = 1000;
+  unsigned mod = 0;
+  unsigned multi = 1e4;
+  double mag = 0.5;
+  double pop = 1.0;
+  double T = 2.0 / log(1.0 + sqrt(2.0));
+  double H = 1.0;
+  double ε = 0.1;
+  unsigned lag = 1e2;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:r:p:l:")) != -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 'r':
+      mag = atof(optarg);
+      break;
+    case 'p':
+      pop = atof(optarg);
+      break;
+    case 'l':
+      lag = (unsigned)atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  isingModel ising(L, isingZ(L), isingBMod(L, mod, H));
+  isingPopulate(ising, L, pop, mag);
+
+  auto g = isingGen(L);
+
+  measurement<signed, D, TorusGroup<signed, D>, signed> A_empty;
+
+  ising.wolff(T, {g}, A_empty, N);
+
+  Autocorrelation A(lag);
+
+  while (true) {
+    ising.wolff(T, {g}, A, N);
+    std::array<double, 2> τ = A.E.τ();
+    std::cout << A.E.num_added() << " " << A.E.avg() << " " << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
+    if (τ[1] / τ[0] < ε && τ[0] * multi < A.E.num_added()) {
+      break;
+    }
+  }
+
+  std::ofstream outfile;
+  outfile.open("out.dat", std::ios::app);
+
+  std::array<double, 2> act = A.E.τ();
+  std::vector<double> ρ = A.E.ρ();
+
+  outfile << L << " " << T << " " << mod << " " << H << " " << A.E.num_added() << " " <<
+  A.E.avg() << " " << A.E.serr() << " " << act[0] << " " << act[1]; for (double ρi : ρ) {
+    outfile << " " << ρi;
+  }
+  outfile << "\n";
+
+  return 0;
+}
diff --git a/octree.hpp b/octree.hpp
index bb13242..d8f0aa3 100644
--- a/octree.hpp
+++ b/octree.hpp
@@ -5,8 +5,6 @@
 #include <set>
 #include <unordered_map>
 
-#include <iostream>
-
 #include "spin.hpp"
 #include "vector.hpp"
 
diff --git a/quantity.hpp b/quantity.hpp
index 55d5a28..d2a7dea 100644
--- a/quantity.hpp
+++ b/quantity.hpp
@@ -1,57 +1,54 @@
 
 #pragma once
 
-#include <vector>
-#include <list>
-#include <cmath>
 #include <array>
+#include <cmath>
+#include <list>
+#include <vector>
 
 class Quantity {
 private:
   double total;
   double total2;
   std::vector<double> C;
-  unsigned wait;
 
 public:
   unsigned n;
   std::list<double> hist;
 
-  Quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
+  Quantity(unsigned lag) : C(lag) {
     n = 0;
     total = 0;
     total2 = 0;
   }
 
   void add(double x) {
-    if (n > wait) {
-      hist.push_front(x);
-      if (hist.size() > C.size()) {
-        hist.pop_back();
-        unsigned t = 0;
-        for (double a : hist) {
-          C[t] += a * x;
-          t++;
-        }
-        total += x;
-        total2 += pow(x, 2);
+    hist.push_front(x);
+    if (hist.size() > C.size()) {
+      hist.pop_back();
+      unsigned t = 0;
+      for (double a : hist) {
+        C[t] += a * x;
+        t++;
       }
+      total += x;
+      total2 += pow(x, 2);
+      n++;
     }
-    n++;
   }
 
-  double avg() const { return total / (n - wait); }
+  double avg() const { return total / n; }
 
-  double avg2() const { return total2 / (n - wait); }
+  double avg2() const { return total2 / n; }
 
   std::vector<double> ρ() const {
-    double C0 = C.front() / (n - wait);
-    double avg2 = pow(total / (n - wait), 2);
+    double C0 = C.front() / n;
+    double avg2 = pow(total / n, 2);
 
     std::vector<double> ρtmp;
 
     for (double Ct : C) {
-      ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
+      ρtmp.push_back((Ct / n - avg2) / (C0 - avg2));
     }
 
     return ρtmp;
@@ -69,15 +66,15 @@ public:
       M++;
     }
 
-    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
   }
 
   double σ() const {
-    return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+    return 2.0 / n * this->τ()[0] * (C[0] / n - pow(this->avg(), 2));
   }
 
   double serr() const { return sqrt(this->σ()); }
 
-  unsigned num_added() const { return n - wait; }
+  unsigned num_added() const { return n; }
 };