added new algorithm methods, specifically metropolis and wolff with a static ghost spin

3 files changed, 151 insertions, 15 deletions

diff --git a/lib/wolff.h b/lib/wolff.h
index 82ccb9e..bd49bf0 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -18,6 +18,12 @@
 
 #include "queue.h"
 
+typedef enum {
+  WOLFF,
+  WOLFF_GHOST,
+  METROPOLIS
+} sim_t;
+
 typedef struct {
   graph_t *g;
   bool *spins;
@@ -29,6 +35,7 @@ typedef struct {
   uint32_t nv;
   int32_t dJb;
   int32_t dHb;
+  bool hit_ghost;
 } cluster_t;
 
 typedef struct {
@@ -44,12 +51,12 @@ typedef struct {
 
 int32_t sign(double x);
 
-cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
+cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost,
                         gsl_rng *r);
 
 graph_t *graph_add_ext(const graph_t *g);
 
-uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
+uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
                     double *ps);
 
 void update_meas(meas_t *m, double x);
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index 59d04fc..bf78fda 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -2,6 +2,14 @@
 #include "queue.h"
 #include "wolff.h"
 
+int32_t spin_to_sign(bool spin) {
+  if (spin) {
+    return -1;
+  } else {
+    return 1;
+  }
+}
+
 int32_t sign(double x) {
   return x > 0 ? 1 : -1;
 }
@@ -50,7 +58,7 @@ graph_t *graph_add_ext(const graph_t *g) {
   return h;
 }
 
-cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
+cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost,
                         gsl_rng *r) {
   uint32_t v0;
   int32_t n_h_bonds, n_bonds;
@@ -73,7 +81,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
     v = stack_pop(&stack);
     nn = g->vei[v + 1] - g->vei[v];
 
-    if (x[v] == x0) { // if the vertex hasn't already been flipped
+    if (x[v] == x0 && !(c->hit_ghost && stop_on_ghost)) { // if the vertex hasn't already been flipped
       x[v] = !x[v];   // flip the vertex
 
       for (uint32_t i = 0; i < nn; i++) {
@@ -88,7 +96,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
 
         vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself
 
-        is_ext = (v1 == g->nv - 1 || v2 == g->nv - 1);
+        is_ext = (v1 == g->nv - 1 || v2 == g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph
 
         bond_counter = is_ext ? &(c->dHb) : &(c->dJb);
         prob = is_ext ? ps[1] : ps[0];
@@ -98,6 +106,9 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
           (*bond_counter)++;
 
           if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
+            if (is_ext && stop_on_ghost) {
+              c->hit_ghost = true;
+            }
             stack_push(&stack, vn); // push the neighboring vertex to the stack
           }
         } else {
@@ -114,8 +125,9 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
   return c;
 }
 
-uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
+uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
                     double *ps) {
+  uint32_t n_flips;
   bool no_r, no_ps;
   no_r = false;
   no_ps = false;
@@ -133,14 +145,74 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
     ps[1] = 1 - exp(-2 * fabs(H) / T);
   }
 
-  cluster_t *c = flip_cluster(s->g, ps, s->spins, r);
+  switch (sim) {
+    case METROPOLIS: {
+      uint32_t v0 = gsl_rng_uniform_int(r, s->g->nv); // pick a random vertex
+      uint32_t nn = s->g->vei[v0 + 1] - s->g->vei[v0];
+
+      double dE = 0;
+      for (uint32_t i = 0; i < nn; i++) {
+        bool is_ext;
+        uint32_t e, v1, v2, vn;
+        int32_t *bond_counter;
+        double prob;
+
+        e = s->g->ve[s->g->vei[v0] + i]; // select the ith bond connected to site
+        v1 = s->g->ev[2 * e];
+        v2 = s->g->ev[2 * e + 1];
 
-  s->M += - sign(H) * 2 * c->dHb;
-  s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
+        vn = v0 == v1 ? v2 : v1; // distinguish neighboring site from site itself
 
-  uint32_t n_flips = c->nv;
+        is_ext = (v1 == s->g->nv - 1 || v2 == s->g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph
 
-  free(c);
+        if (is_ext) {
+          dE += 2 * H * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]);
+        } else {
+          dE += 2 * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]);
+        }
+      }
+
+      double p = exp(-dE / T);
+      if (gsl_rng_uniform(r) < p) {
+        s->M += - sign(H) * 2 * spin_to_sign(s->spins[v0]) * spin_to_sign(s->spins[s->g->nv - 1]);
+        s->H += dE;
+        s->spins[v0] = !s->spins[v0];
+      }
+
+      n_flips = 1;
+                     }
+                     break;
+
+    case WOLFF: {
+      cluster_t *c = flip_cluster(s->g, ps, s->spins, false, r);
+      s->M += - sign(H) * 2 * c->dHb;
+      s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
+      n_flips = c->nv;
+
+      free(c);
+                }
+                break;
+    case WOLFF_GHOST: {
+      bool *spins_bak;
+      spins_bak = (bool *)malloc(s->g->nv * sizeof(bool));
+      memcpy(spins_bak, s->spins, s->g->nv * sizeof(bool));
+
+      cluster_t *c = flip_cluster(s->g, ps, s->spins, true, r);
+
+      if (c->hit_ghost) {
+        memcpy(s->spins, spins_bak, s->g->nv * sizeof(bool));
+      } else {
+        s->M += - sign(H) * 2 * c->dHb;
+        s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
+      }
+      free(spins_bak);
+
+      n_flips = c->nv;
+
+      free(c);
+                      }
+                break;
+  }
 
   if (no_ps) {
     free(ps);
diff --git a/src/wolff.c b/src/wolff.c
index da57057..5091c0c 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -3,25 +3,28 @@
 
 int main(int argc, char *argv[]) {
   int opt;
-  bool output_state, use_dual, M_stop;
+  sim_t sim;
+  bool output_state, use_dual, M_stop, record_autocorrelation;
   lattice_t lat;
   uint16_t L;
-  uint32_t min_runs, lattice_i;
+  uint32_t min_runs, lattice_i, sim_i;
   uint64_t N;
   double T, H, eps;
 
+  sim = WOLFF;
   L = 128;
   N = 1000;
   lat = SQUARE_LATTICE;
   use_dual = false;
   M_stop = false;
+  record_autocorrelation = false;
   T = 2.3;
   H = 0;
   eps = 1e30;
   output_state = false;
   min_runs = 10;
 
-  while ((opt = getopt(argc, argv, "N:L:T:H:m:e:oq:DM")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:H:m:e:oq:DMas:")) != -1) {
     switch (opt) {
     case 'N':
       N = (uint64_t)atof(optarg);
@@ -50,6 +53,27 @@ int main(int argc, char *argv[]) {
     case 'D':
       use_dual = true;
       break;
+    case 'a':
+      record_autocorrelation = true;
+      break;
+    case 's':
+      sim_i = atoi(optarg);
+      switch (sim_i) {
+      case 0:
+        sim = WOLFF;
+        break;
+      case 1:
+        sim = WOLFF_GHOST;
+        break;
+      case 2:
+        sim = METROPOLIS;
+        break;
+      default:
+        printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 "
+               "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
+        exit(EXIT_FAILURE);
+      }
+      break;
     case 'q':
       lattice_i = atoi(optarg);
       switch (lattice_i) {
@@ -97,8 +121,19 @@ int main(int argc, char *argv[]) {
 
   double diff = 1e31;
   uint64_t n_runs = 0;
+  uint64_t n_steps = 0;
   double clust_per_sweep = 0;
 
+  uint64_t auto_n = 10 * h->nv;
+  uint64_t *auto_time;
+  double *auto_energy;
+  if (record_autocorrelation) {
+    auto_time = (uint64_t *)malloc(auto_n * sizeof(uint64_t));
+    auto_energy = (double *)malloc(auto_n * sizeof(double));
+    auto_time[0] = 0;
+    auto_energy[0] = s->H / h->nv;
+  }
+
   meas_t *M, *aM, *eM, *mM, *E, *eE, *mE;
 
   M = calloc(1, sizeof(meas_t));
@@ -119,7 +154,7 @@ int main(int argc, char *argv[]) {
     uint32_t n_clust = 0;
 
     while (n_flips / h->nv < N) {
-      n_flips += wolff_step(T, H, s, r, bond_probs);
+      n_flips += wolff_step(T, H, s, sim, r, bond_probs);
       n_clust++;
     }
 
@@ -149,6 +184,17 @@ int main(int argc, char *argv[]) {
     clust_per_sweep = add_to_avg(clust_per_sweep, n_clust * 1. / N, n_runs);
 
     n_runs++;
+    n_steps += n_flips;
+
+    if (record_autocorrelation) {
+      if (n_runs == auto_n) {
+        auto_n *= 2;
+        auto_time = (uint64_t *)realloc(auto_time, auto_n * sizeof(uint64_t));
+        auto_energy = (double *)realloc(auto_energy, auto_n * sizeof(double));
+      }
+      auto_time[n_runs] = n_steps;
+      auto_energy[n_runs] = s->H / h->nv;
+    }
   }
 
   printf("\033[F\033[JWOLFF: sweep %" PRIu64
@@ -175,6 +221,17 @@ int main(int argc, char *argv[]) {
     fclose(state_file);
   }
 
+  if (record_autocorrelation) {
+    FILE *auto_file = fopen("autocorrelation.dat", "a");
+    for (uint64_t i = 0; i < n_runs + 1; i++) { 
+      fprintf(auto_file, "%" PRIu64 " %.15f ", auto_time[i], auto_energy[i]);
+    }
+    fprintf(auto_file, "\n");
+    fclose(auto_file);
+    free(auto_time);
+    free(auto_energy);
+  }
+
   gsl_rng_free(r);
   graph_free(s->g);
   free(s->spins);