summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2019-09-03 16:50:52 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2019-09-03 16:50:52 -0400
commit243c5b3e6f3761d51acc4097fa2ed2b870290956 (patch)
tree1e49d280f173f7f6fd6a15a337b148bcf866769a
parent1586d6b2247c510c11d9a4847bde67cb0947243d (diff)
downloadspace_wolff-243c5b3e6f3761d51acc4097fa2ed2b870290956.tar.gz
space_wolff-243c5b3e6f3761d51acc4097fa2ed2b870290956.tar.bz2
space_wolff-243c5b3e6f3761d51acc4097fa2ed2b870290956.zip
partial completion of more modular energy and bond probabiliesefficiency-work
-rw-r--r--ising.cpp197
-rw-r--r--space_wolff.hpp696
2 files changed, 443 insertions, 450 deletions
diff --git a/ising.cpp b/ising.cpp
index fffce19..8f39798 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -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);
}
+ }
};