Started splitting up shared functions into another file.

2 files changed, 214 insertions, 211 deletions

diff --git a/distinguishable.cpp b/distinguishable.cpp
index d9f51ea..4b66989 100644
--- a/distinguishable.cpp
+++ b/distinguishable.cpp
@@ -1,192 +1,25 @@
-#include <eigen3/Eigen/Dense>
-#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
 #include <iostream>
-#include <list>
-#include <vector>
 #include <queue>
 
-#include "pcg-cpp/include/pcg_random.hpp"
-#include "randutils/randutils.hpp"
+#include "glass.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 <unsigned 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;
-}
-
-class State {
+class DistinguishableState {
 public:
   unsigned type;
 
-  State() {
-    type = 0;
-  }
-
-  State(unsigned t) {
-    type = t;
-  }
+  DistinguishableState() { type = 0; }
 
-  bool isEmpty() const { return type == 0; }
+  DistinguishableState(unsigned t) { type = t; }
 
+  bool empty() const { return type == 0; }
   void remove() { type = 0; }
 };
 
-template <unsigned 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;
-    }
-
-    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;
-    }
-
-    Transformation<D> inverse() const {
-      return Transformation<D>(L, m.transpose(), -m.transpose() * v);
-    }
-};
-
-template <unsigned D> class HalfEdge;
-
-template <unsigned D> class Vertex {
-public:
-  Vector<D> position;
-  Vector<D> initialPosition;
-  State state;
-  std::vector<HalfEdge<D>> adjacentEdges;
-  bool marked;
-
-  bool isEmpty() const { return state.isEmpty(); }
-};
-
-template <unsigned D> class HalfEdge {
-public:
-  Vertex<D>& neighbor;
-  Vector<D> Δx;
-
-  HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
-};
-
 template <unsigned D> class System {
 public:
   const unsigned L;
   unsigned N;
-  std::vector<Vertex<D>> vertices;
+  std::vector<Vertex<D, DistinguishableState>> vertices;
   std::vector<double> interaction;
 
   unsigned vectorToIndex(const Vector<D>& x) const {
@@ -218,8 +51,8 @@ public:
       Δ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>(vertices[j], Δx));
-        vertices[j].adjacentEdges.push_back(HalfEdge<D>(vertices[i], -Δx));
+        vertices[i].adjacentEdges.push_back(HalfEdge<D, DistinguishableState>(vertices[j], Δx));
+        vertices[j].adjacentEdges.push_back(HalfEdge<D, DistinguishableState>(vertices[i], -Δx));
       }
     }
 
@@ -232,8 +65,8 @@ public:
     }
   }
 
