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