began generalizing the potts code to work with all finite spin spaces

10 files changed, 321 insertions, 20 deletions

diff --git a/lib/cluster_finite.c b/lib/cluster_finite.c
new file mode 100644
index 0000000..f11a3ea
--- /dev/null
+++ b/lib/cluster_finite.c
@@ -0,0 +1,100 @@
+
+#include "cluster_finite.h"
+
+v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) {
+  q_t *rot = s->transformations + s->q * rot_ind;
+  q_t *R_inv = symmetric_invert(s->q, s->R);
+  v_t nv = 0;
+
+  ll_t *stack = NULL;     // create a new stack
+  stack_push(&stack, v0); // push the initial vertex to the stack
+
+  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
+
+  while (stack != NULL) {
+    v_t v = stack_pop(&stack);
+
+    if (!marks[v]) {
+      q_t s_old, s_new;
+      q_t *R_new, *R_inv_new; 
+      bool external_flipped;
+
+      marks[v] = true;
+
+      if (v == s->g->nv - 1) {
+        R_new = symmetric_compose(s->q, rot, s->R);
+        R_inv_new = symmetric_invert(s->q, R_new);
+        external_flipped = true;
+      } else {
+        s_old = s->spins[v];
+        s_new = symmetric_act(rot, s_old);
+        external_flipped = false;
+      }
+
+      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
+
+      for (v_t i = 0; i < nn; i++) {
+        q_t sn;
+        double prob;
+        bool external_neighbor = false;
+
+        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
+
+        if (vn == s->g->nv - 1) {
+          external_neighbor = true;
+        } else {
+          sn = s->spins[vn];
+        }
+
+        if (external_flipped || external_neighbor) {
+          q_t rot_s_old, rot_s_new;
+
+          if (external_neighbor) {
+            rot_s_old = symmetric_act(R_inv, s_old);
+            rot_s_new = symmetric_act(R_inv, s_new);
+          } else {
+            rot_s_old = symmetric_act(R_inv, sn);
+            rot_s_new = symmetric_act(R_inv_new, sn);
+          }
+
+          prob = s->H_probs[rot_s_new * s->q + rot_s_old];
+
+          s->M[rot_s_old]--;
+          s->M[rot_s_new]++;
+
+          s->E += - s->H[rot_s_new] + s->H[rot_s_old];
+        } else {
+          q_t diff_old = (s_old + s->q - sn) % s->q;
+          q_t diff_new = (s_new + s->q - sn) % s->q;
+
+          prob = s->J_probs[diff_new * s->q + diff_old];
+
+          s->E += - s->J[diff_new] + s->J[diff_old];
+        }
+
+        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
+          stack_push(&stack, vn); // push the neighboring vertex to the stack
+        }
+      }
+
+      if (external_flipped) {
+        free(s->R);
+        free(R_inv);
+        s->R = R_new;
+        R_inv = R_inv_new;
+      } else {
+        s->spins[v] = s_new;
+      }
+
+      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
+        nv++;
+      }
+    }
+  }
+
+  free(marks);
+  free(R_inv);
+
+  return nv;
+}
+
diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h
new file mode 100644
index 0000000..abdc8fc
--- /dev/null
+++ b/lib/cluster_finite.h
@@ -0,0 +1,42 @@
+
+#pragma once
+
+#include <assert.h>
+#include <fftw3.h>
+#include <float.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_rng.h>
+#include <inttypes.h>
+#include <math.h>
+#include <stdbool.h>
+#include <string.h>
+#include <sys/types.h>
+
+#include "types.h"
+#include "rand.h"
+#include "stack.h"
+#include "convex.h"
+#include "graph.h"
+#include "tree.h"
+#include "measurement.h"
+#include "symmetric.h"
+#include "yule_walker.h"
+
+typedef struct {
+  graph_t *g;
+  q_t q;
+  R_t n_transformations;
+  q_t *transformations;
+  double T;
+  double *J;
+  double *H;
+  double *J_probs;
+  double *H_probs;
+  q_t *spins;
+  q_t *R;
+  double E;
+  v_t *M;
+} state_finite_t;
+
+v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot, gsl_rng *r);
+
diff --git a/lib/dihedral.c b/lib/dihedral.c
index 623041a..ac74a23 100644
--- a/lib/dihedral.c
+++ b/lib/dihedral.c
@@ -25,4 +25,16 @@ q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) {
   }
 }
 
