diff options
Diffstat (limited to 'space_wolff.hpp')
| -rw-r--r-- | space_wolff.hpp | 484 |
1 files changed, 14 insertions, 470 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp index 7345726..196b3b5 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -1,488 +1,36 @@ -#include <cmath> -#include <eigen3/Eigen/Dense> -#include <fstream> +#pragma once + +#include <array> #include <functional> -#include <iostream> #include <list> #include <queue> -#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" -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}}, - {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}}, - {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}}, - {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}}, - {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}}, - {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}}, - {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}, - {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}, - {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}}, - {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}}, - {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}}, - {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}}, - {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}}, - {{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 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) { - Vector<T, D> v; - - for (unsigned i = 0; i < D; i++) { - v(i) = std::abs(v1(i) - v2(i)); - if (v(i) > L / 2) { - v(i) = L - v(i); - } - } - - return v; -} - -template <class T, int D, class S> class Spin { -public: - Vector<T, D> x; - S s; -}; - -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; - -public: - Vector<T, D> t; - Matrix<T, D> r; - - /** brief TorusGroup - default constructor, constructs the identity - */ - TorusGroup(T L) : L(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; - } - } - } - - TorusGroup(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) { - 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; - - for (unsigned i = 0; i < D; i++) { - s_new.x(i) = fmod(L + s_new.x(i), L); - } - - return s_new; - } - - TorusGroup act(const TorusGroup& x) const { - Vector<T, D> tnew = r * x.t + t; - Matrix<T, D> rnew = r * x.r; - - for (unsigned i = 0; i < D; i++) { - tnew(i) = fmod(L + tnew(i), L); - } - - TorusGroup pnew(this->L, tnew, rnew); - - return pnew; - } - - TorusGroup inverse() const { - Vector<T, D> tnew = -r.transpose() * t; - Matrix<T, D> rnew = r.transpose(); - - TorusGroup pnew(this->L, tnew, rnew); - - return pnew; - } -}; - -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; - -public: - Octree(T L, unsigned depth) : L(L), N(pow(2, depth)), data(pow(N, D)){}; - - unsigned ind(Vector<T, D> x) const { - unsigned pos_ind = 0; - - for (unsigned i = 0; i < D; i++) { - pos_ind += pow(N, i) * (unsigned)std::floor(N * x(i) / L); - } - - assert(pos_ind < data.size()); - - return pos_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); }; - - 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, - std::set<unsigned>& seen) const { - ns.insert(data[ind].begin(), data[ind].end()); - seen.insert(ind); - - 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)); - - if (!seen.contains(nind)) { - seen.insert(nind); - nearest_neighbors_of(nind, depth - 1, ns, seen); - } - } - } - } - } -}; -*/ - -class quantity { -private: - double total; - double total2; - std::vector<double> C; - unsigned wait; - -public: - unsigned n; - std::list<double> hist; - - quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) { - n = 0; - total = 0; - total2 = 0; - } - - void add(double x) { - if (n > wait) { - hist.push_front(x); - if (hist.size() > C.size()) { - hist.pop_back(); - unsigned t = 0; - for (double a : hist) { - C[t] += a * x; - t++; - } - total += x; - total2 += pow(x, 2); - } - } - n++; - } - - double avg() const { return total / (n - wait); } - - double avg2() const { return total2 / (n - wait); } - - std::vector<double> ρ() const { - double C0 = C.front() / (n - wait); - double avg2 = pow(total / (n - wait), 2); - - std::vector<double> ρtmp; - - for (double Ct : C) { - ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); - } - - return ρtmp; - } - - std::array<double, 2> τ() const { - double τtmp = 0.5; - unsigned M = 1; - double c = 8.0; - - std::vector<double> ρ_tmp = this->ρ(); - - while (c * τtmp > M && M < C.size()) { - τtmp += ρ_tmp[M]; - M++; - } - - return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)}; - } - - double σ() const { - return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2)); - } - - double serr() const { return sqrt(this->σ()); } - - unsigned num_added() const { return n - wait; } -}; +#include "euclidean.hpp" +#include "octree.hpp" +#include "quantity.hpp" +#include "spin.hpp" template <class U, int D, class R, class S> class Model; template <class U, int D, class R, class S> class measurement { public: virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const R&){}; - virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*, const Spin<U, D, S>*, - const Spin<U, D, S>&, double){}; + virtual void plain_bond_visited(const Model<U, D, R, 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, R, S>&, const Spin<U, D, S>*, const Spin<U, D, S>&){}; - virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&, const Spin<U, D, S>&, - double){}; + virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&, + const Spin<U, D, S>&, double){}; virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){}; virtual void post_cluster(const Model<U, D, R, S>&){}; }; -template<int D> -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(); - for (std::array<double, D>& sequence : sequences) { - std::array<double, D> new_sequence = sequence; - new_sequence[new_level] = -1; - sequences.push_front(new_sequence); - } - one_sequences<D>(sequences, new_level); - } -} - -template <class U, int D> -std::vector<Vector<U, D>> torus_vecs(U L) { - std::vector<Vector<U, D>> vecs; - std::array<double, D> ini_sequence; - ini_sequence.fill(1); - std::list<std::array<double, D>> sequences; - sequences.push_back(ini_sequence); - one_sequences<D>(sequences, D); - sequences.pop_front(); // don't want the identity matrix! - - for (std::array<double, D> sequence : sequences) { - Vector<U, D> v; - for (unsigned i = 0; i < D; i++) { - if (sequence[i] == 1) { - v(i) = 0; - } else { - v(i) = L / 2; - } - } - - vecs.push_back(v); - } - - return vecs; -} - -template <class U, int D> -std::vector<Matrix<U, D>> torus_mats() { - std::vector<Matrix<U, D>> mats; - - std::array<double, D> ini_sequence; - ini_sequence.fill(1); - std::list<std::array<double, D>> sequences; - sequences.push_back(ini_sequence); - - one_sequences<D>(sequences, D); - - sequences.pop_front(); // don't want the identity matrix! - - for (std::array<double, D> sequence : sequences) { - Matrix<U, D> m; - for (unsigned i = 0; i < D; i++) { - for (unsigned j = 0; j < D; j++) { - if (i == j) { - m(i, j) = sequence[i]; - } else { - m(i, j) = 0; - } - } - } - - mats.push_back(m); - } - - for (unsigned i = 0; i < D; i++) { - for (unsigned j = 0; j < D; j++) { - if (i != j) { - Matrix<U, D> m; - for (unsigned k = 0; k < D; k++) { - for (unsigned l = 0; l < D; l++) { - if ((k == i && l == j) || (k == j && l == i)) { - if (i < j) { - m(k, l) = 1; - } else { - m(k, l) = -1; - } - } else { - m(k, l) = 0; - } - } - } - mats.push_back(m); - } - } - } - - return mats; -} - template <class U, int D, class R, class S> class Model { public: U L; @@ -493,14 +41,11 @@ public: 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) : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {} 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; queue.push(&(s[ind])); @@ -556,17 +101,15 @@ public: dict.remove(si); *si = si_new; dict.insert(si); - cluster_size++; } } } } - void resize(double T, double P) { - } + void resize(double T, double P) {} void wolff(double T, std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)> g, - measurement<U, D, R, S>& A, unsigned N) { + measurement<U, D, R, S>& A, unsigned N) { for (unsigned i = 0; i < N; i++) { R r = g(*this, rng); unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1)); @@ -579,3 +122,4 @@ public: } } }; + |