many changes, including reworking the measurements system

6 files changed, 91 insertions, 113 deletions

diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp
index c10a211..b730c8d 100644
--- a/lib/include/wolff.hpp
+++ b/lib/include/wolff.hpp
@@ -3,7 +3,7 @@
 #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) {
+void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(std::mt19937&, X_t)> gen_R, wolff_measurement<R_t, X_t>& m, std::mt19937& r) {
 
 #ifdef FINITE_STATES
 #ifdef NOFIELD
@@ -15,21 +15,16 @@ void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(std::mt19937&, X
 
   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);
-  }
+    m.pre_cluster(s, steps, N, v0, step);
+
+    flip_cluster<R_t, X_t>(s, v0, step, r, m);
 
-  if (!silent) {
-    printf("\033[F\033[J");
+    m.post_cluster(s, steps, N);
   }
-  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
index e9dff7b..805e2c3 100644
--- a/lib/include/wolff/cluster.hpp
+++ b/lib/include/wolff/cluster.hpp
@@ -9,11 +9,11 @@
 #include "types.h"
 #include "state.hpp"
 #include "graph.hpp"
+#include "meas.h"
 
 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;
+void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand, wolff_measurement<R_t, X_t>& m) {
+  std::uniform_real_distribution<double> dist(0.0, 1.0);
 
   std::stack<v_t> stack;
   stack.push(v0);
@@ -75,8 +75,8 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand
           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);
+          // run measurement hooks for encountering a ghost bond
+          m.ghost_bond_added(non_ghost, rs_old, rs_new, dE);
         } else // this is a perfectly normal bond!
 #endif
         {
@@ -90,9 +90,10 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand
 #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);
+          // run measurement hooks for encountering a plain bond
+          m.plain_bond_added(v, s.spins[v], si_new, vn, s.spins[vn], dE);
+        }
 
 #ifndef FINITE_STATES
         prob = 1.0 - exp(-dE / s.T);
@@ -105,16 +106,15 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand
 
 #ifndef NOFIELD
       if (v_is_ghost) {
+        m.ghost_site_transformed(s.R, R_new);
         s.R = R_new;
       } else
 #endif
       {
+        m.plain_site_transformed(v, s.spins[v], si_new);
         s.spins[v] = si_new;
-        nv++;
       }
     }
   }
-
-  s.last_cluster_size = nv;
 }
 
diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp
index 165544a..0aeb6af 100644
--- a/lib/include/wolff/graph.hpp
+++ b/lib/include/wolff/graph.hpp
@@ -1,13 +1,16 @@
 
 #pragma once
 
-#include <inttypes.h>
 #include <cmath>
-#include <stdlib.h>
 #include <vector>
 
 #include "types.h"
 