+q_t *dihedral_gen_transformations(q_t q) {
+  q_t *transformations = (q_t *)malloc(q * q * sizeof(q_t));
+
+  for (q_t i = 0; i < q; i++) {
+    for (q_t j = 0; j < q; j++) {
+      transformations[q * i + j] = dihedral_act(q, i, j);
+    }
+  }
+
+  return transformations;
+}
+
 
diff --git a/lib/dihedral.h b/lib/dihedral.h
index 813e796..e5e4cbd 100644
--- a/lib/dihedral.h
+++ b/lib/dihedral.h
@@ -15,3 +15,5 @@ q_t dihedral_act(q_t q, q_t gi, q_t s);
 
 q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s);
 
+q_t *dihedral_gen_transformations(q_t q);
+
diff --git a/lib/measurement.h b/lib/measurement.h
index eaa260b..46c034f 100644
--- a/lib/measurement.h
+++ b/lib/measurement.h
@@ -24,6 +24,11 @@ typedef struct {
   double O2;
 } autocorr_t;
 
+typedef struct {
+  void (*f)(state_finite_t *, void *);
+  void *data;
+} measurement_t;
+
 void meas_update(meas_t *m, double x);
 
 double meas_dx(const meas_t *m);
diff --git a/lib/symmetric.c b/lib/symmetric.c
new file mode 100644
index 0000000..729b38c
--- /dev/null
+++ b/lib/symmetric.c
@@ -0,0 +1,27 @@
+
+#include "symmetric.h"
+
+q_t *symmetric_compose(q_t q, const q_t *g1, const q_t *g2) {
+  q_t *g3 = (q_t *)malloc(q * sizeof(q_t));
+
+  for (q_t i = 0; i < q; i++) {
+    g3[i] = g1[g2[i]];
+  }
+
+  return g3;
+}
+
+q_t symmetric_act(const q_t *g, q_t s) {
+  return g[s];
+}
+
+q_t *symmetric_invert(q_t q, const q_t *g) {
+  q_t *g_inv = (q_t *)malloc(q * sizeof(q_t));
+
+  for (q_t i = 0; i < q; i++) {
+    g_inv[g[i]] = i;
+  }
+
+  return g_inv;
+}
+
diff --git a/lib/symmetric.h b/lib/symmetric.h
new file mode 100644
index 0000000..6e00f52
--- /dev/null
+++ b/lib/symmetric.h
@@ -0,0 +1,13 @@
+
+#pragma once
+
+#include <stdlib.h>
+
+#include "types.h"
+
+q_t *symmetric_compose(q_t q, const q_t *g1, const q_t *g2);
+
+q_t symmetric_act(const q_t *g, q_t s);
+
+q_t *symmetric_invert(q_t q, const q_t *g);
+
diff --git a/lib/types.h b/lib/types.h
index fcc2ce7..572ec64 100644
--- a/lib/types.h
+++ b/lib/types.h
@@ -3,6 +3,7 @@
 
 typedef uint_fast32_t v_t;
 typedef uint_fast8_t q_t;
+typedef uint_fast16_t R_t;
 typedef uint_fast8_t D_t;
 typedef uint_fast16_t L_t;
 typedef uint_fast64_t count_t;
@@ -10,6 +11,7 @@ typedef int_fast64_t h_t;
 
 #define MAX_v UINT_FAST32_MAX
 #define MAX_Q UINT_FAST8_MAX
+#define MAX_R UINT_FAST16_MAX
 #define MAX_D UINT_FAST8_MAX
 #define MAX_L UINT_FAST16_MAX
 #define MAX_COUNT UINT_FAST64_MAX
@@ -18,6 +20,7 @@ typedef int_fast64_t h_t;
 
 #define PRIv PRIuFAST32
 #define PRIq PRIuFAST8
+#define PRIR PRIuFAST16
 #define PRID PRIuFAST8
 #define PRIL PRIuFAST16
 #define PRIcount PRIuFAST64
