refactored much of the fracture library and fixed some percolation measurements

6 files changed, 359 insertions, 379 deletions

diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
index 49a9602..2592f3d 100644
--- a/lib/CMakeLists.txt
+++ b/lib/CMakeLists.txt
@@ -1,4 +1,4 @@
 
-add_library(frac SHARED src/graph.cpp src/network.cpp)
+add_library(frac SHARED src/graph.cpp src/network.cpp src/problem.cpp)
 target_include_directories(frac PUBLIC include PRIVATE voronoi/src)
 
diff --git a/lib/include/current_info.hpp b/lib/include/current_info.hpp
index 90dd3c7..a84fb5b 100644
--- a/lib/include/current_info.hpp
+++ b/lib/include/current_info.hpp
@@ -4,7 +4,7 @@
 #include <vector>
 
 typedef struct current_info {
-  double conductivity;
+  std::array<double, 2> conductivity;
   std::vector<double> currents;
 } current_info;
 
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 15cb8ea..1ffe742 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -6,44 +6,14 @@
 #include <utility>
 #include <random>
 #include <limits>
+#include <valarray>
 
-#include <cholmod.h>
-
+#include "problem.hpp"
 #include "graph.hpp"
 #include "hooks.hpp"
 #include "current_info.hpp"
 #include "array_hash.hpp"
 
-#ifdef FRACTURE_LONGINT
-
-#define CHOL_F(x) cholmod_l_##x
-#define CHOL_INT long int
-
-#else
-
-#define CHOL_F(x) cholmod_##x
-#define CHOL_INT int
-
-#endif
-
-class problem {
-  private:
-    const graph& G;
-    cholmod_dense* b;
-    cholmod_factor* factor;
-    cholmod_sparse* voltcurmat;
-    cholmod_common* c;
-    unsigned axis;
-
-  public:
-    problem(const graph&, unsigned axis, cholmod_common*);
-    problem(const graph&, unsigned axis, cholmod_sparse*, cholmod_common*);
-    problem(const problem&);
-    ~problem();
-    current_info solve(std::vector<bool>& fuses);
-    void break_edge(unsigned, bool unbreak = false);
-};
-
 class network {
   public:
     const graph& G;
@@ -54,45 +24,49 @@ class network {
     network(const network&);
 
     void set_thresholds(double, std::mt19937&);
+    void fracture(hooks&, bool one_axis = false);
+
     virtual void break_edge(unsigned e, bool unbreak = false) {fuses[e] = !unbreak;};
     virtual current_info get_current_info() {current_info empty; return empty;};
 };
 
 class fuse_network : public network {
-  public:
-    problem ohm;
+  private:
+    double weight;
+    problem ohm_x;
+    problem ohm_y;
 
-    fuse_network(const graph&, cholmod_common*);
+  public:
+    fuse_network(const graph&, cholmod_common*, double weight = 1.0);
+    fuse_network(const fuse_network&);
 
-    void fracture(hooks&, double cutoff = 1e-13);
+    void break_edge(unsigned, bool unbreak = false) override;
+    current_info get_current_info() override;
 };
 
 class elastic_network : public network {
   private:
     double weight;
-  public:
     problem hook_x;
     problem hook_y;
 
-    elastic_network(const graph&, cholmod_common*);
+  public:
+    elastic_network(const graph&, cholmod_common*, double weight = 0.5);
     elastic_network(const elastic_network&);
 
-    void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10);
     void break_edge(unsigned, bool unbreak = false) override;
     current_info get_current_info() override;
 };
 
 class percolation_network : public network {
   private:
-    double weight;
-  public:
-    problem hook_x;
-    problem hook_y;
+    problem px;
+    problem py;
 
+  public:
     percolation_network(const graph&, cholmod_common*);
     percolation_network(const percolation_network&);
 
-    void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10);
     void break_edge(unsigned, bool unbreak = false) override;
     current_info get_current_info() override;
 };
