#include "perc_meas.hpp"
#include <iostream>
#include <cstdio>

class badcycleException: public std::exception
{
  virtual const char* what() const throw()
  {
    return "Could not find a valid cycle on the broken system.";
  }
} badcycleex;

void update_distribution_file(std::string id, const std::vector<uint64_t>& data, std::string model_string) {
  std::string filename = model_string + id + ".dat";
  std::ifstream file(filename);

  std::vector<uint64_t> data_old(data.size(), 0);

  if (file.is_open()) {
    for (unsigned i = 0; i < data.size(); i++) {
      uint64_t num;
      file >> num;
      data_old[i] = num;
    }

    file.close();
  }

  std::ofstream file_out(filename);

  for (unsigned i = 0; i < data.size(); i++) {
    file_out <<std::fixed<< data_old[i] + data[i] << " ";
  }

  file_out.close();
}

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();
}

pm::pm(unsigned n, double a) :
  G(2 * n),
  rank(2 * n),
  parent(2 * n),
  ds(&rank[0], &parent[0]),
  sn(3 * n),
  sN(3 * n),
  ss(2 * n),
  sm(2 * n),
  sl(2 * n),
  sb(3 * n),
  sd(3 * n)
{
  model_string = "percolation_" + std::to_string(n) + "_" + std::to_string(a) + "_";
  for (std::vector<uint64_t> &x : sn) {
    x.resize(2 * n);
  }
  for (std::vector<uint64_t> &x : sN) {
    x.resize(2 * n);
  }
  for (std::vector<uint64_t> &x : sb) {
    x.resize(n);
  }
}

pm::pm(unsigned Lx, unsigned Ly) :
  G(Lx * Ly / 2),
  rank(Lx * Ly / 2),
  parent(Lx * Ly / 2),
  ds(&rank[0], &parent[0]),
  sn(Lx * Ly),
  sN(Lx * Ly),
  ss(Lx * Ly / 2),
  sm(Lx * Ly / 2),
  sl(Lx * Ly / 2),
  sb(Lx * Ly),
  sd(Lx * Ly)
{
  model_string = "percolation_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_";
  for (std::vector<uint64_t> &x : sn) {
    x.resize(Lx * Ly / 2);
  }
  for (std::vector<uint64_t> &x : sN) {
    x.resize(Lx * Ly / 2);
  }
  for (std::vector<uint64_t> &x : sb) {
    x.resize(Lx * Ly / 2);
  }
}

pm::~pm() {
  update_distribution_file("sn", sn, model_string);
  update_distribution_file("sN", sN, model_string);
  update_distribution_file("ss", ss, model_string);
  update_distribution_file("sm", sm, model_string);
  update_distribution_file("sl", sl, model_string);
  update_distribution_file("sb", sb, model_string);
  update_distribution_file("sd", sd, model_string);
}

void pm::pre_fracture(const network&) {
  boost::remove_edge_if(trivial, G);
  initialize_incremental_components(G, ds);
  incremental_components(G, ds);
  r = 0;
  sN_tmp.clear();
}

void pm::bond_broken(const network& net, const current_info& cur, unsigned i) {
  boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
  ds.union_set(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1]);

  boost::component_index<VertexIndex> components(parent.begin(), parent.end());
  std::vector<unsigned> counts(components.size());

  sN_tmp.push_front({});

  BOOST_FOREACH(VertexIndex current_index, components) {
    unsigned comp_size = 0;
    BOOST_FOREACH(VertexIndex child_index, components[current_index]) {
      comp_size++;
    }

    sn[r][comp_size - 1]++;
    sN_tmp.front().push_back(comp_size - 1);
  }

  std::vector<bool> vertex_in(net.G.vertices.size());

  for (unsigned i = 0; i < net.G.edges.size(); i++) {
    if (cur.currents[i] > CURRENT_CUTOFF) {
      vertex_in[net.G.edges[i].v[0]] = true;
      vertex_in[net.G.edges[i].v[1]] = true;
    }
  }

  unsigned bb_size = 0;

  for (unsigned i = 0; i < net.G.vertices.size(); i++) {
    if (vertex_in[i]) bb_size++;
  }

  sb[r][bb_size - 1]++;

  r++;
  last_cur = cur;
}

void pm::post_fracture(network &n) {
  auto post_cracks = find_minimal_crack(G, n);
  std::vector<unsigned> component(boost::num_vertices(G));
  unsigned num = boost::connected_components(G, &component[0]);
  if (post_cracks.size() > 2 || post_cracks.size() == 0) {
    throw badcycleex;
  }
  for (auto c : post_cracks) {
    sl[c.second.size() - 1]++;
  }
  unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]];

  std::vector<std::list<graph::coordinate>> components(num);

  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
    components[component[i]].push_back(n.G.dual_vertices[i].r);
  }

  for (unsigned i = 0; i < num; i++) {
    if (i != crack_component) {
      sm[components[i].size() - 1]++;
    } else {
      ss[components[i].size() - 1]++;
    }
  }

  unsigned dr = 0;
  for (std::list<unsigned> l : sN_tmp) {
    for (unsigned size : l) {
      sN[dr][size]++;
    }
    dr++;
  }

  sd[r - 1]++;
}