many changes, including reworking the measurements system

10 files changed, 400 insertions, 100 deletions

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index baf0a10..a435b29 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -3,5 +3,8 @@ add_library(wolff_examples INTERFACE)
 
 target_include_directories(wolff_examples INTERFACE include)
 
+add_executable(simple_ising simple_ising.cpp)
+target_link_libraries(simple_ising wolff)
+
 add_subdirectory(src)
 
diff --git a/examples/include/correlation.hpp b/examples/include/correlation.hpp
index 042cff3..da8ab7f 100644
--- a/examples/include/correlation.hpp
+++ b/examples/include/correlation.hpp
@@ -6,18 +6,18 @@
 
 #include <fftw3.h>
 
-template <class R_t, class X_t>
-double correlation_length(const state_t <R_t, X_t>& s) {
+template <class X_t>
+double correlation_length(const std::vector<typename X_t::F_t>& ReF, const std::vector<typename X_t::F_t>& ImF, D_t D) {
   double total = 0;
 
 #ifdef DIMENSION
   for (D_t j = 0; j < DIMENSION; j++) {
 #else
-  for (D_t j = 0; j < s.D; j++) {
+  for (D_t j = 0; j < D; j++) {
 #endif
-    total += norm_squared(s.ReF[j]) + norm_squared(s.ImF[j]);
+    total += norm_squared(ReF[j]) + norm_squared(ImF[j]);
   }
 
-  return total / s.D;
+  return total / D;
 }
 
diff --git a/examples/include/measure.hpp b/examples/include/measure.hpp
index e20353c..bf8baab 100644
--- a/examples/include/measure.hpp
+++ b/examples/include/measure.hpp
@@ -2,6 +2,7 @@
 #pragma once
 
 #include <wolff/state.hpp>
+#include <wolff/meas.h>
 #include "correlation.hpp"
 #include <functional>
 
@@ -13,51 +14,146 @@ const unsigned char measurement_fourierZero    = 1 << 3;
 
 char const *measurement_labels[] = {"E", "S", "M", "F"};
 
-FILE **measure_setup_files(unsigned char flags, unsigned long timestamp) {
-  FILE **files = (FILE **)calloc(POSSIBLE_MEASUREMENTS, sizeof(FILE *));
+template <class R_t, class X_t>
+class wolff_research_measurements : public wolff_measurement<R_t, X_t> {
+  private:
+    std::vector<std::vector<double>> precomputed_sin;
+    std::vector<std::vector<double>> precomputed_cos;
+    FILE **files;
+    D_t D;
+    unsigned char flags;
+    std::function <void(const state_t <R_t, X_t>&, const wolff_research_measurements<R_t, X_t>&)> other_f;
+    bool silent;
+  public:
+    v_t last_cluster_size;
+    double E;
+    typename X_t::M_t M;
+    std::vector<typename X_t::F_t> ReF;
+    std::vector<typename X_t::F_t> ImF;
+
+    wolff_research_measurements(unsigned char flags, unsigned long timestamp, std::function <void(const state_t <R_t, X_t>&, const wolff_research_measurements<R_t, X_t>&)> other_f, state_t<R_t, X_t>& s, bool silent) : flags(flags), D(s.D), other_f(other_f), silent(silent) {
+      FILE **files = (FILE **)calloc(POSSIBLE_MEASUREMENTS, sizeof(FILE *));
+
+      for (uint8_t i = 0; i < POSSIBLE_MEASUREMENTS; i++) {
+        if (flags & (1 << i)) {
+          char *filename = (char *)malloc(255 * sizeof(char));
+          sprintf(filename, "wolff_%lu_%s.dat", timestamp, measurement_labels[i]);
+          files[i] = fopen(filename, "wb");
+          free(filename);
+        }
+      }
+
+#ifdef BOND_DEPENDENCE
+      E = 0;
+      for (v_t v = 0; v < s.nv; v++) {
+        for (const v_t &vn : s.g.v_adj[v]) {
+          if (v < vn) {
+            E -= s.J(v, s.spins[v], vn, s.spins[vn]);
+          }
+        }
+      }
+#else
+      E = - (double)s.ne * s.J(s.spins[0], s.spins[0]);
+#endif
+
+#ifndef NOFIELD
+#ifdef SITE_DEPENDENCE
+      for (v_t i = 0; i < s.nv; i++) {
+        E -= s.H(i, s.spins[i]);
+      }
+#else
+      E -= (double)s.nv * s.H(s.spins[0]);
+#endif
+#endif
+
+      M = s.spins[0] * s.nv;
+
+      ReF.resize(s.D);
+      ImF.resize(s.D);
+      for (D_t i = 0; i < s.D; i++) {
+        ReF[i] = s.spins[0] * 0.0;
+        ImF[i] = s.spins[0] * 0.0;
+      }
+      precomputed_cos.resize(s.nv);
+      precomputed_sin.resize(s.nv);
+      for (v_t i = 0; i < s.nv; i++) {
+        precomputed_cos[i].resize(s.D);
+        precomputed_sin[i].resize(s.D);
+        for (D_t j = 0; j < s.D; j++) {
+          precomputed_cos[i][j] = cos(2 * M_PI * s.g.coordinate[i][j] / (double)s.L);
+          precomputed_sin[i][j] = sin(2 * M_PI * s.g.coordinate[i][j] / (double)s.L);
+        }
+      }
+
+      if (!silent) printf("\n");
 
-  for (uint8_t i = 0; i < POSSIBLE_MEASUREMENTS; i++) {
-    if (flags & (1 << i)) {
-      char *filename = (char *)malloc(255 * sizeof(char));
-      sprintf(filename, "wolff_%lu_%s.dat", timestamp, measurement_labels[i]);
-      files[i] = fopen(filename, "wb");
-      free(filename);
     }
-  }
 
-  return files;
-}
+    void pre_cluster(const state_t<R_t, X_t>& s, count_t step, count_t N, v_t v0, const R_t& R) {
+      last_cluster_size = 0;
+      if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", step, N, E, last_cluster_size);
 
-template <class R_t, class X_t>
-std::function <void(const state_t <R_t, X_t>&)> measure_function_write_files(unsigned char flags, FILE **files, std::function <void(const state_t <R_t, X_t>&)> other_f) {
-  return [=] (const state_t <R_t, X_t>& s) {
-    if (flags & measurement_energy) {
-      float smaller_E = (float)s.E;
-      fwrite(&smaller_E, sizeof(float), 1, files[0]);
     }
-    if (flags & measurement_clusterSize) {
-      fwrite(&(s.last_cluster_size), sizeof(uint32_t), 1, files[1]);
+
+    void plain_bond_added(v_t vi, const X_t& si_old, const X_t& si_new, v_t vn, const X_t& sn, double dE) {
+      E += dE;
     }
-    if (flags & measurement_magnetization) {
-      write_magnetization(s.M, files[2]);
+
+    void ghost_bond_added(v_t v, const X_t& rs_old, const X_t& rs_new, double dE) {
+      E += dE;
+      M += rs_new - rs_old;
+
+#ifdef DIMENSION
+      for (D_t i = 0; i < DIMENSION; i++)
+#else
+      for (D_t i = 0; i < D; i++)
+#endif
+      {
+        ReF[i] += (rs_new - rs_old) * precomputed_cos[v][i];
+        ImF[i] += (rs_new - rs_old) * precomputed_sin[v][i];
+      }
     }
-    if (flags & measurement_fourierZero) {
-      float smaller_X = (float)correlation_length(s);
-      fwrite(&smaller_X, sizeof(float), 1, files[3]);
+
+    void plain_site_transformed(v_t v, const X_t& s_old, const X_t& s_new) {
+      last_cluster_size++;
     }
 
-    other_f(s);
-  };
-}
+    void ghost_site_transformed(const R_t& r_old, const R_t& r_new) {
+    }
+
+    void post_cluster(const state_t<R_t, X_t>& s, count_t step, count_t N) {
+      if (flags & measurement_energy) {
+        float smaller_E = (float)E;
+        fwrite(&smaller_E, sizeof(float), 1, files[0]);
+      }
+      if (flags & measurement_clusterSize) {
+        fwrite(&(last_cluster_size), sizeof(uint32_t), 1, files[1]);
+      }
+      if (flags & measurement_magnetization) {
+        write_magnetization(M, files[2]);
+      }
+      if (flags & measurement_fourierZero) {
+        float smaller_X = (float)correlation_length<X_t>(ReF, ImF, D);
+        fwrite(&smaller_X, sizeof(float), 1, files[3]);
+      }
+
+      other_f(s, *this);
 
-void measure_free_files(unsigned char flags, FILE **files) {
-  for (uint8_t i = 0; i < POSSIBLE_MEASUREMENTS; i++) {
-    if (flags & (1 << i)) {
-      fclose(files[i]);
     }
-  }
 
-  free(files);
-}
+    ~wolff_research_measurements() {
+      for (uint8_t i = 0; i < POSSIBLE_MEASUREMENTS; i++) {
+        if (flags & (1 << i)) {
+          fclose(files[i]);
+        }
+      }
 
+      if (!silent) {
+        printf("\033[F\033[J");
+      }
+      printf("WOLFF COMPLETE: E = %.2f, S = %" PRIv "\n", E, last_cluster_size);
+
+      free(files);
+    }
+};
 
diff --git a/examples/simple_ising.cpp b/examples/simple_ising.cpp
new file mode 100644
index 0000000..24e4ae5
--- /dev/null
+++ b/examples/simple_ising.cpp
@@ -0,0 +1,192 @@
+
+#include <getopt.h>
+#include <iostream>
+
+#include "include/randutils/randutils.hpp"
+
+#include <wolff.hpp>
+
+// define your R_t and X_t. here both are the same, ising_t
+class ising_t {
+  public:
+    bool x;
+
+    // both X_t and R_t need a default constructor (and destructor, if relevant)
+    ising_t() : x(false) {}
+    ising_t(bool x) : x(x) {}
+
+    // R_t needs the member functions
+    //   X_t act(const X_t& s) const {}
+    //   R_t act(const R_t& s) const {}
+    // to define action on both spins and other transformations
+    ising_t act(const ising_t& s) const {
+      if (x) {
+        return ising_t(!(s.x));
+      } else {
+        return ising_t(s.x);
+      }
+    }
+
+    // R_t needs the member functions
+    //   X_t act_inverse(const X_t& s) const {}
+    //   R_t act_inverse(const R_t& s) const {}
+    // to define action of its inverse on both spins and other transformations
+    ising_t act_inverse(const ising_t& s) const {
+      return this->act(s);
+    }
+};
+
+// define how measurements should be taken by importing the interface wolff_measurement<R_t, X_t>
+class ising_measurements : public wolff_measurement<ising_t, ising_t> {
+  private:
+    count_t n;
+
+    double E;
+    int M;
+    v_t S;
+
+    double totalE;
+    double totalM;
+    double totalS;
+
+  public:
+    ising_measurements(D_t D, L_t L, double H) {
+      n = 0;
+      M = (int)pow(L, D);
+      E = -D * pow(L, D) - H * pow(L, D);
+
+      totalE = 0;
+      totalM = 0;
+      totalS = 0;
+    }
+
+    void pre_cluster(const state_t<ising_t, ising_t>& s, count_t step, count_t N, v_t v0, const ising_t& R) {
+      S = 0;
+    }
+
+    void plain_bond_added(v_t v, const ising_t& s_old, const ising_t& s_new, v_t vn, const ising_t& sn, double dE) {
+      E += dE;
+    }
+
+    void ghost_bond_added(v_t v, const ising_t& s_old, const ising_t& s_new, double dE) {
+      E += dE;
+      
+      if (s_old.x) {
+        M++;
+      } else {
+        M--;
+      }
+
+      if (s_new.x) {
+        M--;
+      } else {
+        M++;
+      }
+    }
+
+    void plain_site_transformed(v_t v, const ising_t& s_old, const ising_t& s_new) {
+      S++;
+    }
+
+    void ghost_site_transformed(const ising_t& R_old, const ising_t& R_new) {
+    }
+
+    void post_cluster(const state_t<ising_t, ising_t>& s, count_t step, count_t N) {
+      totalE += E;
+      totalM += M;
+      totalS += S;
+      n++;
+    }
+
+    double avgE() {
+      return totalE / n;
+    }
+
+    double avgM() {
+      return totalM / n;
+    }
+
+    double avgS() {
+      return totalS / n;
+    }
+};
+
+int main(int argc, char *argv[]) {
+
+  // set defaults
+  count_t N = (count_t)1e4;
+  D_t D = 2;
+  L_t L = 128;
+  double T = 2.26918531421;
+  double H = 0.0;
+
+  int opt;
+
+  // take command line arguments
+  while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
+    switch (opt) {
+    case 'N': // number of steps
+      N = (count_t)atof(optarg);
+      break;
+    case 'D': // dimension
+      D = atoi(optarg);
+      break;
+    case 'L': // linear size
+      L = atoi(optarg);
+      break;
+    case 'T': // temperature
+      T = atof(optarg);
+      break;
+    case 'H': // external field
+      H = atof(optarg);
+      break;
+    default:
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  // define the spin-spin coupling
+  std::function <double(const ising_t&, const ising_t&)> Z = [] (const ising_t& s1, const ising_t& s2) -> double {
+    if (s1.x == s2.x) {
+      return 1.0;
+    } else {
+      return -1.0;
+    }
+  };
+
+  // define the spin-field coupling
+  std::function <double(const ising_t&)> B = [=] (const ising_t& s) -> double {
+    if (s.x) {
+      return -H;
+    } else {
+      return H;
+    }
+  };
+
+  // initialize the system
+  state_t<ising_t, ising_t> s(D, L, T, Z, B);
+
+  // initialize the random number generator
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng{seeds};
+
+  // define function that generates self-inverse rotations
+  std::function <ising_t(std::mt19937&, const ising_t&)> gen_R = [] (std::mt19937&, const ising_t& s) -> ising_t {
+    return ising_t(true);
+  };
+
+  // initailze the measurement object
+  ising_measurements m(D, L, H);
+
+  // run wolff N times
+  wolff<ising_t, ising_t>(N, s, gen_R, m, rng);
+
+  // print the result of our measurements
+  std::cout << "Wolff complete!\nThe average energy per site was " << m.avgE() / s.nv
+    << ".\nThe average magnetization per site was " << m.avgM() / s.nv
+    << ".\nThe average cluster size per site was " << m.avgS() / s.nv << ".\n";
+
+  // exit
+  return 0;
+}
+
diff --git a/examples/src/models/On/wolff_On.cpp b/examples/src/models/On/wolff_On.cpp
index ad2ac77..67f28a5 100644
--- a/examples/src/models/On/wolff_On.cpp
+++ b/examples/src/models/On/wolff_On.cpp
@@ -188,14 +188,12 @@ int main(int argc, char *argv[]) {
 
   fclose(outfile_info);
 
-  FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
-
-  std::function <void(const On_t&)> other_f;
+  std::function <void(const On_t&, const wolff_research_measurements<orthogonal_R_t, vector_R_t>&)> other_f;
   uint64_t sum_of_clusterSize = 0;
 
   if (N_is_sweeps) {
-    other_f = [&] (const On_t& s) {
-      sum_of_clusterSize += s.last_cluster_size;
+    other_f = [&] (const On_t& s, const wolff_research_measurements<orthogonal_R_t, vector_R_t>& m) {
+      sum_of_clusterSize += m.last_cluster_size;
     };
   } else if (draw) {
 #ifdef HAVE_GLUT
@@ -209,7 +207,7 @@ int main(int argc, char *argv[]) {
     glLoadIdentity();
     gluOrtho2D(0.0, L, 0.0, L);
 
-    other_f = [&] (const On_t& s) {
+    other_f = [&] (const On_t& s, const wolff_research_measurements<orthogonal_R_t, vector_R_t>& m) {
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(L, 2); i++) {
 #ifdef NOFIELD
@@ -228,11 +226,9 @@ int main(int argc, char *argv[]) {
     };
 #endif
   } else {
-    other_f = [] (const On_t& s) {};
+    other_f = [] (const On_t& s, const wolff_research_measurements<orthogonal_R_t, vector_R_t>& m) {};
   }
 
-  std::function <void(const On_t&)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f);
-
   std::function <double(const vector_R_t&)> H;
 
   if (modulated_field) {
@@ -251,20 +247,21 @@ int main(int argc, char *argv[]) {
   state_t <orthogonal_R_t, vector_R_t> s(D, L, T, dot <N_COMP, double>);
 #endif
 
+  wolff_research_measurements<orthogonal_R_t, vector_R_t> m(measurement_flags, timestamp, other_f, s, silent);
+
   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 <orthogonal_R_t, vector_R_t> (N, s, gen_R, measurements, rng, silent);
+      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, m.E, m.last_cluster_size);
+      wolff <orthogonal_R_t, vector_R_t> (N, s, gen_R, m, rng);
       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);
+    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, m.E, m.last_cluster_size);
   } else {
-    wolff <orthogonal_R_t, vector_R_t> (N, s, gen_R, measurements, rng, silent);
+    wolff <orthogonal_R_t, vector_R_t> (N, s, gen_R, m, rng);
   }
 
-  measure_free_files(measurement_flags, outfiles);
   free(H_vec);
 
   return 0;
diff --git a/examples/src/models/ising/ising.hpp b/examples/src/models/ising/ising.hpp
index ae20840..73b06ed 100644
--- a/examples/src/models/ising/ising.hpp
+++ b/examples/src/models/ising/ising.hpp
@@ -5,17 +5,31 @@
 
 #include <wolff/types.h>
 
+// all that is required to use wolff.hpp is a default constructor
 class ising_t {
   public:
     bool x;
 
-    typedef int M_t;
-    typedef double F_t;
-
     ising_t() : x(false) {}
-    ising_t(bool x) : x(x) {}
+
+    // optional constructors for syntactic sugar
+    ising_t(bool x) : x(x) {} 
     ising_t(int x) : x((bool)x) {}
 
+    /* below this comment is code required only for using measure.hpp in the
+     * examples folder, which provides an interface for measuring several
+     * generic features of models. these require
+     *
+     *  - an M_t, representing the magnetization or sum of all spins
+     *  - an F_t, representing a double-weighted version of the magnetization
+     *  - the overloaded operator *, which takes a v_t (unsigned int) and returns an M_t
+     *  - the overloaded operator *, which takes a double and returns an F_t
+     *  - the overloaded operator -, which takes another X_t and returns an M_t
+     */
+
+    typedef int M_t;
+    typedef double F_t;
+
     inline int operator*(v_t a) const {
       if (x) {
         return -(int)a;
@@ -45,6 +59,12 @@ class ising_t {
     }
 };
 
+/* using measure.hpp additionally requires a norm_squared function which takes
+ * an F_t to a double, and a write_magnetization function, which takes an M_t
+ * and a FILE pointer and appropriately records the contents of the former to
+ * the latter.
+ */
+
 double norm_squared(double s) {
   return pow(s, 2);
 }
@@ -53,6 +73,11 @@ void write_magnetization(int M, FILE *outfile) {
   fwrite(&M, sizeof(int), 1, outfile);
 }
 
+/* these definitions allow wolff/finite_states.hpp to be invoked and provide
+ * much faster performance for models whose number of possible spin
+ * configurations is finite.
+ */
+
 #define N_STATES 2
 const ising_t states[2] = {ising_t(0), ising_t(1)};
 q_t state_to_ind(ising_t state) { return (q_t)state.x; }
diff --git a/examples/src/models/ising/wolff_ising.cpp b/examples/src/models/ising/wolff_ising.cpp
index f7a3ada..de04f32 100644
--- a/examples/src/models/ising/wolff_ising.cpp
+++ b/examples/src/models/ising/wolff_ising.cpp
@@ -127,14 +127,12 @@ int main(int argc, char *argv[]) {
     return z2_t(true);
   };
 
-  FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
-
-  std::function <void(const state_t<z2_t, ising_t>&)> other_f;
+  std::function <void(const state_t<z2_t, ising_t>&, const wolff_research_measurements<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;
+    other_f = [&] (const state_t<z2_t, ising_t>& s, const wolff_research_measurements<z2_t, ising_t>& meas) {
+      sum_of_clusterSize += meas.last_cluster_size;
     };
   } else if (draw) {
 #ifdef HAVE_GLUT
@@ -148,7 +146,7 @@ int main(int argc, char *argv[]) {
     glLoadIdentity();
     gluOrtho2D(0.0, L, 0.0, L);
 
-    other_f = [] (const state_t <z2_t, ising_t>& s) {
+    other_f = [] (const state_t <z2_t, ising_t>& s, const wolff_research_measurements<z2_t, ising_t>& meas) {
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(s.L, 2); i++) {
 #ifdef NOFIELD
@@ -166,10 +164,10 @@ int main(int argc, char *argv[]) {
     };
 #endif
   } else {
-    other_f = [] (const state_t<z2_t, ising_t>& s) {};
+    other_f = [] (const state_t<z2_t, ising_t>& s, const wolff_research_measurements<z2_t, ising_t>& meas) {};
   }
 
-  std::function <void(const state_t<z2_t, ising_t>&)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f);
+  wolff_research_measurements<z2_t, ising_t> m(measurement_flags, timestamp, other_f, s, silent);
 
   // add line to metadata file with run info
   {
@@ -185,18 +183,15 @@ int main(int argc, char *argv[]) {
     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, rng, silent);
+      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, m.E, m.last_cluster_size);
+      wolff(N, s, gen_R, m, rng);
       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);
+    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, m.E, m.last_cluster_size);
   } else {
-    wolff(N, s, gen_R, measurements, rng, silent);
+    wolff(N, s, gen_R, m, rng);
   }
 
-  measure_free_files(measurement_flags, outfiles);
-
   return 0;
-
 }
 
diff --git a/examples/src/models/potts/symmetric.hpp b/examples/src/models/potts/symmetric.hpp
index 8636f15..bc8673f 100644
--- a/examples/src/models/potts/symmetric.hpp
+++ b/examples/src/models/potts/symmetric.hpp
@@ -36,6 +36,7 @@ class symmetric_t : public std::array<q_t, q> {
         }
       }
 
+      printf("Your spin wasn't a valid state!", s.x);
       exit(EXIT_FAILURE);
     }
 
diff --git a/examples/src/models/potts/wolff_clock.cpp b/examples/src/models/potts/wolff_clock.cpp
index 020415d..0706cc5 100644
--- a/examples/src/models/potts/wolff_clock.cpp
+++ b/examples/src/models/potts/wolff_clock.cpp
@@ -9,6 +9,7 @@
 #include "dihedral.hpp"
 #include "potts.hpp"
 #include <colors.h>
+#include <measure.hpp>
 
 // hack to speed things up considerably
 #define N_STATES POTTSQ
@@ -95,7 +96,7 @@ int main(int argc, char *argv[]) {
   std::function <dihedral_t<POTTSQ>(std::mt19937&, potts_t<POTTSQ>)> gen_R = [] (std::mt19937& r, potts_t<POTTSQ> v) -> dihedral_t<POTTSQ> {
     dihedral_t<POTTSQ> rot;
     rot.is_reflection = true;
-    std::uniform_int_distribution<q_t> dist(0, POTTSQ - 1);
+    std::uniform_int_distribution<q_t> dist(0, POTTSQ - 2);
     q_t x = dist(r);
     rot.x = (2 * v.x + x + 1) % POTTSQ;
 
@@ -103,13 +104,11 @@ int main(int argc, char *argv[]) {
   };
 
   // define function that updates any number of measurements
-  std::function <void(const sim_t&)> measurement;
+  std::function <void(const sim_t&, const wolff_research_measurements<dihedral_t<POTTSQ>, potts_t<POTTSQ>>&)> measurement;
 
-  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;
+    measurement = [&] (const sim_t& s, const wolff_research_measurements<dihedral_t<POTTSQ>, potts_t<POTTSQ>>&) {
     };
   } else {
     // a more complex example: measure the average magnetization, and draw the spin configuration to the screen
@@ -125,8 +124,7 @@ 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;
+    measurement = [&] (const sim_t& s, const wolff_research_measurements<dihedral_t<POTTSQ>, potts_t<POTTSQ>>&) {
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(L, 2); i++) {
         potts_t<POTTSQ> tmp_s = s.R.act_inverse(s.spins[i]);
@@ -138,17 +136,13 @@ int main(int argc, char *argv[]) {
 #endif
   }
 
-  // run wolff for N cluster flips
-  wolff(N, s, gen_R, measurement, rng, silent);
+  wolff_research_measurements<dihedral_t<POTTSQ>, potts_t<POTTSQ>> m(0, 0, measurement, s, silent);
 
-  // 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);
+  // run wolff for N cluster flips
+  wolff(N, s, gen_R, m, rng);
 
   // free the random number generator
 
-  if (draw) {
-  }
-
   return 0;
 
 }
diff --git a/examples/src/models/potts/wolff_potts.cpp b/examples/src/models/potts/wolff_potts.cpp
index a1e9284..7b92ac1 100644
--- a/examples/src/models/potts/wolff_potts.cpp
+++ b/examples/src/models/potts/wolff_potts.cpp
@@ -119,7 +119,7 @@ int main(int argc, char *argv[]) {
   std::function <symmetric_t<POTTSQ>(std::mt19937&, potts_t<POTTSQ>)> gen_R = [] (std::mt19937& r, potts_t<POTTSQ> v) -> symmetric_t<POTTSQ> {
     symmetric_t<POTTSQ> rot;
 
-    std::uniform_int_distribution<q_t> dist(0, POTTSQ - 1);
+    std::uniform_int_distribution<q_t> dist(0, POTTSQ - 2);
     q_t j = dist(r);
     q_t swap_v;
     if (j < v.x) {
@@ -134,14 +134,12 @@ int main(int argc, char *argv[]) {
     return rot;
   };
 
-  FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
-
-  std::function <void(const sim_t&)> other_f;
+  std::function <void(const sim_t&, const wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>>&)> 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;
+    other_f = [&] (const sim_t& s, const wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>>& m) {
+      sum_of_clusterSize += m.last_cluster_size;
     };
   } else if (draw) {
 #ifdef HAVE_GLUT
@@ -155,7 +153,7 @@ int main(int argc, char *argv[]) {
     glLoadIdentity();
     gluOrtho2D(0.0, L, 0.0, L);
 
-    other_f = [] (const sim_t& s) {
+    other_f = [] (const sim_t& s, const wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>>& m) {
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(s.L, 2); i++) {
         potts_t<POTTSQ> tmp_s = s.R.act_inverse(s.spins[i]);
@@ -166,10 +164,10 @@ int main(int argc, char *argv[]) {
     };
 #endif
   } else {
-    other_f = [] (const sim_t& s) {};
+    other_f = [] (const sim_t& s, const wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>>& m) {};
   }
 
-  std::function <void(const sim_t&)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f);
+  wolff_research_measurements<symmetric_t<POTTSQ>, potts_t<POTTSQ>> m(measurement_flags, timestamp, other_f, s, silent);
 
   // add line to metadata file with run info
   {
@@ -194,18 +192,17 @@ int main(int argc, char *argv[]) {
     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, rng, silent);
+      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, m.E, m.last_cluster_size);
+      wolff(N, s, gen_R, m, rng);
       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);
+    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, m.E, m.last_cluster_size);
   } else {
-    wolff(N, s, gen_R, measurements, rng, silent);
+    wolff(N, s, gen_R, m, rng);
   }
 
   // free the random number generator
   free(H_vec);
-  measure_free_files(measurement_flags, outfiles);
 
   return 0;