summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ciamarra.cpp79
-rw-r--r--glass.hpp (renamed from glass.cpp)451
2 files changed, 279 insertions, 251 deletions
diff --git a/ciamarra.cpp b/ciamarra.cpp
new file mode 100644
index 0000000..8449cbb
--- /dev/null
+++ b/ciamarra.cpp
@@ -0,0 +1,79 @@
+#include <iostream>
+
+#include "glass.hpp"
+
+
+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) {
+ std::cerr << "◀";
+ } else if (v.state(0) == 0 && v.state(1) == -1) {
+ std::cerr << "▲";
+ } else if (v.state(0) == 0 && v.state(1) == 1) {
+ std::cerr << "▼";
+ } else if (v.state(0) == 0 && v.state(1) == 0) {
+ std::cerr << " ";
+ } else {
+ std::cerr << "X";
+ }
+
+ if (v.position(0) == s.L - 1) {
+ std::cerr << std::endl;
+ }
+ }
+}
+
+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;
+ unsigned Nmin = 2e2;
+ unsigned Nmax = 2e5;
+ double Tmin = 0.04;
+ double Tmax = 0.2;
+ double δT = 0.02;
+
+ CiamarraSystem<D> s(L);
+
+ Rng r;
+
+ double z = exp(1 / 0.08);
+
+ while (s.density() < 0.818) {
+ s.sweepGrandCanonical(z, r);
+ }
+
+ if (!s.compatible()) {
+ std::cerr << "Net compatible!" << std::endl;
+ return 1;
+ }
+
+ std::cerr << "Found state with appropriate density." << std::endl;
+
+ CiamarraSystem<D> s0 = s;
+
+ std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+
+ unsigned n = 1;
+ for (unsigned i = 0; i < 1e5; i++) {
+ if (n < 20 * log(i + 1)) {
+ n++;
+ std::cout << i << " " << s.overlap(s0) << std::endl;
+ }
+ s.sweepLocal(r);
+// s.sweepSwap(r);
+// s.swendsenWang(Transformation<D>(L, ms, r), r);
+ }
+
+ return 0;
+}
+
diff --git a/glass.cpp b/glass.hpp
index f2f44c1..da06542 100644
--- a/glass.cpp
+++ b/glass.hpp
@@ -1,9 +1,8 @@
-#include <eigen3/Eigen/Dense>
-#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>
-#include <iostream>
#include <list>
-#include <vector>
#include <queue>
+#include <vector>
+
+#include <eigen3/Eigen/Dense>
#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"
@@ -26,9 +25,7 @@ int iPow(int x, unsigned p) {
return x * tmp * tmp;
}
-unsigned mod(signed a, unsigned b) {
- return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b;
-}
+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;
@@ -51,7 +48,7 @@ template <int D> void one_sequences(std::list<std::array<double, D>>& sequences,
}
}
-template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
+template <int D> std::vector<Matrix<D>> generateTorusMatrices() {
std::vector<Matrix<D>> mats;
std::array<double, D> ini_sequence;
@@ -105,22 +102,22 @@ template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
return mats;
}
-template <unsigned D> class State : public Vector<D> {
+template <int D> class CiamarraState : public Vector<D> {
public:
- State() : Vector<D>(Vector<D>::Zero()) {}
+ CiamarraState() : Vector<D>(Vector<D>::Zero()) {}
- State(const Vector<D>& v) : Vector<D>(v) {}
+ CiamarraState(const Vector<D>& v) : Vector<D>(v) {}
- State(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; }
+ CiamarraState(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})) {}
+ CiamarraState(Rng& r) : CiamarraState(r.uniform(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);
}
@@ -128,74 +125,92 @@ public:
}
};
-template <unsigned D> class Transformation {
- public:
- unsigned L;
- Matrix<D> m;
- Vector<D> v;
+class BiroliState {
+public:
+ unsigned type;
+ unsigned occupiedNeighbors;
- Transformation(unsigned L) : L(L) {
- m.setIdentity();
- v.setZero();
- }
+ BiroliState() {
+ type = 0;
+ occupiedNeighbors = 0;
+ }
- Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {}
+ BiroliState(unsigned t, unsigned o) {
+ type = t;
+ occupiedNeighbors = o;
+ }
- 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);
- }
+ bool isEmpty() const { return type == 0; }
- v = v - m * v;
- }
+ void remove() { type = 0; };
+};
- Vector<D> apply(const Vector<D>& x) const {
- return mod<D>(v + m * x, L);
- }
+template <int D> class Transformation {
+public:
+ unsigned L;
+ Matrix<D> m;
+ Vector<D> v;
- Transformation<D> apply(const Transformation<D>& t) const {
- Transformation<D> tNew(L);
+ Transformation(unsigned L) : L(L) {
+ m.setIdentity();
+ v.setZero();
+ }
- tNew.m = m * t.m;
- tNew.v = apply(t.v);
+ Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {}
- return tNew;
+ 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);
}
- State<D> apply(const State<D>& x) const {
- return State<D>(m * x);
- }
+ v = v - m * v;
+ }
- Transformation<D> inverse() const {
- return Transformation<D>(L, m.transpose(), -m.transpose() * 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 <unsigned D> class HalfEdge;
+template <int D, typename State> class HalfEdge;
-template <unsigned D> class Vertex {
+template <int D, typename State> class Vertex {
public:
Vector<D> position;
- State<D> state;
- std::vector<HalfEdge<D>> adjacentEdges;
+ std::vector<HalfEdge<D, State>> adjacentEdges;
+ State state;
bool marked;
bool isEmpty() const { return state.isEmpty(); }
};
-template <unsigned D> class HalfEdge {
+template <int D, class State> class HalfEdge {
public:
- Vertex<D>& neighbor;
+ Vertex<D, State>& neighbor;
Vector<D> Δx;
- HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
+ HalfEdge(Vertex<D, State>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
};
-template <unsigned D> class System {
+template <int D, class State> class System {
public:
const unsigned L;
unsigned N;
- std::vector<Vertex<D>> vertices;
+ std::vector<Vertex<D, State>> vertices;
Transformation<D> orientation;
unsigned vectorToIndex(const Vector<D>& x) const {
@@ -226,126 +241,32 @@ public:
Δ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));
+ vertices[i].adjacentEdges.push_back(HalfEdge<D, State>(vertices[j], Δx));
+ vertices[j].adjacentEdges.push_back(HalfEdge<D, State>(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;
+ virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&,
+ const State&, bool = false) { return {}; };
- if (s.isEmpty()) {
- return o;
- }
+ virtual bool insert(Vertex<D, State>& v, const State& s) { return false; };
- if (!v.isEmpty() && !excludeSelf) {
- o.push_back(v);
- }
+ virtual bool remove(Vertex<D, State>& v) { return false;};
- 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);
- }
- }
- }
+ virtual bool tryRandomMove(Rng& r) { return false;};
- 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;
- }
- }
+ virtual bool swap(Vertex<D, State>& v1, Vertex<D, State>& v2) { 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;
+ Vertex<D, State>& v1 = r.pick(vertices);
+ Vertex<D, State>& v2 = r.pick(vertices);
- 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++;
- }
- }
+ return swap(v1, v2);
}
bool compatible() {
- for (Vertex<D>& v : vertices) {
+ for (Vertex<D, State>& v : vertices) {
if (overlaps(v, v.state, true).size() > 0) {
return false;
}
@@ -356,31 +277,28 @@ public:
double density() const { return N / pow(L, D); }
- unsigned maxOccupation() const {
- // return (2 * D * iPow(L, D)) / (2 * D + 1);
- return iPow(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 = maxOccupation() * z / (N + 1);
+ double pIns = size() * z / (N + 1);
if (pIns > r.uniform(0.0, 1.0)) {
while (true) {
- Vertex<D>& v = r.pick(vertices);
+ Vertex<D, State>& v = r.pick(vertices);
if (v.isEmpty()) {
- insert(v, State<D>(r));
+ insert(v, State(r));
break;
}
}
}
} else {
- double pDel = N / (z * maxOccupation());
+ double pDel = N / (z * size());
if (pDel > r.uniform(0.0, 1.0)) {
- tryDeletion(r.pick(vertices));
+ remove(r.pick(vertices));
}
}
@@ -400,39 +318,38 @@ public:
}
}
- unsigned flipCluster(const Transformation<D>& R, Vertex<D>& v0, Rng& r, bool dry = false) {
- std::queue<std::reference_wrapper<Vertex<D>>> q;
+ unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, Rng& r, bool dry = false) {
+ std::queue<std::reference_wrapper<Vertex<D, State>>> q;
q.push(v0);
unsigned n = 0;
while (!q.empty()) {
- Vertex<D>& v = q.front();
+ Vertex<D, State>& v = q.front();
q.pop();
if (!v.marked) {
Vector<D> xNew = R.apply(v.position);
- Vertex<D>& vNew = vertices[vectorToIndex(xNew)];
+ Vertex<D, State>& vNew = vertices[vectorToIndex(xNew)];
v.marked = true;
vNew.marked = true;
- State<D> s = R.apply(v.state);
- State<D> sNew = R.apply(vNew.state);
+ CiamarraState<D> s = R.apply(v.state);
+ CiamarraState<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);
+ std::list<std::reference_wrapper<Vertex<D, State>>> overlaps1 = overlaps(vNew, s, true);
+ std::list<std::reference_wrapper<Vertex<D, State>>> overlaps2 = overlaps(v, sNew, true);
overlaps1.splice(overlaps1.begin(), overlaps2);
- for (Vertex<D>& vn : overlaps1) {
+ for (Vertex<D, State>& vn : overlaps1) {
if (!vn.marked) {
q.push(vn);
}
}
if (!dry) {
- v.state = sNew;
- vNew.state = s;
+ swap(v, vNew);
}
n += 1;
@@ -443,7 +360,7 @@ public:
}
void swendsenWang(const Transformation<D>& R, Rng& r) {
- for (Vertex<D>& v : vertices) {
+ for (Vertex<D, State>& v : vertices) {
if (!v.marked) {
bool dry = 0.5 < r.uniform(0.0, 1.0);
unsigned n = flipCluster(R, v, r, dry);
@@ -453,101 +370,133 @@ public:
}
}
- for (Vertex<D>& v : vertices) {
+ for (Vertex<D, State>& v : vertices) {
v.marked = false;
}
}
- int overlap(const System<D>& s) const {
- int o = 0;
+ 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;
+ }
- 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);
+ 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;
}
-};
-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) {
- std::cerr << "◀";
- } else if (v.state(0) == 0 && v.state(1) == -1) {
- std::cerr << "▲";
- } else if (v.state(0) == 0 && v.state(1) == 1) {
- std::cerr << "▼";
- } else if (v.state(0) == 0 && v.state(1) == 0) {
- std::cerr << " ";
+ bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s) override {
+ if (overlaps(v, s).empty()) {
+ v.state = s;
+ this->N++;
+ return true;
} else {
- std::cerr << "X";
- }
-
- if (v.position(0) == s.L - 1) {
- std::cerr << std::endl;
+ return false;
}
}
-}
-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);
+ bool remove(Vertex<D, CiamarraState<D>>& v) override {
+ if (v.isEmpty()) {
+ return false;
+ } else {
+ v.state.remove();
+ this->N--;
+ return true;
+ }
}
- return x;
-}
-
-int main() {
- const unsigned D = 3;
- unsigned L = 15;
- unsigned Nmin = 2e2;
- unsigned Nmax = 2e5;
- double Tmin = 0.04;
- double Tmax = 0.2;
- double δT = 0.02;
-
- System<D> s(L);
- Rng r;
+ bool tryRandomMove(Rng& r) override {
+ Vertex<D, CiamarraState<D>>& v = r.pick(this->vertices);
+ CiamarraState<D> oldState = v.state;
- double z = exp(1 / 0.08);
+ if (!remove(v)) {
+ return false;
+ }
- while (s.density() < 0.818) {
- s.sweepGrandCanonical(z, r);
+ 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;
}
- if (!s.compatible()) {
- std::cerr << "Net compatible!" << std::endl;
- return 1;
- }
+ bool swap(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) {
+ std::swap(v1.state, v2.state);
+ return true;
+ } else {
+ return false;
+ }
+ }
- std::cerr << "Found state with appropriate density." << std::endl;
+ void setGroundState() {
+ this->N = 0;
- System<D> s0 = s;
+ 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;
- std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+ v.state.setZero() = Vector<D>::Zero();
- unsigned n = 1;
- for (unsigned i = 0; i < 1e5; i++) {
- if (n < 20 * log(i + 1)) {
- 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;
+ if (0 < a && a <= D) {
+ v.state(a - 1) = -1;
+ this->N++;
+ } else if (D < a) {
+ v.state(2 * D - a) = 1;
+ this->N++;
+ }
}
-// s.sweepLocal(r);
-// s.sweepSwap(r);
-// s.swendsenWang(Transformation<D>(L, ms, r), r);
}
- return 0;
-}
+ 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;
+ }
+};