diff options
-rw-r--r-- | space_wolff.hpp | 113 |
1 files changed, 41 insertions, 72 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp index a78116d..02c544d 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -30,16 +30,13 @@ 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++) { @@ -52,19 +49,17 @@ Vector<T, D> diff(T L, Vector<T, D> v1, Vector<T, D> v2) { 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; @@ -85,8 +80,7 @@ class Euclidean { 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; @@ -122,14 +116,13 @@ class Euclidean { } }; -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 { @@ -148,16 +141,14 @@ class Octree { 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); @@ -165,9 +156,8 @@ class Octree { 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); @@ -180,13 +170,13 @@ class Octree { }; class quantity { - private: +private: double total; double total2; std::vector<double> C; unsigned wait; - public: +public: unsigned n; std::list<double> hist; @@ -246,8 +236,7 @@ class quantity { } 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->σ()); } @@ -255,32 +244,25 @@ class quantity { unsigned num_added() const { return n - wait; } }; -template <class U, int D, class S> -class Model; +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>*, +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 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; @@ -294,8 +276,7 @@ class Model { double E; 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(); @@ -308,18 +289,10 @@ class Model { } } - 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; @@ -327,7 +300,7 @@ class Model { 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; @@ -421,13 +394,10 @@ class Model { } 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) { @@ -461,8 +431,7 @@ class Model { 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; |