From 15744ab1864b9ee7b59fbb050afd785bec54f8a8 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 23 Mar 2021 14:59:51 +0100 Subject: Working Biroli-Mezard clusters. --- glass.hpp | 241 ++++++++++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 205 insertions(+), 36 deletions(-) (limited to 'glass.hpp') diff --git a/glass.hpp b/glass.hpp index e67bb51..b35964b 100644 --- a/glass.hpp +++ b/glass.hpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -130,7 +131,7 @@ public: unsigned occupiedNeighbors; BiroliState() { - type = 0; + type = std::numeric_limits::max(); occupiedNeighbors = 0; } @@ -139,9 +140,13 @@ public: occupiedNeighbors = o; } - bool isEmpty() const { return type == 0; } + BiroliState(Rng& r) { + type = r.pick({1,2,2,2,2,2,3,3,3,3}); + } + + bool isEmpty() const { return type == std::numeric_limits::max(); } - void remove() { type = 0; }; + void remove() { type = std::numeric_limits::max(); }; }; template class Transformation { @@ -246,10 +251,9 @@ public: } } - virtual std::list>> overlaps(Vertex&, - const State&, bool = false) { return {}; }; + virtual std::list>> overlaps(Vertex&) { return {}; }; - virtual bool insert(Vertex& v, const State& s) { return false; }; + virtual bool insert(Vertex& v, const State& s, bool force = false) { return false; }; virtual bool remove(Vertex& v) { return false;}; @@ -268,7 +272,7 @@ public: bool compatible() { for (Vertex& v : vertices) { - if (overlaps(v, v.state, true).size() > 0) { + if (overlaps(v).size() > 0) { return false; } } @@ -319,43 +323,78 @@ public: } } - unsigned flipCluster(const Transformation& R, Vertex& v0, Rng& r, bool dry = false) { - std::queue>> q; - q.push(v0); - + unsigned flipCluster(const Transformation& R, Vertex& v0, bool dry = false) { unsigned n = 0; - while (!q.empty()) { - Vertex& v = q.front(); - q.pop(); + Vertex& v0New = vertices[vectorToIndex(R.apply(v0.position))]; - if (!v.marked) { - Vertex& vNew = vertices[vectorToIndex(R.apply(v.position))]; + if (&v0New != &v0) { + std::queue>> q; + + v0.marked = true; + v0New.marked = true; + + swap(v0, v0New); + + for (Vertex& vn : overlaps(v0)) { + if (!vn.marked) { + q.push(vn); + } + } + for (Vertex& vn : overlaps(v0New)) { + if (!vn.marked) { + q.push(vn); + } + } - v.marked = true; - vNew.marked = true; + while (!q.empty()) { + Vertex& v = q.front(); + q.pop(); - State s = R.apply(v.state); - State sNew = R.apply(vNew.state); + if (!v.marked && !overlaps(v).empty()) { + v.marked = true; + Vertex& vNew = vertices[vectorToIndex(R.apply(v.position))]; - std::list>> overlaps1 = overlaps(vNew, s, true); - std::list>> overlaps2 = overlaps(v, sNew, true); - overlaps1.splice(overlaps1.begin(), overlaps2); + if (&vNew != &v) { + vNew.marked = true; + + swap(v, vNew); + + for (Vertex& vn : overlaps(v)) { + if (!vn.marked) { + q.push(vn); + } + } + for (Vertex& vn : overlaps(vNew)) { + if (!vn.marked) { + q.push(vn); + } + } - for (Vertex& vn : overlaps1) { - if (!vn.marked) { - q.push(vn); + n += 2; + } else { + n += 1; } } - - if (!dry) { - swap(v, vNew); - v.state = sNew; - vNew.state = s; + if (q.empty()) { + for (Vertex& vv : vertices) { + if (!vv.marked && !overlaps(vv).empty()) { +// std::cerr << "Found bad state at end" << std::endl; + q.push(vv); + } + } } - - n += 1; } + + + } + + if (n > size() / 4) { + orientation = R.apply(orientation); + } + + for (Vertex& v : vertices) { + v.marked = false; } return n; @@ -365,7 +404,7 @@ public: for (Vertex& v : vertices) { if (!v.marked) { bool dry = 0.5 < r.uniform(0.0, 1.0); - unsigned n = flipCluster(R, v, r, dry); + unsigned n = flipCluster(R, v, dry); if (n > size() / 4 && !dry) { orientation = R.apply(orientation); } @@ -408,8 +447,8 @@ template class CiamarraSystem : public System> { return o; } - bool insert(Vertex>& v, const CiamarraState& s) override { - if (overlaps(v, s).empty()) { + bool insert(Vertex>& v, const CiamarraState& s, bool force = false) override { + if (force || overlaps(v, s).empty()) { v.state = s; this->N++; return true; @@ -506,3 +545,133 @@ template class CiamarraSystem : public System> { return o; } }; + +template class BiroliSystem : public System { +public: + using System::System; + + bool canReplaceWith(const Vertex& 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& e : v.adjacentEdges) { + if (e.neighbor.state.type < e.neighbor.state.occupiedNeighbors + Δn) { + return false; + } + } + } + return true; + } + + std::list>> + overlaps(Vertex& v) override { + std::list>> 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& 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& v, const BiroliState& s, bool force = false) override { + if (force || (v.isEmpty() && canReplaceWith(v, s))) { + v.state.type = s.type; + for (HalfEdge& e : v.adjacentEdges) { + e.neighbor.state.occupiedNeighbors++; + } + this->N++; + return true; + } else { + return false; + } + } + + bool remove(Vertex& v) override { + if (v.isEmpty()) { + return false; + } else { + v.state.remove(); + for (HalfEdge& e : v.adjacentEdges) { + e.neighbor.state.occupiedNeighbors--; + } + this->N--; + return true; + } + } + + bool randomMove(Rng& r) override { + Vertex& v = r.pick(this->vertices); + + return exchange(v, r.pick(v.adjacentEdges).neighbor); + } + + void swap(Vertex& v1, Vertex& 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& v1, Vertex& v2) override { + if (canReplaceWith(v1, v2.state) && canReplaceWith(v2, v1.state)) { + swap(v1, v2); + return true; + } else { + return false; + } + } + + int overlap(const System& 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; + } +}; -- cgit v1.2.3-54-g00ecf