summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--biroli-mezard.cpp214
-rw-r--r--glass.hpp4
2 files changed, 93 insertions, 125 deletions
diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index 5b6213f..52a4360 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -24,41 +24,9 @@ public:
void remove() {
maximumNeighbors = Empty;
}
-}
-
-template <int D> class Vertex {
-public:
- Vector<D> initialPosition;
- Vector<D> position;
- std::vector<std::reference_wrapper<Vertex<D>>> neighbors;
- unsigned occupiedNeighbors;
- unsigned maximumNeighbors;
- bool marked;
- bool visited;
-
- Vertex() {
- marked = false;
- }
-
- 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;
- }
- void setInitial() {
- initialPosition = position;
- }
};
-template <int D> class BiroliSystem : public HardSystem<D, BiroliState> {
+template <int D> class BiroliSystem : public System<D, BiroliState> {
public:
std::vector<unsigned> N;
std::unordered_set<Vertex<D, BiroliState>*> occupants;
@@ -69,11 +37,23 @@ public:
return BiroliSystem::size() - N[0];
}
- double density() const { return (double)occupancy() / size(); }
+ double density() const { return (double)occupancy() / BiroliSystem::size(); }
+
+ bool valid(const Vertex<D, BiroliState>& v) const {
+ unsigned o = 0;
+
+ for (const Vertex<D, BiroliState>& vn : v.neighbors) {
+ if (!vn.empty()) {
+ o++;
+ }
+ }
+
+ return o == v.occupiedNeighbors;
+ }
bool compatible() const {
- for (const Vertex<D>& v : vertices) {
- if (v.frustrated()) {
+ for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
+ if (v.state.frustrated()) {
return false;
}
}
@@ -82,8 +62,8 @@ public:
}
bool valid() const {
- for (const Vertex<D>& v : vertices) {
- if (!v.valid()) {
+ for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
+ if (!valid(v)) {
return false;
}
}
@@ -91,25 +71,12 @@ public:
return true;
}
- BiroliSystem(unsigned L, unsigned T) : L(L), N(T + 1, 0), vertices(iPow(L, D)), orientation(L) {
- N[0] = size();
-
- 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 < 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]);
- }
- }
+ BiroliSystem(unsigned L, unsigned T) : System<D, BiroliState>(L), N(T + 1, 0) {
+ N[0] = BiroliSystem::size();
}
- std::list<std::reference_wrapper<Vertex<D>>> overlaps(Vertex<D>& v) {
- std::list<std::reference_wrapper<Vertex<D>>> o;
+ std::list<std::reference_wrapper<Vertex<D, BiroliState>>> overlaps(Vertex<D, BiroliState>& v) {
+ std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o;
if (v.empty()) { // an empty site cannot be frustrating anyone
return o;
@@ -118,7 +85,7 @@ public:
bool selfFrustrated = v.frustrated();
bool anyNeighborFrustrated = false;
- for (Vertex<D>& vn : v.neighbors) {
+ for (Vertex<D, BiroliState>& vn : v.neighbors) {
if (!vn.empty()) {
bool thisNeighborFrustrated = vn.frustrated();
@@ -139,28 +106,29 @@ public:
return o;
}
- std::list<std::reference_wrapper<Vertex<D>>> overlaps(const std::list<std::reference_wrapper<Vertex<D>>>& vs) {
- std::list<std::reference_wrapper<Vertex<D>>> o;
+ std::list<std::reference_wrapper<Vertex<D, BiroliState>>> overlaps(const std::list<std::reference_wrapper<Vertex<D, BiroliState>>>& vs) {
+ std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o;
- for (Vertex<D>& v : vs) {
+ for (Vertex<D, BiroliState>& v : vs) {
o.splice(o.begin(), overlaps(v));
}
return o;
}
- bool canReplaceWith(const Vertex<D>& v, unsigned t) const {
+ bool canReplaceWith(const Vertex<D, BiroliState>& v, unsigned t) const {
if (t == Empty) {
return true;
}
- if (t < v.occupiedNeighbors) {
+ if (t < v.state.occupiedNeighbors) {
return false;
}
if (v.empty()) {
- for (const Vertex<D>& vn : v.neighbors) {
- if (vn.maximumNeighbors < vn.occupiedNeighbors + 1) {
+ for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+ const Vertex<D, BiroliState>& vn = e.neighbor;
+ if (vn.state.maximumNeighbors < vn.state.occupiedNeighbors + 1) {
return false;
}
}
@@ -170,16 +138,16 @@ public:
}
void setInitial() {
- for (Vertex<D>& v : vertices) {
+ for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
v.setInitial();
}
}
- bool insert(Vertex<D>& v, unsigned t, bool force = false) {
+ bool insert(Vertex<D, BiroliState>& v, unsigned t, bool force = false) {
if (force || (v.empty() && canReplaceWith(v, t))) {
- v.maximumNeighbors = t;
- for (Vertex<D>& vn : v.neighbors) {
- vn.occupiedNeighbors++;
+ v.state.maximumNeighbors = t;
+ for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+ e.neighbor.state.occupiedNeighbors++;
}
N[t]++;
N[0]--;
@@ -192,39 +160,39 @@ public:
}
}
- bool remove(Vertex<D>& v) {
+ bool remove(Vertex<D, BiroliState>& v) {
if (v.empty()) {
return false;
} else {
- if (v.maximumNeighbors > 0) {
+ if (v.state.maximumNeighbors > 0) {
occupants.erase(&v);
}
- N[v.maximumNeighbors]--;
+ N[v.state.maximumNeighbors]--;
N[0]++;
- v.maximumNeighbors = Empty;
- for (Vertex<D>& vn : v.neighbors) {
- vn.occupiedNeighbors--;
+ v.state.maximumNeighbors = Empty;
+ for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+ e.neighbor.state.occupiedNeighbors--;
}
return true;
}
}
bool randomLocalExchange(Rng& r) {
- Vertex<D>& v = r.pick(vertices);
+ Vertex<D, BiroliState>& v = r.pick(BiroliSystem::vertices);
return exchange(v, r.pick(v.neighbors));
}
- void swap(Vertex<D>& v1, Vertex<D>& v2) {
- if (v1.maximumNeighbors != v2.maximumNeighbors) {
+ void swap(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) {
+ if (v1.state.maximumNeighbors != v2.state.maximumNeighbors) {
if (!v1.empty() && !v2.empty()) {
- std::swap(v1.maximumNeighbors, v2.maximumNeighbors);
+ std::swap(v1.state.maximumNeighbors, v2.state.maximumNeighbors);
} else if (v1.empty()) {
- unsigned t = v2.maximumNeighbors;
+ unsigned t = v2.state.maximumNeighbors;
remove(v2);
insert(v1, t, true);
} else {
- unsigned t = v1.maximumNeighbors;
+ unsigned t = v1.state.maximumNeighbors;
remove(v1);
insert(v2, t, true);
}
@@ -232,8 +200,8 @@ public:
}
}
- bool exchange(Vertex<D>& v1, Vertex<D>& v2) {
- if (canReplaceWith(v1, v2.maximumNeighbors) && canReplaceWith(v2, v1.maximumNeighbors)) {
+ bool exchange(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) {
+ if (canReplaceWith(v1, v2.state.maximumNeighbors) && canReplaceWith(v2, v1.state.maximumNeighbors)) {
swap(v1, v2);
return true;
} else {
@@ -246,10 +214,10 @@ public:
using namespace std::complex_literals;
- for (const Vertex<D>& v : vertices) {
+ for (const Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
if (v.maximumNeighbors == t) {
for (unsigned d = 0; d < D; d++) {
- F += exp(1i * 2.0 * M_PI * (double)k * (double)(v.position(d) - v.initialPosition(d)) / (double)L);
+ F += exp(1i * 2.0 * M_PI * (double)k * (double)(v.position(d) - v.initialPosition(d)) / (double)BiroliSystem::L);
}
}
}
@@ -259,11 +227,11 @@ public:
return F;
}
- double overlap(const System<D>& s, unsigned t = 2) const {
+ double overlap(const BiroliSystem<D>& s, unsigned t = 2) const {
int o = 0;
- for (const Vertex<D>& v1 : vertices) {
- const Vertex<D>& v2 = s.vertices[vectorToIndex(orientation.inverse().apply(v1.position))];
+ for (const Vertex<D, BiroliState>& v1 : BiroliSystem::vertices) {
+ const Vertex<D, BiroliState>& v2 = s.BiroliSystem::vertices[vectorToIndex(BiroliSystem::orientation.inverse().apply(v1.position))];
if (v1.maximumNeighbors == v2.maximumNeighbors) {
o++;
}
@@ -275,12 +243,12 @@ public:
nn += iPow(n, 2);
}
- return ((double)o / (double)size() - (double)nn / pow((double)size(), 2)) / (1 - (double)nn / pow((double)size(), 2));
+ return ((double)o / (double)BiroliSystem::size() - (double)nn / pow((double)BiroliSystem::size(), 2)) / (1 - (double)nn / pow((double)BiroliSystem::size(), 2));
}
bool randomExchange(Rng& r) {
- Vertex<D>& v1 = r.pick(vertices);
- Vertex<D>& v2 = r.pick(vertices);
+ Vertex<D, BiroliState>& v1 = r.pick(BiroliSystem::vertices);
+ Vertex<D, BiroliState>& v2 = r.pick(BiroliSystem::vertices);
return exchange(v1, v2);
}
@@ -290,24 +258,24 @@ public:
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);
+ double pIns = BiroliSystem::size() * z / (occupancy() + 1);
if (pIns > r.uniform(0.0, 1.0)) {
- Vertex<D>& v = r.pick(vertices);
+ Vertex<D, BiroliState>& v = r.pick(BiroliSystem::vertices);
insert(v, t);
}
} else {
double pDel = density() / z;
if (pDel > r.uniform(0.0, 1.0)) {
- Vertex<D>* v = r.pick(occupants);
+ Vertex<D, BiroliState>* v = r.pick(occupants);
remove(*v);
}
}
}
void sweepGrandCanonical(double z, Rng& r) {
- for (unsigned i = 0; i < size(); i++) {
+ for (unsigned i = 0; i < BiroliSystem::size(); i++) {
stepGrandCanonical(z, r);
randomExchange(r);
@@ -315,36 +283,36 @@ public:
}
void sweepLocal(Rng& r) {
- for (unsigned i = 0; i < size() / 2; i++) {
+ for (unsigned i = 0; i < BiroliSystem::size() / 2; i++) {
randomLocalExchange(r);
}
}
void sweepSwap(Rng& r) {
- for (unsigned i = 0; i < size() / 2; i++) {
+ for (unsigned i = 0; i < BiroliSystem::size() / 2; i++) {
randomExchange(r);
}
}
- Vertex<D>& apply(const Transformation<D>& R, const Vertex<D>& v) {
- return vertices[vectorToIndex(R.apply(v.position))];
+ Vertex<D, BiroliState>& apply(const Transformation<D>& R, const Vertex<D, BiroliState>& v) {
+ return BiroliSystem::vertices[vectorToIndex(R.apply(v.position))];
}
- bool eventChain(const Transformation<D>& R, Vertex<D>& v0) {
+ bool eventChain(const Transformation<D>& R, Vertex<D, BiroliState>& v0) {
unsigned n = 0;
if (!v0.empty()) {
- Vertex<D>& v0New = apply(R, v0);;
+ Vertex<D, BiroliState>& v0New = apply(R, v0);;
unsigned t0 = v0.maximumNeighbors;
- std::queue<std::pair<std::reference_wrapper<Vertex<D>>, unsigned>> q;
+ std::queue<std::pair<std::reference_wrapper<Vertex<D, BiroliState>>, unsigned>> q;
q.push({v0New, t0});
remove(v0);
while (!q.empty()) {
auto [vr, t] = q.front();
- Vertex<D>& v = vr;
+ Vertex<D, BiroliState>& v = vr;
q.pop();
if (!v.empty()) {
@@ -355,7 +323,7 @@ public:
insert(v, t, true);
v.marked = true;
- for (Vertex<D>& vn : overlaps({v})) {
+ for (Vertex<D, BiroliState>& vn : overlaps({v})) {
if (!vn.marked) {
q.push({apply(R, vn), vn.maximumNeighbors});
remove(vn);
@@ -365,7 +333,7 @@ public:
* frustrated, it may be due to one of its neighbors!
*/
if (&vn != &v) {
- for (Vertex<D>& vnn : overlaps(vn)) {
+ for (Vertex<D, BiroliState>& vnn : overlaps(vn)) {
if (!vnn.marked) {
q.push({apply(R, vnn), vnn.maximumNeighbors});
remove(vnn);
@@ -378,7 +346,7 @@ public:
}
}
- for (Vertex<D>& v : vertices) {
+ for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
v.marked = false;
}
@@ -386,44 +354,44 @@ public:
}
bool randomEventChain(Rng& r) {
- Transformation<D> R(L);
+ Transformation<D> R(BiroliSystem::L);
R.v[r.uniform((unsigned)0, (unsigned)D-1)] = r.pick({-1, 1});
- return eventChain(R, r.pick(vertices));
+ return eventChain(R, r.pick(BiroliSystem::vertices));
}
- unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0) {
+ unsigned flipCluster(const Transformation<D>& R, Vertex<D, BiroliState>& v0) {
unsigned n = 0;
- Vertex<D>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
+ Vertex<D, BiroliState>& v0New = BiroliSystem::vertices[vectorToIndex(R.apply(v0.position))];
if (&v0New != &v0) {
- std::queue<std::reference_wrapper<Vertex<D>>> q;
+ std::queue<std::reference_wrapper<Vertex<D, BiroliState>>> q;
v0.marked = true;
v0New.marked = true;
swap(v0, v0New);
- for (Vertex<D>& vn : overlaps({v0, v0New})) {
+ for (Vertex<D, BiroliState>& vn : overlaps({v0, v0New})) {
if (!vn.marked) {
q.push(vn);
}
}
while (!q.empty()) {
- Vertex<D>& v = q.front();
+ Vertex<D, BiroliState>& v = q.front();
q.pop();
if (!v.marked && !overlaps(v).empty()) {
v.marked = true;
- Vertex<D>& vNew = vertices[vectorToIndex(R.apply(v.position))];
+ Vertex<D, BiroliState>& vNew = BiroliSystem::vertices[vectorToIndex(R.apply(v.position))];
if (&vNew != &v) {
vNew.marked = true;
swap(v, vNew);
- for (Vertex<D>& vn : overlaps({v, vNew})) {
+ for (Vertex<D, BiroliState>& vn : overlaps({v, vNew})) {
if (!vn.marked) {
q.push(vn);
} else {
@@ -431,7 +399,7 @@ public:
* frustrated, it may be due to one of its neighbors!
*/
if (&vn != &v && &vn != &vNew) {
- for (Vertex<D>& vnn : overlaps(vn)) {
+ for (Vertex<D, BiroliState>& vnn : overlaps(vn)) {
if (!vnn.marked) {
q.push(vnn);
}
@@ -449,11 +417,11 @@ public:
n = 1;
}
- if (n > size() / 2) {
- orientation = R.apply(orientation);
+ if (n > BiroliSystem::size() / 2) {
+ BiroliSystem::orientation = R.apply(BiroliSystem::orientation);
}
- for (Vertex<D>& v : vertices) {
+ for (Vertex<D, BiroliState>& v : BiroliSystem::vertices) {
v.marked = false;
}
@@ -461,10 +429,10 @@ public:
}
};
-void print(const System<2>& s) {
- for (const Vertex<2>& v : s.vertices) {
+void print(const BiroliSystem<2>& s) {
+ for (const Vertex<2, BiroliState>& v : s.vertices) {
if (!v.empty()) {
- std::cerr << v.maximumNeighbors;
+ std::cerr << v.state.maximumNeighbors;
} else {
std::cerr << " ";
}
@@ -484,14 +452,14 @@ int main() {
double μmax = 10;
double dμ = 0.05;
- System<D> s(L, 3);
+ BiroliSystem<D> s(L, 3);
Rng r;
for (double μ = μmin; μ <= μmax; μ += dμ) {
double z = exp(μ);
for (unsigned i = 0; i < 1e4; i++) {
- for (unsigned j = 0; j < s.size(); j++) {
+ for (unsigned j = 0; j < s.BiroliSystem::size(); j++) {
s.stepGrandCanonical(z, r);
// s.randomEventChain(r);
s.randomExchange(r);
diff --git a/glass.hpp b/glass.hpp
index d313ced..fc7b186 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -573,7 +573,7 @@ public:
bool dry = 0.5 < r.uniform(0.0, 1.0);
unsigned n = flipCluster(R, v, r, dry);
if (n > pow(HardSystem::L, D) / 4 && !dry) {
- orientation = R.apply(orientation);
+ HardSystem::orientation = R.apply(HardSystem::orientation);
}
}
}
@@ -587,7 +587,7 @@ public:
int o = 0;
for (unsigned i = 0; i < HardSystem::vertices.size(); i++) {
- S s2 = orientation.apply(s.vertices[HardSystem::vectorToIndex(orientation.inverse().apply(HardSystem::indexToVector(i)))].state);
+ S s2 = HardSystem::orientation.apply(s.vertices[HardSystem::vectorToIndex(HardSystem::orientation.inverse().apply(HardSystem::indexToVector(i)))].state);
o += HardSystem::vertices[i].state.dot(s2);
}