diff --git a/lib/wolff_finite.c b/lib/wolff_finite.c
new file mode 100644
index 0000000..64de9ba
--- /dev/null
+++ b/lib/wolff_finite.c
@@ -0,0 +1,70 @@
+
+#include "cluster_finite.h"
+
+void wolff_finite(state_finite_t *s, count_t sweeps, count_t sweeps_per_measurement, count_t n_measurements, measurement_t *measurements) {
+  for (count_t i = 0; i < sweeps; i++) {
+
+    count_t n_flips = 0;
+
+    while (n_flips / h->nv < sweeps_per_measurement) {
+      v_t v0 = gsl_rng_uniform_int(r, h->nv);
+      R_t step;
+     
+      bool changed = false;
+      while (!changed) {
+        step = gsl_rng_uniform_int(r, s->n_transformations);
+        if v(symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) {
+          changed = true;
+        }
+      }
+
+      v_t tmp_flips = flip_cluster_finite(s, v0, step, r);
+      n_flips += tmp_flips;
+
+      if (n_runs > 0) {
+        n_steps++;
+        meas_update(clust, tmp_flips);
+
+        if (record_autocorrelation && n_steps % ac_skip == 0) {
+          update_autocorr(autocorr, s->E);
+        }
+
+      }
+
+    }
+
+    for (q_t i = 0; i < q; i++) {
+      meas_update(M[i], s->M[i]);
+    }
+    meas_update(E, s->E);
+
+    q_t n_at_max = 0;
+    q_t max_M_i = 0;
+    v_t max_M = 0;
+
+    for (q_t i = 0; i < q; i++) {
+      if (s->M[i] > max_M) {
+        n_at_max = 1;
+        max_M_i = i;
+        max_M = s->M[i];
+      } else if (s->M[i] == max_M) {
+        n_at_max++;
+      }
+    }
+
+    if (record_distribution) {
+      mag_dist[s->M[0]]++;
+    }
+
+    if (n_at_max == 1) {
+      for (q_t i = 0; i < q; i++) {
+        meas_update(sM[max_M_i][i], s->M[i]);
+      }
+      meas_update(sE[max_M_i], s->E);
+      freqs[max_M_i]++;
+    }
+
+    diff = fabs(meas_dx(clust) / clust->x);
+  }
+}
+
diff --git a/src/wolff_potts.c b/src/wolff_potts.c
index f40c216..b081bec 100644
--- a/src/wolff_potts.c
+++ b/src/wolff_potts.c
@@ -1,7 +1,8 @@
 
 #include <getopt.h>
 
