made collection of data more standard and efficient

3 files changed, 92 insertions, 74 deletions

diff --git a/lib/wolff.h b/lib/wolff.h
index 428ce85..82ccb9e 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -31,6 +31,17 @@ typedef struct {
   int32_t dHb;
 } cluster_t;
 
+typedef struct {
+  uint64_t n;
+  double x;
+  double dx;
+  double x2;
+  double m2;
+  double m4;
+  double c;
+  double dc;
+} meas_t;
+
 int32_t sign(double x);
 
 cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
@@ -40,3 +51,8 @@ graph_t *graph_add_ext(const graph_t *g);
 
 uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
                     double *ps);
+
+void update_meas(meas_t *m, double x);
+
+double add_to_avg(double mx, double x, uint64_t n);
+
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index 7886bab..59d04fc 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -152,3 +152,26 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
 
   return n_flips;
 }
+
+double add_to_avg(double mx, double x, uint64_t n) {
+  return mx * (n / (n + 1.)) + x * 1. / (n + 1.);
+}
+
+void update_meas(meas_t *m, double x) {
+  uint64_t n = m->n;
+
+  m->x = add_to_avg(m->x, x, n);
+  m->x2 = add_to_avg(m->x2, pow(x, 2), n);
+
+  m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n);
+  m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n);
+
+  if (n > 1) {
+    double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2));
+    m->dx = sqrt(s2 / n);
+    m->c = s2;
+    m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n);
+  }
+
+  (m->n)++;
+}
diff --git a/src/wolff.c b/src/wolff.c
index 3b07a13..4c8d744 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -97,34 +97,23 @@ int main(int argc, char *argv[]) {
 
   double diff = 1e31;
   uint64_t n_runs = 0;
-
-  double E1, E2, dE1, aM1, daM1, M1, M2, dM1, C, dC, aX, daX, X, dX, aMmu2, Mmu2, aMmu4, Mmu4, Emu2, Emu4;
   double clust_per_sweep = 0;
 
-  E1 = 0;
-  E2 = 0;
-  M1 = 0;
-  aM1 = 0;
-  M2 = 0;
-  C = 0;
-  dC = 0;
-  X = 0;
-  aX = 0;
-  dX = 0;
-  dE1 = 0;
-  dM1 = 0;
-  Mmu2 = 0;
-  Mmu4 = 0;
-  aMmu2 = 0;
-  aMmu4 = 0;
-  Emu2 = 0;
-  Emu4 = 0;
+  meas_t *M, *aM, *eM, *mM, *E, *eE, *mE;
+
+  M = calloc(1, sizeof(meas_t));
+  aM = calloc(1, sizeof(meas_t));
+  eM = calloc(1, sizeof(meas_t));
+  mM = calloc(1, sizeof(meas_t));
+  E = calloc(1, sizeof(meas_t));
+  eE = calloc(1, sizeof(meas_t));
+  mE = calloc(1, sizeof(meas_t));
 
   printf("\n");
   while (diff > eps || diff == 0. || 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(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
+           n_runs, fabs(E->dx / E->x), M->dx / M->x, E->dc / E->c, M->dc / M->c, clust_per_sweep);
 
     uint32_t n_flips = 0;
     uint32_t n_clust = 0;
@@ -134,69 +123,52 @@ int main(int argc, char *argv[]) {
       n_clust++;
     }
 
-    E1 = E1 * (n_runs / (n_runs + 1.)) + s->H * 1. / (n_runs + 1.);
-    M1 = M1 * (n_runs / (n_runs + 1.)) + s->M * 1. / (n_runs + 1.);
-    aM1 = aM1 * (n_runs / (n_runs + 1.)) + abs(s->M) * 1. / (n_runs + 1.);
-    E2 = E2 * (n_runs / (n_runs + 1.)) + pow(s->H, 2) * 1. / (n_runs + 1.);
-    M2 = M2 * (n_runs / (n_runs + 1.)) + pow(s->M, 2) * 1. / (n_runs + 1.);
-
-    Mmu2 = Mmu2 * (n_runs / (n_runs + 1.)) +
-           pow(s->M - M1, 2) * 1. / (n_runs + 1.);
-    Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) +
-           pow(s->M - M1, 4) * 1. / (n_runs + 1.);
-    aMmu2 = aMmu2 * (n_runs / (n_runs + 1.)) +
-           pow(abs(s->M) - aM1, 2) * 1. / (n_runs + 1.);
-    aMmu4 = aMmu4 * (n_runs / (n_runs + 1.)) +
-           pow(abs(s->M) - aM1, 4) * 1. / (n_runs + 1.);
-    Emu2 = Emu2 * (n_runs / (n_runs + 1.)) +
-           pow(s->H - E1, 2) * 1. / (n_runs + 1.);
-    Emu4 = Emu4 * (n_runs / (n_runs + 1.)) +
-           pow(s->H - E1, 4) * 1. / (n_runs + 1.);
-
-    if (n_runs > 1) {
-      double Msigma2 = n_runs / (n_runs - 1) * (M2 - pow(M1, 2));
-      double aMsigma2 = n_runs / (n_runs - 1) * (M2 - pow(aM1, 2));
-      X = Msigma2 / T;
-      aX = aMsigma2 / T;
-      dX =
-          sqrt((Mmu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Mmu2, 2)) / n_runs) /
-          T;
-      daX =
-          sqrt((aMmu4 - (n_runs - 3.) / (n_runs - 1.) * pow(aMmu2, 2)) / n_runs) /
-          T;
-      
-
-      double Esigma2 = n_runs / (n_runs - 1) * (E2 - pow(E1, 2));
-      C = Esigma2 / T;
-      dC =
-          sqrt((Emu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Emu2, 2)) / n_runs) /
-          T;
-
-      dE1 = sqrt(Esigma2 / n_runs);
-      dM1 = sqrt(Msigma2 / n_runs);
-      daM1 = sqrt(aMsigma2 / n_runs);
-
-      if (M_stop) {
-        diff = fabs(dM1 / M1);
-      } else {
-        diff = fabs(dX / X);
-      }
+    double ss = 1;
+    if (s->spins[s->g->nv - 1]) {
+      ss = -1;
+    }
+
+    double HH = 1;
+    if (H < 0) {
+      HH = -1;
+    }
+
+    update_meas(M, (double)(s->M));
+    update_meas(aM, HH * fabs((double)(s->M)));
+    update_meas(E, s->H);
+
+    if (HH * s->M * ss > 0) {
+      update_meas(eM, HH * fabs((double)(s->M)));
+      update_meas(eE, s->H);
+    } else {
+      update_meas(mM, - HH * fabs((double)(s->M)));
+      update_meas(mE, s->H);
+    }
+
+    if (M_stop) {
+      diff = fabs(M->dx / M->x);
+    } else {
+      diff = fabs(M->dc / M->c);
     }
 
