From 1cf96564223e5f3550c5f140f50435d4e55dbe0a Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 24 Jun 2021 11:18:52 +0200 Subject: Started splitting up shared functions into another file. --- distinguishable.cpp | 258 ++++++++++------------------------------------------ glass.hpp | 167 ++++++++++++++++++++++++++++++++++ 2 files changed, 214 insertions(+), 211 deletions(-) create mode 100644 glass.hpp diff --git a/distinguishable.cpp b/distinguishable.cpp index d9f51ea..4b66989 100644 --- a/distinguishable.cpp +++ b/distinguishable.cpp @@ -1,192 +1,25 @@ -#include -#include #include -#include -#include #include -#include "pcg-cpp/include/pcg_random.hpp" -#include "randutils/randutils.hpp" +#include "glass.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; -} - -class State { +class DistinguishableState { public: unsigned type; - State() { - type = 0; - } - - State(unsigned t) { - type = t; - } + DistinguishableState() { type = 0; } - bool isEmpty() const { return type == 0; } + DistinguishableState(unsigned t) { type = t; } + bool empty() 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; - } - - 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; - } - - Transformation inverse() const { - return Transformation(L, m.transpose(), -m.transpose() * v); - } -}; - -template class HalfEdge; - -template class Vertex { -public: - Vector position; - Vector initialPosition; - 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; + std::vector> vertices; std::vector interaction; unsigned vectorToIndex(const Vector& x) const { @@ -218,8 +51,8 @@ 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(vertices[j], Δx)); - vertices[j].adjacentEdges.push_back(HalfEdge(vertices[i], -Δx)); + vertices[i].adjacentEdges.push_back(HalfEdge(vertices[j], Δx)); + vertices[j].adjacentEdges.push_back(HalfEdge(vertices[i], -Δx)); } } @@ -232,8 +65,8 @@ public: } } - double pairEnergy(const State& s1, const State& s2) const { - if (s1.isEmpty() || s2.isEmpty()) { + double pairEnergy(const DistinguishableState& s1, const DistinguishableState& s2) const { + if (s1.empty() || s2.empty()) { return 0; } else { if (s1.type < s2.type) { @@ -244,10 +77,10 @@ public: } } - double siteEnergy(const Vertex& v) const { + double siteEnergy(const Vertex& v) const { double E = 0; - for (const HalfEdge& e : v.adjacentEdges) { + for (const HalfEdge& e : v.adjacentEdges) { E += pairEnergy(v.state, e.neighbor.state); } @@ -257,14 +90,15 @@ public: double energy() const { double E = 0; - for (const Vertex& v : vertices) { + for (const Vertex& v : vertices) { E += siteEnergy(v); } return E / 2; } - bool trySwap(Vertex& v1, Vertex& v2, double T, Rng& r) { + bool trySwap(Vertex& v1, Vertex& v2, double T, + Rng& r) { double E0 = siteEnergy(v1) + siteEnergy(v2); std::swap(v1.state, v2.state); double E1 = siteEnergy(v1) + siteEnergy(v2); @@ -282,15 +116,15 @@ public: } bool tryRandomMove(double T, Rng& r) { - Vertex& v1 = r.pick(vertices); + Vertex& v1 = r.pick(vertices); - if (v1.isEmpty()) { + if (v1.empty()) { return false; } - Vertex& v2 = (r.pick(v1.adjacentEdges)).neighbor; + Vertex& v2 = (r.pick(v1.adjacentEdges)).neighbor; - if (!v2.isEmpty()) { + if (!v2.empty()) { return false; } @@ -298,8 +132,8 @@ public: } bool tryRandomSwap(double T, Rng& r) { - Vertex& v1 = r.pick(vertices); - Vertex& v2 = r.pick(vertices); + Vertex& v1 = r.pick(vertices); + Vertex& v2 = r.pick(vertices); if (v1.state.type != v2.state.type) { return trySwap(v1, v2, T, r); @@ -322,10 +156,11 @@ public: } } - unsigned flipCluster(const Transformation& R, Vertex& v0, double T, Rng& r, bool dry = false) { - std::queue>, 2>> q; + unsigned flipCluster(const Transformation& R, Vertex& v0, double T, + Rng& r, bool dry = false) { + std::queue>, 2>> q; Vector x0New = R.apply(v0.position); - Vertex& v0New = vertices[vectorToIndex(x0New)]; + Vertex& v0New = vertices[vectorToIndex(x0New)]; q.push({v0, v0New}); unsigned n = 0; @@ -334,24 +169,24 @@ public: auto [vR, vNewR] = q.front(); q.pop(); - Vertex& v = vR; - Vertex& vNew = vNewR; + Vertex& v = vR; + Vertex& vNew = vNewR; if (!v.marked && !vNew.marked) { v.marked = true; vNew.marked = true; - for (HalfEdge& e : v.adjacentEdges) { - Vertex& vn = e.neighbor; + for (HalfEdge& e : v.adjacentEdges) { + Vertex& vn = e.neighbor; Vector xnNew = R.apply(vn.position); - Vertex& vnNew = vertices[vectorToIndex(xnNew)]; + Vertex& vnNew = vertices[vectorToIndex(xnNew)]; if (!vn.marked && !vnNew.marked) { double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state); double E1 = pairEnergy(vNew.state, vn.state) + pairEnergy(v.state, vnNew.state); double ΔE = E1 - E0; - + if (exp(-ΔE / T) < r.uniform(0.0, 1.0)) { q.push({vn, vnNew}); } @@ -372,21 +207,21 @@ public: unsigned wolff(const Transformation& R, double T, Rng& r) { unsigned n = flipCluster(R, r.pick(vertices), T, r); - for (Vertex& v : vertices) { + for (Vertex& v : vertices) { v.marked = false; } return n; } void swendsenWang(const Transformation& R, double T, Rng& r) { - for (Vertex& v : vertices) { + for (Vertex& v : vertices) { if (!v.marked) { bool dry = 0.5 < r.uniform(0.0, 1.0); flipCluster(R, v, T, r, dry); } } - for (Vertex& v : vertices) { + for (Vertex& v : vertices) { v.marked = false; } } @@ -395,15 +230,16 @@ public: 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); + // DistinguishableState s2 = + // orientation.apply(s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].state); + // o += vertices[i].state.dot(s2); } return o; } void setInitialPosition() { - for (Vertex& v: vertices) { + for (Vertex& v : vertices) { v.initialPosition = v.position; } } @@ -413,8 +249,8 @@ public: Vector k1 = {M_PI, 0}; Vector k2 = {0, M_PI}; - for (const Vertex& v : vertices) { - if (!v.isEmpty()) { + for (const Vertex& v : vertices) { + if (!v.empty()) { F += cos(k1.dot(v.position - v.initialPosition)); F += cos(k2.dot(v.position - v.initialPosition)); } @@ -447,7 +283,7 @@ int main() { std::vector> ms = generateTorusMatrices(); for (unsigned j = 0; j < 10 * s.vertices.size(); j++) { - //unsigned nC = s.wolff(Transformation(L, m, v.position - m * v.position), T, r); + // unsigned nC = s.wolff(Transformation(L, m, v.position - m * v.position), T, r); s.tryRandomSwap(1000, r); } @@ -463,17 +299,17 @@ int main() { /* for (unsigned i = 0; i < 1e5; i++) { Matrix m = r.pick(ms); - Vertex& v = r.pick(s.vertices); + Vertex& v = r.pick(s.vertices); s.wolff(Transformation(L, m, v.position - m * v.position), T, r); } */ /* - */ + */ s.setInitialPosition(); std::cout << T << " "; for (unsigned i = 0; i < 1e4; i++) { Matrix m = r.pick(ms); - Vertex& v = r.pick(s.vertices); + Vertex& v = r.pick(s.vertices); s.wolff(Transformation(L, m, v.position - m * v.position), T, r); std::cout << s.selfIntermediateScattering() << " "; } @@ -487,9 +323,9 @@ int main() { */ std::cout << std::endl; std::cerr << T << " " << s.energy() / N << std::endl; -// s.sweepLocal(r); -// s.sweepSwap(r); -// s.swendsenWang(Transformation(L, ms, r), r); + // 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..fbff7d6 --- /dev/null +++ b/glass.hpp @@ -0,0 +1,167 @@ +#include + +#include +#include + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" + +template using Vector = Eigen::Matrix; +template using Matrix = Eigen::Matrix; + +using Rng = randutils::random_generator; + +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 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; + } + + Transformation inverse() const { + return Transformation(L, m.transpose(), -m.transpose() * v); + } +}; + +template concept State = requires(T v) { + { v.empty() } -> std::same_as; + {v.remove()}; +}; + +template class HalfEdge; + +template class Vertex { +public: + Vector position; + Vector initialPosition; + S state; + std::vector> adjacentEdges; + bool marked; + + bool empty() const { return state.empty(); } +}; + +template class HalfEdge { +public: + Vertex& neighbor; + Vector Δx; + + HalfEdge(Vertex& n, const Vector& d) : neighbor(n), Δx(d) {} +}; + -- cgit v1.2.3-70-g09d2