summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--percolation.cpp169
1 files changed, 87 insertions, 82 deletions
diff --git a/percolation.cpp b/percolation.cpp
index fac04ec..c1237e1 100644
--- a/percolation.cpp
+++ b/percolation.cpp
@@ -2,6 +2,7 @@
#include <list>
#include <cstdlib>
#include <getopt.h>
+#include <fstream>
#include "eigen/Eigen/Dense"
#include "randutils/randutils.hpp"
@@ -9,6 +10,7 @@
using Rng = randutils::random_generator<pcg32>;
+using Bool = short unsigned;
using Real = double;
using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
@@ -24,8 +26,8 @@ public:
class Vertex {
public:
- std::vector<HalfEdge> bonds;
unsigned index;
+ std::vector<HalfEdge> outgoingEdges;
std::vector<unsigned> coordinates;
};
@@ -58,73 +60,72 @@ public:
return D * Nv();
}
- unsigned squaredDistance(unsigned vi, unsigned vj) const {
- unsigned sd = 0;
- for (unsigned i = 0; i < D; i++) {
- int x1 = vertices[vi].coordinates[i];
- int x2 = vertices[vj].coordinates[i];
- unsigned Δx = abs(x1 - x2);
+ unsigned squaredDistance(unsigned i₁, unsigned i₂) const {
+ unsigned Δx² = 0;
+ for (unsigned d = 0; d < D; d++) {
+ int x₁ = vertices[i₁].coordinates[d];
+ int x₂ = vertices[i₂].coordinates[d];
+ unsigned Δx = abs(x₁ - x₂);
if (Δx > L / 2) {
Δx = L - Δx;
}
- sd += Δx * Δx;
+ Δx² += Δx * Δx;
}
- return sd;
+ return Δx²;
}
- /* Initialize a square lattice */
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].coordinates.resize(D);
vertices[i].index = i;
- vertices[i].bonds.reserve(2 * D);
- for (unsigned j = 0; j < D; j++) {
- vertices[i].coordinates[j] = (i / uPow(L, j)) % L;
+ vertices[i].coordinates.reserve(D);
+ vertices[i].outgoingEdges.reserve(2 * D);
+ for (unsigned d = 0; d < D; d++) {
+ vertices[i].coordinates.push_back((i / uPow(L, d)) % L);
- 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 j₁ = uPow(L, d + 1) * (i / uPow(L, d + 1)) + (i + uPow(L, d)) % uPow(L, d + 1);
+ unsigned j₂ = uPow(L, d + 1) * (i / uPow(L, d + 1)) + (uPow(L, d + 1) + i - uPow(L, d)) % uPow(L, d + 1);
- unsigned e1 = j * Nv() + i;
- unsigned e2 = j * Nv() + n2;
+ unsigned e₁ = d * Nv() + i;
+ unsigned e₂ = d * Nv() + j₂;
- HalfEdge he1(e1, vertices[n1]);
- HalfEdge he2(e2, vertices[n2]);
+ HalfEdge he₁(e₁, vertices[j₁]);
+ HalfEdge he₂(e₂, vertices[j₂]);
- vertices[i].bonds.push_back(he1);
- vertices[i].bonds.push_back(he2);
+ vertices[i].outgoingEdges.push_back(he₁);
+ vertices[i].outgoingEdges.push_back(he₂);
}
}
- for (const Vertex& v1 : vertices) {
- for (const Vertex& v2 : vertices) {
- unsigned dist = squaredDistance(v1.index, v2.index);
- if (dist > 0) {
- multiplicities[dist - 1]++;
+ for (const Vertex& v₁ : vertices) {
+ for (const Vertex& v₂ : vertices) {
+ unsigned Δx² = squaredDistance(v₁.index, v₂.index);
+ if (Δx² > 0) {
+ multiplicities[Δx² - 1]++;
}
}
}
}
- std::vector<Cluster> markClusters(const std::vector<short unsigned>& activeEdges) const {
- std::vector<Cluster> clusters;
- std::vector<unsigned> indicies(Nv());
+ std::list<Cluster> findClusters(const std::vector<Bool>& activeEdges) const {
+ std::list<Cluster> clusters;
+ std::vector<unsigned> clusterIndicies(Nv());
unsigned nClusters = 0;
- for (const Vertex& v : vertices) {
- if (indicies[v.index] == 0) {
+ for (const Vertex& v₀ : vertices) {
+ if (clusterIndicies[v₀.index] == 0) {
nClusters++;
clusters.push_back({});
- std::stack<std::reference_wrapper<const Vertex>> inCluster;
- inCluster.push(v);
- while (!inCluster.empty()) {
- const Vertex& vn = inCluster.top();
- inCluster.pop();
- if (indicies[vn.index] == 0) {
- indicies[vn.index] = nClusters;
- clusters[nClusters - 1].push_back(vn);
- }
- for (const HalfEdge& e : vn.bonds) {
- if (activeEdges[e.index] && indicies[e.neighbor.index] == 0) {
- inCluster.push(e.neighbor);
+ std::stack<std::reference_wrapper<const Vertex>> verticiesToProcess;
+ verticiesToProcess.push(v₀);
+ while (!verticiesToProcess.empty()) {
+ const Vertex& v₁ = verticiesToProcess.top();
+ verticiesToProcess.pop();
+ if (clusterIndicies[v₁.index] == 0) {
+ clusterIndicies[v₁.index] = nClusters;
+ clusters.back().push_back(v₁);
+ for (const HalfEdge& e : v₁.outgoingEdges) {
+ if (activeEdges[e.index] && clusterIndicies[e.neighbor.index] == 0) {
+ verticiesToProcess.push(e.neighbor);
+ }
}
}
}
@@ -134,11 +135,11 @@ public:
return clusters;
}
- Matrix laplacian(const std::vector<short unsigned>& activeEdges) const {
+ Matrix laplacian(const std::vector<Bool>& activeEdges) const {
Matrix L = Matrix::Zero(Nv(), Nv());
for (const Vertex& v : vertices) {
- for (const HalfEdge& e : v.bonds) {
+ for (const HalfEdge& e : v.outgoingEdges) {
if (activeEdges[e.index]) {
L(v.index, v.index) += 1;
L(v.index, e.neighbor.index) -= 1;
@@ -153,14 +154,15 @@ public:
int main(int argc, char* argv[]) {
Rng r;
- unsigned D = 2;
- unsigned L = 8;
- double p = 0.5;
- unsigned n = 10;
+ unsigned D = 2; // Dimension of space
+ unsigned L = 8; // Size of lattice
+ Real p = 0.5; // Probability of bond
+ unsigned M = 1; // Moments of conductivity
+ unsigned n = 1; // Number of experiments
int opt;
- while ((opt = getopt(argc, argv, "D:L:p:n:")) != -1) {
+ while ((opt = getopt(argc, argv, "D:L:p:n:M:")) != -1) {
switch (opt) {
case 'D':
D = atoi(optarg);
@@ -171,6 +173,9 @@ int main(int argc, char* argv[]) {
case 'p':
p = atof(optarg);
break;
+ case 'M':
+ M = atoi(optarg);
+ break;
case 'n':
n = (unsigned)atof(optarg);
break;
@@ -181,59 +186,59 @@ int main(int argc, char* argv[]) {
Graph G(D, L);
- for (unsigned trials = 0; trials < n; trials++) {
- std::vector<short unsigned> activeEdges(G.Ne());
- for (short unsigned& edge : activeEdges) {
+ for (unsigned trial = 0; trial < n; trial++) {
+ std::vector<Bool> activeEdges(G.Ne());
+ for (Bool& edge : activeEdges) {
edge = r.uniform(0.0, 1.0) < p;
}
- std::vector<Cluster> clusters = G.markClusters(activeEdges);
+
+ std::list<Cluster> clusters = G.findClusters(activeEdges);
Matrix laplacian = G.laplacian(activeEdges);
- std::vector<Real> conductivities((L * L * D) / 4);
- std::vector<unsigned> probabilities((L * L * D) / 4);
+ std::vector<std::vector<Real>> conductivities(M + 1);
+
+ for (std::vector<Real>& dat : conductivities) {
+ dat.resize((L * L * D) / 4);
+ }
- for (const Cluster& c : clusters) {
- std::vector<unsigned> inds;
- inds.reserve(c.size());
- for (const Vertex& v : c){
- inds.push_back(v.index);
+ for (const Cluster& cluster : clusters) {
+ std::vector<unsigned> clusterIndicies;
+ clusterIndicies.reserve(cluster.size());
+ for (const Vertex& v : cluster){
+ clusterIndicies.push_back(v.index);
}
- Matrix subLaplacian = laplacian(inds, inds);
+ Matrix subLaplacian = laplacian(clusterIndicies, clusterIndicies);
subLaplacian(0,0)++; /* Set voltage of first node to zero */
Matrix subLaplacianInv = subLaplacian.inverse();
- for (unsigned i = 0; i < c.size(); i++) {
+ for (unsigned i = 0; i < cluster.size(); i++) {
for (unsigned j = 0; j < i; j++) {
- Vector input = Vector::Zero(inds.size());
+ Vector input = Vector::Zero(clusterIndicies.size());
input(i) = 1;
input(j) = -1;
Vector output = subLaplacianInv * input;
- double ΔV = std::abs(output(i) - output(j));
+ Real ΔV = std::abs(output(i) - output(j));
- conductivities[G.squaredDistance(inds[i], inds[j]) - 1] += 1 / ΔV;
- probabilities[G.squaredDistance(inds[i], inds[j]) - 1]++;
+ for (unsigned m = 0; m <= M; m++) {
+ conductivities[m][G.squaredDistance(clusterIndicies[i], clusterIndicies[j]) - 1] += pow(1 / ΔV, m);
+ }
}
}
}
- for (unsigned i = 0; i < conductivities.size(); i++) {
- if (G.multiplicities[i] != 0) {
- std::cout << conductivities[i] / G.multiplicities[i] << " ";
- } else {
- std::cout << 0 << " ";
- }
- }
- std::cout << std::endl;
- for (unsigned i = 0; i < probabilities.size(); i++) {
- if (G.multiplicities[i] != 0) {
- std::cout << ((Real)probabilities[i]) / G.multiplicities[i] << " ";
- } else {
- std::cout << 0 << " ";
+ for (unsigned m = 0; m <= M; m++) {
+ std::ofstream file("perc_" + std::to_string(D) + "_" + std::to_string(L) + "_" + std::to_string(p) + "_" + std::to_string(m) + ".dat", std::ios_base::app);
+ for (unsigned i = 0; i < conductivities[m].size(); i++) {
+ if (G.multiplicities[i] != 0) {
+ file << conductivities[m][i] / G.multiplicities[i] << " ";
+ } else {
+ file << 0 << " ";
+ }
}
+ file << std::endl;
}
- std::cout << std::endl;
}
return 0;