From fadb6926d8e07ecacc9e0d1031a567494db4730a Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 22 Mar 2021 16:37:13 +0100 Subject: Split up functions, tired to generize in preparation for implementing Biroli-Mezard. --- ciamarra.cpp | 79 +++++++++ glass.cpp | 553 ----------------------------------------------------------- glass.hpp | 502 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 581 insertions(+), 553 deletions(-) create mode 100644 ciamarra.cpp delete mode 100644 glass.cpp create mode 100644 glass.hpp 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 + +#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 Vector randomVector(unsigned L, Rng& r) { + Vector 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 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 s0 = s; + + std::vector> ms = generateTorusMatrices(); + + 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(L, ms, r), r); + } + + return 0; +} + diff --git a/glass.cpp b/glass.cpp deleted file mode 100644 index f2f44c1..0000000 --- a/glass.cpp +++ /dev/null @@ -1,553 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include "pcg-cpp/include/pcg_random.hpp" -#include "randutils/randutils.hpp" - -using Rng = randutils::random_generator; - -template using Vector = Eigen::Matrix; -template using Matrix = Eigen::Matrix; - -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 Vector mod(const Eigen::MatrixBase& v, unsigned b) { - Vector u; - for (unsigned i = 0; i < D; i++) { - u(i) = mod(v(i), b); - } - return u; -} - -template void one_sequences(std::list>& sequences, unsigned level) { - if (level > 0) { - unsigned new_level = level - 1; - unsigned old_length = sequences.size(); - for (std::array& sequence : sequences) { - std::array new_sequence = sequence; - new_sequence[new_level] = -1; - sequences.push_front(new_sequence); - } - one_sequences(sequences, new_level); - } -} - -template std::vector> generateTorusMatrices() { - std::vector> mats; - - std::array ini_sequence; - ini_sequence.fill(1); - std::list> sequences; - sequences.push_back(ini_sequence); - - one_sequences(sequences, D); - - sequences.pop_back(); // don't want the identity matrix! - - for (std::array sequence : sequences) { - Matrix 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 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 class State : public Vector { -public: - State() : Vector(Vector::Zero()) {} - - State(const Vector& v) : Vector(v) {} - - State(unsigned a, signed b) : Vector(Vector::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 flip() const { - State s; - for (unsigned i = 0; i < D; i++) { - s(i) = -this->operator()(i); - } - return s; - } -}; - -template class Transformation { - public: - unsigned L; - Matrix m; - Vector v; - - Transformation(unsigned L) : L(L) { - m.setIdentity(); - v.setZero(); - } - - Transformation(unsigned L, const Matrix& m, const Vector& v) : L(L), m(m), v(v) {} - - Transformation(unsigned L, const std::vector>& 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 apply(const Vector& x) const { - return mod(v + m * x, L); - } - - Transformation apply(const Transformation& t) const { - Transformation tNew(L); - - tNew.m = m * t.m; - tNew.v = apply(t.v); - - return tNew; - } - - State apply(const State& x) const { - return State(m * x); - } - - Transformation inverse() const { - return Transformation(L, m.transpose(), -m.transpose() * v); - } -}; - -template class HalfEdge; - -template class Vertex { -public: - Vector position; - State state; - std::vector> adjacentEdges; - bool marked; - - bool isEmpty() const { return state.isEmpty(); } -}; - -template class HalfEdge { -public: - Vertex& neighbor; - Vector Δx; - - HalfEdge(Vertex& n, const Vector& d) : neighbor(n), Δx(d) {} -}; - -template class System { -public: - const unsigned L; - unsigned N; - std::vector> vertices; - Transformation orientation; - - unsigned vectorToIndex(const Vector& x) const { - unsigned i = 0; - for (unsigned d = 0; d < D; d++) { - i += x[d] * iPow(L, d); - } - return i; - } - - Vector indexToVector(unsigned i) const { - Vector x; - for (unsigned d = 0; d < D; d++) { - x[d] = (i / iPow(L, d)) % L; - } - return x; - } - - 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; - } - - for (unsigned d = 0; d < D; d++) { - Vector Δx = Vector::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(vertices[j], Δx)); - vertices[j].adjacentEdges.push_back(HalfEdge(vertices[i], -Δx)); - } - } - } - - std::list>> overlaps(Vertex& v, const State& s, - bool excludeSelf = false) { - std::list>> o; - - if (s.isEmpty()) { - return o; - } - - if (!v.isEmpty() && !excludeSelf) { - o.push_back(v); - } - - for (const HalfEdge& 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; - } - - bool insert(Vertex& v, const State& s) { - if (overlaps(v, s).empty()) { - v.state = s; - N++; - return true; - } else { - return false; - } - } - - bool tryDeletion(Vertex& v) { - if (v.isEmpty()) { - return false; - } else { - v.state.remove(); - N--; - return true; - } - } - - bool tryRandomMove(Rng& r) { - Vertex& v = r.pick(vertices); - State oldState = v.state; - - if (!tryDeletion(v)) { - return false; - } - - if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) { - for (HalfEdge& e : v.adjacentEdges) { - if (1 == e.Δx.dot(oldState)) { - if (insert(e.neighbor, oldState.flip())) { - return true; - } - break; - } - } - } else { - State newState(r); - while (newState.dot(oldState) == 1) { - newState = State(r); - } - if (insert(v, newState)) { - return true; - } - } - v.state = oldState; - N++; - return false; - } - - bool trySwap(Vertex& v1, Vertex& 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& v1 = r.pick(vertices); - Vertex& v2 = r.pick(vertices); - - return trySwap(v1, v2); - } - - void setGroundState() { - N = 0; - - for (Vertex& v : vertices) { - unsigned a = 0; - for (unsigned d = 0; d < D; d++) { - a += (d + 1) * v.position(d); - } - a %= 2 * D + 1; - - v.state.setZero() = Vector::Zero(); - - if (0 < a && a <= D) { - v.state(a - 1) = -1; - N++; - } else if (D < a) { - v.state(2 * D - a) = 1; - N++; - } - } - } - - bool compatible() { - for (Vertex& 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& v = r.pick(vertices); - if (v.isEmpty()) { - insert(v, State(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& R, Vertex& v0, Rng& r, bool dry = false) { - std::queue>> q; - q.push(v0); - - unsigned n = 0; - - while (!q.empty()) { - Vertex& v = q.front(); - q.pop(); - - if (!v.marked) { - Vector xNew = R.apply(v.position); - Vertex& vNew = vertices[vectorToIndex(xNew)]; - - v.marked = true; - vNew.marked = true; - - State s = R.apply(v.state); - State sNew = R.apply(vNew.state); - - std::list>> overlaps1 = overlaps(vNew, s, true); - std::list>> overlaps2 = overlaps(v, sNew, true); - overlaps1.splice(overlaps1.begin(), overlaps2); - - for (Vertex& vn : overlaps1) { - if (!vn.marked) { - q.push(vn); - } - } - - if (!dry) { - v.state = sNew; - vNew.state = s; - } - - n += 1; - } - } - - return n; - } - - void swendsenWang(const Transformation& R, Rng& r) { - for (Vertex& 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& v : vertices) { - v.marked = false; - } - } - - int overlap(const System& s) const { - int o = 0; - - for (unsigned i = 0; i < vertices.size(); i++) { - State 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) { - 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 Vector randomVector(unsigned L, Rng& r) { - Vector 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; - - System 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; - - System s0 = s; - - std::vector> ms = generateTorusMatrices(); - - 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 m = r.pick(ms); - Vertex& v = r.pick(s.vertices); - unsigned nC = s.flipCluster(Transformation(L, m, v.position - m * v.position), v, r); - std::cerr << nC << std::endl; - for (Vertex& v : s.vertices) { - v.marked = false; - } -// s.sweepLocal(r); -// s.sweepSwap(r); -// s.swendsenWang(Transformation(L, ms, r), r); - } - - return 0; -} - diff --git a/glass.hpp b/glass.hpp new file mode 100644 index 0000000..da06542 --- /dev/null +++ b/glass.hpp @@ -0,0 +1,502 @@ +#include +#include +#include + +#include + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" + +using Rng = randutils::random_generator; + +template using Vector = Eigen::Matrix; +template using Matrix = Eigen::Matrix; + +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 Vector mod(const Eigen::MatrixBase& v, unsigned b) { + Vector u; + for (unsigned i = 0; i < D; i++) { + u(i) = mod(v(i), b); + } + return u; +} + +template void one_sequences(std::list>& sequences, unsigned level) { + if (level > 0) { + unsigned new_level = level - 1; + unsigned old_length = sequences.size(); + for (std::array& sequence : sequences) { + std::array new_sequence = sequence; + new_sequence[new_level] = -1; + sequences.push_front(new_sequence); + } + one_sequences(sequences, new_level); + } +} + +template std::vector> generateTorusMatrices() { + std::vector> mats; + + std::array ini_sequence; + ini_sequence.fill(1); + std::list> sequences; + sequences.push_back(ini_sequence); + + one_sequences(sequences, D); + + sequences.pop_back(); // don't want the identity matrix! + + for (std::array sequence : sequences) { + Matrix 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 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 class CiamarraState : public Vector { +public: + CiamarraState() : Vector(Vector::Zero()) {} + + CiamarraState(const Vector& v) : Vector(v) {} + + CiamarraState(unsigned a, signed b) : Vector(Vector::Zero()) { this->operator()(a) = b; } + + CiamarraState(Rng& r) : CiamarraState(r.uniform(0, D - 1), r.pick({-1, 1})) {} + + bool isEmpty() const { return this->squaredNorm() == 0; } + + void remove() { this->setZero(); } + + CiamarraState flip() const { + CiamarraState s; + for (unsigned i = 0; i < D; i++) { + s(i) = -this->operator()(i); + } + return s; + } +}; + +class BiroliState { +public: + unsigned type; + unsigned occupiedNeighbors; + + BiroliState() { + type = 0; + occupiedNeighbors = 0; + } + + BiroliState(unsigned t, unsigned o) { + type = t; + occupiedNeighbors = o; + } + + bool isEmpty() const { return type == 0; } + + void remove() { type = 0; }; +}; + +template class Transformation { +public: + unsigned L; + Matrix m; + Vector v; + + Transformation(unsigned L) : L(L) { + m.setIdentity(); + v.setZero(); + } + + Transformation(unsigned L, const Matrix& m, const Vector& v) : L(L), m(m), v(v) {} + + Transformation(unsigned L, const std::vector>& 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; + } + + Transformation inverse() const { + return Transformation(L, m.transpose(), -m.transpose() * v); + } + + Vector apply(const Vector& x) const { return mod(v + m * x, L); } + + Transformation apply(const Transformation& t) const { + Transformation tNew(L); + + tNew.m = m * t.m; + tNew.v = apply(t.v); + + return tNew; + } + + CiamarraState apply(const CiamarraState& s) const { return CiamarraState(m * s); } + + BiroliState apply(const BiroliState& s) const { return s; } +}; + +template class HalfEdge; + +template class Vertex { +public: + Vector position; + std::vector> adjacentEdges; + State state; + bool marked; + + bool isEmpty() const { return state.isEmpty(); } +}; + +template class HalfEdge { +public: + Vertex& neighbor; + Vector Δx; + + HalfEdge(Vertex& n, const Vector& d) : neighbor(n), Δx(d) {} +}; + +template class System { +public: + const unsigned L; + unsigned N; + std::vector> vertices; + Transformation orientation; + + unsigned vectorToIndex(const Vector& x) const { + unsigned i = 0; + for (unsigned d = 0; d < D; d++) { + i += x[d] * iPow(L, d); + } + return i; + } + + Vector indexToVector(unsigned i) const { + Vector x; + for (unsigned d = 0; d < D; d++) { + x[d] = (i / iPow(L, d)) % L; + } + return x; + } + + 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; + } + + for (unsigned d = 0; d < D; d++) { + Vector Δx = Vector::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(vertices[j], Δx)); + vertices[j].adjacentEdges.push_back(HalfEdge(vertices[i], -Δx)); + } + } + } + + virtual std::list>> overlaps(Vertex&, + const State&, bool = false) { return {}; }; + + virtual bool insert(Vertex& v, const State& s) { return false; }; + + virtual bool remove(Vertex& v) { return false;}; + + virtual bool tryRandomMove(Rng& r) { return false;}; + + virtual bool swap(Vertex& v1, Vertex& v2) { return false;}; + + bool tryRandomSwap(Rng& r) { + Vertex& v1 = r.pick(vertices); + Vertex& v2 = r.pick(vertices); + + return swap(v1, v2); + } + + bool compatible() { + for (Vertex& v : vertices) { + if (overlaps(v, v.state, true).size() > 0) { + return false; + } + } + + return true; + } + + double density() const { return N / pow(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 = size() * z / (N + 1); + + if (pIns > r.uniform(0.0, 1.0)) { + while (true) { + Vertex& v = r.pick(vertices); + if (v.isEmpty()) { + insert(v, State(r)); + break; + } + } + } + } else { + + double pDel = N / (z * size()); + + if (pDel > r.uniform(0.0, 1.0)) { + remove(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& R, Vertex& v0, Rng& r, bool dry = false) { + std::queue>> q; + q.push(v0); + + unsigned n = 0; + + while (!q.empty()) { + Vertex& v = q.front(); + q.pop(); + + if (!v.marked) { + Vector xNew = R.apply(v.position); + Vertex& vNew = vertices[vectorToIndex(xNew)]; + + v.marked = true; + vNew.marked = true; + + CiamarraState s = R.apply(v.state); + CiamarraState sNew = R.apply(vNew.state); + + std::list>> overlaps1 = overlaps(vNew, s, true); + std::list>> overlaps2 = overlaps(v, sNew, true); + overlaps1.splice(overlaps1.begin(), overlaps2); + + for (Vertex& vn : overlaps1) { + if (!vn.marked) { + q.push(vn); + } + } + + if (!dry) { + swap(v, vNew); + } + + n += 1; + } + } + + return n; + } + + void swendsenWang(const Transformation& R, Rng& r) { + for (Vertex& 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& v : vertices) { + v.marked = false; + } + } + + virtual int overlap(const System& s) const { return 0; }; +}; + +template class CiamarraSystem : public System> { + public: + + using System>::System; + std::list>>> + overlaps(Vertex>& v, const CiamarraState& s, + bool excludeSelf = false) override { + std::list>>> o; + + if (s.isEmpty()) { + return o; + } + + if (!v.isEmpty() && !excludeSelf) { + o.push_back(v); + } + + for (const HalfEdge>& 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; + } + + bool insert(Vertex>& v, const CiamarraState& s) override { + if (overlaps(v, s).empty()) { + v.state = s; + this->N++; + return true; + } else { + return false; + } + } + + bool remove(Vertex>& v) override { + if (v.isEmpty()) { + return false; + } else { + v.state.remove(); + this->N--; + return true; + } + } + + bool tryRandomMove(Rng& r) override { + Vertex>& v = r.pick(this->vertices); + CiamarraState oldState = v.state; + + if (!remove(v)) { + return false; + } + + if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) { + for (HalfEdge>& e : v.adjacentEdges) { + if (1 == e.Δx.dot(oldState)) { + if (insert(e.neighbor, oldState.flip())) { + return true; + } + break; + } + } + } else { + CiamarraState newState(r); + while (newState.dot(oldState) == 1) { + newState = CiamarraState(r); + } + if (insert(v, newState)) { + return true; + } + } + v.state = oldState; + this->N++; + return false; + } + + bool swap(Vertex>& v1, Vertex>& 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; + } + } + + void setGroundState() { + this->N = 0; + + for (Vertex>& v : this->vertices) { + unsigned a = 0; + for (unsigned d = 0; d < D; d++) { + a += (d + 1) * v.position(d); + } + a %= 2 * D + 1; + + v.state.setZero() = Vector::Zero(); + + if (0 < a && a <= D) { + v.state(a - 1) = -1; + this->N++; + } else if (D < a) { + v.state(2 * D - a) = 1; + this->N++; + } + } + } + + int overlap(const System>& s) const override { + int o = 0; + + for (unsigned i = 0; i < this->vertices.size(); i++) { + CiamarraState 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; + } +}; -- cgit v1.2.3-70-g09d2