summaryrefslogtreecommitdiff
path: root/ising.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'ising.cpp')
-rw-r--r--ising.cpp250
1 files changed, 151 insertions, 99 deletions
diff --git a/ising.cpp b/ising.cpp
index 1987a8f..38dce3b 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -1,14 +1,53 @@
+#include "smiley.hpp"
#include "space_wolff.hpp"
+#include "torus_symmetries.hpp"
+#include <GL/glut.h>
+
+const unsigned D = 2;
+typedef Model<signed, D, TorusGroup<signed, D>, signed> model;
+
+class animation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
+private:
+ uint64_t t1;
+ uint64_t t2;
+ unsigned n;
+ unsigned tmp;
+
+public:
+ animation(double L, unsigned w, int argc, char* argv[]) {
+ t1 = 0;
+ t2 = 0;
+ n = 0;
+
+ glutInit(&argc, argv);
+ glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+ glutInitWindowSize(w, w);
+ glutCreateWindow("wolff");
+ glClearColor(0.0, 0.0, 0.0, 0.0);
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+ gluOrtho2D(0, L, 0, L);
+ }
-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 H * s.s * cos(2 * M_PI * n * s.x(0) / ((double)L));
- };
-}
+ void post_cluster(const model& m) override {
+ glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
+ glClear(GL_COLOR_BUFFER_BIT);
+ glClearColor(0.0, 1.0, 0.0, 0.0);
+ for (const Spin<signed, 2, signed>* s : m.s) {
+ if (s->s == 1) {
+ glColor3f(0.0, 0.0, 0.0);
+ } else if (s->s == -1) {
+ glColor3f(1.0, 1.0, 1.0);
+ }
+ Vector<signed, 2> xx = m.s0.inverse().act(*s).x;
+ glRecti(xx(0), xx(1), xx(0) + 1, xx(1) + 1);
+ }
+ glFlush();
+ }
+};
int main(int argc, char* argv[]) {
- const unsigned D = 2;
unsigned L = 32;
unsigned N = 1000;
@@ -23,120 +62,133 @@ 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>)> Z =
- [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double {
- bool old_one_one = false;
- bool old_many_ones = false;
- bool old_any_two = false;
-
- vector<signed, D> old_diff = diff<signed, D>(L, s1.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;
- }
+ std::function<double(const Spin<signed, D, signed>&, const Spin<signed, D, signed>&)> Z =
+ [L, pZ](const Spin<signed, D, signed>& s1, const Spin<signed, D, signed>& s2) -> double {
+ bool old_one_one = false;
+ bool old_many_ones = false;
+ bool old_any_two = false;
+
+ Vector<signed, D> old_diff = diff<signed, D>(L, s1.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;
}
+ }
- bool were_on_someone = !old_one_one && !old_any_two;
- bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;
-
- if (were_on_someone) {
- return -std::numeric_limits<double>::infinity();;
- } else if (were_nearest_neighbors) {
- return s1.s * s2.s;
- } 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;
+ bool were_on_someone = !old_one_one && !old_any_two;
+ bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;
- if (mod > 0) {
- B = B_sin(L, mod, H);
- } else {
- B = B_face;
- }
+ if (were_on_someone) {
+ return -std::numeric_limits<double>::infinity();
+ ;
+ } else if (were_nearest_neighbors) {
+ return s1.s * s2.s;
+ } else {
+ return 0.0;
+ }
+ };
- 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;
- };
+ std::function<double(const Spin<signed, D, signed>&)> B_face =
+ [L, H](const Spin<signed, D, signed>& s) -> double {
+ return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L];
+ };
- model<signed, D, signed> ising(L, L, Z, B, neighbors);
+ std::function<double(Spin<signed, D, signed>)> B = B_face;
- randutils::auto_seed_128 seeds;
- std::mt19937 rng{seeds};
+ Model<signed, D, TorusGroup<signed, D>, signed> ising(L, Z, B);
- std::uniform_int_distribution<unsigned> coin(0, 1);
+ Rng rng;
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) * mag) || down >= pow(L, 2) * mag) {
- ising.s.push_back({{i, j}, 1});
+ Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
+ ss->x = {i, j};
+ if (rng.uniform<double>(0, 1) < mag) {
+ ss->s = 1;
up++;
} else {
- ising.s.push_back({{i, j}, -1});
+ ss->s = -1;
down++;
}
- ising.dict.record({i, j}, n);
+ ising.s.push_back(ss);
+ ising.dict.insert(ss);
n++;
}
}
+
+ std::vector<Vector<signed, 2>> torusVectors = torus_vecs<signed, 2>(L);
+ std::vector<Matrix<signed, 2>> torusMatrices = torus_mats<signed, 2>();
+
+ Gen<signed, D, TorusGroup<signed, D>, signed> g =
+ [L, &torusVectors,
+ &torusMatrices](Model<signed, D, TorusGroup<signed, D>, signed>& M,
+ Rng& r) -> Transformation<signed, D, TorusGroup<signed, D>, signed>* {
+ Matrix<signed, 2> m;
+ Vector<signed, 2> t;
+
+ m = r.pick(torusMatrices);
+ t(0) = r.uniform<signed>(0, L);
+ t(1) = r.uniform<signed>(0, L);
+ t = t - m * t;
+
+ TorusGroup<signed, 2> g = TorusGroup<signed, 2>({(signed)L, t, m});
+
+ Spin<signed, 2, signed>* ss = r.pick(M.s);
+ Spin<signed, 2, signed> s2 = g.act(*ss);
+
+ while (ss->x(0) == s2.x(0) && ss->x(1) == s2.x(1)) {
+ ss = r.pick(M.s);
+ s2 = g.act(*ss);
+ }
+
+ std::set<Spin<signed, 2, signed>*> s2s = M.dict.at(s2.x);
+
+ return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
+ };
+
+ animation A(L, 750, argc, argv);
+
+ ising.wolff(T, {g}, A, N);
+
/*
for (unsigned i = 0; i < L; i++) {
for (unsigned j = 0; j < L; j++) {
@@ -151,8 +203,7 @@ int main(int argc, char* argv[]) {
}
*/
- ising.update_energy();
-
+ /*
while (true) {
ising.wolff(T, N, rng);
std::array<double, 2> τ = ising.Eq.τ();
@@ -168,8 +219,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];
- for (double ρi : ρ) {
+ outfile << L << " " << T << " " << mod << " " << H << " " << ising.Eq.num_added() << " " <<
+ ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1]; for (double ρi : ρ) {
outfile << " " << ρi;
}
outfile << "\n";
@@ -195,6 +246,7 @@ int main(int argc, char* argv[]) {
}
snapfile.close();
+ */
return 0;
}