various changes

7 files changed, 90 insertions, 49 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index f3271b0..9dcf50d 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -25,78 +25,64 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, gsl_rng *rand) {
     stack.pop();
 
     if (!marks[v]) { // don't reprocess anyone we've already visited!
-      X_t si_old, si_new;
-      R_t R_old, R_new;
+      X_t si_new;
+      R_t R_new;
 
-      R_old = s.R;
       marks[v] = true;
 
       bool v_is_ghost = (v == s.nv); // ghost site has the last index
 
       if (v_is_ghost) {
-        R_new = r.act(R_old); // if we are, then we're moving the transformation
+        R_new = r.act(s.R); // if we are, then we're moving the transformation
       } else {
-        si_old = s.spins[v];
-        si_new = r.act(si_old); // otherwise, we're moving the spin at our site
+        si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site
       }
 
       for (const v_t &vn : s.g.v_adj[v]) {
-        X_t sj;
         bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
 
-        if (!vn_is_ghost) {
-          sj = s.spins[vn];
-        }
-
-        double prob;
+        double dE, prob;
 
-        if (v_is_ghost || vn_is_ghost) { // if this is a ghost-involved bond...
+        if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond
           X_t rs_old, rs_new;
           v_t non_ghost;
+
           if (vn_is_ghost) {
-            rs_old = R_old.act_inverse(si_old);
-            rs_new = R_old.act_inverse(si_new);
+            // if our neighbor is the ghost, the current site is a normal
+            // spin - rotate it back!
+            rs_old = s.R.act_inverse(s.spins[v]);
+            rs_new = s.R.act_inverse(si_new);
             non_ghost = v;
           } else {
-            rs_old = R_old.act_inverse(sj);
-            rs_new = R_new.act_inverse(sj);
+            /* if we're the ghost, we need to rotate our neighbor back in
+               both the old and new ways */
+            rs_old = s.R.act_inverse(s.spins[vn]);
+            rs_new = R_new.act_inverse(s.spins[vn]);
             non_ghost = vn;
           }
 
-          double dE = s.H(rs_old) - s.H(rs_new);
+          dE = s.H(rs_old) - s.H(rs_new);
 
 #ifdef FINITE_STATES
           prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
-#else
-          prob = 1.0 - exp(-dE / s.T);
 #endif
 
-          s.M -= rs_old;
-          s.M += rs_new;
-
-          s.E += dE;
-
-          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[i] -= rs_old * s.precomputed_cos[x];
-            s.ReF[i] += rs_new * s.precomputed_cos[x];
-
-            s.ImF[i] -= rs_old * s.precomputed_sin[x];
-            s.ImF[i] += rs_new * s.precomputed_sin[x];
-          }
-        } else { // otherwise, we're at a perfectly normal bond!
-          double dE = s.J(si_old, sj) - s.J(si_new, sj);
+          s.update_magnetization(rs_old, rs_new);
+          s.update_fourierZero(non_ghost, rs_old, rs_new);
+        } else { // this is a perfectly normal bond!
+          dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]);
 
 #ifdef FINITE_STATES
-          prob = J_probs[state_to_ind(si_old)][state_to_ind(si_new)][state_to_ind(sj)];
-#else
-          prob = 1.0 - exp(-dE / s.T);
+          prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])];
 #endif
-
-          s.E += dE;
         }
 
+        s.update_energy(dE);
+
+#ifndef FINITE_STATES
+        prob = 1.0 - exp(-dE / s.T);
+#endif
+
         if (gsl_rng_uniform(rand) < prob) {
           stack.push(vn); // push the neighboring vertex to the stack
         }
