diff options
| -rw-r--r-- | biroli-mezard.cpp | 161 | ||||
| -rw-r--r-- | ciamarra.cpp | 489 | ||||
| -rw-r--r-- | glass.hpp | 677 |
3 files changed, 600 insertions, 727 deletions
diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp index 7e200b9..e2fb71a 100644 --- a/biroli-mezard.cpp +++ b/biroli-mezard.cpp @@ -5,6 +5,7 @@ #include <queue> #include <vector> #include <chrono> +#include <unordered_set> #include <eigen3/Eigen/Dense> @@ -152,6 +153,7 @@ public: unsigned occupiedNeighbors; unsigned maximumNeighbors; bool marked; + bool visited; Vertex() { occupiedNeighbors = 0; @@ -181,6 +183,7 @@ template <int D> class System { public: const unsigned L; std::vector<unsigned> N; + std::unordered_set<Vertex<D>*> occupants; std::vector<Vertex<D>> vertices; Transformation<D> orientation; @@ -322,6 +325,9 @@ public: } N[t]++; N[0]--; + if (t > 0) { + occupants.insert(&v); + } return true; } else { return false; @@ -332,6 +338,9 @@ public: if (v.empty()) { return false; } else { + if (v.maximumNeighbors > 0) { + occupants.erase(&v); + } N[v.maximumNeighbors]--; N[0]++; v.maximumNeighbors = Empty; @@ -418,30 +427,32 @@ public: return exchange(v1, v2); } - void sweepGrandCanonical(double z, Rng& r) { - for (unsigned i = 0; i < size(); i++) { - if (0.5 < r.uniform(0.0, 1.0)) { - double pIns = size() * z / (occupancy() + 1); + const std::array<double, 3> ratios = {0.1, 0.5, 0.4}; - if (pIns > r.uniform(0.0, 1.0)) { - while (true) { - Vertex<D>& v = r.pick(vertices); - if (v.empty()) { - insert(v, r.pick({1, 2, 2, 2, 2, 2, 3, 3, 3, 3})); - break; - } - } - } - } else { + void stepGrandCanonical(double z, Rng& r) { + if (1.0 / 11.0 < r.uniform(0.0, 1.0)) { + unsigned t = r.pick({1, 2, 2, 2, 2, 2, 3, 3, 3, 3}); + double pIns = size() * z / (occupancy() + 1); - double pDel = density() / z; + if (pIns > r.uniform(0.0, 1.0)) { + Vertex<D>& v = r.pick(vertices); + insert(v, t); + } + } else { + double pDel = density() / z; - if (pDel > r.uniform(0.0, 1.0)) { - remove(r.pick(vertices)); - } + if (pDel > r.uniform(0.0, 1.0)) { + Vertex<D>* v = r.pick(occupants); + remove(*v); } + } + } - randomLocalExchange(r); + void sweepGrandCanonical(double z, Rng& r) { + for (unsigned i = 0; i < size(); i++) { + stepGrandCanonical(z, r); + + randomExchange(r); } } @@ -457,6 +468,71 @@ public: } } + Vertex<D>& apply(const Transformation<D>& R, const Vertex<D>& v) { + return vertices[vectorToIndex(R.apply(v.position))]; + } + + bool eventChain(const Transformation<D>& R, Vertex<D>& v0) { + unsigned n = 0; + + if (!v0.empty()) { + Vertex<D>& v0New = apply(R, v0);; + unsigned t0 = v0.maximumNeighbors; + + std::queue<std::pair<std::reference_wrapper<Vertex<D>>, unsigned>> q; + + q.push({v0New, t0}); + remove(v0); + + while (!q.empty()) { + auto [vr, t] = q.front(); + Vertex<D>& v = vr; + q.pop(); + + if (!v.empty()) { + q.push({apply(R, v), v.maximumNeighbors}); + remove(v); + } + + insert(v, t, true); + v.marked = true; + + for (Vertex<D>& vn : overlaps({v})) { + if (!vn.marked) { + q.push({apply(R, vn), vn.maximumNeighbors}); + remove(vn); + vn.marked = true; + } else { + /* If a neighbor has already been moved but is still + * frustrated, it may be due to one of its neighbors! + */ + if (&vn != &v) { + for (Vertex<D>& vnn : overlaps(vn)) { + if (!vnn.marked) { + q.push({apply(R, vnn), vnn.maximumNeighbors}); + remove(vnn); + vnn.marked = true; + } + } + } + } + } + } + } + + for (Vertex<D>& v : vertices) { + v.marked = false; + } + + return true; + } + + bool randomEventChain(Rng& r) { + Transformation<D> R(L); + R.v[r.uniform((unsigned)0, (unsigned)D-1)] = r.pick({-1, 1}); + return eventChain(R, r.pick(vertices)); + } + unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0) { unsigned n = 0; @@ -546,30 +622,35 @@ int main() { unsigned L = 15; unsigned Nmin = 2e2; unsigned Nmax = 2e5; - double Tmin = 0.04; - double Tmax = 0.2; - double δT = 0.02; + double μmin = -3; + double μmax = 10; + double dμ = 0.05; System<D> s(L, 3); Rng r; - double z = exp(1 / 0.15); + for (double μ = μmin; μ <= μmax; μ += dμ) { + double z = exp(μ); + for (unsigned i = 0; i < 1e4; i++) { + for (unsigned j = 0; j < s.size(); j++) { + s.stepGrandCanonical(z, r); +// s.randomEventChain(r); + s.randomExchange(r); + } - if (!s.compatible()) { - std::cerr << "Storted incompatible!" << std::endl; - return 1; - } + if (!s.compatible()) { + std::cerr << "Not compatible!" << std::endl; + return 1; + } - while (s.density() < 0.56) { - s.sweepGrandCanonical(z, r); - } + } - if (!s.compatible()) { - std::cerr << "Not compatible!" << std::endl; - return 1; + std::cout << μ << " " << s.density() << std::endl; } + + /* std::cerr << "Found state with appropriate density." << std::endl; System<D> s0 = s; @@ -586,22 +667,19 @@ int main() { auto start = std::chrono::high_resolution_clock::now(); while (nC / s.size() < 1e5) { if (n < 20 * log(i + 1)) { + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = duration_cast<std::chrono::microseconds>(stop - start); n++; - std::cout << nC / s.size() << " " << s.overlap(s0) << std::endl; + std::cout << duration.count() << " " << s.overlap(s0) << std::endl; } - //unsigned nn = s.flipCluster(Transformation<D>(L, ms, r), r.pick(s.vertices)); - //nC += nn; +// unsigned nn = s.flipCluster(Transformation<D>(L, ms, r), r.pick(s.vertices)); +// nC += nn; // clusterDist[nn]++; //s.sweepLocal(r); nC += s.size(); s.sweepSwap(r); i++; } - auto stop = std::chrono::high_resolution_clock::now(); - - auto duration = duration_cast<std::chrono::microseconds>(stop - start); - - std::cerr << duration.count() << std::endl; if (!s.compatible()) { std::cerr << "Not compatible!" << std::endl; @@ -618,6 +696,7 @@ int main() { file << i << " "; } file.close(); + */ return 0; } diff --git a/ciamarra.cpp b/ciamarra.cpp index 6b932bf..f2f44c1 100644 --- a/ciamarra.cpp +++ b/ciamarra.cpp @@ -1,10 +1,477 @@ +#include <eigen3/Eigen/Dense> +#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h> #include <iostream> +#include <list> +#include <vector> +#include <queue> -#include "glass.hpp" +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" +using Rng = randutils::random_generator<pcg32>; -void print(const CiamarraSystem<2>& s) { - for (const Vertex<2, CiamarraState<2>>& v : s.vertices) { +template <int D> using Vector = Eigen::Matrix<int, D, 1>; +template <int D> using Matrix = Eigen::Matrix<int, D, D>; + +int iPow(int x, unsigned p) { + if (p == 0) + return 1; + if (p == 1) + return x; + + int tmp = iPow(x, p / 2); + if (p % 2 == 0) + return tmp * tmp; + else + return x * tmp * tmp; +} + +unsigned mod(signed a, unsigned b) { + return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b; +} + +template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) { + Vector<D> u; + for (unsigned i = 0; i < D; i++) { + u(i) = mod(v(i), b); + } + return u; +} + +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 <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() { + std::vector<Matrix<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_back(); // don't want the identity matrix! + + for (std::array<double, D> sequence : sequences) { + Matrix<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<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 if (k == l && (k != i && k != j)) { + m(k, l) = 1; + } else { + m(k, l) = 0; + } + } + } + mats.push_back(m); + } + } + } + + return mats; +} + +template <unsigned D> class State : public Vector<D> { +public: + State() : Vector<D>(Vector<D>::Zero()) {} + + State(const Vector<D>& v) : Vector<D>(v) {} + + State(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; } + + State(Rng& r) : State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {} + + bool isEmpty() const { return this->squaredNorm() == 0; } + + void remove() { this->setZero(); } + + State<D> flip() const { + State<D> s; + for (unsigned i = 0; i < D; i++) { + s(i) = -this->operator()(i); + } + return s; + } +}; + +template <unsigned D> class Transformation { + public: + unsigned L; + Matrix<D> m; + Vector<D> v; + + Transformation(unsigned L) : L(L) { + m.setIdentity(); + v.setZero(); + } + + Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {} + + Transformation(unsigned L, const std::vector<Matrix<D>>& ms, Rng& r) : m(r.pick(ms)), L(L) { + for (unsigned i = 0; i < D; i++) { + v[i] = r.uniform((unsigned)0, L - 1); + } + + v = v - m * v; + } + + Vector<D> apply(const Vector<D>& x) const { + return mod<D>(v + m * x, L); + } + + Transformation<D> apply(const Transformation<D>& t) const { + Transformation<D> tNew(L); + + tNew.m = m * t.m; + tNew.v = apply(t.v); + + return tNew; + } + + State<D> apply(const State<D>& x) const { + return State<D>(m * x); + } + + Transformation<D> inverse() const { + return Transformation<D>(L, m.transpose(), -m.transpose() * v); + } +}; + +template <unsigned D> class HalfEdge; + +template <unsigned D> class Vertex { +public: + Vector<D> position; + State<D> state; + std::vector<HalfEdge<D>> adjacentEdges; + bool marked; + + bool isEmpty() const { return state.isEmpty(); } +}; + +template <unsigned D> class HalfEdge { +public: + Vertex<D>& neighbor; + Vector<D> Δx; + + HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {} +}; + +template <unsigned D> class System { +public: + const unsigned L; + unsigned N; + std::vector<Vertex<D>> vertices; + Transformation<D> orientation; + + unsigned vectorToIndex(const Vector<D>& x) const { + unsigned i = 0; + for (unsigned d = 0; d < D; d++) { + i += x[d] * iPow(L, d); + } + return i; + } + + Vector<D> indexToVector(unsigned i) const { + Vector<D> x; + for (unsigned d = 0; d < D; d++) { + x[d] = (i / iPow(L, d)) % L; + } + return x; + } + + System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) { + for (unsigned i = 0; i < iPow(L, D); i++) { + vertices[i].position = indexToVector(i); + vertices[i].adjacentEdges.reserve(2 * D); + vertices[i].marked = false; + } + + for (unsigned d = 0; d < D; d++) { + Vector<D> Δx = Vector<D>::Zero(); + Δx[d] = 1; + for (signed i = 0; i < iPow(L, D); i++) { + unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1)); + vertices[i].adjacentEdges.push_back(HalfEdge<D>(vertices[j], Δx)); + vertices[j].adjacentEdges.push_back(HalfEdge<D>(vertices[i], -Δx)); + } + } + } + + std::list<std::reference_wrapper<Vertex<D>>> overlaps(Vertex<D>& v, const State<D>& s, + bool excludeSelf = false) { + std::list<std::reference_wrapper<Vertex<D>>> o; + + if (s.isEmpty()) { + return o; + } + + if (!v.isEmpty() && !excludeSelf) { + o.push_back(v); + } + + for (const HalfEdge<D>& e : v.adjacentEdges) { + if (!e.neighbor.isEmpty()) { + if (s.dot(e.Δx) == 1 || e.neighbor.state.dot(e.Δx) == -1) { + o.push_back(e.neighbor); + } + } + } + + return o; + } + + bool insert(Vertex<D>& v, const State<D>& s) { + if (overlaps(v, s).empty()) { + v.state = s; + N++; + return true; + } else { + return false; + } + } + + bool tryDeletion(Vertex<D>& v) { + if (v.isEmpty()) { + return false; + } else { + v.state.remove(); + N--; + return true; + } + } + + bool tryRandomMove(Rng& r) { + Vertex<D>& v = r.pick(vertices); + State<D> oldState = v.state; + + if (!tryDeletion(v)) { + return false; + } + + if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) { + for (HalfEdge<D>& e : v.adjacentEdges) { + if (1 == e.Δx.dot(oldState)) { + if (insert(e.neighbor, oldState.flip())) { + return true; + } + break; + } + } + } else { + State<D> newState(r); + while (newState.dot(oldState) == 1) { + newState = State<D>(r); + } + if (insert(v, newState)) { + return true; + } + } + v.state = oldState; + N++; + return false; + } + + bool trySwap(Vertex<D>& v1, Vertex<D>& v2) { + if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) { + std::swap(v1.state, v2.state); + return true; + } else { + return false; + } + } + + bool tryRandomSwap(Rng& r) { + Vertex<D>& v1 = r.pick(vertices); + Vertex<D>& v2 = r.pick(vertices); + + return trySwap(v1, v2); + } + + void setGroundState() { + N = 0; + + for (Vertex<D>& v : vertices) { + unsigned a = 0; + for (unsigned d = 0; d < D; d++) { + a += (d + 1) * v.position(d); + } + a %= 2 * D + 1; + + v.state.setZero() = Vector<D>::Zero(); + + if (0 < a && a <= D) { + v.state(a - 1) = -1; + N++; + } else if (D < a) { + v.state(2 * D - a) = 1; + N++; + } + } + } + + bool compatible() { + for (Vertex<D>& v : vertices) { + if (overlaps(v, v.state, true).size() > 0) { + return false; + } + } + + return true; + } + + double density() const { return N / pow(L, D); } + + unsigned maxOccupation() const { + // return (2 * D * iPow(L, D)) / (2 * D + 1); + return iPow(L, D); + } + + void sweepGrandCanonical(double z, Rng& r) { + for (unsigned i = 0; i < iPow(L, D); i++) { + if (0.5 < r.uniform(0.0, 1.0)) { + double pIns = maxOccupation() * z / (N + 1); + + if (pIns > r.uniform(0.0, 1.0)) { + while (true) { + Vertex<D>& v = r.pick(vertices); + if (v.isEmpty()) { + insert(v, State<D>(r)); + break; + } + } + } + } else { + + double pDel = N / (z * maxOccupation()); + + if (pDel > r.uniform(0.0, 1.0)) { + tryDeletion(r.pick(vertices)); + } + } + + tryRandomMove(r); + } + } + + void sweepLocal(Rng& r) { + for (unsigned i = 0; i < iPow(L, D); i++) { + tryRandomMove(r); + } + } + + void sweepSwap(Rng& r) { + for (unsigned i = 0; i < iPow(L, D); i++) { + tryRandomSwap(r); + } + } + + unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0, Rng& r, bool dry = false) { + std::queue<std::reference_wrapper<Vertex<D>>> q; + q.push(v0); + + unsigned n = 0; + + while (!q.empty()) { + Vertex<D>& v = q.front(); + q.pop(); + + if (!v.marked) { + Vector<D> xNew = R.apply(v.position); + Vertex<D>& vNew = vertices[vectorToIndex(xNew)]; + + v.marked = true; + vNew.marked = true; + + State<D> s = R.apply(v.state); + State<D> sNew = R.apply(vNew.state); + + std::list<std::reference_wrapper<Vertex<D>>> overlaps1 = overlaps(vNew, s, true); + std::list<std::reference_wrapper<Vertex<D>>> overlaps2 = overlaps(v, sNew, true); + overlaps1.splice(overlaps1.begin(), overlaps2); + + for (Vertex<D>& vn : overlaps1) { + if (!vn.marked) { + q.push(vn); + } + } + + if (!dry) { + v.state = sNew; + vNew.state = s; + } + + n += 1; + } + } + + return n; + } + + void swendsenWang(const Transformation<D>& R, Rng& r) { + for (Vertex<D>& v : vertices) { + if (!v.marked) { + bool dry = 0.5 < r.uniform(0.0, 1.0); + unsigned n = flipCluster(R, v, r, dry); + if (n > pow(L, D) / 4 && !dry) { + orientation = R.apply(orientation); + } + } + } + + for (Vertex<D>& v : vertices) { + v.marked = false; + } + } + + int overlap(const System<D>& s) const { + int o = 0; + + for (unsigned i = 0; i < vertices.size(); i++) { + State<D> s2 = orientation.apply(s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].state); + o += vertices[i].state.dot(s2); + } + + return o; + } +}; + +void print(const System<2>& s) { + for (const Vertex<2>& v : s.vertices) { if (v.state(0) == 1 && v.state(1) == 0) { std::cerr << "▶"; } else if (v.state(0) == -1 && v.state(1) == 0) { @@ -42,7 +509,7 @@ int main() { double Tmax = 0.2; double δT = 0.02; - CiamarraSystem<D> s(L); + System<D> s(L); Rng r; @@ -58,11 +525,8 @@ int main() { } std::cerr << "Found state with appropriate density." << std::endl; - if (!s.compatible()) { - std::cerr <<"whoops!" << std::endl; - } - CiamarraSystem<D> s0 = s; + System<D> s0 = s; std::vector<Matrix<D>> ms = generateTorusMatrices<D>(); @@ -72,9 +536,16 @@ int main() { n++; std::cout << i << " " << s.overlap(s0) << std::endl; } + Matrix<D> m = r.pick(ms); + Vertex<D>& v = r.pick(s.vertices); + unsigned nC = s.flipCluster(Transformation<D>(L, m, v.position - m * v.position), v, r); + std::cerr << nC << std::endl; + for (Vertex<D>& v : s.vertices) { + v.marked = false; + } // s.sweepLocal(r); // s.sweepSwap(r); - s.swendsenWang(Transformation<D>(L, ms, r), r); +// s.swendsenWang(Transformation<D>(L, ms, r), r); } return 0; diff --git a/glass.hpp b/glass.hpp deleted file mode 100644 index b35964b..0000000 --- a/glass.hpp +++ /dev/null @@ -1,677 +0,0 @@ -#include <limits> -#include <list> -#include <queue> -#include <vector> - -#include <eigen3/Eigen/Dense> - -#include "pcg-cpp/include/pcg_random.hpp" -#include "randutils/randutils.hpp" - -using Rng = randutils::random_generator<pcg32>; - -template <int D> using Vector = Eigen::Matrix<int, D, 1>; -template <int D> using Matrix = Eigen::Matrix<int, D, D>; - -int iPow(int x, unsigned p) { - if (p == 0) - return 1; - if (p == 1) - return x; - - int tmp = iPow(x, p / 2); - if (p % 2 == 0) - return tmp * tmp; - else - return x * tmp * tmp; -} - -unsigned mod(signed a, unsigned b) { return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b; } - -template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) { - Vector<D> u; - for (unsigned i = 0; i < D; i++) { - u(i) = mod(v(i), b); - } - return u; -} - -template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) { - if (level > 0) { - unsigned new_level = level - 1; - 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 <int D> std::vector<Matrix<D>> generateTorusMatrices() { - std::vector<Matrix<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_back(); // don't want the identity matrix! - - for (std::array<double, D> sequence : sequences) { - Matrix<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<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 if (k == l && (k != i && k != j)) { - m(k, l) = 1; - } else { - m(k, l) = 0; - } - } - } - mats.push_back(m); - } - } - } - - return mats; -} - -template <int D> class CiamarraState : public Vector<D> { -public: - CiamarraState() : Vector<D>(Vector<D>::Zero()) {} - - CiamarraState(const Vector<D>& v) : Vector<D>(v) {} - - CiamarraState(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; } - - CiamarraState(Rng& r) : CiamarraState(r.uniform(0, D - 1), r.pick({-1, 1})) {} - - bool isEmpty() const { return this->squaredNorm() == 0; } - - void remove() { this->setZero(); } - - CiamarraState<D> flip() const { - CiamarraState<D> s; - for (unsigned i = 0; i < D; i++) { - s(i) = -this->operator()(i); - } - return s; - } -}; - -class BiroliState { -public: - unsigned type; - unsigned occupiedNeighbors; - - BiroliState() { - type = std::numeric_limits<unsigned>::max(); - occupiedNeighbors = 0; - } - - BiroliState(unsigned t, unsigned o) { - type = t; - occupiedNeighbors = o; - } - - BiroliState(Rng& r) { - type = r.pick({1,2,2,2,2,2,3,3,3,3}); - } - - bool isEmpty() const { return type == std::numeric_limits<unsigned>::max(); } - - void remove() { type = std::numeric_limits<unsigned>::max(); }; -}; - -template <int D> class Transformation { -public: - unsigned L; - Matrix<D> m; - Vector<D> v; - - Transformation(unsigned L) : L(L) { - m.setIdentity(); - v.setZero(); - } - - Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {} - - Transformation(unsigned L, const std::vector<Matrix<D>>& ms, Rng& r) : m(r.pick(ms)), L(L) { - for (unsigned i = 0; i < D; i++) { - v[i] = r.uniform((unsigned)0, L - 1); - } - - v = v - m * v; - } - - Transformation<D> inverse() const { - return Transformation<D>(L, m.transpose(), -m.transpose() * v); - } - - Vector<D> apply(const Vector<D>& x) const { return mod<D>(v + m * x, L); } - - Transformation<D> apply(const Transformation<D>& t) const { - Transformation<D> tNew(L); - - tNew.m = m * t.m; - tNew.v = apply(t.v); - - return tNew; - } - - CiamarraState<D> apply(const CiamarraState<D>& s) const { return CiamarraState<D>(m * s); } - - BiroliState apply(const BiroliState& s) const { return s; } -}; - -template <int D, typename State> class HalfEdge; - -template <int D, typename State> class Vertex { -public: - Vector<D> position; - std::vector<HalfEdge<D, State>> adjacentEdges; - State state; - bool marked; - - bool isEmpty() const { return state.isEmpty(); } -}; - -template <int D, class State> class HalfEdge { -public: - Vertex<D, State>& neighbor; - Vector<D> Δx; - - HalfEdge(Vertex<D, State>& n, const Vector<D>& d) : neighbor(n), Δx(d) {} -}; - -template <int D, class State> class System { -public: - const unsigned L; - unsigned N; - std::vector<Vertex<D, State>> vertices; - Transformation<D> orientation; - - unsigned vectorToIndex(const Vector<D>& x) const { - unsigned i = 0; - for (unsigned d = 0; d < D; d++) { - i += x[d] * iPow(L, d); - } - return i; - } - - Vector<D> indexToVector(unsigned i) const { - Vector<D> x; - for (unsigned d = 0; d < D; d++) { - x[d] = (i / iPow(L, d)) % L; - } - return x; - } - - System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) { - for (unsigned i = 0; i < iPow(L, D); i++) { - vertices[i].position = indexToVector(i); - vertices[i].adjacentEdges.reserve(2 * D); - vertices[i].marked = false; - } - - for (unsigned d = 0; d < D; d++) { - Vector<D> Δx = Vector<D>::Zero(); - Δx[d] = 1; - for (signed i = 0; i < iPow(L, D); i++) { - unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1)); - vertices[i].adjacentEdges.push_back(HalfEdge<D, State>(vertices[j], Δx)); - vertices[j].adjacentEdges.push_back(HalfEdge<D, State>(vertices[i], -Δx)); - } - } - } - - virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&) { return {}; }; - - virtual bool insert(Vertex<D, State>& v, const State& s, bool force = false) { return false; }; - - virtual bool remove(Vertex<D, State>& v) { return false;}; - - virtual bool randomMove(Rng& r) { return false;}; - - virtual void swap(Vertex<D, State>& v1, Vertex<D, State>& v2) {}; - - virtual bool exchange(Vertex<D, State>& v1, Vertex<D, State>& v2) { return false;}; - - bool randomExchange(Rng& r) { - Vertex<D, State>& v1 = r.pick(vertices); - Vertex<D, State>& v2 = r.pick(vertices); - - return exchange(v1, v2); - } - - bool compatible() { - for (Vertex<D, State>& v : vertices) { - if (overlaps(v).size() > 0) { - return false; - } - } - - return true; - } - - double density() const { return N / pow(L, D); } - - unsigned size() const { return vertices.size(); } - - void sweepGrandCanonical(double z, Rng& r) { - for (unsigned i = 0; i < iPow(L, D); i++) { - if (0.5 < r.uniform(0.0, 1.0)) { - double pIns = size() * z / (N + 1); - - if (pIns > r.uniform(0.0, 1.0)) { - while (true) { - Vertex<D, State>& v = r.pick(vertices); - if (v.isEmpty()) { - insert(v, State(r)); - break; - } - } - } - } else { - - double pDel = N / (z * size()); - - if (pDel > r.uniform(0.0, 1.0)) { - remove(r.pick(vertices)); - } - } - - randomMove(r); - } - } - - void sweepLocal(Rng& r) { - for (unsigned i = 0; i < iPow(L, D); i++) { - randomMove(r); - } - } - - void sweepSwap(Rng& r) { - for (unsigned i = 0; i < iPow(L, D); i++) { - randomExchange(r); - } - } - - unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, bool dry = false) { - unsigned n = 0; - - Vertex<D, State>& v0New = vertices[vectorToIndex(R.apply(v0.position))]; - - if (&v0New != &v0) { - std::queue<std::reference_wrapper<Vertex<D, State>>> q; - - v0.marked = true; - v0New.marked = true; - - swap(v0, v0New); - - for (Vertex<D, State>& vn : overlaps(v0)) { - if (!vn.marked) { - q.push(vn); - } - } - for (Vertex<D, State>& vn : overlaps(v0New)) { - if (!vn.marked) { - q.push(vn); - } - } - - while (!q.empty()) { - Vertex<D, State>& v = q.front(); - q.pop(); - - if (!v.marked && !overlaps(v).empty()) { - v.marked = true; - Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.position))]; - - if (&vNew != &v) { - vNew.marked = true; - - swap(v, vNew); - - for (Vertex<D, State>& vn : overlaps(v)) { - if (!vn.marked) { - q.push(vn); - } - } - for (Vertex<D, State>& vn : overlaps(vNew)) { - if (!vn.marked) { - q.push(vn); - } - } - - n += 2; - } else { - n += 1; - } - } - if (q.empty()) { - for (Vertex<D, State>& vv : vertices) { - if (!vv.marked && !overlaps(vv).empty()) { -// std::cerr << "Found bad state at end" << std::endl; - q.push(vv); - } - } - } - } - - - } - - if (n > size() / 4) { - orientation = R.apply(orientation); - } - - for (Vertex<D, State>& v : vertices) { - v.marked = false; - } - - return n; - } - - void swendsenWang(const Transformation<D>& R, Rng& r) { - for (Vertex<D, State>& v : vertices) { - if (!v.marked) { - bool dry = 0.5 < r.uniform(0.0, 1.0); - unsigned n = flipCluster(R, v, dry); - if (n > size() / 4 && !dry) { - orientation = R.apply(orientation); - } - } - } - - for (Vertex<D, State>& v : vertices) { - v.marked = false; - } - } - - virtual int overlap(const System<D, State>& s) const { return 0; }; -}; - -template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> { - public: - - using System<D, CiamarraState<D>>::System; - std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>> - overlaps(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s, - bool excludeSelf = false) override { - std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>> o; - - if (s.isEmpty()) { - return o; - } - - if (!v.isEmpty() && !excludeSelf) { - o.push_back(v); - } - - for (const HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) { - if (!e.neighbor.isEmpty()) { - if (s.dot(e.Δx) == 1 || e.neighbor.state.dot(e.Δx) == -1) { - o.push_back(e.neighbor); - } - } - } - - return o; - } - - bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s, bool force = false) override { - if (force || overlaps(v, s).empty()) { - v.state = s; - this->N++; - return true; - } else { - return false; - } - } - - bool remove(Vertex<D, CiamarraState<D>>& v) override { - if (v.isEmpty()) { - return false; - } else { - v.state.remove(); - this->N--; - return true; - } - } - - bool randomMove(Rng& r) override { - Vertex<D, CiamarraState<D>>& v = r.pick(this->vertices); - CiamarraState<D> oldState = v.state; - - if (!remove(v)) { - return false; - } - - if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) { - for (HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) { - if (1 == e.Δx.dot(oldState)) { - if (insert(e.neighbor, oldState.flip())) { - return true; - } - break; - } - } - } else { - CiamarraState<D> newState(r); - while (newState.dot(oldState) == 1) { - newState = CiamarraState<D>(r); - } - if (insert(v, newState)) { - return true; - } - } - v.state = oldState; - this->N++; - return false; - } - - void swap(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override { - std::swap(v1.state, v2.state); - } - - bool exchange(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override { - if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) { - swap(v1, v2); - return true; - } else { - return false; - } - } - - void setGroundState() { - this->N = 0; - - for (Vertex<D, CiamarraState<D>>& v : this->vertices) { - unsigned a = 0; - for (unsigned d = 0; d < D; d++) { - a += (d + 1) * v.position(d); - } - a %= 2 * D + 1; - - v.state.setZero() = Vector<D>::Zero(); - - if (0 < a && a <= D) { - v.state(a - 1) = -1; - this->N++; - } else if (D < a) { - v.state(2 * D - a) = 1; - this->N++; - } - } - } - - int overlap(const System<D, CiamarraState<D>>& s) const override { - int o = 0; - - for (unsigned i = 0; i < this->vertices.size(); i++) { - CiamarraState<D> s2 = this->orientation.apply( - s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state); - o += this->vertices[i].state.dot(s2); - } - - return o; - } -}; - -template <int D> class BiroliSystem : public System<D, BiroliState> { -public: - using System<D, BiroliState>::System; - - bool canReplaceWith(const Vertex<D, BiroliState>& v, const BiroliState& s) { - if (s.isEmpty()) { - return true; - } - if (s.type < v.state.occupiedNeighbors) { - return false; - } else { - int Δn = 0; - if (v.isEmpty()) { - Δn += 1; - } - for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) { - if (e.neighbor.state.type < e.neighbor.state.occupiedNeighbors + Δn) { - return false; - } - } - } - return true; - } - - std::list<std::reference_wrapper<Vertex<D, BiroliState>>> - overlaps(Vertex<D, BiroliState>& v) override { - std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o; - - if (v.isEmpty()) { // an empty site cannot be frustrating anyone - return o; - } - - bool selfFrustrated = v.state.occupiedNeighbors > v.state.type; - bool anyNeighborFrustrated = false; - - for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) { - if (!e.neighbor.isEmpty()) { - bool thisNeighborFrustrated = e.neighbor.state.type < e.neighbor.state.occupiedNeighbors; - - if (thisNeighborFrustrated) { - anyNeighborFrustrated = true; - } - - if (selfFrustrated || thisNeighborFrustrated) { - o.push_back(e.neighbor); - } - } - } - - if (selfFrustrated || anyNeighborFrustrated) { - o.push_back(v); - } - - return o; - } - - bool insert(Vertex<D, BiroliState>& v, const BiroliState& s, bool force = false) override { - if (force || (v.isEmpty() && canReplaceWith(v, s))) { - v.state.type = s.type; - for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) { - e.neighbor.state.occupiedNeighbors++; - } - this->N++; - return true; - } else { - return false; - } - } - - bool remove(Vertex<D, BiroliState>& v) override { - if (v.isEmpty()) { - return false; - } else { - v.state.remove(); - for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) { - e.neighbor.state.occupiedNeighbors--; - } - this->N--; - return true; - } - } - - bool randomMove(Rng& r) override { - Vertex<D, BiroliState>& v = r.pick(this->vertices); - - return exchange(v, r.pick(v.adjacentEdges).neighbor); - } - - void swap(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) override { - if (v1.state.type != v2.state.type) { - if (!v1.isEmpty() && !v2.isEmpty()) { - std::swap(v1.state.type, v2.state.type); - } else if (v1.isEmpty()) { - unsigned t = v2.state.type; - remove(v2); - insert(v1, BiroliState(t, 0), true); - } else { - unsigned t = v1.state.type; - remove(v1); - insert(v2, BiroliState(t, 0), true); - } - } - } - - bool exchange(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) override { - if (canReplaceWith(v1, v2.state) && canReplaceWith(v2, v1.state)) { - swap(v1, v2); - return true; - } else { - return false; - } - } - - int overlap(const System<D, BiroliState>& s) const override { - int o = 0; - - for (unsigned i = 0; i < this->vertices.size(); i++) { - BiroliState s2 = this->orientation.apply( - s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state); - if (s2.type > 0 && s2.type == this->vertices[i].state.type) { - o++; - } else { - o--; - } - } - - return o; - } -}; |