From c6076430a445ba07866eb7fbdb083f016be57894 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 21 Mar 2021 22:51:35 +0100 Subject: Updates. --- glass.cpp | 138 +++++++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 88 insertions(+), 50 deletions(-) (limited to 'glass.cpp') diff --git a/glass.cpp b/glass.cpp index ba983d2..f2f44c1 100644 --- a/glass.cpp +++ b/glass.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -117,6 +118,14 @@ public: bool isEmpty() const { return this->squaredNorm() == 0; } void remove() { this->setZero(); } + + State flip() const { + State s; + for (unsigned i = 0; i < D; i++) { + s(i) = -this->operator()(i); + } + return s; + } }; template class Transformation { @@ -125,21 +134,41 @@ template class Transformation { Matrix m; Vector v; + Transformation(unsigned L) : L(L) { + m.setIdentity(); + v.setZero(); + } + Transformation(unsigned L, const Matrix& m, const Vector& v) : L(L), m(m), v(v) {} Transformation(unsigned L, const std::vector>& 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 apply(const Vector& x) const { - return mod(v + m * (x - v), L); + return mod(v + m * x, L); + } + + Transformation apply(const Transformation& t) const { + Transformation tNew(L); + + tNew.m = m * t.m; + tNew.v = apply(t.v); + + return tNew; } State apply(const State& x) const { return State(m * x); } + + Transformation inverse() const { + return Transformation(L, m.transpose(), -m.transpose() * v); + } }; template class HalfEdge; @@ -167,6 +196,7 @@ public: const unsigned L; unsigned N; std::vector> vertices; + Transformation orientation; unsigned vectorToIndex(const Vector& 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& v = r.pick(vertices); State oldState = v.state; + if (!tryDeletion(v)) { return false; } - if (0.2 > r.uniform(0.0, 1.0)) { - if (insert(v, State(r))) { - return true; + if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) { + for (HalfEdge& 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(r))) { + State newState(r); + while (newState.dot(oldState) == 1) { + newState = State(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& R, Rng& r) { for (Vertex& 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& v : vertices) { v.marked = false; } } + + int overlap(const System& s) const { + int o = 0; + + for (unsigned i = 0; i < vertices.size(); i++) { + State 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 Vector 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> ms = generateTorusMatrices(); - - 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(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 s0 = s; - for (unsigned n = 0; n < N; n++) { - s.sweepGrandCanonical(z, r); - } + std::vector> ms = generateTorusMatrices(); - 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 m = r.pick(ms); + Vertex& v = r.pick(s.vertices); + unsigned nC = s.flipCluster(Transformation(L, m, v.position - m * v.position), v, r); + std::cerr << nC << std::endl; + for (Vertex& v : s.vertices) { + v.marked = false; } +// s.sweepLocal(r); +// s.sweepSwap(r); +// s.swendsenWang(Transformation(L, ms, r), r); } - */ return 0; } -- cgit v1.2.3-70-g09d2