added ising example to cpp collection

2 files changed, 196 insertions, 3 deletions

diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index d718440..76831ba 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -40,6 +40,7 @@ int main(int argc, char *argv[]) {
 
   bool silent = false;
   bool use_pert = false;
+  bool N_is_sweeps = false;
 
   bool modulated_field = false;
   int order = 2;
@@ -51,7 +52,7 @@ int main(int argc, char *argv[]) {
 
   unsigned char measurement_flags = measurement_energy | measurement_clusterSize;
 
-  while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:S")) != -1) {
     switch (opt) {
     case 'N': // number of steps
       N = (count_t)atof(optarg);
@@ -82,11 +83,14 @@ int main(int argc, char *argv[]) {
       modulated_field = true;
       break;
     case 'M':
-      measurement_flags |= 1 << atoi(optarg);
+      measurement_flags ^= 1 << atoi(optarg);
       break;
     case 'o':
       order = atoi(optarg);
       break;
+    case 'S':
+      N_is_sweeps = true;
+      break;
     default:
       exit(EXIT_FAILURE);
     }
@@ -179,6 +183,13 @@ int main(int argc, char *argv[]) {
     n_measurements++;
   }
 
+  meas_t *meas_sweeps;
+  if (N_is_sweeps) {
+    meas_sweeps = (meas_t *)calloc(1, sizeof(meas_t));
+    measurements[n_measurements] = measurement_average_cluster<orthogonal_R_t, vector_R_t> (meas_sweeps);
+    n_measurements++;
+  }
+
   std::function <double(vector_R_t)> H;
 
   if (modulated_field) {
@@ -187,7 +198,24 @@ int main(int argc, char *argv[]) {
     H = std::bind(H_vector <N_COMP, double>, std::placeholders::_1, H_vec);
   }
 
-  wolff <orthogonal_R_t, vector_R_t> (N, D, L, T, dot <N_COMP, double>, H, gen_R, n_measurements, measurements, silent);
+  // initialize random number generator
+  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+  gsl_rng_set(r, rand_seed());
+
+  state_t <orthogonal_R_t, vector_R_t> s(D, L, T, dot <N_COMP, double>, H);
+
+  if (N_is_sweeps) {
+    count_t N_rounds = 0;
+    printf("\n");
+    while (N_rounds * N * meas_sweeps->x < N * s.nv) {
+      printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+      wolff <orthogonal_R_t, vector_R_t> (N, &s, gen_R, n_measurements, 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)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+  } else {
+    wolff <orthogonal_R_t, vector_R_t> (N, &s, gen_R, n_measurements, measurements, r, silent);
+  }
 
   free(measurements);
 
@@ -204,7 +232,12 @@ int main(int argc, char *argv[]) {
     fclose(outfile_F);
   }
 
+  if (N_is_sweeps) {
+    free(meas_sweeps);
+  }
+
   free(H_vec);
+  gsl_rng_free(r);
 
   return 0;
 }
diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp
new file mode 100644
index 0000000..4d4c791
--- /dev/null
+++ b/src/wolff_ising.cpp
@@ -0,0 +1,160 @@
+
+#include <getopt.h>
+
+#include <wolff.h>
+#include <correlation.h>
+#include <measure.h>
+
+int main(int argc, char *argv[]) {
+
+  count_t N = (count_t)1e7;
+
+  D_t D = 2;
+  L_t L = 128;
+  double T = 2.26918531421;
+  double H = 0.0;
+
+  bool silent = false;
+  bool N_is_sweeps = false;
+
+  int opt;
+
+  unsigned char measurement_flags = measurement_energy | measurement_clusterSize;
+
+  while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:sM:S")) != -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. nth call couples to state n
+      H = atof(optarg);
+      break;
+    case 's': // don't print anything during simulation. speeds up slightly
+      silent = true;
+      break;
+    case 'M':
+      measurement_flags ^= 1 << atoi(optarg);
+      break;
+    case 'S':
+      N_is_sweeps = true;
+      break;
+    default:
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  unsigned long timestamp;
+
+  {
+    struct timespec spec;
+    clock_gettime(CLOCK_REALTIME, &spec);
+    timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec;
+  }
+
+  FILE *outfile_info = fopen("wolff_metadata.txt", "a");
+
+  fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"ISING\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> %.15f |>\n", timestamp, D, L, (v_t)pow(L, D), D * (v_t)pow(L, D), T, H);
+
+  fclose(outfile_info);
+
+  unsigned int n_measurements = 0;
+  std::function <void(const state_t <z2_t, ising_t> *)> *measurements = (std::function <void(const state_t <z2_t, ising_t> *)> *)calloc(POSSIBLE_MEASUREMENTS, sizeof(std::function <void(const state_t <z2_t, ising_t> *)>));
+  FILE *outfile_M, *outfile_E, *outfile_S, *outfile_F;
+
+  if (measurement_flags & measurement_energy) {
+    char *filename_E = (char *)malloc(255 * sizeof(char));
+    sprintf(filename_E, "wolff_%lu_E.dat", timestamp);
+    outfile_E = fopen(filename_E, "wb");
+    free(filename_E);
+    measurements[n_measurements] = measurement_energy_file<z2_t, ising_t> (outfile_E);
+    n_measurements++;
+  }
+
+  if (measurement_flags & measurement_clusterSize) {
+    char *filename_S = (char *)malloc(255 * sizeof(char));
+    sprintf(filename_S, "wolff_%lu_S.dat", timestamp);
+    outfile_S = fopen(filename_S, "wb");
+    free(filename_S);
+    measurements[n_measurements] = measurement_cluster_file<z2_t, ising_t> (outfile_S);
+    n_measurements++;
+  }
+
+  if (measurement_flags & measurement_magnetization) {
+    char *filename_M = (char *)malloc(255 * sizeof(char));
+    sprintf(filename_M, "wolff_%lu_M.dat", timestamp);
+    outfile_M = fopen(filename_M, "wb");
+    free(filename_M);
+    measurements[n_measurements] = measurement_magnetization_file<z2_t, ising_t> (outfile_M);
+    n_measurements++;
+  }
+
+  if (measurement_flags & measurement_fourierZero) {
+    char *filename_F = (char *)malloc(255 * sizeof(char));
+    sprintf(filename_F, "wolff_%lu_F.dat", timestamp);
+    outfile_F = fopen(filename_F, "wb");
+    free(filename_F);
+    measurements[n_measurements] = measurement_fourier_file<z2_t, ising_t> (outfile_F);
+    n_measurements++;
+  }
+
+  meas_t *meas_sweeps;
+  if (N_is_sweeps) {
+    meas_sweeps = (meas_t *)calloc(1, sizeof(meas_t));
+    measurements[n_measurements] = measurement_average_cluster<z2_t, ising_t> (meas_sweeps);
+    n_measurements++;
+  }
+
+  // initialize random number generator
+  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+  gsl_rng_set(r, rand_seed());
+
+  state_t <z2_t, ising_t> s(D, L, T, ising_dot, std::bind(scalar_field, std::placeholders::_1, H));
+
+  std::function <z2_t(gsl_rng *, const state_t <z2_t, ising_t> *)> gen = generate_ising_rotation;
+
+  if (N_is_sweeps) {
+    count_t N_rounds = 0;
+    printf("\n");
+    while (N_rounds * N * meas_sweeps->x < N * s.nv) {
+      printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+      wolff(N, &s, gen, n_measurements, 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)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+  } else {
+    wolff(N, &s, gen, n_measurements, measurements, r, silent);
+  }
+
+  free(measurements);
+
+  if (measurement_flags & measurement_energy) {
+    fclose(outfile_E);
+  }
+  if (measurement_flags & measurement_clusterSize) {
+    fclose(outfile_S);
+  }
+  if (measurement_flags & measurement_magnetization) {
+    fclose(outfile_M);
+  }
+  if (measurement_flags & measurement_fourierZero) {
+    fclose(outfile_F);
+  }
+
+  if (N_is_sweeps) {
+    free(meas_sweeps);
+  }
+
+  gsl_rng_free(r);
+
+  return 0;
+}
+