diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/angle.h | 48 | ||||
| -rw-r--r-- | lib/circle_group.h | 46 | ||||
| -rw-r--r-- | lib/cluster.h | 16 | ||||
| -rw-r--r-- | lib/dihedral.h | 1 | ||||
| -rw-r--r-- | lib/state.h | 6 | ||||
| -rw-r--r-- | lib/torus.h | 64 | ||||
| -rw-r--r-- | lib/vector.h | 4 |
7 files changed, 172 insertions, 13 deletions
diff --git a/lib/angle.h b/lib/angle.h new file mode 100644 index 0000000..c3f128e --- /dev/null +++ b/lib/angle.h @@ -0,0 +1,48 @@ +#pragma once + +#include "types.h" + +#include <cmath> +#include "vector.h" + +class angle_t { + public: + double x; + + typedef vector_t<2, double> M_t; + typedef vector_t<2, double> F_t; + + angle_t() : x(0) {} + angle_t(double x) : x(x) {} + + inline vector_t<2, double> operator*(v_t a) const { + vector_t<2, double>M; + M[0] = a * cos(x); + M[1] = a * sin(x); + + return M; + } + + inline vector_t<2, double> operator*(double a) const { + vector_t<2, double>M; + M[0] = a * cos(x); + M[1] = a * sin(x); + + return M; + } +}; + +inline vector_t<2,double>& operator+=(vector_t<2,double>& M, const angle_t& theta) { + M[0] += cos(theta.x); + M[1] += sin(theta.x); + + return M; +} + +inline vector_t<2,double>& operator-=(vector_t<2,double>& M, const angle_t& theta) { + M[0] -= cos(theta.x); + M[1] -= sin(theta.x); + + return M; +} + diff --git a/lib/circle_group.h b/lib/circle_group.h new file mode 100644 index 0000000..cb7cadd --- /dev/null +++ b/lib/circle_group.h @@ -0,0 +1,46 @@ +#pragma once + +#include "angle.h" + +class circle_group_t { + public: + bool is_reflection; + double x; + + circle_group_t() : is_reflection(false), x(0) {} + circle_group_t(bool x, double y) : is_reflection(x), x(y) {} + + angle_t act(const angle_t& theta) const { + if (is_reflection) { + return angle_t(fmod(2 * M_PI + x - theta.x, 2 * M_PI)); + } else { + return angle_t(fmod(x + theta.x, 2 * M_PI)); + } + } + + circle_group_t act(const circle_group_t& r) const { + if (is_reflection) { + return circle_group_t(!r.is_reflection, fmod(2 * M_PI + x - r.x, 2 * M_PI)); + } else { + return circle_group_t(r.is_reflection, fmod(x + r.x, 2 * M_PI)); + } + } + + angle_t act_inverse(const angle_t& theta) const { + if (is_reflection) { + return act(theta); + } else { + return angle_t(fmod(2 * M_PI + theta.x - x, 2 * M_PI)); + } + } + + circle_group_t act_inverse(const circle_group_t& r) const { + if (is_reflection) { + return act(r); + } else { + return circle_group_t(r.is_reflection, fmod(2 * M_PI + r.x - x, 2 * M_PI)); + } + } +}; + + diff --git a/lib/cluster.h b/lib/cluster.h index c5f2be7..f948586 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -31,7 +31,9 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { R_old = state->R; marks[v] = true; - if (v == state->nv) { + bool v_is_ghost = (v == state->nv); + + if (v_is_ghost) { R_new = r.act(R_old); } else { si_old = state->spins[v]; @@ -40,19 +42,18 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { for (const v_t &vn : state->g.v_adj[v]) { X_t sj; + bool vn_is_ghost = (vn == state->nv); - if (vn != state->nv) { + if (!vn_is_ghost) { sj = state->spins[vn]; } double prob; - bool is_ext = (v == state->nv || vn == state->nv); - - if (is_ext) { + if (v_is_ghost || vn_is_ghost) { X_t rs_old, rs_new; v_t non_ghost; - if (vn == state->nv) { + if (vn_is_ghost) { rs_old = R_old.act_inverse(si_old); rs_new = R_old.act_inverse(si_new); non_ghost = v; @@ -61,6 +62,7 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { 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[state_to_ind(rs_old)][state_to_ind(rs_new)]; @@ -97,7 +99,7 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { } } - if (v == state->nv) { + if (v_is_ghost) { state->R = R_new; } else { state->spins[v] = si_new; diff --git a/lib/dihedral.h b/lib/dihedral.h index 5588afa..8d0472b 100644 --- a/lib/dihedral.h +++ b/lib/dihedral.h @@ -21,7 +21,6 @@ class dihedral_t { } } - 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); diff --git a/lib/state.h b/lib/state.h index eefa0cb..5abf65b 100644 --- a/lib/state.h +++ b/lib/state.h @@ -27,10 +27,10 @@ class state_t { std::vector<double> precomputed_cos; std::vector<double> precomputed_sin; - std::function <double(X_t, X_t)> J; - std::function <double(X_t)> H; + std::function <double(const X_t&, const X_t&)> J; + std::function <double(const 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), R(), J(J), H(H) { + state_t(D_t D, L_t L, double T, std::function <double(const X_t&, const X_t&)> J, std::function <double(const 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(); diff --git a/lib/torus.h b/lib/torus.h new file mode 100644 index 0000000..2aead52 --- /dev/null +++ b/lib/torus.h @@ -0,0 +1,64 @@ + +#pragma once + +#include <cmath> +#include <array> +#include "types.h" + +template <q_t n> +class torus_t : public std::array<double, n> { + public: + typedef std::array<double, n> M_t; + typedef std::array<double, n> F_t; + + torus_t() { + this->fill(0); + } + + inline torus_t<n> operator*(v_t a) const { + torus_t<n> x; + for (q_t i = 0; i < n; i++) { + x[i] = a * (*this)[i]; + } + + return x; + } + + inline torus_t<n> operator*(double a) const { + torus_t<n> x; + for (q_t i = 0; i < n; i++) { + x[i] = a * (*this)[i]; + } + + return x; + } + + inline torus_t<n>& operator+=(const torus_t<n>& x) { + for (q_t i = 0; i < n; i++) { + (*this)[i] += x[i]; + } + } + + inline torus_t<n>& operator-=(const torus_t<n>& x) { + for (q_t i = 0; i < n; i++) { + (*this)[i] -= x[i]; + } + } +}; + +template <q_t n> +double norm_squared(const torus_t<n>& x) { + double tmp = 0; + for (const double& xi : x) { + tmp += pow(xi, 2); + } + return tmp; +} + +void write_magnetization(const torus_t<n>& x, FILE *outfile) { + for (const double& xi : x) { + float tmp_xi = (float)xi; + fwrite(&tmp_xi, sizeof(float), 1, outfile); + } +} + diff --git a/lib/vector.h b/lib/vector.h index beee1a7..2f4077a 100644 --- a/lib/vector.h +++ b/lib/vector.h @@ -109,7 +109,7 @@ void write_magnetization(vector_t <q, double> M, FILE *outfile) { } template <q_t q, class T> -T dot(vector_t <q, T> v1, vector_t <q, T> v2) { +T dot(const vector_t <q, T>& v1, const vector_t <q, T>& v2) { T prod = 0; for (q_t i = 0; i < q; i++) { @@ -120,7 +120,7 @@ T dot(vector_t <q, T> v1, vector_t <q, T> v2) { } template <q_t q, class T> -double H_vector(vector_t <q, T> v1, T *H) { +double H_vector(const vector_t <q, T>& v1, T *H) { vector_t <q, T> H_vec(H); return (double)(dot <q, T> (v1, H_vec)); } |