summaryrefslogtreecommitdiff
path: root/lib/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src')
-rw-r--r--lib/src/network.cpp631
1 files changed, 340 insertions, 291 deletions
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index ed23d45..6b3873b 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -1,28 +1,13 @@
#include "network.hpp"
-class nofuseException: public std::exception
-{
- virtual const char* what() const throw()
- {
- return "No valid fuse was available to break.";
- }
+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);
-}
+network::network(const graph& G, cholmod_common* c)
+ : G(G), C(G.dual_vertices.size()), fuses(G.edges.size()), backbone(G.edges.size()),
+ thresholds(G.edges.size()), px(G, 0, c), py(G, 1, c) {}
void network::set_thresholds(double beta, std::mt19937& rng) {
if (beta == 0.0) {
@@ -78,7 +63,7 @@ void network::fracture(hooks& m, bool one_axis) {
}
}
- if (max_pos == UINT_MAX) {
+ if (max_pos == UINT_MAX) {
throw nofuseex;
}
@@ -89,155 +74,199 @@ void network::fracture(hooks& m, bool one_axis) {
m.post_fracture(*this);
}
-void network::update_backbone(const std::vector<double>& c) {
+void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
+ std::set<unsigned>& seen_vertices, unsigned v_prev,
+ unsigned v_cur) const {
+ for (unsigned ei : G.dual_vertices[v_cur].ne) {
+ const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
+ unsigned vn = e[0] == v_cur ? e[1] : e[0];
+
+ if (vn != v_prev) {
+ auto it = seen_vertices.find(vn);
+ if (it == seen_vertices.end()) {
+ // if (seen_vertices.contains()) { // need to wait for better c++20 support...
+ cycle_edges.insert(vn);
+ } else {
+ this->get_cycle_edges_helper(cycle_edges, seen_vertices, v_cur, vn);
+ }
+ }
+ }
+}
- 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;
+std::set<unsigned> network::get_cycle_edges(unsigned v0) const {
+ std::set<unsigned> seen_vertices;
+ std::set<unsigned> cycle_edges;
- get_cycle_edges(std::list<unsigned>& bE) : bE(bE) {};
+ this->get_cycle_edges_helper(cycle_edges, seen_vertices, v0, v0);
- void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
- tE.insert(g[e].index);
- }
+ return cycle_edges;
+}
- 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);
- }
- }
- };
+std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, unsigned v1) const {
+ for (unsigned ei : G.dual_vertices[v0].ne) {
+ const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
+ unsigned vn = e[0] == v0 ? e[1] : e[0];
+ if (vn == v1) {
+ return {ei};
+ } else {
+ std::list<unsigned> edges = this->get_cycle_edges(vn, v1);
+ edges.splice(edges.end(), {ei});
+ return edges;
+ }
+ }
+
+ return {};
+}
+
+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 find_cycle : public boost::default_dfs_visitor {
- public:
- const graph& G;
- std::array<unsigned, 2>& E;
- unsigned end;
- struct done{};
+ 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) {}
+ 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 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 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]--;
- }
+ 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;
+ public:
+ std::set<unsigned>& E;
- get_cluster_edges(std::set<unsigned>& E) : E(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);
- }
+ 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)) {
+ if ((!backbone[i]) && C.same_component(G.dual_edges[i].v[0], G.dual_edges[i].v[1])) {
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&) {
+ std::list<unsigned> cedges = this->get_cycle_edges(G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
+
+ 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);
}
- 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 (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 (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);
+ 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);
+ 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 {
@@ -279,8 +308,6 @@ void network::update_backbone(const std::vector<double>& c) {
RE.insert(e);
}
}
-
-
};
std::set<unsigned> bb_verts;
@@ -297,161 +324,189 @@ void network::update_backbone(const std::vector<double>& c) {
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 l = (i + j) % G.vertices[v].nd.size();
+ if (C.same_component(G.vertices[v].nd[i], G.vertices[v].nd[l])) {
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++;
+ 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 (!C.same_component(G.vertices[v].nd[i], G.vertices[v].nd[k])) {
+ any_intervening1 = 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()]]) {
- ie2++;
+ for (unsigned k = (ig + 1); k % G.vertices[v].nd.size() != il; k++) {
+ if (!C.same_component(G.vertices[v].nd[i],
+ G.vertices[v].nd[k % G.vertices[v].nd.size()])) {
+ any_intervening2 = true;
+ break;
}
}
- }
- 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 (any_intervening2 && !any_intervening1) {
+ for (unsigned k = il + 1; k <= ig; k++) {
+ if (!backbone[G.vertices[v].ne[k]]) {
+ ie1++;
}
}
- 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];
+ }
+ 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<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];
- }
+ 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]++;
}
- 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];
+ 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 (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 (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);
+ }
+ }
+ }
}
- }
- 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);
}
}
-
- }
- }
- for (unsigned e : cedges) {
- boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
- }
- }
- }
+ }
}
}
}
@@ -461,15 +516,13 @@ void network::update_backbone(const std::vector<double>& c) {
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]);
+ C.add_bond(G.dual_edges[e]);
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 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) {}
@@ -482,7 +535,7 @@ current_info fuse_network::get_current_info() {
current_info ctot;
ctot.currents.resize(G.edges.size());
- ctot.conductivity = {px.sol.conductivity[0], py.sol.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++) {
@@ -495,19 +548,17 @@ current_info fuse_network::get_current_info() {
} 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]);
+ 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 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) {}
+elastic_network::elastic_network(const elastic_network& n) : network(n), weight(n.weight) {}
current_info elastic_network::get_current_info() {
px.solve(fuses);
@@ -518,8 +569,8 @@ current_info elastic_network::get_current_info() {
current_info ctot;
ctot.currents.resize(G.edges.size());
- ctot.conductivity[0] = px.sol.conductivity[0];
- ctot.conductivity[1] = py.sol.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++) {
@@ -539,7 +590,6 @@ current_info elastic_network::get_current_info() {
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) {}
@@ -558,7 +608,6 @@ current_info percolation_network::get_current_info() {
ctot.currents[i] = 1.0;
}
}
-
+
return ctot;
}
-