diff options
| -rw-r--r-- | glass.cpp | 138 |
1 files changed, 88 insertions, 50 deletions
@@ -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; } |