summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--space_wolff.hpp215
-rw-r--r--spheres.cpp49
2 files changed, 144 insertions, 120 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp
index 02c544d..1335f07 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -244,125 +244,136 @@ public:
unsigned num_added() const { return n - wait; }
};
-template <class U, int D, class S> class Model;
+template <class U, int D, class R, class S> class Model;
-template <class U, int D, class S> class measurement {
+template <class U, int D, class R, class S> class measurement {
public:
- virtual void pre_cluster(const Model<U, D, S>&, unsigned, const Euclidean<U, D>&){};
- virtual void plain_bond_visited(const Model<U, D, S>&, const Spin<U, D, S>*, const Spin<U, D, S>*,
+ 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_site_transformed(const Model<U, D, S>&, const Spin<U, D, S>*,
+ 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, 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_site_transformed(const Model<U, D, S>&, const Euclidean<U, D>&){};
+ virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){};
- virtual void post_cluster(const Model<U, D, S>&){};
+ virtual void post_cluster(const Model<U, D, R, S>&){};
};
-template <class U, int D, class S> class Model {
-public:
- U L;
- Euclidean<U, D> s0;
- std::vector<Spin<U, D, S>> s;
- Octree<U, D, S> dict;
- unsigned dDepth;
- unsigned nDepth;
- std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z;
- std::function<double(const Spin<U, D, S>&)> B;
- std::vector<Matrix<U, D>> mats;
- std::vector<Vector<U, D>> steps;
- double E;
- measurement<U, D, S>& A;
-
- 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(sequences, new_level);
+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(sequences, new_level);
}
+}
- 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, unsigned dDepth, unsigned nDepth,
- measurement<U, D, S>& A)
- : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), dDepth(dDepth), nDepth(nDepth), A(A) {
- std::array<double, D> ini_sequence;
- ini_sequence.fill(1);
- std::list<std::array<double, D>> sequences;
- sequences.push_back(ini_sequence);
+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);
+ 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);
+ }
- one_sequences(sequences, D);
+ return vecs;
+}
- sequences.pop_front(); // don't want the identity matrix!
+template <class U, int D>
+std::vector<Matrix<U, D>> torus_mats() {
+ std::vector<Matrix<U, D>> mats;
- 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;
- }
- }
- }
+ std::array<double, D> ini_sequence;
+ ini_sequence.fill(1);
+ std::list<std::array<double, D>> sequences;
+ sequences.push_back(ini_sequence);
- mats.push_back(m);
+ one_sequences(sequences, D);
- Vector<U, D> v;
- for (unsigned i = 0; i < D; i++) {
- if (sequence[i] == 1) {
- v(i) = 0;
+ 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 {
- v(i) = L / 2;
+ m(i, j) = 0;
}
}
-
- steps.push_back(v);
}
- 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;
- }
+ 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) = 0;
+ m(k, l) = -1;
}
+ } else {
+ m(k, l) = 0;
}
}
- mats.push_back(m);
}
+ mats.push_back(m);
}
}
}
- void update_energy() {
- E = 0;
- for (unsigned i = 0; i < s.size(); i++) {
- for (unsigned j = 0; j < i; j++) {
- E -= Z(s[i], s[j]);
- }
+ return mats;
+}
- E -= B(s0.inverse().act(s[i]));
- }
- }
+template <class U, int D, class R, class S> class Model {
+public:
+ U L;
+ R s0;
+ std::vector<Spin<U, D, S>> s;
+ Octree<U, D, S> dict;
+ unsigned dDepth;
+ unsigned nDepth;
+ std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z;
+ std::function<double(const Spin<U, D, S>&)> B;
+ std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g;
+ randutils::mt19937_rng rng;
+ measurement<U, D, R, S>& A;
+
+
+ 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, std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g, unsigned dDepth, unsigned nDepth,
+ measurement<U, D, R, S>& A)
+ : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), g(g), dDepth(dDepth), nDepth(nDepth), rng(), A(A) {}
- void step(double T, unsigned ind, Euclidean<U, D> r, std::mt19937& rng) {
+ void step(double T, unsigned ind, Euclidean<U, D> r) {
unsigned cluster_size = 0;
std::uniform_real_distribution<double> dist(0.0, 1.0);
@@ -427,37 +438,11 @@ public:
}
}
- void wolff(double T, unsigned N, std::mt19937& rng) {
- std::uniform_real_distribution<double> t_dist(0, L);
- std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
+ void wolff(double T, unsigned N) {
std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
- std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1);
for (unsigned i = 0; i < N; i++) {
- Vector<U, D> t;
- Matrix<U, D> m;
- unsigned flip = coin(rng);
- if (flip < mats.size()) {
- for (unsigned j = 0; j < D; j++) {
- t(j) = (U)t_dist(rng);
- }
- 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 = steps[flip - mats.size()];
- }
-
- Euclidean<U, D> g(L, t, m);
-
+ R g = g(*this, rng);
unsigned ind = ind_dist(rng);
A.pre_cluster(*this, ind, g);
@@ -465,8 +450,6 @@ public:
this->step(T, ind, g, rng);
A.post_cluster(*this);
-
- this->update_energy();
}
}
};
diff --git a/spheres.cpp b/spheres.cpp
index 25c5a67..9e67b87 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -2,7 +2,10 @@
#include "space_wolff.hpp"
#include <GL/glut.h>
-class animation : public measurement<double, 2, double> {
+const unsigned D = 2;
+typedef Model<double, D, Euclidean<double, D>, double> model;
+
+class animation : public measurement<double, 2, Euclidean<double, D>, double> {
public:
animation(double L, unsigned w, int argc, char *argv[]) {
glutInit(&argc, argv);
@@ -15,7 +18,7 @@ class animation : public measurement<double, 2, double> {
gluOrtho2D(-1, L + 1, -1 , L + 1);
}
- void post_cluster(const Model<double, 2, double>& m) override {
+ 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) {
@@ -32,6 +35,41 @@ class animation : public measurement<double, 2, double> {
}
};
+std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs) {
+ return [&mats, &vecs] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+ std::uniform_real_distribution<double> t_dist(0, M.L);
+ std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
+ std::uniform_int_distribution<unsigned> ind_dist(0, M.s.size() - 1);
+ std::uniform_int_distribution<unsigned> coin(0, mats.size() + vecs.size() - 1);
+
+ Vector<double, D> t;
+ Matrix<double, D> m;
+ unsigned flip = coin(rng);
+ if (flip < mats.size()) {
+ for (unsigned j = 0; j < D; j++) {
+ t(j) = (double)t_dist(rng);
+ }
+ 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()];
+ }
+
+ Euclidean<double, D> g(M.L, t, m);
+ return g;
+ };
+
+}
+
int main(int argc, char* argv[]) {
const unsigned D = 2;
@@ -85,8 +123,11 @@ int main(int argc, char* argv[]) {
return H * sin(2 * M_PI * 3 * s.x(0) / L);
};
+ std::vector<Matrix<double, D>> mats = torus_mats<double, D>();
+ std::vector<Vector<double, D>> vecs = torus_vecs<double, D>(L);
+ auto g = eGen(mats, vecs);
animation A(L, 750, argc, argv);
- Model<double, D, double> sphere(L, Z, B, std::floor(log2(L)), 2, A);
+ model sphere(L, Z, B, g, std::floor(log2(L)), 2, A);
randutils::auto_seed_128 seeds;
std::mt19937 rng{seeds};
@@ -102,7 +143,7 @@ int main(int argc, char* argv[]) {
}
- sphere.wolff(T, N, rng);
+ sphere.wolff(T, N);
std::ofstream snapfile;
snapfile.open("sphere_snap.dat");