diff options
| -rw-r--r-- | space_wolff.hpp | 157 | ||||
| -rw-r--r-- | spheres.cpp | 18 | ||||
| -rw-r--r-- | spheres_infinite.cpp | 189 |
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; +} |