summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ising.cpp280
-rw-r--r--ising.hpp106
-rw-r--r--ising_animate.cpp105
-rw-r--r--ising_autocorrelation.cpp115
-rw-r--r--octree.hpp2
-rw-r--r--quantity.hpp47
6 files changed, 348 insertions, 307 deletions
diff --git a/ising.cpp b/ising.cpp
deleted file mode 100644
index 838ec9e..0000000
--- a/ising.cpp
+++ /dev/null
@@ -1,280 +0,0 @@
-
-#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;
- bool color;
-
-public:
- animation(double L, unsigned w, bool tcolor, int argc, char* argv[]) {
- color = tcolor;
- 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);
- }
-
- void post_cluster(const model& m) override {
- glClearColor(1.0, 1.0, 1.0, 1.0);
- glClear(GL_COLOR_BUFFER_BIT);
- for (const Spin<signed, 2, signed>* s : m.s) {
- if (s->s == 1) {
- if (color)
- glColor3f(1.0, 0.0, 0.0);
- else
- glColor3f(1.0, 1.0, 1.0);
- } else if (s->s == -1) {
- if (color)
- glColor3f(0.0, 0.0, 1.0);
- else
- glColor3f(0.0, 0.0, 0.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[]) {
-
- unsigned L = 32;
- unsigned N = 1000;
- unsigned mod = 0;
- unsigned multi = 1e4;
- double mag = 0.5;
- double pop = 1.0;
- double T = 2.0 / log(1.0 + sqrt(2.0));
- double H = 1.0;
- double ε = 0.1;
- bool color = false;
-
- int opt;
-
- while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:r:p:c")) != -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 'r':
- mag = atof(optarg);
- break;
- case 'p':
- pop = atof(optarg);
- break;
- case 'c':
- color = true;
- break;
- default:
- exit(1);
- }
- }
-
- double pZ = 1.0 - exp(-1.0 / T);
-
- 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;
-
- if (mod == 0) {
- B = [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];
- };
- } else {
- B = [mod, L, H](const Spin<signed, D, signed>& s) -> double {
- return H * s.s * sin(s.x(0) * mod * 2 * M_PI / L);
- };
- }
-
-
- Model<signed, D, TorusGroup<signed, D>, signed> ising(L, Z, B);
-
- 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 (rng.uniform<double>(0, 1) < pop) {
- 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 {
- ss->s = -1;
- down++;
- }
- 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);
-
- if (s2s.empty()) {
- return new SpinFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss);
- } else {
- return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
- }
- };
-
- animation A(L, 750, color, argc, argv);
-
- ising.wolff(T, {g}, A, 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++;
- }
- }
- */
-
- /*
- while (true) {
- ising.wolff(T, N, rng);
- std::array<double, 2> τ = ising.Eq.τ();
- std::cout << ising.Eq.num_added() << " " << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
- if (τ[1] / τ[0] < ε && τ[0] * multi < ising.Eq.num_added()) {
- break;
- }
- }
-
- std::ofstream outfile;
- outfile.open("out.dat", std::ios::app);
-
- 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 << " " << ρi;
- }
- outfile << "\n";
-
- std::ofstream snapfile;
- snapfile.open("snap.dat");
-
- std::vector<std::vector<signed>> snap(L);
- for (std::vector<signed>& line : snap) {
- line.resize(L);
- }
-
- for (spin<signed, D, signed> s : ising.s) {
- spin<signed, D, signed> snew = ising.s0.inverse().act(s);
- snap[snew.x(0)][snew.x(1)] = snew.s;
- }
-
- for (std::vector<signed> row : snap) {
- for (signed s : row) {
- snapfile << s << " ";
- }
- snapfile << "\n";
- }
-
- snapfile.close();
- */
-
- return 0;
-}
-
diff --git a/ising.hpp b/ising.hpp
new file mode 100644
index 0000000..6f140be
--- /dev/null
+++ b/ising.hpp
@@ -0,0 +1,106 @@
+
+#include "smiley.hpp"
+#include "space_wolff.hpp"
+#include "torus_symmetries.hpp"
+
+const unsigned D = 2;
+typedef Model<signed, D, TorusGroup<signed, D>, signed> isingModel;
+typedef Spin<signed, D, signed> isingSpin;
+
+std::function<double(const isingSpin&, const isingSpin&)> isingZ(unsigned L) {
+ return [L](const isingSpin& s1, const isingSpin& 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(isingSpin)> isingBFace(unsigned L, double H) {
+ return [L, H](const isingSpin& s) -> double {
+ return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L];
+ };
+}
+
+std::function<double(isingSpin)> isingBMod(unsigned L, unsigned mod, double H) {
+ return [mod, L, H](const isingSpin& s) -> double {
+ return H * s.s * sin(s.x(0) * mod * 2 * M_PI / L);
+ };
+}
+
+Gen<signed, D, TorusGroup<signed, D>, signed> isingGen(unsigned L) {
+ std::vector<Vector<signed, 2>> torusVectors = torus_vecs<signed, 2>(L);
+ std::vector<Matrix<signed, 2>> torusMatrices = torus_mats<signed, 2>();
+ return [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);
+
+ if (s2s.empty()) {
+ return new SpinFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss);
+ } else {
+ return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
+ }
+ };
+}
+
+void isingPopulate(isingModel& m, unsigned L, double pop, double mag) {
+ Rng rng;
+
+ for (unsigned i = 0; i < L; i++) {
+ for (unsigned j = 0; j < L; j++) {
+ if (rng.uniform<double>(0, 1) < pop) {
+ Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
+ ss->x = {i, j};
+ if (rng.uniform<double>(0, 1) < mag) {
+ ss->s = 1;
+ } else {
+ ss->s = -1;
+ }
+ m.s.push_back(ss);
+ m.dict.insert(ss);
+ }
+ }
+ }
+}
diff --git a/ising_animate.cpp b/ising_animate.cpp
new file mode 100644
index 0000000..e6586d9
--- /dev/null
+++ b/ising_animate.cpp
@@ -0,0 +1,105 @@
+#include "ising.hpp"
+#include <GL/glut.h>
+
+class Animation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
+private:
+ bool color;
+
+public:
+ Animation(double L, unsigned w, bool tcolor, int argc, char* argv[]) {
+ color = tcolor;
+
+ 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);
+ }
+
+ void post_cluster(const isingModel& m) override {
+ glClearColor(1.0, 1.0, 1.0, 1.0);
+ glClear(GL_COLOR_BUFFER_BIT);
+ for (const Spin<signed, 2, signed>* s : m.s) {
+ if (s->s == 1) {
+ if (color)
+ glColor3f(1.0, 0.0, 0.0);
+ else
+ glColor3f(1.0, 1.0, 1.0);
+ } else if (s->s == -1) {
+ if (color)
+ glColor3f(0.0, 0.0, 1.0);
+ else
+ glColor3f(0.0, 0.0, 0.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[]) {
+ unsigned L = 32;
+ unsigned N = 1000;
+ unsigned mod = 0;
+ double mag = 0.5;
+ double pop = 1.0;
+ double T = 2.0 / log(1.0 + sqrt(2.0));
+ double H = 1.0;
+ bool color = false;
+
+ int opt;
+
+ while ((opt = getopt(argc, argv, "N:L:T:H:m:r:p:c")) != -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 'm':
+ mod = atoi(optarg);
+ break;
+ case 'r':
+ mag = atof(optarg);
+ break;
+ case 'p':
+ pop = atof(optarg);
+ break;
+ case 'c':
+ color = true;
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ std::function<double(Spin<signed, D, signed>)> B;
+
+ if (mod == 0) {
+ B = isingBFace(L, H);
+ } else {
+ B = isingBMod(L, mod, H);
+ }
+
+ isingModel ising(L, isingZ(L), B);
+ isingPopulate(ising, L, pop, mag);
+
+ auto g = isingGen(L);
+
+ Animation A(L, 750, color, argc, argv);
+
+ ising.wolff(T, {g}, A, N);
+
+ return 0;
+}
diff --git a/ising_autocorrelation.cpp b/ising_autocorrelation.cpp
new file mode 100644
index 0000000..053be01
--- /dev/null
+++ b/ising_autocorrelation.cpp
@@ -0,0 +1,115 @@
+#include <iostream>
+#include <fstream>
+
+#include "ising.hpp"
+#include "quantity.hpp"
+
+class Autocorrelation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
+public:
+ Quantity E;
+ Autocorrelation(unsigned lag) : E(lag) {}
+
+ void post_cluster(const isingModel& m) override {
+ double E_tmp = 0;
+
+ TorusGroup<signed, D> s0Inv = m.s0.inverse();
+
+ for (const isingSpin* si : m.s) {
+ E_tmp -= m.B(s0Inv.act(*si));
+
+ for (const isingSpin* sj : m.dict.neighbors(si->x)) {
+ if (si != sj) {
+ E_tmp -= m.Z(*si, *sj) / 2;
+ }
+ }
+ }
+
+ E.add(E_tmp);
+ }
+};
+
+int main(int argc, char* argv[]) {
+ unsigned L = 32;
+ unsigned N = 1000;
+ unsigned mod = 0;
+ unsigned multi = 1e4;
+ double mag = 0.5;
+ double pop = 1.0;
+ double T = 2.0 / log(1.0 + sqrt(2.0));
+ double H = 1.0;
+ double ε = 0.1;
+ unsigned lag = 1e2;
+
+ int opt;
+
+ while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:r:p:l:")) != -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 'r':
+ mag = atof(optarg);
+ break;
+ case 'p':
+ pop = atof(optarg);
+ break;
+ case 'l':
+ lag = (unsigned)atof(optarg);
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ isingModel ising(L, isingZ(L), isingBMod(L, mod, H));
+ isingPopulate(ising, L, pop, mag);
+
+ auto g = isingGen(L);
+
+ measurement<signed, D, TorusGroup<signed, D>, signed> A_empty;
+
+ ising.wolff(T, {g}, A_empty, N);
+
+ Autocorrelation A(lag);
+
+ while (true) {
+ ising.wolff(T, {g}, A, N);
+ std::array<double, 2> τ = A.E.τ();
+ std::cout << A.E.num_added() << " " << A.E.avg() << " " << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
+ if (τ[1] / τ[0] < ε && τ[0] * multi < A.E.num_added()) {
+ break;
+ }
+ }
+
+ std::ofstream outfile;
+ outfile.open("out.dat", std::ios::app);
+
+ std::array<double, 2> act = A.E.τ();
+ std::vector<double> ρ = A.E.ρ();
+
+ outfile << L << " " << T << " " << mod << " " << H << " " << A.E.num_added() << " " <<
+ A.E.avg() << " " << A.E.serr() << " " << act[0] << " " << act[1]; for (double ρi : ρ) {
+ outfile << " " << ρi;
+ }
+ outfile << "\n";
+
+ return 0;
+}
diff --git a/octree.hpp b/octree.hpp
index bb13242..d8f0aa3 100644
--- a/octree.hpp
+++ b/octree.hpp
@@ -5,8 +5,6 @@
#include <set>
#include <unordered_map>
-#include <iostream>
-
#include "spin.hpp"
#include "vector.hpp"
diff --git a/quantity.hpp b/quantity.hpp
index 55d5a28..d2a7dea 100644
--- a/quantity.hpp
+++ b/quantity.hpp
@@ -1,57 +1,54 @@
#pragma once
-#include <vector>
-#include <list>
-#include <cmath>
#include <array>
+#include <cmath>
+#include <list>
+#include <vector>
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) {
+ Quantity(unsigned lag) : C(lag) {
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);
+ 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++;
}
- n++;
}
- double avg() const { return total / (n - wait); }
+ double avg() const { return total / n; }
- double avg2() const { return total2 / (n - wait); }
+ double avg2() const { return total2 / n; }
std::vector<double> ρ() const {
- double C0 = C.front() / (n - wait);
- double avg2 = pow(total / (n - wait), 2);
+ double C0 = C.front() / n;
+ double avg2 = pow(total / n, 2);
std::vector<double> ρtmp;
for (double Ct : C) {
- ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
+ ρtmp.push_back((Ct / n - avg2) / (C0 - avg2));
}
return ρtmp;
@@ -69,15 +66,15 @@ public:
M++;
}
- return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+ return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
}
double σ() const {
- return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+ return 2.0 / n * this->τ()[0] * (C[0] / n - pow(this->avg(), 2));
}
double serr() const { return sqrt(this->σ()); }
- unsigned num_added() const { return n - wait; }
+ unsigned num_added() const { return n; }
};