+typedef enum lattice_t {
+  SQUARE_LATTICE,
+  DIAGONAL_LATTICE
+} lattice_t;
+
 class graph_t {
   public:
     v_t ne;
@@ -15,7 +18,7 @@ class graph_t {
     std::vector<std::vector<v_t>> v_adj;
     std::vector<std::vector<double>> coordinate;
 
-    graph_t(D_t D, L_t L);
+    graph_t(D_t D, L_t L, lattice_t lat = SQUARE_LATTICE);
     void add_ext();
 };
 
diff --git a/lib/include/wolff/meas.h b/lib/include/wolff/meas.h
new file mode 100644
index 0000000..4895eac
--- /dev/null
+++ b/lib/include/wolff/meas.h
@@ -0,0 +1,19 @@
+
+#pragma once
+
+#include "types.h"
+
+template <class R_t, class X_t>
+class wolff_measurement {
+  public:
+    virtual void pre_cluster(const state_t<R_t, X_t>&, count_t, count_t, v_t, const R_t&) = 0;
+
+    virtual void plain_bond_added(v_t, const X_t&, const X_t&, v_t, const X_t&, double) = 0;
+    virtual void ghost_bond_added(v_t, const X_t&, const X_t&, double) = 0;
+
+    virtual void plain_site_transformed(v_t, const X_t&, const X_t&) = 0;
+    virtual void ghost_site_transformed(const R_t&, const R_t&) = 0;
+
+    virtual void post_cluster(const state_t<R_t, X_t>&, count_t, count_t) = 0;
+};
+
diff --git a/lib/include/wolff/state.hpp b/lib/include/wolff/state.hpp
index e7c0ac3..4909881 100644
--- a/lib/include/wolff/state.hpp
+++ b/lib/include/wolff/state.hpp
@@ -9,26 +9,17 @@
 
 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; // the number of vertices in the lattice
-    v_t ne; // the number of edges in the lattice
-    graph_t g; // the graph defining the lattice without ghost
+    D_t D; // the dimension of the system
+    L_t L; // the linear size of the lattice
+    v_t nv; // the number of vertices in the original lattice
+    v_t ne; // the number of edges in the original lattice
+    graph_t g; // the graph defining the lattice with ghost
     double T; // the temperature
     std::vector<X_t> spins; // the state of the ordinary spins
 #ifndef NOFIELD
     R_t R; // the current state of the ghost site
 #endif
-    double E; // the system's total energy
-    typename X_t::M_t M; // the "sum" of the spins, like the total magnetization
-    v_t last_cluster_size; // the size of the last cluster
-    std::vector<typename X_t::F_t> ReF;
-    std::vector<typename X_t::F_t> ImF;
 
 #ifdef BOND_DEPENDENCE
     std::function <double(v_t, const X_t&, v_t, const X_t&)> J; // coupling between sites
@@ -57,7 +48,7 @@ class state_t {
         , std::function <double(const X_t&)> H
 #endif
 #endif
-           ) : D(D), L(L), g(D, L), T(T),
+        , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), g(D, L, lat), T(T),
 #ifndef NOFIELD
               R(),
 #endif
@@ -69,69 +60,9 @@ class state_t {
       nv = g.nv;
       ne = g.ne;
       spins.resize(nv);
-#ifdef BOND_DEPENDENCE
-      E = 0;
-      for (v_t v = 0; v < nv; v++) {
-        for (const v_t &vn : g.v_adj[v]) {
-          if (v < vn) {
-            E -= J(v, spins[v], vn, spins[vn]);
-          }
-        }
-      }
-#else
-      E = - (double)ne * J(spins[0], spins[0]);
-#endif
-
 #ifndef NOFIELD
       g.add_ext();
-#ifdef SITE_DEPENDENCE
-      for (v_t i = 0; i < nv; i++) {
-        E -= H(i, spins[i]);
-      }
-#else
-      E -= (double)nv * H(spins[0]);
-#endif
-#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/src/graph.cpp b/lib/src/graph.cpp
index 4043413..c6f0ba6 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -1,25 +1,55 @@
 
 #include <wolff/graph.hpp>
 
-graph_t::graph_t(D_t D, L_t L) {
-  nv = pow(L, D);
-  ne = D * nv;
+graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
+  switch (lat) {
+    case SQUARE_LATTICE: {
+        nv = pow(L, D);
+        ne = D * nv;
 
-  v_adj.resize(nv);
-  coordinate.resize(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 (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;
+        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)));
+          }
+        }
+        break;
+      }
+    case DIAGONAL_LATTICE: {
+        nv = D * pow(L, D);
+        ne = D * nv;
+
+        v_adj.resize(nv);
+        coordinate.resize(nv);
+
+        for (std::vector<v_t> v_adj_i : v_adj) {
+          v_adj_i.reserve(4 * (D - 1));
+        }
+
+        for (D_t i = 0; i < D; i++) {
+          v_t sb = i * pow(L, D);
+
+          for (v_t j = 0; j < pow(L, D); j++) {
+            v_t vc = sb + j;
 
-      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)));
-    }
+            v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + j);
+            v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L);
+            v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
+            v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
+          }
+        }
+        break;
+      }
   }
 }