diff options
| -rw-r--r-- | space_wolff.hpp | 116 |
1 files changed, 61 insertions, 55 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp index 4c02406..f95c68a 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -29,12 +29,13 @@ const std::array<std::array<unsigned, 16>, 16> smiley = { {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}}, {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}}}; -template <class U, unsigned D> using vector = Eigen::Matrix<U, D, 1>; +template <class U, unsigned D> using Vector = Eigen::Matrix<U, D, 1>; -template <class U, unsigned D> using matrix = Eigen::Matrix<U, D, D>; +template <class U, unsigned D> using Matrix = Eigen::Matrix<U, D, D>; -template <class U, unsigned D> vector<U, D> diff(U L, vector<U, D> v1, vector<U, D> v2) { - vector<U, D> v; +template <class U, unsigned D> +Vector<U, D> diff(U L, Vector<U, D> v1, Vector<U, D> v2) { + Vector<U, D> v; for (unsigned i = 0; i < D; i++) { v(i) = std::abs(v1(i) - v2(i)); @@ -46,20 +47,20 @@ template <class U, unsigned D> vector<U, D> diff(U L, vector<U, D> v1, vector<U, return v; } -template <class U, unsigned D, class state> class spin { +template <class U, unsigned D, class state> class Spin { public: - vector<U, D> x; + Vector<U, D> x; state s; }; -template <class U, unsigned D> class euclidean { +template <class U, unsigned D> class Euclidean { private: U L; public: - vector<U, D> t; - matrix<U, D> r; - euclidean(U L) : L(L) { + Vector<U, D> t; + Matrix<U, D> r; + Euclidean(U L) : L(L) { for (unsigned i = 0; i < D; i++) { t(i) = 0; r(i, i) = 1; @@ -69,14 +70,14 @@ public: } } - euclidean(U L, vector<U, D> t0, matrix<U, D> r0) : L(L) { + Euclidean(U L, Vector<U, D> t0, Matrix<U, D> r0) : L(L) { t = t0; r = r0; } template <class state> - spin<U, D, state> act(const spin<U, D, state> &s) const { - spin<U, D, state> s_new; + Spin<U, D, state> act(const Spin<U, D, state> &s) const { + Spin<U, D, state> s_new; s_new.x = t + r * s.x; s_new.s = s.s; @@ -88,39 +89,39 @@ public: return s_new; } - euclidean act(const euclidean& x) const { - vector<U, D> tnew = r * x.t + t; - matrix<U, D> rnew = r * x.r; + Euclidean act(const Euclidean &x) const { + Vector<U, D> tnew = r * x.t + t; + Matrix<U, D> rnew = r * x.r; for (unsigned i = 0; i < D; i++) { tnew(i) = fmod(L + tnew(i), L); } - euclidean pnew(this->L, tnew, rnew); + Euclidean pnew(this->L, tnew, rnew); return pnew; } - euclidean inverse() const { - vector<U, D> tnew = -r.transpose() * t; - matrix<U, D> rnew = r.transpose(); + Euclidean inverse() const { + Vector<U, D> tnew = -r.transpose() * t; + Matrix<U, D> rnew = r.transpose(); - euclidean pnew(this->L, tnew, rnew); + Euclidean pnew(this->L, tnew, rnew); return pnew; } }; -template <class T, unsigned D> class dictionary { +template <class T, unsigned D> class Dictionary { private: unsigned N; T L; std::vector<std::set<unsigned>> d; public: - dictionary(unsigned Ni, double Li) : N(Ni), L(Li), d(pow(Ni, D)) {}; + Dictionary(unsigned Ni, double Li) : N(Ni), L(Li), d(pow(Ni, D)){}; - unsigned dictionary_index(vector<T, D> x) const { + unsigned dictionary_index(Vector<T, D> x) const { unsigned pos_ind = 0; for (unsigned i = 0; i < D; i++) { @@ -130,15 +131,15 @@ public: return pos_ind; } - void record(vector<T, D> x, unsigned ind) { + void record(Vector<T, D> x, unsigned ind) { d[this->dictionary_index(x)].insert(ind); }; - void remove(vector<T, D> x, unsigned ind) { + void remove(Vector<T, D> x, unsigned ind) { d[this->dictionary_index(x)].erase(ind); }; - std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const { + std::set<unsigned> neighbors(Vector<T, D> x, unsigned depth) const { return nearest_neighbors_of(this->dictionary_index(x), depth, {}); }; @@ -169,7 +170,7 @@ public: }; }; -class quantity { +class Quantity { private: double total; double total2; @@ -180,7 +181,7 @@ public: unsigned n; std::list<double> hist; - quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) { + Quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) { n = 0; total = 0; total2 = 0; @@ -244,20 +245,22 @@ public: unsigned num_added() const { return n - wait; } }; -template <class U, unsigned D, class state> class model { +template <class U, unsigned D, class state> class Model { public: U L; - euclidean<U, D> s0; - std::vector<spin<U, D, state>> s; - dictionary<U, D> dict; - std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors; - std::function<double(spin<U, D, state>, spin<U, D, state>)> Z; - std::function<double(spin<U, D, state>)> B; - std::vector<matrix<U, D>> mats; - std::vector<vector<U, D>> steps; + Euclidean<U, D> s0; + std::vector<Spin<U, D, state>> s; + Dictionary<U, D> dict; + std::function<std::set<unsigned>(Model<U, D, state> &, unsigned, + Spin<U, D, state>)> + neighbors; + std::function<double(Spin<U, D, state>, Spin<U, D, state>)> Z; + std::function<double(Spin<U, D, state>)> B; + std::vector<Matrix<U, D>> mats; + std::vector<Vector<U, D>> steps; long double E; - quantity Eq; - quantity Cq; + Quantity Eq; + Quantity Cq; void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) { if (level > 0) { @@ -272,10 +275,14 @@ public: } } - model(U L, unsigned N, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z, - std::function<double(spin<U, D, state>)> B, - std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) - : L(L), s0(L), dict(N, L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) { + Model(U L, unsigned N, + std::function<double(Spin<U, D, state>, Spin<U, D, state>)> Z, + std::function<double(Spin<U, D, state>)> B, + std::function<std::set<unsigned>(Model<U, D, state> &, unsigned, + Spin<U, D, state>)> + ns) + : L(L), s0(L), dict(N, L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), + Cq(1000, 1000) { std::array<double, D> ini_sequence; ini_sequence.fill(1); std::list<std::array<double, D>> sequences; @@ -286,7 +293,7 @@ public: sequences.pop_front(); // don't want the identity matrix! for (std::array<double, D> sequence : sequences) { - matrix<U, D> m; + Matrix<U, D> m; for (unsigned i = 0; i < D; i++) { for (unsigned j = 0; j < D; j++) { if (i == j) { @@ -299,7 +306,7 @@ public: mats.push_back(m); - vector<U, D> v; + Vector<U, D> v; for (unsigned i = 0; i < D; i++) { if (sequence[i] == 1) { v(i) = 0; @@ -314,7 +321,7 @@ public: for (unsigned i = 0; i < D; i++) { for (unsigned j = 0; j < D; j++) { if (i != j) { - matrix<U, D> m; + Matrix<U, D> m; for (unsigned k = 0; k < D; k++) { for (unsigned l = 0; l < D; l++) { if ((k == i && l == j) || (k == j && l == i)) { @@ -345,7 +352,7 @@ public: } } - void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) { + void step(double T, unsigned ind, Euclidean<U, D> r, std::mt19937 &rng) { unsigned cluster_size = 0; std::uniform_real_distribution<double> dist(0.0, 1.0); @@ -363,8 +370,8 @@ public: bool we_are_ghost = i == s.size(); - spin<U, D, state> si_new; - euclidean<U, D> s0_new(L); + Spin<U, D, state> si_new; + Euclidean<U, D> s0_new(L); if (we_are_ghost) { s0_new = r.act(s0); @@ -378,7 +385,7 @@ public: bool neighbor_is_ghost = j == s.size(); if (we_are_ghost || neighbor_is_ghost) { - spin<U, D, state> s0s_old, s0s_new; + Spin<U, D, state> s0s_old, s0s_new; unsigned non_ghost; if (neighbor_is_ghost) { @@ -423,8 +430,8 @@ public: std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1); for (unsigned i = 0; i < N; i++) { - vector<U, D> t; - matrix<U, D> m; + Vector<U, D> t; + Matrix<U, D> m; unsigned flip = coin(rng); if (flip < mats.size()) { for (unsigned j = 0; j < D; j++) { @@ -445,7 +452,7 @@ public: t = steps[flip - mats.size()]; } - euclidean<U, D> g(L, t, m); + Euclidean<U, D> g(L, t, m); this->step(T, ind_dist(rng), g, rng); @@ -454,4 +461,3 @@ public: } } }; - |