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

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

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

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

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

  std::ofstream file_out(filename);

  file_out << old_N + N << "\n";
  for (unsigned j = 0; j < data.size(); j++) {
    file_out << data_old[j] + data[j] << " ";
  }

  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);
    }
  }
}
*/

void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data,
                      const std::list<coordinate>& pos, std::array<unsigned, 2> count) {
  for (std::list<coordinate>::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) {
    for (std::list<coordinate>::const_iterator it2 = it1; it2 != pos.end(); it2++) {
      double Δx_tmp = fabs(it1->x - it2->x);
      double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp;

      double Δy_tmp = fabs(it1->y - it2->y);
      double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp;

      if (count[0] % 2 == 0) {
        data[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δx / Lx)) +
             (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δy / Ly))]++;
      } else {
        data[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δy / Ly)) +
             (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δx / Lx))]++;
      }
    }
  }
}

void correlation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data,
                  const std::list<coordinate>& pos1,
                  const std::list<coordinate>& pos2) {
  for (std::list<coordinate>::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) {
    for (std::list<coordinate>::const_iterator it2 = pos2.begin(); it2 != pos2.end();
         it2++) {
      double Δx_tmp = fabs(it1->x - it2->x);
      double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp;

      double Δy_tmp = fabs(it1->y - it2->y);
      double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp;

      data[floor((Mx * Δx) / Lx) + (Mx / 2) * floor((My * Δy) / Ly)]++;
    }
  }
}
/*
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(coordinate r, double Lx, double Ly, unsigned Mx, unsigned My) {
  return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
}

ma::ma(unsigned n, double a, double beta, double weight)
    : sc(2 * n), sn(2 * n), ss(2 * n), sm(2 * n), sl(2 * n), sb(n + 1), sd(3 * n), sa(3 * n),
      sA(3 * n), sp(3 * n), si(10000), sI(10000), cc(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
      cn(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), cs(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
      cm(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), cl(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
      cb(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), ca(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
      cA(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), cq(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
      cS(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), NS(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), last_clusters(2 * n) {
  if (beta != 0.0) {
    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" +
                   std::to_string(beta) + "_" + std::to_string(weight) + "_";
  } else {
    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_INF_" +
                   std::to_string(weight) + "_";
  }
  Mx = 4 * ceil(sqrt(n)) * a;
  My = 4 * ceil(sqrt(n)) / a;
  Nc = 0;
  Nl = 0;
  Nm = 0;
  Ns = 0;
  Nn = 0;
  Nb = 0;
  Na = 0;
  NA = 0;
  Nq = 0;
}

ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
    : sc(Lx * Ly / 2), sn(Lx * Ly / 2), ss(Lx * Ly / 2), sm(Lx * Ly / 2), sl(Lx * Ly / 2),
      sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), sp(Lx * Ly), si(10000), sI(10000),
      cc((Lx / 2 + 1) * (Ly / 2 + 1)), cn((Lx / 2 + 1) * (Ly / 2 + 1)),
      cs((Lx / 2 + 1) * (Ly / 2 + 1)), cm((Lx / 2 + 1) * (Ly / 2 + 1)),
      cl((Lx / 2 + 1) * (Ly / 2 + 1)), cb((Lx / 2 + 1) * (Ly / 2 + 1)),
      ca((Lx / 2 + 1) * (Ly / 2 + 1)), cA((Lx / 2 + 1) * (Ly / 2 + 1)),
      cq((Lx / 2 + 1) * (Ly / 2 + 1)), cS(Lx * Ly / 2), NS(Lx * Ly / 2), last_clusters(Lx * Ly / 2) {
  if (beta != 0.0) {
    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" +
                   std::to_string(beta) + "_" + std::to_string(weight) + "_";
  } else {
    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_" +
                   std::to_string(weight) + "_";
  }
  Mx = Lx;
  My = Ly;
  Nc = 0;
  Nl = 0;
  Nm = 0;
  Ns = 0;
  Nn = 0;
  Nb = 0;
  Na = 0;
  NA = 0;
  Nq = 0;
}

ma::~ma() {
  update_distribution_file("sc", sc, 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);
  update_distribution_file("sa", sa, model_string);
  update_distribution_file("sA", sA, model_string);
  update_distribution_file("sp", sp, model_string);
  update_distribution_file("si", si, model_string);
  update_distribution_file("sI", sI, model_string);
  update_field_file("cc", Nc, cc, model_string);
  update_field_file("cl", Nl, cl, model_string);
  update_field_file("cm", Nm, cm, model_string);
  update_field_file("cs", Ns, cs, model_string);
  update_field_file("cn", Nn, cn, model_string);
  update_field_file("cb", Nb, cb, model_string);
  update_field_file("ca", Na, ca, model_string);
  update_field_file("cA", NA, cA, model_string);
  update_field_file("cq", Nq, cq, model_string);
  update_field_file("cS", 0, cS, model_string);
  update_field_file("NS", 0, NS, model_string);
}

void ma::pre_fracture(const network&) {
  lv = std::numeric_limits<long double>::lowest();
  avalanches = {};
  num = 0;
}

void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
  long double c = net.thresholds[i] - logl(cur.currents[i]);
  if (c > lv) {
    last_clusters = net.C;
    last_currents = cur;
    last_backbone = net.backbone;
    lv = c;
    avalanches.push_back({i});
  } else {
    avalanches.back().push_back(i);
  }

  for (unsigned j = 0; j < cur.currents.size(); j++) {
    if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 &&
        (!net.backbone[j] || j == i)) {
      si[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++;
    } else if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && net.backbone[j] &&
               j != i) {
      sI[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++;
    }
  }

  num++;
}

void ma::post_fracture(network& n) {
  auto crack = find_minimal_crack(n, avalanches.back().back());

  std::list<coordinate> cl_cs;
  sl[crack.second.size() - 1]++;
  for (unsigned e : crack.second) {
    cl_cs.push_back(n.G.dual_edges[e].r);
  }
  Nl += crack.second.size();
  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cl, cl_cs, crack.first);

  std::vector<std::list<coordinate>> crit_components(n.G.dual_vertices.size());

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

  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
    if (crit_components[i].size() > 0) {
      sc[crit_components[i].size() - 1]++;
    }
    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cc, crit_components[i], crack.first);
  }
  Nc += n.G.dual_vertices.size();

  std::vector<std::list<coordinate>> components(n.G.dual_vertices.size());

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

  unsigned crack_component = n.C.findroot(n.G.dual_edges[avalanches.back().back()].v[0]);

  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
    if (components[i].size() > 0) {
      if (i != crack_component) {
        sm[components[i].size() - 1]++;
        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cm, components[i], crack.first);
      } else {
        ss[components[i].size() - 1]++;
        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cs, components[i], crack.first);
        Ns += components[i].size();
      }
      sn[components[i].size() - 1]++;
      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cn, components[i], crack.first);
    }
  }
  Nn += n.G.dual_vertices.size();

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

  for (unsigned i = 0; i < n.G.edges.size(); i++) {
    if (!last_backbone[i]) {
      vertex_in[n.G.edges[i].v[0]] = true;
      vertex_in[n.G.edges[i].v[1]] = true;
    }
  }

  unsigned bb_size = 0;

  std::list<coordinate> cb_co;

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

  sb[bb_size]++;
  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cb, cb_co, crack.first);
  Nb += bb_size;

  auto av_it = avalanches.rbegin();
  av_it++;

  while (av_it != avalanches.rend()) {
    sa[(*av_it).size() - 1]++;
    Na += (*av_it).size();
    Nq += (*av_it).size();
    std::list<coordinate> ca_co;
    for (unsigned e : (*av_it)) {
      ca_co.push_back(n.G.edges[e].r);
    }
    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ca, ca_co, crack.first);
    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cq, ca_co, {0, 1});
    av_it++;
  }

  sA[avalanches.back().size() - 1]++;
  std::list<coordinate> cA_co;
  for (unsigned e : avalanches.back()) {
    cA_co.push_back(n.G.edges[e].r);
  }
  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cA, cA_co, crack.first);
  NA += avalanches.back().size();

  sd[num - 1]++;

  // Currents relative to tip of largest crack fragment
  std::map<unsigned, std::set<unsigned>> fragments;

  for (unsigned e : crack.second) {
    if (last_backbone[e]) {
      unsigned root = last_clusters.findroot(n.G.dual_edges[e].v[0]);
      if (fragments.contains(root)) {
        fragments[root].insert(e);
      } else {
        fragments[root] = {e};
      }
    }
  }

  if (fragments.size() > 0) {
    std::set<unsigned> biggest_fragment =
        std::max_element(fragments.begin(), fragments.end(), [](const auto& l1, const auto& l2) {
          return l1.second.size() < l2.second.size();
        })->second;

    sp[biggest_fragment.size()]++;

    // We have all the edges in the biggest fragment, in no meaningful order.
    // To find its ends, we will start at a random place and traverse the
    // fragment until we cannot continue. Then, we go back from where we
    // started and do the same in the other direction. Along the way we will
    // compute the (unwrapped) displacement between the two ends so that we
    // will know which is to the left and which is to the right.

    unsigned ce = *biggest_fragment.begin();
    unsigned cv = n.G.dual_edges[ce].v[0];
    coordinate Δr = n.G.Δr(ce); // direction is from v0 to v1

    while (true) {
      signed ne = -1;
      unsigned nv;

      for (unsigned e : n.G.dual_vertices[cv].ne) {
        if (e != ce && biggest_fragment.contains(e)) {
          ne = e;
          unsigned v0 = n.G.dual_edges[ne].v[0];
          unsigned v1 = n.G.dual_edges[ne].v[1];
          Δr += v0 == cv ? n.G.Δr(ne) : -n.G.Δr(ne);
          nv = v0 == cv ? v1 : v0;

          break;
        }
      }

      if (ne < 0) {
        break;
      } else {
        ce = ne;
        cv = nv;
      }
    }

    unsigned end_1e = ce;
    unsigned end_1v = cv;

    ce = *biggest_fragment.begin();
    cv = n.G.dual_edges[ce].v[1];

    while (true) {
      signed ne = -1;
      unsigned nv;

      for (unsigned e : n.G.dual_vertices[cv].ne) {
        if (e != ce && biggest_fragment.contains(e)) {
          ne = e;
          unsigned v0 = n.G.dual_edges[ne].v[0];
          unsigned v1 = n.G.dual_edges[ne].v[1];
          Δr += v0 == cv ? -n.G.Δr(ne) : n.G.Δr(ne);
          nv = v0 == cv ? v1 : v0;

          break;
        }
      }

      if (ne < 0) {
        break;
      } else {
        ce = ne;
        cv = nv;
      }
    }

    unsigned end_0e = ce;
    unsigned end_0v = cv;

    bool comp;

    if (crack.first[0] % 2 == 1) {
      comp = Δr.x > 0;
    } else {
      comp = Δr.y > 0;
    }

    unsigned left_end, right_end;

    if (comp) {
      left_end = end_0e;
      right_end = end_1e;
    } else {
      left_end = end_1e;
      right_end = end_0e;
    }

    for (unsigned i = 0; i < n.G.dual_edges.size(); i++) {
      if (!last_backbone[i]) {
        double Δx_tmpl, Δy_tmpl, Δx_tmpr, Δy_tmpr, Δxr, Δxl, Δyr, Δyl;

        if (crack.first[0] % 2 == 1) {
          Δx_tmpl = n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x;
          Δxl = Δx_tmpl < 0 ? Δx_tmpl + n.G.L.x : Δx_tmpl;

          Δy_tmpl = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
          Δyl = Δy_tmpl > n.G.L.y / 2 ? n.G.L.y - Δy_tmpl : Δy_tmpl;

          Δx_tmpr = n.G.dual_edges[i].r.x - n.G.dual_edges[right_end].r.x;
          Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.x : Δx_tmpr;

          Δy_tmpr = fabs(n.G.dual_edges[right_end].r.y - n.G.dual_edges[i].r.y);
          Δyr = Δy_tmpr > n.G.L.y / 2 ? n.G.L.y - Δy_tmpr : Δy_tmpr;
        } else {
          Δx_tmpl = n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y;
          Δxl = Δx_tmpl < 0 ? Δx_tmpl + n.G.L.y : Δx_tmpl;

          Δy_tmpl = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
          Δyl = Δy_tmpl > n.G.L.x / 2 ? n.G.L.x - Δy_tmpl : Δy_tmpl;

          Δx_tmpr = n.G.dual_edges[i].r.y - n.G.dual_edges[right_end].r.y;
          Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.y : Δx_tmpr;

          Δy_tmpr = fabs(n.G.dual_edges[right_end].r.x - n.G.dual_edges[i].r.x);
          Δyr = Δy_tmpr > n.G.L.x / 2 ? n.G.L.x - Δy_tmpr : Δy_tmpr;
        }

        cS[(unsigned)floor(Mx * (Δyl / n.G.L.y)) +
           (Mx / 2) * (unsigned)floor(My * (Δxl / n.G.L.x))] +=
            last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
        NS[(unsigned)floor(Mx * (Δyl / n.G.L.y)) +
           (Mx / 2) * (unsigned)floor(My * (Δxl / n.G.L.x))]++;
        NS[(unsigned)floor(Mx * (Δyr / n.G.L.y)) +
           (Mx / 2) * (unsigned)floor(My * (Δxr / n.G.L.x))]++;
        cS[(unsigned)floor(Mx * (Δyr / n.G.L.y)) +
           (Mx / 2) * (unsigned)floor(My * (Δxr / n.G.L.x))] +=
            last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
      }
    }

  } else {
    sp[0]++;
  }

  aG.push_back(last_currents.conductivity);
}