diff options
Diffstat (limited to 'biroli-mezard.cpp')
-rw-r--r-- | biroli-mezard.cpp | 178 |
1 files changed, 18 insertions, 160 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++) { |