implemented updating the first fourier moment in wolff_finite, but also make wolff_finite obselete by adding a hacky preprocessor method for making wolff.h as efficient

8 files changed, 119 insertions, 2 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index 5124061..b8c98e5 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -64,7 +64,11 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
             non_ghost = vn;
           }
           double dE = state->H(rs_old) - state->H(rs_new);
+#ifdef FINITE_STATES
+          prob = H_probs[N_STATES * state_to_ind(rs_old) + state_to_ind(rs_new)];
+#else
           prob = 1.0 - exp(-dE / state->T);
+#endif
 
           add(&(state->M), -1, rs_old);
           add(&(state->M),  1, rs_new);
@@ -84,7 +88,11 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
           free_spin (rs_new);
         } else {
           double dE = state->J(si_old, sj) - state->J(si_new, sj);
+#ifdef FINITE_STATES
+          prob = J_probs[N_STATES * N_STATES * state_to_ind(si_old) + N_STATES * state_to_ind(si_new) + state_to_ind(sj)];
+#else
           prob = 1.0 - exp(-dE / state->T);
+#endif
           state->E += dE;
         }
 
diff --git a/lib/cluster_finite.c b/lib/cluster_finite.c
index 9392cf8..2041df0 100644
--- a/lib/cluster_finite.c
+++ b/lib/cluster_finite.c
@@ -34,7 +34,7 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) {
       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;
+        q_t sn, non_ghost;
         double prob;
         bool external_neighbor = false;
 
@@ -42,8 +42,10 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) {
 
         if (vn == s->g->nv - 1) {
           external_neighbor = true;
+          non_ghost = v;
         } else {
           sn = s->spins[vn];
+          non_ghost = vn;
         }
 
         if (external_flipped || external_neighbor) {
@@ -62,6 +64,16 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) {
           s->M[rot_s_old]--;
           s->M[rot_s_new]++;
 
+          for (D_t i = 0; i < s->D; i++) {
+            L_t x = (non_ghost / (v_t)pow(s->L, s->D - i - 1)) % s->L;
+
+            s->ReF[s->D * i + rot_s_old] -= s->precomputed_cos[i];
+            s->ReF[s->D * i + rot_s_new] += s->precomputed_cos[i];
+
+            s->ImF[s->D * i + rot_s_old] -= s->precomputed_sin[i];
+            s->ImF[s->D * i + rot_s_new] += s->precomputed_sin[i];
+          }
+
         } else {
           q_t diff_old = s->bond_with_zero_type[s->transformations[s->q * s->transform_site_to_zero[sn] + s_old]];
           q_t diff_new = s->bond_with_zero_type[s->transformations[s->q * s->transform_site_to_zero[sn] + s_new]];
diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h
index b2d764e..117087e 100644
--- a/lib/cluster_finite.h
+++ b/lib/cluster_finite.h
@@ -45,6 +45,10 @@ typedef struct {
   q_t *R;
   v_t *B;
   v_t *M;
+  double *ReF;
+  double *ImF;
+  double *precomputed_cos;
+  double *precomputed_sin;
 } state_finite_t;
 
 v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot, gsl_rng *r);
diff --git a/lib/finite_states.h b/lib/finite_states.h
new file mode 100644
index 0000000..1e8800c
--- /dev/null
+++ b/lib/finite_states.h
@@ -0,0 +1,33 @@
+
+#pragma once
+
+#include <cmath>
+#include <functional>
+
+#define FINITE_STATES
+
+// must have N_STATES, states[N_STATES], and state_to_ind defined before
+// invoking header
+
+double J_probs[N_STATES * N_STATES * N_STATES];
+double H_probs[N_STATES * N_STATES];
+
+template <class X_t>
+void initialize_probs(std::function <double(X_t, X_t)> J, std::function <double(X_t)> H, double T) {
+  for (q_t i = 0; i < N_STATES; i++) {
+    for (q_t j = 0; j < N_STATES; j++) {
+      for (q_t k = 0; k < N_STATES; k++) {
+        J_probs[i * N_STATES * N_STATES + j * N_STATES +k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T);
+      }
+    }
+  }
+  for (q_t i = 0; i < N_STATES; i++) {
+    for (q_t j = 0; j < N_STATES; j++) {
+      H_probs[i * N_STATES + j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T);
+    }
+  }
+}
+
+
+
+
diff --git a/lib/initial_finite.c b/lib/initial_finite.c
index 6ea76ef..59e7ec4 100644
--- a/lib/initial_finite.c
+++ b/lib/initial_finite.c
@@ -1,6 +1,24 @@
 
 #include "initial_finite.h"
 
+double *precompute_cos(L_t L) {
+  double *x = (double *)malloc(L * sizeof(double));
+  for (L_t i = 0; i < L; i++) {
+    x[i] = cos(2 * M_PI * (double)i / (double)L);
+  }
+
+  return x;
+}
+
+double *precompute_sin(L_t L) {
+  double *x = (double *)malloc(L * sizeof(double));
+  for (L_t i = 0; i < L; i++) {
+    x[i] = sin(2 * M_PI * (double)i / (double)L);
+  }
+
+  return x;
+}
+
 double *Jprobs_from_J(q_t q, double T, double *J) {
   double *J_probs = (double *)calloc(pow(q, 2), sizeof(double));
 
@@ -112,6 +130,13 @@ state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H)
   s->B = (v_t *)calloc(2, sizeof(v_t));
   s->B[0] = s->ne;
 
+  s->ReF = (double *)calloc(D * 2, sizeof(double));
+  s->ImF = (double *)calloc(D * 2, sizeof(double));
+
+
+  s->precomputed_cos = precompute_cos(L);
+  s->precomputed_sin = precompute_sin(L);
+
   return s;
 }
 
