1 files changed, 111 insertions, 0 deletions
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
new file mode 100644
index 0000000..104f3c2
--- /dev/null
+++ b/lib/include/wolff/cluster.hpp
@@ -0,0 +1,111 @@
+
+#pragma once
+
+#include <random>
+#include <cmath>
+#include <vector>
+#include <stack>
+
+#include "types.h"
+#include "state.hpp"
+#include "graph.hpp"
+
+template <class R_t, class X_t>
+void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand) {
+  std::uniform_real_distribution<double> dist(0.0,1.0);
+  v_t nv = 0;
+
+  std::stack<v_t> stack;
+  stack.push(v0);
+
+  std::vector<bool> marks(s.g.nv, false);
+
+  while (!stack.empty()) {
+    v_t v = stack.top();
+    stack.pop();
+
+    if (!marks[v]) { // don't reprocess anyone we've already visited!
+      marks[v] = true;
+
+      X_t si_new;
+#ifndef NOFIELD
+      R_t R_new;
+
+      bool v_is_ghost = (v == s.nv); // ghost site has the last index
+
+      if (v_is_ghost) {
+        R_new = r.act(s.R); // if we are, then we're moving the transformation
+      } else
+#endif
+      {
+        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]) {
+        double dE, prob;
+
+#ifndef NOFIELD
+        bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
+
+        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) {
+            // 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 {
+            /* 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;
+          }
+
+          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)];
+#endif
+
+          s.update_magnetization(rs_old, rs_new);
+          s.update_fourierZero(non_ghost, rs_old, rs_new);
+        } else // this is a perfectly normal bond!
+#endif
+        {
+          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(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])];
+#endif
+        }
+
+        s.update_energy(dE);
+
+#ifndef FINITE_STATES
+        prob = 1.0 - exp(-dE / s.T);
+#endif
+
+        if (dist(rand) < prob) {
+          stack.push(vn); // push the neighboring vertex to the stack
+        }
+      }
+
+#ifndef NOFIELD
+      if (v_is_ghost) {
+        s.R = R_new;
+      } else
+#endif
+      {
+        s.spins[v] = si_new;
+        nv++;
+      }
+    }
+  }
+
+  s.last_cluster_size = nv;
+}
+