diff --git a/lib/include/problem.hpp b/lib/include/problem.hpp
new file mode 100644
index 0000000..c2c9e09
--- /dev/null
+++ b/lib/include/problem.hpp
@@ -0,0 +1,38 @@
+
+#include <vector>
+#include <limits>
+
+#include <cholmod.h>
+#include "graph.hpp"
+#include "current_info.hpp"
+
+#ifdef FRACTURE_LONGINT
+
+#define CHOL_F(x) cholmod_l_##x
+#define CHOL_INT long int
+
+#else
+
+#define CHOL_F(x) cholmod_##x
+#define CHOL_INT int
+
+#endif
+
+class problem {
+  private:
+    const graph& G;
+    unsigned axis;
+    cholmod_dense* b;
+    cholmod_factor* factor;
+    cholmod_sparse* voltcurmat;
+    cholmod_common* c;
+
+  public:
+    problem(const graph&, unsigned, cholmod_common*);
+    problem(const graph&, unsigned, cholmod_sparse*, cholmod_common*);
+    problem(const problem&);
+    ~problem();
+    current_info solve(std::vector<bool>& fuses);
+    void break_edge(unsigned, bool unbreak = false);
+};
+
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 1d9b3f9..b081a3c 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -1,13 +1,6 @@
 
 #include "network.hpp"
-
-class nanException: public std::exception
-{
-  virtual const char* what() const throw()
-  {
-    return "The linear problem returned NaN.";
-  }
-} nanex;
+#include <iostream>
 
 class nofuseException: public std::exception
 {
@@ -17,210 +10,6 @@ class nofuseException: public std::exception
   }
 } nofuseex;
 
-problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), voltcurmat(vcmat), c(c), axis(axis) {
-  b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
-  for (unsigned i = 0; i < G.edges.size(); i++) {
-    graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
-    graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
-
-    if (G.edges[i].crossings[axis]) {
-      bool ind;
-      if (axis == 1) {
-        ind = v0.y < v1.y ? 0 : 1;
-      } else {
-        ind = v0.x < v1.x ? 0 : 1;
-      }
-
-      ((double *)b->x)[G.edges[i].v[ind]] += 1.0;
-      ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0;
-    }
-  }
-
-  unsigned nnz = G.vertices.size() + G.edges.size();
-
-  cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz,  1, CHOLMOD_REAL, c);
-
-  for (unsigned i = 0; i < G.vertices.size(); i++) {
-    ((CHOL_INT *)t->i)[i] = i;
-    ((CHOL_INT *)t->j)[i] = i;
-    ((double *)t->x)[i] = 0.0;
-  }
-
-  unsigned terms = G.vertices.size();
-
-  std::unordered_map<std::array<unsigned, 2>, unsigned> known_edges;
-
-  for (unsigned i = 0; i < G.edges.size(); i++) {
-    unsigned v0 = G.edges[i].v[0];
-    unsigned v1 = G.edges[i].v[1];
-
-    ((double *)t->x)[v0]++;
-    ((double *)t->x)[v1]++;
-
-    unsigned s0 = v0 < v1 ? v0 : v1;
-    unsigned s1 = v0 < v1 ? v1 : v0;
-
-    auto it = known_edges.find({s0, s1});
-
-    if (it == known_edges.end()) {
-      ((CHOL_INT *)t->i)[terms] = s1;
-      ((CHOL_INT *)t->j)[terms] = s0;
-      ((double *)t->x)[terms] = -1.0;
-
-      known_edges[{s0, s1}] = terms;
-      terms++;
-    } else {
-      ((double *)t->x)[it->second] -= 1.0;
-    }
-  }
-
-  ((double *)t->x)[0]++;
-
-  t->nnz = terms;
-
-  cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c);
-  CHOL_F(free_triplet)(&t, c);
-  factor = CHOL_F(analyze)(laplacian, c);
-  CHOL_F(factorize)(laplacian, factor, c);
-  CHOL_F(free_sparse)(&laplacian, c);
-}
-
-problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) {
-  cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
-
-  t->nnz = 2 * G.edges.size();
-
-  for (unsigned i = 0; i < G.edges.size(); i++) {
-    ((CHOL_INT *)t->i)[2 * i] = i;
-    ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0];
-    ((double *)t->x)[2 * i] = 1.0;
-
-    ((CHOL_INT *)t->i)[2 * i + 1] = i;
-    ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1];
-    ((double *)t->x)[2 * i + 1] = -1.0;
-  }
-
-  voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c);
-
-  CHOL_F(free_triplet)(&t, c);
-}
-
-problem::problem(const problem& other) : G(other.G), c(other.c), axis(other.axis) {
-  b = CHOL_F(copy_dense)(other.b, c);
-  factor = CHOL_F(copy_factor)(other.factor, c);
-  voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
-} 
-
-problem::~problem() {
-  CHOL_F(free_dense)(&b, c);
-  CHOL_F(free_factor)(&factor, c);
-  CHOL_F(free_sparse)(&voltcurmat, c);
-}
-
-current_info problem::solve(std::vector<bool>& fuses) {
-  cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
-
-  if (((double *)x->x)[0] != ((double *)x->x)[0]) {
-    throw nanex;
-  }
-
-  cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
-
-  double alpha[2] = {1, 0};
-  double beta[2] = {0, 0};
-  CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
-
-  std::vector<double> currents(G.edges.size());
-
-  double total_current = 0;
-
-  for (int i = 0; i < G.edges.size(); i++) {
-    if (fuses[i]) {
-      currents[i] = 0;
-    } else {
-      currents[i] = ((double *)y->x)[i];
-
-      graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
-      graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
-
-      if (G.edges[i].crossings[axis]) {
-        bool comp;
-        if (axis == 1) {
-          comp = v0.y > v1.y;
-        } else {
-          comp = v0.x > v1.x;
-        }
-
-        if (comp) {
-          currents[i] += 1.0;
-          total_current += currents[i];
-        } else {
-          currents[i] -= 1.0;
-          total_current -= currents[i];
-        }
-      }
-    }
-  }
-
-  CHOL_F(free_dense)(&x, c);
-  CHOL_F(free_dense)(&y, c);
-
-  return {total_current, currents};
-}
-
-void problem::break_edge(unsigned e, bool unbreak) {
-  unsigned v0 = G.edges[e].v[0];
-  unsigned v1 = G.edges[e].v[1];
-
-  unsigned n = factor->n;
-
-  cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
-
-  unsigned s1, s2;
-  s1 = v0 < v1 ? v0 : v1;
-  s2 = v0 < v1 ? v1 : v0;
-
-  CHOL_INT *pp = (CHOL_INT *)update_mat->p;
-  CHOL_INT *ii = (CHOL_INT *)update_mat->i;
-  double *xx = (double *)update_mat->x;
-
-  for (unsigned i = 0; i <= s1; i++) {
-    pp[i] = 0;
-  }
-
-  for (unsigned i = s1 + 1; i <= n; i++) {
-    pp[i] = 2;
-  }
-
-  ii[0] = s1;
-  ii[1] = s2;
-  xx[0] = 1;
-  xx[1] = -1;
-
-  cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
-      update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c);
-
-  CHOL_F(updown)(unbreak, perm_update_mat, factor, c);
-
-  CHOL_F(free_sparse)(&perm_update_mat, c);
-  CHOL_F(free_sparse)(&update_mat, c);
-
-  graph::coordinate r0 = G.vertices[v0].r;
-  graph::coordinate r1 = G.vertices[v1].r;
-
-  if (G.edges[e].crossings[axis]) {
-    bool ind;
-    if (axis == 1) {
-      ind = r0.y < r1.y ? unbreak : !unbreak;
-    } else {
-      ind = r0.x < r1.x ? unbreak : !unbreak;
-    }
-
-    ((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
-    ((double *)b->x)[G.edges[e].v[!ind]] += 1.0;
-  }
-}
-
 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) {}
@@ -244,16 +33,19 @@ void network::set_thresholds(double beta, std::mt19937& rng) {
   }
 }
 