diff --git a/lib/correlation.h b/lib/correlation.h
index d08f713..29357a5 100644
--- a/lib/correlation.h
+++ b/lib/correlation.h
@@ -10,7 +10,11 @@ template <class R_t, class X_t>
 double correlation_length(const state_t <R_t, X_t>& s) {
   double total = 0;
 
+#ifdef DIMENSION
+  for (D_t j = 0; j < DIMENSION; j++) {
+#else
   for (D_t j = 0; j < s.D; j++) {
+#endif
     total += norm_squared(s.ReF[j]) + norm_squared(s.ImF[j]);
   }
 
diff --git a/lib/graph.cpp b/lib/graph.cpp
index 8c97274..ca251f3 100644
--- a/lib/graph.cpp
+++ b/lib/graph.cpp
@@ -6,13 +6,17 @@ graph_t::graph_t(D_t D, L_t L) {
   ne = D * nv;
 
   v_adj.resize(nv);
+  coordinate.resize(nv);
 
   for (std::vector<v_t> v_adj_i : v_adj) {
     v_adj_i.reserve(2 * D);
   }
 
   for (v_t i = 0; i < nv; i++) {
+    coordinate[i].resize(D);
     for (D_t j = 0; j < D; j++) {
+      coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
+
       v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
       v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
     }
@@ -25,12 +29,15 @@ void graph_t::add_ext() {
   }
 
   v_adj.resize(nv + 1);
+  coordinate.resize(nv + 1);
   v_adj[nv].reserve(nv);
 
   for (v_t i = 0; i < nv; i++) {
     v_adj[nv].push_back(i);
   }
 
+  coordinate[nv].resize(coordinate[0].size());
+
   ne += nv;
   nv += 1;
 }
diff --git a/lib/graph.h b/lib/graph.h
index a4732fb..de06924 100644
--- a/lib/graph.h
+++ b/lib/graph.h
@@ -13,6 +13,7 @@ class graph_t {
     v_t ne;
     v_t nv;
     std::vector<std::vector<v_t>> v_adj;
+    std::vector<std::vector<double>> coordinate;
 
     graph_t(D_t D, L_t L);
     void add_ext();
diff --git a/lib/ising.h b/lib/ising.h
index 473a18a..aaa0c52 100644
--- a/lib/ising.h
+++ b/lib/ising.h
@@ -53,6 +53,18 @@ class ising_t {
         return a;
       }
     }
+
+    inline int operator-(const ising_t &s) const {
+      if (x == s.x) {
+        return 0;
+      } else {
+        if (x) {
+          return -2;
+        } else {
+          return 2;
+        }
+      }
+    }
 };
 
 inline int& operator+=(int& M, const ising_t &s) {
diff --git a/lib/state.h b/lib/state.h
index 5abf65b..550d100 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -9,6 +9,10 @@
 
 template <class R_t, class X_t>
 class state_t {
+  private:
+    // updating fourier terms F requires many cos and sin calls, faster to do it beforehand.
+    std::vector<std::vector<double>> precomputed_cos;
+    std::vector<std::vector<double>> precomputed_sin;
   public:
     D_t D;
     L_t L;
@@ -23,9 +27,6 @@ class state_t {
     v_t last_cluster_size;
     std::vector<typename X_t::F_t> ReF;
     std::vector<typename X_t::F_t> ImF;
-    // updating fourier terms F requires many cos and sin calls, faster to do it beforehand.
-    std::vector<double> precomputed_cos;
-    std::vector<double> precomputed_sin;
 
     std::function <double(const X_t&, const X_t&)> J;
     std::function <double(const X_t&)> H;
@@ -44,11 +45,35 @@ class state_t {
         ReF[i] = spins[0] * 0.0;
         ImF[i] = spins[0] * 0.0;
       }
-      precomputed_cos.resize(L);
-      precomputed_sin.resize(L);
-      for (L_t i = 0; i < L; i++) {
-        precomputed_cos[i] = cos(2 * M_PI * (double)i / (double)L);
-        precomputed_sin[i] = sin(2 * M_PI * (double)i / (double)L);
+      precomputed_cos.resize(nv);
+      precomputed_sin.resize(nv);
+      for (v_t i = 0; i < nv; i++) {
+        precomputed_cos[i].resize(D);
+        precomputed_sin[i].resize(D);
+        for (D_t j = 0; j < D; j++) {
+          precomputed_cos[i][j] = cos(2 * M_PI * g.coordinate[i][j] / (double)L);
+          precomputed_sin[i][j] = sin(2 * M_PI * g.coordinate[i][j] / (double)L);
+        }
+      }
+    }
+
+    void update_magnetization(const X_t& s_old, const X_t& s_new) {
+      M -= s_old;
+      M += s_new;
+    }
+
+    void update_energy(const double& dE) {
+      E += dE;
+    }
+
+    void update_fourierZero(v_t v, const X_t& s_old, const X_t& s_new) {
+#ifdef DIMENSION
+      for (D_t i = 0; i < DIMENSION; i++) {
+#else
+      for (D_t i = 0; i < D; i++) {
+#endif
+        ReF[i] += (s_new - s_old) * precomputed_cos[v][i];
+        ImF[i] += (s_new - s_old) * precomputed_sin[v][i];
       }
     }
 };
diff --git a/lib/vector.h b/lib/vector.h
index 72633a8..7d0ee36 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -61,6 +61,12 @@ class vector_t : public std::array<T, q> {
 
       return result;
     }
+
+    inline vector_t<q, T> operator-(const vector_t<q, T>& v) const {
+      vector_t<q, T> diff = *this;
+      diff -= v;
+      return diff;
+    }
 };