implemented measurements with autocorrelation-based error tracking

1 files changed, 151 insertions, 16 deletions

diff --git a/src/new_model.cpp b/src/new_model.cpp
index 4393a3d..4f473da 100644
--- a/src/new_model.cpp
+++ b/src/new_model.cpp
@@ -62,19 +62,129 @@ double vsin(const vector_t<2, double>& v1, const vector_t<2, double>& v2) {
 double cos_xy;
 vector_t<2, double> sin_xy;
 
+std::vector<double> ρ(unsigned N, const std::vector<double>& C, double avg) {
+  double C0 = C.front() / N;
+  double avg2 = pow(avg / N, 2);
+
+  std::vector<double> ρtmp;
+
+  for (double Ct : C) {
+    ρtmp.push_back((Ct / N - avg2) / (C0 - avg2));
+  }
+
+  return ρtmp;
+}
+
+double τ(unsigned n, const std::vector<double>& C, double avg) {
+  double τtmp = 0.5;
+  unsigned M = 1;
+  double c = 8.0;
+
+  std::vector<double> ρ_tmp = ρ(n, C, avg);
+
+  while (c * τtmp > M && M < C.size()) {
+    τtmp += ρ_tmp[M];
+    M++;
+  }
+
+  return τtmp;
+}
+
+class quantity {
+  private:
+    double total;
+    std::vector<double> C;
+
+  public:
+    unsigned n;
+    std::list<double> hist;
+    quantity(unsigned lag) : C(lag) {
+      n = 0;
+      total = 0;
+    }
+
+    void add(double x) {
+      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;
+        n++;
+      }
+    }
+
+    double avg() {
+      return total / n;
+    }
+
+    double σ() {
+      return 2.0 / n * τ(n, C, total) * (C[0] / n - pow(this->avg(), 2));
+    }
+
+    double serr() {
+      return sqrt(this->σ());
+    }
+};
+
+class quantity_pair {
+  private:
+    std::vector<double> Cab;
+    double total;
+    double T;
+    quantity a;
+    quantity b;
+
+  public:
+    quantity_pair(unsigned lag, double T) : a(lag), b(lag), Cab(lag), T(T) {};
+
+    void add(double x, double y) {
+      a.add(x);
+      b.add(y);
+
+      if (a.hist.size() == Cab.size()) {
+        unsigned t = 0;
+        auto a_it = a.hist.begin();
+        auto b_it = b.hist.begin();
+        while (a_it != a.hist.end()) {
+          Cab[t] += (*a_it) * (*b_it);
+          t++;
+          a_it++;
+          b_it++;
+        }
+      }
+    }
+
+    double cov() {
+      return 1.0 / a.n * τ(a.n, Cab, sqrt(a.avg() * b.avg()) * a.n) * (Cab[0] / a.n - a.avg() * b.avg());
+    }
+
+    double avg() {
+      return a.avg() - b.avg() / T;
+    }
+
+    double serr() {
+      return sqrt(a.σ() + b.σ() / pow(T, 2) - 2 * this->cov() / T);
+    }
+};
+
 class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i> {
   private:
     double a;
     double K1;
-
+    unsigned w;
     int A;
+    unsigned N;
 
   public:
-    double totalA2;
-    i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a, double K1) : a(a), K1(K1) {
-      A = S.nv;
+    quantity A2;
 
-      totalA2 = 0;
+    i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a, double K1, unsigned nh, unsigned w) : a(a), w(w), K1(K1), A2(nh) {
+      A = S.nv;
+      N = 0;
     }
 
     void ghost_bond_visited(const ising::wolff::system<ising_t, ising_t, graph_i>& S, const graph_i::vertex& v, const ising_t& s_old, const ising_t& s_new, double dE) override {
@@ -90,7 +200,10 @@ class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i
     }
 
     void post_cluster(unsigned, unsigned, const ising::wolff::system<ising_t, ising_t, graph_i>&) override {
-      totalA2 += pow(A, 2);
+      if (N > w) {
+        A2.add(pow(A, 2));
+      }
+      N++;
     }
 };
 
