summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--space_wolff.hpp123
1 files changed, 79 insertions, 44 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp
index 81ff776..a78116d 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -1,5 +1,4 @@
-#include "randutils/randutils.hpp"
#include <cmath>
#include <eigen3/Eigen/Dense>
#include <fstream>
@@ -11,6 +10,8 @@
#include <set>
#include <vector>
+#include "randutils/randutils.hpp"
+
const std::array<std::array<unsigned, 16>, 16> smiley = {
{{{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}},
{{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
@@ -29,13 +30,16 @@ const std::array<std::array<unsigned, 16>, 16> smiley = {
{{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
{{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}}};
-template <class T, int D> using Vector = Eigen::Matrix<T, D, 1>;
+template <class T, int D>
+using Vector = Eigen::Matrix<T, D, 1>;
-template <class T, int D> using Matrix = Eigen::Matrix<T, D, D>;
+template <class T, int D>
+using Matrix = Eigen::Matrix<T, D, D>;
/** brief diff - periodic subtraction of vectors
*/
-template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) {
+template <class T, int D>
+Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) {
Vector<T, D> v;
for (unsigned i = 0; i < D; i++) {
@@ -48,17 +52,19 @@ template <class T, int D> Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v
return v;
}
-template <class T, int D, class S> class Spin {
-public:
+template <class T, int D, class S>
+class Spin {
+ public:
Vector<T, D> x;
S s;
};
-template <class T, int D> class Euclidean {
-private:
+template <class T, int D>
+class Euclidean {
+ private:
T L;
-public:
+ public:
Vector<T, D> t;
Matrix<T, D> r;
@@ -79,7 +85,8 @@ public:
r = r0;
}
- template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+ 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;
@@ -115,13 +122,14 @@ public:
}
};
-template <class T, int D, class S> class Octree {
-private:
+template <class T, int D, class S>
+class Octree {
+ private:
T L;
unsigned N;
std::vector<std::set<Spin<T, D, S>*>> data;
-public:
+ public:
Octree(T L, unsigned depth) : L(L), N(pow(2, depth)), data(pow(N, D)){};
unsigned ind(Vector<T, D> x) const {
@@ -140,14 +148,16 @@ public:
void remove(Spin<T, D, S>* s) { data[ind(s->x)].erase(s); };
- std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x, unsigned depth) const {
+ std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x,
+ unsigned depth) const {
std::set<Spin<T, D, S>*> ns;
std::set<unsigned> seen;
nearest_neighbors_of(ind(x), depth, ns, seen);
return ns;
};
- void nearest_neighbors_of(unsigned ind, unsigned depth, std::set<Spin<T, D, S>*>& ns,
+ void nearest_neighbors_of(unsigned ind, unsigned depth,
+ std::set<Spin<T, D, S>*>& ns,
std::set<unsigned>& seen) const {
ns.insert(data[ind].begin(), data[ind].end());
seen.insert(ind);
@@ -155,8 +165,9 @@ public:
if (depth > 0) {
for (unsigned i = 0; i < D; i++) {
for (signed j : {-1, 1}) {
- unsigned nind = pow(N, i + 1) * (ind / ((unsigned)pow(N, i + 1))) +
- fmod(pow(N, i + 1) + ind + j * pow(N, i), pow(N, i + 1));
+ unsigned nind =
+ pow(N, i + 1) * (ind / ((unsigned)pow(N, i + 1))) +
+ fmod(pow(N, i + 1) + ind + j * pow(N, i), pow(N, i + 1));
if (!seen.contains(nind)) {
seen.insert(nind);
@@ -169,13 +180,13 @@ public:
};
class quantity {
-private:
+ private:
double total;
double total2;
std::vector<double> C;
unsigned wait;
-public:
+ public:
unsigned n;
std::list<double> hist;
@@ -235,7 +246,8 @@ public:
}
double σ() const {
- return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+ return 2.0 / (n - wait) * this->τ()[0] *
+ (C[0] / (n - wait) - pow(this->avg(), 2));
}
double serr() const { return sqrt(this->σ()); }
@@ -243,23 +255,32 @@ public:
unsigned num_added() const { return n - wait; }
};
-template <class U, int D, class S> class Model;
-
template <class U, int D, 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>*, const Spin<U, D, S>& , double) {};
- virtual void plain_site_transformed(const Model<U, D, 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>&, double) {};
- virtual void ghost_site_transformed(const Model<U, D, S>&, const Euclidean<U, D>&) {};
+class Model;
- virtual void post_cluster(const Model<U, D, S>&) {};
- };
+template <class U, int D, 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>*, const Spin<U, D, S>&,
+ double){};
+ virtual void plain_site_transformed(const Model<U, D, 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>&, double){};
+ virtual void ghost_site_transformed(const Model<U, D, S>&,
+ const Euclidean<U, D>&){};
+
+ virtual void post_cluster(const Model<U, D, S>&){};
+};
-template <class U, int D, class S> class Model {
-public:
+template <class U, int D, class S>
+class Model {
+ public:
U L;
Euclidean<U, D> s0;
std::vector<Spin<U, D, S>> s;
@@ -271,9 +292,10 @@ public:
std::vector<Matrix<U, D>> mats;
std::vector<Vector<U, D>> steps;
double E;
- measurement<U, D, S> &A;
+ measurement<U, D, S>& A;
- void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+ 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();
@@ -286,9 +308,18 @@ public:
}
}
- 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) {
+ 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;
@@ -296,7 +327,7 @@ public:
one_sequences(sequences, D);
- sequences.pop_front(); // don't want the identity matrix!
+ sequences.pop_front(); // don't want the identity matrix!
for (std::array<double, D> sequence : sequences) {
Matrix<U, D> m;
@@ -390,10 +421,13 @@ 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, nDepth);
+ std::set<Spin<U, D, S>*> new_neighbors =
+ dict.neighbors(si_new.x, nDepth);
- all_neighbors.insert(current_neighbors.begin(), current_neighbors.end());
+ all_neighbors.insert(current_neighbors.begin(),
+ current_neighbors.end());
all_neighbors.insert(new_neighbors.begin(), new_neighbors.end());
all_neighbors.insert(NULL);
for (Spin<U, D, S>* sj : all_neighbors) {
@@ -427,7 +461,8 @@ public:
std::uniform_real_distribution<double> t_dist(0, L);
std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
- std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.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;