-#include <cluster.h>
+#include <dihedral.h>
+#include <cluster_finite.h>
 
 int main(int argc, char *argv[]) {
 
@@ -18,6 +19,7 @@ int main(int argc, char *argv[]) {
   double eps = 0;
   bool pretend_ising = false;
   bool planar_potts = false;
+  bool sim_dgm = false;
   bool silent = false;
   bool snapshots = false;
   bool snapshot = false;
@@ -30,7 +32,7 @@ int main(int argc, char *argv[]) {
   q_t J_ind = 0;
   q_t H_ind = 0;
 
-  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:J:H:m:e:IpsSPak:W:d")) != -1) {
+  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:J:H:m:e:IpsSPak:W:dr")) != -1) {
     switch (opt) {
     case 'N':
       N = (count_t)atof(optarg);
@@ -91,6 +93,9 @@ int main(int argc, char *argv[]) {
     case 'd':
       record_distribution = true;
       break;
+    case 'r':
+      sim_dgm = true;
+      break;
     default:
       exit(EXIT_FAILURE);
     }
@@ -111,19 +116,27 @@ int main(int argc, char *argv[]) {
     }
   }
 
-  ising_state_t *s = (ising_state_t *)calloc(1, sizeof(ising_state_t));
+  if (sim_dgm) {
+    for (q_t i = 0; i < q / 2 + 1; i++) {
+      J[i] = -pow(i, 2);
+    }
+    for (q_t i = 1; i < (q + 1) / 2; i++) {
+      J[q - i] = -pow(i, 2);
+    }
+  }
+
+  state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t));
 
   graph_t *h = graph_create_square(D, L);
   s->g = graph_add_ext(h);
 
   s->q = q;
-
-  s->spins = (q_t *)calloc(h->nv, sizeof(q_t));
+  s->n_transformations = q;
+  s->transformations = dihedral_gen_transformations(q);
 
   s->T = T;
-  s->H = H;
   s->J = J;
-  s->R = (dihedral_t *)calloc(1, sizeof(dihedral_t));
+  s->H = H;
 
   s->J_probs = (double *)calloc(pow(q, 2), sizeof(double));
   for (q_t i = 0; i < q; i++) {
@@ -138,9 +151,18 @@ int main(int argc, char *argv[]) {
     }
   }
 
-  s->M = (v_t *)calloc(q, sizeof(v_t));
-  s->M[0] = h->nv;
+  s->spins = (q_t *)calloc(h->nv, sizeof(q_t)); // everyone starts in state 0
+  s->R = (q_t *)malloc(q * sizeof(q_t)); // transformation is the identity, (1 ... q)
+
+  for (q_t i = 0; i < q; i++) {
+    s->R[i] = i;
+  }
+
+  // energy is the number of edges times the energy of an aligned edge minus
+  // the number of vertices times the energy of a 0-aligned vertex 
   s->E = - ((double)h->ne) * s->J[0] - ((double)h->nv) * s->H[0];
+  s->M = (v_t *)calloc(q, sizeof(v_t));
+  s->M[0] = h->nv; // everyone starts in state 0, remember?
 
   double diff = 1e31;
   count_t n_runs = 0;
@@ -192,15 +214,17 @@ int main(int argc, char *argv[]) {
 
     while (n_flips / h->nv < n) {
       v_t v0 = gsl_rng_uniform_int(r, h->nv);
-      q_t step;
+      R_t step;
      
-      if (q == 2) {
-        step = 1;
-      } else {
-        step = gsl_rng_uniform_int(r, q);
+      bool changed = false;
+      while (!changed) {
+        step = gsl_rng_uniform_int(r, s->n_transformations);
+        if (symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) {
+          changed = true;
+        }
       }
 
-      v_t tmp_flips = flip_cluster(s, v0, step, r);
+      v_t tmp_flips = flip_cluster_finite(s, v0, step, r);
       n_flips += tmp_flips;
 
       if (n_runs > 0) {
@@ -211,9 +235,6 @@ int main(int argc, char *argv[]) {
           update_autocorr(autocorr, s->E);
         }
 
-        if (record_distribution) {
-          mag_dist[s->M[0]]++;
-        }
       }
 
     }
@@ -237,6 +258,10 @@ int main(int argc, char *argv[]) {
       }
     }
 
+    if (record_distribution) {
+      mag_dist[s->M[0]]++;
+    }
+
     if (n_at_max == 1) {
       for (q_t i = 0; i < q; i++) {
         meas_update(sM[max_M_i][i], s->M[i]);
@@ -262,12 +287,13 @@ int main(int argc, char *argv[]) {
   }
 
   if (snapshot) {
+    q_t *R_inv = symmetric_invert(q, s->R);
     FILE *snapfile = fopen("snapshot.m", "a");
     fprintf(snapfile, "{{");
     for (L_t i = 0; i < L; i++) {
       fprintf(snapfile, "{");
       for (L_t j = 0; j < L; j++) {
-        fprintf(snapfile, "%" PRIq, dihedral_inverse_act(q, s->R, s->spins[L * i + j]));
+        fprintf(snapfile, "%" PRIq, symmetric_act(R_inv, s->spins[L * i + j]));
         if (j != L - 1) {
           fprintf(snapfile, ",");
         }
@@ -277,7 +303,7 @@ int main(int argc, char *argv[]) {
         fprintf(snapfile, ",");
       }
     }
-    fprintf(snapfile, "},{%" PRIq ",%d}}\n", s->R->i, s->R->r);
+    fprintf(snapfile, "}}\n");
     fclose(snapfile);
   }
 
@@ -446,6 +472,7 @@ int main(int argc, char *argv[]) {
   free(s->M);
   free(s->spins);
   free(s->R);
+  free(s->transformations);
   graph_free(s->g);
   free(s);
   free(H);