summaryrefslogtreecommitdiff
path: root/lib/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src')
-rw-r--r--lib/src/graph.cpp54
-rw-r--r--lib/src/network.cpp479
-rw-r--r--lib/src/problem.cpp26
3 files changed, 490 insertions, 69 deletions
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 974d3b3..bfdb952 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -33,6 +33,8 @@ graph::graph(unsigned Nx, unsigned Ny) {
for (unsigned i = 0; i < nv; i++) {
vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
vertices[i].r.y = (double)(i / (Nx / 2));
+ signed f = (1 + i / (Nx / 2)) % 2 == 1 ? 1 : -1;
+ vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
vertices[i].polygon = {
{vertices[i].r.x - 1.0, vertices[i].r.y},
{vertices[i].r.x, vertices[i].r.y - 1.0},
@@ -42,6 +44,7 @@ graph::graph(unsigned Nx, unsigned Ny) {
dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
dual_vertices[i].r.y = (double)(i / (Nx / 2));
+ dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx * (i / (Nx / 2)) + (i % (Nx / 2)), (nv + i - Nx / 2) % nv};
dual_vertices[i].polygon = {
{dual_vertices[i].r.x - 1.0, vertices[i].r.y},
{dual_vertices[i].r.x, vertices[i].r.y - 1.0},
@@ -67,6 +70,29 @@ graph::graph(unsigned Nx, unsigned Ny) {
}
}
+ for (vertex& v : vertices) {
+ v.ne.resize(v.nd.size());
+ }
+
+ for (unsigned i = 0; i < edges.size(); i++) {
+ for (unsigned vi : edges[i].v) {
+ auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]);
+ auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]);
+
+ unsigned d1 = std::distance(vertices[vi].nd.begin(), it1);
+ unsigned d2 = std::distance(vertices[vi].nd.begin(), it2);
+
+ unsigned dv1 = d1 < d2 ? d1 : d2;
+ unsigned dv2 = d1 < d2 ? d2 : d1;
+
+ if (dv2 - dv1 == 1) {
+ vertices[vi].ne[dv1] = i;
+ } else {
+ vertices[vi].ne[dv2] = i;
+ }
+ }
+ }
+
}
class eulerException: public std::exception
@@ -239,12 +265,15 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
if (it1 == known_vertices.end()) {
vi1 = dual_vertices.size();
- dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{{site->p.x, site->p.y}}});
+ dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{i1},{},{{site->p.x, site->p.y}}});
known_vertices[t1] = vi1;
} else {
vi1 = it1->second;
+ dual_vertices[vi1].nd.push_back(i1);
dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y});
}
+
+ vertices[i1].nd.push_back(vi1);
if (i1 < i2 || (i1 == i2 && !self_bonded)) {
if (i1 == i2) {
@@ -359,6 +388,29 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
}
}
+ for (vertex& v : vertices) {
+ v.ne.resize(v.nd.size());
+ }
+
+ for (unsigned i = 0; i < edges.size(); i++) {
+ for (unsigned vi : edges[i].v) {
+ auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]);
+ auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]);
+
+ unsigned d1 = std::distance(vertices[vi].nd.begin(), it1);
+ unsigned d2 = std::distance(vertices[vi].nd.begin(), it2);
+
+ unsigned dv1 = d1 < d2 ? d1 : d2;
+ unsigned dv2 = d1 < d2 ? d2 : d1;
+
+ if (dv2 - dv1 == 1) {
+ vertices[vi].ne[dv1] = i;
+ } else {
+ vertices[vi].ne[dv2] = i;
+ }
+ }
+ }
+
jcv_diagram_free(&diagram);
}
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 3424b94..dcc5e0e 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -9,9 +9,20 @@ class nofuseException: public std::exception
}
} nofuseex;
-network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {}
+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), fuses(n.fuses), thresholds(n.thresholds) {}
+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) {
@@ -36,10 +47,15 @@ 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();
+ 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;
}
@@ -52,7 +68,7 @@ void network::fracture(hooks& m, bool one_axis) {
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] > CURRENT_CUTOFF) {
+ if (!backbone[i]) {
long double val = logl(c.currents[i]) - thresholds[i];
if (val > max_val) {
@@ -73,43 +89,413 @@ void network::fracture(hooks& m, bool one_axis) {
m.post_fracture(*this);
}
+void network::update_backbone(const std::vector<double>& c) {
+ 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(), [](std::array<unsigned, 2> c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) {
+ 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;
+ }
-fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) :
- network(G), ohm_x(G, 0, c), ohm_y(G, 1, c) {}
+ void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
+ if (!initialized && disc > 0) {
+ start_edge = vd;
+ initialized = true;
+ }
+ disc++;
+ }
-fuse_network::fuse_network(const fuse_network& n) :
- network(n), ohm_x(n.ohm_x), ohm_y(n.ohm_y) {
+ 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;
+
+ 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]]) {
+ any_intervening_e_1 = true;
+ break;
+ }
+ }
+ }
+ 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()]]) {
+ any_intervening_e_2 = true;
+ break;
+ }
+ }
+ }
+ if ((any_intervening1 && any_intervening2) || (any_intervening1 && any_intervening_e_2) || (any_intervening2 && any_intervening_e_1)) {
+ found_pair = true;
+ p1 = il;
+ p2 = ig;
+ d1 = G.vertices[v].nd[il];
+ d2 = G.vertices[v].nd[ig];
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ if (found_pair) {
+ 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(), [](std::array<unsigned, 2> c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) {
+ 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 fuse_network::break_edge(unsigned e, bool unbreak) {
+void network::break_edge(unsigned e, bool unbreak) {
fuses[e] = !unbreak;
- ohm_x.break_edge(e, unbreak);
- ohm_y.break_edge(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() {
- current_info cx = ohm_x.solve(fuses);
- current_info cy = ohm_y.solve(fuses);
+ px.solve(fuses);
+ py.solve(fuses);
- bool done_x = cx.conductivity[0] < 1.0 / G.edges.size();
- bool done_y = cy.conductivity[1] < 1.0 / G.edges.size();
+ 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 = {cx.conductivity[0], cy.conductivity[1]};
+ 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 * cy.currents[i] / cy.conductivity[1]);
+ 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) * cx.currents[i] / cx.conductivity[0]);
+ 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) * cx.currents[i] / cx.conductivity[0] +
- weight * cy.currents[i] / cy.conductivity[1]);
+ ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0] +
+ weight * py.sol.currents[i] / py.sol.conductivity[1]);
}
}
@@ -118,42 +504,35 @@ current_info fuse_network::get_current_info() {
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) {}
+ network(G, c), weight(weight) {}
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);
-}
+ network(n), weight(n.weight) {}
current_info elastic_network::get_current_info() {
- current_info cx = hook_x.solve(fuses);
- current_info cy = hook_y.solve(fuses);
+ px.solve(fuses);
+ py.solve(fuses);
- bool done_x = cx.conductivity[0] < 1.0 / G.edges.size();
- bool done_y = cy.conductivity[1] < 1.0 / G.edges.size();
+ 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] = cx.conductivity[0];
- ctot.conductivity[1] = cy.conductivity[1];
+ 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(cy.currents[i]) / cy.conductivity[1];
+ 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(cx.currents[i]) / cx.conductivity[0];
+ 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) * cx.currents[i] / cx.conductivity[0], 2) +
- pow(weight * cy.currents[i] / cy.conductivity[1], 2));
+ 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));
}
}
@@ -161,23 +540,21 @@ current_info elastic_network::get_current_info() {
}
-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 graph& G, cholmod_common* c) : network(G, c) {}
-percolation_network::percolation_network(const percolation_network& n) :
- network(n), px(n.px), py(n.py) {}
+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);
- current_info cx = px.solve(fuses);
- current_info cy = py.solve(fuses);
+ px.solve(fuses);
+ py.solve(fuses);
- ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]};
+ ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]};
for (unsigned i = 0; i < G.edges.size(); i++) {
- if (fabs(cx.currents[i]) > CURRENT_CUTOFF || fabs(cy.currents[i]) > CURRENT_CUTOFF) {
+ if (fabs(px.sol.currents[i]) > CURRENT_CUTOFF || fabs(py.sol.currents[i]) > CURRENT_CUTOFF) {
ctot.currents[i] = 1.0;
}
}
@@ -185,9 +562,3 @@ current_info percolation_network::get_current_info() {
return ctot;
}
-void percolation_network::break_edge(unsigned e, bool unbreak) {
- fuses[e] = !unbreak;
- px.break_edge(e, unbreak);
- py.break_edge(e, unbreak);
-}
-
diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp
index f66a35c..0645244 100644
--- a/lib/src/problem.cpp
+++ b/lib/src/problem.cpp
@@ -75,6 +75,8 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
factor = CHOL_F(analyze)(laplacian, c);
CHOL_F(factorize)(laplacian, factor, c);
CHOL_F(free_sparse)(&laplacian, c);
+
+ sol.currents.resize(G.edges.size());
}
problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) {
@@ -97,7 +99,7 @@ problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G,
CHOL_F(free_triplet)(&t, c);
}
-problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c) {
+problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c), sol(other.sol) {
b = CHOL_F(copy_dense)(other.b, c);
factor = CHOL_F(copy_factor)(other.factor, c);
voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
@@ -109,7 +111,7 @@ problem::~problem() {
CHOL_F(free_sparse)(&voltcurmat, c);
}
-current_info problem::solve(std::vector<bool>& fuses) {
+void 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]) {
@@ -122,15 +124,13 @@ current_info problem::solve(std::vector<bool>& fuses) {
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};
+ sol.conductivity = {0, 0};
for (int i = 0; i < G.edges.size(); i++) {
if (fuses[i]) {
- currents[i] = 0;
+ sol.currents[i] = 0;
} else {
- currents[i] = ((double *)y->x)[i];
+ sol.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;
@@ -144,22 +144,20 @@ current_info problem::solve(std::vector<bool>& fuses) {
}
if (comp) {
- currents[i] += 1.0;
- total_current[axis] += currents[i];
+ sol.currents[i] += 1.0;
+ sol.conductivity[axis] += sol.currents[i];
} else {
- currents[i] -= 1.0;
- total_current[axis] -= currents[i];
+ sol.currents[i] -= 1.0;
+ sol.conductivity[axis] -= sol.currents[i];
}
}
}
}
- total_current[!axis] = 0.0;
+ sol.conductivity[!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) {