@@ -168,6 +193,12 @@ state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, doub
   s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t));
   s->B[0] = s->ne; // everyone starts in state 0, remember?
 
+  s->ReF = (double *)calloc(D * q, sizeof(double));
+  s->ImF = (double *)calloc(D * q, sizeof(double));
+
+  s->precomputed_cos = precompute_cos(L);
+  s->precomputed_sin = precompute_sin(L);
+
   return s;
 }
 
@@ -230,10 +261,15 @@ state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, doub
   s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t));
   s->B[0] = s->ne; // everyone starts in state 0, remember?
 
+  s->ReF = (double *)calloc(D * q, sizeof(double));
+  s->ImF = (double *)calloc(D * q, sizeof(double));
+
+  s->precomputed_cos = precompute_cos(L);
+  s->precomputed_sin = precompute_sin(L);
+
   return s;
 }
 
-
 state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double *H) {
   state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t));
 
@@ -291,6 +327,12 @@ state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double
   s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t));
   s->B[0] = s->nv; // everyone starts in state 0, remember?
 
+  s->ReF = (double *)calloc(D * q, sizeof(double));
+  s->ImF = (double *)calloc(D * q, sizeof(double));
+
+  s->precomputed_cos = precompute_cos(L);
+  s->precomputed_sin = precompute_sin(L);
+
   return s;
 }
 
@@ -321,6 +363,10 @@ void state_finite_free(state_finite_t *s) {
   free(s->involutions);
   free(s->transform_site_to_zero);
   free(s->bond_with_zero_type);
+  free(s->ReF);
+  free(s->ImF);
+  free(s->precomputed_cos);
+  free(s->precomputed_sin);
   free(s);
 }
 
diff --git a/lib/wolff.h b/lib/wolff.h
index 8286024..da4b7b6 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -5,6 +5,10 @@
 template <class R_t, class X_t>
 void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s0)> gen_R, std::function <void(const state_t <R_t, X_t> *)> measurements, gsl_rng *r, bool silent) {
 
+#ifdef FINITE_STATES
+  initialize_probs(s->J, s->H, s->T);
+#endif
+
   if (!silent) printf("\n");
   for (count_t steps = 0; steps < N; steps++) {
     if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s->E, s->last_cluster_size);
diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp
index 9866248..f614707 100644
--- a/src/wolff_ising.cpp
+++ b/src/wolff_ising.cpp
@@ -8,6 +8,11 @@
 #include <z2.h>
 #include <ising.h>
 
+#define N_STATES 2
+const ising_t states[N_STATES] = {false, true};
+q_t state_to_ind(ising_t state) { return (q_t)state.x; }
+#include <finite_states.h>
+
 // include wolff.h
 #include <rand.h>
 #include <wolff.h>
diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp
index e3259e4..6b6f602 100644
--- a/src/wolff_potts.cpp
+++ b/src/wolff_potts.cpp
@@ -10,6 +10,11 @@
 #include <potts.h>
 #include <colors.h>
 
+#define N_STATES POTTSQ
+const potts_t<POTTSQ> states[8] = {0, 1, 2, 3, 4, 5, 6, 7};
+q_t state_to_ind(potts_t<POTTSQ> state) { return (q_t)state.x; }
+#include <finite_states.h>
+
 // include wolff.h
 #include <rand.h>
 #include <wolff.h>