1 files changed, 126 insertions, 0 deletions
diff --git a/percolation.cpp b/percolation.cpp
new file mode 100644
index 0000000..f0b0610
--- /dev/null
+++ b/percolation.cpp
@@ -0,0 +1,126 @@
+
+#include <fstream>
+#include <random>
+#include <cinttypes>
+#include <vector>
+#include <string>
+#include <algorithm>
+
+#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>
+
+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;
+
+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);
+
+  std::vector<std::vector<uint64_t>> data_old(data.size());
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    data_old[i].resize(data[0].size());
+  }
+
+  if (file.is_open()) {
+    for (unsigned i = 0; i < data.size(); i++) {
+      for (unsigned j = 0; j < data[0].size(); j++) {
+        uint64_t num;
+        file >> num;
+        data_old[i][j] = num;
+      }
+    }
+
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    for (unsigned j = 0; j < data[0].size(); j++) {
+      file_out << std::fixed << data_old[i][j] + data[i][j] << " ";
+    }
+    file_out << "\n";
+  }
+
+  file_out.close();
+}
+
+int main(int argc, char* argv[]) {
+  unsigned L = 16;
+  unsigned N = 100;
+
+  int opt;
+
+  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);
+    }
+  }
+
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng{seeds};
+
+  unsigned n_vertices = pow(L, 2);;
+  unsigned n_bonds = 2 * n_vertices;
+
+  std::vector<std::array<unsigned, 2>> bonds(n_bonds);
+
+  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};
+    bonds[2 * v_i + 1] = {v_i, L * (((v_i / L) + 1) % L) + v_i % L};
+  }
+
+  std::vector<std::vector<uint64_t>> sn(n_bonds);
+
+  for (std::vector<uint64_t>& sni : sn) {
+    sni.resize(n_vertices);
+  }
+
+  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);
+
+    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]++;
+      }
+    }
+  }
+
+  update_distribution_file("sn", sn, "percolation_" + std::to_string(L) + "_");
+
+  return 0;
+}
+
+