summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitmodules3
m---------randutils0
-rw-r--r--space_wolff.cpp446
3 files changed, 449 insertions, 0 deletions
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..0a97926
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "randutils"]
+ path = randutils
+ url = https://gist.github.com/imneme/540829265469e673d045/
diff --git a/randutils b/randutils
new file mode 160000
+Subproject 8486a610a954a8248c12485fb4cfc390a5f5f85
diff --git a/space_wolff.cpp b/space_wolff.cpp
new file mode 100644
index 0000000..b88944e
--- /dev/null
+++ b/space_wolff.cpp
@@ -0,0 +1,446 @@
+
+#include <vector>
+#include <list>
+#include <set>
+#include <iostream>
+#include <functional>
+#include <random>
+#include <queue>
+#include <eigen3/Eigen/Dense>
+#include "randutils/randutils.hpp"
+
+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, class state>
+class spin {
+ public:
+ vector<T, D> x;
+ state s;
+};
+
+template <class T, int D>
+class euclidean {
+ private:
+ unsigned L;
+
+ public:
+ vector<T, D> t;
+ matrix<T, D> r;
+ euclidean(unsigned 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(unsigned L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
+ t = t0;
+ r = r0;
+ }
+
+ template <class state>
+ spin<T, D, state> act(spin<T, D, state> s) {
+ spin<T, D, state> s_new;
+
+ 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;
+ }
+
+ euclidean act(const euclidean& x) {
+ vector<T, D> tnew = r * x.t + t;
+ matrix<T, D> rnew = r * x.r;
+
+ euclidean pnew(this->L, tnew, rnew);
+
+ return pnew;
+ }
+
+ euclidean inverse() {
+ vector<T, D> tnew = - r.transpose() * t;
+ matrix<T, D> rnew = r.transpose();
+
+ euclidean pnew(this->L, tnew, rnew);
+
+ return pnew;
+ }
+};
+
+template <int D>
+class dictionary {
+ private:
+ unsigned L;
+ std::vector<std::set<unsigned>> d;
+
+ public:
+ dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {};
+
+ template <class T>
+ unsigned dictionary_index(vector<T, D> x) {
+ unsigned pos_ind = 0;
+
+ for (unsigned i = 0; i < D; i++) {
+ pos_ind += pow(L, i) * (unsigned)x(i);
+ };
+
+ return pos_ind;
+ }
+
+ template <class T>
+ void record(vector<T, D> x, unsigned ind) {
+ d[dictionary_index<T>(x)].insert(ind);
+ };
+
+ template <class T>
+ void remove(vector<T, D> x, unsigned ind) {
+ d[dictionary_index<T>(x)].erase(ind);
+ };
+
+ template <class T>
+ std::set<unsigned> on_site(vector<T, D> x) {
+ return d[dictionary_index<T>(x)];
+ };
+
+ template <class T>
+ std::set<unsigned> nearest_neighbors(vector<T, D> x) {
+ unsigned ind = dictionary_index<T>(x);
+ std::set<unsigned> ns;
+
+ for (unsigned i = 0; i < D; 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));
+ for (unsigned nii : d[ni]) {
+ ns.insert(nii);
+ }
+ }
+ }
+
+ return ns;
+ };
+
+ template <class T>
+ std::set<unsigned> next_nearest_neighbors(vector<T, D> x) {
+ unsigned ind = dictionary_index<T>(x);
+ std::set<unsigned> ns;
+
+ for (unsigned i = 0; i < D; 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));
+ for (unsigned k = 0; k < D; k++) {
+ if (k != i) {
+ for (signed l : {-1, 1}) {
+ unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1));
+ for (unsigned nnii : d[nni]) {
+ ns.insert(nnii);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return ns;
+ };
+
+ template <class T>
+ std::set<unsigned> next_next_nearest_neighbors(vector<T, D> x) {
+ unsigned ind = dictionary_index<T>(x);
+ std::set<unsigned> ns;
+
+ for (unsigned i = 0; i < D; 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));
+ for (unsigned k = 0; k < D; k++) {
+ if (k != i) {
+ for (signed l : {-1, 1}) {
+ unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1));
+ for (unsigned m = 0; m < D; m++) {
+ if (m != i && m != k) {
+ for (signed n : {-1, 1}) {
+ unsigned nnni = pow(L, m + 1) * (nni / ((unsigned)pow(L, m + 1))) + fmod(pow(L, m + 1) + nni + n * pow(L, m), pow(L, m + 1));
+ for (unsigned nnnii : d[nnni]) {
+ ns.insert(nnnii);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return ns;
+ };
+};
+
+template <class U, int D, class state>
+class model {
+ public:
+ unsigned 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>)> B;
+ double E;
+
+ model(unsigned L, 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) {
+ }
+
+ void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) {
+ std::uniform_real_distribution<double> dist(0.0, 1.0);
+
+ std::queue<unsigned> queue;
+ queue.push(ind);
+
+ std::vector<bool> visited(s.size() + 1, false);
+
+ while (!queue.empty()) {
+ unsigned i = queue.front();
+ queue.pop();
+
+ if (!visited[i]) {
+ visited[i] = true;
+
+ bool we_are_ghost = i == s.size();
+
+ spin<U, D, state> si_new;
+ euclidean<U, D> s0_new(L);
+
+ if (we_are_ghost) {
+ s0_new = r.act(s0);
+ } else {
+ si_new = r.act(s[i]);
+ }
+
+ for (unsigned j : neighbors(*this, i, si_new)) {
+ if (j != i) {
+ double dE;
+ bool neighbor_is_ghost = j == s.size();
+
+ if (we_are_ghost || neighbor_is_ghost) {
+ spin<U, D, state> s0s_old, s0s_new;
+ unsigned non_ghost;
+
+ 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]);
+ }
+
+ dE = B(s0s_old) - B(s0s_new);
+ } else {
+ dE = Z(s[i], s[j]) - Z(si_new, s[j]);
+ }
+
+ double p = 1.0 - exp(-dE / T);
+
+ if (dist(rng) < p) {
+ queue.push(j);
+ }
+ }
+ }
+
+ if (we_are_ghost) {
+ s0 = s0_new;
+ } else {
+ dict.remove(s[i].x, i);
+ s[i] = si_new;
+ dict.record(s[i].x, i);
+ }
+ }
+ }
+ }
+
+ 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);
+
+ for (unsigned i = 0; i < N; i++) {
+ vector<U, D> t;
+ matrix<U, D> 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 {
+ m(j, k) = 1;
+ }
+ } else if ((j == k && j != flip_D1) && j != flip_D2) {
+ m(j, k) = 1;
+ } else {
+ m(j, k) = 0;
+ }
+ }
+ }
+
+ euclidean<U, D> g(L, t, m);
+
+ this->step(T, ind_dist(rng), g, rng);
+ }
+ }
+};
+
+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;
+}
+