added measurement tools

2 files changed, 142 insertions, 70 deletions

diff --git a/ising.cpp b/ising.cpp
index 6939749..645c759 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -8,10 +8,11 @@ int main(int argc, char* argv[]) {
   unsigned N = 1000;
   double T = 2.0 / log(1.0 + sqrt(2.0));
   double H = 1.0;
+  double ε = 0.1;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:L:T:H:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:H:e:")) != -1) {
     switch (opt) {
       case 'N': 
         N = (unsigned)atof(optarg);
@@ -25,6 +26,9 @@ int main(int argc, char* argv[]) {
       case 'H':
         H = atof(optarg);
         break;
+      case 'e':
+        ε = atof(optarg);
+        break;
       default:
         exit(1);
     }
@@ -67,12 +71,10 @@ int main(int argc, char* argv[]) {
     [] (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> os1 = m.dict.on_site(s1.x);
-        std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x);
-        std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x);
+        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(os1.begin(), os1.end());
         nn.insert(m.s.size());
       } else {
         for (unsigned i = 0; i < m.s.size(); i++) {
@@ -119,7 +121,16 @@ int main(int argc, char* argv[]) {
   }
   */
 
-  ising.wolff(T, N, rng);
+  ising.update_energy();
+
+  while (true) {
+    ising.wolff(T, N, rng);
+    std::array<double, 2> τ = ising.Eq.τ();
+    std::cout << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
+    if (τ[1] / τ[0] < ε) {
+      break;
+    }
+  }
 
   std::vector<signed> output(pow(L, D));
 
@@ -128,13 +139,21 @@ int main(int argc, char* argv[]) {
     output[L * rs.x(1) + rs.x(0)] = s.s;
   }
 
+  std::ofstream outfile;
+  outfile.open("snap.dat");
+
   for (unsigned i = 0; i < L; i++) {
     for (unsigned j = 0; j < L; j++) {
       unsigned out = output[L * i + j] == 1 ? 1 : 0;
-      std::cout << out;
+      outfile << out << " ";
     }
-    std::cout << "\n";
+    outfile << "\n";
   }
+  outfile.close();
+
+  std::array<double, 2> act = ising.Eq.τ();
+
+  std::cout << act[0] << " " << act[1] << "\n";
 
   return 0;
 }
diff --git a/space_wolff.hpp b/space_wolff.hpp
index aaef673..b60e274 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -3,6 +3,7 @@
 #include <list>
 #include <set>
 #include <iostream>
+#include <fstream>
 #include <functional>
 #include <random>
 #include <queue>
@@ -66,7 +67,7 @@ class euclidean {
     }
 
     template <class state>
-    spin<T, D, state> act(spin<T, D, state> s) {
+    spin<T, D, state> act(const spin<T, D, state>& s) const {
       spin<T, D, state> s_new;
 
       s_new.x = t + r * s.x;
@@ -79,7 +80,7 @@ class euclidean {
       return s_new;
     }
 
-    euclidean act(const euclidean& x) {
+    euclidean act(const euclidean& x) const {
       vector<T, D> tnew = r * x.t + t;
       matrix<T, D> rnew = r * x.r;
 
@@ -88,7 +89,7 @@ class euclidean {
       return pnew;
     }
 
-    euclidean inverse() {
+    euclidean inverse() const {
       vector<T, D> tnew = - r.transpose() * t;
       matrix<T, D> rnew = r.transpose();
 
@@ -98,6 +99,7 @@ class euclidean {
     }
 };
 
+// to-do: clean up nnn... code
 template <int D>
 class dictionary {
   private:
@@ -108,7 +110,7 @@ class dictionary {
     dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {};
 
     template <class T>
-    unsigned dictionary_index(vector<T, D> x) {
+    unsigned dictionary_index(vector<T, D> x) const {
       unsigned pos_ind = 0;
 
       for (unsigned i = 0; i < D; i++) {
@@ -129,42 +131,25 @@ class dictionary {
     };
 
     template <class T>
-    std::set<unsigned> on_site(vector<T, D> x) {
-      return d[dictionary_index<T>(x)];
+    std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const {
+      return nearest_neighbors_of(dictionary_index<T>(x), depth, {});
     };
 
-    template <class T>
-    std::set<unsigned> nearest_neighbors(vector<T, D> x) {
-      unsigned ind = dictionary_index<T>(x);
-      std::set<unsigned> ns;
+    std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, std::list<unsigned> ignores) const {
+      std::set<unsigned> ns = d[ind];
 
-      for (unsigned i = 0; i < D; i++) {
-        for (signed j : {-1, 1}) {
-          unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));
-          for (unsigned nii : d[ni]) {
-            ns.insert(nii);
-          }
-        }
-      }
-
-      return ns;
-    };
+      if (depth > 0) {
+        for (unsigned i = 0; i < D; i++) {
+          if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k){return i==k;})) {
+            std::list<unsigned> ignores_new = ignores;
+            ignores_new.push_back(i);
+            for (signed j : {-1, 1}) {
+              unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));
 
-    template <class T>
-    std::set<unsigned> next_nearest_neighbors(vector<T, D> x) {
-      unsigned ind = dictionary_index<T>(x);
-      std::set<unsigned> ns;
+              std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new);
 
-      for (unsigned i = 0; i < D; i++) {
-        for (signed j : {-1, 1}) {
-          unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));
-          for (unsigned k = 0; k < D; k++) {
-            if (k != i) {
-              for (signed l : {-1, 1}) {
-                unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1));
-                for (unsigned nnii : d[nni]) {
-                  ns.insert(nnii);
-                }
+              for (unsigned guy : nns) {
+                ns.insert(guy);
               }
             }
           }
@@ -173,37 +158,85 @@ class dictionary {
 
       return ns;
     };
+};
 
-    template <class T>
-    std::set<unsigned> next_next_nearest_neighbors(vector<T, D> x) {
-      unsigned ind = dictionary_index<T>(x);
-      std::set<unsigned> ns;
+class quantity {
+  private:
+    double total;
+    double total2;
+    std::vector<double> C;
+    unsigned wait;
 
-      for (unsigned i = 0; i < D; i++) {
-        for (signed j : {-1, 1}) {
-          unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));
-          for (unsigned k = 0; k < D; k++) {
-            if (k != i) {
-              for (signed l : {-1, 1}) {
-                unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1));
-                for (unsigned m = 0; m < D; m++) {
-                  if (m != i && m != k) {
-                    for (signed n : {-1, 1}) {
-                      unsigned nnni = pow(L, m + 1) * (nni / ((unsigned)pow(L, m + 1))) + fmod(pow(L, m + 1) + nni + n * pow(L, m), pow(L, m + 1));
-                      for (unsigned nnnii : d[nnni]) {
-                        ns.insert(nnnii);
-                      }
-                    }
-                  }
-                }
-              }
-            }
+  public:
+    unsigned n;
+    std::list<double> hist;
+
+    quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
+      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);
         }
       }
+      n++;
+    }
 
-      return ns;
-    };
+    double avg() const {
+      return total / (n - wait);
+    }
+
+    double avg2() const {
+      return total2 / (n - wait);
+    }
+
+    std::vector<double> ρ() const {
+      double C0 = C.front() / (n - wait);
+      double avg2 = pow(total / (n - wait), 2);
+
+      std::vector<double> ρtmp;
+
+      for (double Ct : C) {
+        ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
+      }
+
+      return ρtmp;
+    }
+
+    std::array<double, 2> τ() const {
+      double τtmp = 0.5;
+      unsigned M = 1;
+      double c = 8.0;
+
+      std::vector<double> ρ_tmp = this->ρ();
+
+      while (c * τtmp > M && M < C.size()) {
+        τtmp += ρ_tmp[M];
+        M++;
+      }
+
+      return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+    }
+
+    double σ() const {
+      return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+    }
+
+    double serr() const {
+      return sqrt(this->σ());
+    }
 };
 
 template <class U, int D, class state>
@@ -216,15 +249,29 @@ class model {
     std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors;
     std::function<double(spin<U, D, state>, spin<U, D, state>)> Z;
     std::function<double(spin<U, D, state>)> B;
-    double E;
+    long double E;
+    quantity Eq;
+    quantity Cq;
 
     model(unsigned L, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z,
         std::function<double(spin<U, D, state>)> B,
         std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) : 
-      L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B) {
+      L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) {
         }
 
+    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]);
+        }
+
+        E -= B(s0.inverse().act(s[i]));
+      }
+    }
+
     void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) {
+      unsigned cluster_size = 0;
       std::uniform_real_distribution<double> dist(0.0, 1.0);
 
       std::queue<unsigned> queue;
@@ -288,9 +335,12 @@ class model {
             dict.remove(s[i].x, i);
             s[i] = si_new;
             dict.record(s[i].x, i);
+            cluster_size++;
           }
         }
       }
+
+      Cq.add(cluster_size);
     }
 
     void wolff(double T, unsigned N, std::mt19937& rng) {
@@ -327,6 +377,9 @@ class model {
         euclidean<U, D> g(L, t, m);
 
         this->step(T, ind_dist(rng), g, rng);
+
+        this->update_energy();
+        Eq.add(E);
       }
     }
 };