big rearrangement of files to make libraries and example (research) files clearer, and changed to c++ std lib random numbers

6 files changed, 332 insertions, 0 deletions

diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp
new file mode 100644
index 0000000..c10a211
--- /dev/null
+++ b/lib/include/wolff.hpp
@@ -0,0 +1,35 @@
+
+#include "wolff/cluster.hpp"
+#include "wolff/state.hpp"
+
+template <class R_t, class X_t>
+void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(std::mt19937&, X_t)> gen_R, std::function <void(const state_t <R_t, X_t>&)> measurements, std::mt19937& r, bool silent) {
+
+#ifdef FINITE_STATES
+#ifdef NOFIELD
+  initialize_probs(s.J, s.T);
+#else
+  initialize_probs(s.J, s.H, s.T);
+#endif
+#endif
+
+  std::uniform_int_distribution<v_t> dist(0, s.nv);
+
+  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);
+
+    v_t v0 = dist(r);
+    R_t step = gen_R(r, s.spins[v0]);
+    flip_cluster <R_t, X_t> (s, v0, step, r);
+
+    measurements(s);
+  }
+
+  if (!silent) {
+    printf("\033[F\033[J");
+  }
+  printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size);
+
+}
+
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;
+}
+
diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp
new file mode 100644
index 0000000..426edad
--- /dev/null
+++ b/lib/include/wolff/finite_states.hpp
@@ -0,0 +1,40 @@
+
+#pragma once
+
+#include <cmath>
+#include <functional>
+#include <array>
+
+#define FINITE_STATES
+
+// must have N_STATES, states[N_STATES], and state_to_ind defined before
+// invoking header
+
+std::array<std::array<std::array<double, N_STATES>, N_STATES>, N_STATES> J_probs;
+#ifndef NOFIELD
+std::array<std::array<double, N_STATES>, N_STATES> H_probs;
+#endif
+
+template <class X_t>
+#ifndef NOFIELD
+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++) {
+      H_probs[i][j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T);
+    }
+  }
+#else
+void initialize_probs(std::function <double(X_t, X_t)> J, double T) {
+#endif
+  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][j][k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T);
+      }
+    }
+  }
+}
+
+
+
+
diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp
new file mode 100644
index 0000000..165544a
--- /dev/null
+++ b/lib/include/wolff/graph.hpp
@@ -0,0 +1,21 @@
+
+#pragma once
+
+#include <inttypes.h>
+#include <cmath>
+#include <stdlib.h>
+#include <vector>
+
+#include "types.h"
+
+class graph_t {
+  public:
+    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/include/wolff/state.hpp b/lib/include/wolff/state.hpp
new file mode 100644
index 0000000..1f5e359
--- /dev/null
+++ b/lib/include/wolff/state.hpp
@@ -0,0 +1,89 @@
+
+#pragma once
+
+#include <functional>
+#include <vector>
+
+#include "types.h"
+#include "graph.hpp"
+
+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;
+    v_t nv;
+    v_t ne;
+    graph_t g;
+    double T;
+    std::vector<X_t> spins;
+    R_t R;
+    double E;
+    typename X_t::M_t M; // the "sum" of the spins, like the total magnetization
+    v_t last_cluster_size;
+    std::vector<typename X_t::F_t> ReF;
+    std::vector<typename X_t::F_t> ImF;
+
+    std::function <double(const X_t&, const X_t&)> J;
+#ifndef NOFIELD
+    std::function <double(const X_t&)> H;
+
+    state_t(D_t D, L_t L, double T, std::function <double(const X_t&, const X_t&)> J, std::function <double(const X_t&)> H) : D(D), L(L), g(D, L), T(T), R(), J(J), H(H) {
+#else
+    state_t(D_t D, L_t L, double T, std::function <double(const X_t&, const X_t&)> J) : D(D), L(L), g(D, L), T(T), R(), J(J) {
+#endif
+      nv = g.nv;
+      ne = g.ne;
+      spins.resize(nv);
+#ifndef NOFIELD
+      g.add_ext();
+      E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
+#else
+      E = - (double)ne * J(spins[0], spins[0]);
+#endif
+      M = spins[0] * nv;
+      last_cluster_size = 0;
+      ReF.resize(D);
+      ImF.resize(D);
+      for (D_t i = 0; i < D; i++) {
+        ReF[i] = spins[0] * 0.0;
+        ImF[i] = spins[0] * 0.0;
+      }
+      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_new - s_old;
+    }
+
+    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/include/wolff/types.h b/lib/include/wolff/types.h
new file mode 100644
index 0000000..ec9efaf
--- /dev/null
+++ b/lib/include/wolff/types.h
@@ -0,0 +1,36 @@
+
+#include <inttypes.h>
+
+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;
+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
+#define MAX_H INT_FAST64_MAX
+#define MIN_H INT_FAST64_MIN
+
+#define PRIv PRIuFAST32
+#define PRIq PRIuFAST8
+#define PRIR PRIuFAST16
+#define PRID PRIuFAST8
+#define PRIL PRIuFAST16
+#define PRIcount PRIuFAST64
+#define PRIh PRIdFAST64
+
+#define SCNv SCNuFAST32
+#define SCNq SCNuFAST8
+#define SCNR SCNuFAST16
+#define SCND SCNuFAST8
+#define SCNL SCNuFAST16
+#define SCNcount SCNuFAST64
+#define SCNh SCNdFAST64
+