-fuse_network::fuse_network(const graph& G, cholmod_common *c) : network(G), ohm(G, 1, c) {
-}
-
-void fuse_network::fracture(hooks& m, double cutoff) {
+void network::fracture(hooks& m, bool one_axis) {
   m.pre_fracture(*this);
 
   while (true) {
-    current_info ci = ohm.solve(fuses);
+    current_info c = this->get_current_info();
+
+    double min_cond = 1.0 / G.edges.size();
 
-    if (ci.conductivity < cutoff * G.vertices.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;
     }
 
@@ -261,8 +53,9 @@ void fuse_network::fracture(hooks& m, double cutoff) {
     long double max_val = std::numeric_limits<long double>::lowest();
 
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
-        long double val = logl(fabs(ci.currents[i])) - thresholds[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;
@@ -274,172 +67,130 @@ void fuse_network::fracture(hooks& m, double cutoff) {
       throw nofuseex;
     }
 
-    fuses[max_pos] = true;
-    ohm.break_edge(max_pos);
-
-    m.bond_broken(*this, ci, max_pos);
+    this->break_edge(max_pos);
+    m.bond_broken(*this, c, max_pos);
   }
 
   m.post_fracture(*this);
 }
 
-elastic_network::elastic_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) {
-}
-
-elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) {
-}
-
-void elastic_network::fracture(hooks& m, double weight, double cutoff) {
-  this->weight = weight;
-  m.pre_fracture(*this);
-
-  while (true) {
-    current_info ctot = this->get_current_info();
-
-    if (ctot.conductivity == 0) {
-      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] && ctot.currents[i] > 1.0 / G.edges.size()) {
-        long double val = logl(ctot.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);
+fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) :
+  network(G), ohm_x(G, 0, c), ohm_y(G, 1, c) {}
 
-    m.bond_broken(*this, ctot, max_pos);
-  }
+fuse_network::fuse_network(const fuse_network& n) :
+  network(n), ohm_x(n.ohm_x), ohm_y(n.ohm_y) {
+}
 
-  m.post_fracture(*this);
+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 elastic_network::get_current_info() {
-  current_info cx = hook_x.solve(fuses);
-  current_info cy = hook_y.solve(fuses);
+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 < 1.0 / G.edges.size();
-  bool done_y = cy.conductivity < 1.0 / G.edges.size();
+  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) {
-    ctot.conductivity = 0;
-  } else if (done_x) {
-    ctot.conductivity = cy.conductivity;
+  if (done_x && !done_y) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity;
+      ctot.currents[i] = fabs(weight * cy.currents[i] / cy.conductivity[1]);
     }
-  } else if (done_y) {
-    ctot.conductivity = cx.conductivity;
+  } 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;
+      ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0]);
     }
-  } else {
-    ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); 
+  } 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, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2));
+      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);
 }
 
-percolation_network::percolation_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) {
-}
+current_info elastic_network::get_current_info() {
+  current_info cx = hook_x.solve(fuses);
+  current_info cy = hook_y.solve(fuses);
 
-percolation_network::percolation_network(const percolation_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) {
-}
+  bool done_x = cx.conductivity[0] < 1.0 / G.edges.size();
+  bool done_y = cy.conductivity[1] < 1.0 / G.edges.size();
 
-void percolation_network::fracture(hooks& m, double weight, double cutoff) {
-  this->weight = weight;
-  m.pre_fracture(*this);
-
-  while (true) {
-    current_info ctot = this->get_current_info();
+  current_info ctot;
+  ctot.currents.resize(G.edges.size());
+  ctot.conductivity[0] = cx.conductivity[0]; 
+  ctot.conductivity[1] = cy.conductivity[1]; 
 
-    if (ctot.conductivity == 0) {
-      break;
+  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];
     }
-
-    unsigned max_pos = UINT_MAX;
-    long double max_val = std::numeric_limits<long double>::lowest();
-
+  } else if (done_y && !done_x) {
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      if (!fuses[i] && ctot.currents[i] > 1.0 / G.edges.size()) {
-        long double val = - thresholds[i];
-        if (val > max_val) {
-          max_val = val;
-          max_pos = i;
-        }
-      }
+      ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity[0];
     }
-
-    if (max_pos == UINT_MAX)  {
-      throw nofuseex;
+  } 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));
     }
