diff options
| -rw-r--r-- | biroli-mezard.cpp | 178 | ||||
| -rw-r--r-- | ciamarra.cpp | 468 | ||||
| -rw-r--r-- | compile_flags.txt | 3 | ||||
| -rw-r--r-- | distinguishable.cpp | 6 | ||||
| -rw-r--r-- | glass.hpp | 209 |
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); } @@ -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; + } +}; |