got everyone recording data, and fixed huge bug in the updating of ReF and ImF

7 files changed, 169 insertions, 79 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index b8c98e5..3261969 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -77,11 +77,11 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
           for (D_t i = 0; i < state->D; i++) {
             L_t x = (non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L;
 
-            add(&(state->ReF[i]), -state->precomputed_cos[i], rs_old);
-            add(&(state->ReF[i]),  state->precomputed_cos[i], rs_new);
+            add(&(state->ReF[i]), -state->precomputed_cos[x], rs_old);
+            add(&(state->ReF[i]),  state->precomputed_cos[x], rs_new);
 
-            add(&(state->ImF[i]), -state->precomputed_sin[i], rs_old);
-            add(&(state->ImF[i]),  state->precomputed_sin[i], rs_new);
+            add(&(state->ImF[i]), -state->precomputed_sin[x], rs_old);
+            add(&(state->ImF[i]),  state->precomputed_sin[x], rs_new);
           }
 
           free_spin (rs_old);
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index ce2d300..0ba5eee 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -8,6 +8,7 @@
 
 #include "state.h"
 #include "types.h"
+#include "vector.h"
 
 template <q_t q, class T>
 struct orthogonal_t { bool is_reflection; T *x; };
diff --git a/lib/potts.h b/lib/potts.h
index 732a76f..e61e4e1 100644
--- a/lib/potts.h
+++ b/lib/potts.h
@@ -4,6 +4,7 @@
 #include <stdio.h>
 
 #include "types.h"
+#include "vector.h"
 
 /* The following is the minimum definition of a spin class.
  *
@@ -30,8 +31,8 @@ class potts_t {
   public:
     q_t x;
 
-    typedef int *M_t;
-    typedef double *F_t;
+    typedef vector_t<q, int> M_t;
+    typedef vector_t<q, double> F_t;
 };
 
 template <q_t q>
@@ -44,56 +45,54 @@ void free_spin(potts_t <q> s) {
   // do nothing!
 }
 
-void free_spin(int *s) {
-  free(s);
-}
-
-void free_spin(double *s) {
-  free(s);
-}
-
 template <q_t q>
 potts_t <q> copy(potts_t <q> s) {
   return s;
 }
 
 template <q_t q>
-void add(typename potts_t<q>::M_t *s1, int a, potts_t <q> s2) {
-  (*s1)[s2.x] += a;
+void add(vector_t<q, int> *s1, int a, potts_t <q> s2) {
+  (s1->x)[s2.x] += a;
 }
 
 template <q_t q>
-void add(typename potts_t<q>::F_t *s1, double a, potts_t <q> s2) {
-  (*s1)[s2.x] += a;
+void add(vector_t<q, double> *s1, double a, potts_t <q> s2) {
+  (s1->x)[s2.x] += a;
 }
 
 template <q_t q>
-typename potts_t<q>::M_t scalar_multiple(int factor, potts_t <q> s) {
-  int *M = (int *)calloc(q, sizeof(int));
-  M[s.x] += factor;
+vector_t<q, int> scalar_multiple(int factor, potts_t <q> s) {
+  vector_t<q, int> M;
+  M.x = (int *)calloc(q, sizeof(int));
+  M.x[s.x] += factor;
   return M;
 }
 
 template <q_t q>
-typename potts_t<q>::F_t scalar_multiple(double factor, potts_t <q> s) {
-  double *F = (double *)calloc(q, sizeof(double));
-  F[s.x] += factor;
+vector_t<q, double> scalar_multiple(double factor, potts_t <q> s) {
+  vector_t<q, double> F;
+  F.x = (double *)calloc(q, sizeof(double));
+  F.x[s.x] += factor;
   return F;
 }
 
+// we could inherit norm_squared from vector.h, but convention dictates that
+// potts norms be changed by a constant factor
 template <q_t q>
-double norm_squared(typename potts_t<q>::F_t s) {
+double norm_squared(vector_t<q, double> s) {
   double total = 0;
   for (q_t i = 0; i < q; i++) {
-    total += pow(s[i], 2);
+    total += pow(s.x[i], 2);
   }
 
   return total * (double)q / ((double)q - 1.0);
 }
 
+// we could inherit write_magnetization from vector.h, but since M.x must sum
+// to nv we don't need to write the last element
 template <q_t q>
-void write_magnetization(typename potts_t<q>::M_t M, FILE *outfile) {
-  fwrite(&M, sizeof(int), q, outfile);
+void write_magnetization(vector_t<q, int> M, FILE *outfile) {
+  fwrite(M.x, sizeof(int), q - 1, outfile);
 }
 
 // knock yourself out
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index 997ec09..a277e0f 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -6,11 +6,10 @@
 #include <GL/glut.h>
 #endif
 
-#include <vector.h>
 #include <orthogonal.h>
+#include <vector.h>
 
 #include <wolff.h>
-#include <correlation.h>
 #include <measure.h>
 #include <colors.h>
 #include <rand.h>
@@ -235,7 +234,6 @@ int main(int argc, char *argv[]) {
     wolff <orthogonal_R_t, vector_R_t> (N, &s, gen_R, measurements, r, silent);
   }
 
-
   measure_free_files(measurement_flags, outfiles);
   free(H_vec);
   gsl_rng_free(r);
diff --git a/src/wolff_clock.cpp b/src/wolff_clock.cpp
index e186c44..86badfe 100644
--- a/src/wolff_clock.cpp
+++ b/src/wolff_clock.cpp
@@ -107,7 +107,7 @@ int main(int argc, char *argv[]) {
   if (!draw) {
     // a very simple example: measure the average magnetization
     measurement = [&] (const sim_t *s) {
-      average_M += (double)s->M[0] / (double)N / (double)s->nv;
+      average_M += (double)s->M.x[0] / (double)N / (double)s->nv;
     };
   } else {
     // a more complex example: measure the average magnetization, and draw the spin configuration to the screen
@@ -124,7 +124,7 @@ int main(int argc, char *argv[]) {
     gluOrtho2D(0.0, L, 0.0, L);
 
     measurement = [&] (const sim_t *s) {
-      average_M += (double)s->M[0] / (double)N / (double)s->nv;
+      average_M += (double)s->M.x[0] / (double)N / (double)s->nv;
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(L, 2); i++) {
         potts_t<POTTSQ> tmp_s = act_inverse(s->R, s->spins[i]);
diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp
index 7492ebf..5e44cab 100644
--- a/src/wolff_ising.cpp
+++ b/src/wolff_ising.cpp
@@ -1,5 +1,8 @@
 
 #include <getopt.h>
+#include <stdio.h>
+
+// if you have GLUT installed, you can see graphics!
 #ifdef HAVE_GLUT
 #include <GL/glut.h>
 #endif
@@ -8,10 +11,17 @@
 #include <z2.h>
 #include <ising.h>
 
+// finite_states.h can be included for spin types that have special variables
+// defined, and it causes wolff execution to use precomputed bond probabilities
 #include <finite_states.h>
 
-// include wolff.h
+// rand.h uses a unix-specific way to seed the random number generator
 #include <rand.h>
+
+// measure.h contains useful functions for saving timeseries to files
+#include <measure.h>
+
+// include wolff.h
 #include <wolff.h>
 
 int main(int argc, char *argv[]) {
@@ -25,11 +35,15 @@ int main(int argc, char *argv[]) {
 
   bool silent = false;
   bool draw = false;
+  bool N_is_sweeps = false;
   unsigned int window_size = 512;
 
+  // don't measure anything by default
+  unsigned char measurement_flags = 0;
+
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -1) {
     switch (opt) {
     case 'N': // number of steps
       N = (count_t)atof(optarg);
@@ -49,6 +63,9 @@ int main(int argc, char *argv[]) {
     case 's': // don't print anything during simulation. speeds up slightly
       silent = true;
       break;
+    case 'S':
+      N_is_sweeps = true;
+      break;
     case 'd':
 #ifdef HAVE_GLUT
       draw = true;
@@ -60,11 +77,23 @@ int main(int argc, char *argv[]) {
     case 'w':
       window_size = atoi(optarg);
       break;
+    case 'M':
+      measurement_flags ^= 1 << atoi(optarg);
+      break;
     default:
       exit(EXIT_FAILURE);
     }
   }
 
+  // get nanosecond timestamp for unique run id
+  unsigned long timestamp;
+
+  {
+    struct timespec spec;
+    clock_gettime(CLOCK_REALTIME, &spec);
+    timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec;
+  }
+
   // initialize random number generator
   gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
   gsl_rng_set(r, rand_seed());
@@ -97,18 +126,16 @@ int main(int argc, char *argv[]) {
     return rot;
   };
 
-  // define function that updates any number of measurements
-  std::function <void(const state_t <z2_t, ising_t> *)> measurement;
+  FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
 
-  double average_M = 0;
-  if (!draw) {
-    // a very simple example: measure the average magnetization
-    measurement = [&] (const state_t <z2_t, ising_t> *s) {
-      average_M += (double)s->M / (double)N / (double)s->nv;
-    };
-  } else {
-    // a more complex example: measure the average magnetization, and draw the spin configuration to the screen
+  std::function <void(const state_t<z2_t, ising_t> *)> other_f;
+  uint64_t sum_of_clusterSize = 0;
 
+  if (N_is_sweeps) {
+    other_f = [&] (const state_t<z2_t, ising_t> *s) {
+      sum_of_clusterSize += s->last_cluster_size;
+    };
+  } else if (draw) {
 #ifdef HAVE_GLUT
     // initialize glut
     glutInit(&argc, argv);
@@ -120,34 +147,52 @@ int main(int argc, char *argv[]) {
     glLoadIdentity();
     gluOrtho2D(0.0, L, 0.0, L);
 
-    measurement = [&] (const state_t <z2_t, ising_t> *s) {
-      average_M += (double)s->M / (double)N / (double)s->nv;
+    other_f = [] (const state_t <z2_t, ising_t> *s) {
       glClear(GL_COLOR_BUFFER_BIT);
-      for (v_t i = 0; i < pow(L, 2); i++) {
+      for (v_t i = 0; i < pow(s->L, 2); i++) {
         if (s->spins[i].x == s->R.x) {
           glColor3f(0.0, 0.0, 0.0);
         } else {
           glColor3f(1.0, 1.0, 1.0);
         }
-        glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1);
+        glRecti(i / s->L, i % s->L, (i / s->L) + 1, (i % s->L) + 1);
       }
       glFlush();
     };
 #endif
+  } else {
+    other_f = [] (const state_t<z2_t, ising_t> *s) {};
   }
 
-  // run wolff for N cluster flips
-  wolff(N, &s, gen_R, measurement, r, silent);
+  std::function <void(const state_t<z2_t, ising_t> *)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f);
 
-  // tell us what we found!
-  printf("%" PRIcount " Ising runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, <M> = %g\n", N, D, L, T, H, average_M);
+  // add line to metadata file with run info
+  {
+    FILE *outfile_info = fopen("wolff_metadata.txt", "a");
 
-  // free the random number generator
-  gsl_rng_free(r);
+    fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"ISING\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> %.15f |>\n", timestamp, s.D, s.L, s.nv, s.ne, T, H);
 
-  if (draw) {
+    fclose(outfile_info);
   }
 
+  // run wolff for N cluster flips
+  if (N_is_sweeps) {
+    count_t N_rounds = 0;
+    printf("\n");
+    while (sum_of_clusterSize < N * s.nv) {
+      printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
+      wolff(N, &s, gen_R, measurements, r, silent);
+      N_rounds++;
+    }
+    printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
+  } else {
+    wolff(N, &s, gen_R, measurements, r, silent);
+  }
+
+  // free the random number generator
+  gsl_rng_free(r);
+  measure_free_files(measurement_flags, outfiles);
+
   return 0;
 
 }
diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp
index f8f1523..cdc4c07 100644
--- a/src/wolff_potts.cpp
+++ b/src/wolff_potts.cpp
@@ -1,5 +1,6 @@
 
 #include <getopt.h>
+#include <stdio.h>
 
 #ifdef HAVE_GLUT
 #include <GL/glut.h>
@@ -8,13 +9,14 @@
 // include your group and spin space
 #include <symmetric.h>
 #include <potts.h>
-#include <colors.h>
 
 // hack to speed things up considerably
 #define N_STATES POTTSQ
 #include <finite_states.h>
 
 // include wolff.h
+#include <measure.h>
+#include <colors.h>
 #include <rand.h>
 #include <wolff.h>
 
@@ -31,12 +33,16 @@ int main(int argc, char *argv[]) {
 
   bool silent = false;
   bool draw = false;
+  bool N_is_sweeps = false;
   unsigned int window_size = 512;
 
+  // don't measure anything by default
+  unsigned char measurement_flags = 0;
+
   int opt;
   q_t H_ind = 0;
 
-  while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -1) {
     switch (opt) {
     case 'N': // number of steps
       N = (count_t)atof(optarg);
@@ -57,6 +63,9 @@ int main(int argc, char *argv[]) {
     case 's': // don't print anything during simulation. speeds up slightly
       silent = true;
       break;
+    case 'S':
+      N_is_sweeps = true;
+      break;
     case 'd':
 #ifdef HAVE_GLUT
       draw = true;
@@ -68,11 +77,23 @@ int main(int argc, char *argv[]) {
     case 'w':
       window_size = atoi(optarg);
       break;
+    case 'M':
+      measurement_flags ^= 1 << atoi(optarg);
+      break;
     default:
       exit(EXIT_FAILURE);
     }
   }
 
+  // get nanosecond timestamp for unique run id
+  unsigned long timestamp;
+
+  {
+    struct timespec spec;
+    clock_gettime(CLOCK_REALTIME, &spec);
+    timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec;
+  }
+
   // initialize random number generator
   gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
   gsl_rng_set(r, rand_seed());
@@ -113,18 +134,16 @@ int main(int argc, char *argv[]) {
     return rot;
   };
 
-  // define function that updates any number of measurements
-  std::function <void(const sim_t *)> measurement;
+  FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
 
-  double average_M = 0;
-  if (!draw) {
-    // a very simple example: measure the average magnetization
-    measurement = [&] (const sim_t *s) {
-      average_M += (double)s->M[0] / (double)N / (double)s->nv;
-    };
-  } else {
-    // a more complex example: measure the average magnetization, and draw the spin configuration to the screen
+  std::function <void(const sim_t *)> other_f;
+  uint64_t sum_of_clusterSize = 0;
 
+  if (N_is_sweeps) {
+    other_f = [&] (const sim_t *s) {
+      sum_of_clusterSize += s->last_cluster_size;
+    };
+  } else if (draw) {
 #ifdef HAVE_GLUT
     // initialize glut
     glutInit(&argc, argv);
@@ -136,31 +155,59 @@ int main(int argc, char *argv[]) {
     glLoadIdentity();
     gluOrtho2D(0.0, L, 0.0, L);
 
-    measurement = [&] (const sim_t *s) {
-      average_M += (double)s->M[0] / (double)N / (double)s->nv;
+    other_f = [] (const sim_t *s) {
       glClear(GL_COLOR_BUFFER_BIT);
-      for (v_t i = 0; i < pow(L, 2); i++) {
+      for (v_t i = 0; i < pow(s->L, 2); i++) {
         potts_t<POTTSQ> tmp_s = act_inverse(s->R, s->spins[i]);
         glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ));
-        glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1);
+        glRecti(i / s->L, i % s->L, (i / s->L) + 1, (i % s->L) + 1);
       }
       glFlush();
     };
 #endif
+  } else {
+    other_f = [] (const sim_t *s) {};
   }
 
-  // run wolff for N cluster flips
-  wolff(N, &s, gen_R, measurement, r, silent);
+  std::function <void(const sim_t *)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f);
 
-  // tell us what we found!
-  printf("%" PRIcount " %d-Potts runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, <M> = %g\n", N, POTTSQ, D, L, T, H_vec[0], average_M);
+  // add line to metadata file with run info
+  {
+    FILE *outfile_info = fopen("wolff_metadata.txt", "a");
 
-  // free the random number generator
-  gsl_rng_free(r);
+    fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"POTTS\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, POTTSQ, s.D, s.L, s.nv, s.ne, T);
+
+    for (q_t i = 0; i < POTTSQ; i++) {
+      fprintf(outfile_info, "%.15f", H_vec[i]);
+      if (i < POTTSQ - 1) {
+        fprintf(outfile_info, ", ");
+      }
+    }
 
-  if (draw) {
+    fprintf(outfile_info, "} |>\n");
+
+    fclose(outfile_info);
   }
 
+  // run wolff for N cluster flips
+  if (N_is_sweeps) {
+    count_t N_rounds = 0;
+    printf("\n");
+    while (sum_of_clusterSize < N * s.nv) {
+      printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
+      wolff(N, &s, gen_R, measurements, r, silent);
+      N_rounds++;
+    }
+    printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size);
+  } else {
+    wolff(N, &s, gen_R, measurements, r, silent);
+  }
+
+  // free the random number generator
+  gsl_rng_free(r);
+  free(H_vec);
+  measure_free_files(measurement_flags, outfiles);
+
   return 0;
 
 }