summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/cluster.h22
-rw-r--r--lib/dihedral.h112
-rw-r--r--lib/dihedral_inf.h97
-rw-r--r--lib/finite_states.h9
-rw-r--r--lib/orthogonal.h200
-rw-r--r--lib/state.h7
-rw-r--r--lib/symmetric.h98
-rw-r--r--lib/wolff.h1
-rw-r--r--lib/z2.h60
-rw-r--r--src/analyze_correlations.cpp3
-rw-r--r--src/wolff_On.cpp2
-rw-r--r--src/wolff_cgm.cpp4
-rw-r--r--src/wolff_clock.cpp10
-rw-r--r--src/wolff_dgm.cpp4
-rw-r--r--src/wolff_ising.cpp6
-rw-r--r--src/wolff_potts.cpp7
16 files changed, 247 insertions, 395 deletions
diff --git a/lib/cluster.h b/lib/cluster.h
index 427c3c8..c5f2be7 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -32,13 +32,13 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
marks[v] = true;
if (v == state->nv) {
- R_new = act (r, R_old);
+ R_new = r.act(R_old);
} else {
si_old = state->spins[v];
- si_new = act (r, si_old);
+ si_new = r.act(si_old);
}
- for (v_t vn : state->g.v_adj[v]) {
+ for (const v_t &vn : state->g.v_adj[v]) {
X_t sj;
if (vn != state->nv) {
@@ -53,17 +53,17 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
X_t rs_old, rs_new;
v_t non_ghost;
if (vn == state->nv) {
- rs_old = act_inverse (R_old, si_old);
- rs_new = act_inverse (R_old, si_new);
+ rs_old = R_old.act_inverse(si_old);
+ rs_new = R_old.act_inverse(si_new);
non_ghost = v;
} else {
- rs_old = act_inverse (R_old, sj);
- rs_new = act_inverse (R_new, sj);
+ rs_old = R_old.act_inverse(sj);
+ rs_new = R_new.act_inverse(sj);
non_ghost = vn;
}
double dE = state->H(rs_old) - state->H(rs_new);
#ifdef FINITE_STATES
- prob = H_probs[N_STATES * state_to_ind(rs_old) + state_to_ind(rs_new)];
+ prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
#else
prob = 1.0 - exp(-dE / state->T);
#endif
@@ -85,7 +85,7 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
} else {
double dE = state->J(si_old, sj) - state->J(si_new, sj);
#ifdef FINITE_STATES
- prob = J_probs[N_STATES * N_STATES * state_to_ind(si_old) + N_STATES * state_to_ind(si_new) + state_to_ind(sj)];
+ prob = J_probs[state_to_ind(si_old)][state_to_ind(si_new)][state_to_ind(sj)];
#else
prob = 1.0 - exp(-dE / state->T);
#endif
@@ -98,13 +98,9 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
}
if (v == state->nv) {
- free_spin(state->R);
state->R = R_new;
} else {
state->spins[v] = si_new;
- }
-
- if (v != state->nv) { // count the number of non-external sites that flip
nv++;
}
}
diff --git a/lib/dihedral.h b/lib/dihedral.h
index 8238835..5588afa 100644
--- a/lib/dihedral.h
+++ b/lib/dihedral.h
@@ -4,76 +4,46 @@
#include "types.h"
#include "potts.h"
-template <class T, q_t q>
-struct dihedral_t { bool is_reflection; T x; };
-
-template <class T, q_t q>
-void init(dihedral_t<T, q> *ptr) {
- ptr->is_reflection = false;
- ptr->x = (T)0;
-}
-
-template <class T, q_t q>
-dihedral_t<T, q> copy(dihedral_t<T, q> r) {
- return r;
-}
-
-template <class T, q_t q>
-void free_spin(dihedral_t<T, q> r) {
- // do nothing!
-}
-
-template <q_t q>
-potts_t<q> act(dihedral_t<q_t, q> r, potts_t<q> s) {
- potts_t<q> s2;
- if (r.is_reflection) {
- s2.x = ((q + r.x) - s.x) % q;
- } else {
- s2.x = (r.x + s.x) % q;
- }
-
- return s2;
-}
-
-template <q_t q>
-dihedral_t<q_t,q> act(dihedral_t<q_t,q> r1, dihedral_t<q_t,q> r2) {
- dihedral_t<q_t,q> r3;
-
- if (r1.is_reflection) {
- r3.is_reflection = !(r2.is_reflection);
- r3.x = ((q + r1.x) - r2.x) % q;
- } else {
- r3.is_reflection = r2.is_reflection;
- r3.x = (r1.x + r2.x) % q;
- }
-
- return r3;
-}
-
-template <q_t q>
-potts_t<q> act_inverse(dihedral_t<q_t,q> r, potts_t<q> s) {
- potts_t<q> s2;
- if (r.is_reflection) {
- s2.x = ((r.x + q) - s.x) % q;
- } else {
- s2.x = ((s.x + q) - r.x) % q;
- }
-
- return s2;
-}
-
template <q_t q>
-dihedral_t<q_t, q> act_inverse(dihedral_t<q_t,q> r1, dihedral_t<q_t,q> r2) {
- dihedral_t<q_t,q> r3;
-
- if (r1.is_reflection) {
- r3.is_reflection = !(r2.is_reflection);
- r3.x = ((r1.x + q) - r2.x) % q;
- } else {
- r3.is_reflection = r2.is_reflection;
- r3.x = ((r2.x + q) - r1.x) % q;
- }
-
- return r3;
-}
+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/dihedral_inf.h b/lib/dihedral_inf.h
index f7c4297..a064b48 100644
--- a/lib/dihedral_inf.h
+++ b/lib/dihedral_inf.h
@@ -4,75 +4,44 @@
#include "height.h"
template <class T>
-struct dihedral_inf_t { bool is_reflection; T x; };
-
-template <class T>
-void init(dihedral_inf_t<T> *ptr) {
- ptr->is_reflection = false;
- ptr->x = (T)0;
-}
-
-template <class T>
-dihedral_inf_t<T> copy(dihedral_inf_t<T> r) {
- return r;
-}
-
-template <class T>
-void free_spin(dihedral_inf_t<T> r) {
- // do nothing!
-}
-
-template <class T>
-height_t<T> act(dihedral_inf_t<T> r, height_t<T> h) {
- height_t<T> h2;
- if (r.is_reflection) {
- h2.x = r.x - h.x;
- } else {
- h2.x = r.x + h.x;
+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);
+ }
}
- return h2;
-}
-
-template <class T>
-dihedral_inf_t<T> act(dihedral_inf_t<T> r1, dihedral_inf_t<T> r2) {
- dihedral_inf_t<T> r3;
-
- if (r1.is_reflection) {
- r3.is_reflection = !(r2.is_reflection);
- r3.x = r1.x - r2.x;
- } else {
- r3.is_reflection = r2.is_reflection;
- r3.x = r1.x + r2.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);
+ }
}
- return r3;
-}
-
-template <class T>
-height_t<T> act_inverse(dihedral_inf_t<T> r, height_t<T> h) {
- height_t<T> h2;
- if (r.is_reflection) {
- h2.x = r.x - h.x;
- } else {
- h2.x = h.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);
+ }
}
- return h2;
-}
-
-template <class T>
-dihedral_inf_t<T> act_inverse(dihedral_inf_t<T> r1, dihedral_inf_t<T> r2) {
- dihedral_inf_t<T> r3;
-
- if (r1.is_reflection) {
- r3.is_reflection = !(r2.is_reflection);
- r3.x = r1.x - r2.x;
- } else {
- r3.is_reflection = r2.is_reflection;
- r3.x = r2.x - r1.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);
+ }
}
-
- return r3;
-}
+};
diff --git a/lib/finite_states.h b/lib/finite_states.h
index 1e8800c..08eff30 100644
--- a/lib/finite_states.h
+++ b/lib/finite_states.h
@@ -3,27 +3,28 @@
#include <cmath>
#include <functional>
+#include <array>
#define FINITE_STATES
// must have N_STATES, states[N_STATES], and state_to_ind defined before
// invoking header
-double J_probs[N_STATES * N_STATES * N_STATES];
-double H_probs[N_STATES * N_STATES];
+std::array<std::array<std::array<double, N_STATES>, N_STATES>, N_STATES> J_probs;
+std::array<std::array<double, N_STATES>, N_STATES> H_probs;
template <class X_t>
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++) {
for (q_t k = 0; k < N_STATES; k++) {
- J_probs[i * N_STATES * N_STATES + j * N_STATES +k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T);
+ 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 < N_STATES; i++) {
for (q_t j = 0; j < N_STATES; j++) {
- H_probs[i * N_STATES + j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T);
+ H_probs[i][j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T);
}
}
}
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index c0958f2..34cd44e 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -11,164 +11,134 @@
#include "vector.h"
template <q_t q, class T>
-struct orthogonal_t { bool is_reflection; T *x; };
+class orthogonal_t : public std::array<std::array<T, q>, q> {
+ public :
+ bool is_reflection;
-template <q_t q, class T>
-void init(orthogonal_t <q, T> *ptr) {
- ptr->is_reflection = false;
- ptr->x = (T *)calloc(q * q, sizeof(T));
-
- for (q_t i = 0; i < q; i++) {
- ptr->x[q * i + i] = (T)1;
- }
-}
-
-template <q_t q, class T>
-orthogonal_t <q, T> copy (orthogonal_t <q, T> m) {
- orthogonal_t <q, T> m_copy;
- m_copy.is_reflection = m.is_reflection;
-
- q_t size;
-
- if (m.is_reflection) {
- size = q;
- } else {
- size = q * q;
- }
-
- m_copy.x = (T *)calloc(size, sizeof(T));
-
- for (q_t i = 0; i < size; i++) {
- m_copy.x[i] = m.x[i];
+ orthogonal_t() : is_reflection(false) {
+ for (q_t i = 0; i < q; i++) {
+ (*this)[i].fill(0);
+ (*this)[i][i] = (T)1;
+ }
}
- return m_copy;
-}
-
-template <q_t q, class T>
-void free_spin (orthogonal_t <q, T> m) {
- free(m.x);
-}
-
-template <q_t q, class T>
-vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) {
- vector_t <q, T> v_rot;
+ vector_t<q, T> act(const vector_t <q, T>& v) const {
+ vector_t <q, T> v_rot;
- if (m.is_reflection) {
- double prod = 0;
- for (q_t i = 0; i < q; i++) {
- prod += v[i] * m.x[i];
- }
- for (q_t i = 0; i < q; i++) {
- v_rot[i] = v[i] - 2 * prod * m.x[i];
- }
- } else {
- for (q_t i = 0; i < q; i++) {
- for (q_t j = 0; j < q; j++) {
- v_rot[i] += m.x[q * i + j] * v[j];
+ 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;
-}
+ return v_rot;
+ }
-template <q_t q, class T>
-orthogonal_t <q, T> act (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
- orthogonal_t <q, T> m2_rot;
+ orthogonal_t<q, T> act(const orthogonal_t <q, T>& m) const {
+ orthogonal_t <q, T> m_rot;
- m2_rot.is_reflection = false;
- m2_rot.x = (T *)calloc(q * q, sizeof(T));
+ m_rot.is_reflection = false;
- if (m1.is_reflection) {
- for (q_t i = 0; i < q; i++) {
- double akOki = 0;
+ if (is_reflection) {
+ for (q_t i = 0; i < q; i++) {
+ double akOki = 0;
- for (q_t k = 0; k < q; k++) {
- akOki += m1.x[k] * m2.x[q * k + i];
- }
+ for (q_t k = 0; k < q; k++) {
+ akOki += (*this)[0][k] * m[k][i];
+ }
- for (q_t j = 0; j < q; j++) {
- m2_rot.x[q * j + i] = m2.x[q * j + i] - 2 * akOki * m1.x[j];
+ 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++) {
- for (q_t j = 0; j < q; j++) {
- for (q_t k = 0; k < q; k++) {
- m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k];
+ } 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 m2_rot;
-}
+ return m_rot;
+ }
-template <q_t q, class T>
-vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) {
- if (m.is_reflection) {
- return act(m, v); // reflections are their own inverse
- } else {
- vector_t <q, T> v_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;
- for (q_t i = 0; i < q; i++) {
- for (q_t j = 0; j < q; j++) {
- v_rot[i] += m.x[q * j + i] * v[j];
+ 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;
+ return v_rot;
+ }
}
-}
-
-template <q_t q, class T>
-orthogonal_t <q, T> act_inverse (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
- if (m1.is_reflection) {
- return act(m1, m2); // reflections are their own inverse
- } else {
- orthogonal_t <q, T> m2_rot;
- m2_rot.x = (T *)calloc(q * q, sizeof(T));
- for (q_t i = 0; i < q; i++) {
- for (q_t j = 0; j < q; j++) {
- for (q_t k = 0; k < q; k++) {
- m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k];
+ 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 m2_rot;
+ return m_rot;
+ }
}
-}
+
+};
+
template <q_t q>
-orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, vector_t <q, double> v) {
+orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, const vector_t <q, double>& v) {
orthogonal_t <q, double> ptr;
ptr.is_reflection = true;
- ptr.x = (double *)calloc(q, sizeof(double));
double v2 = 0;
for (q_t i = 0; i < q; i++) {
- ptr.x[i] = gsl_ran_ugaussian(r);
- v2 += ptr.x[i] * ptr.x[i];
+ ptr[0][i] = gsl_ran_ugaussian(r);
+ v2 += ptr[0][i] * ptr[0][i];
}
double mag_v = sqrt(v2);
for (q_t i = 0; i < q; i++) {
- ptr.x[i] /= mag_v;
+ ptr[0][i] /= mag_v;
}
return ptr;
}
template <q_t q>
-orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q, double> v0, double epsilon, unsigned int n) {
+orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const vector_t <q, double>& v0, double epsilon, unsigned int n) {
orthogonal_t <q, double> m;
m.is_reflection = true;
- m.x = (double *)malloc(q * sizeof(double));
vector_t <q, double> v;
@@ -191,22 +161,22 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q
double m_dot_v = 0;
for (q_t i = 0; i < q; i++) {
- m.x[i] = gsl_ran_ugaussian(r);
- m_dot_v += m.x[i] * v[i];
+ m[0][i] = gsl_ran_ugaussian(r);
+ m_dot_v += m[0][i] * v[i];
}
double v2 = 0;
double factor = epsilon * gsl_ran_ugaussian(r);
for (q_t i = 0; i < q; i++) {
- m.x[i] = m.x[i] - m_dot_v * v[i] + factor * v[i];
- v2 += pow(m.x[i], 2);
+ m[0][i] = m[0][i] - m_dot_v * v[i] + factor * v[i];
+ v2 += pow(m[0][i], 2);
}
double mag_v = sqrt(v2);
for (q_t i = 0; i < q; i++) {
- m.x[i] /= mag_v;
+ m[0][i] /= mag_v;
}
return m;
diff --git a/lib/state.h b/lib/state.h
index b57a9e2..eefa0cb 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -30,12 +30,11 @@ class state_t {
std::function <double(X_t, X_t)> J;
std::function <double(X_t)> H;
- state_t(D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H) : D(D), L(L), g(D, L), T(T), J(J), H(H) {
+ state_t(D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H) : D(D), L(L), g(D, L), T(T), R(), J(J), H(H) {
nv = g.nv;
ne = g.ne;
g.add_ext();
spins.resize(nv);
- init (&R);
E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
M = spins[0] * nv;
last_cluster_size = 0;
@@ -52,10 +51,6 @@ class state_t {
precomputed_sin[i] = sin(2 * M_PI * (double)i / (double)L);
}
}
-
- ~state_t() {
- free_spin(R);
- }
};
diff --git a/lib/symmetric.h b/lib/symmetric.h
index 41c8fb6..9e9b6e4 100644
--- a/lib/symmetric.h
+++ b/lib/symmetric.h
@@ -2,84 +2,50 @@
#pragma once
#include <stdlib.h>
+#include <array>
#include "types.h"
#include "potts.h"
template <q_t q>
-class symmetric_t {
+class symmetric_t : public std::array<q_t, q> {
public:
- q_t *perm;
-};
-
-template <q_t q>
-void init(symmetric_t<q> *p) {
- p->perm = (q_t *)malloc(q * sizeof(q_t));
-
- for (q_t i = 0; i < q; i++) {
- p->perm[i] = i;
- }
-}
-
-template <q_t q>
-void free_spin(symmetric_t<q> p) {
- free(p.perm);
-}
-template <q_t q>
-symmetric_t<q> copy(symmetric_t<q> x) {
- symmetric_t<q> x2;
- x2.perm = (q_t *)malloc(q * sizeof(q_t));
-
- for (q_t i = 0; i < q; i++) {
- x2.perm[i] = x.perm[i];
- }
-
- return x2;
-}
-
-template <q_t q>
-potts_t<q> act(symmetric_t<q> r, potts_t<q> s) {
- potts_t<q> s2;
- s2.x = r.perm[s.x];
- return s2;
-}
-
-template <q_t q>
-symmetric_t<q> act(symmetric_t<q> r1, symmetric_t<q> r2) {
- symmetric_t<q> r3;
- r3.perm = (q_t *)malloc(q * sizeof(q_t));
- for (q_t i = 0; i < q; i++) {
- r3.perm[i] = r1.perm[r2.perm[i]];
- }
-
- return r3;
-}
+ symmetric_t() {
+ for (q_t i = 0; i < q; i++) {
+ (*this)[i] = i;
+ }
+ }
-template <q_t q>
-potts_t<q> act_inverse(symmetric_t<q> r, potts_t<q> s) {
- potts_t<q> s2;
+ potts_t<q> act(const potts_t<q> &s) const {
+ return potts_t<q>((*this)[s.x]);
+ }
- q_t i;
+ 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]];
+ }
- for (i = 0; i < q; i++) {
- if (r.perm[i] == s.x) {
- break;
+ return r_rot;
}
- }
- s2.x = i;
+ 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);
+ }
+ }
- return s2;
-}
+ exit(EXIT_FAILURE);
+ }
-template <q_t q>
-symmetric_t<q> act_inverse(symmetric_t<q> r1, symmetric_t<q> r2) {
- symmetric_t<q> r3;
- r3.perm = (q_t *)malloc(q * sizeof(q_t));
- for (q_t i = 0; i < q; i++) {
- r3.perm[r1.perm[i]] = r2.perm[i];
- }
+ 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 r3;
-}
+ return r_rot;
+ }
+};
diff --git a/lib/wolff.h b/lib/wolff.h
index da4b7b6..a4a663c 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -16,7 +16,6 @@ void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s
v_t v0 = gsl_rng_uniform_int(r, s->nv);
R_t step = gen_R(r, s->spins[v0]);
flip_cluster <R_t, X_t> (s, v0, step, r);
- free_spin(step);
measurements(s);
}
diff --git a/lib/z2.h b/lib/z2.h
index dae1eb7..fd4722e 100644
--- a/lib/z2.h
+++ b/lib/z2.h
@@ -17,49 +17,37 @@
*
*/
-struct z2_t { bool x; };
+class z2_t {
+ public:
+ bool x;
-void init(z2_t *p) {
- p->x = false;
-}
+ z2_t() : x(false) {}
-void free_spin(z2_t p) {
- // do nothing!
-}
+ z2_t(bool x) : x(x) {}
-z2_t copy(z2_t x) {
- return x;
-}
-
-ising_t act(z2_t r, ising_t s) {
- ising_t rs;
-
- if (r.x) {
- rs.x = !s.x;
- return rs;
- } else {
- rs.x = s.x;
- return rs;
+ ising_t act(const ising_t& s) {
+ if (x) {
+ return ising_t(!s.x);
+ } else {
+ return ising_t(s.x);
+ }
}
-}
-z2_t act(z2_t r1, z2_t r2) {
- z2_t r3;
+ z2_t act(const z2_t& r) {
+ if (x) {
+ return z2_t(!r.x);
+ } else {
+ return z2_t(r.x);
+ }
+ }
- if (r1.x) {
- r3.x = !r2.x;
- return r3;
- } else {
- r3.x = r2.x;
- return r3;
+ ising_t act_inverse(const ising_t& s) {
+ return this->act(s);
}
-}
-ising_t act_inverse(z2_t r, ising_t s) {
- return act(r, s);
-}
+ z2_t act_inverse(const z2_t& r) {
+ return this->act(r);
+ }
+};
-z2_t act_inverse(z2_t r1, z2_t r2) {
- return act(r1, r2);
-}
diff --git a/src/analyze_correlations.cpp b/src/analyze_correlations.cpp
index e108f98..5f3dda0 100644
--- a/src/analyze_correlations.cpp
+++ b/src/analyze_correlations.cpp
@@ -1,6 +1,7 @@
#include <types.h>
-#include <measurement.h>
+#include <cmath>
+#include <cstring>
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index 15a9553..5458860 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -189,7 +189,7 @@ int main(int argc, char *argv[]) {
other_f = [&] (const On_t *s) {
glClear(GL_COLOR_BUFFER_BIT);
for (v_t i = 0; i < pow(L, 2); i++) {
- vector_R_t v_tmp = act_inverse(s->R, s->spins[i]);
+ vector_R_t v_tmp = s->R.act_inverse(s->spins[i]);
double thetai = fmod(2 * M_PI + theta(v_tmp), 2 * M_PI);
double saturation = 0.7;
double value = 0.9;
diff --git a/src/wolff_cgm.cpp b/src/wolff_cgm.cpp
index 10477bd..7613691 100644
--- a/src/wolff_cgm.cpp
+++ b/src/wolff_cgm.cpp
@@ -128,7 +128,7 @@ int main(int argc, char *argv[]) {
double max_h = INT64_MIN;
double min_h = INT64_MAX;
for (v_t i = 0; i < pow(L, 2); i++) {
- double cur_h = act_inverse(s->R, s->spins[i]).x;
+ double cur_h = (s->R.act_inverse(s->spins[i])).x;
if (cur_h < min_h) {
min_h = cur_h;
}
@@ -138,7 +138,7 @@ int main(int argc, char *argv[]) {
}
for (v_t i = 0; i < pow(L, 2); i++) {
- double cur_h = act_inverse(s->R, s->spins[i]).x;
+ double cur_h = (s->R.act_inverse(s->spins[i])).x;
double mag = ((double)(cur_h - min_h)) / ((double)(max_h - min_h));
glColor3f(mag, mag, mag);
glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1);
diff --git a/src/wolff_clock.cpp b/src/wolff_clock.cpp
index e186c44..bc2c5d1 100644
--- a/src/wolff_clock.cpp
+++ b/src/wolff_clock.cpp
@@ -18,7 +18,7 @@
#include <rand.h>
#include <wolff.h>
-typedef state_t <dihedral_t<q_t,POTTSQ>, potts_t<POTTSQ>> sim_t;
+typedef state_t <dihedral_t<POTTSQ>, potts_t<POTTSQ>> sim_t;
int main(int argc, char *argv[]) {
@@ -88,11 +88,11 @@ int main(int argc, char *argv[]) {
};
// initialize state object
- state_t <dihedral_t<q_t,POTTSQ>, potts_t<POTTSQ>> s(D, L, T, Z, B);
+ state_t <dihedral_t<POTTSQ>, potts_t<POTTSQ>> s(D, L, T, Z, B);
// define function that generates self-inverse rotations
- std::function <dihedral_t<q_t,POTTSQ>(gsl_rng *, potts_t<POTTSQ>)> gen_R = [] (gsl_rng *r, potts_t<POTTSQ> v) -> dihedral_t<q_t,POTTSQ> {
- dihedral_t<q_t,POTTSQ> rot;
+ std::function <dihedral_t<POTTSQ>(gsl_rng *, potts_t<POTTSQ>)> gen_R = [] (gsl_rng *r, potts_t<POTTSQ> v) -> dihedral_t<POTTSQ> {
+ dihedral_t<POTTSQ> rot;
rot.is_reflection = true;
q_t x = gsl_rng_uniform_int(r, POTTSQ - 1);
rot.x = (2 * v.x + x + 1) % POTTSQ;
@@ -127,7 +127,7 @@ int main(int argc, char *argv[]) {
average_M += (double)s->M[0] / (double)N / (double)s->nv;
glClear(GL_COLOR_BUFFER_BIT);
for (v_t i = 0; i < pow(L, 2); i++) {
- potts_t<POTTSQ> tmp_s = act_inverse(s->R, s->spins[i]);
+ potts_t<POTTSQ> tmp_s = s->R.act_inverse(s->spins[i]);
glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ));
glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1);
}
diff --git a/src/wolff_dgm.cpp b/src/wolff_dgm.cpp
index 6271b79..2583704 100644
--- a/src/wolff_dgm.cpp
+++ b/src/wolff_dgm.cpp
@@ -133,7 +133,7 @@ int main(int argc, char *argv[]) {
int64_t max_h = INT64_MIN;
int64_t min_h = INT64_MAX;
for (v_t i = 0; i < pow(L, 2); i++) {
- int64_t cur_h = act_inverse(s->R, s->spins[i]).x;
+ int64_t cur_h = (s->R.act_inverse(s->spins[i])).x;
if (cur_h < min_h) {
min_h = cur_h;
}
@@ -143,7 +143,7 @@ int main(int argc, char *argv[]) {
}
for (v_t i = 0; i < pow(L, 2); i++) {
- int64_t cur_h = act_inverse(s->R, s->spins[i]).x;
+ int64_t cur_h = (s->R.act_inverse(s->spins[i])).x;
double mag = ((double)(cur_h - min_h)) / ((double)(max_h - min_h));
glColor3f(mag, mag, mag);
glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1);
diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp
index 93b7ea4..410e046 100644
--- a/src/wolff_ising.cpp
+++ b/src/wolff_ising.cpp
@@ -120,10 +120,8 @@ int main(int argc, char *argv[]) {
state_t <z2_t, ising_t> s(D, L, T, Z, B);
// define function that generates self-inverse rotations
- std::function <z2_t(gsl_rng *, ising_t)> gen_R = [] (gsl_rng *, ising_t s) -> z2_t {
- z2_t rot;
- rot.x = true;
- return rot;
+ std::function <z2_t(gsl_rng *, ising_t)> gen_R = [] (gsl_rng *, const ising_t& s) -> z2_t {
+ return z2_t(true);
};
FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp
index cdc4c07..07663d5 100644
--- a/src/wolff_potts.cpp
+++ b/src/wolff_potts.cpp
@@ -118,7 +118,6 @@ int main(int argc, char *argv[]) {
// define function that generates self-inverse rotations
std::function <symmetric_t<POTTSQ>(gsl_rng *, potts_t<POTTSQ>)> gen_R = [] (gsl_rng *r, potts_t<POTTSQ> v) -> symmetric_t<POTTSQ> {
symmetric_t<POTTSQ> rot;
- init(&rot);
q_t j = gsl_rng_uniform_int(r, POTTSQ - 1);
q_t swap_v;
@@ -128,8 +127,8 @@ int main(int argc, char *argv[]) {
swap_v = j + 1;
}
- rot.perm[v.x] = swap_v;
- rot.perm[swap_v] = v.x;
+ rot[v.x] = swap_v;
+ rot[swap_v] = v.x;
return rot;
};
@@ -158,7 +157,7 @@ int main(int argc, char *argv[]) {
other_f = [] (const sim_t *s) {
glClear(GL_COLOR_BUFFER_BIT);
for (v_t i = 0; i < pow(s->L, 2); i++) {
- potts_t<POTTSQ> tmp_s = act_inverse(s->R, s->spins[i]);
+ potts_t<POTTSQ> tmp_s = s->R.act_inverse(s->spins[i]);
glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ));
glRecti(i / s->L, i % s->L, (i / s->L) + 1, (i % s->L) + 1);
}