diff options
| -rw-r--r-- | ising.cpp | 39 | ||||
| -rw-r--r-- | space_wolff.hpp | 681 |
2 files changed, 337 insertions, 383 deletions
@@ -54,17 +54,13 @@ int main(int argc, char* argv[]) { double pZ = 1.0 - exp(-1.0 / T); - std::function<double(spin<signed, D, signed>, spin<signed, D, signed>, spin<signed, D, signed>)> Z = - [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2, spin<signed, D, signed> s1_new) -> double { + std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z = + [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double { bool old_one_one = false; bool old_many_ones = false; bool old_any_two = false; - bool new_one_one = false; - bool new_many_ones = false; - bool new_any_two = false; vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x); - vector<signed, D> new_diff = diff<signed, D>(L, s1_new.x, s2.x); for (unsigned i = 0; i < D; i++) { if (old_diff(i) == 1 && !old_one_one) { @@ -74,38 +70,15 @@ int main(int argc, char* argv[]) { } else if (old_diff(i) > 1) { old_any_two = true; } - if (new_diff(i) == 1 && !new_one_one) { - new_one_one = true; - } else if (new_diff(i) == 1 && new_one_one) { - new_many_ones = true; - } else if (new_diff(i) > 1) { - new_any_two = true; - } } bool were_on_someone = !old_one_one && !old_any_two; - bool are_on_someone = !new_one_one && !new_any_two; bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two; - bool are_nearest_neighbors = new_one_one && !new_many_ones && !new_any_two; if (were_on_someone) { - return 0.0; - } else if (are_on_someone) { - return 1.0; - } else if (were_nearest_neighbors && are_nearest_neighbors) { - return 0.0; + return -std::numeric_limits<double>::infinity();; } else if (were_nearest_neighbors) { - if (s1.s * s2.s == 1) { - return pZ; - } else { - return 0.0; - } - } else if (are_nearest_neighbors) { - if (s1_new.s * s2.s == -1) { - return pZ; - } else { - return 0.0; - } + return s1.s * s2.s; } else { return 0.0; } @@ -141,7 +114,7 @@ int main(int argc, char* argv[]) { return nn; }; - model<signed, D, signed> ising(L, Z, B, neighbors); + model<signed, D, signed> ising(L, L, Z, B, neighbors); randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; @@ -160,7 +133,7 @@ int main(int argc, char* argv[]) { ising.s.push_back({{i, j}, -1}); down++; } - ising.dict.record<signed>({i, j}, n); + ising.dict.record({i, j}, n); n++; } } diff --git a/space_wolff.hpp b/space_wolff.hpp index c87697e..6ac9788 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -1,44 +1,39 @@ -#include <vector> -#include <list> -#include <set> -#include <iostream> +#include "randutils/randutils.hpp" +#include <cmath> +#include <eigen3/Eigen/Dense> #include <fstream> #include <functional> -#include <random> +#include <iostream> +#include <list> #include <queue> -#include <cmath> -#include <eigen3/Eigen/Dense> -#include "randutils/randutils.hpp" +#include <random> +#include <set> +#include <vector> -const std::array<std::array<unsigned, 16>, 16> smiley = - {{ - {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}, - {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}}, - {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}}, - {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}}, - {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}}, - {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}}, - {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}}, - {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}, - {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}, - {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}}, - {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}}, - {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}}, - {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}}, - {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}}, - {{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 T, int D> -using vector = Eigen::Matrix<T, D, 1>; - -template <class T, int D> -using matrix = Eigen::Matrix<T, D, D>; - -template <class T, int D> -vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) { +const std::array<std::array<unsigned, 16>, 16> smiley = { + {{{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}, + {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}}, + {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}}, + {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}}, + {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}}, + {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}}, + {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}}, + {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}, + {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}, + {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}}, + {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}}, + {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}}, + {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}}, + {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}}, + {{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 T, int D> using vector = Eigen::Matrix<T, D, 1>; + +template <class T, int D> using matrix = Eigen::Matrix<T, D, D>; + +template <class T, int D> vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) { vector<T, D> v; for (unsigned i = 0; i < D; i++) { @@ -51,425 +46,411 @@ vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) { return v; } -template <class T, int D, class state> -class spin { - public: - vector<T, D> x; - state s; +template <class T, int D, class state> class spin { +public: + vector<T, D> x; + state s; }; -template <class T, int D> -class euclidean { - private: - T L; - - public: - vector<T, D> t; - matrix<T, D> r; - euclidean(T L) : L(L) { - for (unsigned i = 0; i < D; i++) { - t(i) = 0; - r(i, i) = 1; - for (unsigned j = 1; j < D; j++) { - r(i, (i + j) % D) = 0; - } +template <class T, int D> class euclidean { +private: + T L; + +public: + vector<T, D> t; + matrix<T, D> r; + euclidean(T L) : L(L) { + for (unsigned i = 0; i < D; i++) { + t(i) = 0; + r(i, i) = 1; + for (unsigned j = 1; j < D; j++) { + r(i, (i + j) % D) = 0; } } + } - euclidean(T L, vector<T, D> t0, matrix<T, D> r0) : L(L) { - t = t0; - r = r0; - } + euclidean(T L, vector<T, D> t0, matrix<T, D> r0) : L(L) { + t = t0; + r = r0; + } - template <class state> - spin<T, D, state> act(const spin<T, D, state>& s) const { - spin<T, D, state> s_new; + template <class state> spin<T, D, state> act(const spin<T, D, state>& s) const { + spin<T, D, state> s_new; - s_new.x = t + r * s.x; - s_new.s = s.s; + s_new.x = t + r * s.x; + s_new.s = s.s; - for (unsigned i = 0; i < D; i++) { - s_new.x(i) = fmod(L + s_new.x(i), L); - } - - return s_new; + for (unsigned i = 0; i < D; i++) { + s_new.x(i) = fmod(L + s_new.x(i), L); } - euclidean act(const euclidean& x) const { - vector<T, D> tnew = r * x.t + t; - matrix<T, D> rnew = r * x.r; - - for (unsigned i = 0; i < D; i++) { - tnew(i) = fmod(L + tnew(i), L); - } + return s_new; + } - euclidean pnew(this->L, tnew, rnew); + euclidean act(const euclidean& x) const { + vector<T, D> tnew = r * x.t + t; + matrix<T, D> rnew = r * x.r; - return pnew; + for (unsigned i = 0; i < D; i++) { + tnew(i) = fmod(L + tnew(i), L); } - euclidean inverse() const { - vector<T, D> tnew = - r.transpose() * t; - matrix<T, D> rnew = r.transpose(); + euclidean pnew(this->L, tnew, rnew); + + return pnew; + } - euclidean pnew(this->L, tnew, rnew); + euclidean inverse() const { + vector<T, D> tnew = -r.transpose() * t; + matrix<T, D> rnew = r.transpose(); - return pnew; - } + euclidean pnew(this->L, tnew, rnew); + + return pnew; + } }; -template <int D> -class dictionary { - private: - unsigned L; - std::vector<std::set<unsigned>> d; +template <class T, int D> class dictionary { +private: + unsigned N; + T L; + std::vector<std::set<unsigned>> d; - public: - dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {}; +public: + dictionary(unsigned Ni, double Li) : N(Ni), L(Li), d(pow(Ni, D)) {}; - template <class T> - unsigned dictionary_index(vector<T, D> x) const { - unsigned pos_ind = 0; + unsigned dictionary_index(vector<T, D> x) const { + unsigned pos_ind = 0; - for (unsigned i = 0; i < D; i++) { - pos_ind += pow(L, i) * (unsigned)x(i); - }; + for (unsigned i = 0; i < D; i++) { + pos_ind += pow(N, i) * (unsigned)std::floor(x(i) * N / L); + }; - return pos_ind; - } + return pos_ind; + } - template <class T> - void record(vector<T, D> x, unsigned ind) { - d[dictionary_index<T>(x)].insert(ind); - }; + void record(vector<T, D> x, unsigned ind) { + d[this->dictionary_index(x)].insert(ind); + }; - template <class T> - void remove(vector<T, D> x, unsigned ind) { - d[dictionary_index<T>(x)].erase(ind); - }; + void remove(vector<T, D> x, unsigned ind) { + d[this->dictionary_index(x)].erase(ind); + }; - template <class T> - std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const { - return nearest_neighbors_of(dictionary_index<T>(x), depth, {}); - }; + std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const { + return nearest_neighbors_of(this->dictionary_index(x), depth, {}); + }; - std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, std::list<unsigned> ignores) const { - std::set<unsigned> ns = d[ind]; + std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, + std::list<unsigned> ignores) const { + std::set<unsigned> ns = d[ind]; - if (depth > 0) { - for (unsigned i = 0; i < D; i++) { - if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k){return i==k;})) { - std::list<unsigned> ignores_new = ignores; - ignores_new.push_back(i); - for (signed j : {-1, 1}) { - unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1)); + if (depth > 0) { + for (unsigned i = 0; i < D; i++) { + if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k) { return i == k; })) { + std::list<unsigned> ignores_new = ignores; + ignores_new.push_back(i); + for (signed j : {-1, 1}) { + unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1)); - std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new); + std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new); - for (unsigned guy : nns) { - ns.insert(guy); - } + for (unsigned guy : nns) { + ns.insert(guy); } } } } + } - return ns; - }; + return ns; + }; }; class quantity { - private: - double total; - double total2; - std::vector<double> C; - unsigned wait; - - public: - unsigned n; - std::list<double> hist; - - quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) { - n = 0; - total = 0; - total2 = 0; - } +private: + double total; + double total2; + std::vector<double> C; + unsigned wait; + +public: + unsigned n; + std::list<double> hist; + + quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) { + n = 0; + total = 0; + total2 = 0; + } - void add(double x) { - if (n > wait) { - hist.push_front(x); - if (hist.size() > C.size()) { - hist.pop_back(); - unsigned t = 0; - for (double a : hist) { - C[t] += a * x; - t++; - } - total += x; - total2 += pow(x, 2); + void add(double x) { + if (n > wait) { + hist.push_front(x); + if (hist.size() > C.size()) { + hist.pop_back(); + unsigned t = 0; + for (double a : hist) { + C[t] += a * x; + t++; } + total += x; + total2 += pow(x, 2); } - n++; - } - - double avg() const { - return total / (n - wait); } + n++; + } - double avg2() const { - return total2 / (n - wait); - } + double avg() const { return total / (n - wait); } - std::vector<double> ρ() const { - double C0 = C.front() / (n - wait); - double avg2 = pow(total / (n - wait), 2); + double avg2() const { return total2 / (n - wait); } - std::vector<double> ρtmp; + std::vector<double> ρ() const { + double C0 = C.front() / (n - wait); + double avg2 = pow(total / (n - wait), 2); - for (double Ct : C) { - ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); - } + std::vector<double> ρtmp; - return ρtmp; + for (double Ct : C) { + ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); } - std::array<double, 2> τ() const { - double τtmp = 0.5; - unsigned M = 1; - double c = 8.0; + return ρtmp; + } - std::vector<double> ρ_tmp = this->ρ(); + std::array<double, 2> τ() const { + double τtmp = 0.5; + unsigned M = 1; + double c = 8.0; - while (c * τtmp > M && M < C.size()) { - τtmp += ρ_tmp[M]; - M++; - } + std::vector<double> ρ_tmp = this->ρ(); - return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)}; + while (c * τtmp > M && M < C.size()) { + τtmp += ρ_tmp[M]; + M++; } - double σ() const { - return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2)); - } + return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)}; + } - double serr() const { - return sqrt(this->σ()); - } + double σ() const { + return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2)); + } - unsigned num_added() const { - return n - wait; - } + double serr() const { return sqrt(this->σ()); } + + unsigned num_added() const { return n - wait; } }; -template <class U, int D, class state> -class model { - public: - U L; - euclidean<U, D> s0; - std::vector<spin<U, D, state>> s; - dictionary<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>, 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; - - void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) { - if (level > 0) { - unsigned new_level = level - 1; - unsigned old_length = sequences.size(); - for (std::array<double, D>& sequence : sequences) { - std::array<double, D> new_sequence = sequence; - new_sequence[new_level] = -1; - sequences.push_front(new_sequence); - } - one_sequences(sequences, new_level); +template <class U, int 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; + long double E; + quantity Eq; + quantity Cq; + + void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) { + if (level > 0) { + unsigned new_level = level - 1; + unsigned old_length = sequences.size(); + for (std::array<double, D>& sequence : sequences) { + std::array<double, D> new_sequence = sequence; + new_sequence[new_level] = -1; + sequences.push_front(new_sequence); } + one_sequences(sequences, new_level); } + } - model(U L, std::function<double(spin<U, D, state>, spin<U, D, state>, spin<U, D, state>)> Z, + 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(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; - sequences.push_back(ini_sequence); + 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; + sequences.push_back(ini_sequence); - one_sequences(sequences, D); + one_sequences(sequences, D); - sequences.pop_front(); // don't want the identity matrix! + sequences.pop_front(); // don't want the identity matrix! - for (std::array<double, D> sequence : sequences) { - matrix<U, D> m; - for (unsigned i = 0; i < D; i++) { - for (unsigned j = 0; j < D; j++) { - if (i == j) { - m(i, j) = sequence[i]; - } else { - m(i, j) = 0; - } - } + for (std::array<double, D> sequence : sequences) { + matrix<U, D> m; + for (unsigned i = 0; i < D; i++) { + for (unsigned j = 0; j < D; j++) { + if (i == j) { + m(i, j) = sequence[i]; + } else { + m(i, j) = 0; } + } + } - mats.push_back(m); - - vector<U, D> v; - for (unsigned i = 0; i < D; i++) { - if (sequence[i] == 1) { - v(i) = 0; - } else { - v(i) = L / 2; - } - } + mats.push_back(m); - steps.push_back(v); + vector<U, D> v; + for (unsigned i = 0; i < D; i++) { + if (sequence[i] == 1) { + v(i) = 0; + } else { + v(i) = L / 2; } + } - for (unsigned i = 0; i < D; i++) { - for (unsigned j = 0; j < D; j++) { - if (i != j) { - 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)) { - if (i < j) { - m(k, l) = 1; - } else { - m(k, l) = -1; - } - } else { - m(k, l) = 0; - } + steps.push_back(v); + } + + for (unsigned i = 0; i < D; i++) { + for (unsigned j = 0; j < D; j++) { + if (i != j) { + 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)) { + if (i < j) { + m(k, l) = 1; + } else { + m(k, l) = -1; } + } else { + m(k, l) = 0; } - mats.push_back(m); } } + mats.push_back(m); } } + } + } - void update_energy() { - E = 0; - for (unsigned i = 0; i < s.size(); i++) { - for (unsigned j = 0; j < i; j++) { - // E -= Z(s[i], s[j]); - } - - E -= B(s0.inverse().act(s[i])); + void update_energy() { + E = 0; + for (unsigned i = 0; i < s.size(); i++) { + for (unsigned j = 0; j < i; j++) { + E -= Z(s[i], s[j]); } - } - 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); + E -= B(s0.inverse().act(s[i])); + } + } - std::queue<unsigned> queue; - queue.push(ind); + 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); - std::vector<bool> visited(s.size() + 1, false); + std::queue<unsigned> queue; + queue.push(ind); - while (!queue.empty()) { - unsigned i = queue.front(); - queue.pop(); + std::vector<bool> visited(s.size() + 1, false); - if (!visited[i]) { - visited[i] = true; + while (!queue.empty()) { + unsigned i = queue.front(); + queue.pop(); - bool we_are_ghost = i == s.size(); + if (!visited[i]) { + visited[i] = true; - spin<U, D, state> si_new; - euclidean<U, D> s0_new(L); + bool we_are_ghost = i == s.size(); - if (we_are_ghost) { - s0_new = r.act(s0); - } else { - si_new = r.act(s[i]); - } + spin<U, D, state> si_new; + euclidean<U, D> s0_new(L); - for (unsigned j : neighbors(*this, i, si_new)) { - if (j != i) { - double p; - bool neighbor_is_ghost = j == s.size(); + if (we_are_ghost) { + s0_new = r.act(s0); + } else { + si_new = r.act(s[i]); + } - if (we_are_ghost || neighbor_is_ghost) { - spin<U, D, state> s0s_old, s0s_new; - unsigned non_ghost; + for (unsigned j : neighbors(*this, i, si_new)) { + if (j != i) { + double p; + bool neighbor_is_ghost = j == s.size(); - if (neighbor_is_ghost) { - non_ghost = i; - s0s_old = s0.inverse().act(s[i]); - s0s_new = s0.inverse().act(si_new); - } else { - non_ghost = j; - s0s_old = s0.inverse().act(s[j]); - s0s_new = s0_new.inverse().act(s[j]); - } + if (we_are_ghost || neighbor_is_ghost) { + spin<U, D, state> s0s_old, s0s_new; + unsigned non_ghost; - p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T); + if (neighbor_is_ghost) { + non_ghost = i; + s0s_old = s0.inverse().act(s[i]); + s0s_new = s0.inverse().act(si_new); } else { - p = Z(s[i], s[j], si_new); + non_ghost = j; + s0s_old = s0.inverse().act(s[j]); + s0s_new = s0_new.inverse().act(s[j]); } - if (dist(rng) < p) { - queue.push(j); - } + p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T); + } else { + p = 1.0 - exp(-(Z(si_new, s[j]) - Z(s[i], s[j])) / T); } - } - if (we_are_ghost) { - s0 = s0_new; - } else { - dict.remove(s[i].x, i); - s[i] = si_new; - dict.record(s[i].x, i); - cluster_size++; + if (dist(rng) < p) { + queue.push(j); + } } } - } - Cq.add(cluster_size); + if (we_are_ghost) { + s0 = s0_new; + } else { + dict.remove(s[i].x, i); + s[i] = si_new; + dict.record(s[i].x, i); + cluster_size++; + } + } } - void wolff(double T, unsigned N, std::mt19937& rng) { - std::uniform_real_distribution<double> t_dist(0, L); - std::uniform_int_distribution<unsigned> r_dist(0, D - 1); - std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1); - 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; - unsigned flip = coin(rng); - if (flip < mats.size()) { - for (unsigned j = 0; j < D; j++) { - t(j) = (U)t_dist(rng); - } - m = mats[flip]; - } else { - for (unsigned j = 0; j < D; j++) { - for (unsigned k = 0; k < D; k++) { - if (j == k) { - m(j, k) = 1; - } else { - m(j, k) = 0; - } + Cq.add(cluster_size); + } + + void wolff(double T, unsigned N, std::mt19937& rng) { + std::uniform_real_distribution<double> t_dist(0, L); + std::uniform_int_distribution<unsigned> r_dist(0, D - 1); + std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1); + 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; + unsigned flip = coin(rng); + if (flip < mats.size()) { + for (unsigned j = 0; j < D; j++) { + t(j) = (U)t_dist(rng); + } + m = mats[flip]; + } else { + for (unsigned j = 0; j < D; j++) { + for (unsigned k = 0; k < D; k++) { + if (j == k) { + m(j, k) = 1; + } else { + m(j, k) = 0; } } - - t = steps[flip - mats.size()]; } - euclidean<U, D> g(L, t, m); + t = steps[flip - mats.size()]; + } - this->step(T, ind_dist(rng), g, rng); + euclidean<U, D> g(L, t, m); - this->update_energy(); - Eq.add(E); - } + this->step(T, ind_dist(rng), g, rng); + + this->update_energy(); + Eq.add(E); } + } }; |