summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/include/wolff.hpp33
-rw-r--r--lib/include/wolff/cluster.hpp118
-rw-r--r--lib/include/wolff/finite_states.hpp41
-rw-r--r--lib/include/wolff/graph.hpp7
-rw-r--r--lib/include/wolff/meas.h19
-rw-r--r--lib/include/wolff/measurement.hpp24
-rw-r--r--lib/include/wolff/models/dihedral.hpp48
-rw-r--r--lib/include/wolff/models/dihedral_inf.hpp47
-rw-r--r--lib/include/wolff/models/height.hpp43
-rw-r--r--lib/include/wolff/models/ising.hpp81
-rw-r--r--lib/include/wolff/models/orthogonal.hpp199
-rw-r--r--lib/include/wolff/models/potts.hpp50
-rw-r--r--lib/include/wolff/models/symmetric.hpp51
-rw-r--r--lib/include/wolff/models/vector.hpp108
-rw-r--r--lib/include/wolff/state.hpp68
-rw-r--r--lib/include/wolff/system.hpp70
-rw-r--r--lib/include/wolff/types.h25
-rw-r--r--lib/src/graph.cpp36
18 files changed, 845 insertions, 223 deletions
diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp
index b730c8d..b78cd4b 100644
--- a/lib/include/wolff.hpp
+++ b/lib/include/wolff.hpp
@@ -1,30 +1,27 @@
+#ifndef WOLFF_H
+#define WOLFF_H
+
#include "wolff/cluster.hpp"
-#include "wolff/state.hpp"
template <class R_t, class X_t>
-void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(std::mt19937&, X_t)> gen_R, wolff_measurement<R_t, X_t>& m, std::mt19937& r) {
-
-#ifdef FINITE_STATES
-#ifdef NOFIELD
- initialize_probs(s.J, s.T);
-#else
- initialize_probs(s.J, s.H, s.T);
-#endif
-#endif
+void wolff(N_t N, wolff_system<R_t, X_t>& S,
+ std::function <R_t(std::mt19937&, X_t)> r_gen,
+ wolff_measurement<R_t, X_t>& A, std::mt19937& rng) {
- std::uniform_int_distribution<v_t> dist(0, s.nv);
+ std::uniform_int_distribution<v_t> dist(0, S.nv - 1);
- for (count_t steps = 0; steps < N; steps++) {
- v_t v0 = dist(r);
- R_t step = gen_R(r, s.spins[v0]);
+ for (N_t n = 0; n < N; n++) {
+ v_t i0 = dist(rng);
+ R_t r = r_gen(rng, S.s[i0]);
- m.pre_cluster(s, steps, N, v0, step);
+ A.pre_cluster(n, N, S, i0, r);
- flip_cluster<R_t, X_t>(s, v0, step, r, m);
+ wolff_cluster_flip<R_t, X_t>(S, i0, r, rng, A);
- m.post_cluster(s, steps, N);
+ A.post_cluster(n, N, S);
}
-
}
+#endif
+
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
index 805e2c3..055cdf3 100644
--- a/lib/include/wolff/cluster.hpp
+++ b/lib/include/wolff/cluster.hpp
@@ -1,120 +1,122 @@
-#pragma once
+#ifndef WOLFF_CLUSTER_H
+#define WOLFF_CLUSTER_H
#include <random>
#include <cmath>
#include <vector>
-#include <stack>
+#include <queue>
#include "types.h"
-#include "state.hpp"
+#include "system.hpp"
#include "graph.hpp"
-#include "meas.h"
+#include "measurement.hpp"
template <class R_t, class X_t>
-void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, std::mt19937& rand, wolff_measurement<R_t, X_t>& m) {
+void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
+ std::mt19937& rng, wolff_measurement<R_t, X_t>& A) {
std::uniform_real_distribution<double> dist(0.0, 1.0);
- std::stack<v_t> stack;
- stack.push(v0);
+ std::queue<v_t> queue;
+ queue.push(i0);
- std::vector<bool> marks(s.g.nv, false);
+ std::vector<bool> marks(S.G.nv, false);
- while (!stack.empty()) {
- v_t v = stack.top();
- stack.pop();
+ while (!queue.empty()) {
+ v_t i = queue.front();
+ queue.pop();
- if (!marks[v]) { // don't reprocess anyone we've already visited!
- marks[v] = true;
+ if (!marks[i]) { // don't reprocess anyone we've already visited!
+ marks[i] = true;
X_t si_new;
-#ifndef NOFIELD
- R_t R_new;
+#ifndef WOLFF_NO_FIELD
+ R_t s0_new;
- bool v_is_ghost = (v == s.nv); // ghost site has the last index
+ bool we_are_ghost = (i == S.nv);
- if (v_is_ghost) {
- R_new = r.act(s.R); // if we are, then we're moving the transformation
+ if (we_are_ghost) {
+ s0_new = r.act(S.s0);
} else
#endif
{
- si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site
+ si_new = r.act(S.s[i]);
}
- for (const v_t &vn : s.g.v_adj[v]) {
- double dE, prob;
+ for (const v_t &j : S.G.adj[i]) {
+ double dE, p;
-#ifndef NOFIELD
- bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
+#ifndef WOLFF_NO_FIELD
+ bool neighbor_is_ghost = (j == S.nv);
- if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond
- X_t rs_old, rs_new;
+ if (we_are_ghost || neighbor_is_ghost) {
+ X_t s0s_old, s0s_new;
v_t non_ghost;
- if (vn_is_ghost) {
- // if our neighbor is the ghost, the current site is a normal
- // spin - rotate it back!
- rs_old = s.R.act_inverse(s.spins[v]);
- rs_new = s.R.act_inverse(si_new);
- non_ghost = v;
+ if (neighbor_is_ghost) {
+ non_ghost = i;
+ s0s_old = S.s0.act_inverse(S.s[i]);
+ s0s_new = S.s0.act_inverse(si_new);
} else {
- /* if we're the ghost, we need to rotate our neighbor back in
- both the old and new ways */
- rs_old = s.R.act_inverse(s.spins[vn]);
- rs_new = R_new.act_inverse(s.spins[vn]);
- non_ghost = vn;
+ non_ghost = j;
+ s0s_old = S.s0.act_inverse(S.s[j]);
+ s0s_new = s0_new.act_inverse(S.s[j]);
}
-#ifdef SITE_DEPENDENCE
- dE = s.H(non_ghost, rs_old) - s.H(non_ghost, rs_new);
+#ifdef WOLFF_SITE_DEPENDENCE
+ dE = S.B(non_ghost, s0s_old) - S.B(non_ghost, s0s_new);
#else
- dE = s.H(rs_old) - s.H(rs_new);
+ dE = S.B(s0s_old) - S.B(s0s_new);
#endif
-#ifdef FINITE_STATES
- prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
+#ifdef WOLFF_FINITE_STATES
+ p = finite_states_Bp[finite_states_enum(s0s_old)]
+ [finite_states_enum(s0s_new)];
#endif
// run measurement hooks for encountering a ghost bond
- m.ghost_bond_added(non_ghost, rs_old, rs_new, dE);
+ A.ghost_bond_visited(S, non_ghost, s0s_old, s0s_new, dE);
} else // this is a perfectly normal bond!
#endif
{
-#ifdef BOND_DEPENDENCE
- dE = s.J(v, s.spins[v], vn, s.spins[vn]) - s.J(v, si_new, vn, s.spins[vn]);
+#ifdef WOLFF_BOND_DEPENDENCE
+ dE = S.Z(i, S.s[i], j, S.s[j]) - S.Z(i, si_new, j, S.s[j]);
#else
- dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]);
+ dE = S.Z(S.s[i], S.s[j]) - S.Z(si_new, S.s[j]);
#endif
-
-#ifdef FINITE_STATES
- prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])];
+#ifdef WOLFF_FINITE_STATES
+ p = finite_states_Zp[finite_states_enum(S.s[i])]
+ [finite_states_enum(si_new)]
+ [finite_states_enum(S.s[j])];
#endif
// run measurement hooks for encountering a plain bond
- m.plain_bond_added(v, s.spins[v], si_new, vn, s.spins[vn], dE);
+ A.plain_bond_visited(S, i, si_new, j, dE);
}
#ifndef FINITE_STATES
- prob = 1.0 - exp(-dE / s.T);
+ p = 1.0 - exp(-dE / S.T);
#endif
- if (dist(rand) < prob) {
- stack.push(vn); // push the neighboring vertex to the stack
+ if (dist(rng) < p) {
+ queue.push(j); // push the neighboring vertex to the queue
}
}
-#ifndef NOFIELD
- if (v_is_ghost) {
- m.ghost_site_transformed(s.R, R_new);
- s.R = R_new;
+#ifndef WOLFF_NO_FIELD
+ if (we_are_ghost) {
+ A.ghost_site_transformed(S, s0_new);
+ S.s0 = s0_new;
} else
#endif
{
- m.plain_site_transformed(v, s.spins[v], si_new);
- s.spins[v] = si_new;
+ A.plain_site_transformed(S, i, si_new);
+ S.s[i] = si_new;
}
}
}
}
+#endif
+
diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp
index 426edad..4cca69e 100644
--- a/lib/include/wolff/finite_states.hpp
+++ b/lib/include/wolff/finite_states.hpp
@@ -1,40 +1,35 @@
-#pragma once
+
+#ifndef WOLFF_FINITE_STATES
+#define WOLFF_FINITE_STATES
#include <cmath>
-#include <functional>
#include <array>
-#define FINITE_STATES
-
-// must have N_STATES, states[N_STATES], and state_to_ind defined before
-// invoking header
+#include "system.hpp"
-std::array<std::array<std::array<double, N_STATES>, N_STATES>, N_STATES> J_probs;
-#ifndef NOFIELD
-std::array<std::array<double, N_STATES>, N_STATES> H_probs;
+std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Zp;
+#ifndef WOLFF_NO_FIELD
+std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Bp;
#endif
-template <class X_t>
-#ifndef NOFIELD
-void initialize_probs(std::function <double(X_t, X_t)> J, std::function <double(X_t)> H, double T) {
- for (q_t i = 0; i < N_STATES; i++) {
- for (q_t j = 0; j < N_STATES; j++) {
- H_probs[i][j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T);
+template <class R_t, class X_t>
+void finite_states_init(const wolff_system<R_t, X_t>& S) {
+#ifndef WOLFF_NO_FIELD
+ for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+ for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+ finite_states_Bp[i][j] = 1.0 - exp(-(S.B(finite_states_possible[i]) - S.B(finite_states_possible[j])) / S.T);
}
}
-#else
-void initialize_probs(std::function <double(X_t, X_t)> J, double T) {
#endif
- for (q_t i = 0; i < N_STATES; i++) {
- for (q_t j = 0; j < N_STATES; j++) {
- for (q_t k = 0; k < N_STATES; k++) {
- J_probs[i][j][k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T);
+ 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++) {
+ finite_states_Zp[i][j][k] = 1.0 - exp(-(S.Z(finite_states_possible[i], finite_states_possible[k]) - S.Z(finite_states_possible[j], finite_states_possible[k])) / S.T);
}
}
}
}
-
-
+#endif
diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp
index 0aeb6af..fcd73fa 100644
--- a/lib/include/wolff/graph.hpp
+++ b/lib/include/wolff/graph.hpp
@@ -1,5 +1,6 @@
-#pragma once
+#ifndef WOLFF_GRAPH_H
+#define WOLFF_GRAPH_H
#include <cmath>
#include <vector>
@@ -15,10 +16,12 @@ class graph_t {
public:
v_t ne;
v_t nv;
- std::vector<std::vector<v_t>> v_adj;
+ std::vector<std::vector<v_t>> adj;
std::vector<std::vector<double>> coordinate;
graph_t(D_t D, L_t L, lattice_t lat = SQUARE_LATTICE);
void add_ext();
};
+#endif
+
diff --git a/lib/include/wolff/meas.h b/lib/include/wolff/meas.h
deleted file mode 100644
index 4895eac..0000000
--- a/lib/include/wolff/meas.h
+++ /dev/null
@@ -1,19 +0,0 @@
-
-#pragma once
-
-#include "types.h"
-
-template <class R_t, class X_t>
-class wolff_measurement {
- public:
- virtual void pre_cluster(const state_t<R_t, X_t>&, count_t, count_t, v_t, const R_t&) = 0;
-
- virtual void plain_bond_added(v_t, const X_t&, const X_t&, v_t, const X_t&, double) = 0;
- virtual void ghost_bond_added(v_t, const X_t&, const X_t&, double) = 0;
-
- virtual void plain_site_transformed(v_t, const X_t&, const X_t&) = 0;
- virtual void ghost_site_transformed(const R_t&, const R_t&) = 0;
-
- virtual void post_cluster(const state_t<R_t, X_t>&, count_t, count_t) = 0;
-};
-
diff --git a/lib/include/wolff/measurement.hpp b/lib/include/wolff/measurement.hpp
new file mode 100644
index 0000000..a6a5f79
--- /dev/null
+++ b/lib/include/wolff/measurement.hpp
@@ -0,0 +1,24 @@
+
+#ifndef WOLFF_MEASUREMENTS
+#define WOLFF_MEASUREMENTS
+
+#include "system.hpp"
+
+template <class R_t, class X_t>
+class wolff_measurement {
+ public:
+ virtual void pre_cluster(N_t, N_t, const wolff_system<R_t, X_t>&, v_t, const R_t&) = 0;
+
+ virtual void plain_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, v_t, double) = 0;
+ virtual void plain_site_transformed(const wolff_system<R_t, X_t>&, v_t, const X_t&) = 0;
+
+#ifndef WOLFF_NO_FIELD
+ virtual void ghost_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) = 0;
+ virtual void ghost_site_transformed(const wolff_system<R_t, X_t>&, const R_t&) = 0;
+#endif
+
+ virtual void post_cluster(N_t, N_t, const wolff_system<R_t, X_t>&) = 0;
+};
+
+#endif
+
diff --git a/lib/include/wolff/models/dihedral.hpp b/lib/include/wolff/models/dihedral.hpp
new file mode 100644
index 0000000..cbc5687
--- /dev/null
+++ b/lib/include/wolff/models/dihedral.hpp
@@ -0,0 +1,48 @@
+
+#pragma once
+
+#include <wolff/types.h>
+#include "potts.hpp"
+
+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);
+ }
+ }
+};
+
diff --git a/lib/include/wolff/models/dihedral_inf.hpp b/lib/include/wolff/models/dihedral_inf.hpp
new file mode 100644
index 0000000..bded387
--- /dev/null
+++ b/lib/include/wolff/models/dihedral_inf.hpp
@@ -0,0 +1,47 @@
+
+#include "../types.h"
+#include <cmath>
+#include "height.hpp"
+
+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);
+ }
+ }
+};
+
diff --git a/lib/include/wolff/models/height.hpp b/lib/include/wolff/models/height.hpp
new file mode 100644
index 0000000..bd0ceb6
--- /dev/null
+++ b/lib/include/wolff/models/height.hpp
@@ -0,0 +1,43 @@
+
+#pragma once
+
+#include <cmath>
+#include <wolff/types.h>
+
+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 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
new file mode 100644
index 0000000..c6a3a47
--- /dev/null
+++ b/lib/include/wolff/models/ising.hpp
@@ -0,0 +1,81 @@
+
+#ifndef WOLFF_MODELS_ISING
+#define WOLFF_MODELS_ISING
+
+#include "../types.h"
+
+class ising_t {
+ public:
+ bool x;
+
+ ising_t() : x(false) {}
+
+ ising_t(bool x) : x(x) {}
+ ising_t(int 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;
+ }
+ }
+};
+
+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;
+}
+
+#define WOLFF_FINITE_STATES_N 2
+const ising_t finite_states_possible[2] = {ising_t(0), ising_t(1)};
+q_t finite_states_enum(ising_t state) { return (q_t)state.x; }
+
+#endif
+
diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp
new file mode 100644
index 0000000..9dd5ddd
--- /dev/null
+++ b/lib/include/wolff/models/orthogonal.hpp
@@ -0,0 +1,199 @@
+
+#pragma once
+
+#include <random>
+#include <cmath>
+
+#include <wolff/types.h>
+#include "vector.hpp"
+
+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 vector_t <q, double>& v) {
+ 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 vector_t <q, double>& v0, 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] = v0[0] * cosr - v0[1] * sinr;
+ v[1] = v0[1] * cosr + v0[0] * sinr;
+
+ for (q_t i = 2; i < q; i++) {
+ v[i] = v0[i];
+ }
+ } else {
+ v = v0;
+ }
+
+ 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;
+}
+
diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp
new file mode 100644
index 0000000..903f25f
--- /dev/null
+++ b/lib/include/wolff/models/potts.hpp
@@ -0,0 +1,50 @@
+
+#pragma once
+
+#include <cmath>
+
+#include "../types.h"
+#include "vector.hpp"
+
+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;
+ }
+};
+
+const potts_t<POTTSQ> states[256] = {{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {21}, {22}, {23}, {24}, {25}, {26}, {27}, {28}, {29}, {30}, {31}, {32}, {33}, {34}, {35}, {36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46}, {47}, {48}, {49}, {50}, {51}, {52}, {53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}, {61}, {62}, {63}, {64}, {65}, {66}, {67}, {68}, {69}, {70}, {71}, {72}, {73}, {74}, {75}, {76}, {77}, {78}, {79}, {80}, {81}, {82}, {83}, {84}, {85}, {86}, {87}, {88}, {89}, {90}, {91}, {92}, {93}, {94}, {95}, {96}, {97}, {98}, {99}, {100}, {101}, {102}, {103}, {104}, {105}, {106}, {107}, {108}, {109}, {110}, {111}, {112}, {113}, {114}, {115}, {116}, {117}, {118}, {119}, {120}, {121}, {122}, {123}, {124}, {125}, {126}, {127}, {128}, {129}, {130}, {131}, {132}, {133}, {134}, {135}, {136}, {137}, {138}, {139}, {140}, {141}, {142}, {143}, {144}, {145}, {146}, {147}, {148}, {149}, {150}, {151}, {152}, {153}, {154}, {155}, {156}, {157}, {158}, {159}, {160}, {161}, {162}, {163}, {164}, {165}, {166}, {167}, {168}, {169}, {170}, {171}, {172}, {173}, {174}, {175}, {176}, {177}, {178}, {179}, {180}, {181}, {182}, {183}, {184}, {185}, {186}, {187}, {188}, {189}, {190}, {191}, {192}, {193}, {194}, {195}, {196}, {197}, {198}, {199}, {200}, {201}, {202}, {203}, {204}, {205}, {206}, {207}, {208}, {209}, {210}, {211}, {212}, {213}, {214}, {215}, {216}, {217}, {218}, {219}, {220}, {221}, {222}, {223}, {224}, {225}, {226}, {227}, {228}, {229}, {230}, {231}, {232}, {233}, {234}, {235}, {236}, {237}, {238}, {239}, {240}, {241}, {242}, {243}, {244}, {245}, {246}, {247}, {248}, {249}, {250}, {251}, {252}, {253}, {254}, {255}};
+template <q_t q>
+q_t finite_states_enum(potts_t<q> state) { return (q_t)state.x; }
+
diff --git a/lib/include/wolff/models/symmetric.hpp b/lib/include/wolff/models/symmetric.hpp
new file mode 100644
index 0000000..7bcc9a6
--- /dev/null
+++ b/lib/include/wolff/models/symmetric.hpp
@@ -0,0 +1,51 @@
+
+#pragma once
+
+#include <array>
+
+#include "../types.h"
+#include "potts.hpp"
+
+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;
+ }
+};
+
diff --git a/lib/include/wolff/models/vector.hpp b/lib/include/wolff/models/vector.hpp
new file mode 100644
index 0000000..1d635c8
--- /dev/null
+++ b/lib/include/wolff/models/vector.hpp
@@ -0,0 +1,108 @@
+
+#pragma once
+
+#include <cmath>
+#include <array>
+#include <iostream>
+
+#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;
+}
+
diff --git a/lib/include/wolff/state.hpp b/lib/include/wolff/state.hpp
deleted file mode 100644
index 4909881..0000000
--- a/lib/include/wolff/state.hpp
+++ /dev/null
@@ -1,68 +0,0 @@
-
-#pragma once
-
-#include <functional>
-#include <vector>
-
-#include "types.h"
-#include "graph.hpp"
-
-template <class R_t, class X_t>
-class state_t {
- public:
- D_t D; // the dimension of the system
- L_t L; // the linear size of the lattice
- v_t nv; // the number of vertices in the original lattice
- v_t ne; // the number of edges in the original lattice
- graph_t g; // the graph defining the lattice with ghost
- double T; // the temperature
- std::vector<X_t> spins; // the state of the ordinary spins
-#ifndef NOFIELD
- R_t R; // the current state of the ghost site
-#endif
-
-#ifdef BOND_DEPENDENCE
- std::function <double(v_t, const X_t&, v_t, const X_t&)> J; // coupling between sites
-#else
- std::function <double(const X_t&, const X_t&)> J; // coupling between sites
-#endif
-
-#ifndef NOFIELD
-#ifdef SITE_DEPENDENCE
- std::function <double(v_t, const X_t&)> H; // coupling with the external field
-#else
- std::function <double(const X_t&)> H; // coupling with the external field
-#endif
-#endif
-
- state_t(D_t D, L_t L, double T,
-#ifdef BOND_DEPENDENCE
- std::function <double(v_t, const X_t&, v_t, const X_t&)> J
-#else
- std::function <double(const X_t&, const X_t&)> J
-#endif
-#ifndef NOFIELD
-#ifdef SITE_DEPENDENCE
- , std::function <double(v_t, const X_t&)> H
-#else
- , std::function <double(const X_t&)> H
-#endif
-#endif
- , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), g(D, L, lat), T(T),
-#ifndef NOFIELD
- R(),
-#endif
- J(J)
-#ifndef NOFIELD
- , H(H)
-#endif
- {
- nv = g.nv;
- ne = g.ne;
- spins.resize(nv);
-#ifndef NOFIELD
- g.add_ext();
-#endif
- }
-};
-
diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp
new file mode 100644
index 0000000..ac82f44
--- /dev/null
+++ b/lib/include/wolff/system.hpp
@@ -0,0 +1,70 @@
+
+#ifndef WOLFF_STATE_H
+#define WOLFF_STATE_H
+
+#include <functional>
+#include <vector>
+
+#include "types.h"
+#include "graph.hpp"
+
+template <class R_t, class X_t>
+class wolff_system {
+ public:
+ D_t D; // dimension
+ L_t L; // linear size
+ v_t nv; // number of vertices
+ v_t ne; // number of edges
+ graph_t G; // the graph defining the lattice with ghost
+ double T; // the temperature
+ std::vector<X_t> s; // the state of the ordinary spins
+#ifndef WOLFF_NO_FIELD
+ R_t s0; // the current state of the ghost site
+#endif
+
+#ifdef WOLFF_BOND_DEPENDENCE
+ std::function <double(v_t, const X_t&, v_t, const X_t&)> Z; // coupling between sites
+#else
+ std::function <double(const X_t&, const X_t&)> Z; // coupling between sites
+#endif
+
+#ifndef WOLFF_NO_FIELD
+#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
+
+ wolff_system(D_t D, L_t L, 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
+ , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), G(D, L, lat), 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_ext();
+#endif
+#ifdef WOLFF_FINITE_STATES
+ finite_states_init(*this);
+#endif
+ }
+};
+
+#endif
+
diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h
index ec9efaf..68bd09e 100644
--- a/lib/include/wolff/types.h
+++ b/lib/include/wolff/types.h
@@ -1,36 +1,27 @@
#include <inttypes.h>
-typedef uint_fast32_t v_t;
-typedef uint_fast8_t q_t;
-typedef uint_fast16_t R_t;
-typedef uint_fast8_t D_t;
-typedef uint_fast16_t L_t;
-typedef uint_fast64_t count_t;
-typedef int_fast64_t h_t;
+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_R UINT_FAST16_MAX
#define MAX_D UINT_FAST8_MAX
#define MAX_L UINT_FAST16_MAX
-#define MAX_COUNT UINT_FAST64_MAX
-#define MAX_H INT_FAST64_MAX
-#define MIN_H INT_FAST64_MIN
+#define MAX_N UINT_FAST64_MAX
#define PRIv PRIuFAST32
#define PRIq PRIuFAST8
-#define PRIR PRIuFAST16
#define PRID PRIuFAST8
#define PRIL PRIuFAST16
-#define PRIcount PRIuFAST64
-#define PRIh PRIdFAST64
+#define PRIN PRIuFAST64
#define SCNv SCNuFAST32
#define SCNq SCNuFAST8
-#define SCNR SCNuFAST16
#define SCND SCNuFAST8
#define SCNL SCNuFAST16
-#define SCNcount SCNuFAST64
-#define SCNh SCNdFAST64
+#define SCNN SCNuFAST64
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index c6f0ba6..b42aa78 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -7,11 +7,11 @@ graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
nv = pow(L, D);
ne = D * nv;
- v_adj.resize(nv);
+ adj.resize(nv);
coordinate.resize(nv);
- for (std::vector<v_t> v_adj_i : v_adj) {
- v_adj_i.reserve(2 * D);
+ for (std::vector<v_t> adj_i : adj) {
+ adj_i.reserve(2 * D);
}
for (v_t i = 0; i < nv; i++) {
@@ -19,21 +19,21 @@ graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
for (D_t j = 0; j < D; j++) {
coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
- v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
- v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
+ adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
+ adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
}
}
break;
}
case DIAGONAL_LATTICE: {
nv = D * pow(L, D);
- ne = D * nv;
+ ne = 2 * nv;
- v_adj.resize(nv);
+ adj.resize(nv);
coordinate.resize(nv);
- for (std::vector<v_t> v_adj_i : v_adj) {
- v_adj_i.reserve(4 * (D - 1));
+ for (std::vector<v_t> adj_i : adj) {
+ adj_i.reserve(4 * (D - 1));
}
for (D_t i = 0; i < D; i++) {
@@ -42,10 +42,10 @@ graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
for (v_t j = 0; j < pow(L, D); j++) {
v_t vc = sb + j;
- v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + j);
- v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L);
- v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
- v_adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
+ adj[vc].push_back(((i + 1) % D) * pow(L, D) + j);
+ adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L);
+ adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
+ adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
}
}
break;
@@ -54,16 +54,16 @@ graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
}
void graph_t::add_ext() {
- for (std::vector<v_t>& v_adj_i : v_adj) {
- v_adj_i.push_back(nv);
+ for (std::vector<v_t>& adj_i : adj) {
+ adj_i.push_back(nv);
}
- v_adj.resize(nv + 1);
+ adj.resize(nv + 1);
coordinate.resize(nv + 1);
- v_adj[nv].reserve(nv);
+ adj[nv].reserve(nv);
for (v_t i = 0; i < nv; i++) {
- v_adj[nv].push_back(i);
+ adj[nv].push_back(i);
}
coordinate[nv].resize(coordinate[0].size());