1 files changed, 18 insertions, 160 deletions
diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index e2fb71a..5b6213f 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -1,149 +1,30 @@
 #include <fstream>
 #include <iostream>
-#include <limits>
-#include <list>
-#include <queue>
-#include <vector>
-#include <chrono>
 #include <unordered_set>
 
-#include <eigen3/Eigen/Dense>
-
-#include "pcg-cpp/include/pcg_random.hpp"
-#include "randutils/randutils.hpp"
-
-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>;
+#include "glass.hpp"
 
 const unsigned Empty = std::numeric_limits<unsigned>::max();
 
-int iPow(int x, unsigned p) {
-  if (p == 0)
-    return 1;
-  if (p == 1)
-    return x;
-
-  int tmp = iPow(x, p / 2);
-  if (p % 2 == 0)
-    return tmp * tmp;
-  else
-    return x * tmp * tmp;
-}
-
-unsigned mod(signed a, unsigned b) { return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b; }
-
-template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) {
-  Vector<D> u;
-  for (unsigned i = 0; i < D; i++) {
-    u(i) = mod(v(i), b);
-  }
-  return u;
-}
-
-template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
-  if (level > 0) {
-    unsigned new_level = level - 1;
-    for (std::array<double, D>& sequence : sequences) {
-      std::array<double, D> new_sequence = sequence;
-      new_sequence[new_level] = -1;
-      sequences.push_front(new_sequence);
-    }
-    one_sequences<D>(sequences, new_level);
-  }
-}
-
-template <int D> std::vector<Matrix<D>> generateTorusMatrices() {
-  std::vector<Matrix<D>> mats;
-
-  std::array<double, D> ini_sequence;
-  ini_sequence.fill(1);
-  std::list<std::array<double, D>> sequences;
-  sequences.push_back(ini_sequence);
-
-  one_sequences<D>(sequences, D);
-
-  sequences.pop_back(); // don't want the identity matrix!
-
-  for (std::array<double, D> sequence : sequences) {
-    Matrix<D> m;
-    for (unsigned i = 0; i < D; i++) {
-      for (unsigned j = 0; j < D; j++) {
-        if (i == j) {
-          m(i, j) = sequence[i];
-        } else {
-          m(i, j) = 0;
-        }
-      }
-    }
-
-    mats.push_back(m);
-  }
-
-  for (unsigned i = 0; i < D; i++) {
-    for (unsigned j = 0; j < D; j++) {
-      if (i != j) {
-        Matrix<D> m;
-        for (unsigned k = 0; k < D; k++) {
-          for (unsigned l = 0; l < D; l++) {
-            if ((k == i && l == j) || (k == j && l == i)) {
-              if (i < j) {
-                m(k, l) = 1;
-              } else {
-                m(k, l) = -1;
-              }
-            } else if (k == l && (k != i && k != j)) {
-              m(k, l) = 1;
-            } else {
-              m(k, l) = 0;
-            }
-          }
-        }
-        mats.push_back(m);
-      }
-    }
-  }
-
-  return mats;
-}
-
-template <int D> class Transformation {
+class BiroliState {
 public:
-  unsigned L;
-  Matrix<D> m;
-  Vector<D> v;
-
-  Transformation(unsigned L) : L(L) {
-    m.setIdentity();
-    v.setZero();
-  }
-
-  Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {}
-
-  Transformation(unsigned L, const std::vector<Matrix<D>>& ms, Rng& r) : m(r.pick(ms)), L(L) {
-    for (unsigned i = 0; i < D; i++) {
-      v[i] = r.uniform((unsigned)0, L - 1);
-    }
+  unsigned maximumNeighbors;
+  unsigned occupiedNeighbors;
 
-    v = v - m * v;
+  BiroliState() {
+    maximumNeighbors = Empty;
+    occupiedNeighbors = 0;
   }
 
-  Transformation<D> inverse() const {
-    return Transformation<D>(L, m.transpose(), -m.transpose() * v);
+  bool empty() const { return maximumNeighbors == Empty; }
+  bool frustrated() const { return occupiedNeighbors > maximumNeighbors; }
+  bool operator==(const BiroliState& s) const {
+    return s.maximumNeighbors == maximumNeighbors;
   }
-
-  Vector<D> apply(const Vector<D>& x) const { return mod<D>(v + m * x, L); }
-
-  Transformation<D> apply(const Transformation<D>& t) const {
-    Transformation<D> tNew(L);
-
-    tNew.m = m * t.m;
-    tNew.v = apply(t.v);
-
-    return tNew;
+  void remove() {
+    maximumNeighbors = Empty;
   }
-};
+}
 
 template <int D> class Vertex {
 public:
@@ -156,8 +37,6 @@ public:
   bool visited;
 
   Vertex() {
-    occupiedNeighbors = 0;
-    maximumNeighbors = Empty;
     marked = false;
   }
 
@@ -179,20 +58,15 @@ public:
   }
 };
 
-template <int D> class System {
+template <int D> class BiroliSystem : public HardSystem<D, BiroliState> {
 public:
-  const unsigned L;
   std::vector<unsigned> N;
-  std::unordered_set<Vertex<D>*> occupants;
-  std::vector<Vertex<D>> vertices;
-  Transformation<D> orientation;
-
-  unsigned size() const { return vertices.size(); }
+  std::unordered_set<Vertex<D, BiroliState>*> occupants;
 
   unsigned types() const { return N.size() - 1; }
 
   unsigned occupancy() const {
-    return size() - N[0];
+    return BiroliSystem::size() - N[0];
   }
 
   double density() const { return (double)occupancy() / size(); }
@@ -217,23 +91,7 @@ public:
     return true;
   }
 
-  unsigned vectorToIndex(const Vector<D>& x) const {
-    unsigned i = 0;
-    for (unsigned d = 0; d < D; d++) {
-      i += x[d] * iPow(L, d);
-    }
-    return i;
-  }
-
-  Vector<D> indexToVector(unsigned i) const {
-    Vector<D> x;
-    for (unsigned d = 0; d < D; d++) {
-      x[d] = (i / iPow(L, d)) % L;
-    }
-    return x;
-  }
-
-  System(unsigned L, unsigned T) : L(L), N(T + 1, 0), vertices(iPow(L, D)), orientation(L) {
+  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++) {