summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--glass.cpp570
1 files changed, 368 insertions, 202 deletions
diff --git a/glass.cpp b/glass.cpp
index e8739d0..ba983d2 100644
--- a/glass.cpp
+++ b/glass.cpp
@@ -1,294 +1,447 @@
+#include <eigen3/Eigen/Dense>
#include <iostream>
-#include <limits>
-#include <vector>
#include <list>
-#include <set>
-#include <eigen3/Eigen/Dense>
-#include <unordered_map>
+#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>;
+template <int D> using Vector = Eigen::Matrix<int, D, 1>;
+template <int D> using Matrix = Eigen::Matrix<int, D, D>;
-// Integer power function
int iPow(int x, unsigned p) {
- if (p == 0) return 1;
- if (p == 1) return x;
+ 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;
+ if (p % 2 == 0)
+ return tmp * tmp;
+ else
+ return x * tmp * tmp;
}
unsigned mod(signed a, unsigned b) {
- return ((a %= b) < 0) ? a + b : a;
+ return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b;
}
-template<unsigned D> class State {
- public:
- Vector<D> σ;
+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;
+}
- State(unsigned a, signed b) : σ(Vector<D>::Zero()) {
- σ[a] = b;
+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);
+ }
+}
- State(Rng& r) : State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {}
+template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
+ std::vector<Matrix<D>> mats;
- State<D> apply(const Matrix<D>& m) const {
- return {m * σ};
+ 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;
+ }
+ }
}
-};
-template <unsigned D> class Vertex;
+ mats.push_back(m);
+ }
-template<unsigned D> class Particle {
- public:
- Vertex<D>& position;
- State<D> state;
- bool marked;
+ 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);
+ }
+ }
+ }
- Particle(Vertex<D>& p, const State<D>& s) : position(p), state(s) {}
-};
+ return mats;
+}
-template<unsigned D> class HalfEdge {
- public:
- Vertex<D>& neighbor;
- Vector<D> Δx;
+template <unsigned D> class State : public Vector<D> {
+public:
+ State() : Vector<D>(Vector<D>::Zero()) {}
- HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
-};
+ State(const Vector<D>& v) : Vector<D>(v) {}
-template<unsigned D> class Vertex {
- public:
- Vector<D> position;
- std::vector<HalfEdge<D>> adjacentEdges;
- typename std::unordered_map<unsigned, Particle<D>>::iterator occupant;
+ 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(); }
};
-template<unsigned D> class System {
+template <unsigned D> class Transformation {
public:
- const unsigned L;
- std::vector<Vertex<D>> vertices;
- std::unordered_map<unsigned, Particle<D>> particles;
+ unsigned L;
+ Matrix<D> m;
+ Vector<D> v;
- unsigned vectorToIndex(const Vector<D>& x) const {
- unsigned i = 0;
- for (unsigned d = 0; d < D; d++) {
- i += x[d] * iPow(L, d);
+ 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);
}
- 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;
+ Vector<D> apply(const Vector<D>& x) const {
+ return mod<D>(v + m * (x - v), L);
}
- bool isEmpty(const Vertex<D>& v) const {
- return v.occupant == particles.end();
+ State<D> apply(const State<D>& x) const {
+ return State<D>(m * x);
}
+};
+
+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;
- unsigned N() const {
- return particles.size();
+ 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;
+ }
- System(unsigned L) : L(L), vertices(iPow(L, D)) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
- vertices[i].position = indexToVector(i);
- vertices[i].adjacentEdges.reserve(2 * D);
- vertices[i].occupant = particles.end();
- }
+ 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;
+ }
- 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));
- }
+ System(unsigned L) : L(L), N(0), vertices(iPow(L, D)) {
+ 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<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<typename std::unordered_map<unsigned, Particle<D>>::iterator> overlaps(const Particle<D>& p, bool excludeSelf = false) const {
- std::list<typename std::unordered_map<unsigned, Particle<D>>::iterator> o;
+ 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;
- if (!excludeSelf && !isEmpty(p.position)) {
- o.push_back(p.position.occupant);
- }
+ if (s.isEmpty()) {
+ return o;
+ }
+
+ if (!v.isEmpty() && !excludeSelf) {
+ o.push_back(v);
+ }
- for (const HalfEdge<D>& e : p.position.adjacentEdges) {
- if (!isEmpty(e.neighbor)) {
- if (p.state.σ.dot(e.Δx) == 1 || e.neighbor.occupant->second.state.σ.dot(e.Δx) == -1) {
- o.push_back(e.neighbor.occupant);
- }
+ 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);
}
}
+ }
- return o;
+ 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 insert(const Particle<D>& p, Rng& r) {
- if (overlaps(p).empty()) {
- typename std::unordered_map<unsigned, Particle<D>>::iterator it;
- std::tie(it, std::ignore) = particles.insert({r.uniform((unsigned)0, (unsigned)pow(2, 28)), p});
- p.position.occupant = it;
+ bool tryRandomMove(Rng& r) {
+ Vertex<D>& v = r.pick(vertices);
+ State<D> oldState = v.state;
+ if (!tryDeletion(v)) {
+ return false;
+ }
+
+ if (0.2 > r.uniform(0.0, 1.0)) {
+ if (insert(v, State<D>(r))) {
+ return true;
+ }
+ } else {
+ if (insert(r.pick(v.adjacentEdges).neighbor, State<D>(r))) {
return true;
- } else {
- return false;
}
}
+ v.state = oldState;
+ N++;
+ return false;
+ }
- void remove(typename std::unordered_map<unsigned, Particle<D>>::iterator ip) {
- ip->second.position.occupant = particles.end();
- particles.erase(ip);
+ 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 randomMove(Rng& r) {
- Particle<D>& p = *(r.pick(particles));
- Particle<D> pN = p;
+ bool tryRandomSwap(Rng& r) {
+ Vertex<D>& v1 = r.pick(vertices);
+ Vertex<D>& v2 = r.pick(vertices);
- if (0.2 < r.uniform(0.0, 1.0)) {
- pN.position = r.pick(p.position.adjacentEdges).neighbor;
+ 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);
}
- pN.state = randomState<D>(r);
+ a %= 2 * D + 1;
- if (overlaps(pN, true).empty()) {
- remove(&p);
- insert(pN);
- return true;
- } else {
- return false;
+ 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++;
}
}
- */
+ }
- bool swap(Particle<D>& p1, Particle<D>& p2) {
- std::swap(p1.state, p2.state);
- if (overlaps(p1, true).empty() && overlaps(p2, true).empty()) {
- return true;
- } else {
- std::swap(p1.state, p2.state);
+ bool compatible() {
+ for (Vertex<D>& v : vertices) {
+ if (overlaps(v, v.state, true).size() > 0) {
return false;
}
}
- bool randomSwap(Rng& r) {
- return swap(r.pick(particles).second, r.pick(particles).second);
- }
+ return true;
+ }
- void setGroundState(Rng& r) {
- for (typename std::unordered_map<unsigned, Particle<D>>::iterator ip = particles.begin(); ip != particles.end(); ip++) {
- remove(ip);
- }
+ double density() const { return N / pow(L, D); }
- 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;
-
- if (a > 0) {
- if (a <= D) {
- State<D> s(a - 1, -1);
- insert(Particle<D>(v, s), r);
- } else if (D < a) {
- State<D> s(2 * D - a, 1);
- insert(Particle<D>(v, s), r);
+ 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());
- bool compatible() {
- for (Particle<D> p : particles) {
- if (!overlaps(p, true).empty()) {
- return false;
+ if (pDel > r.uniform(0.0, 1.0)) {
+ tryDeletion(r.pick(vertices));
}
}
- return true;
+ tryRandomMove(r);
+ tryRandomSwap(r);
}
+ }
- double density() const {
- return N() / pow(L, D);
+ void sweepLocal(double z, Rng& r) {
+ for (unsigned i = 0; i < iPow(L, D); i++) {
+ tryRandomMove(r);
}
+ }
- unsigned maxOccupation() const {
-// return (2 * D * iPow(L, D)) / (2 * D + 1);
- return iPow(L, D);
+ void sweepSwap(double z, Rng& r) {
+ for (unsigned i = 0; i < iPow(L, D); i++) {
+ tryRandomSwap(r);
}
+ }
- 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 (isEmpty(v)) {
- insert(Particle<D>(v, State<D>(r)), r);
- break;
- }
- }
- }
- } else {
+ 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)];
- double pDel = N() / (z * maxOccupation());
+ v.marked = true;
+ vNew.marked = true;
- if (pDel > r.uniform(0.0, 1.0)) {
- remove(r.choose(particles));
+ 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);
}
}
-// randomMove(r);
- randomSwap(r);
- }
- }
+ if (!dry) {
+ v.state = sNew;
+ vNew.state = s;
+ }
- void sweepLocal(double z, Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
-// randomMove(r);
+ n += 1;
}
}
- void sweepSwap(double z, Rng& r) {
- for (unsigned i = 0; i < iPow(L, D); i++) {
- randomSwap(r);
+ return n;
+ }
+
+ void swendsenWang(const Transformation<D>& R, Rng& r) {
+ for (Vertex<D>& v : vertices) {
+ if (!v.marked) {
+ unsigned n = flipCluster(R, v, r, 0.5 < r.uniform(0.0, 1.0));
+ std::cerr << n << " ";
}
}
- void flipCluster(const Matrix<D>& R, Vertex<D>& v0, Rng& r, bool dry = false) {
- std::list<std::reference_wrapper<Vertex<D>>> q = {v0};
+ std::cerr << std::endl;
- while (!q.empty()) {
- Vertex<D>& v = q.front();
- q.pop_front();
+ for (Vertex<D>& v : vertices) {
+ v.marked = false;
+ }
+ }
+};
- Vector<D> xNew = R * v.position;
- State<D> sNew = v.state.apply(R);
+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";
+ }
- for (Vertex<D>& vn : overlaps(vertices[vectorToIndex(xNew)], sNew)) {
- if (!vn.marked) {
- vn.marked = true;
- q.push_back(vn);
- }
- }
- }
+ if (v.position(0) == s.L - 1) {
+ std::cerr << std::endl;
}
-};
+ }
+}
-template<unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
+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);
@@ -309,12 +462,32 @@ int main() {
Rng r;
+ double z = exp(1 / 0.15);
+
+ s.setGroundState();
+ std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+
+ for (unsigned n = 0; n < 1e3; n++) {
+ s.sweepGrandCanonical(z, r);
+ }
+
+ if (s.compatible()) {
+ std::cerr << "Still compatible!" << std::endl;
+ }
+
+
+ for (unsigned n = 0; n < 10; n++) {
+ s.swendsenWang(Transformation<D>(L, ms, r), r);
+ }
+
+ if (s.compatible()) {
+ std::cerr << "Still compatible!" << std::endl;
+ }
+
/*
+
for (unsigned N = Nmin; N <= Nmax; N *= 10) {
s.setGroundState();
- if (!s.compatible()) {
- return 1;
- }
for (double T = Tmin; T < Tmax; T *= 1 + δT) {
double z = exp(1 / T);
@@ -337,13 +510,6 @@ int main() {
}
*/
- s.setGroundState(r);
-
- for (unsigned n = 0; n < 10; n++) {
- s.sweepGrandCanonical(exp(1/0.15), r);
- }
-
return 0;
}
-