diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/include/wolff.hpp | 21 | ||||
| -rw-r--r-- | lib/include/wolff/cluster.hpp | 57 | ||||
| -rw-r--r-- | lib/include/wolff/finite_states.hpp | 6 | ||||
| -rw-r--r-- | lib/include/wolff/graph.hpp | 21 | ||||
| -rw-r--r-- | lib/include/wolff/measurement.hpp | 18 | ||||
| -rw-r--r-- | lib/include/wolff/models/ising.hpp | 11 | ||||
| -rw-r--r-- | lib/include/wolff/models/orthogonal.hpp | 15 | ||||
| -rw-r--r-- | lib/include/wolff/models/potts.hpp | 2 | ||||
| -rw-r--r-- | lib/include/wolff/system.hpp | 30 | ||||
| -rw-r--r-- | lib/include/wolff/types.h | 2 | ||||
| -rw-r--r-- | lib/src/graph.cpp | 79 |
11 files changed, 139 insertions, 123 deletions
diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp index b78cd4b..7e909c8 100644 --- a/lib/include/wolff.hpp +++ b/lib/include/wolff.hpp @@ -3,25 +3,30 @@ #define WOLFF_H #include "wolff/cluster.hpp" +#include "wolff/measurement.hpp" + +namespace wolff{ template <class R_t, class X_t> -void wolff(N_t N, wolff_system<R_t, X_t>& S, - std::function <R_t(std::mt19937&, X_t)> r_gen, - wolff_measurement<R_t, X_t>& A, std::mt19937& rng) { +void system<R_t, X_t>::run_wolff(N_t N, + std::function <R_t(std::mt19937&, const system<R_t, X_t>&, v_t)> r_gen, + measurement<R_t, X_t>& A, std::mt19937& rng) { - std::uniform_int_distribution<v_t> dist(0, S.nv - 1); + std::uniform_int_distribution<v_t> dist(0, nv - 1); for (N_t n = 0; n < N; n++) { v_t i0 = dist(rng); - R_t r = r_gen(rng, S.s[i0]); + R_t r = r_gen(rng, *this, i0); - A.pre_cluster(n, N, S, i0, r); + A.pre_cluster(n, N, *this, i0, r); - wolff_cluster_flip<R_t, X_t>(S, i0, r, rng, A); + this->flip_cluster(i0, r, rng, A); - A.post_cluster(n, N, S); + A.post_cluster(n, N, *this); } } +} + #endif diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp index 055cdf3..b66f367 100644 --- a/lib/include/wolff/cluster.hpp +++ b/lib/include/wolff/cluster.hpp @@ -7,20 +7,19 @@ #include <vector> #include <queue> -#include "types.h" #include "system.hpp" -#include "graph.hpp" -#include "measurement.hpp" + +namespace wolff { template <class R_t, class X_t> -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) { +void system<R_t, X_t>::flip_cluster(v_t i0, const R_t& r, + std::mt19937& rng, measurement<R_t, X_t>& A) { std::uniform_real_distribution<double> dist(0.0, 1.0); std::queue<v_t> queue; queue.push(i0); - std::vector<bool> marks(S.G.nv, false); + std::vector<bool> marks(G.nv, false); while (!queue.empty()) { v_t i = queue.front(); @@ -33,21 +32,21 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r, #ifndef WOLFF_NO_FIELD R_t s0_new; - bool we_are_ghost = (i == S.nv); + bool we_are_ghost = (i == nv); if (we_are_ghost) { - s0_new = r.act(S.s0); + s0_new = r.act(s0); } else #endif { - si_new = r.act(S.s[i]); + si_new = r.act(s[i]); } - for (const v_t &j : S.G.adj[i]) { + for (const v_t &j : G.adjacency[i]) { double dE, p; #ifndef WOLFF_NO_FIELD - bool neighbor_is_ghost = (j == S.nv); + bool neighbor_is_ghost = (j == nv); if (we_are_ghost || neighbor_is_ghost) { X_t s0s_old, s0s_new; @@ -55,18 +54,18 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r, if (neighbor_is_ghost) { non_ghost = i; - s0s_old = S.s0.act_inverse(S.s[i]); - s0s_new = S.s0.act_inverse(si_new); + s0s_old = s0.act_inverse(s[i]); + s0s_new = s0.act_inverse(si_new); } else { non_ghost = j; - s0s_old = S.s0.act_inverse(S.s[j]); - s0s_new = s0_new.act_inverse(S.s[j]); + s0s_old = s0.act_inverse(s[j]); + s0s_new = s0_new.act_inverse(s[j]); } #ifdef WOLFF_SITE_DEPENDENCE - dE = S.B(non_ghost, s0s_old) - S.B(non_ghost, s0s_new); + dE = B(non_ghost, s0s_old) - B(non_ghost, s0s_new); #else - dE = S.B(s0s_old) - S.B(s0s_new); + dE = B(s0s_old) - B(s0s_new); #endif #ifdef WOLFF_FINITE_STATES @@ -75,28 +74,28 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r, #endif // run measurement hooks for encountering a ghost bond - A.ghost_bond_visited(S, non_ghost, s0s_old, s0s_new, dE); + A.ghost_bond_visited(*this, non_ghost, s0s_old, s0s_new, dE); } else // this is a perfectly normal bond! #endif { #ifdef WOLFF_BOND_DEPENDENCE - dE = S.Z(i, S.s[i], j, S.s[j]) - S.Z(i, si_new, j, S.s[j]); + dE = Z(i, s[i], j, s[j]) - Z(i, si_new, j, s[j]); #else - dE = S.Z(S.s[i], S.s[j]) - S.Z(si_new, S.s[j]); + dE = Z(s[i], s[j]) - Z(si_new, s[j]); #endif #ifdef WOLFF_FINITE_STATES - p = finite_states_Zp[finite_states_enum(S.s[i])] + p = finite_states_Zp[finite_states_enum(s[i])] [finite_states_enum(si_new)] - [finite_states_enum(S.s[j])]; + [finite_states_enum(s[j])]; #endif // run measurement hooks for encountering a plain bond - A.plain_bond_visited(S, i, si_new, j, dE); + A.plain_bond_visited(*this, i, si_new, j, dE); } #ifndef FINITE_STATES - p = 1.0 - exp(-dE / S.T); + p = 1.0 - exp(-dE / T); #endif if (dist(rng) < p) { @@ -106,17 +105,19 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r, #ifndef WOLFF_NO_FIELD if (we_are_ghost) { - A.ghost_site_transformed(S, s0_new); - S.s0 = s0_new; + A.ghost_site_transformed(*this, s0_new); + s0 = s0_new; } else #endif { - A.plain_site_transformed(S, i, si_new); - S.s[i] = si_new; + A.plain_site_transformed(*this, i, si_new); + s[i] = si_new; } } } } +} + #endif diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp index 4cca69e..7e54eff 100644 --- a/lib/include/wolff/finite_states.hpp +++ b/lib/include/wolff/finite_states.hpp @@ -8,13 +8,15 @@ #include "system.hpp" +namespace wolff { + std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Zp; #ifndef WOLFF_NO_FIELD std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Bp; #endif template <class R_t, class X_t> -void finite_states_init(const wolff_system<R_t, X_t>& S) { +void finite_states_init(const system<R_t, X_t>& S) { #ifndef WOLFF_NO_FIELD for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) { for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) { @@ -31,5 +33,7 @@ void finite_states_init(const wolff_system<R_t, X_t>& S) { } } +} + #endif diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp index fcd73fa..65b3941 100644 --- a/lib/include/wolff/graph.hpp +++ b/lib/include/wolff/graph.hpp @@ -5,23 +5,26 @@ #include <cmath> #include <vector> -#include "types.h" +namespace wolff { -typedef enum lattice_t { - SQUARE_LATTICE, - DIAGONAL_LATTICE -} lattice_t; +#include "types.h" -class graph_t { +class graph { public: + D_t D; + L_t L; v_t ne; v_t nv; - std::vector<std::vector<v_t>> adj; + std::vector<std::vector<v_t>> adjacency; std::vector<std::vector<double>> coordinate; - graph_t(D_t D, L_t L, lattice_t lat = SQUARE_LATTICE); - void add_ext(); + graph(); + graph(D_t D, L_t L); + + void add_ghost(); }; +} + #endif diff --git a/lib/include/wolff/measurement.hpp b/lib/include/wolff/measurement.hpp index a6a5f79..604eaf5 100644 --- a/lib/include/wolff/measurement.hpp +++ b/lib/include/wolff/measurement.hpp @@ -4,21 +4,25 @@ #include "system.hpp" +namespace wolff { + template <class R_t, class X_t> -class wolff_measurement { +class measurement { public: - virtual void pre_cluster(N_t, N_t, const wolff_system<R_t, X_t>&, v_t, const R_t&) = 0; + virtual void pre_cluster(N_t, N_t, const system<R_t, X_t>&, v_t, const R_t&) = 0; - virtual void plain_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, v_t, double) = 0; - virtual void plain_site_transformed(const wolff_system<R_t, X_t>&, v_t, const X_t&) = 0; + virtual void plain_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, v_t, double) = 0; + virtual void plain_site_transformed(const system<R_t, X_t>&, v_t, const X_t&) = 0; #ifndef WOLFF_NO_FIELD - virtual void ghost_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) = 0; - virtual void ghost_site_transformed(const wolff_system<R_t, X_t>&, const R_t&) = 0; + virtual void ghost_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) = 0; + virtual void ghost_site_transformed(const system<R_t, X_t>&, const R_t&) = 0; #endif - virtual void post_cluster(N_t, N_t, const wolff_system<R_t, X_t>&) = 0; + virtual void post_cluster(N_t, N_t, const system<R_t, X_t>&) = 0; }; +} + #endif diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp index c6a3a47..824c063 100644 --- a/lib/include/wolff/models/ising.hpp +++ b/lib/include/wolff/models/ising.hpp @@ -3,6 +3,9 @@ #define WOLFF_MODELS_ISING #include "../types.h" +#include "../system.hpp" + +namespace wolff { class ising_t { public: @@ -73,9 +76,15 @@ inline ising_t::F_t operator*(double a, const ising_t& s) { return s * a; } +ising_t gen_ising(std::mt19937&, const system<ising_t, ising_t>&, v_t) { + return ising_t(true); +}; + #define WOLFF_FINITE_STATES_N 2 const ising_t finite_states_possible[2] = {ising_t(0), ising_t(1)}; -q_t finite_states_enum(ising_t state) { return (q_t)state.x; } +q_t finite_states_enum(const ising_t& state) { return (q_t)state.x; } + +} #endif diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp index 9dd5ddd..58203b7 100644 --- a/lib/include/wolff/models/orthogonal.hpp +++ b/lib/include/wolff/models/orthogonal.hpp @@ -4,7 +4,8 @@ #include <random> #include <cmath> -#include <wolff/types.h> +#include "../types.h" +#include "../system.hpp" #include "vector.hpp" template <q_t q, class T> @@ -113,7 +114,7 @@ class orthogonal_t : public std::array<std::array<T, q>, q> { }; template <q_t q> -orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const vector_t <q, double>& v) { +orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const system<orthogonal_t<q, double>, vector_t<q, double>>&, v_t) { std::normal_distribution<double> dist(0.0,1.0); orthogonal_t <q, double> ptr; ptr.is_reflection = true; @@ -135,7 +136,7 @@ orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const vecto } template <q_t q> -orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const vector_t <q, double>& v0, double epsilon, unsigned int n) { +orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const system<orthogonal_t<q, double>, vector_t<q, double>>& S, v_t i0, double epsilon, unsigned int n) { std::normal_distribution<double> dist(0.0,1.0); orthogonal_t <q, double> m; m.is_reflection = true; @@ -149,14 +150,14 @@ orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const double cosr = cos(2 * M_PI * rotation / (double)n / 2.0); double sinr = sin(2 * M_PI * rotation / (double)n / 2.0); - v[0] = v0[0] * cosr - v0[1] * sinr; - v[1] = v0[1] * cosr + v0[0] * sinr; + v[0] = S.s[i0][0] * cosr - S.s[i0][1] * sinr; + v[1] = S.s[i0][1] * cosr + S.s[i0][0] * sinr; for (q_t i = 2; i < q; i++) { - v[i] = v0[i]; + v[i] = S.s[i0][i]; } } else { - v = v0; + v = S.s[i0]; } double m_dot_v = 0; diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp index 903f25f..10b440e 100644 --- a/lib/include/wolff/models/potts.hpp +++ b/lib/include/wolff/models/potts.hpp @@ -46,5 +46,5 @@ class potts_t { const potts_t<POTTSQ> states[256] = {{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {21}, {22}, {23}, {24}, {25}, {26}, {27}, {28}, {29}, {30}, {31}, {32}, {33}, {34}, {35}, {36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46}, {47}, {48}, {49}, {50}, {51}, {52}, {53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}, {61}, {62}, {63}, {64}, {65}, {66}, {67}, {68}, {69}, {70}, {71}, {72}, {73}, {74}, {75}, {76}, {77}, {78}, {79}, {80}, {81}, {82}, {83}, {84}, {85}, {86}, {87}, {88}, {89}, {90}, {91}, {92}, {93}, {94}, {95}, {96}, {97}, {98}, {99}, {100}, {101}, {102}, {103}, {104}, {105}, {106}, {107}, {108}, {109}, {110}, {111}, {112}, {113}, {114}, {115}, {116}, {117}, {118}, {119}, {120}, {121}, {122}, {123}, {124}, {125}, {126}, {127}, {128}, {129}, {130}, {131}, {132}, {133}, {134}, {135}, {136}, {137}, {138}, {139}, {140}, {141}, {142}, {143}, {144}, {145}, {146}, {147}, {148}, {149}, {150}, {151}, {152}, {153}, {154}, {155}, {156}, {157}, {158}, {159}, {160}, {161}, {162}, {163}, {164}, {165}, {166}, {167}, {168}, {169}, {170}, {171}, {172}, {173}, {174}, {175}, {176}, {177}, {178}, {179}, {180}, {181}, {182}, {183}, {184}, {185}, {186}, {187}, {188}, {189}, {190}, {191}, {192}, {193}, {194}, {195}, {196}, {197}, {198}, {199}, {200}, {201}, {202}, {203}, {204}, {205}, {206}, {207}, {208}, {209}, {210}, {211}, {212}, {213}, {214}, {215}, {216}, {217}, {218}, {219}, {220}, {221}, {222}, {223}, {224}, {225}, {226}, {227}, {228}, {229}, {230}, {231}, {232}, {233}, {234}, {235}, {236}, {237}, {238}, {239}, {240}, {241}, {242}, {243}, {244}, {245}, {246}, {247}, {248}, {249}, {250}, {251}, {252}, {253}, {254}, {255}}; template <q_t q> -q_t finite_states_enum(potts_t<q> state) { return (q_t)state.x; } +q_t finite_states_enum(const potts_t<q>& state) { return (q_t)state.x; } diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp index ac82f44..32ad38e 100644 --- a/lib/include/wolff/system.hpp +++ b/lib/include/wolff/system.hpp @@ -4,23 +4,25 @@ #include <functional> #include <vector> +#include <random> -#include "types.h" #include "graph.hpp" +namespace wolff { + +#include "types.h" + +template <class R_t, class X_t> +class measurement; + template <class R_t, class X_t> -class wolff_system { +class system { public: - D_t D; // dimension - L_t L; // linear size v_t nv; // number of vertices v_t ne; // number of edges - graph_t G; // the graph defining the lattice with ghost + graph G; // the graph defining the lattice with ghost double T; // the temperature std::vector<X_t> s; // the state of the ordinary spins -#ifndef WOLFF_NO_FIELD - R_t s0; // the current state of the ghost site -#endif #ifdef WOLFF_BOND_DEPENDENCE std::function <double(v_t, const X_t&, v_t, const X_t&)> Z; // coupling between sites @@ -29,6 +31,7 @@ class wolff_system { #endif #ifndef WOLFF_NO_FIELD + R_t s0; // the current state of the ghost site #ifdef WOLFF_SITE_DEPENDENCE std::function <double(v_t, const X_t&)> B; // coupling with the external field #else @@ -36,7 +39,7 @@ class wolff_system { #endif #endif - wolff_system(D_t D, L_t L, double T, + system(graph g, double T, #ifdef WOLFF_BOND_DEPENDENCE std::function <double(v_t, const X_t&, v_t, const X_t&)> Z #else @@ -49,7 +52,7 @@ class wolff_system { , std::function <double(const X_t&)> B #endif #endif - , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), G(D, L, lat), T(T), Z(Z) + ) : G(g), T(T), Z(Z) #ifndef WOLFF_NO_FIELD , s0(), B(B) #endif @@ -58,13 +61,18 @@ class wolff_system { ne = G.ne; s.resize(nv); #ifndef WOLFF_NO_FIELD - G.add_ext(); + G.add_ghost(); #endif #ifdef WOLFF_FINITE_STATES finite_states_init(*this); #endif } + + void flip_cluster(v_t, const R_t&, std::mt19937&, measurement<R_t, X_t>&); + void run_wolff(N_t, std::function <R_t(std::mt19937&, const system<R_t, X_t>&, v_t)> r_gen, measurement<R_t, X_t>& A, std::mt19937& rng); }; +} + #endif diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h index 68bd09e..1bc8c4d 100644 --- a/lib/include/wolff/types.h +++ b/lib/include/wolff/types.h @@ -8,7 +8,7 @@ typedef uint_fast16_t L_t; // linear size typedef uint_fast64_t N_t; // cycle iterator #define MAX_v UINT_FAST32_MAX -#define MAX_Q UINT_FAST8_MAX +#define MAX_q UINT_FAST8_MAX #define MAX_D UINT_FAST8_MAX #define MAX_L UINT_FAST16_MAX #define MAX_N UINT_FAST64_MAX diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index b42aa78..4f91be4 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -1,74 +1,55 @@ #include <wolff/graph.hpp> -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; +namespace wolff { - adj.resize(nv); - coordinate.resize(nv); - - for (std::vector<v_t> adj_i : adj) { - 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; - - adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1))); - 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 = 2 * nv; +graph::graph() { + D = 0; + L = 0; + nv = 0; + ne = 0; +} - adj.resize(nv); - coordinate.resize(nv); +graph::graph(D_t D, L_t L) : D(D), L(L) { + nv = pow(L, D); + ne = D * nv; - for (std::vector<v_t> adj_i : adj) { - adj_i.reserve(4 * (D - 1)); - } + adjacency.resize(nv); + coordinate.resize(nv); - for (D_t i = 0; i < D; i++) { - v_t sb = i * pow(L, D); + for (std::vector<v_t> adj_i : adjacency) { + adj_i.reserve(2 * D); + } - for (v_t j = 0; j < pow(L, D); j++) { - v_t vc = sb + j; + 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; - adj[vc].push_back(((i + 1) % D) * pow(L, D) + j); - 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); - 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); - 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; - } + adjacency[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1))); + adjacency[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))); + } } } -void graph_t::add_ext() { - for (std::vector<v_t>& adj_i : adj) { +void graph::add_ghost() { + for (std::vector<v_t>& adj_i : adjacency) { adj_i.push_back(nv); } - adj.resize(nv + 1); + adjacency.resize(nv + 1); coordinate.resize(nv + 1); - adj[nv].reserve(nv); + adjacency[nv].reserve(nv); for (v_t i = 0; i < nv; i++) { - adj[nv].push_back(i); + adjacency[nv].push_back(i); } coordinate[nv].resize(coordinate[0].size()); ne += nv; - nv += 1; + nv++; +} + } |