@@ -122,20 +235,24 @@ class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vec
 int main (int argc, char *argv[]) {
 
   // set defaults
-  unsigned N = (unsigned)1e4;
+  unsigned N = (unsigned)1e4; // number of steps between precision checks
+  unsigned w = (unsigned)1e3; // steps to skip before data collection begins
+  unsigned nh = 1e2; // lag to keep in autocorrelation functions
+  double εmag = 1e-1; // precision goal in magnetization
+  double εstiff = 1e-1; // precision goal in spin stiffness
 
   const unsigned D = 2; // we're always in 2D
   unsigned L = 128;     // system size
   double T = 1.0;       // temperature
   double J = 1.0;      // ising nn coupling
   double K1 = 1.0;      // xy nn coupling
-  double K2 = 0.0;     // xy nnn coupling
+  double K2 = 0.0;      // xy nnn coupling
   double α = 1.5;       // ising-xy coupling
 
   int opt;
 
   // take command line arguments
-  while ((opt = getopt(argc, argv, "N:L:T:J:K:a:n:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:J:K:a:n:e:f:")) != -1) {
     switch (opt) {
       case 'N': // number of steps
         N = (unsigned)atof(optarg);
@@ -158,6 +275,12 @@ int main (int argc, char *argv[]) {
       case 'a':
         α = atof(optarg);
         break;
+      case 'e':
+        εmag = atof(optarg);
+        break;
+      case 'f':
+        εstiff = atof(optarg);
+        break;
       default:
         exit(EXIT_FAILURE);
     }
@@ -334,17 +457,29 @@ int main (int argc, char *argv[]) {
   cos_xy = K1 * 2 * pow(L, 2) + K1 * K2 * 4 * pow(L, 2);
   sin_xy.fill(0);
 
-  i_measurement mi(si, α, K1);
+  i_measurement mi(si, α, K1, nh, w);
   x_measurement mxy(sxy, K1, K2, α, &(si.s0));
 
-  double total_cos = 0;
-  double total_sin2 = 0;
+  quantity_pair stiffness(nh, T);
+
+  unsigned nn = 0;
 
-  for (unsigned i = 0; i < N; i++) {
+  while (true) {
     si.run_wolff(1, ising::wolff::gen_ising<graph_i>, mi, rng);
     sxy.run_wolff(1, xy::wolff::generate_rotation_uniform<2, graph_x>, mxy, rng);
-    total_cos += cos_xy;
-    total_sin2 += sin_xy * sin_xy;
+    if (nn > w) {
+      stiffness.add(cos_xy, sin_xy * sin_xy);
+
+      if (nn % N == 0) {
+        double err = mi.A2.serr();
+        double val2 = stiffness.avg();
+        double err2 = stiffness.serr();
+        if (err / mi.A2.avg() <= εmag && err2 / val2 < εstiff) {
+          break;
+        }
+      }
+    }
+    nn++;
   }
 
   std::ifstream checkfile("out.dat");
@@ -358,7 +493,7 @@ int main (int argc, char *argv[]) {
   if (!already_exists) {
     outfile << "N L T J K1 K2 \\[Alpha] M \\[Rho] \\[Sigma]\n";
   }
-  outfile << N << " " << L << " " << T << " " << J << " " << K1 << " " << K2 << " " << α << " " << mi.totalA2 / N / pow(si.nv, 2) << " " << (total_cos - total_sin2 / T) / N / sxy.ne << " " <<(total_cos + total_sin2 / T) / N / sxy.ne << "\n";
+  outfile << N << " " << L << " " << T << " " << J << " " << K1 << " " << K2 << " " << α << " " << mi.A2.avg() / pow(si.nv, 2) << " " << mi.A2.serr() / pow(si.nv, 2) << " " << stiffness.avg() / sxy.ne << " " << stiffness.serr() / sxy.ne << "\n";
   outfile.close();
 }