diff options
Diffstat (limited to 'lib/include/wolff/cluster.hpp')
| -rw-r--r-- | lib/include/wolff/cluster.hpp | 118 |
1 files changed, 60 insertions, 58 deletions
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp index 805e2c3..055cdf3 100644 --- a/lib/include/wolff/cluster.hpp +++ b/lib/include/wolff/cluster.hpp @@ -1,120 +1,122 @@ -#pragma once +#ifndef WOLFF_CLUSTER_H +#define WOLFF_CLUSTER_H #include <random> #include <cmath> #include <vector> -#include <stack> +#include <queue> #include "types.h" -#include "state.hpp" +#include "system.hpp" #include "graph.hpp" -#include "meas.h" +#include "measurement.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, wolff_measurement<R_t, X_t>& m) { +void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r, + std::mt19937& rng, wolff_measurement<R_t, X_t>& A) { std::uniform_real_distribution<double> dist(0.0, 1.0); - std::stack<v_t> stack; - stack.push(v0); + std::queue<v_t> queue; + queue.push(i0); - std::vector<bool> marks(s.g.nv, false); + std::vector<bool> marks(S.G.nv, false); - while (!stack.empty()) { - v_t v = stack.top(); - stack.pop(); + while (!queue.empty()) { + v_t i = queue.front(); + queue.pop(); - if (!marks[v]) { // don't reprocess anyone we've already visited! - marks[v] = true; + if (!marks[i]) { // don't reprocess anyone we've already visited! + marks[i] = true; X_t si_new; -#ifndef NOFIELD - R_t R_new; +#ifndef WOLFF_NO_FIELD + R_t s0_new; - bool v_is_ghost = (v == s.nv); // ghost site has the last index + bool we_are_ghost = (i == S.nv); - if (v_is_ghost) { - R_new = r.act(s.R); // if we are, then we're moving the transformation + if (we_are_ghost) { + s0_new = r.act(S.s0); } else #endif { - si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site + si_new = r.act(S.s[i]); } - for (const v_t &vn : s.g.v_adj[v]) { - double dE, prob; + for (const v_t &j : S.G.adj[i]) { + double dE, p; -#ifndef NOFIELD - bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost +#ifndef WOLFF_NO_FIELD + bool neighbor_is_ghost = (j == S.nv); - if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond - X_t rs_old, rs_new; + if (we_are_ghost || neighbor_is_ghost) { + X_t s0s_old, s0s_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; + if (neighbor_is_ghost) { + non_ghost = i; + s0s_old = S.s0.act_inverse(S.s[i]); + s0s_new = S.s0.act_inverse(si_new); } 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; + non_ghost = j; + s0s_old = S.s0.act_inverse(S.s[j]); + s0s_new = s0_new.act_inverse(S.s[j]); } -#ifdef SITE_DEPENDENCE - dE = s.H(non_ghost, rs_old) - s.H(non_ghost, rs_new); +#ifdef WOLFF_SITE_DEPENDENCE + dE = S.B(non_ghost, s0s_old) - S.B(non_ghost, s0s_new); #else - dE = s.H(rs_old) - s.H(rs_new); + dE = S.B(s0s_old) - S.B(s0s_new); #endif -#ifdef FINITE_STATES - prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)]; +#ifdef WOLFF_FINITE_STATES + p = finite_states_Bp[finite_states_enum(s0s_old)] + [finite_states_enum(s0s_new)]; #endif // run measurement hooks for encountering a ghost bond - m.ghost_bond_added(non_ghost, rs_old, rs_new, dE); + A.ghost_bond_visited(S, non_ghost, s0s_old, s0s_new, dE); } else // this is a perfectly normal bond! #endif { -#ifdef BOND_DEPENDENCE - dE = s.J(v, s.spins[v], vn, s.spins[vn]) - s.J(v, si_new, vn, s.spins[vn]); +#ifdef WOLFF_BOND_DEPENDENCE + dE = S.Z(i, S.s[i], j, S.s[j]) - S.Z(i, si_new, j, S.s[j]); #else - dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]); + dE = S.Z(S.s[i], S.s[j]) - S.Z(si_new, S.s[j]); #endif - -#ifdef FINITE_STATES - prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])]; +#ifdef WOLFF_FINITE_STATES + p = finite_states_Zp[finite_states_enum(S.s[i])] + [finite_states_enum(si_new)] + [finite_states_enum(S.s[j])]; #endif // run measurement hooks for encountering a plain bond - m.plain_bond_added(v, s.spins[v], si_new, vn, s.spins[vn], dE); + A.plain_bond_visited(S, i, si_new, j, dE); } #ifndef FINITE_STATES - prob = 1.0 - exp(-dE / s.T); + p = 1.0 - exp(-dE / S.T); #endif - if (dist(rand) < prob) { - stack.push(vn); // push the neighboring vertex to the stack + if (dist(rng) < p) { + queue.push(j); // push the neighboring vertex to the queue } } -#ifndef NOFIELD - if (v_is_ghost) { - m.ghost_site_transformed(s.R, R_new); - s.R = R_new; +#ifndef WOLFF_NO_FIELD + if (we_are_ghost) { + A.ghost_site_transformed(S, s0_new); + S.s0 = s0_new; } else #endif { - m.plain_site_transformed(v, s.spins[v], si_new); - s.spins[v] = si_new; + A.plain_site_transformed(S, i, si_new); + S.s[i] = si_new; } } } } +#endif + |