-    clust_per_sweep = clust_per_sweep * (n_runs / (n_runs + 1.)) +
-                      (n_clust * 1. / N) * 1. / (n_runs + 1.);
+    clust_per_sweep = add_to_avg(clust_per_sweep, n_clust * 1. / N, n_runs);
 
     n_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(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
+         n_runs, fabs(E->dx / E->x), M->dx / M->x, E->dc / E->c, M->dc / M->c, clust_per_sweep);
 
   FILE *outfile = fopen("out.dat", "a");
+
   fprintf(outfile,
-          "%u %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L,
-          T, H, E1 / h->nv, dE1 / h->nv, M1 / h->nv, dM1 / h->nv, C / h->nv,
-          dC / h->nv, X / h->nv, dX / h->nv, aM1 / h->nv, daM1 / h->nv, aX / h->nv, daX / h->nv);
+          "%u %.15f %.15f %" 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\n", L,
+          T, H, n_runs,
+          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);
   fclose(outfile);
 
   if (output_state) {
@@ -212,6 +184,13 @@ int main(int argc, char *argv[]) {
   graph_free(s->g);
   free(s->spins);
   free(s);
+  free(M);
+  free(aM);
+  free(eM);
+  free(mM);
+  free(E);
+  free(eE);
+  free(mE);
   free(bond_probs);
   graph_free(h);