From 93bdaae2bfd869600d9ea2fbb86dd555c6b7608b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 20 Mar 2021 19:07:07 +0100 Subject: Reverted attempt to make things a little more abstract, because it was very slow for some reason. --- glass.cpp | 181 +++++++++++++++++++++++++++++++------------------------------- 1 file changed, 91 insertions(+), 90 deletions(-) diff --git a/glass.cpp b/glass.cpp index e0c5e25..f5e9b31 100644 --- a/glass.cpp +++ b/glass.cpp @@ -1,7 +1,6 @@ #include #include #include -#include #include #include "pcg-cpp/include/pcg_random.hpp" @@ -33,26 +32,39 @@ template class State { σ[a] = b; } + State() : σ(Vector::Zero()) {} + + bool isEmpty() const { + return σ.squaredNorm() == 0; + } + + void remove() { + σ = Vector::Zero(); + } + State apply(const Matrix& m) const { return {m * σ}; } }; -template class Vertex; +template State randomState(Rng& r) { + return State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})); +} + +template class HalfEdge; -template class Particle { +template class Vertex { public: - Vertex& position; + Vector position; State state; + std::vector> adjacentEdges; bool marked; - Particle(Vertex& p, const State& s) : position(p), state(s) {} + bool isEmpty() const { + return state.isEmpty(); + } }; -template State randomState(Rng& r) { - return State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})); -} - template class HalfEdge { public: Vertex& neighbor; @@ -61,22 +73,11 @@ template class HalfEdge { HalfEdge(Vertex& n, const Vector& d) : neighbor(n), Δx(d) {} }; -template class Vertex { - public: - Vector position; - std::vector> adjacentEdges; - Particle* occupant; - - bool empty() const { - return occupant == NULL; - } -}; - template class System { public: const unsigned L; + unsigned N; std::vector> vertices; - std::set*> particles; unsigned vectorToIndex(const Vector& x) const { unsigned i = 0; @@ -94,15 +95,10 @@ template class System { return x; } - unsigned N() const { - return particles.size(); - } - - System(unsigned L) : L(L), vertices(iPow(L, D)) { + System(unsigned L) : L(L), N(0), vertices(iPow(L, D)) { for (unsigned i = 0; i < iPow(L, D); i++) { vertices[i].position = indexToVector(i); vertices[i].adjacentEdges.reserve(2 * D); - vertices[i].occupant = NULL; } for (unsigned d = 0; d < D; d++) { @@ -116,17 +112,21 @@ template class System { } } - std::list*> overlaps(const Particle& p, bool excludeSelf = false) const { - std::list*> o; + std::list>> overlaps(Vertex& v, const State& s, bool excludeSelf = false) { + std::list>> o; + + if (s.isEmpty()) { + return o; + } - if (!p.position.empty() && !excludeSelf) { - o.push_back(p.position.occupant); + if (!v.isEmpty() && !excludeSelf) { + o.push_back(v); } - for (const HalfEdge& e : p.position.adjacentEdges) { - if (!e.neighbor.empty()) { - if (p.state.σ.dot(e.Δx) == 1 || e.neighbor.occupant->state.σ.dot(e.Δx) == -1) { - o.push_back(e.neighbor.occupant); + for (const HalfEdge& 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); } } } @@ -134,61 +134,65 @@ template class System { return o; } - bool insert(const Particle& p) { - if (overlaps(p).empty()) { - Particle* pN = new Particle(p); - particles.insert(pN); - p.position.occupant = pN; + bool tryInsertion(Vertex& v, const State& s) { + if (overlaps(v, s).empty()) { + v.state = s; + N++; return true; } else { return false; } } - void remove(Particle* p) { - p->position.occupant = NULL; - particles.erase(p); - delete p; + bool tryDeletion(Vertex& v) { + if (v.isEmpty()) { + return false; + } else { + v.state.remove(); + N--; + return true; + } } - /* - bool randomMove(Rng& r) { - Particle& p = *(r.pick(particles)); - Particle pN = p; - - if (0.2 < r.uniform(0.0, 1.0)) { - pN.position = r.pick(p.position.adjacentEdges).neighbor; + bool tryRandomMove(Rng& r) { + Vertex& v = r.pick(vertices); + State oldState = v.state; + if (!tryDeletion(v)) { + return false; } - pN.state = randomState(r); - if (overlaps(pN, true).empty()) { - remove(&p); - insert(pN); - return true; + if (0.2 > r.uniform(0.0, 1.0)) { + if (tryInsertion(v, randomState(r))) { + return true; + } } else { - return false; + if (tryInsertion(r.pick(v.adjacentEdges).neighbor, randomState(r))) { + return true; + } } + v.state = oldState; + N++; + return false; } - */ - bool swap(Particle& p1, Particle& p2) { - std::swap(p1.state, p2.state); - if (overlaps(p1, true).empty() && overlaps(p2, true).empty()) { + bool trySwap(Vertex& v1, Vertex& 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 { - std::swap(p1.state, p2.state); return false; } } - bool randomSwap(Rng& r) { - return swap(*r.pick(particles), *r.pick(particles)); + bool tryRandomSwap(Rng& r) { + Vertex& v1 = r.pick(vertices); + Vertex& v2 = r.pick(vertices); + + return trySwap(v1, v2); } void setGroundState() { - for (Particle* p : particles) { - remove(p); - } + N = 0; for (Vertex& v : vertices) { unsigned a = 0; @@ -197,21 +201,21 @@ template class System { } a %= 2 * D + 1; - if (a > 0) { - if (a <= D) { - State s(a - 1, -1); - insert(Particle(v, s)); - } else if (D < a) { - State s(2 * D - a, 1); - insert(Particle(v, s)); - } + v.state.σ = Vector::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 (Particle* p : particles) { - if (!overlaps(*p, true).empty()) { + for (Vertex& v : vertices) { + if (overlaps(v, v.state, true).size() > 0) { return false; } } @@ -220,7 +224,7 @@ template class System { } double density() const { - return N() / pow(L, D); + return N / pow(L, D); } unsigned maxOccupation() const { @@ -231,40 +235,40 @@ template class System { 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); + double pIns = maxOccupation() * z / (N + 1); if (pIns > r.uniform(0.0, 1.0)) { while (true) { Vertex& v = r.pick(vertices); - if (v.empty()) { - insert(Particle(v, randomState(r))); + if (v.isEmpty()) { + tryInsertion(v, randomState(r)); break; } } } } else { - double pDel = N() / (z * maxOccupation()); + double pDel = N / (z * maxOccupation()); if (pDel > r.uniform(0.0, 1.0)) { - remove(r.pick(particles)); + tryDeletion(r.pick(vertices)); } } -// randomMove(r); - randomSwap(r); + tryRandomMove(r); + tryRandomSwap(r); } } void sweepLocal(double z, Rng& r) { for (unsigned i = 0; i < iPow(L, D); i++) { -// randomMove(r); + tryRandomMove(r); } } void sweepSwap(double z, Rng& r) { for (unsigned i = 0; i < iPow(L, D); i++) { - randomSwap(r); + tryRandomSwap(r); } } @@ -311,9 +315,6 @@ int main() { for (unsigned N = Nmin; N <= Nmax; N *= 10) { s.setGroundState(); - if (!s.compatible()) { - return 1; - } for (double T = Tmin; T < Tmax; T *= 1 + δT) { double z = exp(1 / T); -- cgit v1.2.3-70-g09d2