From 3b470b46f600ee56babb71a337f2acfe7b9ae4b3 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 24 Jun 2021 14:08:57 +0200
Subject: More refactoring.

---
 glass.hpp | 230 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 229 insertions(+), 1 deletion(-)

(limited to 'glass.hpp')

diff --git a/glass.hpp b/glass.hpp
index fbff7d6..64fe8b7 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -1,4 +1,5 @@
 #include <list>
+#include <queue>
 
 #include <eigen3/Eigen/Dense>
 #include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
@@ -11,6 +12,14 @@ template <int D> using Matrix = Eigen::Matrix<int, D, D>;
 
 using Rng = randutils::random_generator<pcg32>;
 
+template <unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
+  Vector<D> x;
+  for (unsigned i = 0; i < D; i++) {
+    x[i] = r.uniform((unsigned)0, L - 1);
+  }
+  return x;
+}
+
 int iPow(int x, unsigned p) {
   if (p == 0)
     return 1;
@@ -139,9 +148,10 @@ public:
   }
 };
 
-template <typename T> concept State = requires(T v) {
+template <typename T> concept State = requires(T v, const T& v2) {
   { v.empty() } -> std::same_as<bool>;
   {v.remove()};
+  { v.operator==(v2) } -> std::same_as<bool>;
 };
 
 template <unsigned D, State S> class HalfEdge;
@@ -165,3 +175,221 @@ public:
   HalfEdge(Vertex<D, S>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
 };
 
+template <unsigned D, State S> class System {
+public:
+  const unsigned L;
+  unsigned N;
+  std::vector<Vertex<D, S>> vertices;
+
+  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 N) : L(L), N(N), vertices(iPow(L, D)) {
+    for (unsigned i = 0; i < iPow(L, D); i++) {
+      vertices[i].position = indexToVector(i);
+      vertices[i].initialPosition = vertices[i].position;
+      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, S>(vertices[j], Δx));
+        vertices[j].adjacentEdges.push_back(HalfEdge<D, S>(vertices[i], -Δx));
+      }
+    }
+  }
+
+  double density() const { return N / pow(L, D); }
+
+  void setInitialPosition() {
+    for (Vertex<D, S>& v : vertices) {
+      v.initialPosition = v.position;
+    }
+  }
+
+  double selfIntermediateScattering() const {
+    double F = 0;
+
+    Vector<D> k1 = {M_PI, 0};
+    Vector<D> k2 = {0, M_PI};
+    for (const Vertex<D, S>& v : vertices) {
+      if (!v.empty()) {
+        F += cos(k1.dot(v.position - v.initialPosition));
+        F += cos(k2.dot(v.position - v.initialPosition));
+      }
+    }
+
+    return F / (2 * this->N);
+  }
+};
+
+template <unsigned D, State S> class FiniteEnergySystem : public System<D, S> {
+public:
+  using System<D, S>::System;
+
+  virtual double pairEnergy(const S&, const S&) const { return 0; }
+
+  double siteEnergy(const Vertex<D, S>& v) const {
+    double E = 0;
+
+    for (const HalfEdge<D, S>& e : v.adjacentEdges) {
+      E += pairEnergy(v.state, e.neighbor.state);
+    }
+
+    return E;
+  }
+
+  double energy() const {
+    double E = 0;
+
+    for (const Vertex<D, S>& v : this->vertices) {
+      E += siteEnergy(v);
+    }
+
+    return E / 2;
+  }
+
+  bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& v2, double T, Rng& r) {
+    double E0 = siteEnergy(v1) + siteEnergy(v2);
+    std::swap(v1.state, v2.state);
+    double E1 = siteEnergy(v1) + siteEnergy(v2);
+    double ΔE = E1 - E0;
+
+    if (exp(-ΔE / T) > r.uniform(0.0, 1.0)) {
+      /* Accept the swap. */
+      std::swap(v1.initialPosition, v2.initialPosition);
+      return true;
+    } else {
+      /* Revert the swap. */
+      std::swap(v1.state, v2.state);
+      return false;
+    }
+  }
+
+  bool tryRandomMove(double T, Rng& r) {
+    Vertex<D, S>& v1 = r.pick(this->vertices);
+
+    if (v1.empty()) {
+      return false;
+    }
+
+    Vertex<D, S>& v2 = (r.pick(v1.adjacentEdges)).neighbor;
+
+    if (!v2.empty()) {
+      return false;
+    }
+
+    return trySwap(v1, v2, T, r);
+  }
+
+  bool tryRandomSwap(double T, Rng& r) {
+    Vertex<D, S>& v1 = r.pick(this->vertices);
+    Vertex<D, S>& v2 = r.pick(this->vertices);
+
+    if (v1.state != v2.state) {
+      return trySwap(v1, v2, T, r);
+    } else {
+      return false;
+    }
+  }
+
+  void sweepLocal(Rng& r) {
+    for (unsigned i = 0; i < iPow(this->L, D); i++) {
+      tryRandomMove(r);
+    }
+  }
+
+  void sweepSwap(Rng& r) {
+    for (unsigned i = 0; i < iPow(this->L, D); i++) {
+      tryRandomSwap(r);
+    }
+  }
+
+  unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, double T, Rng& r,
+                       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)];
+    q.push({v0, v0New});
+
+    unsigned n = 0;
+
+    while (!q.empty()) {
+      auto [vR, vNewR] = q.front();
+      q.pop();
+
+      Vertex<D, S>& v = vR;
+      Vertex<D, S>& vNew = vNewR;
+
+      if (!v.marked && !vNew.marked) {
+        v.marked = true;
+        vNew.marked = true;
+
+        for (HalfEdge<D, S>& e : v.adjacentEdges) {
+          Vertex<D, S>& vn = e.neighbor;
+
+          Vector<D> xnNew = R.apply(vn.position);
+          Vertex<D, S>& vnNew = this->vertices[this->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});
+            }
+          }
+        }
+
+        if (!dry) {
+          std::swap(v.state, vNew.state);
+          std::swap(v.initialPosition, vNew.initialPosition);
+        }
+
+        n += 1;
+      }
+    }
+
+    return n;
+  }
+
+  unsigned wolff(const Transformation<D>& R, double T, Rng& r) {
+    unsigned n = flipCluster(R, r.pick(this->vertices), T, r);
+    for (Vertex<D, S>& v : this->vertices) {
+      v.marked = false;
+    }
+    return n;
+  }
+
+  void swendsenWang(const Transformation<D>& R, double T, Rng& r) {
+    for (Vertex<D, S>& v : this->vertices) {
+      if (!v.marked) {
+        bool dry = 0.5 < r.uniform(0.0, 1.0);
+        flipCluster(R, v, T, r, dry);
+      }
+    }
+
+    for (Vertex<D, S>& v : this->vertices) {
+      v.marked = false;
+    }
+  }
+};
+
-- 
cgit v1.2.3-54-g00ecf