From 2187d5abb39f36fe342acffa94b84a392ccee0ae Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 20 Mar 2021 19:52:48 +0100 Subject: Tried to improve efficiency by switching to list iterators as keys. --- glass.cpp | 184 ++++++++++++++++++++++++++++++++------------------------------ 1 file changed, 94 insertions(+), 90 deletions(-) diff --git a/glass.cpp b/glass.cpp index f5e9b31..b6241a2 100644 --- a/glass.cpp +++ b/glass.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include "pcg-cpp/include/pcg_random.hpp" @@ -11,6 +12,7 @@ using Rng = randutils::random_generator; template using Vector = Eigen::Matrix; template using Matrix = Eigen::Matrix; +// Integer power function int iPow(int x, unsigned p) { if (p == 0) return 1; if (p == 1) return x; @@ -32,37 +34,22 @@ template class State { σ[a] = b; } - State() : σ(Vector::Zero()) {} - - bool isEmpty() const { - return σ.squaredNorm() == 0; - } - - void remove() { - σ = Vector::Zero(); - } + State(Rng& r) : State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {} State apply(const Matrix& m) const { return {m * σ}; } }; -template State randomState(Rng& r) { - return State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})); -} - -template class HalfEdge; +template class Vertex; -template class Vertex { +template class Particle { public: - Vector position; + Vertex& position; State state; - std::vector> adjacentEdges; bool marked; - bool isEmpty() const { - return state.isEmpty(); - } + Particle(Vertex& p, const State& s) : position(p), state(s) {} }; template class HalfEdge { @@ -73,11 +60,18 @@ template class HalfEdge { HalfEdge(Vertex& n, const Vector& d) : neighbor(n), Δx(d) {} }; +template class Vertex { + public: + Vector position; + std::vector> adjacentEdges; + typename std::list>::iterator occupant; +}; + template class System { public: const unsigned L; - unsigned N; std::vector> vertices; + std::list> particles; unsigned vectorToIndex(const Vector& x) const { unsigned i = 0; @@ -95,10 +89,19 @@ template class System { return x; } - System(unsigned L) : L(L), N(0), vertices(iPow(L, D)) { + bool isEmpty(const Vertex& v) const { + return v.occupant == particles.end(); + } + + unsigned N() const { + return particles.size(); + } + + System(unsigned L) : L(L), 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 = particles.end(); } for (unsigned d = 0; d < D; d++) { @@ -112,21 +115,17 @@ template class System { } } - std::list>> overlaps(Vertex& v, const State& s, bool excludeSelf = false) { - std::list>> o; + std::list>::iterator> overlaps(const Particle& p, bool excludeSelf = false) const { + std::list>::iterator> o; - if (s.isEmpty()) { - return o; + if (!excludeSelf && !isEmpty(p.position)) { + o.push_back(p.position.occupant); } - if (!v.isEmpty() && !excludeSelf) { - o.push_back(v); - } - - 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); + for (const HalfEdge& e : p.position.adjacentEdges) { + if (!isEmpty(e.neighbor)) { + if (p.state.σ.dot(e.Δx) == 1 || e.neighbor.occupant->state.σ.dot(e.Δx) == -1) { + o.push_back(e.neighbor.occupant); } } } @@ -134,65 +133,59 @@ template class System { return o; } - bool tryInsertion(Vertex& v, const State& s) { - if (overlaps(v, s).empty()) { - v.state = s; - N++; + bool insert(const Particle& p) { + if (overlaps(p).empty()) { + particles.push_back(p); + p.position.occupant = std::prev(particles.end()); return true; } else { return false; } } - bool tryDeletion(Vertex& v) { - if (v.isEmpty()) { - return false; - } else { - v.state.remove(); - N--; - return true; - } + void remove(typename std::list>::iterator ip) { + ip->position.occupant = particles.end(); + particles.erase(ip); } - bool tryRandomMove(Rng& r) { - Vertex& v = r.pick(vertices); - State oldState = v.state; - if (!tryDeletion(v)) { - return false; + /* + 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; } + pN.state = randomState(r); - if (0.2 > r.uniform(0.0, 1.0)) { - if (tryInsertion(v, randomState(r))) { - return true; - } + if (overlaps(pN, true).empty()) { + remove(&p); + insert(pN); + return true; } else { - if (tryInsertion(r.pick(v.adjacentEdges).neighbor, randomState(r))) { - return true; - } + return false; } - v.state = oldState; - N++; - return false; } + */ - 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); + bool swap(Particle& p1, Particle& p2) { + std::swap(p1.state, p2.state); + if (overlaps(p1, true).empty() && overlaps(p2, true).empty()) { return true; } else { + std::swap(p1.state, p2.state); return false; } } - bool tryRandomSwap(Rng& r) { - Vertex& v1 = r.pick(vertices); - Vertex& v2 = r.pick(vertices); - - return trySwap(v1, v2); + bool randomSwap(Rng& r) { + return swap(r.pick(particles), r.pick(particles)); } void setGroundState() { - N = 0; + for (typename std::list>::iterator ip = particles.begin(); ip != particles.end(); ip++) { + remove(ip); + } for (Vertex& v : vertices) { unsigned a = 0; @@ -201,21 +194,21 @@ template class System { } a %= 2 * D + 1; - 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++; + 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)); + } } } } bool compatible() { - for (Vertex& v : vertices) { - if (overlaps(v, v.state, true).size() > 0) { + for (Particle p : particles) { + if (!overlaps(p, true).empty()) { return false; } } @@ -224,7 +217,7 @@ template class System { } double density() const { - return N / pow(L, D); + return N() / pow(L, D); } unsigned maxOccupation() const { @@ -235,40 +228,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.isEmpty()) { - tryInsertion(v, randomState(r)); + if (isEmpty(v)) { + insert(Particle(v, State(r))); break; } } } } else { - double pDel = N / (z * maxOccupation()); + double pDel = N() / (z * maxOccupation()); if (pDel > r.uniform(0.0, 1.0)) { - tryDeletion(r.pick(vertices)); + remove(r.choose(particles)); } } - tryRandomMove(r); - tryRandomSwap(r); +// randomMove(r); + randomSwap(r); } } void sweepLocal(double z, Rng& r) { for (unsigned i = 0; i < iPow(L, D); i++) { - tryRandomMove(r); +// randomMove(r); } } void sweepSwap(double z, Rng& r) { for (unsigned i = 0; i < iPow(L, D); i++) { - tryRandomSwap(r); + randomSwap(r); } } @@ -313,8 +306,12 @@ int main() { Rng r; + /* 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); @@ -335,6 +332,13 @@ int main() { std::cout << L << " " << T << " " << N << " " << δT << " " << 1 / s.density() << std::endl; } } + */ + + s.setGroundState(); + + for (unsigned n = 0; n < 10; n++) { + s.sweepGrandCanonical(exp(1/0.15), r); + } return 0; } -- cgit v1.2.3-54-g00ecf