some optimizations

8 files changed, 100 insertions, 68 deletions

diff --git a/lib/cluster.c b/lib/cluster.c
index fa417fa..7274eb9 100644
--- a/lib/cluster.c
+++ b/lib/cluster.c
@@ -193,8 +193,13 @@ v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) {
 
 //    if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
     if (!marks[v]) {
+      bool v_is_external = false;
       double *s_old, *s_new, *R_tmp; 
 
+      if (v == s->g->nv - 1) {
+        v_is_external = true;
+      }
+
       //tree_insert(&T, v);
       marks[v] = true;
 
@@ -209,17 +214,24 @@ v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) {
 
       for (v_t i = 0; i < nn; i++) {
         v_t vn = s->g->v_adj[s->g->v_i[v] + i];
+
+        bool vn_is_external = false;
+
+        if (vn == s->g->nv - 1) {
+          vn_is_external = true;
+        }
+
         double *sn;
-        if (vn != s->g->nv - 1) {
+
+        if (!vn_is_external) {
           sn = &(s->spins[s->n * vn]);
         }
-        double prob;
 
-        bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
+        double prob;
 
-        if (is_ext) {
+        if (v_is_external || vn_is_external) {
           double *rs_old, *rs_new;
-          if (vn == s->g->nv - 1) {
+          if (vn_is_external) {
             rs_old = vector_rotate_inverse(s->n, s->R, s_old);
             rs_new = vector_rotate_inverse(s->n, s->R, s_new);
           } else {
@@ -237,8 +249,6 @@ v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) {
         } else {
           double dE = (s->J)(vector_dot(s->n, sn, s_old)) - (s->J)(vector_dot(s->n, sn, s_new));
           prob = 1.0 - exp(-dE / s->T);
-          //printf("(%g %g) (%g %g) (%g %g) %g\n", s_old[0], s_old[1], s_new[0], s_new[1], sn[0], sn[1], dE);
-          //getchar();
           s->E += dE;
         }
 
diff --git a/lib/measurement.c b/lib/measurement.c
index 435c749..ad824f6 100644
--- a/lib/measurement.c
+++ b/lib/measurement.c
@@ -2,10 +2,10 @@
 #include "measurement.h"
 
 double add_to_avg(double mx, double x, count_t n) {
-  return mx * (n / (n + 1.0)) + x * 1.0 / (n + 1.0);
+  return mx * (n / (n + 1.0)) + x / (n + 1.0);
 }
 
-void update_meas(meas_t *m, double x) {
+void meas_update(meas_t *m, double x) {
   count_t n = m->n;
 
   m->x = add_to_avg(m->x, x, n);
@@ -14,16 +14,30 @@ void update_meas(meas_t *m, double x) {
   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)++;
 }
 
+double meas_dx(const meas_t *m) {
+  return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2)));
+}
+
+double meas_c(const meas_t *m) {
+  return m->n / (m->n - 1.) * (m->x2 - pow(m->x, 2));
+}
+
+double meas_dc(const meas_t *m) {
+  return sqrt((m->m4 - (m->n - 3.)/(m->n - 1.) * pow(m->m2, 2)) / m->n);
+}
+
 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);
@@ -56,7 +70,7 @@ void update_autocorr(autocorr_t *OO, double O) {
   OO->n++;
 }
 
-double rho(autocorr_t *o, count_t i) {
+double rho(const autocorr_t *o, count_t i) {
   return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2));
 }
 
diff --git a/lib/measurement.h b/lib/measurement.h
index 4ced4dc..a24df9f 100644
--- a/lib/measurement.h
+++ b/lib/measurement.h
@@ -8,12 +8,9 @@
 typedef struct {
   uint64_t n;
   double x;
-  double dx;
   double x2;
   double m2;
   double m4;
-  double c;
-  double dc;
 } meas_t;
 
 typedef struct {
@@ -25,9 +22,15 @@ typedef struct {
   double O2;
 } autocorr_t;
 
-void update_meas(meas_t *m, double x);
+void meas_update(meas_t *m, double x);
+
+double meas_dx(const meas_t *m);
+
+double meas_c(const meas_t *m);
+
+double meas_dc(const meas_t *m);
 
 void update_autocorr(autocorr_t *OO, double O);
 
-double rho(autocorr_t *o, uint64_t i);
+double rho(const autocorr_t *o, uint64_t i);
 
diff --git a/lib/orthogonal.c b/lib/orthogonal.c
index 7e0a71b..87569ae 100644
--- a/lib/orthogonal.c
+++ b/lib/orthogonal.c
@@ -2,24 +2,28 @@
 #include "orthogonal.h"
 
 void vector_replace(q_t n, double *v1, const double *v2) {
+  // writes vector v2 of length n to memory located at v1
   for (q_t i = 0; i < n; i++) {
     v1[i] = v2[i];
   }
 }
 
 void vector_add(q_t n, double *v1, const double *v2) {
+  // adds vector v2 of length n to vector v1
   for (q_t i = 0; i < n; i++) {
     v1[i] += v2[i];
   }
 }
 
 void vector_subtract(q_t n, double *v1, const double *v2) {
+  // subtracts vector v2 of length n from vector v1
   for (q_t i = 0; i < n; i++) {
     v1[i] -= v2[i];
   }
 }
 
-double *vector_rotate(q_t n, double *rot, double *vec) {
+double *vector_rotate(q_t n, const double *rot, const double *vec) {
+  // multiplies n by n rotation matrix rot to vector vec
   double *rot_vec = (double *)malloc(n * sizeof(double));
 
   double prod = 0.0;
@@ -46,7 +50,7 @@ double *vector_rotate_inverse(q_t n, const double *rot, const double *vec) {
   return rot_vec;
 }
 
-double vector_dot(q_t n, double *v1, double *v2) {
+double vector_dot(q_t n, const double *v1, const double *v2) {
   double dot = 0;
 
   for (q_t i = 0; i < n; i++) {
@@ -56,7 +60,7 @@ double vector_dot(q_t n, double *v1, double *v2) {
   return dot;
 }
 
-double *orthogonal_rotate(q_t n, double *r, double *m) {
+double *orthogonal_rotate(q_t n, const double *r, const double *m) {
   double *mul = (double *)calloc(n * n, sizeof(double));
 
   for (q_t i = 0; i < n; i++) {
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index a763b08..60d5f49 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -12,13 +12,13 @@ void vector_add(q_t n, double *v1, const double *v2);
 
 void vector_subtract(q_t n, double *v1, const double *v2);
 
-double *vector_rotate(q_t n, double *rot, double *vec);
+double *vector_rotate(q_t n, const double *rot, const double *vec);
 
 double *vector_rotate_inverse(q_t n, const double *rot, const double *vec);
 
-double vector_dot(q_t n, double *v1, double *v2);
+double vector_dot(q_t n, const double *v1, const double *v2);
 
-double *orthogonal_rotate(q_t n, double *m1, double *m2);
+double *orthogonal_rotate(q_t n, const double *m1, const double *m2);
 
 double *gen_rot(gsl_rng *r, q_t n);
 
diff --git a/src/wolff_dgm.c b/src/wolff_dgm.c
index a9287f1..f11b296 100644
--- a/src/wolff_dgm.c
+++ b/src/wolff_dgm.c
@@ -116,7 +116,7 @@ int main(int argc, char *argv[]) {
   while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) {
     if (!silent) 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, h->nv / clust->x);
+           n_runs, fabs(meas_dx(E) / E->x), meas_dx(M) / M->x, meas_dc(E) / meas_c(E), meas_dc(M) / meas_c(M), h->nv / clust->x);
 
     count_t n_flips = 0;
 
@@ -129,7 +129,7 @@ int main(int argc, char *argv[]) {
 
       if (n_runs > 0) {
         n_steps++;
-        update_meas(clust, tmp_flips);
+        meas_update(clust, tmp_flips);
       }
 
       if (record_autocorrelation && n_runs > 0) {
@@ -139,7 +139,7 @@ int main(int argc, char *argv[]) {
       }
     }
 
-    update_meas(M, s->M);
+    meas_update(M, s->M);
     h_t min_h, max_h;
     min_h = MAX_H;
     max_h = MIN_H;
@@ -150,10 +150,10 @@ int main(int argc, char *argv[]) {
         max_h = s->spins[i];
       }
     }
-    update_meas(dM, max_h - min_h);
-    update_meas(E, s->E);
+    meas_update(dM, max_h - min_h);
+    meas_update(E, s->E);
 
-    diff = fabs(E->dc / E->c);
+    diff = fabs(meas_dc(E) / meas_c(E));
 
     n_runs++;
   }
@@ -162,7 +162,7 @@ int main(int argc, char *argv[]) {
   }
   printf("WOLFF: 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, h->nv / clust->x);
+         n_runs, fabs(meas_dx(E) / E->x), meas_dx(M) / M->x, meas_dc(E) / meas_c(E), meas_dc(M) / meas_c(M), h->nv / clust->x);
 
   double tau = 0;
   bool tau_failed = false;
@@ -217,11 +217,11 @@ int main(int argc, char *argv[]) {
   FILE *outfile = fopen("out.m", "a");
 
   fprintf(outfile, "<|D->%" PRID ",L->%" PRIL ",T->%.15f", D, L, T);
-  fprintf(outfile, ",E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->%.15f", E->x / h->nv, E->dx / h->nv, E->c / h->nv, E->dc / h->nv, M->x / h->nv);
-  fprintf(outfile, ",\\[Delta]M->%.15f", M->dx / h->nv);
-  fprintf(outfile, ",\\[Chi]->%.15f", M->c / h->nv);
-  fprintf(outfile, ",\\[Delta]\\[Chi]->%.15f", M->dc / h->nv);
-  fprintf(outfile, ",w->%.15f,\\[Delta]w->%.15f,wc->%.15f,\\[Delta]wc->%.15f,Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", dM->x, dM->dx, dM->c, dM->dc, clust->x / h->nv, clust->dx / h->nv, clust->c / h->nv, clust->dc / h->nv,tau);
+  fprintf(outfile, ",E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->%.15f", E->x / h->nv, meas_dx(E) / h->nv, meas_c(E) / h->nv, meas_dc(E) / h->nv, M->x / h->nv);
+  fprintf(outfile, ",\\[Delta]M->%.15f", meas_dx(M) / h->nv);
+  fprintf(outfile, ",\\[Chi]->%.15f", meas_c(M) / h->nv);
+  fprintf(outfile, ",\\[Delta]\\[Chi]->%.15f", meas_dc(M) / h->nv);
+  fprintf(outfile, ",w->%.15f,\\[Delta]w->%.15f,wc->%.15f,\\[Delta]wc->%.15f,Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", dM->x, meas_dx(dM), meas_c(dM), meas_dc(dM), clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau);
 
   fclose(outfile);
 
diff --git a/src/wolff_potts.c b/src/wolff_potts.c
index eea9ed7..1d2d6bf 100644
--- a/src/wolff_potts.c
+++ b/src/wolff_potts.c
@@ -180,7 +180,7 @@ int main(int argc, char *argv[]) {
   while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) {
     if (!silent) 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[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x);
+           n_runs, fabs(meas_dx(E) / E->x), meas_dx(M[0]) / M[0]->x, meas_dc(E) / meas_c(E), meas_dc(M[0]) / meas_c(M[0]), h->nv / clust->x);
 
     count_t n_flips = 0;
 
@@ -220,21 +220,21 @@ int main(int argc, char *argv[]) {
           lifetime_n++;
         } else {
           if (cur_M != MAX_Q) {
-            update_meas(lifetimes[cur_M], lifetime_n);
+            meas_update(lifetimes[cur_M], lifetime_n);
           }
           lifetime_n = 0;
           cur_M = max_M_i;
         }
       } else {
         if (cur_M != MAX_Q) {
-          update_meas(lifetimes[cur_M], lifetime_n);
+          meas_update(lifetimes[cur_M], lifetime_n);
           cur_M = MAX_Q;
         }
       }
 
       if (n_runs > 0) {
         n_steps++;
-        update_meas(clust, (double)tmp_flips);
+        meas_update(clust, (double)tmp_flips);
       }
 
       if (record_autocorrelation && n_runs > 0) {
@@ -245,19 +245,19 @@ int main(int argc, char *argv[]) {
     }
 
     for (q_t i = 0; i < q; i++) {
-      update_meas(M[i], s->M[i]);
+      meas_update(M[i], s->M[i]);
     }
-    update_meas(E, s->E);
+    meas_update(E, s->E);
 
     if (n_at_max == 1) {
       for (q_t i = 0; i < q; i++) {
-        update_meas(sM[max_M_i][i], s->M[i]);
+        meas_update(sM[max_M_i][i], s->M[i]);
       }
-      update_meas(sE[max_M_i], s->E);
+      meas_update(sE[max_M_i], s->E);
       freqs[max_M_i]++;
     }
 
-    diff = fabs(sM[0][0]->dc / sM[0][0]->c);
+    diff = fabs(meas_dc(sM[0][0]) / meas_c(sM[0][0]));
 
     n_runs++;
   }
@@ -266,7 +266,7 @@ int main(int argc, char *argv[]) {
   }
   printf("WOLFF: 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[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x);
+         n_runs, fabs(meas_dx(E) / E->x), meas_dx(M[0]) / M[0]->x, meas_dc(E) / meas_c(E), meas_dc(M[0]) / meas_c(M[0]), h->nv / clust->x);
 
   if (snapshots) {
     FILE *snapfile = fopen("snapshots.m", "a");
@@ -359,7 +359,7 @@ int main(int argc, char *argv[]) {
       fprintf(outfile, ",");
     }
   }
-  fprintf(outfile, "},E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->{", E->x / h->nv, E->dx / h->nv, E->c / h->nv, E->dc / h->nv);
+  fprintf(outfile, "},E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->{", E->x / h->nv, meas_dx(E) / h->nv, meas_c(E) / h->nv, meas_dc(E) / h->nv);
   for (q_t i = 0; i < q; i++) {
     fprintf(outfile, "%.15f", M[i]->x / h->nv);
     if (i != q-1) {
@@ -368,27 +368,27 @@ int main(int argc, char *argv[]) {
   }
   fprintf(outfile, "},\\[Delta]M->{");
   for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", M[i]->dx / h->nv);
+    fprintf(outfile, "%.15f", meas_dx(M[i]) / h->nv);
     if (i != q-1) {
       fprintf(outfile, ",");
     }
   }
   fprintf(outfile, "},\\[Chi]->{");
   for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", M[i]->c / h->nv);
+    fprintf(outfile, "%.15f", meas_c(M[i]) / h->nv);
     if (i != q-1) {
       fprintf(outfile, ",");
     }
   }
   fprintf(outfile, "},\\[Delta]\\[Chi]->{");
   for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", M[i]->dc / h->nv);
+    fprintf(outfile, "%.15f", meas_dc(M[i]) / h->nv);
     if (i != q-1) {
       fprintf(outfile, ",");
     }
   }
   for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "},Subscript[E,%" PRIq "]->%.15f,Subscript[\\[Delta]E,%" PRIq "]->%.15f,Subscript[C,%" PRIq "]->%.15f,Subscript[\\[Delta]C,%" PRIq "]->%.15f,Subscript[M,%" PRIq "]->{", i, sE[i]->x / h->nv, i, sE[i]->dx / h->nv, i, sE[i]->c / h->nv, i, sE[i]->dc / h->nv, i);
+    fprintf(outfile, "},Subscript[E,%" PRIq "]->%.15f,Subscript[\\[Delta]E,%" PRIq "]->%.15f,Subscript[C,%" PRIq "]->%.15f,Subscript[\\[Delta]C,%" PRIq "]->%.15f,Subscript[M,%" PRIq "]->{", i, sE[i]->x / h->nv, i, meas_dx(sE[i]) / h->nv, i, meas_c(sE[i]) / h->nv, i, meas_dc(sE[i]) / h->nv, i);
     for (q_t j = 0; j < q; j++) {
       fprintf(outfile, "%.15f", sM[i][j]->x / h->nv);
       if (j != q-1) {
@@ -397,21 +397,21 @@ int main(int argc, char *argv[]) {
     }
     fprintf(outfile, "},Subscript[\\[Delta]M,%" PRIq "]->{", i);
     for (q_t j = 0; j < q; j++) {
-      fprintf(outfile, "%.15f", sM[i][j]->dx / h->nv);
+      fprintf(outfile, "%.15f", meas_dx(sM[i][j]) / h->nv);
       if (j != q-1) {
         fprintf(outfile, ",");
       }
     }
     fprintf(outfile, "},Subscript[\\[Chi],%" PRIq "]->{", i);
     for (q_t j = 0; j < q; j++) {
-      fprintf(outfile, "%.15f", sM[i][j]->c / h->nv);
+      fprintf(outfile, "%.15f", meas_c(sM[i][j]) / h->nv);
       if (j != q-1) {
         fprintf(outfile, ",");
       }
     }
     fprintf(outfile, "},Subscript[\\[Delta]\\[Chi],%" PRIq "]->{", i);
     for (q_t j = 0; j < q; j++) {
-      fprintf(outfile, "%.15f", sM[i][j]->dc / h->nv);
+      fprintf(outfile, "%.15f", meas_dc(sM[i][j]) / h->nv);
       if (j != q-1) {
         fprintf(outfile, ",");
       }
@@ -422,9 +422,9 @@ int main(int argc, char *argv[]) {
     fprintf(outfile, ",Subscript[f,%" PRIq "]->%.15f,Subscript[\\[Delta]f,%" PRIq "]->%.15f", i, (double)freqs[i] / (double)n_runs, i, sqrt(freqs[i]) / (double)n_runs);
   }
   for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, ",Subscript[t,%" PRIq "]->%.15f,Subscript[\\[Delta]t,%" PRIq "]->%.15f", i, lifetimes[i]->x, i, lifetimes[i]->dx);
+    fprintf(outfile, ",Subscript[t,%" PRIq "]->%.15f,Subscript[\\[Delta]t,%" PRIq "]->%.15f", i, lifetimes[i]->x, i, meas_dx(lifetimes[i]));
   }
-  fprintf(outfile, ",Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", clust->x / h->nv, clust->dx / h->nv, clust->c / h->nv, clust->dc / h->nv,tau);
+  fprintf(outfile, ",Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau);
 
   fclose(outfile);
 
diff --git a/src/wolff_vector.c b/src/wolff_vector.c
index 90ce7fc..c87e60d 100644
--- a/src/wolff_vector.c
+++ b/src/wolff_vector.c
@@ -158,7 +158,7 @@ int main(int argc, char *argv[]) {
   while (((diff > eps || diff != diff) && n_runs < N) || n_runs < min_runs) {
     if (!silent) 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), aM->dx / aM->x, E->dc / E->c, aM->dc / aM->c, h->nv / clust->x);
+           n_runs, fabs(meas_dx(E) / E->x), meas_dx(aM) / aM->x, meas_dc(E) / meas_c(E), meas_dc(aM) / meas_c(aM), h->nv / clust->x);
 
     count_t n_flips = 0;
 
@@ -172,34 +172,35 @@ int main(int argc, char *argv[]) {
 
       if (n_runs > 0) {
         n_steps++;
-        update_meas(clust, tmp_flips);
-      }
+        meas_update(clust, tmp_flips);
 
-      if (record_autocorrelation && n_runs > 0) {
-        if (n_steps % ac_skip == 0) {
+        if (record_autocorrelation && n_steps % ac_skip == 0) {
           update_autocorr(autocorr, s->E);
         }
       }
     }
 
     double aM_val = 0;
+
     for (q_t i = 0; i < q; i++) {
-      update_meas(M[i], s->M[i]);
+      meas_update(M[i], s->M[i]);
       aM_val += s->M[i] * s->M[i];
     }
-    update_meas(aM, sqrt(aM_val));
-    update_meas(E, s->E);
 
-    diff = fabs(aM->dc / aM->c);
+    meas_update(aM, sqrt(aM_val));
+    meas_update(E, s->E);
+
+    diff = fabs(meas_dc(aM) / meas_c(aM));
 
     n_runs++;
   }
+
   if (!silent) {
     printf("\033[F\033[J");
   }
   printf("WOLFF: 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[0]->dx / M[0]->x, E->dc / E->c, M[0]->dc / M[0]->c, h->nv / clust->x);
+         n_runs, fabs(meas_dx(E) / E->x), meas_dx(M[0]) / M[0]->x, meas_dc(E) / meas_c(E), meas_dc(M[0]) / meas_c(M[0]), h->nv / clust->x);
 
   double tau = 0;
   bool tau_failed = false;
@@ -260,7 +261,7 @@ int main(int argc, char *argv[]) {
       fprintf(outfile, ",");
     }
   }
-  fprintf(outfile, "},E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->{", E->x / h->nv, E->dx / h->nv, E->c / h->nv, E->dc / h->nv);
+  fprintf(outfile, "},E->%.15f,\\[Delta]E->%.15f,C->%.15f,\\[Delta]C->%.15f,M->{", E->x / h->nv, meas_dx(E) / h->nv, meas_c(E) / h->nv, meas_dc(E) / h->nv);
   for (q_t i = 0; i < q; i++) {
     fprintf(outfile, "%.15f", M[i]->x / h->nv);
     if (i != q-1) {
@@ -269,26 +270,26 @@ int main(int argc, char *argv[]) {
   }
   fprintf(outfile, "},\\[Delta]M->{");
   for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", M[i]->dx / h->nv);
+    fprintf(outfile, "%.15f", meas_dx(M[i]) / h->nv);
     if (i != q-1) {
       fprintf(outfile, ",");
     }
   }
   fprintf(outfile, "},\\[Chi]->{");
   for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", M[i]->c / h->nv);
+    fprintf(outfile, "%.15f", meas_c(M[i]) / h->nv);
     if (i != q-1) {
       fprintf(outfile, ",");
     }
   }
   fprintf(outfile, "},\\[Delta]\\[Chi]->{");
   for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", M[i]->dc / h->nv);
+    fprintf(outfile, "%.15f", meas_dc(M[i]) / h->nv);
     if (i != q-1) {
       fprintf(outfile, ",");
     }
   }
-  fprintf(outfile, "},aM->%.15f,\\[Delta]aM->%.15f,a\\[Chi]->%.15f,\\[Delta]a\\[Chi]->%.15f,Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv, clust->x / h->nv, clust->dx / h->nv, clust->c / h->nv, clust->dc / h->nv,tau);
+  fprintf(outfile, "},aM->%.15f,\\[Delta]aM->%.15f,a\\[Chi]->%.15f,\\[Delta]a\\[Chi]->%.15f,Subscript[n,\"clust\"]->%.15f,Subscript[\\[Delta]n,\"clust\"]->%.15f,Subscript[m,\"clust\"]->%.15f,Subscript[\\[Delta]m,\"clust\"]->%.15f,\\[Tau]->%.15f|>\n", aM->x / h->nv, meas_dx(aM) / h->nv, meas_c(aM) / h->nv, meas_dc(aM) / h->nv, clust->x / h->nv, meas_dx(clust) / h->nv, meas_c(clust) / h->nv, meas_dc(clust) / h->nv,tau);
 
   fclose(outfile);