#include "network.hpp"

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, cholmod_common* c) :
  G(G), bG(G.vertices.size()), dG(G.dual_vertices.size()),
  rank(G.dual_vertices.size()), parent(G.dual_vertices.size()), ds(&rank[0], &parent[0]),
  fuses(G.edges.size()), backbone(G.edges.size()), thresholds(G.edges.size()),
  px(G, 0, c), py(G, 1, c) {
  for (unsigned i = 0; i < G.edges.size(); i++) {
    boost::add_edge(G.edges[i].v[0], G.edges[i].v[1], {i}, bG);
  }
  initialize_incremental_components(dG, ds);
}

network::network(const network& n) : G(n.G), bG(n.bG), dG(n.dG), rank(n.rank), parent(n.parent), ds(&rank[0], &parent[0]), fuses(n.fuses), backbone(n.backbone), thresholds(n.thresholds), px(n.px), py(n.py) {
  initialize_incremental_components(dG, ds);
}

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) {
    double min_cond = 1.0 / G.edges.size();

    current_info c = this->get_current_info();
    if (px.sol.conductivity[0] > min_cond) {
      this->update_backbone(px.sol.currents);
    } else {
      this->update_backbone(py.sol.currents);
    }

    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 (!backbone[i]) {
        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);
}

void network::update_backbone(const std::vector<double>& c) {

        bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
        bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
  class get_cycle_edges : public boost::default_dfs_visitor {
    public:
      std::set<unsigned> tE;
      std::list<unsigned>& bE;

      get_cycle_edges(std::list<unsigned>& bE) : bE(bE) {};

      void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
        tE.insert(g[e].index);
      }

      void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
        if (tE.find(g[e].index) == tE.end()) {
          bE.push_back(g[e].index);
        }
      }
  };

  class find_cycle : public boost::default_dfs_visitor {
    public:
      const graph& G;
      std::array<unsigned, 2>& E;
      unsigned end;
      struct done{};

      find_cycle(const graph& G, std::array<unsigned, 2>& E, unsigned end) : G(G), E(E), end(end) {}

      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
        if (v == end) {
          throw done{};
        }
      }

      void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
        unsigned e = g[ed].index;
        if (G.dual_edges[e].crossings[0]) E[0]++;
        if (G.dual_edges[e].crossings[1]) E[1]++;
      }

      void finish_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
        unsigned e = g[ed].index;
        if (G.dual_edges[e].crossings[0]) E[0]--;
        if (G.dual_edges[e].crossings[1]) E[1]--;
      }
  };

  class get_cluster_edges : public boost::default_dfs_visitor {
    public:
      std::set<unsigned>& E;

      get_cluster_edges(std::set<unsigned>& E) : E(E) {}

      void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
        unsigned e = g[ed].index;
        E.insert(e);
      }
  };

  for (unsigned i = 0; i < G.edges.size(); i++) {
    if ((!backbone[i]) && boost::same_component(G.dual_edges[i].v[0], G.dual_edges[i].v[1], ds)) {
      auto it_search = seen_pairs.find({G.dual_edges[i].v[0], G.dual_edges[i].v[1]});
      if (it_search == seen_pairs.end()) {
      std::list<unsigned> cedges = {};
      get_cycle_edges vis(cedges);
      std::vector<boost::default_color_type> cm1(G.dual_vertices.size());
      boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis, boost::make_iterator_property_map(cm1.begin(), boost::get(boost::vertex_index, dG), cm1[0]));

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

      std::array<unsigned, 2> cycle = {0, 0};
      if (G.dual_edges[i].crossings[0]) cycle[0]++;
      if (G.dual_edges[i].crossings[1]) cycle[1]++;
      find_cycle vis2(G, cycle, G.dual_edges[i].v[1]);
      std::vector<boost::default_color_type> cm2(G.dual_vertices.size());
      try {
        boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis2, boost::make_iterator_property_map(cm2.begin(), boost::get(boost::vertex_index, dG), cm2[0]));
      } catch(find_cycle::done const&) {
        std::list<std::array<unsigned, 2>> other_cycles = {cycle};
        for (unsigned e : cedges) {
          std::array<unsigned, 2> tcycle = {0, 0};
          if (G.dual_edges[i].crossings[0]) tcycle[0]++;
          if (G.dual_edges[i].crossings[1]) tcycle[1]++;
          if (G.dual_edges[e].crossings[0]) tcycle[0]++;
          if (G.dual_edges[e].crossings[1]) tcycle[1]++;
          find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]);
          std::vector<boost::default_color_type> cm3(G.dual_vertices.size());
          try {
            boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis3, boost::make_iterator_property_map(cm3.begin(), boost::get(boost::vertex_index, dG), cm3[0]));
          } catch(find_cycle::done const&) {
          }
          find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]);
          std::vector<boost::default_color_type> cm4(G.dual_vertices.size());
          try {
            boost::depth_first_visit(dG, G.dual_edges[i].v[1], vis4, boost::make_iterator_property_map(cm4.begin(), boost::get(boost::vertex_index, dG), cm4[0]));
          } catch(find_cycle::done const&) {
          }
          other_cycles.push_back(tcycle);
        }

        if (std::any_of(other_cycles.begin(), other_cycles.end(), [done_x, done_y](std::array<unsigned, 2> c){return (c[0] % 2 == 0 && c[1] % 2 == 0) || ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y));})) {
          seen_pairs[{G.edges[i].v[0], G.edges[i].v[1]}] = true;
          backbone[i] = true;
          boost::remove_edge(G.edges[i].v[0], G.edges[i].v[1], bG);

          for (unsigned j = 0; j < 2; j++) {
            std::set<unsigned> component_edges = {};
            get_cluster_edges visc(component_edges);
            std::vector<boost::default_color_type> cm2(G.vertices.size());
            boost::depth_first_visit(bG, G.edges[i].v[j], visc, boost::make_iterator_property_map(cm2.begin(), boost::get(boost::vertex_index, bG), cm2[0]));

            std::array<double, 2> total_cur = {0.0, 0.0};
            for (unsigned e : component_edges) {
              if (G.edges[e].crossings[0]) {
                if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
                  total_cur[0] += c[e];
                } else {
                  total_cur[0] -= c[e];
                }
              }
              if (G.edges[e].crossings[1]) {
                if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
                  total_cur[1] += c[e];
                } else {
                  total_cur[1] -= c[e];
                }
              }
            }


            if (fabs(total_cur[0]) < 1.0 / G.edges.size() && fabs(total_cur[1]) < 1.0 / G.edges.size()) {
              for (unsigned e : component_edges) {
                backbone[e] = true;
                boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
              }
            }
          }
        }
      }
      for (unsigned e : cedges) {
        boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
      }
    }
    }
  }

  class get_tie_edges : public boost::default_dfs_visitor {
    public:
      std::set<unsigned>& LE;
      std::set<unsigned>& RE;
      unsigned start_edge;
      bool initialized;
      bool done_left;
      unsigned disc;
      unsigned dr;

      get_tie_edges(std::set<unsigned>& LE, std::set<unsigned>& RE) : LE(LE), RE(RE) {
        initialized = false;
        done_left = false;
        disc = 0;
        dr = 0;
      }

      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
        if (!initialized && disc > 0) {
          start_edge = vd;
          initialized = true;
        }
        disc++;
      }

      void finish_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
        if (vd == start_edge) {
          done_left = true;
        }
      }

      void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
        unsigned e = g[ed].index;
        if (!done_left) {
          LE.insert(e);
        } else if (LE.find(e) == LE.end()) {
          RE.insert(e);
        }
      }


  };

  std::set<unsigned> bb_verts;

  for (unsigned i = 0; i < G.edges.size(); i++) {
    if (!backbone[i]) {
      bb_verts.insert(G.edges[i].v[0]);
      bb_verts.insert(G.edges[i].v[1]);
    }
  }

  for (unsigned v : bb_verts) {
    bool found_pair = false;
    unsigned d1, d2, p1, p2;
    for (unsigned i = 0; i < G.vertices[v].nd.size(); i++) {
      for (unsigned j = 2; j < G.vertices[v].nd.size() - 1; j++) {
        unsigned l = (i+j)%G.vertices[v].nd.size();
        if (boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[l], ds)) {
          unsigned il = i < l ? i : l;
          unsigned ig = i < l ? l : i;
          auto it_search = seen_pairs.find({G.vertices[v].nd[il], G.vertices[v].nd[ig]});
          if (it_search == seen_pairs.end()) {
          bool any_intervening1 = false;
          bool any_intervening_e_1 = false;
          bool any_intervening2 = false;
          bool any_intervening_e_2 = false;
          unsigned ie1 = 0;
          unsigned ie2 = 0;

          for (unsigned k = il + 1; k < ig; k++) {
            if (!boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[k], ds)) {
              any_intervening1 = true;
              break;
            }
          }
          for (unsigned k = (ig + 1); k%G.vertices[v].nd.size() != il; k++) {
            if (!boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[k%G.vertices[v].nd.size()], ds)) {
              any_intervening2 = true;
              break;
            }
          }
          if (any_intervening2 && !any_intervening1) {
            for (unsigned k = il + 1; k <= ig; k++) {
              if (!backbone[G.vertices[v].ne[k]]) {
                ie1++;
              }
            }
          }
          if (any_intervening1 && !any_intervening2) {
            for (unsigned k = ig + 1; k%G.vertices[v].nd.size() != il + 1; k++) {
              if (!backbone[G.vertices[v].ne[k%G.vertices[v].nd.size()]]) {
                ie2++;
              }
            }
          }
          if ((any_intervening1 && any_intervening2) || (any_intervening1 && ie2 > 1) || (any_intervening2 && ie1 > 1)) {
            found_pair = true;
            p1 = il;
            p2 = ig;
            d1 = G.vertices[v].nd[il];
            d2 = G.vertices[v].nd[ig];

      std::list<unsigned> cedges = {};
      get_cycle_edges vis(cedges);
      std::vector<boost::default_color_type> cm1(G.dual_vertices.size());
      boost::depth_first_visit(dG, d1, vis, boost::make_iterator_property_map(cm1.begin(), boost::get(boost::vertex_index, dG), cm1[0]));

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

      std::array<unsigned, 2> cycle = {0, 0};
      for (unsigned k = p1; k < p2; k++) {
        if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) cycle[0]++;
        if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) cycle[1]++;
      }
      find_cycle vis2(G, cycle, d2);
      std::vector<boost::default_color_type> cm2(G.dual_vertices.size());
      try {
        boost::depth_first_visit(dG, d1, vis2, boost::make_iterator_property_map(cm2.begin(), boost::get(boost::vertex_index, dG), cm2[0]));
      } catch(find_cycle::done const&) {
        std::list<std::array<unsigned, 2>> other_cycles = {cycle};
        for (unsigned e : cedges) {
          std::array<unsigned, 2> tcycle = {0, 0};
          for (unsigned k = p1; k < p2; k++) {
            if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) tcycle[0]++;
            if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) tcycle[1]++;
          }
          if (G.dual_edges[e].crossings[0]) tcycle[0]++;
          if (G.dual_edges[e].crossings[1]) tcycle[1]++;
          find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]);
          std::vector<boost::default_color_type> cm3(G.dual_vertices.size());
          try {
            boost::depth_first_visit(dG, d1, vis3, boost::make_iterator_property_map(cm3.begin(), boost::get(boost::vertex_index, dG), cm3[0]));
          } catch(find_cycle::done const&) {
          }
          find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]);
          std::vector<boost::default_color_type> cm4(G.dual_vertices.size());
          try {
            boost::depth_first_visit(dG, d2, vis4, boost::make_iterator_property_map(cm4.begin(), boost::get(boost::vertex_index, dG), cm4[0]));
          } catch(find_cycle::done const&) {
          }
          other_cycles.push_back(tcycle);
        }

        if (std::any_of(other_cycles.begin(), other_cycles.end(), [done_x, done_y](std::array<unsigned, 2> c){return ((c[0] % 2 == 0 && c[1] % 2 == 0) || (c[0] % 2 == 0 && done_x)) || (c[1] % 2 == 0 && done_y);})) {
          seen_pairs[{d1, d2}] = true;

      std::set<unsigned> left_edges = {};
      std::set<unsigned> right_edges = {};
      get_tie_edges visc(left_edges, right_edges);
      std::vector<boost::default_color_type> cm2(G.vertices.size());
      boost::depth_first_visit(bG, v, visc, boost::make_iterator_property_map(cm2.begin(), boost::get(boost::vertex_index, bG), cm2[0]));

            std::array<double, 2> total_cur_left = {0.0, 0.0};
            for (unsigned e : left_edges) {
              if (G.edges[e].crossings[0]) {
                if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
                  total_cur_left[0] += c[e];
                } else {
                  total_cur_left[0] -= c[e];
                }
              }
              if (G.edges[e].crossings[1]) {
                if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
                  total_cur_left[1] += c[e];
                } else {
                  total_cur_left[1] -= c[e];
                }
              }
            }

            std::array<double, 2> total_cur_right = {0.0, 0.0};
            for (unsigned e : right_edges) {
              if (G.edges[e].crossings[0]) {
                if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
                  total_cur_right[0] += c[e];
                } else {
                  total_cur_right[0] -= c[e];
                }
              }
              if (G.edges[e].crossings[1]) {
                if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
                  total_cur_right[1] += c[e];
                } else {
                  total_cur_right[1] -= c[e];
                }
              }
            }


            if (fabs(total_cur_left[0]) < 1.0 / G.edges.size() && fabs(total_cur_left[1]) < 1.0 / G.edges.size()) {
              for (unsigned e : left_edges) {
                backbone[e] = true;
                boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
              }
            }
            if (fabs(total_cur_right[0]) < 1.0 / G.edges.size() && fabs(total_cur_right[1]) < 1.0 / G.edges.size()) {
              for (unsigned e : right_edges) {
                backbone[e] = true;
                boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
              }
            }

      }
    }
      for (unsigned e : cedges) {
        boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
      }
  }
  }
        }
      }
    }
  }
}

