diff options
Diffstat (limited to 'lib/include/wolff/models')
-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 |
8 files changed, 0 insertions, 713 deletions
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 - |