summaryrefslogtreecommitdiff
path: root/lib/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src')
-rw-r--r--lib/src/network.cpp417
-rw-r--r--lib/src/problem.cpp217
2 files changed, 301 insertions, 333 deletions
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;
+ }
+}
+