added ising example to cpp collection

9 files changed, 243 insertions, 228 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index 2225e2f..4b3c2f6 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -23,6 +23,8 @@
 #include "measurement.h"
 #include "vector.h"
 #include "orthogonal.h"
+#include "ising.h"
+#include "z2.h"
 #include "dihedral.h"
 #include "dihinf.h"
 #include "yule_walker.h"
@@ -107,17 +109,17 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
           double dE = state->H(rs_old) - state->H(rs_new);
           prob = 1.0 - exp(-dE / state->T);
 
-          subtract (&(state->M), rs_old);
-          add (&(state->M), rs_new);
+          add(&(state->M), -1, rs_old);
+          add(&(state->M),  1, rs_new);
           state->E += dE;
 
           for (D_t i = 0; i < state->D; i++) {
             double x = (double)((non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L) / (double)state->L; 
-            scalar_subtract(&(state->ReF[i]), cos(2 * M_PI * x), rs_old);
-            scalar_add(&(state->ReF[i]), cos(2 * M_PI * x), rs_new);
+            add(&(state->ReF[i]), -cos(2 * M_PI * x), rs_old);
+            add(&(state->ReF[i]),  cos(2 * M_PI * x), rs_new);
 
-            scalar_subtract(&(state->ImF[i]), sin(2 * M_PI * x), rs_old);
-            scalar_add(&(state->ImF[i]), sin(2 * M_PI * x), rs_new);
+            add(&(state->ImF[i]), -sin(2 * M_PI * x), rs_old);
+            add(&(state->ImF[i]),  sin(2 * M_PI * x), rs_new);
           }
 
           free_spin (rs_old);
diff --git a/lib/ising.h b/lib/ising.h
new file mode 100644
index 0000000..4e5164b
--- /dev/null
+++ b/lib/ising.h
@@ -0,0 +1,80 @@
+#pragma once
+
+#include <cmath>
+#include <stdlib.h>
+
+#include "types.h"
+
+class ising_t {
+  public:
+    bool x;
+
+    typedef int M_t;
+    typedef double F_t;
+};
+
+void init(ising_t *p) {
+  p->x = false;
+}
+
+void free_spin(ising_t s) {
+  // do nothing!
+}
+
+void free_spin(int s) {
+  // do nothing
+}
+
+void free_spin(double s) {
+  // do nothing
+}
+
+ising_t copy(ising_t s) {
+  return s;
+}
+
+template <class T>
+void add(T *s1, T a, ising_t s2) {
+  if (s2.x) {
+    *s1 -= a;
+  } else {
+    *s1 += a;
+  }
+}
+
+int scalar_multiple(int factor, ising_t s) {
+  if (s.x) {
+    return -factor;
+  } else {
+    return factor;
+  }
+}
+
+
+double norm_squared(double s) {
+  return pow(s, 2);
+}
+
+template <class T>
+void write_magnetization(T M, FILE *outfile) {
+  fwrite(&M, sizeof(T), 1, outfile);
+}
+
+// below this line is unnecessary, but convenient
+
+double ising_dot(ising_t s1, ising_t s2) {
+  if (s1.x == s2.x) {
+    return 1.0;
+  } else {
+    return -1.0;
+  }
+}
+
+double scalar_field(ising_t s, double H) {
+  if (s.x) {
+    return -H;
+  } else {
+    return H;
+  }
+}
+
diff --git a/lib/measure.h b/lib/measure.h
index d20c081..52e43af 100644
--- a/lib/measure.h
+++ b/lib/measure.h
@@ -1,7 +1,9 @@
 
 #pragma once
 
-#define POSSIBLE_MEASUREMENTS 4
+#include "measurement.h"
+
+#define POSSIBLE_MEASUREMENTS 5
 const unsigned char measurement_energy        = 1 << 0;
 const unsigned char measurement_clusterSize   = 1 << 1;
 const unsigned char measurement_magnetization = 1 << 2;
@@ -43,6 +45,13 @@ std::function <void(const state_t <R_t, X_t> *)> measurement_fourier_file(FILE *
   };
 }
 
