From 49ac78a6c04e215950bc9c0f97368605e63da15b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 14 Jan 2019 15:47:59 -0500 Subject: Large refactoring around changes in the graph class. - Graphs now use lists of references instead of vectors of indicies. - Vertices and edges have associated classes that can be given arbitrary properties via template specification. - All essential library headers have been combined into one, wolff.hpp. --- lib/CMakeLists.txt | 16 +- lib/include/wolff.hpp | 32 --- lib/include/wolff/cluster.hpp | 120 ----------- lib/include/wolff/graph.hpp | 30 --- lib/include/wolff/measurement.hpp | 28 --- lib/include/wolff/models/dihedral.hpp | 56 ------ lib/include/wolff/models/dihedral_inf.hpp | 57 ------ lib/include/wolff/models/height.hpp | 57 ------ lib/include/wolff/models/ising.hpp | 93 --------- lib/include/wolff/models/orthogonal.hpp | 208 ------------------- lib/include/wolff/models/potts.hpp | 68 ------- lib/include/wolff/models/symmetric.hpp | 59 ------ lib/include/wolff/models/vector.hpp | 115 ----------- lib/include/wolff/system.hpp | 105 ---------- lib/include/wolff/types.h | 32 --- lib/src/graph.cpp | 55 ----- lib/wolff.hpp | 320 ++++++++++++++++++++++++++++++ lib/wolff_models/dihedral.hpp | 54 +++++ lib/wolff_models/dihedral_inf.hpp | 55 +++++ lib/wolff_models/height.hpp | 56 ++++++ lib/wolff_models/ising.hpp | 92 +++++++++ lib/wolff_models/orthogonal.hpp | 206 +++++++++++++++++++ lib/wolff_models/potts.hpp | 66 ++++++ lib/wolff_models/symmetric.hpp | 57 ++++++ lib/wolff_models/vector.hpp | 113 +++++++++++ 25 files changed, 1025 insertions(+), 1125 deletions(-) delete mode 100644 lib/include/wolff.hpp delete mode 100644 lib/include/wolff/cluster.hpp delete mode 100644 lib/include/wolff/graph.hpp delete mode 100644 lib/include/wolff/measurement.hpp delete mode 100644 lib/include/wolff/models/dihedral.hpp delete mode 100644 lib/include/wolff/models/dihedral_inf.hpp delete mode 100644 lib/include/wolff/models/height.hpp delete mode 100644 lib/include/wolff/models/ising.hpp delete mode 100644 lib/include/wolff/models/orthogonal.hpp delete mode 100644 lib/include/wolff/models/potts.hpp delete mode 100644 lib/include/wolff/models/symmetric.hpp delete mode 100644 lib/include/wolff/models/vector.hpp delete mode 100644 lib/include/wolff/system.hpp delete mode 100644 lib/include/wolff/types.h delete mode 100644 lib/src/graph.cpp create mode 100644 lib/wolff.hpp create mode 100644 lib/wolff_models/dihedral.hpp create mode 100644 lib/wolff_models/dihedral_inf.hpp create mode 100644 lib/wolff_models/height.hpp create mode 100644 lib/wolff_models/ising.hpp create mode 100644 lib/wolff_models/orthogonal.hpp create mode 100644 lib/wolff_models/potts.hpp create mode 100644 lib/wolff_models/symmetric.hpp create mode 100644 lib/wolff_models/vector.hpp (limited to 'lib') diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 1df9807..445e68f 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -1,14 +1,10 @@ -project(libwolff LANGUAGES C CXX) +project(libwolff LANGUAGES CXX) -add_library(wolff SHARED src/graph.cpp) +add_library(wolff INTERFACE) +target_include_directories(wolff INTERFACE .) +target_compile_features(wolff INTERFACE cxx_std_17) -target_include_directories(wolff PUBLIC $ - $) - -install(TARGETS wolff EXPORT wolffConfig - ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} - LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} - RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) -install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) +install(FILES "wolff.hpp" DESTINATION include) +install(DIRECTORY "wolff_models" DESTINATION include) diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp deleted file mode 100644 index 7e909c8..0000000 --- a/lib/include/wolff.hpp +++ /dev/null @@ -1,32 +0,0 @@ - -#ifndef WOLFF_H -#define WOLFF_H - -#include "wolff/cluster.hpp" -#include "wolff/measurement.hpp" - -namespace wolff{ - -template -void system::run_wolff(N_t N, - std::function &, v_t)> r_gen, - measurement& A, std::mt19937& rng) { - - std::uniform_int_distribution dist(0, nv - 1); - - for (N_t n = 0; n < N; n++) { - v_t i0 = dist(rng); - R_t r = r_gen(rng, *this, i0); - - A.pre_cluster(n, N, *this, i0, r); - - this->flip_cluster(i0, r, rng, A); - - A.post_cluster(n, N, *this); - } -} - -} - -#endif - diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp deleted file mode 100644 index 3912446..0000000 --- a/lib/include/wolff/cluster.hpp +++ /dev/null @@ -1,120 +0,0 @@ - -#ifndef WOLFF_CLUSTER_H -#define WOLFF_CLUSTER_H - -#include -#include -#include -#include - -#include "system.hpp" - -namespace wolff { - -template -void system::flip_cluster(v_t i0, const R_t& r, - std::mt19937& rng, measurement& A) { - std::uniform_real_distribution dist(0.0, 1.0); - - std::queue queue; - queue.push(i0); - - std::vector marks(G.nv, false); - - while (!queue.empty()) { - v_t i = queue.front(); - queue.pop(); - - if (!marks[i]) { // don't reprocess anyone we've already visited! - marks[i] = true; - - X_t si_new; -#ifndef WOLFF_NO_FIELD - R_t s0_new; - - bool we_are_ghost = (i == nv); - - if (we_are_ghost) { - s0_new = r.act(s0); - } else -#endif - { - si_new = r.act(s[i]); - } - - for (const v_t &j : G.adjacency[i]) { - double dE, p; - -#ifndef WOLFF_NO_FIELD - bool neighbor_is_ghost = (j == nv); - - if (we_are_ghost || neighbor_is_ghost) { - X_t s0s_old, s0s_new; - v_t non_ghost; - - if (neighbor_is_ghost) { - non_ghost = i; - s0s_old = s0.act_inverse(s[i]); - s0s_new = s0.act_inverse(si_new); - } else { - non_ghost = j; - s0s_old = s0.act_inverse(s[j]); - s0s_new = s0_new.act_inverse(s[j]); - } - -#ifdef WOLFF_SITE_DEPENDENCE - dE = B(non_ghost, s0s_old) - B(non_ghost, s0s_new); -#else - dE = B(s0s_old) - B(s0s_new); -#endif - -#ifdef WOLFF_USE_FINITE_STATES - p = Bp[s0s_old.enumerate()][s0s_new.enumerate()]; -#endif - - // run measurement hooks for encountering a ghost bond - 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 = Z(i, s[i], j, s[j]) - Z(i, si_new, j, s[j]); -#else - dE = Z(s[i], s[j]) - Z(si_new, s[j]); -#endif - -#ifdef WOLFF_USE_FINITE_STATES - p = Zp[s[i].enumerate()][si_new.enumerate()][s[j].enumerate()]; -#endif - - // run measurement hooks for encountering a plain bond - A.plain_bond_visited(*this, i, si_new, j, dE); - } - -#ifndef WOLFF_USE_FINITE_STATES - p = 1.0 - exp(-dE / T); -#endif - - if (dist(rng) < p) { - queue.push(j); // push the neighboring vertex to the queue - } - } - -#ifndef WOLFF_NO_FIELD - if (we_are_ghost) { - A.ghost_site_transformed(*this, s0_new); - s0 = s0_new; - } else -#endif - { - A.plain_site_transformed(*this, i, si_new); - s[i] = si_new; - } - } - } -} - -} - -#endif - diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp deleted file mode 100644 index 65b3941..0000000 --- a/lib/include/wolff/graph.hpp +++ /dev/null @@ -1,30 +0,0 @@ - -#ifndef WOLFF_GRAPH_H -#define WOLFF_GRAPH_H - -#include -#include - -namespace wolff { - -#include "types.h" - -class graph { - public: - D_t D; - L_t L; - v_t ne; - v_t nv; - std::vector> adjacency; - std::vector> coordinate; - - 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 deleted file mode 100644 index 6e352a5..0000000 --- a/lib/include/wolff/measurement.hpp +++ /dev/null @@ -1,28 +0,0 @@ - -#ifndef WOLFF_MEASUREMENTS_H -#define WOLFF_MEASUREMENTS_H - -#include "system.hpp" - -namespace wolff { - -template -class measurement { - public: - virtual void pre_cluster(N_t, N_t, const system&, v_t, const R_t&) {}; - - virtual void plain_bond_visited(const system&, v_t, const X_t&, v_t, double) {}; - virtual void plain_site_transformed(const system&, v_t, const X_t&) {}; - -#ifndef WOLFF_NO_FIELD - virtual void ghost_bond_visited(const system&, v_t, const X_t&, const X_t&, double) {}; - virtual void ghost_site_transformed(const system&, const R_t&) {}; -#endif - - virtual void post_cluster(N_t, N_t, const system&) {}; -}; - -} - -#endif - diff --git a/lib/include/wolff/models/dihedral.hpp b/lib/include/wolff/models/dihedral.hpp deleted file mode 100644 index cf4b067..0000000 --- a/lib/include/wolff/models/dihedral.hpp +++ /dev/null @@ -1,56 +0,0 @@ - -#ifndef WOLFF_MODELS_DIHEDRAL_H -#define WOLFF_MODELS_DIHEDRAL_H - -#include "potts.hpp" - -namespace wolff { - -#include "../types.h" - -template -class dihedral_t { - public: - bool is_reflection; - q_t x; - - dihedral_t() : is_reflection(false), x(0) {} - dihedral_t(bool x, q_t y) : is_reflection(x), x(y) {} - - potts_t act(const potts_t& s) const { - if (this->is_reflection) { - return potts_t(((q + this->x) - s.x) % q); - } else { - return potts_t((this->x + s.x) % q); - } - } - - dihedral_t act(dihedral_t r) const { - if (this->is_reflection) { - return dihedral_t(!(r.is_reflection), ((q + this->x) - r.x) % q); - } else { - return dihedral_t(r.is_reflection, (this->x + r.x) % q); - } - } - - potts_t act_inverse(potts_t s) const { - if (this->is_reflection) { - return this->act(s); - } else { - return potts_t(((s.x + q) - this->x) % q); - } - } - - dihedral_t act_inverse(dihedral_t r) const { - if (this->is_reflection) { - return this->act(r); - } else { - return dihedral_t(r.is_reflection, ((r.x + q) - this->x) % q); - } - } -}; - -} - -#endif - diff --git a/lib/include/wolff/models/dihedral_inf.hpp b/lib/include/wolff/models/dihedral_inf.hpp deleted file mode 100644 index 33d5cec..0000000 --- a/lib/include/wolff/models/dihedral_inf.hpp +++ /dev/null @@ -1,57 +0,0 @@ - -#ifndef WOLFF_MODELS_DIHEDRAL_INF_H -#define WOLFF_MODELS_DIHEDRAL_INF_H - -#include -#include "height.hpp" - -namespace wolff { - -#include "../types.h" - -template -class dihedral_inf_t { - public: - bool is_reflection; - T x; - - dihedral_inf_t() : is_reflection(false), x(0) {} - dihedral_inf_t(bool x, T y) : is_reflection(x), x(y) {} - - height_t act(const height_t& h) const { - if (this->is_reflection) { - return height_t(this->x - h.x); - } else { - return height_t(this->x + h.x); - } - } - - dihedral_inf_t act(const dihedral_inf_t& r) const { - if (this->is_reflection) { - return dihedral_inf_t(!r.is_reflection, this->x - r.x); - } else { - return dihedral_inf_t(r.is_reflection, this->x + r.x); - } - } - - height_t act_inverse(const height_t& h) const { - if (this->is_reflection) { - return this->act(h); - } else { - return height_t(h.x - this->x); - } - } - - dihedral_inf_t act_inverse(const dihedral_inf_t& r) const { - if (this->is_reflection) { - return this->act(r); - } else { - return dihedral_inf_t(r.is_reflection, r.x - this->x); - } - } -}; - -} - -#endif - diff --git a/lib/include/wolff/models/height.hpp b/lib/include/wolff/models/height.hpp deleted file mode 100644 index 9d7e87c..0000000 --- a/lib/include/wolff/models/height.hpp +++ /dev/null @@ -1,57 +0,0 @@ - -#pragma once - -#include -#include - -namespace wolff { - -template -struct height_t { - T x; - - height_t() : x(0) {} - height_t(T x) : x(x) {} - - typedef T M_t; - typedef double F_t; - - inline T operator*(v_t a) const { - return x * a; - } - - inline double operator*(double a) const { - return x * a; - } - - inline T operator-(const height_t& h) const { - return x - h.x; - } -}; - -template -inline typename height_t::M_t operator*(v_t a, height_t h) { - return h * a; -} - -template -inline typename height_t::F_t operator*(double a, height_t h) { - return h * a; -} - -template -inline T& operator+=(T& M, const height_t &h) { - M += h.x; - - return M; -} - -template -inline T& operator-=(T& M, const height_t &h) { - M -= h.x; - - return M; -} - -} - diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp deleted file mode 100644 index d2f6f91..0000000 --- a/lib/include/wolff/models/ising.hpp +++ /dev/null @@ -1,93 +0,0 @@ - -#ifndef WOLFF_MODELS_ISING_H -#define WOLFF_MODELS_ISING_H - -#define WOLFF_FINITE_STATES_N 2 - -#include "../system.hpp" - -namespace wolff { - -#include "../types.h" - -class ising_t { - public: - bool x; - - ising_t() : x(false) {} - - ising_t(bool x) : x(x) {} - ising_t(q_t x) : x((bool)x) {} - - ising_t act(const ising_t& s) const { - if (x) { - return ising_t(!s.x); - } else { - return ising_t(s.x); - } - } - - ising_t act_inverse(const ising_t& s) const { - return this->act(s); - } - - typedef int M_t; - typedef double F_t; - - inline M_t operator*(v_t a) const { - if (x) { - return -(M_t)a; - } else { - return (M_t)a; - } - } - - inline F_t operator*(double a) const { - if (x) { - return -a; - } else { - return a; - } - } - - inline M_t operator-(const ising_t &s) const { - if (x == s.x) { - return 0; - } else { - if (x) { - return -2; - } else { - return 2; - } - } - } - - inline int operator*(const ising_t& s) const { - if (x == s.x) { - return 1; - } else { - return -1; - } - } - - q_t enumerate() const { - return (q_t)x; - } -}; - -inline ising_t::M_t operator*(v_t a, const ising_t& s) { - return s * a; -} - -inline ising_t::F_t operator*(double a, const ising_t& s) { - return s * a; -} - -ising_t gen_ising(std::mt19937&, const system&, v_t) { - return ising_t(true); -}; - -} - -#endif - diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp deleted file mode 100644 index 514c88a..0000000 --- a/lib/include/wolff/models/orthogonal.hpp +++ /dev/null @@ -1,208 +0,0 @@ - -#ifndef WOLFF_MODELS_ORTHOGONAL_H -#define WOLFF_MODELS_ORTHOGONAL_H - -#include -#include - -#include "../system.hpp" -#include "vector.hpp" - -namespace wolff { - -#include "../types.h" - -template -class orthogonal_t : public std::array, q> { - public : - bool is_reflection; - - orthogonal_t() : is_reflection(false) { - for (q_t i = 0; i < q; i++) { - (*this)[i].fill(0); - (*this)[i][i] = (T)1; - } - } - - vector_t act(const vector_t & v) const { - vector_t v_rot; - v_rot.fill(0); - - if (is_reflection) { - double prod = 0; - for (q_t i = 0; i < q; i++) { - prod += v[i] * (*this)[0][i]; - } - for (q_t i = 0; i < q; i++) { - v_rot[i] = v[i] - 2 * prod * (*this)[0][i]; - } - } else { - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - v_rot[i] += (*this)[i][j] * v[j]; - } - } - } - - return v_rot; - } - - orthogonal_t act(const orthogonal_t & m) const { - orthogonal_t m_rot; - - m_rot.is_reflection = false; - - if (is_reflection) { - for (q_t i = 0; i < q; i++) { - double akOki = 0; - - for (q_t k = 0; k < q; k++) { - akOki += (*this)[0][k] * m[k][i]; - } - - for (q_t j = 0; j < q; j++) { - m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][j]; - } - } - } else { - for (q_t i = 0; i < q; i++) { - m_rot[i].fill(0); - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m_rot[i][j] += (*this)[i][j] * m[j][k]; - } - } - } - } - - return m_rot; - } - - vector_t act_inverse(const vector_t & v) const { - if (is_reflection) { - return this->act(v); // reflections are their own inverse - } else { - vector_t v_rot; - v_rot.fill(0); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - v_rot[i] += (*this)[j][i] * v[j]; - } - } - - return v_rot; - } - } - - vector_t act_inverse(const orthogonal_t & m) const { - if (is_reflection) { - return this->act(m); // reflections are their own inverse - } else { - orthogonal_t m_rot; - m_rot.is_reflection = false; - - for (q_t i = 0; i < q; i++) { - m_rot[i].fill(0); - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m_rot[i][j] += (*this)[j][i] * m[j][k]; - } - } - } - - return m_rot; - } - } - -}; - -template -orthogonal_t generate_rotation_uniform (std::mt19937& r, const system, vector_t>&, v_t) { - std::normal_distribution dist(0.0,1.0); - orthogonal_t ptr; - ptr.is_reflection = true; - - double v2 = 0; - - for (q_t i = 0; i < q; i++) { - ptr[0][i] = dist(r); - v2 += ptr[0][i] * ptr[0][i]; - } - - double mag_v = sqrt(v2); - - for (q_t i = 0; i < q; i++) { - ptr[0][i] /= mag_v; - } - - return ptr; -} - -template -orthogonal_t generate_rotation_perturbation (std::mt19937& r, const system, vector_t>& S, v_t i0, double epsilon, unsigned int n) { - std::normal_distribution dist(0.0,1.0); - orthogonal_t m; - m.is_reflection = true; - - vector_t v; - - if (n > 1) { - std::uniform_int_distribution udist(0, n); - unsigned int rotation = udist(r); - - double cosr = cos(2 * M_PI * rotation / (double)n / 2.0); - double sinr = sin(2 * M_PI * rotation / (double)n / 2.0); - - 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] = S.s[i0][i]; - } - } else { - v = S.s[i0]; - } - - double m_dot_v = 0; - - for (q_t i = 0; i < q; i++) { - m[0][i] = dist(r); // create a random vector - m_dot_v += m[0][i] * v[i]; - } - - double v2 = 0; - - for (q_t i = 0; i < q; i++) { - m[0][i] = m[0][i] - m_dot_v * v[i]; // find the component orthogonal to v - v2 += pow(m[0][i], 2); - } - - double mag_v = sqrt(v2); - - for (q_t i = 0; i < q; i++) { - m[0][i] /= mag_v; // normalize - } - - v2 = 0; - - double factor = epsilon * dist(r); - - for (q_t i = 0; i < q; i++) { - m[0][i] += factor * v[i]; // perturb orthogonal vector in original direction - v2 += pow(m[0][i], 2); - } - - mag_v = sqrt(v2); - - for (q_t i = 0; i < q; i++) { - m[0][i] /= mag_v; // normalize - } - - return m; -} - -} - -#endif - diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp deleted file mode 100644 index a0cb368..0000000 --- a/lib/include/wolff/models/potts.hpp +++ /dev/null @@ -1,68 +0,0 @@ - -#ifndef WOLFF_MODELS_POTTS_H -#define WOLFF_MODELS_POTTS_H - -#include - -#include "vector.hpp" - -namespace wolff { - -#include "../types.h" - -template -class potts_t { - public: - q_t x; - - potts_t() : x(0) {} - potts_t(q_t x) : x(x) {} - - typedef vector_t M_t; - typedef vector_t F_t; - - inline vector_t operator*(v_t a) const { - vector_t result; - result.fill(0); - result[x] = (int)a; - - return result; - } - - inline vector_t operator*(double a) const { - vector_t result; - result.fill(0.0); - result[x] = a; - - return result; - } - - inline vector_t operator-(const potts_t &s) const { - vector_t result; - result.fill(0); - - result[x]++; - result[s.x]--; - - return result; - } - - q_t enumerate() const { - return x; - } -}; - -template -inline typename potts_t::M_t operator*(v_t a, const potts_t& s) { - return s * a; -} - -template -inline typename potts_t::F_t operator*(double a, const potts_t& s) { - return s * a; -} - -} - -#endif - diff --git a/lib/include/wolff/models/symmetric.hpp b/lib/include/wolff/models/symmetric.hpp deleted file mode 100644 index d15a7ba..0000000 --- a/lib/include/wolff/models/symmetric.hpp +++ /dev/null @@ -1,59 +0,0 @@ - -#ifndef WOLFF_MODELS_SYMMETRIC_H -#define WOLFF_MODELS_SYMMETRIC_H - -#include - -#include "potts.hpp" - -namespace wolff { - -#include "../types.h" - -template -class symmetric_t : public std::array { - public: - - symmetric_t() { - for (q_t i = 0; i < q; i++) { - (*this)[i] = i; - } - } - - potts_t act(const potts_t &s) const { - return potts_t((*this)[s.x]); - } - - symmetric_t act(const symmetric_t& r) const { - symmetric_t r_rot; - for (q_t i = 0; i < q; i++) { - r_rot[i] = (*this)[r[i]]; - } - - return r_rot; - } - - potts_t act_inverse(const potts_t& s) const { - for (q_t i = 0; i < q; i++) { - if ((*this)[i] == s.x) { - return potts_t(i); - } - } - - exit(EXIT_FAILURE); - } - - symmetric_t act_inverse(const symmetric_t& r) const { - symmetric_t r_rot; - for (q_t i = 0; i < q; i++) { - r_rot[(*this)[i]] = r[i]; - } - - return r_rot; - } -}; - -} - -#endif - diff --git a/lib/include/wolff/models/vector.hpp b/lib/include/wolff/models/vector.hpp deleted file mode 100644 index e4c4a1c..0000000 --- a/lib/include/wolff/models/vector.hpp +++ /dev/null @@ -1,115 +0,0 @@ - -#ifndef WOLFF_MODELS_VECTOR_H -#define WOLFF_MODELS_VECTOR_H - -#include -#include -#include - -namespace wolff { - -#include - -template -class vector_t : public std::array { - public: - - vector_t() { - this->fill((T)0); - (*this)[0] = (T)1; - } - - vector_t(const T *x) { - for (q_t i = 0; i < q; i++) { - (*this)[i] = x[i]; - } - } - - typedef vector_t M_t; - typedef vector_t F_t; - - template - inline vector_t& operator+=(const vector_t &v) { - for (q_t i = 0; i < q; i++) { - (*this)[i] += (T)v[i]; - } - return *this; - } - - template - inline vector_t& operator-=(const vector_t &v) { - for (q_t i = 0; i < q; i++) { - (*this)[i] -= (T)v[i]; - } - return *this; - } - - inline vector_t operator*(v_t x) const { - vector_t result; - for (q_t i = 0; i < q; i++) { - result[i] = x * (*this)[i]; - } - - return result; - } - - inline vector_t operator*(double x) const { - vector_t result; - for (q_t i = 0; i < q; i++) { - result[i] = x * (*this)[i]; - } - - return result; - } - - inline vector_t operator-(const vector_t& v) const { - vector_t diff = *this; - diff -= v; - return diff; - } - - inline T operator*(const vector_t& v) const { - double prod = 0; - - for (q_t i = 0; i < q; i++) { - prod += v[i] * (*this)[i]; - } - - return prod; - } - - template - inline vector_t operator/(U a) const { - vector_t result; - for (q_t i = 0; i < q; i++) { - result[i] = (*this)[i] / a; - } - - return result; - } -}; - -template -inline vector_t operator*(v_t a, const vector_t&v) { - return v * a; -} - -template -inline vector_t operator*(double a, const vector_t&v) { - return v * a; -} - -template -std::ostream& operator<<(std::ostream& os, const vector_t&v) { - os << "( "; - for (T vi : v) { - os << vi << " "; - } - os << ")"; - return os; -} - -} - -#endif - diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp deleted file mode 100644 index fa9a1b0..0000000 --- a/lib/include/wolff/system.hpp +++ /dev/null @@ -1,105 +0,0 @@ - -#ifndef WOLFF_SYSTEM_H -#define WOLFF_SYSTEM_H - -#include -#include -#include - -#include "graph.hpp" - -namespace wolff { - -#include "types.h" - -template -class measurement; - -template -class system { - public: - v_t nv; // number of vertices - v_t ne; // number of edges - graph G; // the graph defining the lattice with ghost - double T; // the temperature - std::vector s; // the state of the ordinary spins - -#ifdef WOLFF_BOND_DEPENDENCE - std::function Z; // coupling between sites -#else - std::function Z; // coupling between sites -#endif - -#ifndef WOLFF_NO_FIELD - R_t s0; // the current state of the ghost site -#ifdef WOLFF_SITE_DEPENDENCE - std::function B; // coupling with the external field -#else - std::function B; // coupling with the external field -#endif -#endif - -#ifdef WOLFF_USE_FINITE_STATES - std::array, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp; -#ifndef WOLFF_NO_FIELD - std::array, WOLFF_FINITE_STATES_N> Bp; -#endif -#endif - - system(graph g, double T, -#ifdef WOLFF_BOND_DEPENDENCE - std::function Z -#else - std::function Z -#endif -#ifndef WOLFF_NO_FIELD -#ifdef WOLFF_SITE_DEPENDENCE - , std::function B -#else - , std::function B -#endif -#endif - ) : G(g), T(T), Z(Z) -#ifndef WOLFF_NO_FIELD - , s0(), B(B) -#endif - { - nv = G.nv; - ne = G.ne; - s.resize(nv); -#ifndef WOLFF_NO_FIELD - G.add_ghost(); -#endif -#ifdef WOLFF_USE_FINITE_STATES - this->finite_states_init(); -#endif - } - - void flip_cluster(v_t, const R_t&, std::mt19937&, measurement&); - void run_wolff(N_t, std::function &, v_t)> r_gen, measurement& A, std::mt19937& rng); - -#ifdef WOLFF_USE_FINITE_STATES - void finite_states_init() { -#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++) { - Bp[i][j] = 1.0 - exp(-(B(X_t(i)) - B(X_t(j))) / T); - } - } -#endif - for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) { - for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) { - for (q_t k = 0; k < WOLFF_FINITE_STATES_N; k++) { - Zp[i][j][k] = 1.0 - exp(-(Z(X_t(i), X_t(k)) - Z(X_t(j), X_t(k))) / T); - } - } - } - } -#endif - -}; - -} - -#endif - diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h deleted file mode 100644 index 5bbaa6d..0000000 --- a/lib/include/wolff/types.h +++ /dev/null @@ -1,32 +0,0 @@ - -#ifndef WOLFF_TYPES_H -#define WOLFF_TYPES_H - -#include - -typedef uint_fast32_t v_t; // vertex and edge indices -typedef uint_fast8_t q_t; // enumerated states -typedef uint_fast8_t D_t; // dimension -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_D UINT_FAST8_MAX -#define MAX_L UINT_FAST16_MAX -#define MAX_N UINT_FAST64_MAX - -#define PRIv PRIuFAST32 -#define PRIq PRIuFAST8 -#define PRID PRIuFAST8 -#define PRIL PRIuFAST16 -#define PRIN PRIuFAST64 - -#define SCNv SCNuFAST32 -#define SCNq SCNuFAST8 -#define SCND SCNuFAST8 -#define SCNL SCNuFAST16 -#define SCNN SCNuFAST64 - -#endif - diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp deleted file mode 100644 index 4f91be4..0000000 --- a/lib/src/graph.cpp +++ /dev/null @@ -1,55 +0,0 @@ - -#include - -namespace wolff { - -graph::graph() { - D = 0; - L = 0; - nv = 0; - ne = 0; -} - -graph::graph(D_t D, L_t L) : D(D), L(L) { - nv = pow(L, D); - ne = D * nv; - - adjacency.resize(nv); - coordinate.resize(nv); - - for (std::vector adj_i : adjacency) { - 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; - - 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::add_ghost() { - for (std::vector& adj_i : adjacency) { - adj_i.push_back(nv); - } - - adjacency.resize(nv + 1); - coordinate.resize(nv + 1); - adjacency[nv].reserve(nv); - - for (v_t i = 0; i < nv; i++) { - adjacency[nv].push_back(i); - } - - coordinate[nv].resize(coordinate[0].size()); - - ne += nv; - nv++; -} - -} - diff --git a/lib/wolff.hpp b/lib/wolff.hpp new file mode 100644 index 0000000..d21d788 --- /dev/null +++ b/lib/wolff.hpp @@ -0,0 +1,320 @@ + +#ifndef WOLFF_H +#define WOLFF_H + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace wolff{ + + template , class edge_prop = std::tuple<>> + class graph { + /* the graph class describes a lattice on which wolff is run. for most + * purposes, using the default square lattice constructor is sufficient, + * but arbitrary lattices can be constructed + */ + public: + unsigned D; // dimension of space + unsigned L; // linear size + unsigned ne; // number of edges + unsigned nv; // number of vertices + + struct _vertex; + + typedef struct _halfedge { + struct _vertex &self; // reference to the vertex this edge comes from + struct _vertex &neighbor; // reference to the vertex this edge goes to + edge_prop prop; // optional properties + + _halfedge(struct _vertex &v1, struct _vertex &v2) : self(v1), neighbor(v2) {}; + } halfedge; + + typedef struct _vertex { + unsigned ind; // index of the vertex + std::list edges; // list of edges incident on this vertex + vertex_prop prop; // optional properties + } vertex; + + std::vector vertices; /* vector of vertices, length nv, with + * vertices[i].ind = i for all 0 <= i < nv + */ + graph() { + // default constructor for empty graph. use it to build your own! + D = 0; + L = 0; + nv = 0; + ne = 0; + }; + + graph(unsigned D, unsigned L) : D(D), L(L) { + // default constructor for square lattice graph + nv = pow(L, D); + ne = D * nv; + + vertices.resize(nv); + + for (unsigned i = 0; i < nv; i++) { + vertices[i].ind = i; + for (unsigned j = 0; j < D; j++) { + unsigned n1 = pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) + + fmod(i + pow(L, j), pow(L, j + 1)); + unsigned n2 = pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) + + fmod(pow(L, j + 1) + i - pow(L, j), pow(L, j + 1)); + + halfedge f(vertices[i], vertices[n1]); + halfedge b(vertices[i], vertices[n2]); + + vertices[i].edges.push_back(f); + vertices[i].edges.push_back(b); + } + } + }; + + void add_ghost() { + // adds a ghost site to any graph + vertices.resize(nv + 1); + + for (auto it = vertices.begin(); it != std::prev(vertices.end()); it++) { + halfedge e1(*it, vertices[nv]); + it->edges.push_back(e1); + + halfedge e2(vertices[nv], *it); + vertices[nv].edges.push_back(e2); + } + + vertices[nv].ind = nv; + + ne += nv; + nv++; + }; + }; + + template > + class measurement; + + template > + class system { + public: + unsigned nv; // number of vertices + unsigned ne; // number of edges + G_t G; // the graph defining the lattice with ghost + double T; // the temperature + std::vector s; // the state of the ordinary spins + +#ifdef WOLFF_BOND_DEPENDENCE + std::function Z; // coupling between sites +#else + std::function Z; // coupling between sites +#endif + +#ifndef WOLFF_NO_FIELD + R_t s0; // the current state of the ghost site +#ifdef WOLFF_SITE_DEPENDENCE + std::function B; // coupling with the external field +#else + std::function B; // coupling with the external field +#endif +#endif + +#ifdef WOLFF_USE_FINITE_STATES + std::array, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp; +#ifndef WOLFF_NO_FIELD + std::array, WOLFF_FINITE_STATES_N> Bp; +#endif +#endif + + system(G_t g, double T, +#ifdef WOLFF_BOND_DEPENDENCE + std::function Z +#else + std::function Z +#endif +#ifndef WOLFF_NO_FIELD +#ifdef WOLFF_SITE_DEPENDENCE + , std::function B +#else + , std::function B +#endif +#endif + ) : G(g), T(T), Z(Z) +#ifndef WOLFF_NO_FIELD + , s0(), B(B) +#endif + { + nv = G.nv; + ne = G.ne; + s.resize(nv); +#ifndef WOLFF_NO_FIELD + G.add_ghost(); +#endif +#ifdef WOLFF_USE_FINITE_STATES + this->finite_states_init(); +#endif + } + + void flip_cluster(typename G_t::vertex& v0, const R_t& r, std::mt19937& rng, + measurement& A) { + std::uniform_real_distribution dist(0.0, 1.0); + + std::queue queue; + queue.push(v0.ind); + + std::vector visited(G.nv, false); + + while (!queue.empty()) { + unsigned i = queue.front(); + queue.pop(); + + if (!visited[i]) { // don't reprocess anyone we've already visited! + visited[i] = true; + + X_t si_new; +#ifndef WOLFF_NO_FIELD + R_t s0_new; + + bool we_are_ghost = (i == nv); + + if (we_are_ghost) { + s0_new = r.act(s0); + } else +#endif + { + si_new = r.act(s[i]); + } + + for (const typename G_t::halfedge &e : G.vertices[i].edges) { + double dE, p; + unsigned j = e.neighbor.ind; + +#ifndef WOLFF_NO_FIELD + bool neighbor_is_ghost = (j == nv); + + if (we_are_ghost || neighbor_is_ghost) { + X_t s0s_old, s0s_new; + unsigned non_ghost; + + if (neighbor_is_ghost) { + non_ghost = i; + s0s_old = s0.act_inverse(s[i]); + s0s_new = s0.act_inverse(si_new); + } else { + non_ghost = j; + s0s_old = s0.act_inverse(s[j]); + s0s_new = s0_new.act_inverse(s[j]); + } + +#ifdef WOLFF_SITE_DEPENDENCE + dE = B(G.vertices[non_ghost], s0s_old) - B(G.vertices[non_ghost], s0s_new); +#else + dE = B(s0s_old) - B(s0s_new); +#endif + +#ifdef WOLFF_USE_FINITE_STATES + p = Bp[s0s_old.enumerate()][s0s_new.enumerate()]; +#endif + + // run measurement hooks for encountering a ghost bond + A.ghost_bond_visited(*this, G.vertices[non_ghost], s0s_old, s0s_new, dE); + } else // this is a perfectly normal bond! +#endif + { +#ifdef WOLFF_BOND_DEPENDENCE + dE = Z(e, s[i], s[j]) - Z(e, si_new, s[j]); +#else + dE = Z(s[i], s[j]) - Z(si_new, s[j]); +#endif + +#ifdef WOLFF_USE_FINITE_STATES + p = Zp[s[i].enumerate()][si_new.enumerate()][s[j].enumerate()]; +#endif + + // run measurement hooks for encountering a plain bond + A.plain_bond_visited(*this, e, si_new, dE); + } + +#ifndef WOLFF_USE_FINITE_STATES + p = 1.0 - exp(-dE / T); +#endif + + if (dist(rng) < p) { + queue.push(j); // push the neighboring vertex to the queue + } + } + +#ifndef WOLFF_NO_FIELD + if (we_are_ghost) { + A.ghost_site_transformed(*this, s0_new); + s0 = s0_new; + } else +#endif + { + A.plain_site_transformed(*this, G.vertices[i], si_new); + s[i] = si_new; + } + } + } + } + + void run_wolff(unsigned N, + std::function &, const typename G_t::vertex&)> r_gen, measurement& A, std::mt19937& rng) { + std::uniform_int_distribution dist(0, nv - 1); + + for (unsigned n = 0; n < N; n++) { + unsigned i0 = dist(rng); + R_t r = r_gen(rng, *this, G.vertices[i0]); + + A.pre_cluster(n, N, *this, G.vertices[i0], r); + + this->flip_cluster(G.vertices[i0], r, rng, A); + + A.post_cluster(n, N, *this); + } + } + +#ifdef WOLFF_USE_FINITE_STATES + void finite_states_init() { +#ifndef WOLFF_NO_FIELD + for (unsigned i = 0; i < WOLFF_FINITE_STATES_N; i++) { + for (unsigned j = 0; j < WOLFF_FINITE_STATES_N; j++) { + Bp[i][j] = 1.0 - exp(-(B(X_t(i)) - B(X_t(j))) / T); + } + } +#endif + for (unsigned i = 0; i < WOLFF_FINITE_STATES_N; i++) { + for (unsigned j = 0; j < WOLFF_FINITE_STATES_N; j++) { + for (unsigned k = 0; k < WOLFF_FINITE_STATES_N; k++) { + Zp[i][j][k] = 1.0 - exp(-(Z(X_t(i), X_t(k)) - Z(X_t(j), X_t(k))) / T); + } + } + } + } +#endif + + }; + + template + class measurement { + public: + virtual void pre_cluster(unsigned, unsigned, const system&, const typename G_t::vertex& v, const R_t&) {}; + + virtual void plain_bond_visited(const system&, const typename G_t::halfedge& e, const X_t&, double) {}; + virtual void plain_site_transformed(const system&, const typename G_t::vertex& v, const X_t&) {}; + +#ifndef WOLFF_NO_FIELD + virtual void ghost_bond_visited(const system&, const typename G_t::vertex& v, const X_t&, const X_t&, double) {}; + virtual void ghost_site_transformed(const system&, const R_t&) {}; +#endif + + virtual void post_cluster(unsigned, unsigned, const system&) {}; + }; + +} + +#endif + diff --git a/lib/wolff_models/dihedral.hpp b/lib/wolff_models/dihedral.hpp new file mode 100644 index 0000000..aef3be6 --- /dev/null +++ b/lib/wolff_models/dihedral.hpp @@ -0,0 +1,54 @@ + +#ifndef WOLFF_MODELS_DIHEDRAL_H +#define WOLFF_MODELS_DIHEDRAL_H + +#include "potts.hpp" + +namespace wolff { + + template + class dihedral_t { + public: + bool is_reflection; + unsigned x; + + dihedral_t() : is_reflection(false), x(0) {} + dihedral_t(bool x, unsigned y) : is_reflection(x), x(y) {} + + potts_t act(const potts_t& s) const { + if (this->is_reflection) { + return potts_t(((q + this->x) - s.x) % q); + } else { + return potts_t((this->x + s.x) % q); + } + } + + dihedral_t act(dihedral_t r) const { + if (this->is_reflection) { + return dihedral_t(!(r.is_reflection), ((q + this->x) - r.x) % q); + } else { + return dihedral_t(r.is_reflection, (this->x + r.x) % q); + } + } + + potts_t act_inverse(potts_t s) const { + if (this->is_reflection) { + return this->act(s); + } else { + return potts_t(((s.x + q) - this->x) % q); + } + } + + dihedral_t act_inverse(dihedral_t r) const { + if (this->is_reflection) { + return this->act(r); + } else { + return dihedral_t(r.is_reflection, ((r.x + q) - this->x) % q); + } + } + }; + +} + +#endif + diff --git a/lib/wolff_models/dihedral_inf.hpp b/lib/wolff_models/dihedral_inf.hpp new file mode 100644 index 0000000..0b09c59 --- /dev/null +++ b/lib/wolff_models/dihedral_inf.hpp @@ -0,0 +1,55 @@ + +#ifndef WOLFF_MODELS_DIHEDRAL_INF_H +#define WOLFF_MODELS_DIHEDRAL_INF_H + +#include +#include "height.hpp" + +namespace wolff { + + template + class dihedral_inf_t { + public: + bool is_reflection; + T x; + + dihedral_inf_t() : is_reflection(false), x(0) {} + dihedral_inf_t(bool x, T y) : is_reflection(x), x(y) {} + + height_t act(const height_t& h) const { + if (this->is_reflection) { + return height_t(this->x - h.x); + } else { + return height_t(this->x + h.x); + } + } + + dihedral_inf_t act(const dihedral_inf_t& r) const { + if (this->is_reflection) { + return dihedral_inf_t(!r.is_reflection, this->x - r.x); + } else { + return dihedral_inf_t(r.is_reflection, this->x + r.x); + } + } + + height_t act_inverse(const height_t& h) const { + if (this->is_reflection) { + return this->act(h); + } else { + return height_t(h.x - this->x); + } + } + + dihedral_inf_t act_inverse(const dihedral_inf_t& r) const { + if (this->is_reflection) { + return this->act(r); + } else { + return dihedral_inf_t(r.is_reflection, r.x - this->x); + } + } + }; + +} + +#endif + diff --git a/lib/wolff_models/height.hpp b/lib/wolff_models/height.hpp new file mode 100644 index 0000000..ae3ad82 --- /dev/null +++ b/lib/wolff_models/height.hpp @@ -0,0 +1,56 @@ + +#pragma once + +#include + +namespace wolff { + + template + struct height_t { + T x; + + height_t() : x(0) {} + height_t(T x) : x(x) {} + + typedef T M_t; + typedef double F_t; + + inline T operator*(unsigned a) const { + return x * a; + } + + inline double operator*(double a) const { + return x * a; + } + + inline T operator-(const height_t& h) const { + return x - h.x; + } + }; + + template + inline typename height_t::M_t operator*(unsigned a, height_t h) { + return h * a; + } + + template + inline typename height_t::F_t operator*(double a, height_t h) { + return h * a; + } + + template + inline T& operator+=(T& M, const height_t &h) { + M += h.x; + + return M; + } + + template + inline T& operator-=(T& M, const height_t &h) { + M -= h.x; + + return M; + } + +} + diff --git a/lib/wolff_models/ising.hpp b/lib/wolff_models/ising.hpp new file mode 100644 index 0000000..345eeff --- /dev/null +++ b/lib/wolff_models/ising.hpp @@ -0,0 +1,92 @@ + +#ifndef WOLFF_MODELS_ISING_H +#define WOLFF_MODELS_ISING_H + +#define WOLFF_FINITE_STATES_N 2 + +#include "wolff.hpp" + +namespace wolff { + + class ising_t { + public: + bool x; + + ising_t() : x(false) {} + + ising_t(bool x) : x(x) {} + ising_t(unsigned x) : x((bool)x) {} + + ising_t act(const ising_t& s) const { + if (x) { + return ising_t(!s.x); + } else { + return ising_t(s.x); + } + } + + ising_t act_inverse(const ising_t& s) const { + return this->act(s); + } + + typedef int M_t; + typedef double F_t; + + inline M_t operator*(unsigned a) const { + if (x) { + return -(M_t)a; + } else { + return (M_t)a; + } + } + + inline F_t operator*(double a) const { + if (x) { + return -a; + } else { + return a; + } + } + + inline M_t operator-(const ising_t &s) const { + if (x == s.x) { + return 0; + } else { + if (x) { + return -2; + } else { + return 2; + } + } + } + + inline int operator*(const ising_t& s) const { + if (x == s.x) { + return 1; + } else { + return -1; + } + } + + unsigned enumerate() const { + return (unsigned)x; + } + }; + + inline ising_t::M_t operator*(unsigned a, const ising_t& s) { + return s * a; + } + + inline ising_t::F_t operator*(double a, const ising_t& s) { + return s * a; + } + + template + ising_t gen_ising(std::mt19937&, const system&, const typename G_t::vertex&) { + return ising_t(true); + }; + +} + +#endif + diff --git a/lib/wolff_models/orthogonal.hpp b/lib/wolff_models/orthogonal.hpp new file mode 100644 index 0000000..ea9af14 --- /dev/null +++ b/lib/wolff_models/orthogonal.hpp @@ -0,0 +1,206 @@ + +#ifndef WOLFF_MODELS_ORTHOGONAL_H +#define WOLFF_MODELS_ORTHOGONAL_H + +#include +#include + +#include +#include "vector.hpp" + +namespace wolff { + + template + class orthogonal_t : public std::array, q> { + public : + bool is_reflection; + + orthogonal_t() : is_reflection(false) { + for (unsigned i = 0; i < q; i++) { + (*this)[i].fill(0); + (*this)[i][i] = (T)1; + } + } + + vector_t act(const vector_t & v) const { + vector_t v_rot; + v_rot.fill(0); + + if (is_reflection) { + double prod = 0; + for (unsigned i = 0; i < q; i++) { + prod += v[i] * (*this)[0][i]; + } + for (unsigned i = 0; i < q; i++) { + v_rot[i] = v[i] - 2 * prod * (*this)[0][i]; + } + } else { + for (unsigned i = 0; i < q; i++) { + for (unsigned j = 0; j < q; j++) { + v_rot[i] += (*this)[i][j] * v[j]; + } + } + } + + return v_rot; + } + + orthogonal_t act(const orthogonal_t & m) const { + orthogonal_t m_rot; + + m_rot.is_reflection = false; + + if (is_reflection) { + for (unsigned i = 0; i < q; i++) { + double akOki = 0; + + for (unsigned k = 0; k < q; k++) { + akOki += (*this)[0][k] * m[k][i]; + } + + for (unsigned j = 0; j < q; j++) { + m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][j]; + } + } + } else { + for (unsigned i = 0; i < q; i++) { + m_rot[i].fill(0); + for (unsigned j = 0; j < q; j++) { + for (unsigned k = 0; k < q; k++) { + m_rot[i][j] += (*this)[i][j] * m[j][k]; + } + } + } + } + + return m_rot; + } + + vector_t act_inverse(const vector_t & v) const { + if (is_reflection) { + return this->act(v); // reflections are their own inverse + } else { + vector_t v_rot; + v_rot.fill(0); + + for (unsigned i = 0; i < q; i++) { + for (unsigned j = 0; j < q; j++) { + v_rot[i] += (*this)[j][i] * v[j]; + } + } + + return v_rot; + } + } + + vector_t act_inverse(const orthogonal_t & m) const { + if (is_reflection) { + return this->act(m); // reflections are their own inverse + } else { + orthogonal_t m_rot; + m_rot.is_reflection = false; + + for (unsigned i = 0; i < q; i++) { + m_rot[i].fill(0); + for (unsigned j = 0; j < q; j++) { + for (unsigned k = 0; k < q; k++) { + m_rot[i][j] += (*this)[j][i] * m[j][k]; + } + } + } + + return m_rot; + } + } + + }; + + template + orthogonal_t generate_rotation_uniform (std::mt19937& r, const system, vector_t, G_t>&, const typename G_t::vertex&) { + std::normal_distribution dist(0.0,1.0); + orthogonal_t ptr; + ptr.is_reflection = true; + + double v2 = 0; + + for (unsigned i = 0; i < q; i++) { + ptr[0][i] = dist(r); + v2 += ptr[0][i] * ptr[0][i]; + } + + double mag_v = sqrt(v2); + + for (unsigned i = 0; i < q; i++) { + ptr[0][i] /= mag_v; + } + + return ptr; + } + + template + orthogonal_t generate_rotation_perturbation (std::mt19937& r, const system, vector_t, G_t>& S, const typename G_t::vertex& v0, double epsilon, unsigned int n) { + std::normal_distribution dist(0.0,1.0); + orthogonal_t m; + m.is_reflection = true; + + vector_t v; + + if (n > 1) { + std::uniform_int_distribution udist(0, n); + unsigned int rotation = udist(r); + + double cosr = cos(2 * M_PI * rotation / (double)n / 2.0); + double sinr = sin(2 * M_PI * rotation / (double)n / 2.0); + + v[0] = S.s[v0.ind][0] * cosr - S.s[v0.ind][1] * sinr; + v[1] = S.s[v0.ind][1] * cosr + S.s[v0.ind][0] * sinr; + + for (unsigned i = 2; i < q; i++) { + v[i] = S.s[v0.ind][i]; + } + } else { + v = S.s[v0.ind]; + } + + double m_dot_v = 0; + + for (unsigned i = 0; i < q; i++) { + m[0][i] = dist(r); // create a random vector + m_dot_v += m[0][i] * v[i]; + } + + double v2 = 0; + + for (unsigned i = 0; i < q; i++) { + m[0][i] = m[0][i] - m_dot_v * v[i]; // find the component orthogonal to v + v2 += pow(m[0][i], 2); + } + + double mag_v = sqrt(v2); + + for (unsigned i = 0; i < q; i++) { + m[0][i] /= mag_v; // normalize + } + + v2 = 0; + + double factor = epsilon * dist(r); + + for (unsigned i = 0; i < q; i++) { + m[0][i] += factor * v[i]; // perturb orthogonal vector in original direction + v2 += pow(m[0][i], 2); + } + + mag_v = sqrt(v2); + + for (unsigned i = 0; i < q; i++) { + m[0][i] /= mag_v; // normalize + } + + return m; + } + +} + +#endif + diff --git a/lib/wolff_models/potts.hpp b/lib/wolff_models/potts.hpp new file mode 100644 index 0000000..dd56134 --- /dev/null +++ b/lib/wolff_models/potts.hpp @@ -0,0 +1,66 @@ + +#ifndef WOLFF_MODELS_POTTS_H +#define WOLFF_MODELS_POTTS_H + +#include + +#include "vector.hpp" + +namespace wolff { + + template + class potts_t { + public: + unsigned x; + + potts_t() : x(0) {} + potts_t(unsigned x) : x(x) {} + + typedef vector_t M_t; + typedef vector_t F_t; + + inline vector_t operator*(unsigned a) const { + vector_t result; + result.fill(0); + result[x] = (int)a; + + return result; + } + + inline vector_t operator*(double a) const { + vector_t result; + result.fill(0.0); + result[x] = a; + + return result; + } + + inline vector_t operator-(const potts_t &s) const { + vector_t result; + result.fill(0); + + result[x]++; + result[s.x]--; + + return result; + } + + unsigned enumerate() const { + return x; + } + }; + + template + inline typename potts_t::M_t operator*(unsigned a, const potts_t& s) { + return s * a; + } + + template + inline typename potts_t::F_t operator*(double a, const potts_t& s) { + return s * a; + } + +} + +#endif + diff --git a/lib/wolff_models/symmetric.hpp b/lib/wolff_models/symmetric.hpp new file mode 100644 index 0000000..98e53e3 --- /dev/null +++ b/lib/wolff_models/symmetric.hpp @@ -0,0 +1,57 @@ + +#ifndef WOLFF_MODELS_SYMMETRIC_H +#define WOLFF_MODELS_SYMMETRIC_H + +#include + +#include "potts.hpp" + +namespace wolff { + + template + class symmetric_t : public std::array { + public: + + symmetric_t() { + for (unsigned i = 0; i < q; i++) { + (*this)[i] = i; + } + } + + potts_t act(const potts_t &s) const { + return potts_t((*this)[s.x]); + } + + symmetric_t act(const symmetric_t& r) const { + symmetric_t r_rot; + for (unsigned i = 0; i < q; i++) { + r_rot[i] = (*this)[r[i]]; + } + + return r_rot; + } + + potts_t act_inverse(const potts_t& s) const { + for (unsigned i = 0; i < q; i++) { + if ((*this)[i] == s.x) { + return potts_t(i); + } + } + + exit(EXIT_FAILURE); + } + + symmetric_t act_inverse(const symmetric_t& r) const { + symmetric_t r_rot; + for (unsigned i = 0; i < q; i++) { + r_rot[(*this)[i]] = r[i]; + } + + return r_rot; + } + }; + +} + +#endif + diff --git a/lib/wolff_models/vector.hpp b/lib/wolff_models/vector.hpp new file mode 100644 index 0000000..4106e62 --- /dev/null +++ b/lib/wolff_models/vector.hpp @@ -0,0 +1,113 @@ + +#ifndef WOLFF_MODELS_VECTOR_H +#define WOLFF_MODELS_VECTOR_H + +#include +#include +#include + +namespace wolff { + + template + class vector_t : public std::array { + public: + + vector_t() { + this->fill((T)0); + (*this)[0] = (T)1; + } + + vector_t(const T *x) { + for (unsigned i = 0; i < q; i++) { + (*this)[i] = x[i]; + } + } + + typedef vector_t M_t; + typedef vector_t F_t; + + template + inline vector_t& operator+=(const vector_t &v) { + for (unsigned i = 0; i < q; i++) { + (*this)[i] += (T)v[i]; + } + return *this; + } + + template + inline vector_t& operator-=(const vector_t &v) { + for (unsigned i = 0; i < q; i++) { + (*this)[i] -= (T)v[i]; + } + return *this; + } + + inline vector_t operator*(unsigned x) const { + vector_t result; + for (unsigned i = 0; i < q; i++) { + result[i] = x * (*this)[i]; + } + + return result; + } + + inline vector_t operator*(double x) const { + vector_t result; + for (unsigned i = 0; i < q; i++) { + result[i] = x * (*this)[i]; + } + + return result; + } + + inline vector_t operator-(const vector_t& v) const { + vector_t diff = *this; + diff -= v; + return diff; + } + + inline T operator*(const vector_t& v) const { + double prod = 0; + + for (unsigned i = 0; i < q; i++) { + prod += v[i] * (*this)[i]; + } + + return prod; + } + + template + inline vector_t operator/(U a) const { + vector_t result; + for (unsigned i = 0; i < q; i++) { + result[i] = (*this)[i] / a; + } + + return result; + } + }; + + template + inline vector_t operator*(unsigned a, const vector_t&v) { + return v * a; + } + + template + inline vector_t operator*(double a, const vector_t&v) { + return v * a; + } + + template + std::ostream& operator<<(std::ostream& os, const vector_t&v) { + os << "( "; + for (T vi : v) { + os << vi << " "; + } + os << ")"; + return os; + } + +} + +#endif + -- cgit v1.2.3-70-g09d2