diff options
| -rw-r--r-- | biroli-mezard.cpp | 214 | ||||
| -rw-r--r-- | glass.hpp | 4 |
2 files changed, 93 insertions, 125 deletions
diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp index 5b6213f..52a4360 100644 --- a/biroli-mezard.cpp +++ b/biroli-mezard.cpp @@ -24,41 +24,9 @@ public: void remove() { maximumNeighbors = Empty; } -} - -template <int D> class Vertex { -public: - Vector<D> initialPosition; - Vector<D> position; - std::vector<std::reference_wrapper<Vertex<D>>> neighbors; - unsigned occupiedNeighbors; - unsigned maximumNeighbors; - bool marked; - bool visited; - - Vertex() { - marked = false; - } - - bool empty() const { return maximumNeighbors == Empty; } - bool frustrated() const { return occupiedNeighbors > maximumNeighbors; } - bool valid() const { - unsigned o = 0; - - for (const Vertex<D>& v : neighbors) { - if (!v.empty()) { - o++; - } - } - - return o == occupiedNeighbors; - } - void setInitial() { - initialPosition = position; - } }; -template <int D> class BiroliSystem : public HardSystem<D, BiroliState> { +template <int D> class BiroliSystem : public System<D, BiroliState> { public: std::vector<unsigned> N; std::unordered_set<Vertex<D, BiroliState>*> occupants; @@ -69,11 +37,23 @@ public: return BiroliSystem::size() - N[0]; } - double density() const { return (double)occupancy() / size(); } + double density() const { return (double)occupancy() / BiroliSystem::size(); } + + bool valid(const Vertex<D, BiroliState>& v) const { + unsigned o = 0; + + for (const Vertex<D, BiroliState>& vn : v.neighbors) { + if (!vn.empty()) { + o++; + } + } + + return o == v.occupiedNeighbors; + } bool compatible() const { - for (const Vertex<D>& v : vertices) { - if (v.frustrated()) { + for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) { + if (v.state.frustrated()) { return false; } } @@ -82,8 +62,8 @@ public: } bool valid() const { - for (const Vertex<D>& v : vertices) { - if (!v.valid()) { + for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) { + if (!valid(v)) { return false; } } @@ -91,25 +71,12 @@ public: return true; } - BiroliSystem(unsigned L, unsigned T) : L(L), N(T + 1, 0), vertices(iPow(L, D)), orientation(L) { - N[0] = size(); - - for (unsigned i = 0; i < size(); i++) { - vertices[i].position = indexToVector(i); - vertices[i].neighbors.reserve(2 * D); - } - - for (unsigned d = 0; d < D; d++) { - for (signed i = 0; i < size(); i++) { - unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1)); - vertices[i].neighbors.push_back(vertices[j]); - vertices[j].neighbors.push_back(vertices[i]); - } - } + BiroliSystem(unsigned L, unsigned T) : System<D, BiroliState>(L), N(T + 1, 0) { + N[0] = BiroliSystem::size(); } - std::list<std::reference_wrapper<Vertex<D>>> overlaps(Vertex<D>& v) { - std::list<std::reference_wrapper<Vertex<D>>> o; + std::list<std::reference_wrapper<Vertex<D, BiroliState>>> overlaps(Vertex<D, BiroliState>& v) { + std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o; if (v.empty()) { // an empty site cannot be frustrating anyone return o; @@ -118,7 +85,7 @@ public: bool selfFrustrated = v.frustrated(); bool anyNeighborFrustrated = false; - for (Vertex<D>& vn : v.neighbors) { + for (Vertex<D, BiroliState>& vn : v.neighbors) { if (!vn.empty()) { bool thisNeighborFrustrated = vn.frustrated(); @@ -139,28 +106,29 @@ public: return o; } - std::list<std::reference_wrapper<Vertex<D>>> overlaps(const std::list<std::reference_wrapper<Vertex<D>>>& vs) { - std::list<std::reference_wrapper<Vertex<D>>> o; + std::list<std::reference_wrapper<Vertex<D, BiroliState>>> overlaps(const std::list<std::reference_wrapper<Vertex<D, BiroliState>>>& vs) { + std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o; - for (Vertex<D>& v : vs) { + for (Vertex<D, BiroliState>& v : vs) { o.splice(o.begin(), overlaps(v)); } return o; } - bool canReplaceWith(const Vertex<D>& v, unsigned t) const { + bool canReplaceWith(const Vertex<D, BiroliState>& v, unsigned t) const { if (t == Empty) { return true; } - if (t < v.occupiedNeighbors) { + if (t < v.state.occupiedNeighbors) { return false; } if (v.empty()) { - for (const Vertex<D>& vn : v.neighbors) { - if (vn.maximumNeighbors < vn.occupiedNeighbors + 1) { + for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) { + const Vertex<D, BiroliState>& vn = e.neighbor; + if (vn.state.maximumNeighbors < vn.state.occupiedNeighbors + 1) { return false; } } @@ -170,16 +138,16 @@ public: } void setInitial() { - for (Vertex<D>& v : vertices) { + for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) { v.setInitial(); } } - bool insert(Vertex<D>& v, unsigned t, bool force = false) { + bool insert(Vertex<D, BiroliState>& v, unsigned t, bool force = false) { if (force || (v.empty() && canReplaceWith(v, t))) { - v.maximumNeighbors = t; - for (Vertex<D>& vn : v.neighbors) { - vn.occupiedNeighbors++; + v.state.maximumNeighbors = t; + for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) { + e.neighbor.state.occupiedNeighbors++; } N[t]++; N[0]--; @@ -192,39 +160,39 @@ public: } } - bool remove(Vertex<D>& v) { + bool remove(Vertex<D, BiroliState>& v) { if (v.empty()) { return false; } else { - if (v.maximumNeighbors > 0) { + if (v.state.maximumNeighbors > 0) { occupants.erase(&v); } - N[v.maximumNeighbors]--; + N[v.state.maximumNeighbors]--; N[0]++; - v.maximumNeighbors = Empty; - for (Vertex<D>& vn : v.neighbors) { - vn.occupiedNeighbors--; + v.state.maximumNeighbors = Empty; + for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) { + e.neighbor.state.occupiedNeighbors--; } return true; } } bool randomLocalExchange(Rng& r) { - Vertex<D>& v = r.pick(vertices); + Vertex<D, BiroliState>& v = r.pick(BiroliSystem::vertices); return exchange(v, r.pick(v.neighbors)); } - void swap(Vertex<D>& v1, Vertex<D>& v2) { - if (v1.maximumNeighbors != v2.maximumNeighbors) { + void swap(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) { + if (v1.state.maximumNeighbors != v2.state.maximumNeighbors) { if (!v1.empty() && !v2.empty()) { - std::swap(v1.maximumNeighbors, v2.maximumNeighbors); + std::swap(v1.state.maximumNeighbors, v2.state.maximumNeighbors); } else if (v1.empty()) { - unsigned t = v2.maximumNeighbors; + unsigned t = v2.state.maximumNeighbors; remove(v2); insert(v1, t, true); } else { - unsigned t = v1.maximumNeighbors; + unsigned t = v1.state.maximumNeighbors; remove(v1); insert(v2, t, true); } @@ -232,8 +200,8 @@ public: } } - bool exchange(Vertex<D>& v1, Vertex<D>& v2) { - if (canReplaceWith(v1, v2.maximumNeighbors) && canReplaceWith(v2, v1.maximumNeighbors)) { + bool exchange(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) { + if (canReplaceWith(v1, v2.state.maximumNeighbors) && canReplaceWith(v2, v1.state.maximumNeighbors)) { swap(v1, v2); return true; } else { @@ -246,10 +214,10 @@ public: using namespace std::complex_literals; - for (const Vertex<D>& v : vertices) { + for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) { if (v.maximumNeighbors == t) { for (unsigned d = 0; d < D; d++) { - F += exp(1i * 2.0 * M_PI * (double)k * (double)(v.position(d) - v.initialPosition(d)) / (double)L); + F += exp(1i * 2.0 * M_PI * (double)k * (double)(v.position(d) - v.initialPosition(d)) / (double)BiroliSystem::L); } } } @@ -259,11 +227,11 @@ public: return F; } - double overlap(const System<D>& s, unsigned t = 2) const { + double overlap(const BiroliSystem<D>& s, unsigned t = 2) const { int o = 0; - for (const Vertex<D>& v1 : vertices) { - const Vertex<D>& v2 = s.vertices[vectorToIndex(orientation.inverse().apply(v1.position))]; + for (const Vertex<D, BiroliState>& v1 : BiroliSystem::vertices) { + const Vertex<D, BiroliState>& v2 = s.BiroliSystem::vertices[vectorToIndex(BiroliSystem::orientation.inverse().apply(v1.position))]; if (v1.maximumNeighbors == v2.maximumNeighbors) { o++; } @@ -275,12 +243,12 @@ public: nn += iPow(n, 2); } - return ((double)o / (double)size() - (double)nn / pow((double)size(), 2)) / (1 - (double)nn / pow((double)size(), 2)); + return ((double)o / (double)BiroliSystem::size() - (double)nn / pow((double)BiroliSystem::size(), 2)) / (1 - (double)nn / pow((double)BiroliSystem::size(), 2)); } bool randomExchange(Rng& r) { - Vertex<D>& v1 = r.pick(vertices); - Vertex<D>& v2 = r.pick(vertices); + Vertex<D, BiroliState>& v1 = r.pick(BiroliSystem::vertices); + Vertex<D, BiroliState>& v2 = r.pick(BiroliSystem::vertices); return exchange(v1, v2); } @@ -290,24 +258,24 @@ public: 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 pIns = BiroliSystem::size() * z / (occupancy() + 1); if (pIns > r.uniform(0.0, 1.0)) { - Vertex<D>& v = r.pick(vertices); + Vertex<D, BiroliState>& v = r.pick(BiroliSystem::vertices); insert(v, t); } } else { double pDel = density() / z; if (pDel > r.uniform(0.0, 1.0)) { - Vertex<D>* v = r.pick(occupants); + Vertex<D, BiroliState>* v = r.pick(occupants); remove(*v); } } } void sweepGrandCanonical(double z, Rng& r) { - for (unsigned i = 0; i < size(); i++) { + for (unsigned i = 0; i < BiroliSystem::size(); i++) { stepGrandCanonical(z, r); randomExchange(r); @@ -315,36 +283,36 @@ public: } void sweepLocal(Rng& r) { - for (unsigned i = 0; i < size() / 2; i++) { + for (unsigned i = 0; i < BiroliSystem::size() / 2; i++) { randomLocalExchange(r); } } void sweepSwap(Rng& r) { - for (unsigned i = 0; i < size() / 2; i++) { + for (unsigned i = 0; i < BiroliSystem::size() / 2; i++) { randomExchange(r); } } - Vertex<D>& apply(const Transformation<D>& R, const Vertex<D>& v) { - return vertices[vectorToIndex(R.apply(v.position))]; + Vertex<D, BiroliState>& apply(const Transformation<D>& R, const Vertex<D, BiroliState>& v) { + return BiroliSystem::vertices[vectorToIndex(R.apply(v.position))]; } - bool eventChain(const Transformation<D>& R, Vertex<D>& v0) { + bool eventChain(const Transformation<D>& R, Vertex<D, BiroliState>& v0) { unsigned n = 0; if (!v0.empty()) { - Vertex<D>& v0New = apply(R, v0);; + Vertex<D, BiroliState>& v0New = apply(R, v0);; unsigned t0 = v0.maximumNeighbors; - std::queue<std::pair<std::reference_wrapper<Vertex<D>>, unsigned>> q; + std::queue<std::pair<std::reference_wrapper<Vertex<D, BiroliState>>, unsigned>> q; q.push({v0New, t0}); remove(v0); while (!q.empty()) { auto [vr, t] = q.front(); - Vertex<D>& v = vr; + Vertex<D, BiroliState>& v = vr; q.pop(); if (!v.empty()) { @@ -355,7 +323,7 @@ public: insert(v, t, true); v.marked = true; - for (Vertex<D>& vn : overlaps({v})) { + for (Vertex<D, BiroliState>& vn : overlaps({v})) { if (!vn.marked) { q.push({apply(R, vn), vn.maximumNeighbors}); remove(vn); @@ -365,7 +333,7 @@ public: * frustrated, it may be due to one of its neighbors! */ if (&vn != &v) { - for (Vertex<D>& vnn : overlaps(vn)) { + for (Vertex<D, BiroliState>& vnn : overlaps(vn)) { if (!vnn.marked) { q.push({apply(R, vnn), vnn.maximumNeighbors}); remove(vnn); @@ -378,7 +346,7 @@ public: } } - for (Vertex<D>& v : vertices) { + for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) { v.marked = false; } @@ -386,44 +354,44 @@ public: } bool randomEventChain(Rng& r) { - Transformation<D> R(L); + Transformation<D> R(BiroliSystem::L); R.v[r.uniform((unsigned)0, (unsigned)D-1)] = r.pick({-1, 1}); - return eventChain(R, r.pick(vertices)); + return eventChain(R, r.pick(BiroliSystem::vertices)); } - unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0) { + unsigned flipCluster(const Transformation<D>& R, Vertex<D, BiroliState>& v0) { unsigned n = 0; - Vertex<D>& v0New = vertices[vectorToIndex(R.apply(v0.position))]; + Vertex<D, BiroliState>& v0New = BiroliSystem::vertices[vectorToIndex(R.apply(v0.position))]; if (&v0New != &v0) { - std::queue<std::reference_wrapper<Vertex<D>>> q; + std::queue<std::reference_wrapper<Vertex<D, BiroliState>>> q; v0.marked = true; v0New.marked = true; swap(v0, v0New); - for (Vertex<D>& vn : overlaps({v0, v0New})) { + for (Vertex<D, BiroliState>& vn : overlaps({v0, v0New})) { if (!vn.marked) { q.push(vn); } } while (!q.empty()) { - Vertex<D>& v = q.front(); + Vertex<D, BiroliState>& v = q.front(); q.pop(); if (!v.marked && !overlaps(v).empty()) { v.marked = true; - Vertex<D>& vNew = vertices[vectorToIndex(R.apply(v.position))]; + Vertex<D, BiroliState>& vNew = BiroliSystem::vertices[vectorToIndex(R.apply(v.position))]; if (&vNew != &v) { vNew.marked = true; swap(v, vNew); - for (Vertex<D>& vn : overlaps({v, vNew})) { + for (Vertex<D, BiroliState>& vn : overlaps({v, vNew})) { if (!vn.marked) { q.push(vn); } else { @@ -431,7 +399,7 @@ public: * frustrated, it may be due to one of its neighbors! */ if (&vn != &v && &vn != &vNew) { - for (Vertex<D>& vnn : overlaps(vn)) { + for (Vertex<D, BiroliState>& vnn : overlaps(vn)) { if (!vnn.marked) { q.push(vnn); } @@ -449,11 +417,11 @@ public: n = 1; } - if (n > size() / 2) { - orientation = R.apply(orientation); + if (n > BiroliSystem::size() / 2) { + BiroliSystem::orientation = R.apply(BiroliSystem::orientation); } - for (Vertex<D>& v : vertices) { + for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) { v.marked = false; } @@ -461,10 +429,10 @@ public: } }; -void print(const System<2>& s) { - for (const Vertex<2>& v : s.vertices) { +void print(const BiroliSystem<2>& s) { + for (const Vertex<2, BiroliState>& v : s.vertices) { if (!v.empty()) { - std::cerr << v.maximumNeighbors; + std::cerr << v.state.maximumNeighbors; } else { std::cerr << " "; } @@ -484,14 +452,14 @@ int main() { double μmax = 10; double dμ = 0.05; - System<D> s(L, 3); + BiroliSystem<D> s(L, 3); Rng r; for (double μ = μmin; μ <= μmax; μ += dμ) { double z = exp(μ); for (unsigned i = 0; i < 1e4; i++) { - for (unsigned j = 0; j < s.size(); j++) { + for (unsigned j = 0; j < s.BiroliSystem::size(); j++) { s.stepGrandCanonical(z, r); // s.randomEventChain(r); s.randomExchange(r); @@ -573,7 +573,7 @@ public: bool dry = 0.5 < r.uniform(0.0, 1.0); unsigned n = flipCluster(R, v, r, dry); if (n > pow(HardSystem::L, D) / 4 && !dry) { - orientation = R.apply(orientation); + HardSystem::orientation = R.apply(HardSystem::orientation); } } } @@ -587,7 +587,7 @@ public: int o = 0; for (unsigned i = 0; i < HardSystem::vertices.size(); i++) { - S s2 = orientation.apply(s.vertices[HardSystem::vectorToIndex(orientation.inverse().apply(HardSystem::indexToVector(i)))].state); + S s2 = HardSystem::orientation.apply(s.vertices[HardSystem::vectorToIndex(HardSystem::orientation.inverse().apply(HardSystem::indexToVector(i)))].state); o += HardSystem::vertices[i].state.dot(s2); } |