- added support for computing spatial fourier transforms

 - measurements now are functions passed to wolff in array

9 files changed, 217 insertions, 50 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index 5f770ad..8061d34 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -52,7 +52,7 @@ template <class T>
 T subtract(T, T);
 
 template <class R_t, class X_t>
-v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
+void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
   v_t nv = 0;
 
   ll_t *stack = NULL;     // create a new stack
@@ -104,8 +104,8 @@ v_t 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);
+          subtract (&(state->M), rs_old);
+          add (&(state->M), rs_new);
           state->E += dE;
 
           free_spin (rs_old);
@@ -137,6 +137,6 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
 
   free(marks);
 
-  return nv;
+  state->last_cluster_size = nv;
 }
 
diff --git a/lib/correlation.h b/lib/correlation.h
new file mode 100644
index 0000000..5722aab
--- /dev/null
+++ b/lib/correlation.h
@@ -0,0 +1,36 @@
+
+#pragma once
+
+#include "types.h"
+#include "state.h"
+
+#include <fftw3.h>
+
+template <class R_t, class X_t>
+double correlation_length(const state_t <R_t, X_t> *s) {
+  double *data = (double *)fftw_malloc(s->nv * sizeof(double));
+  int rank = s->D;
+  int *n = (int *)malloc(rank * sizeof(int));
+  fftw_r2r_kind *kind = (fftw_r2r_kind *)malloc(rank * sizeof(fftw_r2r_kind));
+  for (D_t i = 0; i < rank; i++) {
+    n[i] = s->L;
+    kind[i] = FFTW_R2HC;
+  }
+  fftw_plan plan = fftw_plan_r2r(rank, n, data, data, kind, 0);
+
+  for (v_t i = 0; i < s->nv; i++) {
+    data[i] = correlation_component(s->spins[i]);
+  }
+
+  fftw_execute(plan);
+
+  double length = pow(data[0], 2);
+
+  fftw_destroy_plan(plan);
+  fftw_free(data);
+  free(n);
+  free(kind);
+
+  return length;
+}
+
diff --git a/lib/height.h b/lib/height.h
new file mode 100644
index 0000000..16c4dc5
--- /dev/null
+++ b/lib/height.h
@@ -0,0 +1,62 @@
+
+#pragma once
+
+#include <cmath>
+
+#include "types.h"
+
+// object definition
+template <class T>
+struct height_t { T x; };
+
+// init, copy, add, subtract, scalar_multiple, free_spin, and
+// write_magnetization are necessary for the operation of wolff.h
+template <class T>
+void init(height_t *ptr) {
+  ptr->x = (T)0;
+}
+
+template <class T>
+height_t <T> copy (height_t h) {
+  return h;
+}
+
+template <class T>
+void add (height_t <T> *h1, height_t <T> h2) {
+  h1->x += h2.x;
+}
+
+template <class T>
+void subtract (height_t <T> *h1, height_T <T> h2) {
+  h1->x -= h2.x;
+}
+
+template <class T>
+height_t <T> scalar_multiple(v_t a, height_t <T> h) {
+  height_t <T> hm;
+  hm.x = a * h.x;
+
+  return hm;
+}
+
+template <class T>
+void free_spin (height_t <T> h) {
+}
+
+template <class T>
+void write_magnetization(height_t <T> M, FILE *outfile) {
+  fwrite(&(M.x), sizeof(T), 1, outfile);
+}
+
+template <class T>
+double correlation_component(height_t <T> h) {
+  return (double)h.x;
+}
+
+// below here are not necessary for operation
+
+template <class T>
+T dot(height_t <T> h1, height_t <T> h2) {
+  return (h1.x - h2.x) * (h1.x - h2.x);
+}
+
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index d8ad33d..340ee2c 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -184,9 +184,10 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const state
   }
 
   double v2 = 0;
+  double factor = gsl_ran_ugaussian(r);
 
   for (q_t i = 0; i < q; i++) {
-    tmp_v.x[i] = (tmp_v.x[i] - tmpM * s->M.x[i] / M2) + epsilon * gsl_ran_ugaussian(r);
+    tmp_v.x[i] = (tmp_v.x[i] - tmpM * s->M.x[i] / M2) + epsilon * factor * s->M.x[i] / sqrt(M2);
     v2 += tmp_v.x[i] * tmp_v.x[i];
   }
 
diff --git a/lib/state.h b/lib/state.h
index 8d4265b..1491938 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -19,6 +19,7 @@ class state_t {
     R_t R;
     double E;
     X_t M; // the "sum" of the spins, like the total magnetization
+    v_t last_cluster_size;
 
     std::function <double(X_t, X_t)> J;
     std::function <double(X_t)> H;
@@ -36,6 +37,7 @@ class state_t {
       init (&R);
       E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
       M = scalar_multiple (nv, spins[0]);
+      last_cluster_size = 0;
     }
 
     ~state_t() {
diff --git a/lib/vector.h b/lib/vector.h
index 62ce59e..88934c1 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -30,16 +30,16 @@ vector_t <q, T> copy (vector_t <q, T> v) {
 }
 
 template <q_t q, class T>
-void add (vector_t <q, T> v1, vector_t <q, T> v2) {
+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];
+    v1->x[i] += v2.x[i];
   }
 }
 
 template <q_t q, class T>
