modified algorthim to use Newman and Ziff's, started haskell implementation

2 files changed, 99 insertions, 43 deletions

diff --git a/percolation.cpp b/percolation.cpp
index f0b0610..3f82cc4 100644
--- a/percolation.cpp
+++ b/percolation.cpp
@@ -1,25 +1,67 @@
 
+#include <algorithm>
+#include <cinttypes>
 #include <fstream>
 #include <random>
-#include <cinttypes>
-#include <vector>
 #include <string>
-#include <algorithm>
+#include <vector>
 
 #include "randutils/randutils.hpp"
 
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/incremental_components.hpp>
-#include <boost/pending/disjoint_sets.hpp>
-#include <boost/foreach.hpp>
+class Tree {
+private:
+  std::vector<signed> p;
+  std::vector<unsigned> o;
+
+public:
+  unsigned largest;
+  std::vector<unsigned> c;
 
-typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS> Graph;
-typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
-typedef boost::graph_traits<Graph>::vertices_size_type VertexIndex;
-typedef VertexIndex* Rank;
-typedef Vertex* Parent;
+  Tree(unsigned n) : p(n, -1), o(n, 1), c(n, 0) {
+    largest = 1;
+    c[0] = n;
+  }
+
+  unsigned findroot(unsigned i) {
+    if (p[i] < 0)
+      return i;
+    else
+      return p[i] = this->findroot(p[i]);
+  }
+
+  void join(unsigned i, unsigned j) {
+    c[o[i] - 1]--;
+    c[o[j] - 1]--;
 
-void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data, std::string model_string) {
+    if (o[i] < o[j]) {
+      p[i] = j;
+      o[j] += o[i];
+      c[o[j] - 1]++;
+      if (o[j] > largest) {
+        largest = o[j];
+      }
+    } else {
+      p[j] = i;
+      o[i] += o[j];
+      c[o[i] - 1]++;
+      if (o[i] > largest) {
+        largest = o[i];
+      }
+    }
+  }
+
+  void add_bond(const std::array<unsigned, 2>& b) {
+    unsigned i = this->findroot(b[0]);
+    unsigned j = this->findroot(b[1]);
+
+    if (i != j) {
+      this->join(i, j);
+    }
+  }
+};
+
+void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data,
+                              std::string model_string) {
   std::string filename = model_string + id + ".dat";
   std::ifstream file(filename);
 
@@ -61,24 +103,29 @@ int main(int argc, char* argv[]) {
 
   while ((opt = getopt(argc, argv, "L:N:")) != -1) {
     switch (opt) {
-      case 'L':
-        L = (unsigned)atof(optarg);
-        break;
-      case 'N':
-        N = (unsigned)atof(optarg);
-        break;
-      default:
-        exit(1);
+    case 'L':
+      L = (unsigned)atof(optarg);
+      break;
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    default:
+      exit(1);
     }
   }
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
-  unsigned n_vertices = pow(L, 2);;
+  unsigned n_vertices = pow(L, 2);
   unsigned n_bonds = 2 * n_vertices;
 
   std::vector<std::array<unsigned, 2>> bonds(n_bonds);
+  std::vector<unsigned> bond_indices(n_bonds);
+
+  for (unsigned i = 0; i < n_bonds; i++) {
+    bond_indices[i] = i;
+  }
 
   for (unsigned v_i = 0; v_i < n_vertices; v_i++) {
     bonds[2 * v_i + 0] = {v_i, L * (v_i / L) + (v_i + 1) % L};
@@ -92,28 +139,14 @@ int main(int argc, char* argv[]) {
   }
 
   for (unsigned trial = 0; trial < N; trial++) {
-    Graph G(n_vertices);
-    std::vector<VertexIndex> rank(n_vertices);
-    std::vector<Vertex> parent(n_vertices);
-
-    boost::disjoint_sets<VertexIndex*, Vertex*> ds(&rank[0], &parent[0]);
-
-    initialize_incremental_components(G, ds);
-    incremental_components(G, ds);
-
-    std::shuffle(bonds.begin(), bonds.end(), rng);
+    Tree t(n_vertices);
+    std::shuffle(bond_indices.begin(), bond_indices.end(), rng);
 
     for (unsigned p = 0; p < n_bonds; p++) {
-      boost::add_edge(bonds[p][0], bonds[p][1], G);
-      ds.union_set(bonds[p][0], bonds[p][1]);
-
-      boost::component_index<VertexIndex> components(parent.begin(), parent.end());
-      BOOST_FOREACH(VertexIndex current_index, components) {
-        unsigned comp_size = 0;
-        BOOST_FOREACH(VertexIndex child_index, components[current_index]) {
-          comp_size++;
-        }
-        sn[p][comp_size - 1]++;
+      t.add_bond(bonds[bond_indices[p]]);
+
+      for (unsigned i = 0; i < t.largest; i++) {
+        sn[p][i] += t.c[i];
       }
     }
   }
@@ -123,4 +156,3 @@ int main(int argc, char* argv[]) {
   return 0;
 }
 
-
diff --git a/percolation.hs b/percolation.hs
new file mode 100644
index 0000000..1f66375
--- /dev/null
+++ b/percolation.hs
@@ -0,0 +1,24 @@
+{-# LANGUAGE BangPatterns #-}
+import Control.Monad.State (State, get, put, runState, state)
+import Data.Array (Array)
+
+data Tree = Root Int Int | Node Int Tree
+
+data Forest = Array Int Tree
+
+findRoot :: Tree -> (Int, [(Int, Tree)])
+findRoot (Root i n) = (i, [])
+findRoot (Node i t) = (r, [(i, Node i ] ++ [])
+  
+
+
+squareLatticeEdges :: Int -> Int -> [(Int, Int)]
+squareLatticeEdges d l =
+  [ (i, (l ^ j) * quot i (l ^ j) + mod (l ^ j + i + l ^ (j - 1)) (l ^ j))
+  | j <- [1 .. d]
+  , i <- [0 .. (l ^ d) - 1]
+  ]
+
+
+
+