diff options
| -rw-r--r-- | space_wolff.hpp | 701 | ||||
| -rw-r--r-- | spheres.cpp | 50 |
2 files changed, 359 insertions, 392 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp index c87697e..cc4cd4f 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -1,45 +1,42 @@ -#include <vector> -#include <list> -#include <set> -#include <iostream> +#include "randutils/randutils.hpp" +#include <cmath> +#include <eigen3/Eigen/Dense> #include <fstream> #include <functional> -#include <random> +#include <iostream> +#include <list> #include <queue> -#include <cmath> -#include <eigen3/Eigen/Dense> -#include "randutils/randutils.hpp" +#include <random> +#include <set> +#include <vector> -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>; - -template <class T, int D> -vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) { - vector<T, D> v; +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)); @@ -51,425 +48,411 @@ vector<T, D> diff(T L, vector<T, D> v1, vector<T, D> v2) { return v; } -template <class T, int D, class state> -class spin { - public: - vector<T, D> x; - state s; +template <class T, int D, class S> class Spin { +public: + Vector<T, D> x; + S s; }; -template <class T, int D> -class euclidean { - private: - T L; - - public: - vector<T, D> t; - matrix<T, D> r; - euclidean(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; - } +template <class T, int D> class Euclidean { +private: + T L; + +public: + Vector<T, D> t; + Matrix<T, D> r; + + /** brief Euclidean - default constructor, constructs the identity + */ + Euclidean(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; } } + } - euclidean(T L, vector<T, D> t0, matrix<T, D> r0) : L(L) { - t = t0; - r = r0; - } + Euclidean(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) { + t = t0; + r = r0; + } - template <class state> - spin<T, D, state> act(const spin<T, D, state>& s) const { - spin<T, D, state> s_new; + 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; + 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; + for (unsigned i = 0; i < D; i++) { + s_new.x(i) = fmod(L + s_new.x(i), L); } - euclidean act(const euclidean& 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); - } + return s_new; + } - euclidean pnew(this->L, tnew, rnew); + Euclidean act(const Euclidean& x) const { + Vector<T, D> tnew = r * x.t + t; + Matrix<T, D> rnew = r * x.r; - return pnew; + for (unsigned i = 0; i < D; i++) { + tnew(i) = fmod(L + tnew(i), L); } - euclidean inverse() const { - vector<T, D> tnew = - r.transpose() * t; - matrix<T, D> rnew = r.transpose(); + Euclidean pnew(this->L, tnew, rnew); + + return pnew; + } - euclidean pnew(this->L, tnew, rnew); + Euclidean inverse() const { + Vector<T, D> tnew = -r.transpose() * t; + Matrix<T, D> rnew = r.transpose(); - return pnew; - } -}; + Euclidean pnew(this->L, tnew, rnew); -template <int D> -class dictionary { - private: - unsigned L; - std::vector<std::set<unsigned>> d; + return pnew; + } +}; - public: - dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {}; +template <class T, int D, class S> class Octree { +private: + T L; + unsigned N; + std::vector<std::set<Spin<T, D, S>*>> data; - template <class T> - unsigned dictionary_index(vector<T, D> x) const { - unsigned pos_ind = 0; +public: + Octree(T L, unsigned depth) : L(L), N(pow(2, depth)), data(pow(N, D)){}; - for (unsigned i = 0; i < D; i++) { - pos_ind += pow(L, i) * (unsigned)x(i); - }; + unsigned ind(Vector<T, D> x) const { + unsigned pos_ind = 0; - return pos_ind; + for (unsigned i = 0; i < D; i++) { + pos_ind += pow(N, i) * (unsigned)std::floor(N * x(i) / L); } - template <class T> - void record(vector<T, D> x, unsigned ind) { - d[dictionary_index<T>(x)].insert(ind); - }; + assert(pos_ind < data.size()); - template <class T> - void remove(vector<T, D> x, unsigned ind) { - d[dictionary_index<T>(x)].erase(ind); - }; + return pos_ind; + } - template <class T> - std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const { - return nearest_neighbors_of(dictionary_index<T>(x), depth, {}); - }; + void insert(Spin<T, D, S>* s) { data[ind(s->x)].insert(s); }; - std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, std::list<unsigned> ignores) const { - std::set<unsigned> ns = d[ind]; + void remove( Spin<T, D, S>* s) { data[ind(s->x)].erase(s); }; - if (depth > 0) { - for (unsigned i = 0; i < D; i++) { - if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k){return i==k;})) { - std::list<unsigned> ignores_new = ignores; - ignores_new.push_back(i); - for (signed j : {-1, 1}) { - unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1)); + 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; + }; - std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new); + 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); - for (unsigned guy : nns) { - ns.insert(guy); - } - } + 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); } } } - - return ns; - }; + } + } }; 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; - } +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); + 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); } + n++; + } - double avg2() const { - return total2 / (n - wait); - } + double avg() const { return total / (n - wait); } - std::vector<double> ρ() const { - double C0 = C.front() / (n - wait); - double avg2 = pow(total / (n - wait), 2); + double avg2() const { return total2 / (n - wait); } - std::vector<double> ρtmp; + std::vector<double> ρ() const { + double C0 = C.front() / (n - wait); + double avg2 = pow(total / (n - wait), 2); - for (double Ct : C) { - ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); - } + std::vector<double> ρtmp; - return ρtmp; + for (double Ct : C) { + ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); } - std::array<double, 2> τ() const { - double τtmp = 0.5; - unsigned M = 1; - double c = 8.0; + return ρtmp; + } - std::vector<double> ρ_tmp = this->ρ(); + std::array<double, 2> τ() const { + double τtmp = 0.5; + unsigned M = 1; + double c = 8.0; - while (c * τtmp > M && M < C.size()) { - τtmp += ρ_tmp[M]; - M++; - } + std::vector<double> ρ_tmp = this->ρ(); - return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)}; + while (c * τtmp > M && M < C.size()) { + τtmp += ρ_tmp[M]; + M++; } - double σ() const { - return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2)); - } + return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)}; + } - double serr() const { - return sqrt(this->σ()); - } + double σ() const { + return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2)); + } - unsigned num_added() const { - return n - wait; - } + double serr() const { return sqrt(this->σ()); } + + unsigned num_added() const { return n - wait; } }; -template <class U, int D, class state> -class model { - public: - U L; - euclidean<U, D> s0; - std::vector<spin<U, D, state>> s; - dictionary<D> dict; - std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors; - std::function<double(spin<U, D, state>, spin<U, D, state>, spin<U, D, state>)> Z; - std::function<double(spin<U, D, state>)> B; - std::vector<matrix<U, D>> mats; - std::vector<vector<U, D>> steps; - long double E; - quantity Eq; - quantity Cq; - - 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(sequences, new_level); +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 X_t&, double) {}; + virtual void plain_site_transformed(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const X_t&) {}; + + virtual void ghost_bond_visited(const system<R_t, X_t, G_t>&, const typename G_t::vertex& v, const X_t&, const X_t&, double) {}; + virtual void ghost_site_transformed(const system<R_t, X_t, G_t>&, const R_t&) {}; + + virtual void post_cluster(unsigned, unsigned, const system<R_t, X_t, G_t>&) {}; + }; + */ + +template <class U, int D, class S> class Model { +public: + U L; + Euclidean<U, D> 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; + std::vector<Matrix<U, D>> mats; + std::vector<Vector<U, D>> steps; + double E; + + 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(sequences, new_level); } + } - model(U L, std::function<double(spin<U, D, state>, spin<U, D, state>, spin<U, D, state>)> Z, - std::function<double(spin<U, D, state>)> B, - std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) : - L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) { - std::array<double, D> ini_sequence; - ini_sequence.fill(1); - std::list<std::array<double, D>> sequences; - sequences.push_back(ini_sequence); - - one_sequences(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; - } - } - } + 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) { + std::array<double, D> ini_sequence; + ini_sequence.fill(1); + std::list<std::array<double, D>> sequences; + sequences.push_back(ini_sequence); - mats.push_back(m); + one_sequences(sequences, D); - vector<U, D> v; - for (unsigned i = 0; i < D; i++) { - if (sequence[i] == 1) { - v(i) = 0; - } else { - v(i) = L / 2; - } + 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); - steps.push_back(v); + Vector<U, D> v; + for (unsigned i = 0; i < D; i++) { + if (sequence[i] == 1) { + v(i) = 0; + } else { + v(i) = L / 2; } + } + + steps.push_back(v); + } - 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; - } + 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); } } + mats.push_back(m); } } - - void update_energy() { - E = 0; - for (unsigned i = 0; i < s.size(); i++) { - for (unsigned j = 0; j < i; j++) { - // E -= Z(s[i], s[j]); - } - - E -= B(s0.inverse().act(s[i])); - } } + } - void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) { - unsigned cluster_size = 0; - std::uniform_real_distribution<double> dist(0.0, 1.0); + void update_energy() { + E = 0; + for (unsigned i = 0; i < s.size(); i++) { + for (unsigned j = 0; j < i; j++) { + E -= Z(s[i], s[j]); + } - std::queue<unsigned> queue; - queue.push(ind); + E -= B(s0.inverse().act(s[i])); + } + } - std::vector<bool> visited(s.size() + 1, false); + void step(double T, unsigned ind, Euclidean<U, D> r, std::mt19937& rng) { + unsigned cluster_size = 0; + std::uniform_real_distribution<double> dist(0.0, 1.0); - while (!queue.empty()) { - unsigned i = queue.front(); - queue.pop(); + std::queue<Spin<U, D, S>*> queue; + queue.push(&(s[ind])); - if (!visited[i]) { - visited[i] = true; + std::set<Spin<U, D, S>*> visited; - bool we_are_ghost = i == s.size(); + while (!queue.empty()) { + Spin<U, D, S>* si = queue.front(); + queue.pop(); - spin<U, D, state> si_new; - euclidean<U, D> s0_new(L); + if (!visited.contains(si)) { + visited.insert(si); - if (we_are_ghost) { - s0_new = r.act(s0); - } else { - si_new = r.act(s[i]); + if (si == NULL) { + Euclidean<U, D> 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); + double p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T); + if (dist(rng) < p) { + queue.push(&ss); + } + s0 = s0_new; } - - for (unsigned j : neighbors(*this, i, si_new)) { - if (j != i) { + } 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); + + 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) { + if (sj != si) { double p; - bool neighbor_is_ghost = j == s.size(); - - if (we_are_ghost || neighbor_is_ghost) { - spin<U, D, state> s0s_old, s0s_new; - unsigned non_ghost; - - if (neighbor_is_ghost) { - non_ghost = i; - s0s_old = s0.inverse().act(s[i]); - s0s_new = s0.inverse().act(si_new); - } else { - non_ghost = j; - s0s_old = s0.inverse().act(s[j]); - s0s_new = s0_new.inverse().act(s[j]); - } - + if (sj == NULL) { + Spin<U, D, S> s0s_old = s0.inverse().act(*si); + Spin<U, D, S> s0s_new = s0.inverse().act(si_new); p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T); } else { - p = Z(s[i], s[j], si_new); + p = 1.0 - exp(-(Z(si_new, *sj) - Z(*si, *sj)) / T); } - if (dist(rng) < p) { - queue.push(j); + queue.push(sj); } } } - - if (we_are_ghost) { - s0 = s0_new; - } else { - dict.remove(s[i].x, i); - s[i] = si_new; - dict.record(s[i].x, i); - cluster_size++; - } + dict.remove(si); + *si = si_new; + dict.insert(si); + cluster_size++; } } - - Cq.add(cluster_size); } + } - void wolff(double T, unsigned N, std::mt19937& rng) { - 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); - - for (unsigned i = 0; i < N; i++) { - vector<U, D> t; - matrix<U, D> m; - unsigned flip = coin(rng); - if (flip < mats.size()) { - for (unsigned j = 0; j < D; j++) { - t(j) = (U)t_dist(rng); - } - m = mats[flip]; - } else { - for (unsigned j = 0; j < D; j++) { - for (unsigned k = 0; k < D; k++) { - if (j == k) { - m(j, k) = 1; - } else { - m(j, k) = 0; - } + void wolff(double T, unsigned N, std::mt19937& rng) { + 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); + + for (unsigned i = 0; i < N; i++) { + Vector<U, D> t; + Matrix<U, D> m; + unsigned flip = coin(rng); + if (flip < mats.size()) { + for (unsigned j = 0; j < D; j++) { + t(j) = (U)t_dist(rng); + } + m = mats[flip]; + } else { + for (unsigned j = 0; j < D; j++) { + for (unsigned k = 0; k < D; k++) { + if (j == k) { + m(j, k) = 1; + } else { + m(j, k) = 0; } } - - t = steps[flip - mats.size()]; } - euclidean<U, D> g(L, t, m); + t = steps[flip - mats.size()]; + } - this->step(T, ind_dist(rng), g, rng); + Euclidean<U, D> g(L, t, m); - this->update_energy(); - Eq.add(E); - } + this->step(T, ind_dist(rng), g, rng); + + this->update_energy(); } + } }; - diff --git a/spheres.cpp b/spheres.cpp index b0f4fb4..633a26a 100644 --- a/spheres.cpp +++ b/spheres.cpp @@ -34,65 +34,49 @@ int main(int argc, char* argv[]) { } } - std::function<double(spin<double, D, double>, spin<double, D, double>, spin<double, D, double>)> Z = - [L] (spin<double, D, double> s1, spin<double, D, double> s2, spin<double, D, double> s1_new) -> double { - vector<double, D> diff_old = diff(L, s1.x, s2.x); - vector<double, D> diff_new = diff(L, s1_new.x, s2.x); + std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z = + [L] (const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double { + Vector<double, D> d = diff(L, s1.x, s2.x); double rad_sum = pow(s1.s + s2.s, 2); - bool old_overlap = diff_old.transpose() * diff_old < rad_sum; - bool new_overlap = diff_new.transpose() * diff_new < rad_sum; + bool overlap = d.transpose() * d < rad_sum; - if (new_overlap) { - return 1.0; + if (overlap) { + return -1e8; } else { - return 0.0; + return 0; } }; - std::function<double(spin<double, D, double>)> B = - [L, H] (spin<double, D, double> s) -> double { + std::function<double(Spin<double, D, double>)> B = + [L, H] (Spin<double, D, double> s) -> double { return H * sin(2 * M_PI * 3 * s.x(0) / L); }; - std::function<std::set<unsigned>(model<double, D, double>&, unsigned, spin<double, D, double>)> neighbors = - [] (model<double, D, double>& m, unsigned i0, spin<double, D, double> s1) -> std::set<unsigned> { - std::set<unsigned> nn; - if (i0 < m.s.size()) { - std::set<unsigned> n_old = m.dict.neighbors(m.s[i0].x, 1); - std::set<unsigned> n_new = m.dict.neighbors(s1.x, 1);; - nn.insert(n_old.begin(), n_old.end()); - nn.insert(n_new.begin(), n_new.end()); - nn.insert(m.s.size()); - } else { - for (unsigned i = 0; i < m.s.size(); i++) { - nn.insert(i); - } - } - return nn; - }; - - model<double, D, double> sphere(L, Z, B, neighbors); + Model<double, D, double> sphere(L, Z, B, std::floor(log2(L)), 2); randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; std::uniform_real_distribution<double> dist(0.0, L); + sphere.s.reserve(n); + for (unsigned i = 0; i < n; i++) { - vector<double, D> pos = {dist(rng), dist(rng)}; + Vector<double, D> pos = {dist(rng), dist(rng)}; sphere.s.push_back({pos, 0.5}); - sphere.dict.record<double>(pos, i); + sphere.dict.insert(&sphere.s.back()); } + sphere.wolff(T, N, rng); std::ofstream snapfile; snapfile.open("sphere_snap.dat"); - for (spin<double, D, double> s : sphere.s) { - spin<double, D, double> rs = sphere.s0.inverse().act(s); + for (Spin<double, D, double> s : sphere.s) { + Spin<double, D, double> rs = sphere.s0.inverse().act(s); snapfile << rs.x.transpose() << "\n"; } |