diff options
Diffstat (limited to 'lib/include')
| -rw-r--r-- | lib/include/wolff.hpp | 32 | ||||
| -rw-r--r-- | lib/include/wolff/cluster.hpp | 120 | ||||
| -rw-r--r-- | lib/include/wolff/graph.hpp | 30 | ||||
| -rw-r--r-- | lib/include/wolff/measurement.hpp | 28 | ||||
| -rw-r--r-- | lib/include/wolff/models/dihedral.hpp | 56 | ||||
| -rw-r--r-- | lib/include/wolff/models/dihedral_inf.hpp | 57 | ||||
| -rw-r--r-- | lib/include/wolff/models/height.hpp | 57 | ||||
| -rw-r--r-- | lib/include/wolff/models/ising.hpp | 93 | ||||
| -rw-r--r-- | lib/include/wolff/models/orthogonal.hpp | 208 | ||||
| -rw-r--r-- | lib/include/wolff/models/potts.hpp | 68 | ||||
| -rw-r--r-- | lib/include/wolff/models/symmetric.hpp | 59 | ||||
| -rw-r--r-- | lib/include/wolff/models/vector.hpp | 115 | ||||
| -rw-r--r-- | lib/include/wolff/system.hpp | 105 | ||||
| -rw-r--r-- | lib/include/wolff/types.h | 32 |
14 files changed, 0 insertions, 1060 deletions
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 <class R_t, class X_t> -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, 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 <random> -#include <cmath> -#include <vector> -#include <queue> - -#include "system.hpp" - -namespace wolff { - -template <class R_t, class X_t> -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(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 <cmath> -#include <vector> - -namespace wolff { - -#include "types.h" - -class graph { - public: - D_t D; - L_t L; - v_t ne; - v_t nv; - std::vector<std::vector<v_t>> adjacency; - std::vector<std::vector<double>> 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 R_t, class X_t> -class measurement { - public: - virtual void pre_cluster(N_t, N_t, const system<R_t, X_t>&, v_t, const R_t&) {}; - - virtual void plain_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, v_t, double) {}; - virtual void plain_site_transformed(const system<R_t, X_t>&, v_t, const X_t&) {}; - -#ifndef WOLFF_NO_FIELD - virtual void ghost_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) {}; - virtual void ghost_site_transformed(const system<R_t, X_t>&, const R_t&) {}; -#endif - - virtual void post_cluster(N_t, N_t, const system<R_t, X_t>&) {}; -}; - -} - -#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 <q_t q> -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<q> act(const potts_t<q>& s) const { - if (this->is_reflection) { - return potts_t<q>(((q + this->x) - s.x) % q); - } else { - return potts_t<q>((this->x + s.x) % q); - } - } - - dihedral_t<q> act(dihedral_t<q> r) const { - if (this->is_reflection) { - return dihedral_t<q>(!(r.is_reflection), ((q + this->x) - r.x) % q); - } else { - return dihedral_t<q>(r.is_reflection, (this->x + r.x) % q); - } - } - - potts_t<q> act_inverse(potts_t<q> s) const { - if (this->is_reflection) { - return this->act(s); - } else { - return potts_t<q>(((s.x + q) - this->x) % q); - } - } - - dihedral_t<q> act_inverse(dihedral_t<q> r) const { - if (this->is_reflection) { - return this->act(r); - } else { - return dihedral_t<q>(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 <cmath> -#include "height.hpp" - -namespace wolff { - -#include "../types.h" - -template <class T> -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<T> act(const height_t<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<T> act(const dihedral_inf_t<T>& r) const { - if (this->is_reflection) { - return dihedral_inf_t<T>(!r.is_reflection, this->x - r.x); - } else { - return dihedral_inf_t<T>(r.is_reflection, this->x + r.x); - } - } - - height_t<T> act_inverse(const height_t<T>& h) const { - if (this->is_reflection) { - return this->act(h); - } else { - return height_t(h.x - this->x); - } - } - - dihedral_inf_t<T> act_inverse(const dihedral_inf_t<T>& r) const { - if (this->is_reflection) { - return this->act(r); - } else { - return dihedral_inf_t<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 <cmath> -#include <wolff/types.h> - -namespace wolff { - -template <class T> -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 <class T> -inline typename height_t<T>::M_t operator*(v_t a, height_t<T> h) { - return h * a; -} - -template <class T> -inline typename height_t<T>::F_t operator*(double a, height_t<T> h) { - return h * a; -} - -template <class T> -inline T& operator+=(T& M, const height_t<T> &h) { - M += h.x; - - return M; -} - -template <class T> -inline T& operator-=(T& M, const height_t<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<ising_t, ising_t>&, 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 <random> -#include <cmath> - -#include "../system.hpp" -#include "vector.hpp" - -namespace wolff { - -#include "../types.h" - -template <q_t q, class T> -class orthogonal_t : public std::array<std::array<T, q>, 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<q, T> act(const vector_t <q, T>& v) const { - vector_t <q, 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<q, T> act(const orthogonal_t <q, T>& m) const { - orthogonal_t <q, 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 <q, T> act_inverse(const vector_t <q, T>& v) const { - if (is_reflection) { - return this->act(v); // reflections are their own inverse - } else { - vector_t <q, 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 <q, T> act_inverse(const orthogonal_t <q, T>& m) const { - if (is_reflection) { - return this->act(m); // reflections are their own inverse - } else { - orthogonal_t <q, 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 <q_t q> -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; - - 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 <q_t q> -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; - - vector_t <q, double> v; - - if (n > 1) { - std::uniform_int_distribution<unsigned int> 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 <cmath> - -#include "vector.hpp" - -namespace wolff { - -#include "../types.h" - -template <q_t q> -class potts_t { - public: - q_t x; - - potts_t() : x(0) {} - potts_t(q_t x) : x(x) {} - - typedef vector_t<q, int> M_t; - typedef vector_t<q, double> F_t; - - inline vector_t<q, int> operator*(v_t a) const { - vector_t<q, int> result; - result.fill(0); - result[x] = (int)a; - - return result; - } - - inline vector_t<q, double> operator*(double a) const { - vector_t<q, double> result; - result.fill(0.0); - result[x] = a; - - return result; - } - - inline vector_t<q, int> operator-(const potts_t<q> &s) const { - vector_t<q, int> result; - result.fill(0); - - result[x]++; - result[s.x]--; - - return result; - } - - q_t enumerate() const { - return x; - } -}; - -template<q_t q> -inline typename potts_t<q>::M_t operator*(v_t a, const potts_t<q>& s) { - return s * a; -} - -template<q_t q> -inline typename potts_t<q>::F_t operator*(double a, const potts_t<q>& 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 <array> - -#include "potts.hpp" - -namespace wolff { - -#include "../types.h" - -template <q_t q> -class symmetric_t : public std::array<q_t, q> { - public: - - symmetric_t() { - for (q_t i = 0; i < q; i++) { - (*this)[i] = i; - } - } - - potts_t<q> act(const potts_t<q> &s) const { - return potts_t<q>((*this)[s.x]); - } - - symmetric_t<q> act(const symmetric_t<q>& r) const { - symmetric_t<q> r_rot; - for (q_t i = 0; i < q; i++) { - r_rot[i] = (*this)[r[i]]; - } - - return r_rot; - } - - potts_t<q> act_inverse(const potts_t<q>& s) const { - for (q_t i = 0; i < q; i++) { - if ((*this)[i] == s.x) { - return potts_t<q>(i); - } - } - - exit(EXIT_FAILURE); - } - - symmetric_t<q> act_inverse(const symmetric_t<q>& r) const { - symmetric_t<q> 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 <cmath> -#include <array> -#include <iostream> - -namespace wolff { - -#include <wolff/types.h> - -template <q_t q, class T> -class vector_t : public std::array<T, q> { - 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 <q, T> M_t; - typedef vector_t <q, double> F_t; - - template <class U> - inline vector_t<q, T>& operator+=(const vector_t<q, U> &v) { - for (q_t i = 0; i < q; i++) { - (*this)[i] += (T)v[i]; - } - return *this; - } - - template <class U> - inline vector_t<q, T>& operator-=(const vector_t<q, U> &v) { - for (q_t i = 0; i < q; i++) { - (*this)[i] -= (T)v[i]; - } - return *this; - } - - inline vector_t<q, T> operator*(v_t x) const { - vector_t<q, T> result; - for (q_t i = 0; i < q; i++) { - result[i] = x * (*this)[i]; - } - - return result; - } - - inline vector_t<q, double> operator*(double x) const { - vector_t<q, double> result; - for (q_t i = 0; i < q; i++) { - result[i] = x * (*this)[i]; - } - - return result; - } - - inline vector_t<q, T> operator-(const vector_t<q, T>& v) const { - vector_t<q, T> diff = *this; - diff -= v; - return diff; - } - - inline T operator*(const vector_t<q, T>& v) const { - double prod = 0; - - for (q_t i = 0; i < q; i++) { - prod += v[i] * (*this)[i]; - } - - return prod; - } - - template <class U> - inline vector_t<q, T> operator/(U a) const { - vector_t<q, T> result; - for (q_t i = 0; i < q; i++) { - result[i] = (*this)[i] / a; - } - - return result; - } -}; - -template<q_t q, class T> -inline vector_t<q, T> operator*(v_t a, const vector_t<q, T>&v) { - return v * a; -} - -template<q_t q, class T> -inline vector_t<q, double> operator*(double a, const vector_t<q, T>&v) { - return v * a; -} - -template<q_t q, class T> -std::ostream& operator<<(std::ostream& os, const vector_t<q, 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 <functional> -#include <vector> -#include <random> - -#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 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<X_t> s; // the state of the ordinary spins - -#ifdef WOLFF_BOND_DEPENDENCE - std::function <double(v_t, const X_t&, v_t, const X_t&)> Z; // coupling between sites -#else - std::function <double(const X_t&, const X_t&)> 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 <double(v_t, const X_t&)> B; // coupling with the external field -#else - std::function <double(const X_t&)> B; // coupling with the external field -#endif -#endif - -#ifdef WOLFF_USE_FINITE_STATES - std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp; -#ifndef WOLFF_NO_FIELD - std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Bp; -#endif -#endif - - system(graph g, double T, -#ifdef WOLFF_BOND_DEPENDENCE - std::function <double(v_t, const X_t&, v_t, const X_t&)> Z -#else - std::function <double(const X_t&, const X_t&)> Z -#endif -#ifndef WOLFF_NO_FIELD -#ifdef WOLFF_SITE_DEPENDENCE - , std::function <double(v_t, const X_t&)> B -#else - , std::function <double(const X_t&)> 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<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); - -#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 <inttypes.h> - -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 - |