Tried to improve efficiency by switching to list iterators as keys.

1 files 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;
 }