void network::break_edge(unsigned e, bool unbreak) {
  fuses[e] = !unbreak;
  backbone[e] = !unbreak;
  boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
  boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
  ds.union_set(G.dual_edges[e].v[0], G.dual_edges[e].v[1]);
  px.break_edge(e, unbreak);
  py.break_edge(e, unbreak);
}


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

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

current_info fuse_network::get_current_info() {
  px.solve(fuses);
  py.solve(fuses);

  bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
  bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();

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

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

  return ctot;
}


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

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

current_info elastic_network::get_current_info() {
  px.solve(fuses);
  py.solve(fuses);

  bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
  bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();

  current_info ctot;
  ctot.currents.resize(G.edges.size());
  ctot.conductivity[0] = px.sol.conductivity[0]; 
  ctot.conductivity[1] = py.sol.conductivity[1]; 

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

  return ctot;
}


percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, c) {}

percolation_network::percolation_network(const percolation_network& n) : network(n) {}

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

  px.solve(fuses);
  py.solve(fuses);

  ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]};

  for (unsigned i = 0; i < G.edges.size(); i++) {
    if (fabs(px.sol.currents[i]) > CURRENT_CUTOFF || fabs(py.sol.currents[i]) > CURRENT_CUTOFF) {
      ctot.currents[i] = 1.0;
    }
  }
  
  return ctot;
}