-void subtract (vector_t <q, T> v1, vector_t <q, T> v2) {
+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];
+    v1->x[i] -= v2.x[i];
   }
 }
 
@@ -72,7 +72,7 @@ void free_spin (vector_t <q, T> v) {
 
 template <q_t q, class T>
 void write_magnetization(vector_t <q, T> M, FILE *outfile) {
-  fwrite(M.x, sizeof(double), q, outfile);
+  fwrite(M.x, sizeof(T), q, outfile);
 }
 
 template <q_t q> // save some space and don't write whole doubles
@@ -83,3 +83,8 @@ 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];
+}
+
diff --git a/lib/wolff.h b/lib/wolff.h
index 90470d5..2d7e8b2 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -1,7 +1,4 @@
 
-#include <time.h>
-#include <getopt.h>
-
 #include "cluster.h"
 #include "state.h"
 
@@ -13,7 +10,7 @@ double H_vector(vector_t <q, T> v1, T *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 long timestamp, bool silent) {
+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);
 
@@ -21,52 +18,26 @@ void wolff(count_t N, D_t D, L_t L, double T, std::function <double(X_t, X_t)> J
   gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
   gsl_rng_set(r, rand_seed());
 
-  char *filename_M = (char *)malloc(255 * sizeof(char));
-  char *filename_E = (char *)malloc(255 * sizeof(char));
-  char *filename_S = (char *)malloc(255 * sizeof(char));
-
-  sprintf(filename_M, "wolff_%lu_M.dat", timestamp);
-  sprintf(filename_E, "wolff_%lu_E.dat", timestamp);
-  sprintf(filename_S, "wolff_%lu_S.dat", timestamp);
-
-  FILE *outfile_M = fopen(filename_M, "wb");
-  FILE *outfile_E = fopen(filename_E, "wb");
-  FILE *outfile_S = fopen(filename_S, "wb");
-
-  free(filename_M);
-  free(filename_E);
-  free(filename_S);
-
-  v_t cluster_size = 0;
-
   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, cluster_size);
+    if (!silent) printf("\033[F\033[JWOLFF: sweep %" 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);
-
-    cluster_size = flip_cluster <R_t, X_t> (&s, v0, step, r);
-
+    flip_cluster <R_t, X_t> (&s, v0, step, r);
     free_spin(step);
 
-    {
-      float smaller_E = (float)s.E;
-      fwrite(&smaller_E, sizeof(float), 1, outfile_E);
+    for (unsigned int i = 0; i < n_measurements; i++) {
+      measurements[i](&s);
     }
-    write_magnetization(s.M, outfile_M);
-    fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S);
 
   }
   if (!silent) {
     printf("\033[F\033[J");
   }
-  printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, cluster_size);
+  printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size);
 
-  fclose(outfile_M);
-  fclose(outfile_E);
-  fclose(outfile_S);
+  fftw_cleanup();
 
   gsl_rng_free(r);
 }
diff --git a/src/wolff_heisenberg.cpp b/src/wolff_heisenberg.cpp
index 83ea503..f38fa96 100644
--- a/src/wolff_heisenberg.cpp
+++ b/src/wolff_heisenberg.cpp
@@ -1,6 +1,7 @@
 
 #include <getopt.h>
 
