From 53f05e5f0bc0b79b4422ecfbb3dde7e346fdddd4 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 15 Jan 2020 19:17:50 -0500 Subject: refactor --- euclidean.hpp | 116 ++++++++++++ matrix.hpp | 7 + octree.hpp | 87 +++++++++ quantity.hpp | 83 +++++++++ smiley.hpp | 23 +++ space_wolff.hpp | 484 ++------------------------------------------------- spheres.cpp | 230 ++++++++++++------------ spheres_infinite.cpp | 181 +++++++++---------- spin.hpp | 11 ++ torus_symmetries.hpp | 100 +++++++++++ vector.hpp | 20 +++ 11 files changed, 669 insertions(+), 673 deletions(-) create mode 100644 euclidean.hpp create mode 100644 matrix.hpp create mode 100644 octree.hpp create mode 100644 quantity.hpp create mode 100644 smiley.hpp create mode 100644 spin.hpp create mode 100644 torus_symmetries.hpp create mode 100644 vector.hpp diff --git a/euclidean.hpp b/euclidean.hpp new file mode 100644 index 0000000..b1a5e80 --- /dev/null +++ b/euclidean.hpp @@ -0,0 +1,116 @@ + +#pragma once + +#include "matrix.hpp" +#include "spin.hpp" +#include "vector.hpp" + +template class Euclidean { +public: + Vector t; + Matrix r; + + Euclidean(T 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(Vector t0, Matrix r0) { + t = t0; + r = r0; + } + + template Spin act(const Spin& s) const { + Spin s_new; + + s_new.x = t + r * s.x; + s_new.s = s.s; + + return s_new; + } + + Euclidean act(const Euclidean& x) const { + Vector tnew = r * x.t + t; + Matrix rnew = r * x.r; + + Euclidean pnew(tnew, rnew); + + return pnew; + } + + Euclidean inverse() const { + Vector tnew = -r.transpose() * t; + Matrix rnew = r.transpose(); + + Euclidean pnew(tnew, rnew); + + return pnew; + } +}; + +template class TorusGroup { +private: + T L; + +public: + Vector t; + Matrix r; + + /** brief TorusGroup - default constructor, constructs the identity + */ + TorusGroup(T 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; + } + } + } + + TorusGroup(T L, Vector t0, Matrix r0) : L(L) { + t = t0; + r = r0; + } + + template Spin act(const Spin& s) const { + Spin 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; + } + + TorusGroup act(const TorusGroup& x) const { + Vector tnew = r * x.t + t; + Matrix rnew = r * x.r; + + for (unsigned i = 0; i < D; i++) { + tnew(i) = fmod(L + tnew(i), L); + } + + TorusGroup pnew(this->L, tnew, rnew); + + return pnew; + } + + TorusGroup inverse() const { + Vector tnew = -r.transpose() * t; + Matrix rnew = r.transpose(); + + TorusGroup pnew(this->L, tnew, rnew); + + return pnew; + } +}; + diff --git a/matrix.hpp b/matrix.hpp new file mode 100644 index 0000000..96b63fb --- /dev/null +++ b/matrix.hpp @@ -0,0 +1,7 @@ + +#pragma once + +#include + +template using Matrix = Eigen::Matrix; + diff --git a/octree.hpp b/octree.hpp new file mode 100644 index 0000000..1be6e0a --- /dev/null +++ b/octree.hpp @@ -0,0 +1,87 @@ + +#pragma once + +#include +#include +#include + +#include "spin.hpp" +#include "vector.hpp" + +namespace std { +template struct hash> { + typedef array argument_type; + typedef size_t result_type; + + result_type operator()(const argument_type& a) const { + hash hasher; + result_type h = 0; + for (result_type i = 0; i < N; ++i) { + h = h * 31 + hasher(a[i]); + } + return h; + } +}; +} // namespace std + +template class Octree { +private: + unsigned N; + T L; + std::unordered_map, std::set*>> data; + +public: + Octree(T L, unsigned N) { + L = L; + N = N; + } + + std::array ind(Vector x) const { + std::array ind; + + for (unsigned i = 0; i < D; i++) { + ind[i] = (signed)std::floor(N * x(i) / L); + } + + return ind; + } + + void insert(Spin* s) { data[ind(s->x)].insert(s); }; + + void remove(Spin* s) { + data[ind(s->x)].erase(s); + if (data[ind(s->x)].empty()) { + data.erase(ind(s->x)); + } + }; + + std::set*> neighbors(const Vector& x) const { + std::array i0 = ind(x); + std::set*> ns; + + nearest_neighbors_of(i0, D + 1, ns); + + return ns; + }; + + void nearest_neighbors_of(std::array i0, unsigned depth, + std::set*>& ns) const { + if (depth == 0) { + auto it = data.find(i0); + if (it != data.end()) { + ns.insert(it->second.begin(), it->second.end()); + } + } else { + for (signed j : {-1, 0, 1}) { + std::array i1 = i0; + if (N < 2) { + i1[depth - 1] += j; + } else { + i1[depth - 1] = (N + i1[depth - 1] + j) % N; + } + nearest_neighbors_of(i1, depth - 1, ns); + } + } + }; +}; + diff --git a/quantity.hpp b/quantity.hpp new file mode 100644 index 0000000..55d5a28 --- /dev/null +++ b/quantity.hpp @@ -0,0 +1,83 @@ + +#pragma once + +#include +#include +#include +#include + +class Quantity { +private: + double total; + double total2; + std::vector C; + unsigned wait; + +public: + unsigned n; + std::list hist; + + Quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) { + 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); + } + } + n++; + } + + double avg() const { return total / (n - wait); } + + double avg2() const { return total2 / (n - wait); } + + std::vector ρ() const { + double C0 = C.front() / (n - wait); + double avg2 = pow(total / (n - wait), 2); + + std::vector ρtmp; + + for (double Ct : C) { + ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); + } + + return ρtmp; + } + + std::array τ() const { + double τtmp = 0.5; + unsigned M = 1; + double c = 8.0; + + std::vector ρ_tmp = this->ρ(); + + while (c * τtmp > M && M < C.size()) { + τtmp += ρ_tmp[M]; + M++; + } + + return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)}; + } + + double σ() const { + return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2)); + } + + double serr() const { return sqrt(this->σ()); } + + unsigned num_added() const { return n - wait; } +}; + diff --git a/smiley.hpp b/smiley.hpp new file mode 100644 index 0000000..bd58fc1 --- /dev/null +++ b/smiley.hpp @@ -0,0 +1,23 @@ + +#pragma once + +#include + +const std::array, 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}}}}; + diff --git a/space_wolff.hpp b/space_wolff.hpp index 7345726..196b3b5 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -1,488 +1,36 @@ -#include -#include -#include +#pragma once + +#include #include -#include #include #include -#include -#include #include -#include -#include -#include - -namespace std { -template struct hash> { - typedef array argument_type; - typedef size_t result_type; - - result_type operator()(const argument_type& a) const { - hash hasher; - result_type h = 0; - for (result_type i = 0; i < N; ++i) { - h = h * 31 + hasher(a[i]); - } - return h; - } -}; -} // namespace std #include "randutils/randutils.hpp" -const std::array, 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 using Vector = Eigen::Matrix; - -template using Matrix = Eigen::Matrix; - -/** brief diff - periodic subtraction of vectors - */ -template Vector diff(T L, Vector v1, Vector v2) { - Vector v; - - for (unsigned i = 0; i < D; i++) { - v(i) = std::abs(v1(i) - v2(i)); - if (v(i) > L / 2) { - v(i) = L - v(i); - } - } - - return v; -} - -template class Spin { -public: - Vector x; - S s; -}; - -template class Euclidean { -public: - Vector t; - Matrix r; - - Euclidean(T 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(Vector t0, Matrix r0) { - t = t0; - r = r0; - } - - template Spin act(const Spin& s) const { - Spin s_new; - - s_new.x = t + r * s.x; - s_new.s = s.s; - - return s_new; - } - - Euclidean act(const Euclidean& x) const { - Vector tnew = r * x.t + t; - Matrix rnew = r * x.r; - - Euclidean pnew(tnew, rnew); - - return pnew; - } - - Euclidean inverse() const { - Vector tnew = -r.transpose() * t; - Matrix rnew = r.transpose(); - - Euclidean pnew(tnew, rnew); - - return pnew; - } -}; - -template class TorusGroup { -private: - T L; - -public: - Vector t; - Matrix r; - - /** brief TorusGroup - default constructor, constructs the identity - */ - TorusGroup(T 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; - } - } - } - - TorusGroup(T L, Vector t0, Matrix r0) : L(L) { - t = t0; - r = r0; - } - - template Spin act(const Spin& s) const { - Spin 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; - } - - TorusGroup act(const TorusGroup& x) const { - Vector tnew = r * x.t + t; - Matrix rnew = r * x.r; - - for (unsigned i = 0; i < D; i++) { - tnew(i) = fmod(L + tnew(i), L); - } - - TorusGroup pnew(this->L, tnew, rnew); - - return pnew; - } - - TorusGroup inverse() const { - Vector tnew = -r.transpose() * t; - Matrix rnew = r.transpose(); - - TorusGroup pnew(this->L, tnew, rnew); - - return pnew; - } -}; - -template class Octree { -private: - unsigned N; - T L; - std::unordered_map, std::set*>> data; - -public: - Octree(T L, unsigned N) { - L = L; - N = N; - } - - std::array ind(Vector x) const { - std::array ind; - - for (unsigned i = 0; i < D; i++) { - ind[i] = (signed)std::floor(N * x(i) / L); - } - - return ind; - } - - void insert(Spin* s) { data[ind(s->x)].insert(s); }; - - void remove(Spin* s) { - data[ind(s->x)].erase(s); - if (data[ind(s->x)].empty()) { - data.erase(ind(s->x)); - } - }; - - std::set*> neighbors(const Vector& x) const { - std::array i0 = ind(x); - std::set*> ns; - - nearest_neighbors_of(i0, D + 1, ns); - - return ns; - }; - - void nearest_neighbors_of(std::array i0, unsigned depth, std::set*>& ns) const { - if (depth == 0) { - auto it = data.find(i0); - if (it != data.end()) { - ns.insert(it->second.begin(), it->second.end()); - } - } else { - for (signed j : {-1, 0, 1}) { - std::array i1 = i0; - if (N < 2) { - i1[depth - 1] += j; - } else { - i1[depth - 1] = (N + i1[depth - 1] + j) % N; - } - nearest_neighbors_of(i1, depth - 1, ns); - } - } - }; -}; - -/* -template class Octree { -private: - T L; - unsigned N; - std::vector*>> data; - -public: - Octree(T L, unsigned depth) : L(L), N(pow(2, depth)), data(pow(N, D)){}; - - unsigned ind(Vector x) const { - unsigned pos_ind = 0; - - for (unsigned i = 0; i < D; i++) { - pos_ind += pow(N, i) * (unsigned)std::floor(N * x(i) / L); - } - - assert(pos_ind < data.size()); - - return pos_ind; - } - - void insert(Spin* s) { data[ind(s->x)].insert(s); }; - - void remove(Spin* s) { data[ind(s->x)].erase(s); }; - - std::set*> neighbors(const Vector& x, unsigned depth) const { - std::set*> ns; - std::set seen; - nearest_neighbors_of(ind(x), depth, ns, seen); - return ns; - }; - - void nearest_neighbors_of(unsigned ind, unsigned depth, std::set*>& ns, - std::set& seen) const { - ns.insert(data[ind].begin(), data[ind].end()); - seen.insert(ind); - - if (depth > 0) { - for (unsigned i = 0; i < D; i++) { - for (signed j : {-1, 1}) { - unsigned nind = pow(N, i + 1) * (ind / ((unsigned)pow(N, i + 1))) + - fmod(pow(N, i + 1) + ind + j * pow(N, i), pow(N, i + 1)); - - if (!seen.contains(nind)) { - seen.insert(nind); - nearest_neighbors_of(nind, depth - 1, ns, seen); - } - } - } - } - } -}; -*/ - -class quantity { -private: - double total; - double total2; - std::vector C; - unsigned wait; - -public: - unsigned n; - std::list hist; - - quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) { - 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); - } - } - n++; - } - - double avg() const { return total / (n - wait); } - - double avg2() const { return total2 / (n - wait); } - - std::vector ρ() const { - double C0 = C.front() / (n - wait); - double avg2 = pow(total / (n - wait), 2); - - std::vector ρtmp; - - for (double Ct : C) { - ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); - } - - return ρtmp; - } - - std::array τ() const { - double τtmp = 0.5; - unsigned M = 1; - double c = 8.0; - - std::vector ρ_tmp = this->ρ(); - - while (c * τtmp > M && M < C.size()) { - τtmp += ρ_tmp[M]; - M++; - } - - return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)}; - } - - double σ() const { - return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2)); - } - - double serr() const { return sqrt(this->σ()); } - - unsigned num_added() const { return n - wait; } -}; +#include "euclidean.hpp" +#include "octree.hpp" +#include "quantity.hpp" +#include "spin.hpp" template class Model; template class measurement { public: virtual void pre_cluster(const Model&, unsigned, const R&){}; - virtual void plain_bond_visited(const Model&, const Spin*, const Spin*, - const Spin&, double){}; + virtual void plain_bond_visited(const Model&, const Spin*, + const Spin*, const Spin&, double){}; virtual void plain_site_transformed(const Model&, const Spin*, const Spin&){}; - virtual void ghost_bond_visited(const Model&, const Spin&, const Spin&, - double){}; + virtual void ghost_bond_visited(const Model&, const Spin&, + const Spin&, double){}; virtual void ghost_site_transformed(const Model&, const R&){}; virtual void post_cluster(const Model&){}; }; -template -void one_sequences(std::list>& sequences, unsigned level) { - if (level > 0) { - unsigned new_level = level - 1; - unsigned old_length = sequences.size(); - for (std::array& sequence : sequences) { - std::array new_sequence = sequence; - new_sequence[new_level] = -1; - sequences.push_front(new_sequence); - } - one_sequences(sequences, new_level); - } -} - -template -std::vector> torus_vecs(U L) { - std::vector> vecs; - std::array ini_sequence; - ini_sequence.fill(1); - std::list> sequences; - sequences.push_back(ini_sequence); - one_sequences(sequences, D); - sequences.pop_front(); // don't want the identity matrix! - - for (std::array sequence : sequences) { - Vector v; - for (unsigned i = 0; i < D; i++) { - if (sequence[i] == 1) { - v(i) = 0; - } else { - v(i) = L / 2; - } - } - - vecs.push_back(v); - } - - return vecs; -} - -template -std::vector> torus_mats() { - std::vector> mats; - - std::array ini_sequence; - ini_sequence.fill(1); - std::list> sequences; - sequences.push_back(ini_sequence); - - one_sequences(sequences, D); - - sequences.pop_front(); // don't want the identity matrix! - - for (std::array sequence : sequences) { - Matrix m; - for (unsigned i = 0; i < D; i++) { - for (unsigned j = 0; j < D; j++) { - if (i == j) { - m(i, j) = sequence[i]; - } else { - m(i, j) = 0; - } - } - } - - mats.push_back(m); - } - - for (unsigned i = 0; i < D; i++) { - for (unsigned j = 0; j < D; j++) { - if (i != j) { - Matrix m; - for (unsigned k = 0; k < D; k++) { - for (unsigned l = 0; l < D; l++) { - if ((k == i && l == j) || (k == j && l == i)) { - if (i < j) { - m(k, l) = 1; - } else { - m(k, l) = -1; - } - } else { - m(k, l) = 0; - } - } - } - mats.push_back(m); - } - } - } - - return mats; -} - template class Model { public: U L; @@ -493,14 +41,11 @@ public: std::function&)> B; randutils::mt19937_rng rng; - Model(U L, std::function&, const Spin&)> Z, std::function&)> B) : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {} void step(double T, unsigned ind, R& r, measurement& A) { - unsigned cluster_size = 0; - std::queue*> queue; queue.push(&(s[ind])); @@ -556,17 +101,15 @@ public: dict.remove(si); *si = si_new; dict.insert(si); - cluster_size++; } } } } - void resize(double T, double P) { - } + void resize(double T, double P) {} void wolff(double T, std::function&, randutils::mt19937_rng&)> g, - measurement& A, unsigned N) { + measurement& A, unsigned N) { for (unsigned i = 0; i < N; i++) { R r = g(*this, rng); unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1)); @@ -579,3 +122,4 @@ public: } } }; + diff --git a/spheres.cpp b/spheres.cpp index 97c1e18..b7a9d94 100644 --- a/spheres.cpp +++ b/spheres.cpp @@ -1,100 +1,105 @@ -#include "space_wolff.hpp" #include +#include +#include + +#include "space_wolff.hpp" +#include "torus_symmetries.hpp" const unsigned D = 2; typedef Model, double> model; class animation : public measurement, double> { - 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(-1, L + 1, -1 , L + 1); - } +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(-1, L + 1, -1, L + 1); + } - void pre_cluster(const model&, unsigned, const TorusGroup&) override { - tmp = 0; - } + void pre_cluster(const model&, unsigned, const TorusGroup&) override { tmp = 0; } - void plain_site_transformed(const model&, const Spin*, const Spin&) override { - tmp++; - } + void plain_site_transformed(const model&, const Spin*, + const Spin&) override { + tmp++; + } - void post_cluster(const model& m) override { - glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); - glClear(GL_COLOR_BUFFER_BIT); - for (const Spin& s : m.s) { - glBegin(GL_POLYGON); - unsigned n_points = 50; - glColor3f(0.0f, 0.0f, 0.0f); - for (unsigned i = 0; i < n_points; i++) { - glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points), m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points)); - } - glEnd(); + void post_cluster(const model& m) override { + glClearColor(1.0f, 1.0f, 1.0f, 1.0f); + glClear(GL_COLOR_BUFFER_BIT); + for (const Spin& s : m.s) { + glBegin(GL_POLYGON); + unsigned n_points = 50; + glColor3f(0.0f, 0.0f, 0.0f); + for (unsigned i = 0; i < n_points; i++) { + glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points), + m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points)); } - glFlush(); - - t1 += tmp; - t2 += tmp * tmp; - n++; + glEnd(); } + glFlush(); - void clear() { - t1 = 0; - t2 = 0; - n = 0; - } + t1 += tmp; + t2 += tmp * tmp; + n++; + } - double var() { - return (t2 - t1 * t1 / (double)n) / (double)n; - } + void clear() { + t1 = 0; + t2 = 0; + n = 0; + } + + double var() { return (t2 - t1 * t1 / (double)n) / (double)n; } }; -std::function(const model&, randutils::mt19937_rng&)> eGen(const std::vector>& mats, const std::vector>& vecs, double ε, double L) { - return [&mats, &vecs, L, ε] (const model& M, randutils::mt19937_rng& rng) -> TorusGroup { - Vector t; - Matrix m; - unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1)); - if (flip < mats.size()) { - unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size()); - t = M.s[f_ind].x; - for (unsigned j = 0; j < D; j++) { - t(j) = fmod(10 * L + t(j) + rng.variate(0.0, ε), L); - } - m = mats[flip]; - } else { - for (unsigned j = 0; j < D; j++) { - for (unsigned k = 0; k < D; k++) { - if (j == k) { - m(j, k) = 1; - } else { - m(j, k) = 0; +std::function(const model&, randutils::mt19937_rng&)> +eGen(const std::vector>& mats, const std::vector>& vecs, + double ε, double L) { + return + [&mats, &vecs, L, ε](const model& M, randutils::mt19937_rng& rng) -> TorusGroup { + Vector t; + Matrix m; + unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1)); + if (flip < mats.size()) { + unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size()); + t = M.s[f_ind].x; + for (unsigned j = 0; j < D; j++) { + t(j) = fmod(10 * L + t(j) + rng.variate(0.0, ε), L); + } + m = mats[flip]; + } else { + for (unsigned j = 0; j < D; j++) { + for (unsigned k = 0; k < D; k++) { + if (j == k) { + m(j, k) = 1; + } else { + m(j, k) = 0; + } + } } - } - } - - t = vecs[flip - mats.size()]; - } - TorusGroup g(M.L, t, m); - return g; - }; + t = vecs[flip - mats.size()]; + } + TorusGroup g(M.L, t, m); + return g; + }; } int main(int argc, char* argv[]) { @@ -110,23 +115,23 @@ int main(int argc, char* argv[]) { 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); + 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); } } @@ -134,25 +139,24 @@ int main(int argc, char* argv[]) { double a = 0.05; std::function&, const Spin&)> Z = - [L, a, k] (const Spin& s1, const Spin& s2) -> double { - Vector d = diff(L, s1.x, s2.x); - - double σ = s1.s + s2.s; - double δ = σ - sqrt(d.transpose() * d); - - if (δ > - a * σ) { - return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2)); - } else if (δ > - 2 * a * σ) { - return 0.5 * k * pow(δ + 2 * a * σ, 2); - } else { - return 0; - } - }; + [L, a, k](const Spin& s1, const Spin& s2) -> double { + Vector d = diff(L, s1.x, s2.x); + + double σ = s1.s + s2.s; + double δ = σ - sqrt(d.transpose() * d); - std::function)> B = - [L, H] (Spin s) -> double { - return H * s.x(1); - }; + if (δ > -a * σ) { + return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2)); + } else if (δ > -2 * a * σ) { + return 0.5 * k * pow(δ + 2 * a * σ, 2); + } else { + return 0; + } + }; + + std::function)> B = [L, H](Spin s) -> double { + return H * s.x(1); + }; std::vector> mats = torus_mats(); std::vector> vecs = torus_vecs(L); diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index f2f998a..dba7ab3 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -1,4 +1,7 @@ +#include +#include + #include "space_wolff.hpp" #include @@ -6,67 +9,66 @@ const unsigned D = 2; typedef Model, double> model; class animation : public measurement, double> { - 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(-L, L, -L , L); - } +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(-L, L, -L, L); + } - void pre_cluster(const model&, unsigned, const Euclidean&) override { - tmp = 0; - } + void pre_cluster(const model&, unsigned, const Euclidean&) override { tmp = 0; } - void plain_site_transformed(const model&, const Spin*, const Spin&) override { - tmp++; - } + void plain_site_transformed(const model&, const Spin*, + const Spin&) override { + tmp++; + } - void post_cluster(const model& m) override { - glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); - glClear(GL_COLOR_BUFFER_BIT); - for (const Spin& s : m.s) { - glBegin(GL_POLYGON); - unsigned n_points = 50; - glColor3f(0.0f, 0.0f, 0.0f); - for (unsigned i = 0; i < n_points; i++) { - glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points), m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points)); - } - glEnd(); + void post_cluster(const model& m) override { + glClearColor(1.0f, 1.0f, 1.0f, 1.0f); + glClear(GL_COLOR_BUFFER_BIT); + for (const Spin& s : m.s) { + glBegin(GL_POLYGON); + unsigned n_points = 50; + glColor3f(0.0f, 0.0f, 0.0f); + for (unsigned i = 0; i < n_points; i++) { + glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points), + m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points)); } - glFlush(); - - t1 += tmp; - t2 += tmp * tmp; - n++; + glEnd(); } + glFlush(); - void clear() { - t1 = 0; - t2 = 0; - n = 0; - } + t1 += tmp; + t2 += tmp * tmp; + n++; + } - double var() { - return (t2 - t1 * t1 / (double)n) / (double)n; - } + void clear() { + t1 = 0; + t2 = 0; + n = 0; + } + + double var() { return (t2 - t1 * t1 / (double)n) / (double)n; } }; std::function(const model&, randutils::mt19937_rng&)> eGen(double ε) { - return [ε] (const model& M, randutils::mt19937_rng& rng) -> Euclidean { + return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean { Vector t; Matrix m; @@ -85,7 +87,6 @@ std::function(const model&, randutils::mt19937_rng&)> eGen( Euclidean g(t, m); return g; }; - } int main(int argc, char* argv[]) { @@ -101,51 +102,50 @@ int main(int argc, char* argv[]) { 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); + 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); } } - double k = 1e2; - double a = 0.05; + double k = 1e8; + double a = 0.0; std::function&, const Spin&)> Z = - [L, a, k] (const Spin& s1, const Spin& s2) -> double { - Vector d = s1.x - s2.x; - - double σ = s1.s + s2.s; - double δ = σ - sqrt(d.transpose() * d); - - if (δ > - a * σ) { - return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2)); - } else if (δ > - 2 * a * σ) { - return 0.5 * k * pow(δ + 2 * a * σ, 2); - } else { - return 0; - } - }; + [L, a, k](const Spin& s1, const Spin& s2) -> double { + Vector d = s1.x - s2.x; + + double σ = s1.s + s2.s; + double δ = σ - sqrt(d.transpose() * d); + + if (δ > -a * σ) { + return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2)); + } else if (δ > -2 * a * σ) { + return 0.5 * k * pow(δ + 2 * a * σ, 2); + } else { + return 0; + } + }; - std::function)> B = - [L, H] (Spin s) -> double { - return H * s.x.norm(); - }; + std::function)> B = [L, H](Spin s) -> double { + return H * s.x.norm(); + }; - auto g = eGen(1); + auto g = eGen(0.25); animation A(L, 750, argc, argv); model sphere(1, Z, B); @@ -187,3 +187,4 @@ int main(int argc, char* argv[]) { return 0; } + diff --git a/spin.hpp b/spin.hpp new file mode 100644 index 0000000..d4d39cc --- /dev/null +++ b/spin.hpp @@ -0,0 +1,11 @@ + +#pragma once + +#include "vector.hpp" + +template class Spin { +public: + Vector x; + S s; +}; + diff --git a/torus_symmetries.hpp b/torus_symmetries.hpp new file mode 100644 index 0000000..17772de --- /dev/null +++ b/torus_symmetries.hpp @@ -0,0 +1,100 @@ + +#pragma once + +#include +#include +#include + +#include "matrix.hpp" +#include "vector.hpp" + +template void one_sequences(std::list>& sequences, unsigned level) { + if (level > 0) { + unsigned new_level = level - 1; + unsigned old_length = sequences.size(); + for (std::array& sequence : sequences) { + std::array new_sequence = sequence; + new_sequence[new_level] = -1; + sequences.push_front(new_sequence); + } + one_sequences(sequences, new_level); + } +} + +template std::vector> torus_vecs(U L) { + std::vector> vecs; + std::array ini_sequence; + ini_sequence.fill(1); + std::list> sequences; + sequences.push_back(ini_sequence); + one_sequences(sequences, D); + sequences.pop_front(); // don't want the identity matrix! + + for (std::array sequence : sequences) { + Vector v; + for (unsigned i = 0; i < D; i++) { + if (sequence[i] == 1) { + v(i) = 0; + } else { + v(i) = L / 2; + } + } + + vecs.push_back(v); + } + + return vecs; +} + +template std::vector> torus_mats() { + std::vector> mats; + + std::array ini_sequence; + ini_sequence.fill(1); + std::list> sequences; + sequences.push_back(ini_sequence); + + one_sequences(sequences, D); + + sequences.pop_front(); // don't want the identity matrix! + + for (std::array sequence : sequences) { + Matrix m; + for (unsigned i = 0; i < D; i++) { + for (unsigned j = 0; j < D; j++) { + if (i == j) { + m(i, j) = sequence[i]; + } else { + m(i, j) = 0; + } + } + } + + mats.push_back(m); + } + + for (unsigned i = 0; i < D; i++) { + for (unsigned j = 0; j < D; j++) { + if (i != j) { + Matrix m; + for (unsigned k = 0; k < D; k++) { + for (unsigned l = 0; l < D; l++) { + if ((k == i && l == j) || (k == j && l == i)) { + if (i < j) { + m(k, l) = 1; + } else { + m(k, l) = -1; + } + } else { + m(k, l) = 0; + } + } + } + mats.push_back(m); + } + } + } + + return mats; +} + diff --git a/vector.hpp b/vector.hpp new file mode 100644 index 0000000..2e87acd --- /dev/null +++ b/vector.hpp @@ -0,0 +1,20 @@ + +#pragma once + +#include + +template using Vector = Eigen::Matrix; + +template Vector diff(T L, Vector v1, Vector v2) { + Vector v; + + for (unsigned i = 0; i < D; i++) { + v(i) = std::abs(v1(i) - v2(i)); + if (v(i) > L / 2) { + v(i) = L - v(i); + } + } + + return v; +} + -- cgit v1.2.3-54-g00ecf