Updates.

1 files changed, 88 insertions, 50 deletions

diff --git a/glass.cpp b/glass.cpp
index ba983d2..f2f44c1 100644
--- a/glass.cpp
+++ b/glass.cpp
@@ -1,4 +1,5 @@
 #include <eigen3/Eigen/Dense>
+#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
 #include <iostream>
 #include <list>
 #include <vector>
@@ -117,6 +118,14 @@ public:
   bool isEmpty() const { return this->squaredNorm() == 0; }
 
   void remove() { this->setZero(); }
+
+  State<D> flip() const {
+    State<D> s;
+    for (unsigned i = 0; i < D; i++) {
+      s(i) = -this->operator()(i);
+    }
+    return s;
+  }
 };
 
 template <unsigned D> class Transformation {
@@ -125,21 +134,41 @@ template <unsigned D> class Transformation {
     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 - v), L);
+      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;
     }
 
     State<D> apply(const State<D>& x) const {
       return State<D>(m * x);
     }
+
+    Transformation<D> inverse() const {
+      return Transformation<D>(L, m.transpose(), -m.transpose() * v);
+    }
 };
 
 template <unsigned D> class HalfEdge;
@@ -167,6 +196,7 @@ public:
   const unsigned L;
   unsigned N;
   std::vector<Vertex<D>> vertices;
+  Transformation<D> orientation;
 
   unsigned vectorToIndex(const Vector<D>& x) const {
     unsigned i = 0;
@@ -184,7 +214,7 @@ public:
     return x;
   }
 
-  System(unsigned L) : L(L), N(0), vertices(iPow(L, D)) {
+  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);
@@ -248,16 +278,26 @@ public:
   bool tryRandomMove(Rng& r) {
     Vertex<D>& v = r.pick(vertices);
     State<D> oldState = v.state;
+
     if (!tryDeletion(v)) {
       return false;
     }
 
-    if (0.2 > r.uniform(0.0, 1.0)) {
-      if (insert(v, State<D>(r))) {
-        return true;
+    if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) {
+      for (HalfEdge<D>& e : v.adjacentEdges) {
+        if (1 == e.Δx.dot(oldState)) {
+          if (insert(e.neighbor, oldState.flip())) {
+            return true;
+          }
+          break;
+        }
       }
     } else {
-      if (insert(r.pick(v.adjacentEdges).neighbor, State<D>(r))) {
+      State<D> newState(r);
+      while (newState.dot(oldState) == 1) {
+        newState = State<D>(r);
+      }
+      if (insert(v, newState)) {
         return true;
       }
     }
@@ -345,17 +385,16 @@ public:
       }
 
       tryRandomMove(r);
-      tryRandomSwap(r);
     }
   }
 
-  void sweepLocal(double z, Rng& r) {
+  void sweepLocal(Rng& r) {
     for (unsigned i = 0; i < iPow(L, D); i++) {
       tryRandomMove(r);
     }
   }
 
-  void sweepSwap(double z, Rng& r) {
+  void sweepSwap(Rng& r) {
     for (unsigned i = 0; i < iPow(L, D); i++) {
       tryRandomSwap(r);
     }
@@ -406,17 +445,29 @@ public:
   void swendsenWang(const Transformation<D>& R, Rng& r) {
     for (Vertex<D>& v : vertices) {
       if (!v.marked) {
-        unsigned n = flipCluster(R, v, r, 0.5 < r.uniform(0.0, 1.0));
-        std::cerr << n << " ";
+        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);
+        }
       }
     }
 
-    std::cerr << std::endl;
-
     for (Vertex<D>& v : vertices) {
       v.marked = false;
     }
   }
+
+  int overlap(const System<D>& s) const {
+    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);
+    }
+
+    return o;
+  }
 };
 
 void print(const System<2>& s) {
@@ -451,7 +502,7 @@ template <unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
 
 int main() {
   const unsigned D = 3;
-  unsigned L = 28;
+  unsigned L = 15;
   unsigned Nmin = 2e2;
   unsigned Nmax = 2e5;
   double Tmin = 0.04;
@@ -462,53 +513,40 @@ int main() {
 
   Rng r;
 
-  double z = exp(1 / 0.15);
+  double z = exp(1 / 0.08);
 
-  s.setGroundState();
-  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
-
-  for (unsigned n = 0; n < 1e3; n++) {
+  while (s.density() < 0.818) {
     s.sweepGrandCanonical(z, r);
   }
 
-  if (s.compatible()) {
-    std::cerr << "Still compatible!" << std::endl;
-  }
-
-
-  for (unsigned n = 0; n < 10; n++) {
-    s.swendsenWang(Transformation<D>(L, ms, r), r);
-  }
-
-  if (s.compatible()) {
-    std::cerr << "Still compatible!" << std::endl;
-  }
+  if (!s.compatible()) {
+    std::cerr << "Net compatible!" << std::endl;
+    return 1;
+  } 
 
-  /*
+  std::cerr << "Found state with appropriate density." << std::endl;
 
-  for (unsigned N = Nmin; N <= Nmax; N *= 10) {
-    s.setGroundState();
-    for (double T = Tmin; T < Tmax; T *= 1 + δT) {
-      double z = exp(1 / T);
+  System<D> s0 = s;
 
-      for (unsigned n = 0; n < N; n++) {
-        s.sweepGrandCanonical(z, r);
-      }
+  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
 
-      std::cout << L << " " << T << " " << N << " " << δT << " " << 1 / s.density() << std::endl;
+  unsigned n = 1;
+  for (unsigned i = 0; i < 1e5; i++) {
+    if (n < 20 * log(i + 1)) {
+      n++;
+      std::cout << i << " " << s.overlap(s0) << std::endl;
     }
-
-    for (double T = Tmax; T > Tmin; T /= 1 + δT) {
-      double z = exp(1 / T);
-
-      for (unsigned n = 0; n < N; n++) {
-        s.sweepGrandCanonical(z, r);
-      }
-
-      std::cout << L << " " << T << " " << N << " " << δT << " " << 1 / s.density() << std::endl;
+    Matrix<D> m = r.pick(ms);
+    Vertex<D>& v = r.pick(s.vertices);
+    unsigned nC = s.flipCluster(Transformation<D>(L, m, v.position - m * v.position), v, r);
+    std::cerr << nC << std::endl;
+    for (Vertex<D>& v : s.vertices) {
+      v.marked = false;
     }
+//    s.sweepLocal(r);
+//    s.sweepSwap(r);
+//    s.swendsenWang(Transformation<D>(L, ms, r), r);
   }
-  */
 
   return 0;
 }