-
-    this->break_edge(max_pos);
-
-    m.bond_broken(*this, ctot, max_pos);
   }
 
-  m.post_fracture(*this);
+  return ctot;
 }
 
-current_info percolation_network::get_current_info() {
-  current_info cx = hook_x.solve(fuses);
-  current_info cy = hook_y.solve(fuses);
 
-  bool done_x = cx.conductivity < 1.0 / G.edges.size();
-  bool done_y = cy.conductivity < 1.0 / G.edges.size();
+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());
+  ctot.currents.resize(G.edges.size(), 0);
 
-  if (done_x && done_y) {
-    ctot.conductivity = 0;
-  } else if (done_x) {
-    ctot.conductivity = cy.conductivity;
-    for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity;
-    }
-  } else if (done_y) {
-    ctot.conductivity = cx.conductivity;
-    for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity;
-    }
-  } else {
-    ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); 
-    for (unsigned i = 0; i < G.edges.size(); i++) {
-      ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2));
+  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;
-  hook_x.break_edge(e, unbreak);
-  hook_y.break_edge(e, unbreak);
+  px.break_edge(e, unbreak);
+  py.break_edge(e, unbreak);
 }
 
diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp
new file mode 100644
index 0000000..f66a35c
--- /dev/null
+++ b/lib/src/problem.cpp
@@ -0,0 +1,217 @@
+
+#include "problem.hpp"
+
+class nanException: public std::exception
+{
+  virtual const char* what() const throw()
+  {
+    return "The linear problem returned NaN.";
+  }
+} nanex;
+
+problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), axis(axis), voltcurmat(vcmat), c(c) {
+  b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
+  for (unsigned i = 0; i < G.edges.size(); i++) {
+    graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+    graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+
+    if (G.edges[i].crossings[axis]) {
+      bool ind;
+      if (axis == 1) {
+        ind = v0.y < v1.y ? 0 : 1;
+      } else {
+        ind = v0.x < v1.x ? 0 : 1;
+      }
+
+      ((double *)b->x)[G.edges[i].v[ind]]  += 1.0;
+      ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0;
+    }
+  }
+
+  unsigned nnz = G.vertices.size() + G.edges.size();
+
+  cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz,  1, CHOLMOD_REAL, c);
+
+  for (unsigned i = 0; i < G.vertices.size(); i++) {
+    ((CHOL_INT *)t->i)[i] = i;
+    ((CHOL_INT *)t->j)[i] = i;
+    ((double *)t->x)[i] = 0.0;
+  }
+
+  unsigned terms = G.vertices.size();
+
+  std::unordered_map<std::array<unsigned, 2>, unsigned> known_edges;
+
+  for (unsigned i = 0; i < G.edges.size(); i++) {
+    unsigned v0 = G.edges[i].v[0];
+    unsigned v1 = G.edges[i].v[1];
+
+    ((double *)t->x)[v0]++;
+    ((double *)t->x)[v1]++;
+
+    unsigned s0 = v0 < v1 ? v0 : v1;
+    unsigned s1 = v0 < v1 ? v1 : v0;
+
+    auto it = known_edges.find({s0, s1});
+
+    if (it == known_edges.end()) {
+      ((CHOL_INT *)t->i)[terms] = s1;
+      ((CHOL_INT *)t->j)[terms] = s0;
+      ((double *)t->x)[terms] = -1.0;
+
+      known_edges[{s0, s1}] = terms;
+      terms++;
+    } else {
+      ((double *)t->x)[it->second] -= 1.0;
+    }
+  }
+
+  ((double *)t->x)[0]++;
+
+  t->nnz = terms;
+
+  cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c);
+  CHOL_F(free_triplet)(&t, c);
+  factor = CHOL_F(analyze)(laplacian, c);
+  CHOL_F(factorize)(laplacian, factor, c);
+  CHOL_F(free_sparse)(&laplacian, c);
+}
+
+problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) {
+  cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
+
+  t->nnz = 2 * G.edges.size();
+
+  for (unsigned i = 0; i < G.edges.size(); i++) {
+    ((CHOL_INT *)t->i)[2 * i] = i;
+    ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0];
+    ((double *)t->x)[2 * i] = 1.0;
+
+    ((CHOL_INT *)t->i)[2 * i + 1] = i;
+    ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1];
+    ((double *)t->x)[2 * i + 1] = -1.0;
+  }
+
+  voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c);
+
+  CHOL_F(free_triplet)(&t, c);
+}
+
+problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c) {
+  b = CHOL_F(copy_dense)(other.b, c);
+  factor = CHOL_F(copy_factor)(other.factor, c);
+  voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
+} 
+
+problem::~problem() {
+  CHOL_F(free_dense)(&b, c);
+  CHOL_F(free_factor)(&factor, c);
+  CHOL_F(free_sparse)(&voltcurmat, c);
+}
+
+current_info problem::solve(std::vector<bool>& fuses) {
+  cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+
+  if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+    throw nanex;
+  }
+
+  cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+
+  double alpha[2] = {1, 0};
+  double beta[2] = {0, 0};
+  CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+  std::vector<double> currents(G.edges.size());
+
+  std::array<double, 2> total_current = {0, 0};
+
+  for (int i = 0; i < G.edges.size(); i++) {
+    if (fuses[i]) {
+      currents[i] = 0;
+    } else {
+      currents[i] = ((double *)y->x)[i];
+
+      graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+      graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+
+      if (G.edges[i].crossings[axis]) {
+        bool comp;
+        if (axis == 1) {
+          comp = v0.y > v1.y;
+        } else {
+          comp = v0.x > v1.x;
+        }
+
+        if (comp) {
+          currents[i] += 1.0;
+          total_current[axis] += currents[i];
+        } else {
+          currents[i] -= 1.0;
+          total_current[axis] -= currents[i];
+        }
+      }
+    }
+  }
+
+  total_current[!axis] = 0.0;
+
+  CHOL_F(free_dense)(&x, c);
+  CHOL_F(free_dense)(&y, c);
+
+  return {total_current, currents};
+}
+
+void problem::break_edge(unsigned e, bool unbreak) {
+  unsigned v0 = G.edges[e].v[0];
+  unsigned v1 = G.edges[e].v[1];
+
+  unsigned n = factor->n;
+
+  cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
+
+  unsigned s1, s2;
+  s1 = v0 < v1 ? v0 : v1;
+  s2 = v0 < v1 ? v1 : v0;
+
+  CHOL_INT *pp = (CHOL_INT *)update_mat->p;
+  CHOL_INT *ii = (CHOL_INT *)update_mat->i;
+  double *xx = (double *)update_mat->x;
+
+  for (unsigned i = 0; i <= s1; i++) {
+    pp[i] = 0;
+  }
+
+  for (unsigned i = s1 + 1; i <= n; i++) {
+    pp[i] = 2;
+  }
+
+  ii[0] = s1;
+  ii[1] = s2;
+  xx[0] = 1;
+  xx[1] = -1;
+
+  cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
+      update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c);
+
+  CHOL_F(updown)(unbreak, perm_update_mat, factor, c);
+
+  CHOL_F(free_sparse)(&perm_update_mat, c);
+  CHOL_F(free_sparse)(&update_mat, c);
+
+  graph::coordinate r0 = G.vertices[v0].r;
+  graph::coordinate r1 = G.vertices[v1].r;
+
+  if (G.edges[e].crossings[axis]) {
+    bool ind;
+    if (axis == 1) {
+      ind = r0.y < r1.y ? unbreak : !unbreak;
+    } else {
+      ind = r0.x < r1.x ? unbreak : !unbreak;
+    }
+
+    ((double *)b->x)[G.edges[e].v[ind]]  -= 1.0;
+    ((double *)b->x)[G.edges[e].v[!ind]] += 1.0;
+  }
+}
+