summaryrefslogtreecommitdiff
path: root/lib/include
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2019-01-14 15:47:59 -0500
committerJaron Kent-Dobias <jaron@kent-dobias.com>2019-01-14 15:47:59 -0500
commit49ac78a6c04e215950bc9c0f97368605e63da15b (patch)
tree64b770c543a0c90bc7dcbc06ceaaa31e96e541ce /lib/include
parent994cbf1a3b611ff4c94ced3b1630e51fd249d7ed (diff)
downloadc++-49ac78a6c04e215950bc9c0f97368605e63da15b.tar.gz
c++-49ac78a6c04e215950bc9c0f97368605e63da15b.tar.bz2
c++-49ac78a6c04e215950bc9c0f97368605e63da15b.zip
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.
Diffstat (limited to 'lib/include')
-rw-r--r--lib/include/wolff.hpp32
-rw-r--r--lib/include/wolff/cluster.hpp120
-rw-r--r--lib/include/wolff/graph.hpp30
-rw-r--r--lib/include/wolff/measurement.hpp28
-rw-r--r--lib/include/wolff/models/dihedral.hpp56
-rw-r--r--lib/include/wolff/models/dihedral_inf.hpp57
-rw-r--r--lib/include/wolff/models/height.hpp57
-rw-r--r--lib/include/wolff/models/ising.hpp93
-rw-r--r--lib/include/wolff/models/orthogonal.hpp208
-rw-r--r--lib/include/wolff/models/potts.hpp68
-rw-r--r--lib/include/wolff/models/symmetric.hpp59
-rw-r--r--lib/include/wolff/models/vector.hpp115
-rw-r--r--lib/include/wolff/system.hpp105
-rw-r--r--lib/include/wolff/types.h32
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
-