From 272627e1e0936e3adfd168e4e6ce31cc8954d719 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 31 Jul 2019 16:48:45 -0400 Subject: added another qualification for finishing to minimze bias --- ising.cpp | 37 ++++++++++++-------- space_wolff.hpp | 107 +++++++++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 113 insertions(+), 31 deletions(-) diff --git a/ising.cpp b/ising.cpp index 645c759..e4ab65c 100644 --- a/ising.cpp +++ b/ising.cpp @@ -1,18 +1,25 @@ #include "space_wolff.hpp" +std::function)> B_sin(unsigned L, unsigned n, double H) { + return [n, H, L] (spin s) -> double { + return H * cos(2 * M_PI * n * s.x[0] / ((double)L)); + }; +} + int main(int argc, char* argv[]) { const unsigned D = 2; unsigned L = 32; unsigned N = 1000; + unsigned mod = 0; double T = 2.0 / log(1.0 + sqrt(2.0)); double H = 1.0; double ε = 0.1; int opt; - while ((opt = getopt(argc, argv, "N:L:T:H:e:")) != -1) { + while ((opt = getopt(argc, argv, "N:L:T:H:e:m:")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); @@ -29,6 +36,9 @@ int main(int argc, char* argv[]) { case 'e': ε = atof(optarg); break; + case 'm': + mod = atoi(optarg); + break; default: exit(1); } @@ -62,11 +72,19 @@ int main(int argc, char* argv[]) { } }; - std::function)> B = + std::function)> B_face = [L, H] (spin s) -> double { return H * s.s * smiley[s.x(1) * 16 / L][s.x(0) * 16 / L]; }; + std::function)> B; + + if (mod > 0) { + B = B_sin(L, mod, H); + } else { + B = B_face; + } + std::function(model&, unsigned, spin)> neighbors = [] (model& m, unsigned i0, spin s1) -> std::set { std::set nn; @@ -127,7 +145,7 @@ int main(int argc, char* argv[]) { ising.wolff(T, N, rng); std::array τ = ising.Eq.τ(); std::cout << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n"; - if (τ[1] / τ[0] < ε) { + if (τ[1] / τ[0] < ε && τ[0] * 1e4 < ising.Eq.num_added()) { break; } } @@ -140,20 +158,11 @@ int main(int argc, char* argv[]) { } std::ofstream outfile; - outfile.open("snap.dat"); - - for (unsigned i = 0; i < L; i++) { - for (unsigned j = 0; j < L; j++) { - unsigned out = output[L * i + j] == 1 ? 1 : 0; - outfile << out << " "; - } - outfile << "\n"; - } - outfile.close(); + outfile.open("out.dat", std::ios::app); std::array act = ising.Eq.τ(); - std::cout << act[0] << " " << act[1] << "\n"; + outfile << L << " " << T << " " << mod << " " << H << " " << ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1] << "\n"; return 0; } diff --git a/space_wolff.hpp b/space_wolff.hpp index b60e274..9a4b4b4 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -237,6 +237,10 @@ class quantity { double serr() const { return sqrt(this->σ()); } + + unsigned num_added() const { + return n - wait; + } }; template @@ -249,16 +253,87 @@ class model { std::function(model&, unsigned, spin)> neighbors; std::function, spin)> Z; std::function)> B; + std::vector> mats; + std::vector> steps; long double E; quantity Eq; quantity Cq; + void one_sequences(std::list>& sequences, unsigned level) { + if (level > 0) { + unsigned new_level = level - 1; + unsigned old_length = sequences.size(); + for (std::array& sequence : sequences) { + std::array new_sequence = sequence; + new_sequence[new_level] = -1; + sequences.push_front(new_sequence); + } + one_sequences(sequences, new_level); + } + } + model(unsigned L, std::function, spin)> Z, std::function)> B, std::function(model&, unsigned, spin)> ns) : L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) { + std::array ini_sequence; + ini_sequence.fill(1); + std::list> sequences; + sequences.push_back(ini_sequence); + + one_sequences(sequences, D); + + sequences.pop_front(); // don't want the identity matrix! + + for (std::array sequence : sequences) { + matrix 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 v; + for (unsigned i = 0; i < D; i++) { + if (sequence[i] == 1) { + v(i) = 0; + } else { + v(i) = L / 2; + } + } + + steps.push_back(v); } + for (unsigned i = 0; i < D; i++) { + for (unsigned j = 0; j < D; j++) { + if (i != j) { + matrix 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); + } + } + } + } + void update_energy() { E = 0; for (unsigned i = 0; i < s.size(); i++) { @@ -347,31 +422,29 @@ class model { std::uniform_real_distribution t_dist(0, L); std::uniform_int_distribution r_dist(0, D - 1); std::uniform_int_distribution ind_dist(0, s.size() - 1); + std::uniform_int_distribution coin(0, mats.size() + steps.size() - 1); for (unsigned i = 0; i < N; i++) { vector t; matrix m; - for (unsigned j = 0; j < D; j++) { - t(j) = (U)t_dist(rng); - } - - unsigned flip_D1 = r_dist(rng); - unsigned flip_D2 = r_dist(rng); - - for (unsigned j = 0; j < D; j++) { - for (unsigned k = 0; k < D; k++) { - if ((j == flip_D1 && k == flip_D2) || (j == flip_D2 && k == flip_D1)) { - if (flip_D1 <= flip_D2) { - m(j, k) = -1; - } else { + 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; } - } else if ((j == k && j != flip_D1) && j != flip_D2) { - m(j, k) = 1; - } else { - m(j, k) = 0; } } + + t = steps[flip - mats.size()]; } euclidean g(L, t, m); -- cgit v1.2.3-70-g09d2