diff options
| -rw-r--r-- | ising.cpp | 280 | ||||
| -rw-r--r-- | ising.hpp | 106 | ||||
| -rw-r--r-- | ising_animate.cpp | 105 | ||||
| -rw-r--r-- | ising_autocorrelation.cpp | 115 | ||||
| -rw-r--r-- | octree.hpp | 2 | ||||
| -rw-r--r-- | quantity.hpp | 47 |
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..cd65b6b --- /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("wolffWindow"); + 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; +} @@ -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; } }; |