Removed percolation-specific data from graph class.

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 << " ";
       }