diff options
| -rw-r--r-- | ising.cpp | 197 | ||||
| -rw-r--r-- | space_wolff.hpp | 696 |
2 files changed, 443 insertions, 450 deletions
@@ -2,7 +2,7 @@ #include "space_wolff.hpp" std::function<double(spin<signed, 2, signed>)> B_sin(unsigned L, unsigned n, double H) { - return [n, H, L] (spin<signed, 2, signed> s) -> double { + return [n, H, L](spin<signed, 2, signed> s) -> double { return H * s.s * cos(2 * M_PI * n * s.x(0) / ((double)L)); }; } @@ -23,98 +23,99 @@ int main(int argc, char* argv[]) { while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:n:")) != -1) { switch (opt) { - case 'N': - N = (unsigned)atof(optarg); - break; - case 'L': - L = atoi(optarg); - break; - case 'T': - T = atof(optarg); - break; - case 'H': - H = atof(optarg); - break; - case 'e': - ε = atof(optarg); - break; - case 'm': - mod = atoi(optarg); - break; - case 'M': - multi = atoi(optarg); - break; - case 'n': - mag = atof(optarg); - break; - default: - exit(1); + case 'N': + N = (unsigned)atof(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'T': + T = atof(optarg); + break; + case 'H': + H = atof(optarg); + break; + case 'e': + ε = atof(optarg); + break; + case 'm': + mod = atoi(optarg); + break; + case 'M': + multi = atoi(optarg); + break; + case 'n': + mag = atof(optarg); + break; + default: + exit(1); } } 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 { - 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) { - old_one_one = true; - } else if (old_diff(i) == 1 && old_one_one) { - old_many_ones = true; - } 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; - } + 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 { + 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) { + old_one_one = true; + } else if (old_diff(i) == 1 && old_one_one) { + old_many_ones = true; + } 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) { + 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; + } else if (were_nearest_neighbors) { + if (s1.s * s2.s == 1) { + return pZ; + } else { return 0.0; - } 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; - } + } + } else if (are_nearest_neighbors) { + if (s1_new.s * s2.s == -1) { + return pZ; } else { return 0.0; } - }; + } else { + return 0.0; + } + }; - std::function<double(spin<signed, D, signed>)> B_face = - [L, H] (spin<signed, D, signed> s) -> double { - return H * s.s * smiley[s.x(0) * 16 / L][s.x(1) * 16 / L]; - }; + std::function<double(spin<signed, D, signed>)> B_face = [L, + H](spin<signed, D, signed> s) -> double { + return H * s.s * smiley[s.x(0) * 16 / L][s.x(1) * 16 / L]; + }; std::function<double(spin<signed, D, signed>)> B; @@ -124,22 +125,23 @@ int main(int argc, char* argv[]) { B = B_face; } - std::function<std::set<unsigned>(model<signed, D, signed>&, unsigned, spin<signed, D, signed>)> neighbors = - [] (model<signed, D, signed>& m, unsigned i0, spin<signed, D, signed> s1) -> std::set<unsigned> { - std::set<unsigned> nn; - if (i0 < m.s.size()) { - std::set<unsigned> nn0 = m.dict.neighbors(m.s[i0].x, 1); - std::set<unsigned> nn1 = m.dict.neighbors(s1.x, 1); - nn.insert(nn0.begin(), nn0.end()); - nn.insert(nn1.begin(), nn1.end()); - nn.insert(m.s.size()); - } else { - for (unsigned i = 0; i < m.s.size(); i++) { - nn.insert(i); - } + std::function<std::set<unsigned>(model<signed, D, signed>&, unsigned, spin<signed, D, signed>)> + neighbors = [](model<signed, D, signed>& m, unsigned i0, + spin<signed, D, signed> s1) -> std::set<unsigned> { + std::set<unsigned> nn; + if (i0 < m.s.size()) { + std::set<unsigned> nn0 = m.dict.neighbors(m.s[i0].x, 1); + std::set<unsigned> nn1 = m.dict.neighbors(s1.x, 1); + nn.insert(nn0.begin(), nn0.end()); + nn.insert(nn1.begin(), nn1.end()); + nn.insert(m.s.size()); + } else { + for (unsigned i = 0; i < m.s.size(); i++) { + nn.insert(i); } - return nn; - }; + } + return nn; + }; model<signed, D, signed> ising(L, Z, B, neighbors); @@ -195,7 +197,8 @@ int main(int argc, char* argv[]) { std::array<double, 2> act = ising.Eq.τ(); std::vector<double> ρ = ising.Eq.ρ(); - outfile << L << " " << T << " " << mod << " " << H << " " << ising.Eq.num_added() << " " << ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1]; + outfile << L << " " << T << " " << mod << " " << H << " " << ising.Eq.num_added() << " " + << ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1]; for (double ρi : ρ) { outfile << " " << ρi; } diff --git a/space_wolff.hpp b/space_wolff.hpp index c87697e..7603ec2 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -1,44 +1,42 @@ -#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) { +#define AUTOLAG 1000 +#define AUTOWAIT 10000 + +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 +49,417 @@ 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); - euclidean pnew(this->L, tnew, rnew); + return pnew; + } - return pnew; - } -}; + euclidean inverse() const { + vector<T, D> tnew = -r.transpose() * t; + matrix<T, D> rnew = r.transpose(); -template <int D> -class dictionary { - private: - unsigned L; - std::vector<std::set<unsigned>> d; + euclidean pnew(this->L, tnew, rnew); - public: - dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {}; + return pnew; + } +}; - template <class T> - unsigned dictionary_index(vector<T, D> x) const { - unsigned pos_ind = 0; +template <int D> class dictionary { +private: + unsigned L; + std::vector<std::set<unsigned>> d; - for (unsigned i = 0; i < D; i++) { - pos_ind += pow(L, i) * (unsigned)x(i); - }; +public: + dictionary(unsigned Li) : L(Li), d(pow(Li, D)){}; - return pos_ind; - } + template <class T> unsigned dictionary_index(vector<T, D> x) const { + unsigned pos_ind = 0; - template <class T> - void record(vector<T, D> x, unsigned ind) { - d[dictionary_index<T>(x)].insert(ind); + for (unsigned i = 0; i < D; i++) { + pos_ind += pow(L, i) * (unsigned)x(i); }; - template <class T> - void remove(vector<T, D> x, unsigned ind) { - d[dictionary_index<T>(x)].erase(ind); - }; + return pos_ind; + } - template <class T> - std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const { - return nearest_neighbors_of(dictionary_index<T>(x), depth, {}); - }; + template <class T> void record(vector<T, D> x, unsigned ind) { + d[dictionary_index<T>(x)].insert(ind); + }; - std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, std::list<unsigned> ignores) const { - std::set<unsigned> ns = d[ind]; + template <class T> void remove(vector<T, D> x, unsigned ind) { + d[dictionary_index<T>(x)].erase(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)); + 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> nns = nearest_neighbors_of(ni, depth - 1, ignores_new); + std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, + std::list<unsigned> ignores) const { + std::set<unsigned> ns = d[ind]; - for (unsigned guy : nns) { - ns.insert(guy); - } + 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); + + 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<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>, spin<U, D, state>, spin<U, D, state>)> pZ; + std::function<double(spin<U, D, state>)> B; + std::function<double(spin<U, D, state>, spin<U, D, state>)> pB; + std::vector<matrix<U, D>> mats; + std::vector<vector<U, D>> steps; + double T; + 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, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z, + std::function<double(spin<U, D, state>, spin<U, D, state>, spin<U, D, state>)> pZ, 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); - - one_sequences(sequences, D); - - 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; - } - } + std::function<double(spin<U, D, state>), double(spin<U, D, state>)> pB, + std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns, + double T) + : L(L), s0(L), dict(L), neighbors(ns), Z(Z), pZ(pZ), B(B), pB(pB), T(T), + Eq(AUTOLAG, AUTOWAIT), Cq(AUTOLAG, AUTOWAIT) { + 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); + + 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; } + } + } - mats.push_back(m); + 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; - } - } - - 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(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 = pB(s0s_old, s0s_new); + } else { + p = pZ(s[i], s[j], si_new); } - } - 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); } + } }; |