From 70b13f405bf861dd16e0ddb963293e5546599b3d Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 25 Jun 2021 10:57:46 +0200
Subject: Partial progress towards merging BM model into global objects.

---
 glass.hpp | 209 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 205 insertions(+), 4 deletions(-)

(limited to 'glass.hpp')

diff --git a/glass.hpp b/glass.hpp
index 64fe8b7..d313ced 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -180,6 +180,7 @@ public:
   const unsigned L;
   unsigned N;
   std::vector<Vertex<D, S>> vertices;
+  Transformation<D> orientation;
 
   unsigned vectorToIndex(const Vector<D>& x) const {
     unsigned i = 0;
@@ -197,7 +198,7 @@ public:
     return x;
   }
 
-  System(unsigned L, unsigned N) : L(L), N(N), vertices(iPow(L, D)) {
+  System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) {
     for (unsigned i = 0; i < iPow(L, D); i++) {
       vertices[i].position = indexToVector(i);
       vertices[i].initialPosition = vertices[i].position;
@@ -216,8 +217,14 @@ public:
     }
   }
 
+  unsigned size() const { return vertices.size(); }
+
   double density() const { return N / pow(L, D); }
 
+  unsigned maxOccupation() const {
+    return iPow(L, D);
+  }
+
   void setInitialPosition() {
     for (Vertex<D, S>& v : vertices) {
       v.initialPosition = v.position;
@@ -240,7 +247,7 @@ public:
   }
 };
 
-template <unsigned D, State S> class FiniteEnergySystem : public System<D, S> {
+template <unsigned D, State S> class SoftSystem : public System<D, S> {
 public:
   using System<D, S>::System;
 
@@ -326,7 +333,7 @@ public:
                        bool dry = false) {
     std::queue<std::array<std::reference_wrapper<Vertex<D, S>>, 2>> q;
     Vector<D> x0New = R.apply(v0.position);
-    Vertex<D, S>& v0New = this->vertices[this->vectorToIndex(x0New)];
+    Vertex<D, S>& v0New = SoftSystem::vertices[SoftSystem::vectorToIndex(x0New)];
     q.push({v0, v0New});
 
     unsigned n = 0;
@@ -346,7 +353,7 @@ public:
           Vertex<D, S>& vn = e.neighbor;
 
           Vector<D> xnNew = R.apply(vn.position);
-          Vertex<D, S>& vnNew = this->vertices[this->vectorToIndex(xnNew)];
+          Vertex<D, S>& vnNew = SoftSystem::vertices[SoftSystem::vectorToIndex(xnNew)];
 
           if (!vn.marked && !vnNew.marked) {
             double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state);
@@ -393,3 +400,197 @@ public:
   }
 };
 
