Everything compiles.

2 files changed, 93 insertions, 125 deletions

diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index 5b6213f..52a4360 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -24,41 +24,9 @@ public:
   void remove() {
     maximumNeighbors = Empty;
   }
-}
-
-template <int D> class Vertex {
-public:
-  Vector<D> initialPosition;
-  Vector<D> position;
-  std::vector<std::reference_wrapper<Vertex<D>>> neighbors;
-  unsigned occupiedNeighbors;
-  unsigned maximumNeighbors;
-  bool marked;
-  bool visited;
-
-  Vertex() {
-    marked = false;
-  }
-
-  bool empty() const { return maximumNeighbors == Empty; }
-  bool frustrated() const { return occupiedNeighbors > maximumNeighbors; }
-  bool valid() const {
-    unsigned o = 0;
-
-    for (const Vertex<D>& v : neighbors) {
-      if (!v.empty()) {
-        o++;
-      }
-    }
-
-    return o == occupiedNeighbors;
-  }
-  void setInitial() {
-    initialPosition = position;
-  }
 };
 
-template <int D> class BiroliSystem : public HardSystem<D, BiroliState> {
+template <int D> class BiroliSystem : public System<D, BiroliState> {
 public:
   std::vector<unsigned> N;
   std::unordered_set<Vertex<D, BiroliState>*> occupants;
@@ -69,11 +37,23 @@ public:
     return BiroliSystem::size() - N[0];
   }
 
-  double density() const { return (double)occupancy() / size(); }
+  double density() const { return (double)occupancy() / BiroliSystem::size(); }
+
+  bool valid(const Vertex<D, BiroliState>& v) const {
+    unsigned o = 0;
+
+    for (const Vertex<D, BiroliState>& vn : v.neighbors) {
+      if (!vn.empty()) {
+        o++;
+      }
+    }
+
+    return o == v.occupiedNeighbors;
+  }
 
   bool compatible() const {
-    for (const Vertex<D>& v : vertices) {
-      if (v.frustrated()) {
+    for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
+      if (v.state.frustrated()) {
         return false;
       }
     }
@@ -82,8 +62,8 @@ public:
   }
 
   bool valid() const {
-    for (const Vertex<D>& v : vertices) {
-      if (!v.valid()) {
+    for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
+      if (!valid(v)) {
         return false;
       }
     }
@@ -91,25 +71,12 @@ public:
     return true;
   }
 
-  BiroliSystem(unsigned L, unsigned T) : L(L), N(T + 1, 0), vertices(iPow(L, D)), orientation(L) {
-    N[0] = size();
-
-    for (unsigned i = 0; i < size(); i++) {
-      vertices[i].position = indexToVector(i);
-      vertices[i].neighbors.reserve(2 * D);
-    }
-
-    for (unsigned d = 0; d < D; d++) {
-      for (signed i = 0; i < size(); i++) {
-        unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1));
-        vertices[i].neighbors.push_back(vertices[j]);
-        vertices[j].neighbors.push_back(vertices[i]);
-      }
-    }
+  BiroliSystem(unsigned L, unsigned T) : System<D, BiroliState>(L), N(T + 1, 0) {
+    N[0] = BiroliSystem::size();
   }
 
-  std::list<std::reference_wrapper<Vertex<D>>> overlaps(Vertex<D>& v) {
-    std::list<std::reference_wrapper<Vertex<D>>> o;
+  std::list<std::reference_wrapper<Vertex<D, BiroliState>>> overlaps(Vertex<D, BiroliState>& v) {
+    std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o;
 
     if (v.empty()) { // an empty site cannot be frustrating anyone
       return o;
@@ -118,7 +85,7 @@ public:
     bool selfFrustrated = v.frustrated();
     bool anyNeighborFrustrated = false;
 
-    for (Vertex<D>& vn : v.neighbors) {
+    for (Vertex<D, BiroliState>& vn : v.neighbors) {
       if (!vn.empty()) {
         bool thisNeighborFrustrated = vn.frustrated();
 
@@ -139,28 +106,29 @@ public:
     return o;
   }
 
-  std::list<std::reference_wrapper<Vertex<D>>> overlaps(const std::list<std::reference_wrapper<Vertex<D>>>& vs) {
-    std::list<std::reference_wrapper<Vertex<D>>> o;
+  std::list<std::reference_wrapper<Vertex<D, BiroliState>>> overlaps(const std::list<std::reference_wrapper<Vertex<D, BiroliState>>>& vs) {
+    std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o;
 
-    for (Vertex<D>& v : vs) {
+    for (Vertex<D, BiroliState>& v : vs) {
       o.splice(o.begin(), overlaps(v));
     }
 
     return o;
   }
 
-  bool canReplaceWith(const Vertex<D>& v, unsigned t) const {
+  bool canReplaceWith(const Vertex<D, BiroliState>& v, unsigned t) const {
     if (t == Empty) {
       return true;
     }
 
-    if (t < v.occupiedNeighbors) {
+    if (t < v.state.occupiedNeighbors) {
       return false;
     }
 
     if (v.empty()) {
-      for (const Vertex<D>& vn : v.neighbors) {
-        if (vn.maximumNeighbors < vn.occupiedNeighbors + 1) {
+      for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+        const Vertex<D, BiroliState>& vn = e.neighbor;
+        if (vn.state.maximumNeighbors < vn.state.occupiedNeighbors + 1) {
           return false;
         }
       }
@@ -170,16 +138,16 @@ public:
   }
 
   void setInitial() {
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
       v.setInitial();
     }
   }
 
-  bool insert(Vertex<D>& v, unsigned t, bool force = false) {
+  bool insert(Vertex<D, BiroliState>& v, unsigned t, bool force = false) {
     if (force || (v.empty() && canReplaceWith(v, t))) {
-      v.maximumNeighbors = t;
-      for (Vertex<D>& vn : v.neighbors) {
-        vn.occupiedNeighbors++;
+      v.state.maximumNeighbors = t;
+      for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+        e.neighbor.state.occupiedNeighbors++;
       }
       N[t]++;
       N[0]--;
@@ -192,39 +160,39 @@ public:
     }
   }
 
-  bool remove(Vertex<D>& v) {
+  bool remove(Vertex<D, BiroliState>& v) {
     if (v.empty()) {
       return false;
     } else {
-      if (v.maximumNeighbors > 0) {
+      if (v.state.maximumNeighbors > 0) {
         occupants.erase(&v);
       }
-      N[v.maximumNeighbors]--;
+      N[v.state.maximumNeighbors]--;
       N[0]++;
-      v.maximumNeighbors = Empty;
-      for (Vertex<D>& vn : v.neighbors) {
-        vn.occupiedNeighbors--;
+      v.state.maximumNeighbors = Empty;
+      for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+        e.neighbor.state.occupiedNeighbors--;
       }
       return true;
     }
   }
 
   bool randomLocalExchange(Rng& r) {
-    Vertex<D>& v = r.pick(vertices);
+    Vertex<D, BiroliState>& v = r.pick(BiroliSystem::vertices);
 
     return exchange(v, r.pick(v.neighbors));
   }
 
-  void swap(Vertex<D>& v1, Vertex<D>& v2) {
-    if (v1.maximumNeighbors != v2.maximumNeighbors) {
+  void swap(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) {
+    if (v1.state.maximumNeighbors != v2.state.maximumNeighbors) {
       if (!v1.empty() && !v2.empty()) {
-        std::swap(v1.maximumNeighbors, v2.maximumNeighbors);
+        std::swap(v1.state.maximumNeighbors, v2.state.maximumNeighbors);
       } else if (v1.empty()) {
-        unsigned t = v2.maximumNeighbors;
+        unsigned t = v2.state.maximumNeighbors;
         remove(v2);
         insert(v1, t, true);
       } else {
-        unsigned t = v1.maximumNeighbors;
+        unsigned t = v1.state.maximumNeighbors;
         remove(v1);
         insert(v2, t, true);
       }
@@ -232,8 +200,8 @@ public:
     }
   }
 
-  bool exchange(Vertex<D>& v1, Vertex<D>& v2) {
-    if (canReplaceWith(v1, v2.maximumNeighbors) && canReplaceWith(v2, v1.maximumNeighbors)) {
+  bool exchange(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) {
+    if (canReplaceWith(v1, v2.state.maximumNeighbors) && canReplaceWith(v2, v1.state.maximumNeighbors)) {
       swap(v1, v2);
       return true;
     } else {
@@ -246,10 +214,10 @@ public:
 
     using namespace std::complex_literals;
 
-    for (const Vertex<D>& v : vertices) {
+    for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
       if (v.maximumNeighbors == t) {
         for (unsigned d = 0; d < D; d++) {
-          F += exp(1i * 2.0 * M_PI * (double)k * (double)(v.position(d) - v.initialPosition(d)) / (double)L);
+          F += exp(1i * 2.0 * M_PI * (double)k * (double)(v.position(d) - v.initialPosition(d)) / (double)BiroliSystem::L);
         }
       }
     }
@@ -259,11 +227,11 @@ public:
     return F;
   }
 
-  double overlap(const System<D>& s, unsigned t = 2) const {
+  double overlap(const BiroliSystem<D>& s, unsigned t = 2) const {
     int o = 0;
 
-    for (const Vertex<D>& v1 : vertices) {
-      const Vertex<D>& v2 = s.vertices[vectorToIndex(orientation.inverse().apply(v1.position))];
+    for (const Vertex<D, BiroliState>& v1 : BiroliSystem::vertices) {
+      const Vertex<D, BiroliState>& v2 = s.BiroliSystem::vertices[vectorToIndex(BiroliSystem::orientation.inverse().apply(v1.position))];
       if (v1.maximumNeighbors == v2.maximumNeighbors) {
         o++;
       }
@@ -275,12 +243,12 @@ public:
       nn += iPow(n, 2);
     }
 
-    return ((double)o / (double)size() - (double)nn / pow((double)size(), 2)) / (1 - (double)nn / pow((double)size(), 2));
+    return ((double)o / (double)BiroliSystem::size() - (double)nn / pow((double)BiroliSystem::size(), 2)) / (1 - (double)nn / pow((double)BiroliSystem::size(), 2));
   }
 
   bool randomExchange(Rng& r) {
-    Vertex<D>& v1 = r.pick(vertices);
-    Vertex<D>& v2 = r.pick(vertices);
+    Vertex<D, BiroliState>& v1 = r.pick(BiroliSystem::vertices);
+    Vertex<D, BiroliState>& v2 = r.pick(BiroliSystem::vertices);
 
     return exchange(v1, v2);
   }
@@ -290,24 +258,24 @@ public:
   void stepGrandCanonical(double z, Rng& r) {
     if (1.0 / 11.0 < r.uniform(0.0, 1.0)) {
       unsigned t = r.pick({1, 2, 2, 2, 2, 2, 3, 3, 3, 3});
-      double pIns = size() * z / (occupancy() + 1);
+      double pIns = BiroliSystem::size() * z / (occupancy() + 1);
 
       if (pIns > r.uniform(0.0, 1.0)) {
-        Vertex<D>& v = r.pick(vertices);
+        Vertex<D, BiroliState>& v = r.pick(BiroliSystem::vertices);
         insert(v, t);
       }
     } else {
       double pDel = density() / z;
 
       if (pDel > r.uniform(0.0, 1.0)) {
-        Vertex<D>* v = r.pick(occupants);
+        Vertex<D, BiroliState>* v = r.pick(occupants);
         remove(*v);
       }
     }
   }
 
   void sweepGrandCanonical(double z, Rng& r) {
-    for (unsigned i = 0; i < size(); i++) {
+    for (unsigned i = 0; i < BiroliSystem::size(); i++) {
       stepGrandCanonical(z, r);
 
       randomExchange(r);
@@ -315,36 +283,36 @@ public:
   }
 
   void sweepLocal(Rng& r) {
-    for (unsigned i = 0; i < size() / 2; i++) {
+    for (unsigned i = 0; i < BiroliSystem::size() / 2; i++) {
       randomLocalExchange(r);
     }
   }
 
   void sweepSwap(Rng& r) {
-    for (unsigned i = 0; i < size() / 2; i++) {
+    for (unsigned i = 0; i < BiroliSystem::size() / 2; i++) {
       randomExchange(r);
     }
   }
 
-  Vertex<D>& apply(const Transformation<D>& R, const Vertex<D>& v) {
-    return vertices[vectorToIndex(R.apply(v.position))];
+  Vertex<D, BiroliState>& apply(const Transformation<D>& R, const Vertex<D, BiroliState>& v) {
+    return BiroliSystem::vertices[vectorToIndex(R.apply(v.position))];
   }
 
-  bool eventChain(const Transformation<D>& R, Vertex<D>& v0) {
+  bool eventChain(const Transformation<D>& R, Vertex<D, BiroliState>& v0) {
     unsigned n = 0;
 
     if (!v0.empty()) {
-      Vertex<D>& v0New = apply(R, v0);;
+      Vertex<D, BiroliState>& v0New = apply(R, v0);;
       unsigned t0 = v0.maximumNeighbors;
 
-      std::queue<std::pair<std::reference_wrapper<Vertex<D>>, unsigned>> q;
+      std::queue<std::pair<std::reference_wrapper<Vertex<D, BiroliState>>, unsigned>> q;
 
       q.push({v0New, t0});
       remove(v0);
 
       while (!q.empty()) {
         auto [vr, t] = q.front();
-        Vertex<D>& v = vr;
+        Vertex<D, BiroliState>& v = vr;
         q.pop();
 
         if (!v.empty()) {
@@ -355,7 +323,7 @@ public:
         insert(v, t, true);
         v.marked = true;
 
-        for (Vertex<D>& vn : overlaps({v})) {
+        for (Vertex<D, BiroliState>& vn : overlaps({v})) {
           if (!vn.marked) {
             q.push({apply(R, vn), vn.maximumNeighbors});
             remove(vn);
@@ -365,7 +333,7 @@ public:
              * frustrated, it may be due to one of its neighbors!
              */
             if (&vn != &v) {
-              for (Vertex<D>& vnn : overlaps(vn)) {
+              for (Vertex<D, BiroliState>& vnn : overlaps(vn)) {
                 if (!vnn.marked) {
                   q.push({apply(R, vnn), vnn.maximumNeighbors});
                   remove(vnn);
@@ -378,7 +346,7 @@ public:
       }
     }
 
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
       v.marked = false;
     }
 
@@ -386,44 +354,44 @@ public:
   }
 
   bool randomEventChain(Rng& r) {
-    Transformation<D> R(L);
+    Transformation<D> R(BiroliSystem::L);
     R.v[r.uniform((unsigned)0, (unsigned)D-1)] = r.pick({-1, 1});
-    return eventChain(R, r.pick(vertices));
+    return eventChain(R, r.pick(BiroliSystem::vertices));
   }
 
-  unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0) {
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D, BiroliState>& v0) {
     unsigned n = 0;
 
-    Vertex<D>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
+    Vertex<D, BiroliState>& v0New = BiroliSystem::vertices[vectorToIndex(R.apply(v0.position))];
 
     if (&v0New != &v0) {
-      std::queue<std::reference_wrapper<Vertex<D>>> q;
+      std::queue<std::reference_wrapper<Vertex<D, BiroliState>>> q;
 
       v0.marked = true;
       v0New.marked = true;
 
       swap(v0, v0New);
 
-      for (Vertex<D>& vn : overlaps({v0, v0New})) {
+      for (Vertex<D, BiroliState>& vn : overlaps({v0, v0New})) {
         if (!vn.marked) {
           q.push(vn);
         }
       }
 
       while (!q.empty()) {
-        Vertex<D>& v = q.front();
+        Vertex<D, BiroliState>& v = q.front();
         q.pop();
 
         if (!v.marked && !overlaps(v).empty()) {
           v.marked = true;
-          Vertex<D>& vNew = vertices[vectorToIndex(R.apply(v.position))];
+          Vertex<D, BiroliState>& vNew = BiroliSystem::vertices[vectorToIndex(R.apply(v.position))];
 
           if (&vNew != &v) {
             vNew.marked = true;
 
             swap(v, vNew);
 
-            for (Vertex<D>& vn : overlaps({v, vNew})) {
+            for (Vertex<D, BiroliState>& vn : overlaps({v, vNew})) {
               if (!vn.marked) {
                 q.push(vn);
               } else {
@@ -431,7 +399,7 @@ public:
                  * frustrated, it may be due to one of its neighbors!
                  */
                 if (&vn != &v && &vn != &vNew) {
-                  for (Vertex<D>& vnn : overlaps(vn)) {
+                  for (Vertex<D, BiroliState>& vnn : overlaps(vn)) {
                     if (!vnn.marked) {
                       q.push(vnn);
                     }
@@ -449,11 +417,11 @@ public:
       n = 1;
     }
 
-    if (n > size() / 2) {
-      orientation = R.apply(orientation);
+    if (n > BiroliSystem::size() / 2) {
+      BiroliSystem::orientation = R.apply(BiroliSystem::orientation);
     }
 
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
       v.marked = false;
     }
 
@@ -461,10 +429,10 @@ public:
   }
 };
 
-void print(const System<2>& s) {
-  for (const Vertex<2>& v : s.vertices) {
+void print(const BiroliSystem<2>& s) {
+  for (const Vertex<2, BiroliState>& v : s.vertices) {
     if (!v.empty()) {
-      std::cerr << v.maximumNeighbors;
+      std::cerr << v.state.maximumNeighbors;
     } else {
       std::cerr << " ";
     }
@@ -484,14 +452,14 @@ int main() {
   double μmax = 10;
   double dμ = 0.05;
 
-  System<D> s(L, 3);
+  BiroliSystem<D> s(L, 3);
 
   Rng r;
 
   for (double μ = μmin; μ <= μmax; μ += dμ) {
     double z = exp(μ);
     for (unsigned i = 0; i < 1e4; i++) {
-      for (unsigned j = 0; j < s.size(); j++) {
+      for (unsigned j = 0; j < s.BiroliSystem::size(); j++) {
         s.stepGrandCanonical(z, r);
 //        s.randomEventChain(r);
         s.randomExchange(r);
diff --git a/glass.hpp b/glass.hpp
index d313ced..fc7b186 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -573,7 +573,7 @@ public:
         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);
+          HardSystem::orientation = R.apply(HardSystem::orientation);
         }
       }
     }
@@ -587,7 +587,7 @@ public:
     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);
+      S s2 = HardSystem::orientation.apply(s.vertices[HardSystem::vectorToIndex(HardSystem::orientation.inverse().apply(HardSystem::indexToVector(i)))].state);
       o += HardSystem::vertices[i].state.dot(s2);
     }