1 files changed, 320 insertions, 0 deletions
diff --git a/lib/wolff.hpp b/lib/wolff.hpp
new file mode 100644
index 0000000..d21d788
--- /dev/null
+++ b/lib/wolff.hpp
@@ -0,0 +1,320 @@
+
+#ifndef WOLFF_H
+#define WOLFF_H
+
+#include <functional>
+#include <vector>
+#include <random>
+#include <cmath>
+#include <iterator>
+#include <list>
+#include <tuple>
+#include <queue>
+
+namespace wolff{
+
+  template <class vertex_prop = std::tuple<>, class edge_prop = std::tuple<>>
+  class graph {
+    /* the graph class describes a lattice on which wolff is run. for most
+     * purposes, using the default square lattice constructor is sufficient,
+     * but arbitrary lattices can be constructed
+     */
+    public:
+      unsigned D;  // dimension of space
+      unsigned L;  // linear size
+      unsigned ne; // number of edges
+      unsigned nv; // number of vertices
+
+      struct _vertex;
+
+      typedef struct _halfedge {
+        struct _vertex &self;     // reference to the vertex this edge comes from
+        struct _vertex &neighbor; // reference to the vertex this edge goes to
+        edge_prop prop;           // optional properties
+
+        _halfedge(struct _vertex &v1, struct _vertex &v2) : self(v1), neighbor(v2) {};
+      } halfedge;
+
+      typedef struct _vertex {
+        unsigned ind;              // index of the vertex
+        std::list<halfedge> edges; // list of edges incident on this vertex
+        vertex_prop prop;          // optional properties
+      } vertex;
+
+      std::vector<vertex> vertices; /* vector of vertices, length nv, with 
+                                     * vertices[i].ind = i for all 0 <= i < nv
+                                     */
+      graph() {
+        // default constructor for empty graph. use it to build your own!
+        D = 0;
+        L = 0;
+        nv = 0;
+        ne = 0;
+      };
+
+      graph(unsigned D, unsigned L) : D(D), L(L) {
+        // default constructor for square lattice graph
+        nv = pow(L, D);
+        ne = D * nv;
+
+        vertices.resize(nv);
+
+        for (unsigned i = 0; i < nv; i++) {
+          vertices[i].ind = i;
+          for (unsigned j = 0; j < D; j++) {
+            unsigned n1 = pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) +
+                                              fmod(i + pow(L, j), pow(L, j + 1));
+            unsigned n2 = pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) +
+                              fmod(pow(L, j + 1) + i - pow(L, j), pow(L, j + 1));
+
+            halfedge f(vertices[i], vertices[n1]);
+            halfedge b(vertices[i], vertices[n2]);
+
+            vertices[i].edges.push_back(f);
+            vertices[i].edges.push_back(b);
+          }
+        }
+      };
+
+      void add_ghost() {
+        // adds a ghost site to any graph
+        vertices.resize(nv + 1);
+
+        for (auto it = vertices.begin(); it != std::prev(vertices.end()); it++) {
+          halfedge e1(*it, vertices[nv]);
+          it->edges.push_back(e1);
+
+          halfedge e2(vertices[nv], *it);
+          vertices[nv].edges.push_back(e2);
+        }
+
+        vertices[nv].ind = nv;
+
+        ne += nv;
+        nv++;
+      };
+  };
+
+  template <class R_t, class X_t, class G_t = graph<>>
+  class measurement;
+
+  template <class R_t, class X_t, class G_t = graph<>>
+  class system {
+    public:
+      unsigned nv; // number of vertices
+      unsigned ne; // number of edges
+      G_t G; // the graph defining the lattice with ghost
+      double T; // the temperature
+      std::vector<X_t> s; // the state of the ordinary spins
+
+#ifdef WOLFF_BOND_DEPENDENCE
+      std::function <double(const typename G_t::halfedge&, const X_t&, const X_t&)> Z; // coupling between sites
+#else
+      std::function <double(const X_t&, const X_t&)> Z; // coupling between sites
+#endif
+
+#ifndef WOLFF_NO_FIELD
+      R_t s0; // the current state of the ghost site
+#ifdef WOLFF_SITE_DEPENDENCE
+      std::function <double(const typename G_t::vertex&, const X_t&)> B; // coupling with the external field
+#else
+      std::function <double(const X_t&)> B; // coupling with the external field
+#endif
+#endif
+
+#ifdef WOLFF_USE_FINITE_STATES
+      std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp;
+#ifndef WOLFF_NO_FIELD
+      std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Bp;
+#endif
+#endif
+
+      system(G_t g, double T,
+#ifdef WOLFF_BOND_DEPENDENCE
+          std::function <double(const typename G_t::halfedge&, const X_t&, const X_t&)> Z
+#else
+          std::function <double(const X_t&, const X_t&)> Z
+#endif
+#ifndef WOLFF_NO_FIELD
+#ifdef WOLFF_SITE_DEPENDENCE
+          , std::function <double(const typename G_t::vertex&, const X_t&)> B
+#else
+          , std::function <double(const X_t&)> B
+#endif
+#endif
+          ) : G(g), T(T), Z(Z)
+#ifndef WOLFF_NO_FIELD
+               , s0(), B(B)
+#endif
+      {
+        nv = G.nv;
+        ne = G.ne;
+        s.resize(nv);
+#ifndef WOLFF_NO_FIELD
+        G.add_ghost();
+#endif
+#ifdef WOLFF_USE_FINITE_STATES
+        this->finite_states_init();
+#endif
+      }
+
+      void flip_cluster(typename G_t::vertex& v0, const R_t& r, std::mt19937& rng,
+                          measurement<R_t, X_t, G_t>& A) {
+        std::uniform_real_distribution<double> dist(0.0, 1.0);
+
+        std::queue<unsigned> queue;
+        queue.push(v0.ind);
+
+        std::vector<bool> visited(G.nv, false);
+
+        while (!queue.empty()) {
+          unsigned i = queue.front();
+          queue.pop();
+
+          if (!visited[i]) { // don't reprocess anyone we've already visited!
+            visited[i] = true;
+
+            X_t si_new;
+#ifndef WOLFF_NO_FIELD
+            R_t s0_new;
+
+            bool we_are_ghost = (i == nv);
+
+            if (we_are_ghost) {
+              s0_new = r.act(s0);
+            } else
+#endif
+            {
+              si_new = r.act(s[i]);
+            }
+
+            for (const typename G_t::halfedge &e : G.vertices[i].edges) {
+              double dE, p;
+              unsigned j = e.neighbor.ind;
+
+#ifndef WOLFF_NO_FIELD
+              bool neighbor_is_ghost = (j == nv);
+
+              if (we_are_ghost || neighbor_is_ghost) {
+                X_t s0s_old, s0s_new;
+                unsigned non_ghost;
+
+                if (neighbor_is_ghost) {
+                  non_ghost = i;
+                  s0s_old = s0.act_inverse(s[i]);
+                  s0s_new = s0.act_inverse(si_new);
+                } else {
+                  non_ghost = j;
+                  s0s_old = s0.act_inverse(s[j]);
+                  s0s_new = s0_new.act_inverse(s[j]);
+                }
+
+#ifdef WOLFF_SITE_DEPENDENCE
+                dE = B(G.vertices[non_ghost], s0s_old) - B(G.vertices[non_ghost], s0s_new);
+#else
+                dE = B(s0s_old) - B(s0s_new);
+#endif
+
+#ifdef WOLFF_USE_FINITE_STATES
+                p = Bp[s0s_old.enumerate()][s0s_new.enumerate()];
+#endif
+
+                // run measurement hooks for encountering a ghost bond
+                A.ghost_bond_visited(*this, G.vertices[non_ghost], s0s_old, s0s_new, dE);
+              } else // this is a perfectly normal bond!
+#endif
+              {
+#ifdef WOLFF_BOND_DEPENDENCE
+                dE = Z(e, s[i], s[j]) - Z(e, si_new, s[j]);
+#else
+                dE = Z(s[i], s[j]) - Z(si_new, s[j]);
+#endif
+
+#ifdef WOLFF_USE_FINITE_STATES
+                p = Zp[s[i].enumerate()][si_new.enumerate()][s[j].enumerate()];
+#endif
+
+                // run measurement hooks for encountering a plain bond
+                A.plain_bond_visited(*this, e, si_new, dE);
+              }
+
+#ifndef WOLFF_USE_FINITE_STATES
+              p = 1.0 - exp(-dE / T);
+#endif
+
+              if (dist(rng) < p) {
+                queue.push(j); // push the neighboring vertex to the queue 
+              }
+            }
+
+#ifndef WOLFF_NO_FIELD
+            if (we_are_ghost) {
+              A.ghost_site_transformed(*this, s0_new);
+              s0 = s0_new;
+            } else
+#endif
+            {
+              A.plain_site_transformed(*this, G.vertices[i], si_new);
+              s[i] = si_new;
+            }
+          }
+        }
+      }
+
+      void run_wolff(unsigned N,
+          std::function <R_t(std::mt19937&, const system<R_t, X_t, G_t>&, const typename G_t::vertex&)> r_gen, measurement<R_t, X_t, G_t>& A, std::mt19937& rng) {
+        std::uniform_int_distribution<unsigned> dist(0, nv - 1);
+
+        for (unsigned n = 0; n < N; n++) {
+          unsigned i0 = dist(rng);
+          R_t r = r_gen(rng, *this, G.vertices[i0]);
+
+          A.pre_cluster(n, N, *this, G.vertices[i0], r);
+
+          this->flip_cluster(G.vertices[i0], r, rng, A);
+
+          A.post_cluster(n, N, *this);
+        }
+      }
+
+#ifdef WOLFF_USE_FINITE_STATES
+      void finite_states_init() {
+#ifndef WOLFF_NO_FIELD
+        for (unsigned i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+          for (unsigned j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+            Bp[i][j] = 1.0 - exp(-(B(X_t(i)) - B(X_t(j))) / T);
+          }
+        }
+#endif
+        for (unsigned i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+          for (unsigned j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+            for (unsigned k = 0; k < WOLFF_FINITE_STATES_N; k++) {
+              Zp[i][j][k] = 1.0 - exp(-(Z(X_t(i), X_t(k)) - Z(X_t(j), X_t(k))) / T);
+            }
+          }
+        }
+      }
+#endif
+
+  };
+
+  template <class R_t, class X_t, class G_t>
+  class measurement {
+    public:
+      virtual void pre_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const R_t&) {};
+
+      virtual void plain_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::halfedge& e, const X_t&, double) {};
+      virtual void plain_site_transformed(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const X_t&) {};
+
+#ifndef WOLFF_NO_FIELD
+      virtual void ghost_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const X_t&, const X_t&, double) {};
+      virtual void ghost_site_transformed(const system<R_t, X_t, G_t>&, const R_t&) {};
+#endif
+
+      virtual void post_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&) {};
+  };
+
+}
+
+#endif
+