1 files changed, 234 insertions, 0 deletions
diff --git a/percolation.cpp b/percolation.cpp
new file mode 100644
index 0000000..08c6970
--- /dev/null
+++ b/percolation.cpp
@@ -0,0 +1,234 @@
+#include <stack>
+#include <list>
+#include <cstdlib>
+#include <getopt.h>
+
+#include "eigen/Eigen/Dense"
+#include "randutils/randutils.hpp"
+#include "pcg-cpp/include/pcg_random.hpp"
+
+using Rng = randutils::random_generator<pcg32>;
+
+using Real = double;
+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;
+  Vertex& neighbor;
+  HalfEdge(Edge& e, Vertex& v) : parent(e), 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<Vertex>> {};
+class Clusters {
+public:
+  std::vector<Cluster> lists;
+  std::vector<unsigned> indices;
+};
+
+class Graph {
+public:
+  unsigned D;
+  unsigned L;
+  std::vector<Vertex> vertices;
+  std::vector<Edge> edges;
+  std::vector<Cluster> clusters;
+
+  unsigned Nv() const {
+    return vertices.size();
+  }
+
+  unsigned Ne() const {
+    return edges.size();
+  }
+
+  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));
+      unsigned Δx = abs(x1 - x2);
+      if (Δx > L / 2) {
+        Δx = L - Δx;
+      }
+      sd += pow(Δx, 2);
+    }
+    return sd;
+  }
+
+  /* Initialize a square lattice */
+  Graph(unsigned D, unsigned L) : D(D), L(L), vertices(pow(L, D)), edges(D * pow(L, D)) {
+    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 e1 = j * Nv() + i;
+        unsigned e2 = j * Nv() + n2;
+
+        HalfEdge he1(edges[e1], vertices[n1]);
+        HalfEdge he2(edges[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);
+    }
+  }
+
+  void markClusters() {
+    clusters = {};
+    for (Vertex& v : vertices) {
+      v.cluster = 0;
+    };
+    unsigned nClusters = 0;
+    for (Vertex& v : vertices) {
+      if (v.cluster == 0) {
+        nClusters++;
+        clusters.push_back({});
+        std::stack<std::reference_wrapper<Vertex>> inCluster;
+        inCluster.push(v);
+        while (!inCluster.empty()) {
+          Vertex& vn = inCluster.top();
+          inCluster.pop();
+          if (vn.cluster == 0) {
+            vn.cluster = nClusters;
+            clusters[nClusters - 1].push_back(vn);
+          }
+          for (HalfEdge& e : vn.bonds) {
+            if (e.parent.active && e.neighbor.cluster == 0) {
+              inCluster.push(e.neighbor);
+            }
+          }
+        }
+      }
+    }
+  }
+
+  Matrix laplacian() const {
+    Matrix L = Matrix::Zero(Nv(), Nv());
+
+    for (const Vertex& v : vertices) {
+      for (const HalfEdge& e : v.bonds) {
+        if (e.parent.active) {
+          L(v.index, v.index) += 1;
+          L(v.index, e.neighbor.index) -= 1;
+        }
+      }
+    }
+
+    return L;
+  }
+};
+
+int main(int argc, char* argv[]) {
+  Rng r;
+
+  unsigned D = 2;
+  unsigned L = 8;
+  double p = 0.5;
+  unsigned n = 10;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "D:L:p:n:")) != -1) {
+    switch (opt) {
+    case 'D':
+      D = atoi(optarg);
+      break;
+    case 'L':
+      L = atoi(optarg);
+      break;
+    case 'p':
+      p = atof(optarg);
+      break;
+    case 'n':
+      n = (unsigned)atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  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<Real> conductivities(L * L * D / 4);
+
+    for (const Cluster& c : G.clusters) {
+      std::vector<unsigned> inds;
+      inds.reserve(c.size());
+      for (const Vertex& v : c){
+        inds.push_back(v.index);
+      }
+
+      Matrix subLaplacian = laplacian(inds, inds);
+      subLaplacian(0,0)++; /* Set voltage of first node to zero */
+      Matrix subLaplacianInv = subLaplacian.inverse();
+
+      for (unsigned i = 0; i < c.size(); i++) {
+        for (unsigned j = 0; j < i; j++) {
+          Vector input = Vector::Zero(inds.size());
+          input(i) = 1;
+          input(j) = -1;
+
+          Vector output = subLaplacianInv * input;
+          double ΔV = std::abs(output(i) - output(j));
+
+          conductivities[G.squaredDistance(inds[i], inds[j]) - 1] += 1 / ΔV;
+        }
+      }
+    }
+
+    for (unsigned i = 0; i < conductivities.size(); i++) {
+      if (multiplicities[i] != 0) {
+        std::cout << conductivities[i] / multiplicities[i] << " ";
+      } else {
+        std::cout << 0 << " ";
+      }
+    }
+    std::cout << std::endl;
+  }
+
+  return 0;
+}