summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/include/wolff.hpp15
-rw-r--r--lib/include/wolff/cluster.hpp20
-rw-r--r--lib/include/wolff/graph.hpp9
-rw-r--r--lib/include/wolff/meas.h19
-rw-r--r--lib/include/wolff/state.hpp81
-rw-r--r--lib/src/graph.cpp60
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;
+ }
}
}