summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--biroli-mezard.cpp103
1 files changed, 64 insertions, 39 deletions
diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index 35360d0..c1ed606 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -159,6 +159,17 @@ public:
bool empty() const { return maximumNeighbors == Empty; }
bool frustrated() const { return occupiedNeighbors > maximumNeighbors; }
+ bool valid() const {
+ unsigned o = 0;
+
+ for (const Vertex<D>& v : neighbors) {
+ if (!v.empty()) {
+ o++;
+ }
+ }
+
+ return o == occupiedNeighbors;
+ }
};
template <int D> class System {
@@ -168,6 +179,36 @@ public:
std::vector<Vertex<D>> vertices;
Transformation<D> orientation;
+ unsigned size() const { return vertices.size(); }
+
+ unsigned types() const { return N.size() - 1; }
+
+ unsigned occupancy() const {
+ return size() - N[0];
+ }
+
+ double density() const { return (double)occupancy() / size(); }
+
+ bool compatible() const {
+ for (const Vertex<D>& v : vertices) {
+ if (v.frustrated()) {
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+ bool valid() const {
+ for (const Vertex<D>& v : vertices) {
+ if (!v.valid()) {
+ return false;
+ }
+ }
+
+ return true;
+ }
+
unsigned vectorToIndex(const Vector<D>& x) const {
unsigned i = 0;
for (unsigned d = 0; d < D; d++) {
@@ -187,13 +228,13 @@ public:
System(unsigned L, unsigned T) : L(L), N(T + 1, 0), vertices(iPow(L, D)), orientation(L) {
N[0] = size();
- for (unsigned i = 0; i < iPow(L, D); i++) {
+ for (unsigned i = 0; i < size(); i++) {
vertices[i].position = indexToVector(i);
vertices[i].neighbors.reserve(2 * D);
}
for (unsigned d = 0; d < D; d++) {
- for (signed i = 0; i < iPow(L, D); i++) {
+ for (signed i = 0; i < size(); i++) {
unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1));
vertices[i].neighbors.push_back(vertices[j]);
vertices[j].neighbors.push_back(vertices[i]);
@@ -334,34 +375,8 @@ public:
return exchange(v1, v2);
}
- bool compatible() const {
- for (const Vertex<D>& v : vertices) {
- if (v.frustrated()) {
- return false;
- }
- }
-
- return true;
- }
-
- unsigned occupancy() const {
- unsigned m = 0;
-
- for (unsigned i = 1; i < N.size(); i++) {
- m += N[i];
- }
-
- return m;
- }
-
- unsigned types() const { return N.size() - 1; }
-
- double density() const { return (double)occupancy() / size(); }
-
- unsigned size() const { return vertices.size(); }
-
void sweepGrandCanonical(double z, Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
+ for (unsigned i = 0; i < size(); i++) {
if (0.5 < r.uniform(0.0, 1.0)) {
double pIns = size() * z / (occupancy() + 1);
@@ -388,13 +403,13 @@ public:
}
void sweepLocal(Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
+ for (unsigned i = 0; i < size(); i++) {
randomLocalExchange(r);
}
}
void sweepSwap(Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
+ for (unsigned i = 0; i < size(); i++) {
randomExchange(r);
}
}
@@ -441,6 +456,7 @@ public:
q.push(vn);
}
}
+
for (Vertex<D>& vn : overlaps(vNew)) {
if (!vn.marked) {
q.push(vn);
@@ -452,18 +468,21 @@ public:
n += 1;
}
}
+
+ // For some reason, the process above misses frustrated states. This is an inelegant stopgap.
if (q.empty()) {
for (Vertex<D>& vv : vertices) {
if (!vv.marked && !overlaps(vv).empty()) {
- // std::cerr << "Found bad state at end" << std::endl;
q.push(vv);
}
}
}
}
+ } else {
+ n = 1;
}
- if (n > size() / 4) {
+ if (n > occupancy() / 2) {
orientation = R.apply(orientation);
}
@@ -502,14 +521,14 @@ int main() {
Rng r;
- double z = exp(1 / 0.2);
+ double z = exp(1 / 0.15);
if (!s.compatible()) {
std::cerr << "Storted incompatible!" << std::endl;
return 1;
}
- while (s.density() < 0.57) {
+ while (s.density() < 0.58) {
s.sweepGrandCanonical(z, r);
}
@@ -537,10 +556,11 @@ int main() {
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);
- // s.swendsenWang(Transformation<D>(L, ms, r), r);
+ /*
+ s.sweepLocal(r);
+ nC += s.size();
+ s.sweepSwap(r);
+ */
i++;
}
@@ -549,6 +569,11 @@ int main() {
return 1;
}
+ if (!s.valid()) {
+ std::cerr << "Not valid!" << std::endl;
+ return 1;
+ }
+
std::ofstream file("dist.dat");
for (unsigned i : clusterDist) {
file << i << " ";