+template <class R_t, class X_t>
+std::function <void(const state_t <R_t, X_t> *)> measurement_average_cluster(meas_t *x) {
+  return [=](const state_t <R_t, X_t> *s) {
+    meas_update(x, s->last_cluster_size);
+  };
+}
+
 #endif
 
 
diff --git a/lib/measurement.c b/lib/measurement.c
index b30cf6b..ad824f6 100644
--- a/lib/measurement.c
+++ b/lib/measurement.c
@@ -1,15 +1,6 @@
 
-#include "convex.h"
 #include "measurement.h"
 
-meas_t *meas_initialize(count_t W) {
-  meas_t *m = (meas_t *)calloc(1, sizeof(meas_t));
-  m->W = W;
-  m->xx = (double *)calloc(2 * W + 1, sizeof(double));
-
-  return m;
-}
-
 double add_to_avg(double mx, double x, count_t n) {
   return mx * (n / (n + 1.0)) + x / (n + 1.0);
 }
@@ -19,42 +10,24 @@ void meas_update(meas_t *m, double x) {
 
   m->x = add_to_avg(m->x, x, n);
   m->x2 = add_to_avg(m->x2, pow(x, 2), n);
-  m->x4 = add_to_avg(m->x4, pow(x, 4), 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);
 
-  dll_t *tmp_window = m->x_window;
-  dll_t *pos_save;
-  count_t t = 0;
-
-  while (tmp_window != NULL) {
-    m->xx[t] = add_to_avg(m->xx[t], x * (tmp_window->x), m->n - t - 1);
-    t++;
-    if (t == 2 * (m->W)) {
-      pos_save = tmp_window;
-    }
-    tmp_window = tmp_window->next;
+  /*
+  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);
   }
-
-  if (t == 2 * (m->W) + 1) {
-    if (2 * (m->W) + 1 == 1) {
-      free(m->x_window);
-      m->x_window = NULL;
-    } else {
-      free(pos_save->next);
-      pos_save->next = NULL;
-    }
-  }
-
-  stack_push_d(&(m->x_window), x);
+  */
 
   (m->n)++;
 }
 
 double meas_dx(const meas_t *m) {
-  return 2 * get_tau(m) * Cxx(m, 0) / m->n;
-//  return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2)));
+  return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2)));
 }
 
 double meas_c(const meas_t *m) {
@@ -101,105 +74,3 @@ double rho(const autocorr_t *o, count_t i) {
   return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2));
 }
 
