summaryrefslogtreecommitdiff
path: root/space_wolff.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'space_wolff.hpp')
-rw-r--r--space_wolff.hpp157
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());