diff options
| -rw-r--r-- | ciamarra.cpp | 79 | ||||
| -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; +} + @@ -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; + } +}; |