path: root/src/measurements.cpp

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

void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, std::string model_string) {
  std::string filename = model_string + 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, std::string model_string, unsigned Mx, unsigned My) {
  std::string filename = model_string + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
  std::ifstream file(filename);

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

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

  std::ofstream file_out(filename);

  file_out <<std::fixed<< N_old + 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<graph::coordinate>& pos) {
  for (std::list<graph::coordinate>::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) {
    for (std::list<graph::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;

      data[(unsigned)(Mx * (Δx / Lx)) + (Mx / 2) * (unsigned)(My * (Δy / Ly))]++;
    }
  }
}

void correlation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data, const std::list<graph::coordinate>& pos1, const std::list<graph::coordinate>& pos2) {
  for (std::list<graph::coordinate>::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) {
    for (std::list<graph::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(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(unsigned n, double a, unsigned Mx, unsigned My, double beta) :
  Mx(Mx), My(My), G(2 * n),
  sc(3 * n, 0),
  ss(3 * n, 0),
  sm(3 * n, 0),
  sa(3 * n, 0),
  sl(3 * n, 0),
  sb(3 * n, 0),
  sD(3 * n, 0),
  ccc((Mx / 2) * (My / 2), 0),
  css((Mx / 2) * (My / 2), 0),
  cmm((Mx / 2) * (My / 2), 0),
  caa((Mx / 2) * (My / 2), 0),
  cll((Mx / 2) * (My / 2), 0),
  cbb((Mx / 2) * (My / 2), 0),
  cDD((Mx / 2) * (My / 2), 0),
  csD((Mx / 2) * (My / 2), 0)
{
  N = 0;
  Nc = 0;
  Na = 0;
  Nb = 0;

  if (beta != 0.0) {
    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_";
  } else {
    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_INF_";
  }

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

ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) :
  Mx(Mx), My(My), 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),
  sb(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  sD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
  ccc((Mx / 2) * (My / 2), 0),
  css((Mx / 2) * (My / 2), 0),
  cmm((Mx / 2) * (My / 2), 0),
  caa((Mx / 2) * (My / 2), 0),
  cll((Mx / 2) * (My / 2), 0),
  cbb((Mx / 2) * (My / 2), 0),
  cDD((Mx / 2) * (My / 2), 0),
  csD((Mx / 2) * (My / 2), 0)
{
  N = 0;
  Nc = 0;
  Na = 0;
  Nb = 0;

  if (beta != 0.0) {
    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_";
  } else {
    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_";
  }

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

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, model_string);
  update_distribution_file("ss", ss, N, model_string);
  update_distribution_file("sm", sm, N, model_string);
  update_distribution_file("sa", sa, Na, model_string);
  update_distribution_file("sl", sl, N, model_string);
  update_distribution_file("sb", sb, Nb, model_string);
  update_distribution_file("sD", sD, N, model_string);

  update_field_file("ccc", ccc, Nc, model_string, Mx, My);
  update_field_file("css", css, N, model_string, Mx, My);
  update_field_file("cmm", cmm, N, model_string, Mx, My);
  update_field_file("caa", caa, Na, model_string, Mx, My);
  update_field_file("cll", cll, N, model_string, Mx, My);
  update_field_file("cbb", cbb, Nb, model_string, Mx, My);
  update_field_file("cDD", cDD, N, model_string, Mx, My);
  update_field_file("csD", csD, N, model_string, 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++;

    autocorrelation2(net.G.L.x, net.G.L.y, Mx, My, caa, avalanches.back());

    lv = c;
    avalanches.push_back({net.G.edges[i].r});
    last_avalanche = {i};
  } else {
    avalanches.back().push_back(net.G.edges[i].r);
    last_avalanche.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);

  ss[crack.size() - 1]++;
  std::list<graph::coordinate> sr;

  for (auto edge : crack) {
    sr.push_back(n.G.dual_edges[edge].r);
  }

  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, css, sr);

  unsigned crack_component = component[n.G.dual_edges[crack.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 && components[i].size() > 0) {
      sc[components[i].size() - 1]++;
      Nc++;

      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ccc, components[i]);
    }
  }

  // spanning cluster

  sm[components[crack_component].size() - 1]++;

  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cmm, components[crack_component]);

  /// damage at end
  std::list<graph::coordinate> Dr;

  for (unsigned i = 0; i < n.G.edges.size(); i++) {
    if (n.fuses[i]) { 
      Dr.push_back(n.G.edges[i].r);
    }
  }

  sD[Dr.size() - 1]++;

  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cDD, Dr);
  correlation2(n.G.L.x, n.G.L.y, Mx, My, csD, sr, Dr);

  // ********************** LAST AVALANCHE *********************

  // rewind the last avalanche
  sl[avalanches.back().size() - 1]++;

  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cll, avalanches.back());

  for (unsigned e : last_avalanche) {
    boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
  }

  num = boost::connected_components(G, &component[0]);

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

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

  for (unsigned i = 0; i < num; i++) {
    if (new_components[i].size() > 0) {
      sb[new_components[i].size() - 1]++;
      Nb++;

      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cbb, new_components[i]);
    }
  }

  N++;
}