#include "network.hpp"
#include <iostream>

class nofuseException: public std::exception
{
  virtual const char* what() const throw()
  {
    return "No valid fuse was available to break.";
  }
} nofuseex;

network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {}

network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {}

void network::set_thresholds(double beta, std::mt19937& rng) {
  if (beta == 0.0) {
    /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
    for (long double& threshold : thresholds) {
      threshold = 1.0;
    }
  } else {
    std::uniform_real_distribution<long double> d(0.0, 1.0);

    for (long double& threshold : thresholds) {
      threshold = std::numeric_limits<long double>::lowest();

      while (threshold == std::numeric_limits<long double>::lowest()) {
        threshold = logl(d(rng)) / (long double)beta;
      }
    }
  }
}

void network::fracture(hooks& m, bool one_axis) {
  m.pre_fracture(*this);

  while (true) {
    current_info c = this->get_current_info();

    double min_cond = 1.0 / G.edges.size();

    if (c.conductivity[0] < min_cond && c.conductivity[1] < min_cond) {
      break;
    }

    if (one_axis && (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond)) {
      break;
    }

    unsigned max_pos = UINT_MAX;
    long double max_val = std::numeric_limits<long double>::lowest();

    for (unsigned i = 0; i < G.edges.size(); i++) {
      if (!fuses[i] && c.currents[i] > min_cond) {
        long double val = logl(c.currents[i]) - thresholds[i];

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

    if (max_pos == UINT_MAX)  {
      throw nofuseex;
    }

    this->break_edge(max_pos);
    m.bond_broken(*this, c, max_pos);
  }

  m.post_fracture(*this);
}


fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) :
  network(G), ohm_x(G, 0, c), ohm_y(G, 1, c) {}

fuse_network::fuse_network(const fuse_network& n) :
  network(n), ohm_x(n.ohm_x), ohm_y(n.ohm_y) {
}

void fuse_network::break_edge(unsigned e, bool unbreak) {
  fuses[e] = !unbreak;
  ohm_x.break_edge(e, unbreak);
  ohm_y.break_edge(e, unbreak);
}

current_info fuse_network::get_current_info() {
  current_info cx = ohm_x.solve(fuses);
  current_info cy = ohm_y.solve(fuses);

  bool done_x = cx.conductivity[0] < 1.0 / G.edges.size();
  bool done_y = cy.conductivity[1] < 1.0 / G.edges.size();

  current_info ctot;
  ctot.currents.resize(G.edges.size());
  ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]}; 

  if (done_x && !done_y) {
    for (unsigned i = 0; i < G.edges.size(); i++) {
      ctot.currents[i] = fabs(weight * cy.currents[i] / cy.conductivity[1]);
    }
  } else if (done_y && !done_x) {
    for (unsigned i = 0; i < G.edges.size(); i++) {
      ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0]);
    }
  } else if (!done_x && !done_y) {
    for (unsigned i = 0; i < G.edges.size(); i++) {
      ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0] +
                                    weight * cy.currents[i] / cy.conductivity[1]);
    }
  }

  return ctot;
}


elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight) :
  network(G), weight(weight),
  hook_x(G, 0, c), hook_y(G, 1, c) {}

elastic_network::elastic_network(const elastic_network& n) :
  network(n), weight(n.weight), hook_x(n.hook_x), hook_y(n.hook_y) {}

void elastic_network::break_edge(unsigned e, bool unbreak) {
  fuses[e] = !unbreak;
  hook_x.break_edge(e, unbreak);
  hook_y.break_edge(e, unbreak);
}

current_info elastic_network::get_current_info() {
  current_info cx = hook_x.solve(fuses);
  current_info cy = hook_y.solve(fuses);

  bool done_x = cx.conductivity[0] < 1.0 / G.edges.size();
  bool done_y = cy.conductivity[1] < 1.0 / G.edges.size();

  current_info ctot;
  ctot.currents.resize(G.edges.size());
  ctot.conductivity[0] = cx.conductivity[0]; 
  ctot.conductivity[1] = cy.conductivity[1]; 

  if (done_x && !done_y) {
    for (unsigned i = 0; i < G.edges.size(); i++) {
      ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity[1];
    }
  } else if (done_y && !done_x) {
    for (unsigned i = 0; i < G.edges.size(); i++) {
      ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity[0];
    }
  } else if (!done_x && !done_y) {
    for (unsigned i = 0; i < G.edges.size(); i++) {
      ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity[0], 2) +
                              pow(weight * cy.currents[i] / cy.conductivity[1], 2));
    }
  }

  return ctot;
}


percolation_network::percolation_network(const graph& G, cholmod_common* c) :
  network(G), px(G, 0, c), py(G, 1, c) {}

percolation_network::percolation_network(const percolation_network& n) :
  network(n), px(n.px), py(n.py) {}

current_info percolation_network::get_current_info() {
  current_info ctot;
  ctot.currents.resize(G.edges.size(), 0);

  current_info cx = px.solve(fuses);
  current_info cy = py.solve(fuses);

  ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]};

  double min_cur = 1.0 / G.edges.size();

  for (unsigned i = 0; i < G.edges.size(); i++) {
    if (fabs(cx.currents[i]) > min_cur || fabs(cy.currents[i]) > min_cur) {
      ctot.currents[i] = 1.0;
    }
  }
  
  return ctot;
}

void percolation_network::break_edge(unsigned e, bool unbreak) {
  fuses[e] = !unbreak;
  px.break_edge(e, unbreak);
  py.break_edge(e, unbreak);
}