+#include <correlation.h>
 #include <wolff.h>
 
 int main(int argc, char *argv[]) {
@@ -87,8 +88,51 @@ int main(int argc, char *argv[]) {
 
   fclose(outfile_info);
 
-
-  wolff <orthogonal_t <3, double>, vector_t <3, double>> (N, D, L, T, dot <3, double>, std::bind(H_vector <3, double>, std::placeholders::_1, H), gen_R, timestamp, silent);
+  char *filename_M = (char *)malloc(255 * sizeof(char));
+  char *filename_E = (char *)malloc(255 * sizeof(char));
+  char *filename_S = (char *)malloc(255 * sizeof(char));
+  char *filename_X = (char *)malloc(255 * sizeof(char));
+
+  sprintf(filename_M, "wolff_%lu_M.dat", timestamp);
+  sprintf(filename_E, "wolff_%lu_E.dat", timestamp);
+  sprintf(filename_S, "wolff_%lu_S.dat", timestamp);
+  sprintf(filename_X, "wolff_%lu_X.dat", timestamp);
+
+  FILE *outfile_M = fopen(filename_M, "wb");
+  FILE *outfile_E = fopen(filename_E, "wb");
+  FILE *outfile_S = fopen(filename_S, "wb");
+  FILE *outfile_X = fopen(filename_X, "wb");
+
+  free(filename_M);
+  free(filename_E);
+  free(filename_S);
+  free(filename_X);
+
+  std::function <void(const state_t <orthogonal_t <3, double>, vector_t <3, double>> *)> *measurements = (std::function <void(const state_t <orthogonal_t <3, double>, vector_t <3, double>> *)> *)malloc(4 * sizeof(std::function <void(const state_t <orthogonal_t <3, double>, vector_t <3, double>> *)>));
+
+  measurements[0] = [&](const state_t <orthogonal_t <3, double>, vector_t <3, double>> *s) {
+    float smaller_E = (float)s->E;
+    fwrite(&smaller_E, sizeof(float), 1, outfile_E);
+  };
+  measurements[1] = [&](const state_t <orthogonal_t <3, double>, vector_t <3, double>> *s) {
+    float smaller_X = (float)correlation_length(s);
+    fwrite(&smaller_X, sizeof(float), 1, outfile_X);
+  };
+  measurements[2] = [&](const state_t <orthogonal_t <3, double>, vector_t <3, double>> *s) {
+    write_magnetization(s->M, outfile_M);
+  };
+  measurements[3] = [&](const state_t <orthogonal_t <3, double>, vector_t <3, double>> *s) {
+    fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, outfile_S);
+  };
+
+  wolff <orthogonal_t <3, double>, vector_t <3, double>> (N, D, L, T, dot <3, double>, std::bind(H_vector <3, double>, std::placeholders::_1, H), gen_R, 4, measurements, silent);
+
+  free(measurements);
+
+  fclose(outfile_M);
+  fclose(outfile_E);
+  fclose(outfile_S);
+  fclose(outfile_X);
 
   free(H);
 
diff --git a/src/wolff_planar.cpp b/src/wolff_planar.cpp
index 2a3d02b..4f26471 100644
--- a/src/wolff_planar.cpp
+++ b/src/wolff_planar.cpp
@@ -2,6 +2,7 @@
 #include <getopt.h>
 
 #include <wolff.h>
+#include <correlation.h>
 
 int main(int argc, char *argv[]) {
 
@@ -88,8 +89,53 @@ int main(int argc, char *argv[]) {
 
   fclose(outfile_info);
 
-
-  wolff <orthogonal_t <2, double>, vector_t <2, double>> (N, D, L, T, dot <2, double>, std::bind(H_vector <2, double>, std::placeholders::_1, H), gen_R, timestamp, silent);
+  char *filename_M = (char *)malloc(255 * sizeof(char));
+  char *filename_E = (char *)malloc(255 * sizeof(char));
+  char *filename_S = (char *)malloc(255 * sizeof(char));
+  char *filename_X = (char *)malloc(255 * sizeof(char));
+
+  sprintf(filename_M, "wolff_%lu_M.dat", timestamp);
+  sprintf(filename_E, "wolff_%lu_E.dat", timestamp);
+  sprintf(filename_S, "wolff_%lu_S.dat", timestamp);
+  sprintf(filename_X, "wolff_%lu_X.dat", timestamp);
+
+  FILE *outfile_M = fopen(filename_M, "wb");
+  FILE *outfile_E = fopen(filename_E, "wb");
+  FILE *outfile_S = fopen(filename_S, "wb");
+  FILE *outfile_X = fopen(filename_X, "wb");
+
+  free(filename_M);
+  free(filename_E);
+  free(filename_S);
+  free(filename_X);
+
+  std::function <void(const state_t <orthogonal_t <2, double>, vector_t <2, double>> *)> *measurements = (std::function <void(const state_t <orthogonal_t <2, double>, vector_t <2, double>> *)> *)calloc(4, sizeof(std::function <void(const state_t <orthogonal_t <2, double>, vector_t <2, double>> *)>));
+
+  measurements[0] = (std::function <void(const state_t <orthogonal_t <2, double>, vector_t <2, double>> *)>)[&](const state_t <orthogonal_t <2, double>, vector_t <2, double>> *s) {
+    float smaller_E = (float)s->E;
+    fwrite(&smaller_E, sizeof(float), 1, outfile_E);
+  };
+  measurements[1] = [&](const state_t <orthogonal_t <2, double>, vector_t <2, double>> *s) {
+    float smaller_X = (float)correlation_length(s);
+    fwrite(&smaller_X, sizeof(float), 1, outfile_X);
+  };
+  measurements[2] = [&](const state_t <orthogonal_t <2, double>, vector_t <2, double>> *s) {
+    write_magnetization(s->M, outfile_M);
+  };
+  measurements[3] = [&](const state_t <orthogonal_t <2, double>, vector_t <2, double>> *s) {
+    fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, outfile_S);
+  };
+
+
+
+  wolff <orthogonal_t <2, double>, vector_t <2, double>> (N, D, L, T, dot <2, double>, std::bind(H_vector <2, double>, std::placeholders::_1, H), gen_R, 4, measurements, silent);
+
+  free(measurements);
+
+  fclose(outfile_M);
+  fclose(outfile_E);
+  fclose(outfile_S);
+  fclose(outfile_X);
 
   free(H);