From 7548abd44c37d4495dd498193ea15b0df757b904 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 9 Jun 2021 16:32:27 +0200 Subject: Lots of work. --- biroli-mezard.cpp | 163 ++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 121 insertions(+), 42 deletions(-) (limited to 'biroli-mezard.cpp') diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp index 7e200b9..e2fb71a 100644 --- a/biroli-mezard.cpp +++ b/biroli-mezard.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include @@ -152,6 +153,7 @@ public: unsigned occupiedNeighbors; unsigned maximumNeighbors; bool marked; + bool visited; Vertex() { occupiedNeighbors = 0; @@ -181,6 +183,7 @@ template class System { public: const unsigned L; std::vector N; + std::unordered_set*> occupants; std::vector> vertices; Transformation orientation; @@ -322,6 +325,9 @@ public: } N[t]++; N[0]--; + if (t > 0) { + occupants.insert(&v); + } return true; } else { return false; @@ -332,6 +338,9 @@ public: if (v.empty()) { return false; } else { + if (v.maximumNeighbors > 0) { + occupants.erase(&v); + } N[v.maximumNeighbors]--; N[0]++; v.maximumNeighbors = Empty; @@ -418,30 +427,32 @@ public: return exchange(v1, v2); } - void sweepGrandCanonical(double z, Rng& r) { - for (unsigned i = 0; i < size(); i++) { - if (0.5 < r.uniform(0.0, 1.0)) { - double pIns = size() * z / (occupancy() + 1); - - if (pIns > r.uniform(0.0, 1.0)) { - while (true) { - Vertex& v = r.pick(vertices); - if (v.empty()) { - insert(v, r.pick({1, 2, 2, 2, 2, 2, 3, 3, 3, 3})); - break; - } - } - } - } else { + const std::array ratios = {0.1, 0.5, 0.4}; - double pDel = density() / z; + void stepGrandCanonical(double z, Rng& r) { + if (1.0 / 11.0 < r.uniform(0.0, 1.0)) { + unsigned t = r.pick({1, 2, 2, 2, 2, 2, 3, 3, 3, 3}); + double pIns = size() * z / (occupancy() + 1); - if (pDel > r.uniform(0.0, 1.0)) { - remove(r.pick(vertices)); - } + if (pIns > r.uniform(0.0, 1.0)) { + Vertex& v = r.pick(vertices); + insert(v, t); } + } else { + double pDel = density() / z; - randomLocalExchange(r); + if (pDel > r.uniform(0.0, 1.0)) { + Vertex* v = r.pick(occupants); + remove(*v); + } + } + } + + void sweepGrandCanonical(double z, Rng& r) { + for (unsigned i = 0; i < size(); i++) { + stepGrandCanonical(z, r); + + randomExchange(r); } } @@ -457,6 +468,71 @@ public: } } + Vertex& apply(const Transformation& R, const Vertex& v) { + return vertices[vectorToIndex(R.apply(v.position))]; + } + + bool eventChain(const Transformation& R, Vertex& v0) { + unsigned n = 0; + + if (!v0.empty()) { + Vertex& v0New = apply(R, v0);; + unsigned t0 = v0.maximumNeighbors; + + std::queue>, unsigned>> q; + + q.push({v0New, t0}); + remove(v0); + + while (!q.empty()) { + auto [vr, t] = q.front(); + Vertex& v = vr; + q.pop(); + + if (!v.empty()) { + q.push({apply(R, v), v.maximumNeighbors}); + remove(v); + } + + insert(v, t, true); + v.marked = true; + + for (Vertex& vn : overlaps({v})) { + if (!vn.marked) { + q.push({apply(R, vn), vn.maximumNeighbors}); + remove(vn); + vn.marked = true; + } else { + /* If a neighbor has already been moved but is still + * frustrated, it may be due to one of its neighbors! + */ + if (&vn != &v) { + for (Vertex& vnn : overlaps(vn)) { + if (!vnn.marked) { + q.push({apply(R, vnn), vnn.maximumNeighbors}); + remove(vnn); + vnn.marked = true; + } + } + } + } + } + } + } + + for (Vertex& v : vertices) { + v.marked = false; + } + + return true; + } + + bool randomEventChain(Rng& r) { + Transformation R(L); + R.v[r.uniform((unsigned)0, (unsigned)D-1)] = r.pick({-1, 1}); + return eventChain(R, r.pick(vertices)); + } + unsigned flipCluster(const Transformation& R, Vertex& v0) { unsigned n = 0; @@ -546,30 +622,35 @@ int main() { unsigned L = 15; unsigned Nmin = 2e2; unsigned Nmax = 2e5; - double Tmin = 0.04; - double Tmax = 0.2; - double δT = 0.02; + double μmin = -3; + double μmax = 10; + double dμ = 0.05; System s(L, 3); Rng r; - double z = exp(1 / 0.15); + for (double μ = μmin; μ <= μmax; μ += dμ) { + double z = exp(μ); + for (unsigned i = 0; i < 1e4; i++) { + for (unsigned j = 0; j < s.size(); j++) { + s.stepGrandCanonical(z, r); +// s.randomEventChain(r); + s.randomExchange(r); + } - if (!s.compatible()) { - std::cerr << "Storted incompatible!" << std::endl; - return 1; - } + if (!s.compatible()) { + std::cerr << "Not compatible!" << std::endl; + return 1; + } - while (s.density() < 0.56) { - s.sweepGrandCanonical(z, r); - } + } - if (!s.compatible()) { - std::cerr << "Not compatible!" << std::endl; - return 1; + std::cout << μ << " " << s.density() << std::endl; } + + /* std::cerr << "Found state with appropriate density." << std::endl; System s0 = s; @@ -586,22 +667,19 @@ int main() { auto start = std::chrono::high_resolution_clock::now(); while (nC / s.size() < 1e5) { if (n < 20 * log(i + 1)) { + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = duration_cast(stop - start); n++; - std::cout << nC / s.size() << " " << s.overlap(s0) << std::endl; + std::cout << duration.count() << " " << s.overlap(s0) << std::endl; } - //unsigned nn = s.flipCluster(Transformation(L, ms, r), r.pick(s.vertices)); - //nC += nn; +// unsigned nn = s.flipCluster(Transformation(L, ms, r), r.pick(s.vertices)); +// nC += nn; // clusterDist[nn]++; //s.sweepLocal(r); nC += s.size(); s.sweepSwap(r); i++; } - auto stop = std::chrono::high_resolution_clock::now(); - - auto duration = duration_cast(stop - start); - - std::cerr << duration.count() << std::endl; if (!s.compatible()) { std::cerr << "Not compatible!" << std::endl; @@ -618,6 +696,7 @@ int main() { file << i << " "; } file.close(); + */ return 0; } -- cgit v1.2.3-54-g00ecf