summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--space_wolff.hpp157
-rw-r--r--spheres.cpp18
-rw-r--r--spheres_infinite.cpp189
3 files changed, 340 insertions, 24 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp
index 615fb93..7345726 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -9,6 +9,25 @@
#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"
@@ -56,6 +75,54 @@ public:
};
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;
@@ -63,9 +130,9 @@ public:
Vector<T, D> t;
Matrix<T, D> r;
- /** brief Euclidean - default constructor, constructs the identity
+ /** brief TorusGroup - default constructor, constructs the identity
*/
- Euclidean(T L) : L(L) {
+ TorusGroup(T L) : L(L) {
for (unsigned i = 0; i < D; i++) {
t(i) = 0;
r(i, i) = 1;
@@ -75,7 +142,7 @@ public:
}
}
- Euclidean(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
+ TorusGroup(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
t = t0;
r = r0;
}
@@ -93,7 +160,7 @@ public:
return s_new;
}
- Euclidean act(const Euclidean& x) const {
+ TorusGroup act(const TorusGroup& x) const {
Vector<T, D> tnew = r * x.t + t;
Matrix<T, D> rnew = r * x.r;
@@ -101,16 +168,16 @@ public:
tnew(i) = fmod(L + tnew(i), L);
}
- Euclidean pnew(this->L, tnew, rnew);
+ TorusGroup pnew(this->L, tnew, rnew);
return pnew;
}
- Euclidean inverse() const {
+ TorusGroup inverse() const {
Vector<T, D> tnew = -r.transpose() * t;
Matrix<T, D> rnew = r.transpose();
- Euclidean pnew(this->L, tnew, rnew);
+ TorusGroup pnew(this->L, tnew, rnew);
return pnew;
}
@@ -118,6 +185,67 @@ public:
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;
@@ -168,6 +296,7 @@ public:
}
}
};
+*/
class quantity {
private:
@@ -360,18 +489,16 @@ public:
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;
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, unsigned dDepth, unsigned nDepth)
- : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), dDepth(dDepth), nDepth(nDepth), rng() {}
+ 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, Euclidean<U, D> r, measurement<U, D, R, S>& A) {
+ 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;
@@ -387,7 +514,7 @@ public:
visited.insert(si);
if (si == NULL) {
- Euclidean<U, D> s0_new = r.act(s0);
+ R s0_new = r.act(s0);
for (Spin<U, D, S>& ss : s) {
Spin<U, D, S> s0s_old = s0.inverse().act(ss);
Spin<U, D, S> s0s_new = s0_new.inverse().act(ss);
@@ -402,8 +529,8 @@ public:
} else {
Spin<U, D, S> si_new = r.act(*si);
std::set<Spin<U, D, S>*> all_neighbors;
- std::set<Spin<U, D, S>*> current_neighbors = dict.neighbors(si->x, nDepth);
- std::set<Spin<U, D, S>*> new_neighbors = dict.neighbors(si_new.x, nDepth);
+ std::set<Spin<U, D, S>*> current_neighbors = dict.neighbors(si->x);
+ std::set<Spin<U, D, S>*> new_neighbors = dict.neighbors(si_new.x);
all_neighbors.insert(current_neighbors.begin(), current_neighbors.end());
all_neighbors.insert(new_neighbors.begin(), new_neighbors.end());
diff --git a/spheres.cpp b/spheres.cpp
index f3dde42..97c1e18 100644
--- a/spheres.cpp
+++ b/spheres.cpp
@@ -3,9 +3,9 @@
#include <GL/glut.h>
const unsigned D = 2;
-typedef Model<double, D, Euclidean<double, D>, double> model;
+typedef Model<double, D, TorusGroup<double, D>, double> model;
-class animation : public measurement<double, 2, Euclidean<double, D>, double> {
+class animation : public measurement<double, 2, TorusGroup<double, D>, double> {
private:
uint64_t t1;
uint64_t t2;
@@ -27,7 +27,7 @@ class animation : public measurement<double, 2, Euclidean<double, D>, double> {
gluOrtho2D(-1, L + 1, -1 , L + 1);
}
- void pre_cluster(const model&, unsigned, const Euclidean<double, D>&) override {
+ void pre_cluster(const model&, unsigned, const TorusGroup<double, D>&) override {
tmp = 0;
}
@@ -65,8 +65,8 @@ class animation : public measurement<double, 2, Euclidean<double, D>, 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, double ε, double L) {
- return [&mats, &vecs, L, ε] (const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+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));
@@ -91,7 +91,7 @@ std::function<Euclidean<double, D>(const model&, randutils::mt19937_rng&)> eGen(
t = vecs[flip - mats.size()];
}
- Euclidean<double, D> g(M.L, t, m);
+ TorusGroup<double, D> g(M.L, t, m);
return g;
};
@@ -130,8 +130,8 @@ int main(int argc, char* argv[]) {
}
}
- double k = 100.0;
- double a = 0.2;
+ double k = 1e2;
+ 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 {
@@ -158,7 +158,7 @@ int main(int argc, char* argv[]) {
std::vector<Vector<double, D>> vecs = torus_vecs<double, D>(L);
auto g = eGen(mats, vecs, L, L);
animation A(L, 750, argc, argv);
- model sphere(L, Z, B, std::floor(log2(L)), 2);
+ model sphere(L, Z, B);
randutils::mt19937_rng rng;
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
new file mode 100644
index 0000000..f2f998a
--- /dev/null
+++ b/spheres_infinite.cpp
@@ -0,0 +1,189 @@
+
+#include "space_wolff.hpp"
+#include <GL/glut.h>
+
+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);
+ }
+
+ 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 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();
+ }
+ glFlush();
+
+ t1 += tmp;
+ t2 += tmp * tmp;
+ 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> {
+ Vector<double, D> t;
+ Matrix<double, D> m;
+
+ double θ = rng.uniform((double)0.0, 2 * M_PI);
+ m(0, 0) = cos(θ);
+ m(1, 1) = cos(θ);
+ m(0, 1) = sin(θ);
+ m(1, 0) = -sin(θ);
+
+ 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) = rng.variate<double, std::normal_distribution>(0.0, ε);
+ }
+
+ Euclidean<double, D> g(t, m);
+ return g;
+ };
+
+}
+
+int main(int argc, char* argv[]) {
+ const unsigned D = 2;
+
+ double L = 32;
+ unsigned N = 1000;
+ double T = 2.0 / log(1.0 + sqrt(2.0));
+ double H = 1.0;
+ unsigned n = 25;
+
+ int opt;
+
+ 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);
+ }
+ }
+
+ double k = 1e2;
+ 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 = 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();
+ };
+
+ auto g = eGen(1);
+ animation A(L, 750, argc, argv);
+ model sphere(1, Z, B);
+
+ randutils::mt19937_rng rng;
+
+ sphere.s.reserve(n);
+
+ unsigned nx = floor(sqrt(n));
+ for (unsigned i = 0; i < n; i++) {
+ Vector<double, D> pos = {(i / nx) * L / nx, (i % nx) * L / nx};
+ sphere.s.push_back({pos, 0.5});
+ sphere.dict.insert(&sphere.s.back());
+ }
+
+ sphere.wolff(T, g, A, N);
+ std::ofstream outfile;
+ outfile.open("test.dat");
+
+ for (signed i = -10; i <= 10; i++) {
+ A.clear();
+ double ε = pow(2, -4 + i / 2.0);
+ auto gn = eGen(ε);
+ sphere.wolff(T, gn, A, N);
+ outfile << ε << " " << A.var() / sphere.s.size() << std::endl;
+ std::cout << ε << " " << A.var() / sphere.s.size() << std::endl;
+ }
+
+ outfile.close();
+
+ std::ofstream snapfile;
+ snapfile.open("sphere_snap.dat");
+
+ for (Spin<double, D, double> s : sphere.s) {
+ Spin<double, D, double> rs = sphere.s0.inverse().act(s);
+ snapfile << rs.x.transpose() << "\n";
+ }
+
+ snapfile.close();
+
+ return 0;
+}