completely changed how autocorrelation is computed

3 files changed, 155 insertions, 74 deletions

diff --git a/lib/wolff.h b/lib/wolff.h
index 84b3d22..cf24427 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -50,7 +50,16 @@ typedef struct {
   double dc;
 } meas_t;
 
-int32_t sign(double x);
+typedef struct {
+  uint64_t n;
+  uint64_t W;
+  double *OO;
+  double *Op;
+  double O;
+  double O2;
+} autocorr_t;
+
+int8_t sign(double x);
 
 cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost,
                         gsl_rng *r);
@@ -62,5 +71,9 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
 
 void update_meas(meas_t *m, double x);
 
+void update_autocorr(autocorr_t *OO, double O);
+
 double add_to_avg(double mx, double x, uint64_t n);
 
+double rho(autocorr_t *o, uint64_t i);
+
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index bab0908..92e8302 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -2,7 +2,11 @@
 #include "queue.h"
 #include "wolff.h"
 
-int32_t spin_to_sign(bool spin) {
+int8_t spin_to_sign(bool spin) {
+  /* takes a spin (represented by a bool) and converts it to an integer (plus
+   * or minus one). our convention takes true to be negative spin and false to
+   * be positive spin
+   */
   if (spin) {
     return -1;
   } else {
@@ -10,55 +14,63 @@ int32_t spin_to_sign(bool spin) {
   }
 }
 
-int32_t sign(double x) {
-  return x > 0 ? 1 : -1;
+int8_t sign(double x) {
+  // the sign function. zero returns positive sign.
+  return x >= 0 ? 1 : -1;
 }
 
-graph_t *graph_add_ext(const graph_t *g) {
-  graph_t *h = (graph_t *)calloc(1, sizeof(graph_t));
+graph_t *graph_add_ext(const graph_t *G) {
+  /* takes a graph object G and returns tG, the same graph with an extra vertex
+   * that is connected to every other vertex.
+   */
+  graph_t *tG = (graph_t *)calloc(1, sizeof(graph_t));
 
-  h->nv = g->nv + 1;
-  h->ne = g->ne + g->nv;
+  tG->nv = G->nv + 1;
+  tG->ne = G->ne + G->nv;
 
-  h->ev = (uint32_t *)malloc(2 * h->ne * sizeof(uint32_t));
-  h->vei = (uint32_t *)malloc((h->nv + 1) * sizeof(uint32_t));
-  h->ve = (uint32_t *)malloc(2 * h->ne * sizeof(uint32_t));
-  h->vx = (double *)malloc(2 * h->nv * sizeof(double));
-  h->bq = (bool *)malloc(h->nv * sizeof(bool));
+  tG->ev = (uint32_t *)malloc(2 * tG->ne * sizeof(uint32_t));
+  tG->vei = (uint32_t *)malloc((tG->nv + 1) * sizeof(uint32_t));
+  tG->ve = (uint32_t *)malloc(2 * tG->ne * sizeof(uint32_t));
+  tG->vx = (double *)malloc(2 * tG->nv * sizeof(double));
+  tG->bq = (bool *)malloc(tG->nv * sizeof(bool));
 
-  memcpy(h->ev, g->ev, 2 * g->ne * sizeof(uint32_t));
-  memcpy(h->vx, g->vx, 2 * g->nv * sizeof(double));
-  memcpy(h->bq, g->bq, g->nv * sizeof(bool));
-  h->vx[2 * g->nv] = -1;
-  h->vx[2 * g->nv + 1] = -0.5;
-  h->bq[g->nv] = false;
+  memcpy(tG->ev, G->ev, 2 * G->ne * sizeof(uint32_t));
+  memcpy(tG->vx, G->vx, 2 * G->nv * sizeof(double));
+  memcpy(tG->bq, G->bq, G->nv * sizeof(bool));
 
-  for (uint32_t i = 0; i < g->nv; i++) {
-    h->ev[2 * g->ne + 2 * i] = i;
-    h->ev[2 * g->ne + 2 * i + 1] = g->nv;
+  tG->vx[2 * G->nv] = -1;
+  tG->vx[2 * G->nv + 1] = -0.5;
+  tG->bq[G->nv] = false;
+
+  for (uint32_t i = 0; i < G->nv; i++) {
+    tG->ev[2 * G->ne + 2 * i] = i;
+    tG->ev[2 * G->ne + 2 * i + 1] = G->nv;
   }
 
-  for (uint32_t i = 0; i < g->nv; i++) {
-    h->vei[i] = g->vei[i] + i;
+  for (uint32_t i = 0; i < G->nv; i++) {
+    tG->vei[i] = G->vei[i] + i;
 
-    for (uint32_t j = 0; j < g->vei[i + 1] - g->vei[i]; j++) {
-      h->ve[h->vei[i] + j] = g->ve[g->vei[i] + j];
+    for (uint32_t j = 0; j < G->vei[i + 1] - G->vei[i]; j++) {
+      tG->ve[tG->vei[i] + j] = G->ve[G->vei[i] + j];
     }
 
-    h->ve[h->vei[i] + g->vei[i + 1] - g->vei[i]] = g->ne + i;
+    tG->ve[tG->vei[i] + G->vei[i + 1] - G->vei[i]] = G->ne + i;
   }
 
-  h->vei[g->nv] = g->vei[g->nv] + g->nv;
-  h->vei[g->nv + 1] = h->vei[g->nv] + g->nv;
+  tG->vei[G->nv] = G->vei[G->nv] + G->nv;
+  tG->vei[G->nv + 1] = tG->vei[G->nv] + G->nv;
 
-  for (uint32_t i = 0; i < g->nv; i++) {
-    h->ve[h->vei[g->nv] + i] = g->ne + i;
+  for (uint32_t i = 0; i < G->nv; i++) {
+    tG->ve[tG->vei[G->nv] + i] = G->ne + i;
   }
 
-  return h;
+  return tG;
 }
 
 uint32_t get_neighbor(const graph_t *g, uint32_t v, uint32_t i) {
+  // returns the index of the ith neighbor of vertex v on graph g.
+  assert(i < g->vei[v + 1] - g->vei[v]); // don't request a neighbor over the total number of neighbors!
+  
   uint32_t e, v1, v2;
 
   e = g->ve[g->vei[v] + i]; // select the ith bond connected to site
@@ -68,15 +80,23 @@ uint32_t get_neighbor(const graph_t *g, uint32_t v, uint32_t i) {
   return v == v1 ? v2 : v1; // distinguish neighboring site from site itself
 }
 
-cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost,
+cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *s, bool stop_on_ghost,
                         gsl_rng *r) {
+  /* flips a wolff cluster on the graph g with initial state s. s is always
+   * changed by the routine. the probability of adding a normal bond to the
+   * cluster is passed by ps[0], and the probability of adding a bond with the
+   * external field spin to the cluster is passed by ps[1]. if stop_on_ghost is
+   * true, adding the external field spin to the cluster stops execution of the
+   * routine. flip_cluster returns an object of type cluster_t, which encodes
+   * information about the flipped cluster.
+   */
   uint32_t v0;
   int32_t n_h_bonds, n_bonds;
-  bool x0;
+  bool s0;
   cluster_t *c;
 
   v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex
-  x0 = x[v0];                         // record its orientation
+  s0 = s[v0];                         // record its orientation
 
   ll_t *stack = NULL;     // create a new stack
   stack_push(&stack, v0); // push the initial vertex to the stack
@@ -91,8 +111,8 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_o
     v = stack_pop(&stack);
     nn = g->vei[v + 1] - g->vei[v];
 
-    if (x[v] == x0 && !(c->hit_ghost && stop_on_ghost)) { // if the vertex hasn't already been flipped
-      x[v] = !x[v];   // flip the vertex
+    if (s[v] == s0 && !(c->hit_ghost && stop_on_ghost)) { // if the vertex hasn't already been flipped
+      s[v] = !s[v];   // flip the vertex
       if (stop_on_ghost) {
         stack_push(&(c->spins), v);
       }
@@ -110,8 +130,8 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_o
         bond_counter = is_ext ? &(c->dHb) : &(c->dJb);
         prob = is_ext ? ps[1] : ps[0];
 
-        if (x[vn] ==
-            x0) { // if the neighboring site matches the flipping cluster...
+        if (s[vn] ==
+            s0) { // if the neighboring site matches the flipping cluster...
           (*bond_counter)++;
 
           if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
@@ -147,7 +167,7 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
     gsl_rng_set(r, jst_rand_seed());
   }
 
-  if (ps == NULL) {
+  if (ps == NULL) { // computing exponentials is relatively expensive, so calling functions are given the option to supply values calculated once in the entire runtime
     no_ps = true;
     ps = (double *)malloc(2 * sizeof(double));
     ps[0] = 1 - exp(-2 / T);
@@ -185,8 +205,7 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
       }
 
       n_flips = 1;
-                     }
-                     break;
+      } break;
 
     case WOLFF: {
       cluster_t *c = flip_cluster(s->g, ps, s->spins, false, r);
@@ -195,8 +214,7 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
       n_flips = c->nv;
 
       free(c);
-                }
-                break;
+      } break;
     case WOLFF_GHOST: {
       cluster_t *c = flip_cluster(s->g, ps, s->spins, true, r);
 
@@ -204,10 +222,12 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
         while (c->spins != NULL) {
           uint32_t v = stack_pop(&(c->spins));
           s->spins[v] = !s->spins[v];
+          // if we hit the external spin, undo the cluster flip
         }
       } else {
         while (c->spins != NULL) {
           stack_pop(&(c->spins));
+          // we have to clear the memory on the stack anyway...
         }
         s->M += - sign(H) * 2 * c->dHb;
         s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
@@ -216,8 +236,7 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
       n_flips = c->nv;
 
       free(c);
-                      }
-                break;
+      } break;
   }
 
   if (no_ps) {
@@ -235,6 +254,32 @@ double add_to_avg(double mx, double x, uint64_t n) {
   return mx * (n / (n + 1.)) + x * 1. / (n + 1.);
 }
 
+void update_autocorr(autocorr_t *OO, double O) {
+  OO->O = add_to_avg(OO->O, O, OO->n);
+  OO->O2 = add_to_avg(OO->O2, pow(O, 2), OO->n);
+
+  uint64_t lim = OO->W;
+
+  if (OO->n < OO->W) {
+    lim = OO->n;
+  }
+
+  for (uint64_t t = 0; t < lim; t++) {
+    OO->OO[t] = add_to_avg(OO->OO[t], O * OO->Op[t], OO->n - t - 1);
+  }
+
+  for (uint64_t t = 0; t < OO->W - 1; t++) {
+    OO->Op[OO->W - 1 - t] = OO->Op[OO->W - 2 - t];
+  }
+
+  OO->Op[0] = O;
+  OO->n++;
+}
+
+double rho(autocorr_t *o, uint64_t i) {
+  return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2));
+}
+
 void update_meas(meas_t *m, double x) {
   uint64_t n = m->n;
 
diff --git a/src/wolff.c b/src/wolff.c
index 97d1625..40a88e6 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -8,12 +8,13 @@ int main(int argc, char *argv[]) {
   lattice_t lat;
   uint16_t L;
   uint32_t min_runs, lattice_i, sim_i;
-  uint64_t N, n;
+  uint64_t N, n, W;
   double T, H, eps;
 
   sim = WOLFF;
   L = 128;
   N = 100000000000000;
+  W = 10;
   n = 3;
   lat = SQUARE_LATTICE;
   use_dual = false;
@@ -25,11 +26,14 @@ int main(int argc, char *argv[]) {
   output_state = false;
   min_runs = 10;
 
-  while ((opt = getopt(argc, argv, "n:N:L:T:H:m:e:oq:DMas:")) != -1) {
+  while ((opt = getopt(argc, argv, "n:N:L:T:H:m:e:oq:DMas:W:")) != -1) {
     switch (opt) {
     case 'N':
       N = (uint64_t)atof(optarg);
       break;
+    case 'W':
+      W = (uint64_t)atof(optarg);
+      break;
     case 'n':
       n = (uint64_t)atof(optarg);
       break;
@@ -128,7 +132,7 @@ int main(int argc, char *argv[]) {
   uint64_t n_steps = 0;
   double clust_per_sweep = 0;
 
-  meas_t *M, *aM, *eM, *mM, *E, *eE, *mE, *corrE, *corrmE;
+  meas_t *M, *aM, *eM, *mM, *E, *eE, *mE, *clust;
 
   M = calloc(1, sizeof(meas_t));
   aM = calloc(1, sizeof(meas_t));
@@ -137,16 +141,18 @@ int main(int argc, char *argv[]) {
   E = calloc(1, sizeof(meas_t));
   eE = calloc(1, sizeof(meas_t));
   mE = calloc(1, sizeof(meas_t));
+  clust = calloc(1, sizeof(meas_t));
 
-  corrE = (meas_t *)calloc(1, sizeof(meas_t));
-  corrmE = (meas_t *)calloc(1, sizeof(meas_t));
-  double batch_mean_energy;
-  uint64_t batch_size = exp(0.667 * log(N));
-  uint64_t n_batch = 0;
-  uint64_t batch_flips = 0;
+  autocorr_t *autocorr;
+  if (record_autocorrelation) {
+    autocorr = (autocorr_t *)calloc(1, sizeof(autocorr_t));
+    autocorr->W = W;
+    autocorr->OO = (double *)calloc(W, sizeof(double));
+    autocorr->Op = (double *)calloc(W, sizeof(double));
+  }
 
   printf("\n");
-  while (((diff > eps || diff != diff) && n_steps < N) || n_runs < min_runs) {
+  while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) {
     printf("\033[F\033[JWOLFF: sweep %" PRIu64
            ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n",
            n_runs, fabs(E->dx / E->x), M->dx / M->x, E->dc / E->c, M->dc / M->c, clust_per_sweep);
@@ -158,20 +164,11 @@ int main(int argc, char *argv[]) {
       uint32_t tmp_flips = wolff_step(T, H, s, sim, r, bond_probs);
       n_flips += tmp_flips;
       n_clust++;
+      n_steps++;
 
-      n_steps += tmp_flips;
-
-      if (n_runs > 0 && record_autocorrelation) {
-        update_meas(corrE, s->H);
-        if (batch_flips <= batch_size && batch_flips + tmp_flips > batch_size) {
-          update_meas(corrmE, batch_mean_energy / n_batch);
-          batch_mean_energy = 0;
-          n_batch = 0;
-          batch_flips = (int64_t)(batch_flips + tmp_flips) - (int64_t)batch_size;
-        }
-        batch_mean_energy += s->H;
-        n_batch++;
-        batch_flips += tmp_flips;
+      if (record_autocorrelation && n_runs > 0) {
+        update_autocorr(autocorr, s->H);
+        update_meas(clust, tmp_flips);
       }
     }
 
@@ -209,15 +206,42 @@ int main(int argc, char *argv[]) {
 
   FILE *outfile = fopen("out.dat", "a");
 
-  double tau = batch_size * corrmE->c / corrE->c;
+  double tau = 0;
+  double dtau = 0;
+  if (record_autocorrelation) {
+    double kappa = rho(autocorr, W - 1) / rho(autocorr, W - 2);
+    double R = rho(autocorr, W - 1) / (1-kappa);
+
+    double ttau = 0.5 + R;
+    double X = 0.5;
+    double Y = 0;
+    for (uint64_t i = 0; i < W - 1; i++) {
+      ttau += rho(autocorr, i);
+      X += pow(rho(autocorr, i), 2);
+      if (i > 0) {
+        Y += 0.5 * pow(rho(autocorr, i) - rho(autocorr, i-1), 2);
+      } else {
+        Y += 0.5 * pow(rho(autocorr, i) - 1, 2);
+      }
+    }
+
+    double dttau = sqrt(4.0 / autocorr->n * (pow(ttau, 2) * (W +(1+kappa)/(1-kappa)-0.5) + kappa * X / pow(1-kappa, 2) + pow(kappa, 2) * Y / pow(1-kappa, 4)));
+    
+    free(autocorr->OO);
+    free(autocorr->Op);
+    free(autocorr);
+    
+    tau = ttau * clust->x / h->nv;
+    dtau = tau * sqrt(pow(dttau / ttau, 2) + pow(clust->dx / clust->x, 2));
+  }
 
   fprintf(outfile,
-          "%u %.15f %.15f %" PRIu64 " %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n",
+          "%u %.15f %.15f %" PRIu64 " %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n",
           L, T, H, n_runs, n_steps,
           E->x / h->nv, E->dx / h->nv, M->x / h->nv, M->dx / h->nv, E->c / h->nv, E->dc / h->nv, M->c / h->nv, M->dc / h->nv,
           eE->n, eE->x / h->nv, eE->dx / h->nv, eM->x / h->nv, eM->dx / h->nv, eE->c / h->nv, eE->dc / h->nv, eM->c / h->nv, eM->dc / h->nv,
           mE->n, mE->x / h->nv, mE->dx / h->nv, mM->x / h->nv, mM->dx / h->nv, mE->c / h->nv, mE->dc / h->nv, mM->c / h->nv, mM->dc / h->nv,
-          aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv, tau); 
+          aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv, tau, dtau); 
   fclose(outfile);
 
   if (output_state) {
@@ -240,8 +264,7 @@ int main(int argc, char *argv[]) {
   free(E);
   free(eE);
   free(mE);
-  free(corrE);
-  free(corrmE);
+  free(clust);
   free(bond_probs);
   graph_free(h);