#include <fftw3.h>

#include "fracture.h"

int main(int argc, char *argv[]) {
  int opt;

  fftw_set_timelimit(1);

  // defining variables to be (potentially) set by command line flags
  uint8_t filename_len;
  uint32_t N;
  uint_t L;
  double beta, inf, cutoff, crack_len;
  bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual,
      save_network, save_crit_stress, save_energy, save_conductivity,
      save_damage, save_threshold, save_current_load, save_correlations;
  bound_t boundary;
  lattice_t lattice;

  // assume filenames will be less than 100 characters

  filename_len = 100;

  // set default values

  N = 100;
  L = 16;
  crack_len = 0.;
  beta = .3;
  inf = 1e10;
  cutoff = 1e-9;
  boundary = FREE_BOUND;
  lattice = VORONOI_LATTICE;
  save_data = false;
  save_cluster_dist = false;
  use_voltage_boundaries = false;
  use_dual = false;
  save_network = false;
  save_crit_stress = false;
  save_damage = false;
  save_conductivity = false;
  save_energy = false;
  save_threshold = false;
  save_current_load = false;
  save_correlations = false;

  uint8_t bound_i;
  char boundc2 = 'f';
  uint8_t lattice_i;
  char lattice_c = 'v';
  char dual_c = 'o';

  // get commandline options

  while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TEz")) != -1) {
    switch (opt) {
    case 'n':
      N = atoi(optarg);
      break;
    case 'L':
      L = atoi(optarg);
      break;
    case 'b':
      beta = atof(optarg);
      break;
    case 'l':
      crack_len = atof(optarg);
      break;
    case 'B':
      bound_i = atoi(optarg);
      switch (bound_i) {
      case 0:
        boundary = FREE_BOUND;
        boundc2 = 'f';
        break;
      case 1:
        boundary = CYLINDER_BOUND;
        boundc2 = 'c';
        break;
      case 2:
        boundary = TORUS_BOUND;
        use_voltage_boundaries = true;
        boundc2 = 't';
        break;
      case 3:
        boundary = EMBEDDED_BOUND;
        boundc2 = 'e';
        use_dual = true;
        use_voltage_boundaries = true;
        break;
      default:
        printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), "
               "or 2 (TORUS_BOUND).\n");
        exit(EXIT_FAILURE);
      }
      break;
    case 'q':
      lattice_i = atoi(optarg);
      switch (lattice_i) {
      case 0:
        lattice = VORONOI_LATTICE;
        lattice_c = 'v';
        break;
      case 1:
        lattice = DIAGONAL_LATTICE;
        lattice_c = 'd';
        break;
      case 2:
        lattice = VORONOI_HYPERUNIFORM_LATTICE;
        lattice_c = 'h';
        break;
      case 3:
        lattice = TRIANGULAR_LATTICE;
        lattice_c = 't';
        break;
      case 4:
        lattice = SQUARE_LATTICE;
        lattice_c = 's';
        break;
      default:
        printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 "
               "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
        exit(EXIT_FAILURE);
      }
      break;
    case 'd':
      save_damage = true;
      break;
    case 'V':
      use_voltage_boundaries = true;
      break;
    case 'D':
      use_dual = true;
      dual_c = 'd';
      break;
    case 'c':
      save_cluster_dist = true;
      break;
    case 'o':
      save_data = true;
      break;
    case 'N':
      save_network = true;
      break;
    case 's':
      save_crit_stress = true;
      break;
    case 'r':
      save_conductivity = true;
      break;
    case 'E':
      save_energy = true;
      break;
    case 'T':
      save_threshold = true;
      break;
    case 'C':
      save_current_load = true;
      break;
    case 'z':
      save_correlations = true;
      break;
    default: /* '?' */
      exit(EXIT_FAILURE);
    }
  }

  char boundc;
  if (use_voltage_boundaries) {
    boundc = 'v';
  } else {
    boundc = 'c';
  }

  double *dd_correlations; // damage-damage
  double *dc_correlations; // damage-crack
  double *db_correlations; // damage-backbone
  double *ds_correlations; // damage-stress
  double *dA_correlations; // damage-avalanche
  double *cc_correlations; // crack-crack
  double *cb_correlations; // crack-backbone
  double *cs_correlations; // crack-stress
  double *cA_correlations; // crack-avalanche
  double *bb_correlations; // backbone-backbone
  double *bs_correlations; // backbone-stress
  double *bA_correlations; // backbone-avalanche
  double *ss_correlations; // stress-stress
  double *AA_correlations; // avalanche-avalanche
  double *DD_correlations; // after-crack damage-damage
  double *DC_correlations; // after-crack damage-crack
  double *DB_correlations; // after-crack damage-backbone
  double *CC_correlations; // after-crack crack-crack
  double *CB_correlations; // after-crack crack-backbone
  double *BB_correlations; // after-crack backbone-backbone
  double *fftw_forward_in;
  fftw_complex *fftw_forward_out;
  fftw_complex *fftw_reverse_in;
  double *fftw_reverse_out;
  fftw_plan forward_plan;
  fftw_plan reverse_plan;
  uint64_t N_averaged = 0;
  double mean_D = 0;
  double mean_A = 0;
  double mean_B = 0;
  double mean_C = 0;
  double mean_d = 0;
  double mean_c = 0;
  double mean_b = 0;
  double mean_s = 0;
  char *correlations_filename;
  if (save_correlations) {
    assert(lattice == DIAGONAL_LATTICE);
    dd_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    dc_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    db_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    ds_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    dA_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    cc_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    cb_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    cs_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    cA_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    bb_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    bs_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    bA_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    ss_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    AA_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    DD_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    DC_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    DB_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    CC_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    CB_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    BB_correlations = (double *)calloc(pow(L, 2), sizeof(double));
    fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
    fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
    fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
    fftw_reverse_out = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
    forward_plan = fftw_plan_dft_r2c_2d(L, L, fftw_forward_in, fftw_forward_out, 0);
    reverse_plan = fftw_plan_dft_c2r_2d(L, L, fftw_reverse_in, fftw_reverse_out, 0);

    correlations_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(correlations_filename, filename_len, "corr_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);

    FILE *correlations_out = fopen(correlations_filename, "rb");

    if (correlations_out != NULL) {
      fread(&N_averaged, sizeof(uint64_t), 1, correlations_out);
      fread(&mean_d, sizeof(double), 1, correlations_out);
      fread(&mean_c, sizeof(double), 1, correlations_out);
      fread(&mean_b, sizeof(double), 1, correlations_out);
      fread(&mean_s, sizeof(double), 1, correlations_out);
      fread(&mean_A, sizeof(double), 1, correlations_out);
      fread(&mean_D, sizeof(double), 1, correlations_out);
      fread(&mean_C, sizeof(double), 1, correlations_out);
      fread(&mean_B, sizeof(double), 1, correlations_out);
      fread(dd_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(dc_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(db_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(ds_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(dA_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(cc_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(cb_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(cs_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(cA_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(bb_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(bs_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(bA_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(ss_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(AA_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(DD_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(DC_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(DB_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(CC_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(CB_correlations, sizeof(double), pow(L, 2), correlations_out);
      fread(BB_correlations, sizeof(double), pow(L, 2), correlations_out);
      fclose(correlations_out);
    }
  }
  
  FILE *data_out;

  if (save_data) {
    char *data_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
    data_out = fopen(data_filename, "a");
    free(data_filename);
  }

  uint_t max_verts, max_edges;

  // these are very liberal estimates
  max_verts = 4 * pow(L, 2);
  max_edges = 4 * pow(L, 2);

  if (max_verts > CINT_MAX) {
    exit(EXIT_FAILURE);
  }

  // define arrays for saving cluster and avalanche distributions
  uint32_t *cluster_size_dist;
  uint32_t *avalanche_size_dist;
  char *c_filename;
  char *a_filename;
  if (save_cluster_dist) {
    cluster_size_dist = (uint32_t *)calloc(max_verts, sizeof(uint32_t));
    avalanche_size_dist = (uint32_t *)calloc(max_edges, sizeof(uint32_t));

    c_filename = (char *)malloc(filename_len * sizeof(char));
    a_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
    snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);

    FILE *cluster_out = fopen(c_filename, "rb");
    FILE *avalanche_out = fopen(a_filename, "rb");

    if (cluster_out != NULL) {
      fread(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
      fclose(cluster_out);
    }
    if (avalanche_out != NULL) {
      fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
      fclose(avalanche_out);
    }
  }

  long double *crit_stress;
  if (save_crit_stress) {
    crit_stress = (long double *)malloc(N * sizeof(long double));
  }

  double *conductivity;
  if (save_conductivity) {
    conductivity = (double *)malloc(N * sizeof(double));
  }

  // define arrays for saving damage distributions
  uint32_t *damage;
  char *d_filename;
  if (save_damage) {
    damage = (uint32_t *)calloc(max_edges, sizeof(uint32_t));

    d_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);

    FILE *damage_out = fopen(d_filename, "rb");

    if (damage_out != NULL) {
      fread(damage, sizeof(uint32_t), max_edges, damage_out);
      fclose(damage_out);
    }
  }

  long double *energy;
  if (save_energy) {
    energy = (long double *)malloc(N * sizeof(long double));
  }

  long double *thresholds;
  if (save_threshold) {
    thresholds = (long double *)malloc(N * sizeof(long double));
  }

  long double *loads;
  if (save_current_load) {
    loads = (long double *)malloc(N * sizeof(long double));
  }

  // start cholmod
  cholmod_common c;
  CHOL_F(start)(&c);

  /* if we use voltage boundary conditions, the laplacian matrix is positive
   * definite and we can use a supernodal LL decomposition.  otherwise we need
   * to use the simplicial LDL decomposition
   */
  if (use_voltage_boundaries) {
    //(&c)->supernodal = CHOLMOD_SUPERNODAL;
    (&c)->supernodal = CHOLMOD_SIMPLICIAL;
  } else {
    (&c)->supernodal = CHOLMOD_SIMPLICIAL;
  }

  printf("\n");
  for (uint32_t i = 0; i < N; i++) {
    printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1,
           N);

    graph_t *g = graph_create(lattice, boundary, L, use_dual);
    net_t *net =
        net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
    net_t *tmp_net = net_copy(net, &c);
    data_t *data = net_fracture(tmp_net, &c, cutoff);

    uint_t max_pos = 0;
    long double max_val = 0;

    double cond0;
    {
      double *tmp_voltages = net_voltages(net, &c);
      cond0 = net_conductivity(net, tmp_voltages);
      free(tmp_voltages);
    }

    for (uint_t j = 0; j < data->num_broken; j++) {
      long double val = data->extern_field[j];

      if (val > max_val) {
        max_pos = j;
        max_val = val;
      }
    }

    uint_t av_size = 0;
    long double cur_val = 0;

    for (uint_t j = 0; j < max_pos; j++) {
      uint_t next_broken = data->break_list[j];

      break_edge(net, next_broken, &c);

      long double val = data->extern_field[j];
      if (save_cluster_dist) {
        if (val < cur_val) {
          av_size++;
        }

        if (val > cur_val) {
          avalanche_size_dist[av_size]++;
          av_size = 0;
          cur_val = val;
        }
      }
    }

    if (save_correlations) {
      uint32_t damage1 = 0;
      for (uint32_t j = 0; j < g->ne; j++) {
        if (tmp_net->fuses[j]) {
          fftw_forward_in[j] = 1.0;
          damage1++;
        } else {
          fftw_forward_in[j] = 0.0;
        }
      }

      fftw_execute(forward_plan);
      fftw_complex *D_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
      memcpy(D_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));

      dll_t *cycle = find_cycles(g, tmp_net->fuses);
      components_t *comp = get_clusters(tmp_net);

      uint32_t in_crack1 = 0;
      for (uint32_t j = 0; j < g->ne; j++) {
        if (tmp_net->fuses[j] && comp->labels[g->dev[2 * cycle->e]] == comp->labels[g->dev[2 * j]]) {
          fftw_forward_in[j] = 1.0;
          in_crack1++;
        } else {
          fftw_forward_in[j] = 0.0;
        }
      }

      fftw_execute(forward_plan);
      fftw_complex *C_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
      memcpy(C_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));

      for (uint32_t j = 0; j < g->ne; j++) {
        fftw_forward_in[j] = 0.0;
      }

      uint32_t in_backbone1 = 0;
      dll_t *tmp_cycle = cycle;
      while (tmp_cycle != NULL) {
        fftw_forward_in[tmp_cycle->e] = 1.0;
        in_backbone1++;
        tmp_cycle = tmp_cycle->right;
      }

      fftw_execute(forward_plan);
      fftw_complex *B_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex));
      memcpy(B_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));

      uint32_t in_avalanche = 0;
      for (uint32_t j = 0; j < g->ne; j++) {
        if (tmp_net->fuses[j] && !(net->fuses[j])) {
          fftw_forward_in[j] = 1.0;
        } else {
          fftw_forward_in[j] = 0.0;
        }
      }

      fftw_execute(forward_plan);
      fftw_complex *A_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex));
      memcpy(A_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));

      uint32_t damage2 = 0;
      for (uint32_t j = 0; j < g->ne; j++) {
        if (net->fuses[j]) {
          fftw_forward_in[j] = 1.0;
          damage2++;
        } else {
          fftw_forward_in[j] = 0.0;
        }
      }

      fftw_execute(forward_plan);
      fftw_complex *d_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
      memcpy(d_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));

      uint32_t in_crack2 = 0;
      for (uint32_t j = 0; j < g->ne; j++) {
        if (net->fuses[j] && comp->labels[2 * g->dev[cycle->e]] == comp->labels[g->dev[2 * j]]) {
          fftw_forward_in[j] = 1.0;
          in_crack2++;
        } else {
          fftw_forward_in[j] = 0.0;
        }
      }

      fftw_execute(forward_plan);
      fftw_complex *c_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex));
      memcpy(c_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));

      for (uint32_t j = 0; j < g->ne; j++) {
        fftw_forward_in[j] = 0.0;
      }

      uint32_t in_backbone2 = 0;
      tmp_cycle = cycle;
      while (tmp_cycle != NULL) {
        if (net->fuses[tmp_cycle->e]) {
          fftw_forward_in[tmp_cycle->e] = 1.0;
          in_backbone2++;
        }
        tmp_cycle = tmp_cycle->right;
      }

      fftw_execute(forward_plan);
      fftw_complex *b_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex));
      memcpy(b_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));

      graph_components_free(comp);
      while (cycle != NULL) {
        dll_t *old = cycle;
        cycle = cycle->right;
        free(old);
      }

      double *tmp_voltage = net_voltages(net, &c);
      double *tmp_current = net_currents(net, tmp_voltage, &c);
      free(tmp_voltage);

      double t_stress = 0;
      for (uint32_t j = 0; j < g->ne; j++) {
        fftw_forward_in[j] = tmp_current[j];
        t_stress += tmp_current[j];
      }

      free(tmp_current);

      fftw_execute(forward_plan);
      fftw_complex *s_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
      memcpy(s_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));

      mean_D = (double)damage1 / (1.0 + N_averaged) + (double)N_averaged * mean_D / (N_averaged + 1.0);
      mean_C = (double)in_crack1 / (1.0 + N_averaged) + (double)N_averaged * mean_C / (N_averaged + 1.0);
      mean_B = (double)in_backbone1 / (1.0 + N_averaged) + (double)N_averaged * mean_B / (N_averaged + 1.0);
      mean_A = (double)in_avalanche / (1.0 + N_averaged) + (double)N_averaged * mean_A / (N_averaged + 1.0);
      mean_d = (double)damage2 / (1.0 + N_averaged) + (double)N_averaged * mean_d / (N_averaged + 1.0);
      mean_c = (double)in_crack2 / (1.0 + N_averaged) + (double)N_averaged * mean_c / (N_averaged + 1.0);
      mean_b = (double)in_backbone2 / (1.0 + N_averaged) + (double)N_averaged * mean_b / (N_averaged + 1.0);
      mean_s = (double)t_stress / (1.0 + N_averaged) + (double)N_averaged * mean_s / (N_averaged + 1.0);

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = D_transform[j][0] * D_transform[j][0] + D_transform[j][1] * D_transform[j][1];
        fftw_reverse_in[j][1] = 0.0;
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        DD_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DD_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = D_transform[j][0] * C_transform[j][0] + D_transform[j][1] * C_transform[j][1];
        fftw_reverse_in[j][1] = D_transform[j][0] * C_transform[j][1] - D_transform[j][1] * C_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        DC_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DC_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = D_transform[j][0] * B_transform[j][0] + D_transform[j][1] * B_transform[j][1];
        fftw_reverse_in[j][1] = D_transform[j][0] * B_transform[j][1] - D_transform[j][1] * B_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        DB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DB_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = C_transform[j][0] * C_transform[j][0] + C_transform[j][1] * C_transform[j][1];
        fftw_reverse_in[j][1] = 0.0;
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        CC_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * CC_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = C_transform[j][0] * B_transform[j][0] + C_transform[j][1] * B_transform[j][1];
        fftw_reverse_in[j][1] = C_transform[j][0] * B_transform[j][1] - C_transform[j][1] * B_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        CB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * CB_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = B_transform[j][0] * B_transform[j][0] + B_transform[j][1] * B_transform[j][1];
        fftw_reverse_in[j][1] = 0.0;
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        BB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * BB_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = A_transform[j][0] * A_transform[j][0] + A_transform[j][1] * A_transform[j][1];
        fftw_reverse_in[j][1] = 0.0;
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        AA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * AA_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = d_transform[j][0] * d_transform[j][0] + d_transform[j][1] * d_transform[j][1];
        fftw_reverse_in[j][1] = 0.0;
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        dd_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dd_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = d_transform[j][0] * c_transform[j][0] + d_transform[j][1] * c_transform[j][1];
        fftw_reverse_in[j][1] = d_transform[j][0] * c_transform[j][1] - d_transform[j][1] * c_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        dc_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dc_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = d_transform[j][0] * s_transform[j][0] + d_transform[j][1] * s_transform[j][1];
        fftw_reverse_in[j][1] = d_transform[j][0] * s_transform[j][1] - d_transform[j][1] * s_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        ds_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ds_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = d_transform[j][0] * b_transform[j][0] + d_transform[j][1] * b_transform[j][1];
        fftw_reverse_in[j][1] = d_transform[j][0] * b_transform[j][1] - d_transform[j][1] * b_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        db_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * db_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = d_transform[j][0] * A_transform[j][0] + d_transform[j][1] * A_transform[j][1];
        fftw_reverse_in[j][1] = d_transform[j][0] * A_transform[j][1] - d_transform[j][1] * A_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        dA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dA_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = c_transform[j][0] * c_transform[j][0] + c_transform[j][1] * c_transform[j][1];
        fftw_reverse_in[j][1] = 0.0;
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        cc_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cc_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = c_transform[j][0] * s_transform[j][0] + c_transform[j][1] * s_transform[j][1];
        fftw_reverse_in[j][1] = c_transform[j][0] * s_transform[j][1] - c_transform[j][1] * s_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        cs_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cs_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = c_transform[j][0] * b_transform[j][0] + c_transform[j][1] * b_transform[j][1];
        fftw_reverse_in[j][1] = c_transform[j][0] * b_transform[j][1] - c_transform[j][1] * b_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        cb_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cb_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = c_transform[j][0] * A_transform[j][0] + c_transform[j][1] * A_transform[j][1];
        fftw_reverse_in[j][1] = c_transform[j][0] * A_transform[j][1] - c_transform[j][1] * A_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        cA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cA_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = b_transform[j][0] * b_transform[j][0] + b_transform[j][1] * b_transform[j][1];
        fftw_reverse_in[j][1] = 0.0;
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        bb_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bb_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = b_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1];
        fftw_reverse_in[j][1] = b_transform[j][0] * s_transform[j][1] - s_transform[j][1] * s_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        bs_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bs_correlations[j] / (N_averaged + 1.0);
      }

      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = A_transform[j][0] * b_transform[j][0] + A_transform[j][1] * b_transform[j][1];
        fftw_reverse_in[j][1] = A_transform[j][0] * b_transform[j][1] - A_transform[j][1] * b_transform[j][0];
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        bA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bA_correlations[j] / (N_averaged + 1.0);
      }


      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
        fftw_reverse_in[j][0] = s_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1];
        fftw_reverse_in[j][1] = 0.0;
      }

      fftw_execute(reverse_plan);

      for (uint32_t j = 0; j < g->ne; j++) {
        ss_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ss_correlations[j] / (N_averaged + 1.0);
      }

      free(D_transform);
      free(C_transform);
      free(A_transform);
      free(B_transform);
      free(b_transform);
      free(d_transform);
      free(c_transform);
      free(s_transform);
      N_averaged++;
    }

    net_free(tmp_net, &c);

    if (save_crit_stress) {
      crit_stress[i] = data->extern_field[max_pos];
    }

    if (save_conductivity) {
      if (max_pos > 0) {
        conductivity[i] = data->conductivity[max_pos - 1];
      } else {
        conductivity[i] = cond0;
      }
    }

    if (save_damage) {
      uint_t would_break = 0;
      double *tmp_voltage = net_voltages(net, &c);
      double *tmp_current = net_currents(net, tmp_voltage, &c);
      free(tmp_voltage);
      for (uint_t j = 0; j < g->ne; j++) {
        bool broken = net->fuses[j];
        bool under_thres =
            net->thres[j] < net->thres[data->break_list[max_pos]];
        bool zero_field = fabs(tmp_current[j]) < cutoff;
        if (!broken && under_thres && zero_field) {
          break_edge(net, j, &c);
        }
      }
      damage[net->num_broken]++;
      free(tmp_current);
    }

    if (save_energy) {
      long double tmp_energy = 0;
      if (max_pos > 0) {
        long double sigma1 = data->extern_field[0];
        double cond1 = cond0;
        for (uint_t j = 0; j < max_pos - 1; j++) {
          long double sigma2 = data->extern_field[j + 1];
          double cond2 = data->conductivity[j];
          if (sigma2 > sigma1) {
            tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1;
            sigma1 = sigma2;
            cond1 = cond2;
          }
        }
      }
      energy[i] = tmp_energy;
    }

    if (save_threshold) {
      thresholds[i] = net->thres[data->break_list[max_pos]];
    }

    if (save_current_load) {
      loads[i] =
          data->extern_field[max_pos] / net->thres[data->break_list[max_pos]];
    }

    if (save_data) {
      for (uint_t j = 0; j < data->num_broken; j++) {
        fprintf(data_out, "%u %Lg %g ", data->break_list[j],
                data->extern_field[j], data->conductivity[j]);
      }
      fprintf(data_out, "\n");
    }

    data_free(data);
    if (save_network) {
      FILE *net_out = fopen("network.txt", "w");
      for (uint_t j = 0; j < g->nv; j++) {
        fprintf(net_out, "%f %f ", g->vx[2 * j], g->vx[2 * j + 1]);
      }
      fprintf(net_out, "\n");
      for (uint_t j = 0; j < g->ne; j++) {
        fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]);
      }
      fprintf(net_out, "\n");
      for (uint_t j = 0; j < g->dnv; j++) {
        fprintf(net_out, "%f %f ", g->dvx[2 * j], g->dvx[2 * j + 1]);
      }
      fprintf(net_out, "\n");
      for (uint_t j = 0; j < g->ne; j++) {
        fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]);
      }
      fprintf(net_out, "\n");
      for (uint_t j = 0; j < g->ne; j++) {
        fprintf(net_out, "%d ", net->fuses[j]);
      }
      fclose(net_out);
    }

    if (save_cluster_dist) {
      uint_t *tmp_cluster_dist = get_cluster_dist(net);
      for (uint_t j = 0; j < g->dnv; j++) {
        cluster_size_dist[j] += tmp_cluster_dist[j];
      }
      free(tmp_cluster_dist);
    }

    net_free(net, &c);
    graph_free(g);
  }

  printf("\033[F\033[JFRACTURE: COMPLETE\n");

  if (save_cluster_dist) {
    FILE *cluster_out = fopen(c_filename, "wb");
    FILE *avalanche_out = fopen(a_filename, "wb");

    fwrite(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
    fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);

    fclose(cluster_out);
    fclose(avalanche_out);

    free(c_filename);
    free(a_filename);
    free(cluster_size_dist);
    free(avalanche_size_dist);
  }

  if (save_conductivity) {
    char *cond_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
    FILE *cond_file = fopen(cond_filename, "ab");
    fwrite(conductivity, sizeof(double), N, cond_file);
    fclose(cond_file);
    free(cond_filename);
    free(conductivity);
  }

  if (save_energy) {
    char *tough_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
    FILE *tough_file = fopen(tough_filename, "ab");
    fwrite(energy, sizeof(long double), N, tough_file);
    fclose(tough_file);
    free(tough_filename);
    free(energy);
  }

  if (save_threshold) {
    char *thres_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
    FILE *thres_file = fopen(thres_filename, "ab");
    fwrite(thresholds, sizeof(long double), N, thres_file);
    fclose(thres_file);
    free(thres_filename);
    free(thresholds);
  }

  if (save_current_load) {
    char *load_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
    FILE *load_file = fopen(load_filename, "ab");
    fwrite(loads, sizeof(long double), N, load_file);
    fclose(load_file);
    free(load_filename);
    free(loads);
  }

  if (save_damage) {
    FILE *hdam_file = fopen(d_filename, "wb");
    fwrite(damage, sizeof(uint32_t), max_edges, hdam_file);
    fclose(hdam_file);
    free(d_filename);
    free(damage);
  }

  if (save_data) {
    fclose(data_out);
  }

  if (save_crit_stress) {
    char *str_filename = (char *)malloc(filename_len * sizeof(char));
    snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat",
             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
    FILE *str_file = fopen(str_filename, "ab");
    fwrite(crit_stress, sizeof(long double), N, str_file);
    fclose(str_file);
    free(str_filename);
    free(crit_stress);
  }

  if (save_correlations) {
    FILE *correlations_out = fopen(correlations_filename, "wb");
    fwrite(&N_averaged, sizeof(uint64_t), 1, correlations_out);
    fwrite(&mean_d, sizeof(double), 1, correlations_out);
    fwrite(&mean_c, sizeof(double), 1, correlations_out);
    fwrite(&mean_b, sizeof(double), 1, correlations_out);
    fwrite(&mean_s, sizeof(double), 1, correlations_out);
    fwrite(&mean_A, sizeof(double), 1, correlations_out);
    fwrite(&mean_D, sizeof(double), 1, correlations_out);
    fwrite(&mean_C, sizeof(double), 1, correlations_out);
    fwrite(&mean_B, sizeof(double), 1, correlations_out);
    fwrite(dd_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(dc_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(db_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(ds_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(dA_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(cc_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(cb_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(cs_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(cA_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(bb_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(bs_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(bA_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(ss_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(AA_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(DD_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(DC_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(DB_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(CC_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(CB_correlations, sizeof(double), pow(L, 2), correlations_out);
    fwrite(BB_correlations, sizeof(double), pow(L, 2), correlations_out);
    fclose(correlations_out);

    free(dd_correlations);
    free(dc_correlations);
    free(db_correlations);
    free(ds_correlations);
    free(dA_correlations);
    free(cc_correlations);
    free(cb_correlations);
    free(cs_correlations);
    free(cA_correlations);
    free(bb_correlations);
    free(bs_correlations);
    free(bA_correlations);
    free(ss_correlations);
    free(AA_correlations);
    free(DD_correlations);
    free(DC_correlations);
    free(DB_correlations);
    free(CC_correlations);
    free(CB_correlations);
    free(BB_correlations);
    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);
    free(correlations_filename);
  }

  fftw_cleanup();

  CHOL_F(finish)(&c);

  return 0;
}