Split up functions, tired to generize in preparation for implementing Biroli-Mezard.

1 files changed, 502 insertions, 0 deletions

diff --git a/glass.hpp b/glass.hpp
new file mode 100644
index 0000000..da06542
--- /dev/null
+++ b/glass.hpp
@@ -0,0 +1,502 @@
+#include <list>
+#include <queue>
+#include <vector>
+
+#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>;
+
+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;
+    unsigned old_length = sequences.size();
+    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 CiamarraState : public Vector<D> {
+public:
+  CiamarraState() : Vector<D>(Vector<D>::Zero()) {}
+
+  CiamarraState(const Vector<D>& v) : Vector<D>(v) {}
+
+  CiamarraState(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; }
+
+  CiamarraState(Rng& r) : CiamarraState(r.uniform(0, D - 1), r.pick({-1, 1})) {}
+
+  bool isEmpty() const { return this->squaredNorm() == 0; }
+
+  void remove() { this->setZero(); }
+
+  CiamarraState<D> flip() const {
+    CiamarraState<D> s;
+    for (unsigned i = 0; i < D; i++) {
+      s(i) = -this->operator()(i);
+    }
+    return s;
+  }
+};
+
+class BiroliState {
+public:
+  unsigned type;
+  unsigned occupiedNeighbors;
+
+  BiroliState() {
+    type = 0;
+    occupiedNeighbors = 0;
+  }
+
+  BiroliState(unsigned t, unsigned o) {
+    type = t;
+    occupiedNeighbors = o;
+  }
+
+  bool isEmpty() const { return type == 0; }
+
+  void remove() { type = 0; };
+};
+
+template <int D> class Transformation {
+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);
+    }
+
+    v = v - m * v;
+  }
+
+  Transformation<D> inverse() const {
+    return Transformation<D>(L, m.transpose(), -m.transpose() * v);
+  }
+
+  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;
+  }
+
+  CiamarraState<D> apply(const CiamarraState<D>& s) const { return CiamarraState<D>(m * s); }
+
+  BiroliState apply(const BiroliState& s) const { return s; }
+};
+
+template <int D, typename State> class HalfEdge;
+
+template <int D, typename State> class Vertex {
+public:
+  Vector<D> position;
+  std::vector<HalfEdge<D, State>> adjacentEdges;
+  State state;
+  bool marked;
+
+  bool isEmpty() const { return state.isEmpty(); }
+};
+
+template <int D, class State> class HalfEdge {
+public:
+  Vertex<D, State>& neighbor;
+  Vector<D> Δx;
+
+  HalfEdge(Vertex<D, State>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
+};
+
+template <int D, class State> class System {
+public:
+  const unsigned L;
+  unsigned N;
+  std::vector<Vertex<D, State>> vertices;
+  Transformation<D> orientation;
+
+  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) : 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].adjacentEdges.reserve(2 * D);
+      vertices[i].marked = false;
+    }
+
+    for (unsigned d = 0; d < D; d++) {
+      Vector<D> Δx = Vector<D>::Zero();
+      Δx[d] = 1;
+      for (signed i = 0; i < iPow(L, D); i++) {
+        unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1));
+        vertices[i].adjacentEdges.push_back(HalfEdge<D, State>(vertices[j], Δx));
+        vertices[j].adjacentEdges.push_back(HalfEdge<D, State>(vertices[i], -Δx));
+      }
+    }
+  }
+
+  virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&,
+                                                                       const State&, bool = false) { return {}; };
+
+  virtual bool insert(Vertex<D, State>& v, const State& s) { return false; };
+
+  virtual bool remove(Vertex<D, State>& v) { return false;};
+
+  virtual bool tryRandomMove(Rng& r) { return false;};
+
+  virtual bool swap(Vertex<D, State>& v1, Vertex<D, State>& v2) { return false;};
+
+  bool tryRandomSwap(Rng& r) {
+    Vertex<D, State>& v1 = r.pick(vertices);
+    Vertex<D, State>& v2 = r.pick(vertices);
+
+    return swap(v1, v2);
+  }
+
+  bool compatible() {
+    for (Vertex<D, State>& v : vertices) {
+      if (overlaps(v, v.state, true).size() > 0) {
+        return false;
+      }
+    }
+
+    return true;
+  }
+
+  double density() const { return N / pow(L, D); }
+
+  unsigned size() const { return vertices.size(); }
+
+  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 = size() * z / (N + 1);
+
+        if (pIns > r.uniform(0.0, 1.0)) {
+          while (true) {
+            Vertex<D, State>& v = r.pick(vertices);
+            if (v.isEmpty()) {
+              insert(v, State(r));
+              break;
+            }
+          }
+        }
+      } else {
+
+        double pDel = N / (z * size());
+
+        if (pDel > r.uniform(0.0, 1.0)) {
+          remove(r.pick(vertices));
+        }
+      }
+
+      tryRandomMove(r);
+    }
+  }
+
+  void sweepLocal(Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      tryRandomMove(r);
+    }
+  }
+
+  void sweepSwap(Rng& r) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      tryRandomSwap(r);
+    }
+  }
+
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, Rng& r, bool dry = false) {
+    std::queue<std::reference_wrapper<Vertex<D, State>>> q;
+    q.push(v0);
+
+    unsigned n = 0;
+
+    while (!q.empty()) {
+      Vertex<D, State>& v = q.front();
+      q.pop();
+
+      if (!v.marked) {
+        Vector<D> xNew = R.apply(v.position);
+        Vertex<D, State>& vNew = vertices[vectorToIndex(xNew)];
+
+        v.marked = true;
+        vNew.marked = true;
+
+        CiamarraState<D> s = R.apply(v.state);
+        CiamarraState<D> sNew = R.apply(vNew.state);
+
+        std::list<std::reference_wrapper<Vertex<D, State>>> overlaps1 = overlaps(vNew, s, true);
+        std::list<std::reference_wrapper<Vertex<D, State>>> overlaps2 = overlaps(v, sNew, true);
+        overlaps1.splice(overlaps1.begin(), overlaps2);
+
+        for (Vertex<D, State>& vn : overlaps1) {
+          if (!vn.marked) {
+            q.push(vn);
+          }
+        }
+
+        if (!dry) {
+          swap(v, vNew);
+        }
+
+        n += 1;
+      }
+    }
+
+    return n;
+  }
+
+  void swendsenWang(const Transformation<D>& R, Rng& r) {
+    for (Vertex<D, State>& v : vertices) {
+      if (!v.marked) {
+        bool dry = 0.5 < r.uniform(0.0, 1.0);
+        unsigned n = flipCluster(R, v, r, dry);
+        if (n > pow(L, D) / 4 && !dry) {
+          orientation = R.apply(orientation);
+        }
+      }
+    }
+
+    for (Vertex<D, State>& v : vertices) {
+      v.marked = false;
+    }
+  }
+
+  virtual int overlap(const System<D, State>& s) const { return 0; };
+};
+
+template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
+  public:
+
+    using System<D, CiamarraState<D>>::System;
+  std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>>
+  overlaps(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s,
+           bool excludeSelf = false) override {
+    std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>> o;
+
+    if (s.isEmpty()) {
+      return o;
+    }
+
+    if (!v.isEmpty() && !excludeSelf) {
+      o.push_back(v);
+    }
+
+    for (const HalfEdge<D, CiamarraState<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);
+        }
+      }
+    }
+
+    return o;
+  }
+
+  bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s) override {
+    if (overlaps(v, s).empty()) {
+      v.state = s;
+      this->N++;
+      return true;
+    } else {
+      return false;
+    }
+  }
+
+  bool remove(Vertex<D, CiamarraState<D>>& v) override {
+    if (v.isEmpty()) {
+      return false;
+    } else {
+      v.state.remove();
+      this->N--;
+      return true;
+    }
+  }
+
+  bool tryRandomMove(Rng& r) override {
+    Vertex<D, CiamarraState<D>>& v = r.pick(this->vertices);
+    CiamarraState<D> oldState = v.state;
+
+    if (!remove(v)) {
+      return false;
+    }
+
+    if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) {
+      for (HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) {
+        if (1 == e.Δx.dot(oldState)) {
+          if (insert(e.neighbor, oldState.flip())) {
+            return true;
+          }
+          break;
+        }
+      }
+    } else {
+      CiamarraState<D> newState(r);
+      while (newState.dot(oldState) == 1) {
+        newState = CiamarraState<D>(r);
+      }
+      if (insert(v, newState)) {
+        return true;
+      }
+    }
+    v.state = oldState;
+    this->N++;
+    return false;
+  }
+
+  bool swap(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override {
+    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;
+    }
+  }
+
+  void setGroundState() {
+    this->N = 0;
+
+    for (Vertex<D, CiamarraState<D>>& v : this->vertices) {
+      unsigned a = 0;
+      for (unsigned d = 0; d < D; d++) {
+        a += (d + 1) * v.position(d);
+      }
+      a %= 2 * D + 1;
+
+      v.state.setZero() = Vector<D>::Zero();
+
+      if (0 < a && a <= D) {
+        v.state(a - 1) = -1;
+        this->N++;
+      } else if (D < a) {
+        v.state(2 * D - a) = 1;
+        this->N++;
+      }
+    }
+  }
+
+  int overlap(const System<D, CiamarraState<D>>& s) const override {
+    int o = 0;
+
+    for (unsigned i = 0; i < this->vertices.size(); i++) {
+      CiamarraState<D> s2 = this->orientation.apply(
+          s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state);
+      o += this->vertices[i].state.dot(s2);
+    }
+
+    return o;
+  }
+};