From 2187d5abb39f36fe342acffa94b84a392ccee0ae Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
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 <iostream>
#include <vector>
#include <list>
+#include <set>
#include <eigen3/Eigen/Dense>
#include "pcg-cpp/include/pcg_random.hpp"
@@ -11,6 +12,7 @@ using Rng = randutils::random_generator<pcg32>;
template<int D> using Vector = Eigen::Matrix<int, D, 1>;
template<int D> using Matrix = Eigen::Matrix<int, D, D>;
+// Integer power function
int iPow(int x, unsigned p) {
if (p == 0) return 1;
if (p == 1) return x;
@@ -32,37 +34,22 @@ template<unsigned D> class State {
σ[a] = b;
}
- State() : σ(Vector<D>::Zero()) {}
-
- bool isEmpty() const {
- return σ.squaredNorm() == 0;
- }
-
- void remove() {
- σ = Vector<D>::Zero();
- }
+ State(Rng& r) : State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {}
State<D> apply(const Matrix<D>& m) const {
return {m * σ};
}
};
-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 Vertex;
-template<unsigned D> class Vertex {
+template<unsigned D> class Particle {
public:
- Vector<D> position;
+ Vertex<D>& position;
State<D> state;
- std::vector<HalfEdge<D>> adjacentEdges;
bool marked;
- bool isEmpty() const {
- return state.isEmpty();
- }
+ Particle(Vertex<D>& p, const State<D>& s) : position(p), state(s) {}
};
template<unsigned D> class HalfEdge {
@@ -73,11 +60,18 @@ 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;
+ typename std::list<Particle<D>>::iterator occupant;
+};
+
template<unsigned D> class System {
public:
const unsigned L;
- unsigned N;
std::vector<Vertex<D>> vertices;
+ std::list<Particle<D>> particles;
unsigned vectorToIndex(const Vector<D>& x) const {
unsigned i = 0;
@@ -95,10 +89,19 @@ template<unsigned D> class System {
return x;
}
- System(unsigned L) : L(L), N(0), vertices(iPow(L, D)) {
+ bool isEmpty(const Vertex<D>& 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<unsigned D> class System {
}
}
- 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;
+ std::list<typename std::list<Particle<D>>::iterator> overlaps(const Particle<D>& p, bool excludeSelf = false) const {
+ std::list<typename std::list<Particle<D>>::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<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);
+ for (const HalfEdge<D>& 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<unsigned D> class System {
return o;
}
- bool tryInsertion(Vertex<D>& v, const State<D>& s) {
- if (overlaps(v, s).empty()) {
- v.state = s;
- N++;
+ bool insert(const Particle<D>& p) {
+ if (overlaps(p).empty()) {
+ particles.push_back(p);
+ p.position.occupant = std::prev(particles.end());
return true;
} else {
return false;
}
}
- bool tryDeletion(Vertex<D>& v) {
- if (v.isEmpty()) {
- return false;
- } else {
- v.state.remove();
- N--;
- return true;
- }
+ void remove(typename std::list<Particle<D>>::iterator ip) {
+ ip->position.occupant = particles.end();
+ particles.erase(ip);
}
- bool tryRandomMove(Rng& r) {
- Vertex<D>& v = r.pick(vertices);
- State<D> oldState = v.state;
- if (!tryDeletion(v)) {
- return false;
+ /*
+ 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;
}
+ pN.state = randomState<D>(r);
- if (0.2 > r.uniform(0.0, 1.0)) {
- if (tryInsertion(v, randomState<D>(r))) {
- return true;
- }
+ if (overlaps(pN, true).empty()) {
+ remove(&p);
+ insert(pN);
+ return true;
} else {
- if (tryInsertion(r.pick(v.adjacentEdges).neighbor, randomState<D>(r))) {
- return true;
- }
+ return false;
}
- v.state = oldState;
- N++;
- return false;
}
+ */
- 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);
+ bool swap(Particle<D>& p1, Particle<D>& 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<D>& v1 = r.pick(vertices);
- Vertex<D>& 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<Particle<D>>::iterator ip = particles.begin(); ip != particles.end(); ip++) {
+ remove(ip);
+ }
for (Vertex<D>& v : vertices) {
unsigned a = 0;
@@ -201,21 +194,21 @@ template<unsigned D> class System {
}
a %= 2 * D + 1;
- 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++;
+ 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));
+ }
}
}
}
bool compatible() {
- for (Vertex<D>& v : vertices) {
- if (overlaps(v, v.state, true).size() > 0) {
+ for (Particle<D> p : particles) {
+ if (!overlaps(p, true).empty()) {
return false;
}
}
@@ -224,7 +217,7 @@ template<unsigned D> class System {
}
double density() const {
- return N / pow(L, D);
+ return N() / pow(L, D);
}
unsigned maxOccupation() const {
@@ -235,40 +228,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.isEmpty()) {
- tryInsertion(v, randomState<D>(r));
+ if (isEmpty(v)) {
+ insert(Particle<D>(v, State<D>(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-70-g09d2