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

void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, double Lx, double Ly, double beta) {
  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
  std::ifstream file(filename);

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

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

    file.close();
  }

  std::ofstream file_out(filename);

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

  file_out.close();
}

template <class T>
void update_field_file(std::string id, const std::vector<T>& data, unsigned N, double Lx, double Ly, double beta, unsigned Mx, unsigned My) {
  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
  std::ifstream file(filename);

  uint64_t N_old = 0;
  std::vector<std::vector<uint64_t>> data_old(data.size());
  for (unsigned j = 0; j < data.size(); j++) {
    data_old[j].resize(data[j].size(), 0);
  }

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

  std::ofstream file_out(filename);

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

  file_out.close();
}

fftw_complex* data_transform(unsigned Mx, unsigned My, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
  fftw_execute(forward_plan);

  fftw_complex* output = (fftw_complex*)malloc(Mx * (My / 2 + 1) * sizeof(fftw_complex));

  for (unsigned i = 0; i < Mx * (My / 2 + 1); i++) {
    output[i][0] = fftw_forward_out[i][0];
    output[i][1] = fftw_forward_out[i][1];
  }

  return output;
}

template <class T>
void correlation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& data, const fftw_complex* tx1, const fftw_complex* tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
  for (unsigned i = 0; i < Mx * (My / 2 + 1); i++) {
    fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1];
    fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0];
  }

  fftw_execute(reverse_plan);

  for (unsigned j = 0; j < data.size(); j++) {
    for (unsigned i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) {
      data[j][i] += (T)pow(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My), j + 1);
    }
  }
}

template <class T>
void autocorrelation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
  fftw_execute(forward_plan);

  for (unsigned i = 0; i < My * (Mx / 2 + 1); i++) {
    fftw_reverse_in[i][0] = pow(fftw_forward_out[i][0], 2) + pow(fftw_forward_out[i][1], 2);
    fftw_reverse_in[i][1] = 0.0;
  }

  fftw_execute(reverse_plan);

  for (unsigned j = 0; j < out_data.size(); j++) {
    for (unsigned i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) {
      out_data[j][i] += (T)pow(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My), j + 1);
    }
  }
}

unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, unsigned My) {
  return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
}

ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncum) :
  Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned)ceil(Lx * Ly / 2)),
  sc(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  ss(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  sm(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  sa(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  sl(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  sd(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  sD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  Sc(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  Ss(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  Sm(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  Sa(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  Sl(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  Sd(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  SD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  Ccc(Ncum),
  Css(Ncum),
  Cmm(Ncum),
  Caa(Ncum),
  Cll(Ncum),
  Cdd(Ncum),
  CDD(Ncum),
  CsD(Ncum),
  ccc(Ncum),
  css(Ncum),
  cmm(Ncum),
  caa(Ncum),
  cll(Ncum),
  cdd(Ncum),
  cDD(Ncum),
  csD(Ncum)
{
  N = 0;
  Nc = 0;
  Na = 0;
  for (unsigned i = 0; i < Ncum; i++) {
    Ccc[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    Css[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    Cmm[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    Caa[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    Cll[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    Cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    CDD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    CsD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    ccc[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    css[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    cmm[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    caa[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    cll[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    cDD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
    csD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
  }

  // FFTW setup for correlation functions
  fftw_set_timelimit(1);

  fftw_forward_in = (double *)fftw_malloc(Mx * My * sizeof(double));
  fftw_forward_out = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex));
  fftw_reverse_in = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex));
  fftw_reverse_out = (double *)fftw_malloc(Mx * My * sizeof(double));

  forward_plan = fftw_plan_dft_r2c_2d(My, Mx, fftw_forward_in, fftw_forward_out, 0);
  reverse_plan = fftw_plan_dft_c2r_2d(My, Mx, fftw_reverse_in, fftw_reverse_out, 0);
  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_current.dat";
  //stress_file.open(filename, std::ifstream::app);
}

ma::~ma() {
  // clean up FFTW objects
  fftw_free(fftw_forward_in);
  fftw_free(fftw_forward_out);
  fftw_free(fftw_reverse_in);
  fftw_free(fftw_reverse_out);
  fftw_destroy_plan(forward_plan);
  fftw_destroy_plan(reverse_plan);
  fftw_cleanup();

  update_distribution_file("sc", sc, Nc, Lx, Ly, beta);
  update_distribution_file("ss", ss, N, Lx, Ly, beta);
  update_distribution_file("sm", sm, N, Lx, Ly, beta);
  update_distribution_file("sa", sa, Na, Lx, Ly, beta);
  update_distribution_file("sl", sl, N, Lx, Ly, beta);
  update_distribution_file("sd", sd, N, Lx, Ly, beta);
  update_distribution_file("sD", sD, N, Lx, Ly, beta);

  update_distribution_file("Sc", Sc, Nc, Lx, Ly, beta);
  update_distribution_file("Ss", Ss, N, Lx, Ly, beta);
  update_distribution_file("Sm", Sm, N, Lx, Ly, beta);
  update_distribution_file("Sa", Sa, Na, Lx, Ly, beta);
  update_distribution_file("Sl", Sl, N, Lx, Ly, beta);
  update_distribution_file("Sd", Sd, N, Lx, Ly, beta);
  update_distribution_file("SD", SD, N, Lx, Ly, beta);

  update_field_file("Ccc", Ccc, Nc, Lx, Ly, beta, Mx, My);
  update_field_file("Css", Css, N, Lx, Ly, beta, Mx, My);
  update_field_file("Cmm", Cmm, N, Lx, Ly, beta, Mx, My);
  update_field_file("Caa", Caa, Na, Lx, Ly, beta, Mx, My);
  update_field_file("Cll", Cll, N, Lx, Ly, beta, Mx, My);
  update_field_file("Cdd", Cdd, N, Lx, Ly, beta, Mx, My);
  update_field_file("CDD", CDD, N, Lx, Ly, beta, Mx, My);
  update_field_file("CsD", CsD, N, Lx, Ly, beta, Mx, My);

  update_field_file("ccc", ccc, Nc, Lx, Ly, beta, Mx, My);
  update_field_file("css", css, N, Lx, Ly, beta, Mx, My);
  update_field_file("cmm", cmm, N, Lx, Ly, beta, Mx, My);
  update_field_file("caa", caa, Na, Lx, Ly, beta, Mx, My);
  update_field_file("cll", cll, N, Lx, Ly, beta, Mx, My);
  update_field_file("cdd", cdd, N, Lx, Ly, beta, Mx, My);
  update_field_file("cDD", cDD, N, Lx, Ly, beta, Mx, My);
  update_field_file("csD", csD, N, Lx, Ly, beta, Mx, My);

  //stress_file.close();
}

void ma::pre_fracture(const network&) {
  lv = std::numeric_limits<long double>::lowest();
  avalanches = {{}};
  boost::remove_edge_if(trivial, G);
}

void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
  long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i];
  if (c > lv && avalanches.back().size() > 0) {
    sa[avalanches.back().size() - 1]++;
    Na++;

    std::fill_n(fftw_forward_in, Mx * My, 0.0);

    unsigned avalanches_trunc = 0;

    for (auto e : avalanches.back()) {
      unsigned ind = edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My);
      if (fftw_forward_in[ind] == 0.0) {
        fftw_forward_in[ind] = 1.0;
        avalanches_trunc++;
      }
    }
    Sa[avalanches_trunc - 1]++;

    autocorrelation(Mx, My, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

    std::fill_n(fftw_forward_in, Mx * My, 0.0);

    for (auto e : avalanches.back()) {
      unsigned ind = edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My);
      fftw_forward_in[ind] += 1.0;
    }

    autocorrelation(Mx, My, caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

    lv = c;
    avalanches.push_back({i});
  } else {
    avalanches.back().push_back(i);
  }

  boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
}

void ma::post_fracture(network &n) {
  std::vector<unsigned> component(boost::num_vertices(G));
  unsigned num = boost::connected_components(G, &component[0]);

  std::list<unsigned> crack = find_minimal_crack(G, n);

  // crack surface correlations
  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  unsigned surface_size_trunc = 0;

  for (auto edge : crack) {
    unsigned ind1 = edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My);
    unsigned ind2 = edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My);
    if (fftw_forward_in[ind1] == 0.0) {
      fftw_forward_in[ind1] = 1.0;
      surface_size_trunc++;
    }
    if (fftw_forward_in[ind2] == 0.0) {
      fftw_forward_in[ind2] = 1.0;
      surface_size_trunc++;
    }
  }

  ss[crack.size()]++;

  fftw_complex *Tss = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);

  correlation(Mx, My, Css, Tss, Tss, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  for (auto edge : crack) {
    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] += 0.5;
    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] += 0.5;
  }

  Ss[surface_size_trunc]++;

  fftw_complex *tss = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);

  correlation(Mx, My, css, tss, tss, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];

  std::vector<std::list<unsigned>> components(num);

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

  // non-spanning clusters
  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  for (unsigned i = 0; i < num; i++) {
    if (i != crack_component && components[i].size() > 0) {
      unsigned cluster_size_trunc = 0;
      for (auto it = components[i].begin(); it != components[i].end(); it++) {
        unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My);
        if (fftw_forward_in[ind] == 0.0) {
          cluster_size_trunc++;
          fftw_forward_in[ind] = 1.0;
        }
      }

      Sc[cluster_size_trunc - 1]++;
      sc[components[i].size() - 1]++;

      autocorrelation(Mx, My, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
      Nc++;

      for (auto it = components[i].begin(); it != components[i].end(); it++) {
        fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 0.0;
      }

      for (auto it = components[i].begin(); it != components[i].end(); it++) {
        unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My);
        fftw_forward_in[ind] += 1.0;
      }

      autocorrelation(Mx, My, ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

      for (auto it = components[i].begin(); it != components[i].end(); it++) {
        fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 0.0;
      }
    }
  }

  // spanning cluster
  // std::fill_n(fftw_forward_in, Mx * My, 0.0); we already reset in the last loop

  sm[components[crack_component].size() - 1]++;
  unsigned spanning_cluster_trunc = 0;
  for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
    unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My);
    if (fftw_forward_in[ind] == 0.0) {
      fftw_forward_in[ind] = 1.0;
      spanning_cluster_trunc++;
    }
  }

  autocorrelation(Mx, My, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  Sm[spanning_cluster_trunc - 1]++;
  for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0;
  }

  autocorrelation(Mx, My, cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  /// damage at end
  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  unsigned final_broken = 0;
  unsigned final_trunc = 0;

  for (unsigned i = 0; i < n.G.edges.size(); i++) {
    if (n.fuses[i]) { 
      final_broken++;
      unsigned ind = edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My);
      if (fftw_forward_in[ind] == 0.0) {
        final_trunc++;
        fftw_forward_in[ind] = 1.0;
      }
    }
  }

  fftw_complex *TDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);

  SD[final_trunc]++;
  sD[final_broken]++;

  correlation(Mx, My, CDD, TDD, TDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
  correlation(Mx, My, CsD, Tss, TDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  free(Tss);
  free(TDD);

  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  for (unsigned i = 0; i < n.G.edges.size(); i++) {
    if (n.fuses[i]) { 
      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
    }
  }

  fftw_complex *tDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);

  correlation(Mx, My, cDD, tDD, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
  correlation(Mx, My, csD, tss, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  free(tss);
  free(tDD);

  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  // rewind the last avalanche
  sl[avalanches.back().size() - 1]++;
  unsigned last_trunc = 0;
  for (auto e : avalanches.back()) {
    boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
    n.break_edge(e, true);
    unsigned ind = edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My);
    if (fftw_forward_in[ind] == 0.0) {
      fftw_forward_in[ind] = 1.0;
      last_trunc++;
    }
  }
  Sl[last_trunc - 1]++;

  autocorrelation(Mx, My, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  // rewind the last avalanche
  for (auto e : avalanches.back()) {
    fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
  }

  autocorrelation(Mx, My, cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  // damage size distribution
  unsigned total_broken = 0;

  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  unsigned damage_trunc = 0;

  for (unsigned i = 0; i < n.G.edges.size(); i++) {
    if (n.fuses[i]) { 
      total_broken++;
      unsigned ind = edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My);
      if (fftw_forward_in[ind] == 0.0) {
        fftw_forward_in[ind] = 1.0;
        damage_trunc++;
      }
    }
  }
  Sd[damage_trunc]++;

  autocorrelation(Mx, My, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  std::fill_n(fftw_forward_in, Mx * My, 0.0);

  for (unsigned i = 0; i < n.G.edges.size(); i++) {
    if (n.fuses[i]) { 
      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
    }
  }

  autocorrelation(Mx, My, cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);

  sd[total_broken]++;

  N++;
}