+template <unsigned D, State S> class HardSystem : public System<D, S> {
+public:
+
+  HardSystem(unsigned L) : System<D, S>(L) {}
+
+  virtual std::list<std::reference_wrapper<Vertex<D, S>>> overlaps(Vertex<D, S>&, const S&,
+                                                        bool = false) { return {}; }
+
+  bool insert(Vertex<D, S>& v, const S& s) {
+    if (overlaps(v, s).empty()) {
+      v.state = s;
+      HardSystem::N++;
+      return true;
+    } else {
+      return false;
+    }
+  }
+
+  bool tryDeletion(Vertex<D, S>& v) {
+    if (v.empty()) {
+      return false;
+    } else {
+      v.state.remove();
+      HardSystem::N--;
+      return true;
+    }
+  }
+
+  bool tryRandomMove(Rng& r) {
+    Vertex<D, S>& v = r.pick(HardSystem::vertices);
+    S oldState = v.state;
+
+    if (!tryDeletion(v)) {
+      return false;
+    }
+
+    if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) {
+      for (HalfEdge<D, S>& e : v.adjacentEdges) {
+        if (1 == e.Δx.dot(oldState)) {
+          if (insert(e.neighbor, oldState.flip())) {
+            return true;
+          }
+          break;
+        }
+      }
+    } else {
+      S newState(r);
+      while (newState == oldState) {
+        newState = S(r);
+      }
+      if (insert(v, newState)) {
+        return true;
+      }
+    }
+    v.state = oldState;
+    HardSystem::N++;
+    return false;
+  }
+
+  bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& 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 {
+      return false;
+    }
+  }
+
+  bool tryRandomSwap(Rng& r) {
+    Vertex<D, S>& v1 = r.pick(HardSystem::vertices);
+    Vertex<D, S>& v2 = r.pick(HardSystem::vertices);
+
+    return trySwap(v1, v2);
+  }
+
+
+  bool compatible() {
+    for (Vertex<D, S>& v : HardSystem::vertices) {
+      if (overlaps(v, v.state, true).size() > 0) {
+        return false;
+      }
+    }
+
+    return true;
+  }
+
+  void sweepGrandCanonical(double z, Rng& r) {
+    for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
+      if (0.5 < r.uniform(0.0, 1.0)) {
+        double pIns = HardSystem::maxOccupation() * z / (HardSystem::N + 1);
+
+        if (pIns > r.uniform(0.0, 1.0)) {
+          while (true) {
+            Vertex<D, S>& v = r.pick(HardSystem::vertices);
+            if (v.empty()) {
+              insert(v, S(r));
+              break;
+            }
+          }
+        }
+      } else {
+
+        double pDel = HardSystem::N / (z * HardSystem::maxOccupation());
+
+        if (pDel > r.uniform(0.0, 1.0)) {
+          tryDeletion(r.pick(HardSystem::vertices));
+        }
+      }
+
+      tryRandomMove(r);
+    }
+  }
+
+  void sweepLocal(Rng& r) {
+    for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
+      tryRandomMove(r);
+    }
+  }
+
+  void sweepSwap(Rng& r) {
+    for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
+      tryRandomSwap(r);
+    }
+  }
+
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, Rng& r, bool dry = false) {
+    std::queue<std::reference_wrapper<Vertex<D, S>>> q;
+    q.push(v0);
+
+    unsigned n = 0;
+
+    while (!q.empty()) {
+      Vertex<D, S>& v = q.front();
+      q.pop();
+
+      if (!v.marked) {
+        Vector<D> xNew = R.apply(v.position);
+        Vertex<D, S>& vNew = HardSystem::vertices[HardSystem::vectorToIndex(xNew)];
+
+        v.marked = true;
+        vNew.marked = true;
+
+        S s = R.apply(v.state);
+        S sNew = R.apply(vNew.state);
+
+        std::list<std::reference_wrapper<Vertex<D, S>>> overlaps1 = overlaps(vNew, s, true);
+        std::list<std::reference_wrapper<Vertex<D, S>>> overlaps2 = overlaps(v, sNew, true);
+        overlaps1.splice(overlaps1.begin(), overlaps2);
+
+        for (Vertex<D, S>& vn : overlaps1) {
+          if (!vn.marked) {
+            q.push(vn);
+          }
+        }
+
+        if (!dry) {
+          v.state = sNew;
+          vNew.state = s;
+        }
+
+        n += 1;
+      }
+    }
+
+    return n;
+  }
+
+  void swendsenWang(const Transformation<D>& R, Rng& r) {
+    for (Vertex<D, S>& v : HardSystem::vertices) {
+      if (!v.marked) {
+        bool dry = 0.5 < r.uniform(0.0, 1.0);
+        unsigned n = flipCluster(R, v, r, dry);
+        if (n > pow(HardSystem::L, D) / 4 && !dry) {
+          orientation = R.apply(orientation);
+        }
+      }
+    }
+
+    for (Vertex<D, S>& v : HardSystem::vertices) {
+      v.marked = false;
+    }
+  }
+
+  int overlap(const System<D, S>& s) const {
+    int o = 0;
+
+    for (unsigned i = 0; i < HardSystem::vertices.size(); i++) {
+      S s2 = orientation.apply(s.vertices[HardSystem::vectorToIndex(orientation.inverse().apply(HardSystem::indexToVector(i)))].state);
+      o += HardSystem::vertices[i].state.dot(s2);
+    }
+
+    return o;
+  }
+};
-- 
cgit v1.2.3-54-g00ecf