diff options
Diffstat (limited to 'space_wolff.hpp')
-rw-r--r-- | space_wolff.hpp | 157 |
1 files changed, 142 insertions, 15 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()); |