summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--biroli-mezard.cpp178
-rw-r--r--ciamarra.cpp468
-rw-r--r--compile_flags.txt3
-rw-r--r--distinguishable.cpp6
-rw-r--r--glass.hpp209
5 files changed, 268 insertions, 596 deletions
diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index e2fb71a..5b6213f 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -1,149 +1,30 @@
#include <fstream>
#include <iostream>
-#include <limits>
-#include <list>
-#include <queue>
-#include <vector>
-#include <chrono>
#include <unordered_set>
-#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>;
+#include "glass.hpp"
const unsigned Empty = std::numeric_limits<unsigned>::max();
-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 Transformation {
+class BiroliState {
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);
- }
+ unsigned maximumNeighbors;
+ unsigned occupiedNeighbors;
- v = v - m * v;
+ BiroliState() {
+ maximumNeighbors = Empty;
+ occupiedNeighbors = 0;
}
- Transformation<D> inverse() const {
- return Transformation<D>(L, m.transpose(), -m.transpose() * v);
+ bool empty() const { return maximumNeighbors == Empty; }
+ bool frustrated() const { return occupiedNeighbors > maximumNeighbors; }
+ bool operator==(const BiroliState& s) const {
+ return s.maximumNeighbors == maximumNeighbors;
}
-
- 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;
+ void remove() {
+ maximumNeighbors = Empty;
}
-};
+}
template <int D> class Vertex {
public:
@@ -156,8 +37,6 @@ public:
bool visited;
Vertex() {
- occupiedNeighbors = 0;
- maximumNeighbors = Empty;
marked = false;
}
@@ -179,20 +58,15 @@ public:
}
};
-template <int D> class System {
+template <int D> class BiroliSystem : public HardSystem<D, BiroliState> {
public:
- const unsigned L;
std::vector<unsigned> N;
- std::unordered_set<Vertex<D>*> occupants;
- std::vector<Vertex<D>> vertices;
- Transformation<D> orientation;
-
- unsigned size() const { return vertices.size(); }
+ std::unordered_set<Vertex<D, BiroliState>*> occupants;
unsigned types() const { return N.size() - 1; }
unsigned occupancy() const {
- return size() - N[0];
+ return BiroliSystem::size() - N[0];
}
double density() const { return (double)occupancy() / size(); }
@@ -217,23 +91,7 @@ public:
return true;
}
- 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, unsigned T) : L(L), N(T + 1, 0), vertices(iPow(L, D)), orientation(L) {
+ 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++) {
diff --git a/ciamarra.cpp b/ciamarra.cpp
index f2f44c1..7d7c8af 100644
--- a/ciamarra.cpp
+++ b/ciamarra.cpp
@@ -1,251 +1,58 @@
-#include <eigen3/Eigen/Dense>
-#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
-#include <iostream>
-#include <list>
-#include <vector>
-#include <queue>
-
#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;
-}
+#include "glass.hpp"
-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> class CiamarraState : public Vector<D> {
+public:
+ CiamarraState() : Vector<D>(Vector<D>::Zero()) {}
-template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
- std::vector<Matrix<D>> mats;
+ CiamarraState(const Vector<D>& v) : Vector<D>(v) {}
- std::array<double, D> ini_sequence;
- ini_sequence.fill(1);
- std::list<std::array<double, D>> sequences;
- sequences.push_back(ini_sequence);
+ CiamarraState(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { CiamarraState::operator()(a) = b; }
- one_sequences<D>(sequences, D);
+ CiamarraState(Rng& r) : CiamarraState(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {}
- sequences.pop_back(); // don't want the identity matrix!
+ bool empty() const { return CiamarraState::squaredNorm() == 0; }
- 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;
- }
- }
- }
+ void remove() { CiamarraState::setZero(); }
- mats.push_back(m);
+ bool operator==(const CiamarraState<D>& s) const {
+ return CiamarraState::dot(s) == 1;
}
- 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;
+ CiamarraState<D> flip() const {
+ CiamarraState<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;
+ CiamarraState<D> transform(const Transformation<D>& m) const {
+ return m * *this;
}
+};
- 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;
- }
+template <unsigned D>
+class CiamarraSystem : public HardSystem<D, CiamarraState<D>> {
+ using HardSystem<D, CiamarraState<D>>::HardSystem;
- 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;
+ 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()) {
+ if (s.empty()) {
return o;
}
- if (!v.isEmpty() && !excludeSelf) {
+ if (!v.empty() && !excludeSelf) {
o.push_back(v);
}
- for (const HalfEdge<D>& e : v.adjacentEdges) {
- if (!e.neighbor.isEmpty()) {
+ for (const HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) {
+ if (!e.neighbor.empty()) {
if (s.dot(e.Δx) == 1 || e.neighbor.state.dot(e.Δx) == -1) {
o.push_back(e.neighbor);
}
@@ -255,77 +62,10 @@ public:
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;
+ CiamarraSystem::N = 0;
- for (Vertex<D>& v : vertices) {
+ for (Vertex<D, CiamarraState<D>>& v : CiamarraSystem::vertices) {
unsigned a = 0;
for (unsigned d = 0; d < D; d++) {
a += (d + 1) * v.position(d);
@@ -336,142 +76,18 @@ public:
if (0 < a && a <= D) {
v.state(a - 1) = -1;
- N++;
+ CiamarraSystem::N++;
} else if (D < a) {
v.state(2 * D - a) = 1;
- N++;
+ CiamarraSystem::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) {
+
+void print(const CiamarraSystem<2>& s) {
+ for (const Vertex<2, CiamarraState<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) {
@@ -492,14 +108,6 @@ void print(const System<2>& s) {
}
}
-template <unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
- Vector<D> x;
- for (unsigned i = 0; i < D; i++) {
- x[i] = r.uniform((unsigned)0, L - 1);
- }
- return x;
-}
-
int main() {
const unsigned D = 3;
unsigned L = 15;
@@ -509,7 +117,7 @@ int main() {
double Tmax = 0.2;
double δT = 0.02;
- System<D> s(L);
+ CiamarraSystem<D> s(L);
Rng r;
@@ -526,7 +134,7 @@ int main() {
std::cerr << "Found state with appropriate density." << std::endl;
- System<D> s0 = s;
+ CiamarraSystem<D> s0 = s;
std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
@@ -537,10 +145,10 @@ int main() {
std::cout << i << " " << s.overlap(s0) << std::endl;
}
Matrix<D> m = r.pick(ms);
- Vertex<D>& v = r.pick(s.vertices);
+ Vertex<D, CiamarraState<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) {
+ for (Vertex<D, CiamarraState<D>>& v : s.vertices) {
v.marked = false;
}
// s.sweepLocal(r);
diff --git a/compile_flags.txt b/compile_flags.txt
new file mode 100644
index 0000000..02d9f07
--- /dev/null
+++ b/compile_flags.txt
@@ -0,0 +1,3 @@
+-std=c++20
+-O3
+-I/usr/include/eigen3
diff --git a/distinguishable.cpp b/distinguishable.cpp
index 0b70cf2..22f228a 100644
--- a/distinguishable.cpp
+++ b/distinguishable.cpp
@@ -16,12 +16,14 @@ public:
};
template <unsigned D>
-class DistinguishableSystem : public FiniteEnergySystem<D, DistinguishableState> {
+class DistinguishableSystem : public SoftSystem<D, DistinguishableState> {
public:
std::vector<double> interaction;
DistinguishableSystem(unsigned L, unsigned N, Rng& r)
- : FiniteEnergySystem<D, DistinguishableState>(L, N), interaction(N * N) {
+ : SoftSystem<D, DistinguishableState>(L), interaction(N * N) {
+ DistinguishableSystem::N = N;
+
for (double& V : interaction) {
V = r.uniform(-0.5, 0.5);
}
diff --git a/glass.hpp b/glass.hpp
index 64fe8b7..d313ced 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -180,6 +180,7 @@ public:
const unsigned L;
unsigned N;
std::vector<Vertex<D, S>> vertices;
+ Transformation<D> orientation;
unsigned vectorToIndex(const Vector<D>& x) const {
unsigned i = 0;
@@ -197,7 +198,7 @@ public:
return x;
}
- System(unsigned L, unsigned N) : L(L), N(N), 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].initialPosition = vertices[i].position;
@@ -216,8 +217,14 @@ public:
}
}
+ unsigned size() const { return vertices.size(); }
+
double density() const { return N / pow(L, D); }
+ unsigned maxOccupation() const {
+ return iPow(L, D);
+ }
+
void setInitialPosition() {
for (Vertex<D, S>& v : vertices) {
v.initialPosition = v.position;
@@ -240,7 +247,7 @@ public:
}
};
-template <unsigned D, State S> class FiniteEnergySystem : public System<D, S> {
+template <unsigned D, State S> class SoftSystem : public System<D, S> {
public:
using System<D, S>::System;
@@ -326,7 +333,7 @@ public:
bool dry = false) {
std::queue<std::array<std::reference_wrapper<Vertex<D, S>>, 2>> q;
Vector<D> x0New = R.apply(v0.position);
- Vertex<D, S>& v0New = this->vertices[this->vectorToIndex(x0New)];
+ Vertex<D, S>& v0New = SoftSystem::vertices[SoftSystem::vectorToIndex(x0New)];
q.push({v0, v0New});
unsigned n = 0;
@@ -346,7 +353,7 @@ public:
Vertex<D, S>& vn = e.neighbor;
Vector<D> xnNew = R.apply(vn.position);
- Vertex<D, S>& vnNew = this->vertices[this->vectorToIndex(xnNew)];
+ Vertex<D, S>& vnNew = SoftSystem::vertices[SoftSystem::vectorToIndex(xnNew)];
if (!vn.marked && !vnNew.marked) {
double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state);
@@ -393,3 +400,197 @@ public:
}
};
+template <unsigned D, State S> class HardSystem : public System<D, S> {
+public:
+
+ HardSystem(unsigned L) : System<D, S>(L) {}
+
+ virtual std::list<std::reference_wrapper<Vertex<D, S>>> overlaps(Vertex<D, S>&, const S&,
+ bool = false) { return {}; }
+
+ bool insert(Vertex<D, S>& v, const S& s) {
+ if (overlaps(v, s).empty()) {
+ v.state = s;
+ HardSystem::N++;
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ bool tryDeletion(Vertex<D, S>& v) {
+ if (v.empty()) {
+ return false;
+ } else {
+ v.state.remove();
+ HardSystem::N--;
+ return true;
+ }
+ }
+
+ bool tryRandomMove(Rng& r) {
+ Vertex<D, S>& v = r.pick(HardSystem::vertices);
+ S oldState = v.state;
+
+ if (!tryDeletion(v)) {
+ return false;
+ }
+
+ if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) {
+ for (HalfEdge<D, S>& e : v.adjacentEdges) {
+ if (1 == e.Δx.dot(oldState)) {
+ if (insert(e.neighbor, oldState.flip())) {
+ return true;
+ }
+ break;
+ }
+ }
+ } else {
+ S newState(r);
+ while (newState == oldState) {
+ newState = S(r);
+ }
+ if (insert(v, newState)) {
+ return true;
+ }
+ }
+ v.state = oldState;
+ HardSystem::N++;
+ return false;
+ }
+
+ bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& 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, S>& v1 = r.pick(HardSystem::vertices);
+ Vertex<D, S>& v2 = r.pick(HardSystem::vertices);
+
+ return trySwap(v1, v2);
+ }
+
+
+ bool compatible() {
+ for (Vertex<D, S>& v : HardSystem::vertices) {
+ if (overlaps(v, v.state, true).size() > 0) {
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+ void sweepGrandCanonical(double z, Rng& r) {
+ for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
+ if (0.5 < r.uniform(0.0, 1.0)) {
+ double pIns = HardSystem::maxOccupation() * z / (HardSystem::N + 1);
+
+ if (pIns > r.uniform(0.0, 1.0)) {
+ while (true) {
+ Vertex<D, S>& v = r.pick(HardSystem::vertices);
+ if (v.empty()) {
+ insert(v, S(r));
+ break;
+ }
+ }
+ }
+ } else {
+
+ double pDel = HardSystem::N / (z * HardSystem::maxOccupation());
+
+ if (pDel > r.uniform(0.0, 1.0)) {
+ tryDeletion(r.pick(HardSystem::vertices));
+ }
+ }
+
+ tryRandomMove(r);
+ }
+ }
+
+ void sweepLocal(Rng& r) {
+ for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
+ tryRandomMove(r);
+ }
+ }
+
+ void sweepSwap(Rng& r) {
+ for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
+ tryRandomSwap(r);
+ }
+ }
+
+ unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, Rng& r, bool dry = false) {
+ std::queue<std::reference_wrapper<Vertex<D, S>>> q;
+ q.push(v0);
+
+ unsigned n = 0;
+
+ while (!q.empty()) {
+ Vertex<D, S>& v = q.front();
+ q.pop();
+
+ if (!v.marked) {
+ Vector<D> xNew = R.apply(v.position);
+ Vertex<D, S>& vNew = HardSystem::vertices[HardSystem::vectorToIndex(xNew)];
+
+ v.marked = true;
+ vNew.marked = true;
+
+ S s = R.apply(v.state);
+ S sNew = R.apply(vNew.state);
+
+ std::list<std::reference_wrapper<Vertex<D, S>>> overlaps1 = overlaps(vNew, s, true);
+ std::list<std::reference_wrapper<Vertex<D, S>>> overlaps2 = overlaps(v, sNew, true);
+ overlaps1.splice(overlaps1.begin(), overlaps2);
+
+ for (Vertex<D, S>& 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, S>& v : HardSystem::vertices) {
+ if (!v.marked) {
+ 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);
+ }
+ }
+ }
+
+ for (Vertex<D, S>& v : HardSystem::vertices) {
+ v.marked = false;
+ }
+ }
+
+ int overlap(const System<D, S>& s) const {
+ 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);
+ o += HardSystem::vertices[i].state.dot(s2);
+ }
+
+ return o;
+ }
+};