From 7548abd44c37d4495dd498193ea15b0df757b904 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 9 Jun 2021 16:32:27 +0200
Subject: Lots of work.
---
biroli-mezard.cpp | 163 +++++++++----
ciamarra.cpp | 489 ++++++++++++++++++++++++++++++++++++++-
glass.hpp | 677 ------------------------------------------------------
3 files changed, 601 insertions(+), 728 deletions(-)
delete mode 100644 glass.hpp
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 <queue>
#include <vector>
#include <chrono>
+#include <unordered_set>
#include <eigen3/Eigen/Dense>
@@ -152,6 +153,7 @@ public:
unsigned occupiedNeighbors;
unsigned maximumNeighbors;
bool marked;
+ bool visited;
Vertex() {
occupiedNeighbors = 0;
@@ -181,6 +183,7 @@ template <int D> class System {
public:
const unsigned L;
std::vector<unsigned> N;
+ std::unordered_set<Vertex<D>*> occupants;
std::vector<Vertex<D>> vertices;
Transformation<D> 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<D>& 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<double, 3> 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<D>& v = r.pick(vertices);
+ insert(v, t);
}
+ } else {
+ double pDel = density() / z;
- randomLocalExchange(r);
+ if (pDel > r.uniform(0.0, 1.0)) {
+ Vertex<D>* 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<D>& apply(const Transformation<D>& R, const Vertex<D>& v) {
+ return vertices[vectorToIndex(R.apply(v.position))];
+ }
+
+ bool eventChain(const Transformation<D>& R, Vertex<D>& v0) {
+ unsigned n = 0;
+
+ if (!v0.empty()) {
+ Vertex<D>& v0New = apply(R, v0);;
+ unsigned t0 = v0.maximumNeighbors;
+
+ std::queue<std::pair<std::reference_wrapper<Vertex<D>>, unsigned>> q;
+
+ q.push({v0New, t0});
+ remove(v0);
+
+ while (!q.empty()) {
+ auto [vr, t] = q.front();
+ Vertex<D>& 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<D>& 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<D>& vnn : overlaps(vn)) {
+ if (!vnn.marked) {
+ q.push({apply(R, vnn), vnn.maximumNeighbors});
+ remove(vnn);
+ vnn.marked = true;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ for (Vertex<D>& v : vertices) {
+ v.marked = false;
+ }
+
+ return true;
+ }
+
+ bool randomEventChain(Rng& r) {
+ Transformation<D> 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<D>& R, Vertex<D>& 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<D> 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<D> 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<std::chrono::microseconds>(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<D>(L, ms, r), r.pick(s.vertices));
- //nC += nn;
+// unsigned nn = s.flipCluster(Transformation<D>(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<std::chrono::microseconds>(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;
}
diff --git a/ciamarra.cpp b/ciamarra.cpp
index 6b932bf..f2f44c1 100644
--- a/ciamarra.cpp
+++ b/ciamarra.cpp
@@ -1,10 +1,477 @@
+#include <eigen3/Eigen/Dense>
+#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
#include <iostream>
+#include <list>
+#include <vector>
+#include <queue>
-#include "glass.hpp"
+#include "pcg-cpp/include/pcg_random.hpp"
+#include "randutils/randutils.hpp"
+using Rng = randutils::random_generator<pcg32>;
-void print(const CiamarraSystem<2>& s) {
- for (const Vertex<2, CiamarraState<2>>& v : s.vertices) {
+template <int D> using Vector = Eigen::Matrix<int, D, 1>;
+template <int D> using Matrix = Eigen::Matrix<int, D, D>;
+
+int iPow(int x, unsigned p) {
+ if (p == 0)
+ return 1;
+ if (p == 1)
+ return x;
+
+ int tmp = iPow(x, p / 2);
+ if (p % 2 == 0)
+ return tmp * tmp;
+ else
+ return x * tmp * tmp;
+}
+
+unsigned mod(signed a, unsigned b) {
+ return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b;
+}
+
+template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) {
+ Vector<D> u;
+ for (unsigned i = 0; i < D; i++) {
+ u(i) = mod(v(i), b);
+ }
+ return u;
+}
+
+template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+ if (level > 0) {
+ unsigned new_level = level - 1;
+ unsigned old_length = sequences.size();
+ for (std::array<double, D>& sequence : sequences) {
+ std::array<double, D> new_sequence = sequence;
+ new_sequence[new_level] = -1;
+ sequences.push_front(new_sequence);
+ }
+ one_sequences<D>(sequences, new_level);
+ }
+}
+
+template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
+ std::vector<Matrix<D>> mats;
+
+ std::array<double, D> ini_sequence;
+ ini_sequence.fill(1);
+ std::list<std::array<double, D>> sequences;
+ sequences.push_back(ini_sequence);
+
+ one_sequences<D>(sequences, D);
+
+ sequences.pop_back(); // don't want the identity matrix!
+
+ for (std::array<double, D> sequence : sequences) {
+ Matrix<D> m;
+ for (unsigned i = 0; i < D; i++) {
+ for (unsigned j = 0; j < D; j++) {
+ if (i == j) {
+ m(i, j) = sequence[i];
+ } else {
+ m(i, j) = 0;
+ }
+ }
+ }
+
+ mats.push_back(m);
+ }
+
+ for (unsigned i = 0; i < D; i++) {
+ for (unsigned j = 0; j < D; j++) {
+ if (i != j) {
+ Matrix<D> m;
+ for (unsigned k = 0; k < D; k++) {
+ for (unsigned l = 0; l < D; l++) {
+ if ((k == i && l == j) || (k == j && l == i)) {
+ if (i < j) {
+ m(k, l) = 1;
+ } else {
+ m(k, l) = -1;
+ }
+ } else if (k == l && (k != i && k != j)) {
+ m(k, l) = 1;
+ } else {
+ m(k, l) = 0;
+ }
+ }
+ }
+ mats.push_back(m);
+ }
+ }
+ }
+
+ return mats;
+}
+
+template <unsigned D> class State : public Vector<D> {
+public:
+ State() : Vector<D>(Vector<D>::Zero()) {}
+
+ State(const Vector<D>& v) : Vector<D>(v) {}
+
+ State(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; }
+
+ State(Rng& r) : State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {}
+
+ 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 {
+ public:
+ unsigned L;
+ 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, 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;
+
+template <unsigned D> class Vertex {
+public:
+ Vector<D> position;
+ State<D> state;
+ std::vector<HalfEdge<D>> adjacentEdges;
+ bool marked;
+
+ bool isEmpty() const { return state.isEmpty(); }
+};
+
+template <unsigned D> class HalfEdge {
+public:
+ Vertex<D>& neighbor;
+ Vector<D> Δx;
+
+ HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
+};
+
+template <unsigned D> class System {
+public:
+ const unsigned L;
+ unsigned N;
+ std::vector<Vertex<D>> vertices;
+ Transformation<D> orientation;
+
+ 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) : 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);
+ 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>(vertices[j], Δx));
+ vertices[j].adjacentEdges.push_back(HalfEdge<D>(vertices[i], -Δx));
+ }
+ }
+ }
+
+ std::list<std::reference_wrapper<Vertex<D>>> overlaps(Vertex<D>& v, const State<D>& s,
+ bool excludeSelf = false) {
+ std::list<std::reference_wrapper<Vertex<D>>> o;
+
+ if (s.isEmpty()) {
+ return o;
+ }
+
+ if (!v.isEmpty() && !excludeSelf) {
+ o.push_back(v);
+ }
+
+ for (const HalfEdge<D>& e : v.adjacentEdges) {
+ if (!e.neighbor.isEmpty()) {
+ if (s.dot(e.Δx) == 1 || e.neighbor.state.dot(e.Δx) == -1) {
+ o.push_back(e.neighbor);
+ }
+ }
+ }
+
+ return o;
+ }
+
+ bool insert(Vertex<D>& v, const State<D>& s) {
+ if (overlaps(v, s).empty()) {
+ v.state = s;
+ N++;
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ bool tryDeletion(Vertex<D>& v) {
+ if (v.isEmpty()) {
+ return false;
+ } else {
+ v.state.remove();
+ N--;
+ return true;
+ }
+ }
+
+ bool tryRandomMove(Rng& r) {
+ Vertex<D>& v = r.pick(vertices);
+ State<D> oldState = v.state;
+
+ if (!tryDeletion(v)) {
+ return false;
+ }
+
+ 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 {
+ State<D> newState(r);
+ while (newState.dot(oldState) == 1) {
+ newState = State<D>(r);
+ }
+ if (insert(v, newState)) {
+ return true;
+ }
+ }
+ v.state = oldState;
+ N++;
+ return false;
+ }
+
+ bool trySwap(Vertex<D>& v1, Vertex<D>& v2) {
+ if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) {
+ std::swap(v1.state, v2.state);
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ bool tryRandomSwap(Rng& r) {
+ Vertex<D>& v1 = r.pick(vertices);
+ Vertex<D>& v2 = r.pick(vertices);
+
+ return trySwap(v1, v2);
+ }
+
+ void setGroundState() {
+ N = 0;
+
+ for (Vertex<D>& v : vertices) {
+ unsigned a = 0;
+ for (unsigned d = 0; d < D; d++) {
+ a += (d + 1) * v.position(d);
+ }
+ a %= 2 * D + 1;
+
+ v.state.setZero() = Vector<D>::Zero();
+
+ if (0 < a && a <= D) {
+ v.state(a - 1) = -1;
+ N++;
+ } else if (D < a) {
+ v.state(2 * D - a) = 1;
+ N++;
+ }
+ }
+ }
+
+ bool compatible() {
+ for (Vertex<D>& v : vertices) {
+ if (overlaps(v, v.state, true).size() > 0) {
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+ double density() const { return N / pow(L, D); }
+
+ unsigned maxOccupation() const {
+ // return (2 * D * iPow(L, D)) / (2 * D + 1);
+ return iPow(L, D);
+ }
+
+ void sweepGrandCanonical(double z, Rng& r) {
+ for (unsigned i = 0; i < iPow(L, D); i++) {
+ if (0.5 < r.uniform(0.0, 1.0)) {
+ double pIns = maxOccupation() * z / (N + 1);
+
+ if (pIns > r.uniform(0.0, 1.0)) {
+ while (true) {
+ Vertex<D>& v = r.pick(vertices);
+ if (v.isEmpty()) {
+ insert(v, State<D>(r));
+ break;
+ }
+ }
+ }
+ } else {
+
+ double pDel = N / (z * maxOccupation());
+
+ if (pDel > r.uniform(0.0, 1.0)) {
+ tryDeletion(r.pick(vertices));
+ }
+ }
+
+ tryRandomMove(r);
+ }
+ }
+
+ void sweepLocal(Rng& r) {
+ for (unsigned i = 0; i < iPow(L, D); i++) {
+ tryRandomMove(r);
+ }
+ }
+
+ void sweepSwap(Rng& r) {
+ for (unsigned i = 0; i < iPow(L, D); i++) {
+ tryRandomSwap(r);
+ }
+ }
+
+ unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0, Rng& r, bool dry = false) {
+ std::queue<std::reference_wrapper<Vertex<D>>> q;
+ q.push(v0);
+
+ unsigned n = 0;
+
+ while (!q.empty()) {
+ Vertex<D>& v = q.front();
+ q.pop();
+
+ if (!v.marked) {
+ Vector<D> xNew = R.apply(v.position);
+ Vertex<D>& vNew = vertices[vectorToIndex(xNew)];
+
+ v.marked = true;
+ vNew.marked = true;
+
+ State<D> s = R.apply(v.state);
+ State<D> sNew = R.apply(vNew.state);
+
+ std::list<std::reference_wrapper<Vertex<D>>> overlaps1 = overlaps(vNew, s, true);
+ std::list<std::reference_wrapper<Vertex<D>>> overlaps2 = overlaps(v, sNew, true);
+ overlaps1.splice(overlaps1.begin(), overlaps2);
+
+ for (Vertex<D>& vn : overlaps1) {
+ if (!vn.marked) {
+ q.push(vn);
+ }
+ }
+
+ if (!dry) {
+ v.state = sNew;
+ vNew.state = s;
+ }
+
+ n += 1;
+ }
+ }
+
+ return n;
+ }
+
+ void swendsenWang(const Transformation<D>& R, Rng& r) {
+ for (Vertex<D>& v : vertices) {
+ if (!v.marked) {
+ 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);
+ }
+ }
+ }
+
+ 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) {
+ for (const Vertex<2>& v : s.vertices) {
if (v.state(0) == 1 && v.state(1) == 0) {
std::cerr << "▶";
} else if (v.state(0) == -1 && v.state(1) == 0) {
@@ -42,7 +509,7 @@ int main() {
double Tmax = 0.2;
double δT = 0.02;
- CiamarraSystem<D> s(L);
+ System<D> s(L);
Rng r;
@@ -58,11 +525,8 @@ int main() {
}
std::cerr << "Found state with appropriate density." << std::endl;
- if (!s.compatible()) {
- std::cerr <<"whoops!" << std::endl;
- }
- CiamarraSystem<D> s0 = s;
+ System<D> s0 = s;
std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
@@ -72,9 +536,16 @@ int main() {
n++;
std::cout << i << " " << s.overlap(s0) << 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);
+// s.swendsenWang(Transformation<D>(L, ms, r), r);
}
return 0;
diff --git a/glass.hpp b/glass.hpp
deleted file mode 100644
index b35964b..0000000
--- a/glass.hpp
+++ /dev/null
@@ -1,677 +0,0 @@
-#include <limits>
-#include <list>
-#include <queue>
-#include <vector>
-
-#include <eigen3/Eigen/Dense>
-
-#include "pcg-cpp/include/pcg_random.hpp"
-#include "randutils/randutils.hpp"
-
-using Rng = randutils::random_generator<pcg32>;
-
-template <int D> using Vector = Eigen::Matrix<int, D, 1>;
-template <int D> using Matrix = Eigen::Matrix<int, D, D>;
-
-int iPow(int x, unsigned p) {
- if (p == 0)
- return 1;
- if (p == 1)
- return x;
-
- int tmp = iPow(x, p / 2);
- if (p % 2 == 0)
- return tmp * tmp;
- else
- return x * tmp * tmp;
-}
-
-unsigned mod(signed a, unsigned b) { return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b; }
-
-template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) {
- Vector<D> u;
- for (unsigned i = 0; i < D; i++) {
- u(i) = mod(v(i), b);
- }
- return u;
-}
-
-template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
- if (level > 0) {
- unsigned new_level = level - 1;
- for (std::array<double, D>& sequence : sequences) {
- std::array<double, D> new_sequence = sequence;
- new_sequence[new_level] = -1;
- sequences.push_front(new_sequence);
- }
- one_sequences<D>(sequences, new_level);
- }
-}
-
-template <int D> std::vector<Matrix<D>> generateTorusMatrices() {
- std::vector<Matrix<D>> mats;
-
- std::array<double, D> ini_sequence;
- ini_sequence.fill(1);
- std::list<std::array<double, D>> sequences;
- sequences.push_back(ini_sequence);
-
- one_sequences<D>(sequences, D);
-
- sequences.pop_back(); // don't want the identity matrix!
-
- for (std::array<double, D> sequence : sequences) {
- Matrix<D> m;
- for (unsigned i = 0; i < D; i++) {
- for (unsigned j = 0; j < D; j++) {
- if (i == j) {
- m(i, j) = sequence[i];
- } else {
- m(i, j) = 0;
- }
- }
- }
-
- mats.push_back(m);
- }
-
- for (unsigned i = 0; i < D; i++) {
- for (unsigned j = 0; j < D; j++) {
- if (i != j) {
- Matrix<D> m;
- for (unsigned k = 0; k < D; k++) {
- for (unsigned l = 0; l < D; l++) {
- if ((k == i && l == j) || (k == j && l == i)) {
- if (i < j) {
- m(k, l) = 1;
- } else {
- m(k, l) = -1;
- }
- } else if (k == l && (k != i && k != j)) {
- m(k, l) = 1;
- } else {
- m(k, l) = 0;
- }
- }
- }
- mats.push_back(m);
- }
- }
- }
-
- return mats;
-}
-
-template <int D> class CiamarraState : public Vector<D> {
-public:
- CiamarraState() : Vector<D>(Vector<D>::Zero()) {}
-
- CiamarraState(const Vector<D>& v) : Vector<D>(v) {}
-
- CiamarraState(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; }
-
- CiamarraState(Rng& r) : CiamarraState(r.uniform(0, D - 1), r.pick({-1, 1})) {}
-
- bool isEmpty() const { return this->squaredNorm() == 0; }
-
- void remove() { this->setZero(); }
-
- CiamarraState<D> flip() const {
- CiamarraState<D> s;
- for (unsigned i = 0; i < D; i++) {
- s(i) = -this->operator()(i);
- }
- return s;
- }
-};
-
-class BiroliState {
-public:
- unsigned type;
- unsigned occupiedNeighbors;
-
- BiroliState() {
- type = std::numeric_limits<unsigned>::max();
- occupiedNeighbors = 0;
- }
-
- BiroliState(unsigned t, unsigned o) {
- type = t;
- occupiedNeighbors = o;
- }
-
- BiroliState(Rng& r) {
- type = r.pick({1,2,2,2,2,2,3,3,3,3});
- }
-
- bool isEmpty() const { return type == std::numeric_limits<unsigned>::max(); }
-
- void remove() { type = std::numeric_limits<unsigned>::max(); };
-};
-
-template <int D> class Transformation {
-public:
- unsigned L;
- 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;
- }
-
- Transformation<D> inverse() const {
- return Transformation<D>(L, m.transpose(), -m.transpose() * v);
- }
-
- Vector<D> apply(const Vector<D>& x) const { 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;
- }
-
- CiamarraState<D> apply(const CiamarraState<D>& s) const { return CiamarraState<D>(m * s); }
-
- BiroliState apply(const BiroliState& s) const { return s; }
-};
-
-template <int D, typename State> class HalfEdge;
-
-template <int D, typename State> class Vertex {
-public:
- Vector<D> position;
- std::vector<HalfEdge<D, State>> adjacentEdges;
- State state;
- bool marked;
-
- bool isEmpty() const { return state.isEmpty(); }
-};
-
-template <int D, class State> class HalfEdge {
-public:
- Vertex<D, State>& neighbor;
- Vector<D> Δx;
-
- HalfEdge(Vertex<D, State>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
-};
-
-template <int D, class State> class System {
-public:
- const unsigned L;
- unsigned N;
- std::vector<Vertex<D, State>> vertices;
- Transformation<D> orientation;
-
- 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) : 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);
- 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, State>(vertices[j], Δx));
- vertices[j].adjacentEdges.push_back(HalfEdge<D, State>(vertices[i], -Δx));
- }
- }
- }
-
- virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&) { return {}; };
-
- virtual bool insert(Vertex<D, State>& v, const State& s, bool force = false) { return false; };
-
- virtual bool remove(Vertex<D, State>& v) { return false;};
-
- virtual bool randomMove(Rng& r) { return false;};
-
- virtual void swap(Vertex<D, State>& v1, Vertex<D, State>& v2) {};
-
- virtual bool exchange(Vertex<D, State>& v1, Vertex<D, State>& v2) { return false;};
-
- bool randomExchange(Rng& r) {
- Vertex<D, State>& v1 = r.pick(vertices);
- Vertex<D, State>& v2 = r.pick(vertices);
-
- return exchange(v1, v2);
- }
-
- bool compatible() {
- for (Vertex<D, State>& v : vertices) {
- if (overlaps(v).size() > 0) {
- return false;
- }
- }
-
- return true;
- }
-
- double density() const { return N / pow(L, D); }
-
- unsigned size() const { return vertices.size(); }
-
- void sweepGrandCanonical(double z, Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
- if (0.5 < r.uniform(0.0, 1.0)) {
- double pIns = size() * z / (N + 1);
-
- if (pIns > r.uniform(0.0, 1.0)) {
- while (true) {
- Vertex<D, State>& v = r.pick(vertices);
- if (v.isEmpty()) {
- insert(v, State(r));
- break;
- }
- }
- }
- } else {
-
- double pDel = N / (z * size());
-
- if (pDel > r.uniform(0.0, 1.0)) {
- remove(r.pick(vertices));
- }
- }
-
- randomMove(r);
- }
- }
-
- void sweepLocal(Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
- randomMove(r);
- }
- }
-
- void sweepSwap(Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
- randomExchange(r);
- }
- }
-
- unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, bool dry = false) {
- unsigned n = 0;
-
- Vertex<D, State>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
-
- if (&v0New != &v0) {
- std::queue<std::reference_wrapper<Vertex<D, State>>> q;
-
- v0.marked = true;
- v0New.marked = true;
-
- swap(v0, v0New);
-
- for (Vertex<D, State>& vn : overlaps(v0)) {
- if (!vn.marked) {
- q.push(vn);
- }
- }
- for (Vertex<D, State>& vn : overlaps(v0New)) {
- if (!vn.marked) {
- q.push(vn);
- }
- }
-
- while (!q.empty()) {
- Vertex<D, State>& v = q.front();
- q.pop();
-
- if (!v.marked && !overlaps(v).empty()) {
- v.marked = true;
- Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.position))];
-
- if (&vNew != &v) {
- vNew.marked = true;
-
- swap(v, vNew);
-
- for (Vertex<D, State>& vn : overlaps(v)) {
- if (!vn.marked) {
- q.push(vn);
- }
- }
- for (Vertex<D, State>& vn : overlaps(vNew)) {
- if (!vn.marked) {
- q.push(vn);
- }
- }
-
- n += 2;
- } else {
- n += 1;
- }
- }
- if (q.empty()) {
- for (Vertex<D, State>& vv : vertices) {
- if (!vv.marked && !overlaps(vv).empty()) {
-// std::cerr << "Found bad state at end" << std::endl;
- q.push(vv);
- }
- }
- }
- }
-
-
- }
-
- if (n > size() / 4) {
- orientation = R.apply(orientation);
- }
-
- for (Vertex<D, State>& v : vertices) {
- v.marked = false;
- }
-
- return n;
- }
-
- void swendsenWang(const Transformation<D>& R, Rng& r) {
- for (Vertex<D, State>& v : vertices) {
- if (!v.marked) {
- bool dry = 0.5 < r.uniform(0.0, 1.0);
- unsigned n = flipCluster(R, v, dry);
- if (n > size() / 4 && !dry) {
- orientation = R.apply(orientation);
- }
- }
- }
-
- for (Vertex<D, State>& v : vertices) {
- v.marked = false;
- }
- }
-
- virtual int overlap(const System<D, State>& s) const { return 0; };
-};
-
-template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
- public:
-
- using System<D, CiamarraState<D>>::System;
- std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>>
- overlaps(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s,
- bool excludeSelf = false) override {
- std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>> o;
-
- if (s.isEmpty()) {
- return o;
- }
-
- if (!v.isEmpty() && !excludeSelf) {
- o.push_back(v);
- }
-
- for (const HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) {
- if (!e.neighbor.isEmpty()) {
- if (s.dot(e.Δx) == 1 || e.neighbor.state.dot(e.Δx) == -1) {
- o.push_back(e.neighbor);
- }
- }
- }
-
- return o;
- }
-
- bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s, bool force = false) override {
- if (force || overlaps(v, s).empty()) {
- v.state = s;
- this->N++;
- return true;
- } else {
- return false;
- }
- }
-
- bool remove(Vertex<D, CiamarraState<D>>& v) override {
- if (v.isEmpty()) {
- return false;
- } else {
- v.state.remove();
- this->N--;
- return true;
- }
- }
-
- bool randomMove(Rng& r) override {
- Vertex<D, CiamarraState<D>>& v = r.pick(this->vertices);
- CiamarraState<D> oldState = v.state;
-
- if (!remove(v)) {
- return false;
- }
-
- if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) {
- for (HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) {
- if (1 == e.Δx.dot(oldState)) {
- if (insert(e.neighbor, oldState.flip())) {
- return true;
- }
- break;
- }
- }
- } else {
- CiamarraState<D> newState(r);
- while (newState.dot(oldState) == 1) {
- newState = CiamarraState<D>(r);
- }
- if (insert(v, newState)) {
- return true;
- }
- }
- v.state = oldState;
- this->N++;
- return false;
- }
-
- void swap(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override {
- std::swap(v1.state, v2.state);
- }
-
- bool exchange(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override {
- if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) {
- swap(v1, v2);
- return true;
- } else {
- return false;
- }
- }
-
- void setGroundState() {
- this->N = 0;
-
- for (Vertex<D, CiamarraState<D>>& v : this->vertices) {
- unsigned a = 0;
- for (unsigned d = 0; d < D; d++) {
- a += (d + 1) * v.position(d);
- }
- a %= 2 * D + 1;
-
- v.state.setZero() = Vector<D>::Zero();
-
- if (0 < a && a <= D) {
- v.state(a - 1) = -1;
- this->N++;
- } else if (D < a) {
- v.state(2 * D - a) = 1;
- this->N++;
- }
- }
- }
-
- int overlap(const System<D, CiamarraState<D>>& s) const override {
- int o = 0;
-
- for (unsigned i = 0; i < this->vertices.size(); i++) {
- CiamarraState<D> s2 = this->orientation.apply(
- s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state);
- o += this->vertices[i].state.dot(s2);
- }
-
- return o;
- }
-};
-
-template <int D> class BiroliSystem : public System<D, BiroliState> {
-public:
- using System<D, BiroliState>::System;
-
- bool canReplaceWith(const Vertex<D, BiroliState>& v, const BiroliState& s) {
- if (s.isEmpty()) {
- return true;
- }
- if (s.type < v.state.occupiedNeighbors) {
- return false;
- } else {
- int Δn = 0;
- if (v.isEmpty()) {
- Δn += 1;
- }
- for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
- if (e.neighbor.state.type < e.neighbor.state.occupiedNeighbors + Δn) {
- return false;
- }
- }
- }
- return true;
- }
-
- std::list<std::reference_wrapper<Vertex<D, BiroliState>>>
- overlaps(Vertex<D, BiroliState>& v) override {
- std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o;
-
- if (v.isEmpty()) { // an empty site cannot be frustrating anyone
- return o;
- }
-
- bool selfFrustrated = v.state.occupiedNeighbors > v.state.type;
- bool anyNeighborFrustrated = false;
-
- for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
- if (!e.neighbor.isEmpty()) {
- bool thisNeighborFrustrated = e.neighbor.state.type < e.neighbor.state.occupiedNeighbors;
-
- if (thisNeighborFrustrated) {
- anyNeighborFrustrated = true;
- }
-
- if (selfFrustrated || thisNeighborFrustrated) {
- o.push_back(e.neighbor);
- }
- }
- }
-
- if (selfFrustrated || anyNeighborFrustrated) {
- o.push_back(v);
- }
-
- return o;
- }
-
- bool insert(Vertex<D, BiroliState>& v, const BiroliState& s, bool force = false) override {
- if (force || (v.isEmpty() && canReplaceWith(v, s))) {
- v.state.type = s.type;
- for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
- e.neighbor.state.occupiedNeighbors++;
- }
- this->N++;
- return true;
- } else {
- return false;
- }
- }
-
- bool remove(Vertex<D, BiroliState>& v) override {
- if (v.isEmpty()) {
- return false;
- } else {
- v.state.remove();
- for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
- e.neighbor.state.occupiedNeighbors--;
- }
- this->N--;
- return true;
- }
- }
-
- bool randomMove(Rng& r) override {
- Vertex<D, BiroliState>& v = r.pick(this->vertices);
-
- return exchange(v, r.pick(v.adjacentEdges).neighbor);
- }
-
- void swap(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) override {
- if (v1.state.type != v2.state.type) {
- if (!v1.isEmpty() && !v2.isEmpty()) {
- std::swap(v1.state.type, v2.state.type);
- } else if (v1.isEmpty()) {
- unsigned t = v2.state.type;
- remove(v2);
- insert(v1, BiroliState(t, 0), true);
- } else {
- unsigned t = v1.state.type;
- remove(v1);
- insert(v2, BiroliState(t, 0), true);
- }
- }
- }
-
- bool exchange(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) override {
- if (canReplaceWith(v1, v2.state) && canReplaceWith(v2, v1.state)) {
- swap(v1, v2);
- return true;
- } else {
- return false;
- }
- }
-
- int overlap(const System<D, BiroliState>& s) const override {
- int o = 0;
-
- for (unsigned i = 0; i < this->vertices.size(); i++) {
- BiroliState s2 = this->orientation.apply(
- s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state);
- if (s2.type > 0 && s2.type == this->vertices[i].state.type) {
- o++;
- } else {
- o--;
- }
- }
-
- return o;
- }
-};
--
cgit v1.2.3-70-g09d2