diff options
| -rw-r--r-- | glass.cpp | 181 |
1 files changed, 91 insertions, 90 deletions
@@ -1,7 +1,6 @@ #include <iostream> #include <vector> #include <list> -#include <set> #include <eigen3/Eigen/Dense> #include "pcg-cpp/include/pcg_random.hpp" @@ -33,26 +32,39 @@ template<unsigned D> class State { σ[a] = b; } + State() : σ(Vector<D>::Zero()) {} + + bool isEmpty() const { + return σ.squaredNorm() == 0; + } + + void remove() { + σ = Vector<D>::Zero(); + } + State<D> apply(const Matrix<D>& m) const { return {m * σ}; } }; -template <unsigned D> class Vertex; +template<unsigned D> State<D> randomState(Rng& r) { + return State<D>(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})); +} + +template<unsigned D> class HalfEdge; -template<unsigned D> class Particle { +template<unsigned D> class Vertex { public: - Vertex<D>& position; + Vector<D> position; State<D> state; + std::vector<HalfEdge<D>> adjacentEdges; bool marked; - Particle(Vertex<D>& p, const State<D>& s) : position(p), state(s) {} + bool isEmpty() const { + return state.isEmpty(); + } }; -template<unsigned D> State<D> randomState(Rng& r) { - return State<D>(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})); -} - template<unsigned D> class HalfEdge { public: Vertex<D>& neighbor; @@ -61,22 +73,11 @@ template<unsigned D> class HalfEdge { HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {} }; -template<unsigned D> class Vertex { - public: - Vector<D> position; - std::vector<HalfEdge<D>> adjacentEdges; - Particle<D>* occupant; - - bool empty() const { - return occupant == NULL; - } -}; - template<unsigned D> class System { public: const unsigned L; + unsigned N; std::vector<Vertex<D>> vertices; - std::set<Particle<D>*> particles; unsigned vectorToIndex(const Vector<D>& x) const { unsigned i = 0; @@ -94,15 +95,10 @@ template<unsigned D> 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<unsigned D> class System { } } - std::list<Particle<D>*> overlaps(const Particle<D>& p, bool excludeSelf = false) const { - std::list<Particle<D>*> o; + 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 (!p.position.empty() && !excludeSelf) { - o.push_back(p.position.occupant); + if (!v.isEmpty() && !excludeSelf) { + o.push_back(v); } - for (const HalfEdge<D>& 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<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); } } } @@ -134,61 +134,65 @@ template<unsigned D> class System { return o; } - bool insert(const Particle<D>& p) { - if (overlaps(p).empty()) { - Particle<D>* pN = new Particle<D>(p); - particles.insert(pN); - p.position.occupant = pN; + bool tryInsertion(Vertex<D>& v, const State<D>& s) { + if (overlaps(v, s).empty()) { + v.state = s; + N++; return true; } else { return false; } } - void remove(Particle<D>* p) { - p->position.occupant = NULL; - particles.erase(p); - delete p; + bool tryDeletion(Vertex<D>& v) { + if (v.isEmpty()) { + return false; + } else { + v.state.remove(); + N--; + return true; + } } - /* - bool randomMove(Rng& r) { - Particle<D>& p = *(r.pick(particles)); - Particle<D> pN = p; - - if (0.2 < r.uniform(0.0, 1.0)) { - pN.position = r.pick(p.position.adjacentEdges).neighbor; + bool tryRandomMove(Rng& r) { + Vertex<D>& v = r.pick(vertices); + State<D> oldState = v.state; + if (!tryDeletion(v)) { + return false; } - pN.state = randomState<D>(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<D>(r))) { + return true; + } } else { - return false; + if (tryInsertion(r.pick(v.adjacentEdges).neighbor, randomState<D>(r))) { + return true; + } } + v.state = oldState; + N++; + return false; } - */ - bool swap(Particle<D>& p1, Particle<D>& p2) { - std::swap(p1.state, p2.state); - if (overlaps(p1, true).empty() && overlaps(p2, true).empty()) { + 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 { - 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<D>& v1 = r.pick(vertices); + Vertex<D>& v2 = r.pick(vertices); + + return trySwap(v1, v2); } void setGroundState() { - for (Particle<D>* p : particles) { - remove(p); - } + N = 0; for (Vertex<D>& v : vertices) { unsigned a = 0; @@ -197,21 +201,21 @@ template<unsigned D> class System { } a %= 2 * D + 1; - if (a > 0) { - if (a <= D) { - State<D> s(a - 1, -1); - insert(Particle<D>(v, s)); - } else if (D < a) { - State<D> s(2 * D - a, 1); - insert(Particle<D>(v, s)); - } + v.state.σ = 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 (Particle<D>* p : particles) { - if (!overlaps(*p, true).empty()) { + for (Vertex<D>& v : vertices) { + if (overlaps(v, v.state, true).size() > 0) { return false; } } @@ -220,7 +224,7 @@ template<unsigned D> class System { } double density() const { - return N() / pow(L, D); + return N / pow(L, D); } unsigned maxOccupation() const { @@ -231,40 +235,40 @@ template<unsigned D> 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<D>& v = r.pick(vertices); - if (v.empty()) { - insert(Particle<D>(v, randomState<D>(r))); + if (v.isEmpty()) { + tryInsertion(v, randomState<D>(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); |