diff options
Diffstat (limited to 'lib/src')
| -rw-r--r-- | lib/src/network.cpp | 631 |
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; } - |