Changed many variable names, and now save multiple moments to separate files.HEADmaster

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;