added new method for computing autocorrelation

4 files changed, 159 insertions, 19 deletions

diff --git a/lib/convex.c b/lib/convex.c
new file mode 100644
index 0000000..33af1b4
--- /dev/null
+++ b/lib/convex.c
@@ -0,0 +1,93 @@
+
+#include "convex.h"
+
+double slope(point_t *P, point_t *Q) {
+  return (Q->y - P->y) / ((double)(Q->x) - (double)(P->x));
+}
+
+double *get_convex_minorant(uint64_t n, double *Gammas) {
+  list_t *L = (list_t *)calloc(1, sizeof(list_t));
+  L->p = (point_t *)calloc(1, sizeof(point_t));
+  L->p->x = 0;
+  L->p->y = Gammas[0];
+
+  list_t *pos = L;
+
+  for (uint64_t i = 1; i < n; i++) {
+    pos->next = (list_t *)calloc(1, sizeof(list_t));
+    pos->next->p = (point_t *)calloc(1, sizeof(point_t));
+    pos->next->p->x = i;
+    pos->next->p->y = Gammas[i];
+    pos->next->prev = pos;
+    pos = pos->next;
+  }
+
+  pos->next = (list_t *)calloc(1, sizeof(list_t));
+  pos->next->p = (point_t *)calloc(1, sizeof(point_t));
+  pos->next->p->x = n;
+  pos->next->p->y = 0;
+  pos->next->prev = pos;
+
+  list_t *X = L;
+  list_t *Y = L->next;
+  list_t *Z = Y->next;
+
+  while (true) {
+    if (slope(X->p, Y->p) <= slope(Y->p, Z->p)) {
+      X = Y;
+      Y = Z;
+      if (Z->next == NULL) {
+        break;
+      } else {
+        Z = Z->next;
+      }
+    } else {
+      Y->prev->next = Y->next;
+      Y->next->prev = Y->prev;
+      free(Y->p);
+      free(Y);
+      if (X->prev != NULL) {
+        Y = X;
+        X = X->prev;
+      } else {
+        if (Z->next != NULL) {
+          Y = Z;
+          Z = Z->next;
+        } else {
+          break;
+        }
+      }
+    }
+  }
+
+  pos = L;
+
+  double *g = (double *)calloc(n, sizeof(double));
+  double rho = 0;
+
+  for (uint64_t i = 0; i < n; i++) {
+    if (i > pos->next->p->x) {
+      pos = pos->next;
+    }
+
+    g[i] = pos->p->y + ((double)i - (double)(pos->p->x)) *  (pos->next->p->y - pos->p->y) / ((double)(pos->next->p->x) - (double)(pos->p->x));
+
+    if (Gammas[i] - g[i] > rho) {
+      rho = Gammas[i] - g[i];
+    }
+  }
+
+  for (uint64_t i = 0; i < n; i++) {
+    g[i] += rho / 2;
+  }
+
+  while (L != NULL) {
+    free(L->p);
+    list_t *L_save = L;
+    L = L->next;
+    free(L);
+  }
+
+  return g;
+}
+
diff --git a/lib/convex.h b/lib/convex.h
new file mode 100644
index 0000000..29f4ae0
--- /dev/null
+++ b/lib/convex.h
@@ -0,0 +1,21 @@
+
+#pragma once
+
+#include <inttypes.h>
+#include <stdbool.h>
+#include <stdlib.h>
+#include <string.h>
+
+typedef struct {
+  uint64_t x;
+  double y;
+} point_t;
+
+typedef struct list_tag {
+  struct list_tag *prev;
+  struct list_tag *next;
+  point_t *p;
+} list_t;
+
+double *get_convex_minorant(uint64_t n, double *Gammas);
+
diff --git a/lib/wolff.h b/lib/wolff.h
index 66ce493..dfe7b1c 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -17,6 +17,7 @@
 #include <jst/rand.h>
 
 #include "queue.h"
+#include "convex.h"
 
 typedef enum {
   WOLFF,
diff --git a/src/wolff.c b/src/wolff.c
index 1e62c8c..a806681 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -136,7 +136,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, *clust;
+  meas_t *M, *aM, *eM, *mM, *E, *eE, *mE, *clust, *everyE, *blockE;
 
   M = calloc(1, sizeof(meas_t));
   aM = calloc(1, sizeof(meas_t));
@@ -146,12 +146,18 @@ int main(int argc, char *argv[]) {
   eE = calloc(1, sizeof(meas_t));
   mE = calloc(1, sizeof(meas_t));
   clust = calloc(1, sizeof(meas_t));
+  everyE = calloc(1, sizeof(meas_t));
+
+  blockE = calloc(1,sizeof(meas_t));
+
+  uint64_t blocksize = 1000;
+  double Etot = 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->W = 2 * W + 1;
+    autocorr->OO = (double *)calloc(2 * W + 1, sizeof(double));
   }
 
   printf("\n");
@@ -167,9 +173,17 @@ 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++;
+      if (n_runs > 0){
+        n_steps++;
+      }
 
       if (record_autocorrelation && n_runs > 0) {
+        update_meas(everyE, s->H);
+        if (n_steps % blocksize == 0) {
+          update_meas(blockE, Etot / blocksize);
+          Etot = 0;
+        }
+        Etot += s->H;
         if (n_steps % ac_skip == 0) {
           update_autocorr(autocorr, s->H);
         }
@@ -213,25 +227,36 @@ int main(int argc, char *argv[]) {
 
   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 *Gammas = (double *)malloc((W + 1) * sizeof(double));
+
+    Gammas[0] = 1 + rho(autocorr, 0);
+    printf("%g ", Gammas[0]);
+    for (uint64_t i = 0; i < W; i++) {
+      Gammas[1 + i] = rho(autocorr, 2 * i + 1) + rho(autocorr, 2 * i + 2);
+    } 
+
+    uint64_t n;
+    for (n = 0; n < W + 1; n++) {
+      if (Gammas[n] <= 0) {
+        break;
       }
     }
 
-    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)));
+    if (n == W) {
+      printf("WARNING: correlation function never hit the noise floor.\n");
+    }
+
+    double *conv_Gamma = get_convex_minorant(n, Gammas);
+
+    double ttau = - 0.5;
+
+    for (uint64_t i = 0; i < n; i++) {
+      ttau += conv_Gamma[i];
+    }
     
+    free(Gammas);
     free(autocorr->OO);
     while (autocorr->Op != NULL) {
       stack_pop_d(&(autocorr->Op));
@@ -239,7 +264,7 @@ int main(int argc, char *argv[]) {
     free(autocorr);
     
     tau = ttau * ac_skip * clust->x / h->nv;
-    dtau = tau * sqrt(pow(dttau / ttau, 2) + pow(clust->dx / clust->x, 2));
+    dtau = 0;
   }
 
   fprintf(outfile,