-double Cxx(const meas_t *m, count_t t) {
-  return m->xx[t] - pow(m->x, 2);
-}
-
-double rho_m(const meas_t *m, count_t t) {
-  return Cxx(m, t) / Cxx(m, 0);
-}
-
-double get_tau(const meas_t *m) {
-  double *Gammas = (double *)malloc((m->W + 1) * sizeof(double));
-
-  Gammas[0] = 1 + rho_m(m, 0);
-  for (uint64_t i = 0; i < m->W; i++) {
-    Gammas[1 + i] = rho_m(m, 2 * i + 1) + rho_m(m, 2 * i + 2);
-  }
-
-  uint64_t n;
-  for (n = 0; n < m->W + 1; n++) {
-    if (Gammas[n] <= 0) {
-      break;
-    }
-  }
-
-  double *conv_Gamma = get_convex_minorant(n, Gammas);
-
-  double tau = - 0.5;
-
-  for (uint64_t i = 0; i < n + 1; i++) {
-    tau += conv_Gamma[i];
-  }
-
-  free(Gammas);
-
-  return tau;
-}
-
-void print_meas(const meas_t *m, const char *sym, FILE *outfile) {
-  fprintf(outfile, "%s-><|n->%" PRIcount ",x->%.15f,x^2->%.15f,x^4->%.15f,xx->{", sym, m->n, m->x, m->x2, m->x4);
-  for (count_t i = 0; i < 2 * (m->W) + 1; i++) {
-    fprintf(outfile, "%.15f", m->xx[i]);
-    if (i < 2 * (m->W)) {
-      fprintf(outfile, ",");
-    }
-  }
-  fprintf(outfile, "}|>");
-}
-
-void print_vec_meas(q_t q, const meas_t **m, const char *sym, FILE *outfile) {
-  fprintf(outfile, "%s-><|n->{", sym);
-  for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%" PRIcount, m[i]->n);
-    if (i < q - 1) {
-      fprintf(outfile, ",");
-    }
-  }
-  fprintf(outfile, "},x->{");
-  for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", m[i]->x);
-    if (i < q - 1) {
-      fprintf(outfile, ",");
-    }
-  }
-  fprintf(outfile, "},x^2->{");
-  for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", m[i]->x2);
-    if (i < q - 1) {
-      fprintf(outfile, ",");
-    }
-  }
-  fprintf(outfile, "},x^4->{");
-  for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "%.15f", m[i]->x4);
-    if (i < q - 1) {
-      fprintf(outfile, ",");
-    }
-  }
-  fprintf(outfile, "},xx->{");
-  for (q_t i = 0; i < q; i++) {
-    fprintf(outfile, "{");
-    for (count_t j = 0; j < 2 * (m[i]->W) + 1; j++) {
-      fprintf(outfile, "%.15f", m[i]->xx[j]);
-      if (j < 2 * (m[i]->W)) {
-        fprintf(outfile, ",");
-      }
-    }
-    fprintf(outfile, "}");
-    if (i < q - 1) {
-      fprintf(outfile, ",");
-    }
-  }
-  fprintf(outfile, "}|>");
-}
-
-void free_meas(meas_t *m) {
-  free(m->xx);
-  while (m->x_window != NULL) {
-    stack_pop_d(&(m->x_window));
-  }
-  free(m);
-}
-
-
diff --git a/lib/measurement.h b/lib/measurement.h
index d9bd52e..78fa51b 100644
--- a/lib/measurement.h
+++ b/lib/measurement.h
@@ -3,21 +3,20 @@
 
 #include <math.h>
 #include <stdlib.h>
-#include <stdio.h>
 
 #include "types.h"
 #include "stack.h"
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 typedef struct {
-  count_t n;
+  uint64_t n;
   double x;
   double x2;
-  double x4;
   double m2;
   double m4;
-  count_t W;
-  double *xx;
-  dll_t *x_window;
 } meas_t;
 
 typedef struct {
@@ -41,14 +40,6 @@ void update_autocorr(autocorr_t *OO, double O);
 
 double rho(const autocorr_t *o, uint64_t i);
 
-void print_meas(const meas_t *m, const char *sym, FILE *outfile);
-void print_vec_meas(q_t q, const meas_t **m, const char *sym, FILE *outfile);
-
-void free_meas(meas_t *m);
-
-meas_t *meas_initialize(count_t W);
-
-double get_tau(const meas_t *m);
-
-double Cxx(const meas_t *m, count_t t);
-
+#ifdef __cplusplus
+}
+#endif
diff --git a/lib/state.h b/lib/state.h
index c70459f..b40b85c 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -5,6 +5,7 @@
 
 #include "types.h"
 #include "graph.h"
