diff options
| -rw-r--r-- | ising.cpp | 141 | ||||
| -rw-r--r-- | space_wolff.hpp (renamed from space_wolff.cpp) | 113 | ||||
| -rw-r--r-- | spheres.cpp | 104 |
3 files changed, 245 insertions, 113 deletions
diff --git a/ising.cpp b/ising.cpp new file mode 100644 index 0000000..6939749 --- /dev/null +++ b/ising.cpp @@ -0,0 +1,141 @@ + +#include "space_wolff.hpp" + +int main(int argc, char* argv[]) { + const unsigned D = 2; + + unsigned L = 32; + unsigned N = 1000; + double T = 2.0 / log(1.0 + sqrt(2.0)); + double H = 1.0; + + int opt; + + while ((opt = getopt(argc, argv, "N:L:T:H:")) != -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; + default: + exit(1); + } + } + + std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z = + [] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double { + bool one_one = false; + bool many_ones = false; + bool any_two = false; + + for (unsigned i = 0; i < D; i++) { + unsigned diff = abs(s1.x(i) - s2.x(i)); + if (diff == 1 && !one_one) { + one_one = true; + } else if (diff == 1 && one_one) { + many_ones = true; + break; + } else if (diff > 1) { + any_two = true; + break; + } + } + + if (!one_one && !any_two) { + return -std::numeric_limits<double>::infinity(); + } else if (one_one && !many_ones && !any_two) { + return s1.s * s2.s; + } else { + return 0; + } + }; + + std::function<double(spin<signed, D, signed>)> B = + [L, H] (spin<signed, D, signed> s) -> double { + return H * s.s * smiley[s.x(1) * 16 / L][s.x(0) * 16 / L]; + }; + + 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> os1 = m.dict.on_site(s1.x); + std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x); + std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x); + nn.insert(nn0.begin(), nn0.end()); + nn.insert(nn1.begin(), nn1.end()); + nn.insert(os1.begin(), os1.end()); + nn.insert(m.s.size()); + } else { + for (unsigned i = 0; i < m.s.size(); i++) { + nn.insert(i); + } + } + return nn; + }; + + model<signed, D, signed> ising(L, Z, B, neighbors); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + std::uniform_int_distribution<unsigned> coin(0, 1); + + unsigned n = 0; + unsigned up = 0; + unsigned down = 0; + for (unsigned i = 0; i < L; i++) { + for (unsigned j = 0; j < L; j++) { + if ((coin(rng) && up < pow(L, 2) / 2) || down >= pow(L, 2) / 2) { + ising.s.push_back({{i, j}, 1}); + up++; + } else { + ising.s.push_back({{i, j}, -1}); + down++; + } + ising.dict.record<signed>({i, j}, n); + n++; + } + } + /* + for (unsigned i = 0; i < L; i++) { + for (unsigned j = 0; j < L; j++) { + if (i < L / 2) { + ising.s.push_back({{i, j}, 1}); + } else { + ising.s.push_back({{i, j}, -1}); + } + ising.dict.record<signed>({i, j}, n); + n++; + } + } + */ + + ising.wolff(T, N, rng); + + std::vector<signed> output(pow(L, D)); + + for (spin<signed, D, signed> s : ising.s) { + spin<signed, D, signed> rs = ising.s0.inverse().act(s); + output[L * rs.x(1) + rs.x(0)] = s.s; + } + + for (unsigned i = 0; i < L; i++) { + for (unsigned j = 0; j < L; j++) { + unsigned out = output[L * i + j] == 1 ? 1 : 0; + std::cout << out; + } + std::cout << "\n"; + } + + return 0; +} + diff --git a/space_wolff.cpp b/space_wolff.hpp index b88944e..aaef673 100644 --- a/space_wolff.cpp +++ b/space_wolff.hpp @@ -331,116 +331,3 @@ class model { } }; -int main(int argc, char* argv[]) { - unsigned L = 32; - const unsigned D = 2; - - std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z = - [] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double { - bool one_one = false; - bool many_ones = false; - bool any_two = false; - - for (unsigned i = 0; i < D; i++) { - unsigned diff = abs(s1.x(i) - s2.x(i)); - if (diff == 1 && !one_one) { - one_one = true; - } else if (diff == 1 && one_one) { - many_ones = true; - break; - } else if (diff > 1) { - any_two = true; - break; - } - } - - if (!one_one && !any_two) { - return -std::numeric_limits<double>::infinity(); - } else if (one_one && !many_ones && !any_two) { - return s1.s * s2.s; - } else { - return 0; - } - }; - - std::function<double(spin<signed, D, signed>)> B = - [] (spin<signed, D, signed> s) -> double { - return 100 * s.s * smiley[s.x(0) / 2][s.x(1) / 2]; - }; - - 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> os1 = m.dict.on_site(s1.x); - std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x); - std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x); - nn.insert(nn0.begin(), nn0.end()); - nn.insert(nn1.begin(), nn1.end()); - nn.insert(os1.begin(), os1.end()); - nn.insert(m.s.size()); - } else { - for (unsigned i = 0; i < m.s.size(); i++) { - nn.insert(i); - } - } - return nn; - }; - - model<signed, D, signed> ising(L, Z, B, neighbors); - - randutils::auto_seed_128 seeds; - std::mt19937 rng{seeds}; - - std::uniform_int_distribution<unsigned> coin(0, 1); - - unsigned n = 0; - unsigned up = 0; - unsigned down = 0; - for (unsigned i = 0; i < L; i++) { - for (unsigned j = 0; j < L; j++) { - if ((coin(rng) && up < pow(L, 2) / 2) || down >= pow(L, 2) / 2) { - ising.s.push_back({{i, j}, 1}); - up++; - } else { - ising.s.push_back({{i, j}, -1}); - down++; - } - ising.dict.record<signed>({i, j}, n); - n++; - } - } - /* - for (unsigned i = 0; i < L; i++) { - for (unsigned j = 0; j < L; j++) { - if (i < L / 2) { - ising.s.push_back({{i, j}, 1}); - } else { - ising.s.push_back({{i, j}, -1}); - } - ising.dict.record<signed>({i, j}, n); - n++; - } - } - */ - - ising.wolff(2.0 / log(1.0 + sqrt(2.0)), 5000, rng); - - std::vector<signed> output(pow(L, D)); - - for (spin<signed, D, signed> s : ising.s) { - spin<signed, D, signed> rs = ising.s0.inverse().act(s); - output[L * rs.x(1) + rs.x(0)] = s.s; - } - - for (unsigned i = 0; i < L; i++) { - for (unsigned j = 0; j < L; j++) { - unsigned out = output[L * i + j] == 1 ? 1 : 0; - std::cout << out; - } - std::cout << "\n"; - } - - return 0; -} - diff --git a/spheres.cpp b/spheres.cpp new file mode 100644 index 0000000..106ce4e --- /dev/null +++ b/spheres.cpp @@ -0,0 +1,104 @@ + +#include "space_wolff.hpp" + +int main(int argc, char* argv[]) { + const unsigned D = 2; + + double L = 32; + unsigned N = 1000; + double T = 2.0 / log(1.0 + sqrt(2.0)); + double H = 1.0; + unsigned n = 25; + + int opt; + + while ((opt = getopt(argc, argv, "n:N:L:T:H:")) != -1) { + switch (opt) { + case 'n': + n = (unsigned)atof(optarg); + break; + case 'N': + N = (unsigned)atof(optarg); + break; + case 'L': + L = atof(optarg); + break; + case 'T': + T = atof(optarg); + break; + case 'H': + H = atof(optarg); + break; + default: + exit(1); + } + } + + std::function<double(spin<double, D, double>, spin<double, D, double>)> Z = + [L] (spin<double, D, double> s1, spin<double, D, double> s2) -> double { + vector<double, D> diff = s1.x - s2.x; + for (unsigned i = 0; i < D; i++) { + if (fabs(diff(i)) > L / 2) { + diff(i) = L - fabs(diff(i)); + } else { + diff(i) = fabs(diff(i)); + } + } + if (diff.transpose() * diff < pow(s1.s + s2.s, 2)) { + return -std::numeric_limits<double>::infinity(); + } else { + return 0; + } + }; + + std::function<double(spin<double, D, double>)> B = + [L, H] (spin<double, D, double> s) -> double { + return H * sin(2 * M_PI * 3 * s.x(0) / L); + }; + + std::function<std::set<unsigned>(model<double, D, double>&, unsigned, spin<double, D, double>)> neighbors = + [] (model<double, D, double>& m, unsigned i0, spin<double, D, double> s1) -> std::set<unsigned> { + std::set<unsigned> nn; + if (i0 < m.s.size()) { + std::set<unsigned> os1 = m.dict.on_site(s1.x); + std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x); + std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x); + std::set<unsigned> nnn0 = m.dict.next_nearest_neighbors(m.s[i0].x); + std::set<unsigned> nnn1 = m.dict.next_nearest_neighbors(s1.x); + nn.insert(nn0.begin(), nn0.end()); + nn.insert(nn1.begin(), nn1.end()); + nn.insert(nnn0.begin(), nnn0.end()); + nn.insert(nnn1.begin(), nnn1.end()); + nn.insert(os1.begin(), os1.end()); + nn.insert(m.s.size()); + } else { + for (unsigned i = 0; i < m.s.size(); i++) { + nn.insert(i); + } + } + return nn; + }; + + model<double, D, double> sphere(L, Z, B, neighbors); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + std::uniform_real_distribution<double> dist(0.0, L); + + for (unsigned i = 0; i < n; i++) { + vector<double, D> pos = {dist(rng), dist(rng)}; + sphere.s.push_back({pos, 0.5}); + sphere.dict.record<double>(pos, i); + } + + sphere.wolff(T, N, rng); + + for (spin<double, D, double> s : sphere.s) { + spin<double, D, double> rs = sphere.s0.inverse().act(s); + std::cout << s.x.transpose() << "\n"; + } + + return 0; +} + |