-  double pairEnergy(const State& s1, const State& s2) const {
-    if (s1.isEmpty() || s2.isEmpty()) {
+  double pairEnergy(const DistinguishableState& s1, const DistinguishableState& s2) const {
+    if (s1.empty() || s2.empty()) {
       return 0;
     } else {
       if (s1.type < s2.type) {
@@ -244,10 +77,10 @@ public:
     }
   }
 
-  double siteEnergy(const Vertex<D>& v) const {
+  double siteEnergy(const Vertex<D, DistinguishableState>& v) const {
     double E = 0;
 
-    for (const HalfEdge<D>& e : v.adjacentEdges) {
+    for (const HalfEdge<D, DistinguishableState>& e : v.adjacentEdges) {
       E += pairEnergy(v.state, e.neighbor.state);
     }
 
@@ -257,14 +90,15 @@ public:
   double energy() const {
     double E = 0;
 
-    for (const Vertex<D>& v : vertices) {
+    for (const Vertex<D, DistinguishableState>& v : vertices) {
       E += siteEnergy(v);
     }
 
     return E / 2;
   }
 
-  bool trySwap(Vertex<D>& v1, Vertex<D>& v2, double T, Rng& r) {
+  bool trySwap(Vertex<D, DistinguishableState>& v1, Vertex<D, DistinguishableState>& v2, double T,
+               Rng& r) {
     double E0 = siteEnergy(v1) + siteEnergy(v2);
     std::swap(v1.state, v2.state);
     double E1 = siteEnergy(v1) + siteEnergy(v2);
@@ -282,15 +116,15 @@ public:
   }
 
   bool tryRandomMove(double T, Rng& r) {
-    Vertex<D>& v1 = r.pick(vertices);
+    Vertex<D, DistinguishableState>& v1 = r.pick(vertices);
 
-    if (v1.isEmpty()) {
+    if (v1.empty()) {
       return false;
     }
 
-    Vertex<D>& v2 = (r.pick(v1.adjacentEdges)).neighbor;
+    Vertex<D, DistinguishableState>& v2 = (r.pick(v1.adjacentEdges)).neighbor;
 
-    if (!v2.isEmpty()) {
+    if (!v2.empty()) {
       return false;
     }
 
@@ -298,8 +132,8 @@ public:
   }
 
   bool tryRandomSwap(double T, Rng& r) {
-    Vertex<D>& v1 = r.pick(vertices);
-    Vertex<D>& v2 = r.pick(vertices);
+    Vertex<D, DistinguishableState>& v1 = r.pick(vertices);
+    Vertex<D, DistinguishableState>& v2 = r.pick(vertices);
 
     if (v1.state.type != v2.state.type) {
       return trySwap(v1, v2, T, r);
@@ -322,10 +156,11 @@ public:
     }
   }
 
-  unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0, double T, Rng& r, bool dry = false) {
-    std::queue<std::array<std::reference_wrapper<Vertex<D>>, 2>> q;
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D, DistinguishableState>& v0, double T,
+                       Rng& r, bool dry = false) {
+    std::queue<std::array<std::reference_wrapper<Vertex<D, DistinguishableState>>, 2>> q;
     Vector<D> x0New = R.apply(v0.position);
-    Vertex<D>& v0New = vertices[vectorToIndex(x0New)];
+    Vertex<D, DistinguishableState>& v0New = vertices[vectorToIndex(x0New)];
     q.push({v0, v0New});
 
     unsigned n = 0;
@@ -334,24 +169,24 @@ public:
       auto [vR, vNewR] = q.front();
       q.pop();
 
-      Vertex<D>& v = vR;
-      Vertex<D>& vNew = vNewR;
+      Vertex<D, DistinguishableState>& v = vR;
+      Vertex<D, DistinguishableState>& vNew = vNewR;
 
       if (!v.marked && !vNew.marked) {
         v.marked = true;
         vNew.marked = true;
 
-        for (HalfEdge<D>& e : v.adjacentEdges) {
-          Vertex<D>& vn = e.neighbor;
+        for (HalfEdge<D, DistinguishableState>& e : v.adjacentEdges) {
+          Vertex<D, DistinguishableState>& vn = e.neighbor;
 
           Vector<D> xnNew = R.apply(vn.position);
-          Vertex<D>& vnNew = vertices[vectorToIndex(xnNew)];
+          Vertex<D, DistinguishableState>& vnNew = vertices[vectorToIndex(xnNew)];
 
           if (!vn.marked && !vnNew.marked) {
             double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state);
             double E1 = pairEnergy(vNew.state, vn.state) + pairEnergy(v.state, vnNew.state);
             double ΔE = E1 - E0;
-            
+
             if (exp(-ΔE / T) < r.uniform(0.0, 1.0)) {
               q.push({vn, vnNew});
             }
@@ -372,21 +207,21 @@ public:
 
   unsigned wolff(const Transformation<D>& R, double T, Rng& r) {
     unsigned n = flipCluster(R, r.pick(vertices), T, r);
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, DistinguishableState>& v : vertices) {
       v.marked = false;
     }
     return n;
   }
 
   void swendsenWang(const Transformation<D>& R, double T, Rng& r) {
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, DistinguishableState>& v : vertices) {
       if (!v.marked) {
         bool dry = 0.5 < r.uniform(0.0, 1.0);
         flipCluster(R, v, T, r, dry);
       }
     }
 
-    for (Vertex<D>& v : vertices) {
+    for (Vertex<D, DistinguishableState>& v : vertices) {
       v.marked = false;
     }
   }
@@ -395,15 +230,16 @@ public:
     int o = 0;
 
     for (unsigned i = 0; i < vertices.size(); i++) {
-//      State<D> s2 = orientation.apply(s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].state);
-  //    o += vertices[i].state.dot(s2);
+      //      DistinguishableState<D> s2 =
+      //      orientation.apply(s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].state);
+      //    o += vertices[i].state.dot(s2);
     }
 
     return o;
   }
 
   void setInitialPosition() {
-    for (Vertex<D>& v: vertices) {
+    for (Vertex<D, DistinguishableState>& v : vertices) {
       v.initialPosition = v.position;
     }
   }
@@ -413,8 +249,8 @@ public:
 
     Vector<D> k1 = {M_PI, 0};
     Vector<D> k2 = {0, M_PI};
-    for (const Vertex<D>& v : vertices) {
-      if (!v.isEmpty()) {
+    for (const Vertex<D, DistinguishableState>& v : vertices) {
+      if (!v.empty()) {
         F += cos(k1.dot(v.position - v.initialPosition));
         F += cos(k2.dot(v.position - v.initialPosition));
       }
@@ -447,7 +283,7 @@ int main() {
   std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
 
   for (unsigned j = 0; j < 10 * s.vertices.size(); j++) {
-    //unsigned nC = s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
+    // unsigned nC = s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
     s.tryRandomSwap(1000, r);
   }
 
@@ -463,17 +299,17 @@ int main() {
     /*
     for (unsigned i = 0; i < 1e5; i++) {
       Matrix<D> m = r.pick(ms);
-      Vertex<D>& v = r.pick(s.vertices);
+      Vertex<D, DistinguishableState>& v = r.pick(s.vertices);
       s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
     }
     */
     /*
-    */
+     */
     s.setInitialPosition();
     std::cout << T << " ";
     for (unsigned i = 0; i < 1e4; i++) {
       Matrix<D> m = r.pick(ms);
-      Vertex<D>& v = r.pick(s.vertices);
+      Vertex<D, DistinguishableState>& v = r.pick(s.vertices);
       s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
       std::cout << s.selfIntermediateScattering() << " ";
     }
@@ -487,9 +323,9 @@ int main() {
     */
     std::cout << std::endl;
     std::cerr << T << " " << s.energy() / N << std::endl;
-//    s.sweepLocal(r);
-//    s.sweepSwap(r);
-//    s.swendsenWang(Transformation<D>(L, ms, r), r);
+    //    s.sweepLocal(r);
+    //    s.sweepSwap(r);
+    //    s.swendsenWang(Transformation<D>(L, ms, r), r);
   }
 
   return 0;
diff --git a/glass.hpp b/glass.hpp
new file mode 100644
index 0000000..fbff7d6
--- /dev/null
+++ b/glass.hpp
@@ -0,0 +1,167 @@
+#include <list>
+
+#include <eigen3/Eigen/Dense>
+#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
+
+#include "pcg-cpp/include/pcg_random.hpp"
+#include "randutils/randutils.hpp"
+
+template <int D> using Vector = Eigen::Matrix<int, D, 1>;
+template <int D> using Matrix = Eigen::Matrix<int, D, D>;
+
+using Rng = randutils::random_generator<pcg32>;
+
+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 <unsigned 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 <unsigned 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;
+  }
+
+  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;
+  }
+
+  Transformation<D> inverse() const {
+    return Transformation<D>(L, m.transpose(), -m.transpose() * v);
+  }
+};
+
+template <typename T> concept State = requires(T v) {
+  { v.empty() } -> std::same_as<bool>;
+  {v.remove()};
+};
+
+template <unsigned D, State S> class HalfEdge;
+
+template <unsigned D, State S> class Vertex {
+public:
+  Vector<D> position;
+  Vector<D> initialPosition;
+  S state;
+  std::vector<HalfEdge<D, S>> adjacentEdges;
+  bool marked;
+
+  bool empty() const { return state.empty(); }
+};
+
+template <unsigned D, State S> class HalfEdge {
+public:
+  Vertex<D, S>& neighbor;
+  Vector<D> Δx;
+
+  HalfEdge(Vertex<D, S>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
+};
+