summaryrefslogtreecommitdiff
path: root/percolation.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'percolation.cpp')
-rw-r--r--percolation.cpp118
1 files changed, 55 insertions, 63 deletions
diff --git a/percolation.cpp b/percolation.cpp
index 08c6970..a5c2a36 100644
--- a/percolation.cpp
+++ b/percolation.cpp
@@ -14,128 +14,128 @@ using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
class Vertex;
-class Edge;
class HalfEdge {
public:
- Edge& parent;
+ unsigned index;
Vertex& neighbor;
- HalfEdge(Edge& e, Vertex& v) : parent(e), neighbor(v) {};
+ HalfEdge(unsigned index, Vertex& v) : index(index), neighbor(v) {};
};
class Vertex {
public:
std::vector<HalfEdge> bonds;
unsigned index;
- unsigned cluster;
};
-class Edge {
-public:
- bool active;
-};
+class Cluster : public std::list<std::reference_wrapper<const Vertex>> {};
-class Cluster : public std::list<std::reference_wrapper<Vertex>> {};
-class Clusters {
-public:
- std::vector<Cluster> lists;
- std::vector<unsigned> indices;
-};
+unsigned uPow(unsigned x, unsigned a) {
+ unsigned result = 1;
+ while (a) {
+ if (a & 1) {
+ result *= x;
+ }
+ a >>= 1;
+ x *= x;
+ }
+ return result;
+}
class Graph {
public:
unsigned D;
unsigned L;
std::vector<Vertex> vertices;
- std::vector<Edge> edges;
- std::vector<Cluster> clusters;
+ std::vector<unsigned> multiplicities;
unsigned Nv() const {
return vertices.size();
}
unsigned Ne() const {
- return edges.size();
+ return D * Nv();
}
unsigned squaredDistance(unsigned vi, unsigned vj) const {
unsigned sd = 0;
for (unsigned i = 0; i < D; i++) {
- int x1 = (unsigned)fmod((unsigned)(vi / (unsigned)pow(L, i)), pow(L, i + 1));
- int x2 = (unsigned)fmod((unsigned)(vj / (unsigned)pow(L, i)), pow(L, i + 1));
+ int x1 = (vi / uPow(L, i)) % uPow(L, i + 1);
+ int x2 = (vj / uPow(L, i)) % uPow(L, i + 1);
unsigned Δx = abs(x1 - x2);
if (Δx > L / 2) {
Δx = L - Δx;
}
- sd += pow(Δx, 2);
+ sd += Δx * Δx;
}
return sd;
}
/* Initialize a square lattice */
- Graph(unsigned D, unsigned L) : D(D), L(L), vertices(pow(L, D)), edges(D * pow(L, D)) {
+ Graph(unsigned D, unsigned L) : D(D), L(L), vertices(uPow(L, D)), multiplicities(D * L * L / 4) {
for (unsigned i = 0; i < Nv(); i++) {
vertices[i].index = i;
vertices[i].bonds.reserve(2 * D);
for (unsigned j = 0; j < D; j++) {
- unsigned n1 = (unsigned)pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) +
- (unsigned)fmod(i + pow(L, j), pow(L, j + 1));
- unsigned n2 = (unsigned)pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) +
- (unsigned)fmod(pow(L, j + 1) + i - pow(L, j), pow(L, j + 1));
+ unsigned n1 = uPow(L, j + 1) * (i / uPow(L, j + 1)) + (i + uPow(L, j)) % uPow(L, j + 1);
+ unsigned n2 = uPow(L, j + 1) * (i / uPow(L, j + 1)) + (uPow(L, j + 1) + i - uPow(L, j)) % uPow(L, j + 1);
unsigned e1 = j * Nv() + i;
unsigned e2 = j * Nv() + n2;
- HalfEdge he1(edges[e1], vertices[n1]);
- HalfEdge he2(edges[e2], vertices[n2]);
+ HalfEdge he1(e1, vertices[n1]);
+ HalfEdge he2(e2, vertices[n2]);
vertices[i].bonds.push_back(he1);
vertices[i].bonds.push_back(he2);
}
}
- }
- void dilute(Real p, Rng& r) {
- for (Edge& e : edges) {
- e.active = p > r.variate<Real, std::uniform_real_distribution>(0.0, 1.0);
+ for (const Vertex& v1 : vertices) {
+ for (const Vertex& v2 : vertices) {
+ unsigned dist = squaredDistance(v1.index, v2.index);
+ if (dist > 0) {
+ multiplicities[dist - 1]++;
+ }
+ }
}
}
- void markClusters() {
- clusters = {};
- for (Vertex& v : vertices) {
- v.cluster = 0;
- };
+ std::vector<Cluster> markClusters(const std::vector<short unsigned>& activeEdges) const {
+ std::vector<Cluster> clusters;
+ std::vector<unsigned> indicies(Nv());
unsigned nClusters = 0;
- for (Vertex& v : vertices) {
- if (v.cluster == 0) {
+ for (const Vertex& v : vertices) {
+ if (indicies[v.index] == 0) {
nClusters++;
clusters.push_back({});
- std::stack<std::reference_wrapper<Vertex>> inCluster;
+ std::stack<std::reference_wrapper<const Vertex>> inCluster;
inCluster.push(v);
while (!inCluster.empty()) {
- Vertex& vn = inCluster.top();
+ const Vertex& vn = inCluster.top();
inCluster.pop();
- if (vn.cluster == 0) {
- vn.cluster = nClusters;
+ if (indicies[vn.index] == 0) {
+ indicies[vn.index] = nClusters;
clusters[nClusters - 1].push_back(vn);
}
- for (HalfEdge& e : vn.bonds) {
- if (e.parent.active && e.neighbor.cluster == 0) {
+ for (const HalfEdge& e : vn.bonds) {
+ if (activeEdges[e.index] && indicies[e.neighbor.index] == 0) {
inCluster.push(e.neighbor);
}
}
}
}
}
+
+ return clusters;
}
- Matrix laplacian() const {
+ Matrix laplacian(const std::vector<short unsigned>& activeEdges) const {
Matrix L = Matrix::Zero(Nv(), Nv());
for (const Vertex& v : vertices) {
for (const HalfEdge& e : v.bonds) {
- if (e.parent.active) {
+ if (activeEdges[e.index]) {
L(v.index, v.index) += 1;
L(v.index, e.neighbor.index) -= 1;
}
@@ -175,27 +175,19 @@ int main(int argc, char* argv[]) {
}
}
- std::vector<unsigned> multiplicities(D * L * L / 4);
-
Graph G(D, L);
- for (const Vertex& v1 : G.vertices) {
- for (const Vertex& v2 : G.vertices) {
- unsigned dist = G.squaredDistance(v1.index, v2.index);
- if (dist > 0) {
- multiplicities[dist - 1]++;
- }
- }
- }
-
for (unsigned trials = 0; trials < n; trials++) {
- G.dilute(p, r);
- G.markClusters();
- Matrix laplacian = G.laplacian();
+ std::vector<short unsigned> activeEdges(G.Ne());
+ for (short unsigned& edge : activeEdges) {
+ edge = r.uniform(0.0, 1.0) < p;
+ }
+ std::vector<Cluster> clusters = G.markClusters(activeEdges);
+ Matrix laplacian = G.laplacian(activeEdges);
std::vector<Real> conductivities(L * L * D / 4);
- for (const Cluster& c : G.clusters) {
+ for (const Cluster& c : clusters) {
std::vector<unsigned> inds;
inds.reserve(c.size());
for (const Vertex& v : c){
@@ -221,8 +213,8 @@ int main(int argc, char* argv[]) {
}
for (unsigned i = 0; i < conductivities.size(); i++) {
- if (multiplicities[i] != 0) {
- std::cout << conductivities[i] / multiplicities[i] << " ";
+ if (G.multiplicities[i] != 0) {
+ std::cout << conductivities[i] / G.multiplicities[i] << " ";
} else {
std::cout << 0 << " ";
}