summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--euclidean.hpp116
-rw-r--r--matrix.hpp7
-rw-r--r--octree.hpp87
-rw-r--r--quantity.hpp83
-rw-r--r--smiley.hpp23
-rw-r--r--space_wolff.hpp484
-rw-r--r--spheres.cpp230
-rw-r--r--spheres_infinite.cpp181
-rw-r--r--spin.hpp11
-rw-r--r--torus_symmetries.hpp100
-rw-r--r--vector.hpp20
11 files changed, 669 insertions, 673 deletions
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 T, int D> class Euclidean {
+public:
+ Vector<T, D> t;
+ Matrix<T, D> 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<T, D> t0, Matrix<T, D> r0) {
+ t = t0;
+ r = r0;
+ }
+
+ template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+ Spin<T, D, S> s_new;
+
+ s_new.x = t + r * s.x;
+ s_new.s = s.s;
+
+ return s_new;
+ }
+
+ Euclidean act(const Euclidean& x) const {
+ Vector<T, D> tnew = r * x.t + t;
+ Matrix<T, D> rnew = r * x.r;
+
+ Euclidean pnew(tnew, rnew);
+
+ return pnew;
+ }
+
+ Euclidean inverse() const {
+ Vector<T, D> tnew = -r.transpose() * t;
+ Matrix<T, D> rnew = r.transpose();
+
+ Euclidean pnew(tnew, rnew);
+
+ return pnew;
+ }
+};
+
+template <class T, int D> class TorusGroup {
+private:
+ T L;
+
+public:
+ Vector<T, D> t;
+ Matrix<T, D> 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<T, D> t0, Matrix<T, D> r0) : L(L) {
+ t = t0;
+ r = r0;
+ }
+
+ template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+ Spin<T, D, S> 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<T, D> tnew = r * x.t + t;
+ Matrix<T, D> 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<T, D> tnew = -r.transpose() * t;
+ Matrix<T, D> 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 <eigen3/Eigen/Dense>
+
+template <class T, int D> using Matrix = Eigen::Matrix<T, D, D>;
+
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 <array>
+#include <set>
+#include <unordered_map>
+
+#include "spin.hpp"
+#include "vector.hpp"
+
+namespace std {
+template <typename T, size_t N> struct hash<array<T, N>> {
+ typedef array<T, N> argument_type;
+ typedef size_t result_type;
+
+ result_type operator()(const argument_type& a) const {
+ hash<T> 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 T, int D, class S> class Octree {
+private:
+ unsigned N;
+ T L;
+ std::unordered_map<std::array<signed, D>, std::set<Spin<T, D, S>*>> data;
+
+public:
+ Octree(T L, unsigned N) {
+ L = L;
+ N = N;
+ }
+
+ std::array<signed, D> ind(Vector<T, D> x) const {
+ std::array<signed, D> ind;
+
+ for (unsigned i = 0; i < D; i++) {
+ ind[i] = (signed)std::floor(N * x(i) / L);
+ }
+
+ return ind;
+ }
+
+ void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); };
+
+ void remove(Spin<T, D, S>* s) {
+ data[ind(s->x)].erase(s);
+ if (data[ind(s->x)].empty()) {
+ data.erase(ind(s->x));
+ }
+ };
+
+ std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x) const {
+ std::array<signed, D> i0 = ind(x);
+ std::set<Spin<T, D, S>*> ns;
+
+ nearest_neighbors_of(i0, D + 1, ns);
+
+ return ns;
+ };
+
+ void nearest_neighbors_of(std::array<signed, D> i0, unsigned depth,
+ std::set<Spin<T, D, S>*>& 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<signed, D> 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 <vector>
+#include <list>
+#include <cmath>
+#include <array>
+
+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) {
+ 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<double> ρ() const {
+ double C0 = C.front() / (n - wait);
+ double avg2 = pow(total / (n - wait), 2);
+
+ std::vector<double> ρtmp;
+
+ for (double Ct : C) {
+ ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
+ }
+
+ return ρtmp;
+ }
+
+ std::array<double, 2> τ() const {
+ double τtmp = 0.5;
+ unsigned M = 1;
+ double c = 8.0;
+
+ std::vector<double> ρ_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 <array>
+
+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}}}};
+
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 <cmath>
-#include <eigen3/Eigen/Dense>
-#include <fstream>
+#pragma once
+
+#include <array>
#include <functional>
-#include <iostream>
#include <list>
#include <queue>
-#include <random>
-#include <set>
#include <vector>
-#include <unordered_map>
-#include <array>
-#include <unordered_map>
-
-namespace std {
-template <typename T, size_t N> struct hash<array<T, N>> {
- typedef array<T, N> argument_type;
- typedef size_t result_type;
-
- result_type operator()(const argument_type& a) const {
- hash<T> 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<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>;
-
-/** brief diff - periodic subtraction of vectors
- */
-template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) {
- Vector<T, D> 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 T, int D, class S> class Spin {
-public:
- Vector<T, D> x;
- S s;
-};
-
-template <class T, int D> class Euclidean {
-public:
- Vector<T, D> t;
- Matrix<T, D> 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<T, D> t0, Matrix<T, D> r0) {
- t = t0;
- r = r0;
- }
-
- template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
- Spin<T, D, S> s_new;
-
- s_new.x = t + r * s.x;
- s_new.s = s.s;
-
- return s_new;
- }
-
- Euclidean act(const Euclidean& x) const {
- Vector<T, D> tnew = r * x.t + t;
- Matrix<T, D> rnew = r * x.r;
-
- Euclidean pnew(tnew, rnew);
-
- return pnew;
- }
-
- Euclidean inverse() const {
- Vector<T, D> tnew = -r.transpose() * t;
- Matrix<T, D> rnew = r.transpose();
-
- Euclidean pnew(tnew, rnew);
-
- return pnew;
- }
-};
-
-template <class T, int D> class TorusGroup {
-private:
- T L;
-
-public:
- Vector<T, D> t;
- Matrix<T, D> 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<T, D> t0, Matrix<T, D> r0) : L(L) {
- t = t0;
- r = r0;
- }
-
- template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
- Spin<T, D, S> 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<T, D> tnew = r * x.t + t;
- Matrix<T, D> 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<T, D> tnew = -r.transpose() * t;
- Matrix<T, D> rnew = r.transpose();
-
- TorusGroup pnew(this->L, tnew, rnew);
-
- return pnew;
- }
-};
-
-template <class T, int D, class S> class Octree {
-private:
- unsigned N;
- T L;
- std::unordered_map<std::array<signed, D>, std::set<Spin<T, D, S>*>> data;
-
-public:
- Octree(T L, unsigned N) {
- L = L;
- N = N;
- }
-
- std::array<signed, D> ind(Vector<T, D> x) const {
- std::array<signed, D> ind;
-
- for (unsigned i = 0; i < D; i++) {
- ind[i] = (signed)std::floor(N * x(i) / L);
- }
-
- return ind;
- }
-
- void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); };
-
- void remove(Spin<T, D, S>* s) {
- data[ind(s->x)].erase(s);
- if (data[ind(s->x)].empty()) {
- data.erase(ind(s->x));
- }
- };
-
- std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x) const {
- std::array<signed, D> i0 = ind(x);
- std::set<Spin<T, D, S>*> ns;
-
- nearest_neighbors_of(i0, D + 1, ns);
-
- return ns;
- };
-
- void nearest_neighbors_of(std::array<signed, D> i0, unsigned depth, std::set<Spin<T, D, S>*>& 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<signed, D> 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 T, int D, class S> class Octree {
-private:
- T L;
- unsigned N;
- std::vector<std::set<Spin<T, D, S>*>> data;
-
-public:
- Octree(T L, unsigned depth) : L(L), N(pow(2, depth)), data(pow(N, D)){};
-
- unsigned ind(Vector<T, D> 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<T, D, S>* s) { data[ind(s->x)].insert(s); };
-
- void remove(Spin<T, D, S>* s) { data[ind(s->x)].erase(s); };
-
- std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x, unsigned depth) const {
- std::set<Spin<T, D, S>*> ns;
- std::set<unsigned> seen;
- nearest_neighbors_of(ind(x), depth, ns, seen);
- return ns;
- };
-
- void nearest_neighbors_of(unsigned ind, unsigned depth, std::set<Spin<T, D, S>*>& ns,
- std::set<unsigned>& 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<double> C;
- unsigned wait;
-
-public:
- unsigned n;
- std::list<double> 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<double> ρ() const {
- double C0 = C.front() / (n - wait);
- double avg2 = pow(total / (n - wait), 2);
-
- std::vector<double> ρtmp;
-
- for (double Ct : C) {
- ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
- }
-
- return ρtmp;
- }
-
- std::array<double, 2> τ() const {
- double τtmp = 0.5;
- unsigned M = 1;
- double c = 8.0;
-
- std::vector<double> ρ_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 U, int D, class R, class S> class Model;
template <class U, int D, class R, class S> class measurement {
public:
virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const R&){};
- virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*, const Spin<U, D, S>*,
- const Spin<U, D, S>&, double){};
+ virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*,
+ const Spin<U, D, S>*, const Spin<U, D, S>&, double){};
virtual void plain_site_transformed(const Model<U, D, R, S>&, const Spin<U, D, S>*,
const Spin<U, D, S>&){};
- virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&, const Spin<U, D, S>&,
- double){};
+ virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&,
+ const Spin<U, D, S>&, double){};
virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){};
virtual void post_cluster(const Model<U, D, R, S>&){};
};
-template<int D>
-void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
- if (level > 0) {
- unsigned new_level = level - 1;
- unsigned old_length = sequences.size();
- for (std::array<double, D>& sequence : sequences) {
- std::array<double, D> new_sequence = sequence;
- new_sequence[new_level] = -1;
- sequences.push_front(new_sequence);
- }
- one_sequences<D>(sequences, new_level);
- }
-}
-
-template <class U, int D>
-std::vector<Vector<U, D>> torus_vecs(U L) {
- std::vector<Vector<U, D>> vecs;
- std::array<double, D> ini_sequence;
- ini_sequence.fill(1);
- std::list<std::array<double, D>> sequences;
- sequences.push_back(ini_sequence);
- one_sequences<D>(sequences, D);
- sequences.pop_front(); // don't want the identity matrix!
-
- for (std::array<double, D> sequence : sequences) {
- Vector<U, D> 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 <class U, int D>
-std::vector<Matrix<U, D>> torus_mats() {
- std::vector<Matrix<U, D>> mats;
-
- std::array<double, D> ini_sequence;
- ini_sequence.fill(1);
- std::list<std::array<double, D>> sequences;
- sequences.push_back(ini_sequence);
-
- one_sequences<D>(sequences, D);
-
- sequences.pop_front(); // don't want the identity matrix!
-
- for (std::array<double, D> sequence : sequences) {
- Matrix<U, D> 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<U, D> 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 U, int D, class R, class S> class Model {
public:
U L;
@@ -493,14 +41,11 @@ public:
std::function<double(const Spin<U, D, S>&)> B;
randutils::mt19937_rng rng;
-
Model(U L, std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z,
std::function<double(const Spin<U, D, S>&)> B)
: L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {}
void step(double T, unsigned ind, R& r, measurement<U, D, R, S>& A) {
- unsigned cluster_size = 0;
-
std::queue<Spin<U, D, S>*> 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<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g,
- measurement<U, D, R, S>& A, unsigned N) {
+ measurement<U, D, R, S>& 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 <GL/glut.h>
+#include <fstream>
+#include <iostream>
+
+#include "space_wolff.hpp"
+#include "torus_symmetries.hpp"
const unsigned D = 2;
typedef Model<double, D, TorusGroup<double, D>, double> model;
class animation : public measurement<double, 2, TorusGroup<double, D>, 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<double, D>&) override {
- tmp = 0;
- }
+ void pre_cluster(const model&, unsigned, const TorusGroup<double, D>&) override { tmp = 0; }
- void plain_site_transformed(const model&, const Spin<double, D, double>*, const Spin<double, D, double>&) override {
- tmp++;
- }
+ void plain_site_transformed(const model&, const Spin<double, D, double>*,
+ const Spin<double, D, double>&) 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<double, 2, double>& 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<double, 2, double>& 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<TorusGroup<double, D>(const model&, randutils::mt19937_rng&)> eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs, double ε, double L) {
- return [&mats, &vecs, L, ε] (const model& M, randutils::mt19937_rng& rng) -> TorusGroup<double, 2> {
- Vector<double, D> t;
- Matrix<double, D> 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<double, std::normal_distribution>(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<TorusGroup<double, D>(const model&, randutils::mt19937_rng&)>
+eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs,
+ double ε, double L) {
+ return
+ [&mats, &vecs, L, ε](const model& M, randutils::mt19937_rng& rng) -> TorusGroup<double, 2> {
+ Vector<double, D> t;
+ Matrix<double, D> 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<double, std::normal_distribution>(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<double, D> g(M.L, t, m);
- return g;
- };
+ t = vecs[flip - mats.size()];
+ }
+ TorusGroup<double, D> 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<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
- [L, a, k] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
- Vector<double, D> 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<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
+ Vector<double, D> d = diff(L, s1.x, s2.x);
+
+ double σ = s1.s + s2.s;
+ double δ = σ - sqrt(d.transpose() * d);
- std::function<double(Spin<double, D, double>)> B =
- [L, H] (Spin<double, D, double> 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<double(Spin<double, D, double>)> B = [L, H](Spin<double, D, double> s) -> double {
+ return H * s.x(1);
+ };
std::vector<Matrix<double, D>> mats = torus_mats<double, D>();
std::vector<Vector<double, D>> vecs = torus_vecs<double, D>(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 <fstream>
+#include <iostream>
+
#include "space_wolff.hpp"
#include <GL/glut.h>
@@ -6,67 +9,66 @@ const unsigned D = 2;
typedef Model<double, D, Euclidean<double, D>, double> model;
class animation : public measurement<double, 2, Euclidean<double, D>, 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<double, D>&) override {
- tmp = 0;
- }
+ void pre_cluster(const model&, unsigned, const Euclidean<double, D>&) override { tmp = 0; }
- void plain_site_transformed(const model&, const Spin<double, D, double>*, const Spin<double, D, double>&) override {
- tmp++;
- }
+ void plain_site_transformed(const model&, const Spin<double, D, double>*,
+ const Spin<double, D, double>&) 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<double, 2, double>& 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<double, 2, double>& 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<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(double ε) {
- return [ε] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+ return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
Vector<double, D> t;
Matrix<double, D> m;
@@ -85,7 +87,6 @@ std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(
Euclidean<double, D> 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<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
- [L, a, k] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
- Vector<double, D> 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<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
+ Vector<double, D> 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<double(Spin<double, D, double>)> B =
- [L, H] (Spin<double, D, double> s) -> double {
- return H * s.x.norm();
- };
+ std::function<double(Spin<double, D, double>)> B = [L, H](Spin<double, D, double> 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 T, int D, class S> class Spin {
+public:
+ Vector<T, D> 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 <array>
+#include <list>
+#include <vector>
+
+#include "matrix.hpp"
+#include "vector.hpp"
+
+template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+ if (level > 0) {
+ unsigned new_level = level - 1;
+ unsigned old_length = sequences.size();
+ for (std::array<double, D>& sequence : sequences) {
+ std::array<double, D> new_sequence = sequence;
+ new_sequence[new_level] = -1;
+ sequences.push_front(new_sequence);
+ }
+ one_sequences<D>(sequences, new_level);
+ }
+}
+
+template <class U, int D> std::vector<Vector<U, D>> torus_vecs(U L) {
+ std::vector<Vector<U, D>> vecs;
+ std::array<double, D> ini_sequence;
+ ini_sequence.fill(1);
+ std::list<std::array<double, D>> sequences;
+ sequences.push_back(ini_sequence);
+ one_sequences<D>(sequences, D);
+ sequences.pop_front(); // don't want the identity matrix!
+
+ for (std::array<double, D> sequence : sequences) {
+ Vector<U, D> 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 <class U, int D> std::vector<Matrix<U, D>> torus_mats() {
+ std::vector<Matrix<U, D>> mats;
+
+ std::array<double, D> ini_sequence;
+ ini_sequence.fill(1);
+ std::list<std::array<double, D>> sequences;
+ sequences.push_back(ini_sequence);
+
+ one_sequences<D>(sequences, D);
+
+ sequences.pop_front(); // don't want the identity matrix!
+
+ for (std::array<double, D> sequence : sequences) {
+ Matrix<U, D> 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<U, D> 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 <eigen3/Eigen/Dense>
+
+template <class T, int D> using Vector = Eigen::Matrix<T, D, 1>;
+
+template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) {
+ Vector<T, D> 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;
+}
+