1 files changed, 19 insertions, 16 deletions

diff --git a/glass.cpp b/glass.cpp
index b6241a2..7623765 100644
--- a/glass.cpp
+++ b/glass.cpp
@@ -1,8 +1,10 @@
 #include <iostream>
+#include <limits>
 #include <vector>
 #include <list>
 #include <set>
 #include <eigen3/Eigen/Dense>
+#include <unordered_map>
 
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
@@ -64,14 +66,14 @@ template<unsigned D> class Vertex {
   public:
     Vector<D> position;
     std::vector<HalfEdge<D>> adjacentEdges;
-    typename std::list<Particle<D>>::iterator occupant;
+    typename std::unordered_map<unsigned, Particle<D>>::iterator occupant;
 };
 
 template<unsigned D> class System {
   public:
     const unsigned L;
     std::vector<Vertex<D>> vertices;
-    std::list<Particle<D>> particles;
+    std::unordered_map<unsigned, Particle<D>> particles;
 
     unsigned vectorToIndex(const Vector<D>& x) const {
       unsigned i = 0;
@@ -115,8 +117,8 @@ template<unsigned D> class System {
       }
     }
 
-    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;
+    std::list<typename std::unordered_map<unsigned, Particle<D>>::iterator> overlaps(const Particle<D>& p, bool excludeSelf = false) const {
+      std::list<typename std::unordered_map<unsigned, Particle<D>>::iterator> o;
 
       if (!excludeSelf && !isEmpty(p.position)) {
         o.push_back(p.position.occupant);
@@ -124,7 +126,7 @@ template<unsigned D> class System {
 
       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) {
+          if (p.state.σ.dot(e.Δx) == 1 || e.neighbor.occupant->second.state.σ.dot(e.Δx) == -1) {
             o.push_back(e.neighbor.occupant);
           } 
         }
@@ -133,9 +135,10 @@ template<unsigned D> class System {
       return o;
     }
 
-    bool insert(const Particle<D>& p) {
+    bool insert(const Particle<D>& p, Rng& r) {
       if (overlaps(p).empty()) {
-        particles.push_back(p);
+        typename std::unordered_map<unsigned, Particle<D>>::iterator it;
+        std::tie(it, std::ignore) = particles.insert({r.uniform((unsigned)0, (unsigned)pow(2, 28)), p});
         p.position.occupant = std::prev(particles.end());
         return true;
       } else {
@@ -143,8 +146,8 @@ template<unsigned D> class System {
       }
     }
 
-    void remove(typename std::list<Particle<D>>::iterator ip) {
-      ip->position.occupant = particles.end();
+    void remove(typename std::unordered_map<unsigned, Particle<D>>::iterator ip) {
+      ip->second.position.occupant = particles.end();
       particles.erase(ip);
     }
 
@@ -179,11 +182,11 @@ template<unsigned D> class System {
     }
 
     bool randomSwap(Rng& r) {
-      return swap(r.pick(particles), r.pick(particles));
+      return swap(r.pick(particles).second, r.pick(particles).second);
     }
 
-    void setGroundState() {
-      for (typename std::list<Particle<D>>::iterator ip = particles.begin(); ip != particles.end(); ip++) {
+    void setGroundState(Rng& r) {
+      for (typename std::unordered_map<unsigned, Particle<D>>::iterator ip = particles.begin(); ip != particles.end(); ip++) {
         remove(ip);
       }
 
@@ -197,10 +200,10 @@ template<unsigned D> class System {
         if (a > 0) {
           if (a <= D) {
             State<D> s(a - 1, -1);
-            insert(Particle<D>(v, s));
+            insert(Particle<D>(v, s), r);
           } else if (D < a) {
             State<D> s(2 * D - a, 1);
-            insert(Particle<D>(v, s));
+            insert(Particle<D>(v, s), r);
           }
         }
       }
@@ -234,7 +237,7 @@ template<unsigned D> class System {
             while (true) {
               Vertex<D>& v = r.pick(vertices);
               if (isEmpty(v)) {
-                insert(Particle<D>(v, State<D>(r)));
+                insert(Particle<D>(v, State<D>(r)), r);
                 break;
               }
             }
@@ -334,7 +337,7 @@ int main() {
   }
   */
 
-  s.setGroundState();
+  s.setGroundState(r);
 
   for (unsigned n = 0; n < 10; n++) {
     s.sweepGrandCanonical(exp(1/0.15), r);