+#include "ising.h"
 
 template <class R_t, class X_t>
 class state_t {
@@ -18,10 +19,10 @@ class state_t {
     X_t *spins;
     R_t R;
     double E;
-    X_t M; // the "sum" of the spins, like the total magnetization
+    typename X_t::M_t M; // the "sum" of the spins, like the total magnetization
     v_t last_cluster_size;
-    X_t *ReF;
-    X_t *ImF;
+    typename X_t::F_t *ReF;
+    typename X_t::F_t *ImF;
 
     std::function <double(X_t, X_t)> J;
     std::function <double(X_t)> H;
@@ -38,10 +39,10 @@ class state_t {
       }
       init (&R);
       E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
-      M = scalar_multiple(nv, spins[0]);
+      M = scalar_multiple((int)nv, spins[0]);
       last_cluster_size = 0;
-      ReF = (X_t *)malloc(D * sizeof(X_t));
-      ImF = (X_t *)malloc(D * sizeof(X_t));
+      ReF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t));
+      ImF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t));
       for (D_t i = 0; i < D; i++) {
         ReF[i] = scalar_multiple(0, spins[0]);
         ImF[i] = scalar_multiple(0, spins[0]);
diff --git a/lib/vector.h b/lib/vector.h
index a1084e2..c478618 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -6,17 +6,49 @@
 
 #include "types.h"
 
+/* The following is the minimum definition of a spin class.
+ *
+ * The class must contain an M_t and an F_t for holding the sum of an
+ * integer number of spins and a double-weighted number of spins,
+ * respectively.
+ *
+ * void init(X_t *p);
+ * void free_spin(X_t p);
+ * X_t copy(X_t x);
+ * void add(M_t *x1, int factor, X_t x2);
+ * void add(F_t *x1, double factor, X_t x2);
+ * M_t scalar_multiple(int factor, X_t x);
+ * double norm_squared(F_t x);
+ * void write_magnetization(M_t M, FILE *outfile);
+ *
+ */
+
 template <q_t q, class T>
-struct vector_t { T *x; };
+class vector_t {
+  public: 
+    T *x;
+
+    // M_t needs to hold the sum of nv spins
+    typedef vector_t <q, T> M_t;
+
+    // F_t needs to hold the double-weighted sum of spins
+    typedef vector_t <q, double> F_t;
+};
 
 template <q_t q, class T>
 void init(vector_t <q, T> *ptr) {
   ptr->x = (T *)calloc(q, sizeof(T));
 
+  // initialize a unit vector
   ptr->x[0] = (T)1;
 }
 
 template <q_t q, class T>
+void free_spin (vector_t <q, T> v) {
+  free(v.x);
+}
+
+template <q_t q, class T>
 vector_t <q, T> copy (vector_t <q, T> v) {
   vector_t <q, T> v_copy;
  
@@ -29,46 +61,15 @@ vector_t <q, T> copy (vector_t <q, T> v) {
   return v_copy;
 }
 
-template <q_t q, class T>
-void add (vector_t <q, T> *v1, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] += v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-void scalar_add (vector_t <q, T> *v1, T a, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] += a * v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-void scalar_subtract (vector_t <q, T> *v1, T a, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] -= a * v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-void subtract (vector_t <q, T> *v1, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] -= v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-T norm_squared (vector_t <q, T> v) {
-  T tmp = 0;
+template <q_t q, class T, class U, class V>
+void add(vector_t<q, U> *v1, V a, vector_t <q, T> v2) {
   for (q_t i = 0; i < q; i++) {
-    tmp += v.x[i] * v.x[i];
+    v1->x[i] += (U)(a * v2.x[i]);
   }
-
-  return tmp;
 }
 
 template <q_t q, class T>
-vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
+vector_t <q, T> scalar_multiple(int a, vector_t <q, T> v) {
   vector_t <q, T> multiple;
   multiple.x = (T *)malloc(q * sizeof(T));
   for (q_t i = 0; i < q; i++) {
@@ -79,19 +80,14 @@ vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
 }
 
 template <q_t q, class T>
-T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
-  T prod = 0;
+double norm_squared (vector_t <q, T> v) {
+  double tmp = 0;
 
   for (q_t i = 0; i < q; i++) {
-    prod += v1.x[i] * v2.x[i];
+    tmp += pow(v.x[i], 2);
   }
 
-  return prod;
-}
-
-template <q_t q, class T>
-void free_spin (vector_t <q, T> v) {
-  free(v.x);
+  return tmp;
 }
 
 template <q_t q, class T>
@@ -99,6 +95,8 @@ void write_magnetization(vector_t <q, T> M, FILE *outfile) {
   fwrite(M.x, sizeof(T), q, outfile);
 }
 
+// below functions and definitions are unnecessary for wolff.h but useful.
+
 template <q_t q> // save some space and don't write whole doubles
 void write_magnetization(vector_t <q, double> M, FILE *outfile) {
   for (q_t i = 0; i < q; i++) {
@@ -108,8 +106,14 @@ void write_magnetization(vector_t <q, double> M, FILE *outfile) {
 }
 
 template <q_t q, class T>
-double correlation_component(vector_t <q, T> v) {
-  return (double)v.x[0];
+T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
+  T prod = 0;
+
+  for (q_t i = 0; i < q; i++) {
+    prod += v1.x[i] * v2.x[i];
+  }
+
+  return prod;
 }
 
 template <q_t q, class T>
diff --git a/lib/wolff.h b/lib/wolff.h
index eb357c7..21e88d9 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -3,33 +3,26 @@
 #include "state.h"
 
 template <class R_t, class X_t>
-void wolff(count_t N, D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H, std::function <R_t(gsl_rng *, const state_t <R_t, X_t> *)> gen_R, unsigned int n_measurements, std::function <void(const state_t <R_t, X_t> *)> *measurements, bool silent) {
-
-  state_t <R_t, X_t> s(D, L, T, J, H);
-
-  // initialize random number generator
-  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
-  gsl_rng_set(r, rand_seed());
+void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, const state_t <R_t, X_t> *)> gen_R, unsigned int n_measurements, std::function <void(const state_t <R_t, X_t> *)> *measurements, gsl_rng *r, bool silent) {
 
   if (!silent) printf("\n");
   for (count_t steps = 0; steps < N; steps++) {
-    if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, s.last_cluster_size);
+    if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s->E, s->last_cluster_size);
 
-    v_t v0 = gsl_rng_uniform_int(r, s.nv);
-    R_t step = gen_R(r, &s);
-    flip_cluster <R_t, X_t> (&s, v0, step, r);
+    v_t v0 = gsl_rng_uniform_int(r, s->nv);
+    R_t step = gen_R(r, s);
+    flip_cluster <R_t, X_t> (s, v0, step, r);
     free_spin(step);
 
     for (unsigned int i = 0; i < n_measurements; i++) {
-      measurements[i](&s);
+      measurements[i](s);
     }
   }
 
   if (!silent) {
     printf("\033[F\033[J");
   }
-  printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size);
+  printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s->E, s->last_cluster_size);
 
-  gsl_rng_free(r);
 }
 
diff --git a/lib/z2.h b/lib/z2.h
new file mode 100644
index 0000000..599a6a5
--- /dev/null
+++ b/lib/z2.h
@@ -0,0 +1,64 @@
+
+#pragma once
+
+#include <gsl/gsl_rng.h>
+#include "types.h"
+#include "ising.h"
+#include "state.h"
+
+struct z2_t { bool x; };
+
+void init(z2_t *p) {
+  p->x = false;
+}
+
+void free_spin(z2_t p) {
+  // do nothing!
+}
+
+z2_t copy(z2_t x) {
+  return x;
+}
+
+ising_t act(z2_t r, ising_t s) {
+  ising_t rs;
+
+  if (r.x) {
+    rs.x = !s.x;
+    return rs;
+  } else {
+    rs.x = s.x;
+    return rs;
+  }
+}
+
+z2_t act(z2_t r1, z2_t r2) {
+  z2_t r3;
+
+  if (r1.x) {
+    r3.x = !r2.x;
+    return r3;
+  } else {
+    r3.x = r2.x;
+    return r3;
+  }
+}
+
+ising_t act_inverse(z2_t r, ising_t s) {
+  return act(r, s);
+}
+
+z2_t act_inverse(z2_t r1, z2_t r2) {
+  return act(r1, r2);
+}
+
+// these are all functions necessary for wolff.h
+
+z2_t generate_ising_rotation(gsl_rng *r, const state_t <z2_t, ising_t> *s) {
+  z2_t rot;
+  rot.x = true;
+  return rot;
+}
+
+
+