From d00d4896955c38818685802de42b4e79e6ef933c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 23 Sep 2019 15:43:17 -0400 Subject: started migration away from boost algorithms is the backbone trimming code --- CMakeLists.txt | 4 +- lib/include/clusters.hpp | 52 ++++ lib/include/network.hpp | 49 ++-- lib/src/network.cpp | 631 +++++++++++++++++++++++++---------------------- 4 files changed, 411 insertions(+), 325 deletions(-) create mode 100644 lib/include/clusters.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f5b6d1a..f7477fa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,8 +3,8 @@ cmake_minimum_required(VERSION 3.8) project(fracture) -set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall") -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -march=native") +set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -std=c++2a") +set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -march=native -std=c++2a") set (CMAKE_CXX_STANDARD 17) set (CMAKE_C_STANDARD 11) diff --git a/lib/include/clusters.hpp b/lib/include/clusters.hpp new file mode 100644 index 0000000..4cc9073 --- /dev/null +++ b/lib/include/clusters.hpp @@ -0,0 +1,52 @@ + +#include +#include "graph.hpp" + +class ClusterTree { +private: + std::vector p; + std::vector o; + +public: + std::vector c; + + ClusterTree(unsigned n) : p(n, -1), o(n, 1), c(n, 0) { c[0] = n; } + + unsigned findroot(unsigned i) { + if (p[i] < 0) + return i; + else + return p[i] = this->findroot(p[i]); + } + + void join(unsigned i, unsigned j) { + c[o[i] - 1]--; + c[o[j] - 1]--; + + if (o[i] < o[j]) { + p[i] = j; + o[j] += o[i]; + c[o[j] - 1]++; + } else { + p[j] = i; + o[i] += o[j]; + c[o[i] - 1]++; + } + } + + void add_bond(const graph::edge& b) { + unsigned i = this->findroot(b.v[0]); + unsigned j = this->findroot(b.v[1]); + + if (i != j) { + this->join(i, j); + } + } + + bool same_component(unsigned v0, unsigned v1) { + unsigned i = this->findroot(v0); + unsigned j = this->findroot(v1); + + return i == j; + } +}; diff --git a/lib/include/network.hpp b/lib/include/network.hpp index 8211cfc..bf9015c 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -14,53 +14,38 @@ #include "hooks.hpp" #include "current_info.hpp" #include "array_hash.hpp" +#include "clusters.hpp" #define CURRENT_CUTOFF 1e-10 -#include -#include -#include -#include -#include -#include -#include -#include - -struct EdgeProperties { - unsigned index; -}; +class network { + private: + problem px; + problem py; + std::unordered_map, bool> seen_pairs; -typedef boost::adjacency_list Graph; -typedef boost::graph_traits::vertex_descriptor Vertex; -typedef boost::graph_traits::vertices_size_type VertexIndex; + void update_backbone(const std::vector& c); + void break_edge(unsigned, bool unbreak = false); + void get_cycle_edges_helper(std::set&, std::set&, unsigned, unsigned) const; + std::set get_cycle_edges(unsigned) const; + std::array find_cycle(const std::set&, unsigned, unsigned) const; -class network { public: const graph& G; - Graph bG; - Graph dG; - std::vector rank; - std::vector parent; - boost::disjoint_sets ds; - std::unordered_map, bool> seen_pairs; + ClusterTree C; std::vector fuses; std::vector backbone; std::vector thresholds; - problem px; - problem py; - network(const graph&, cholmod_common*); network(const network&); void set_thresholds(double, std::mt19937&); - void fracture(hooks&, bool one_axis = false); - void update_backbone(const std::vector& c); + void fracture(hooks&, bool one_axis = true); - void break_edge(unsigned, bool unbreak = false); - virtual current_info get_current_info() {current_info empty; return empty;}; + virtual current_info get_current_info() const {current_info empty; return empty;}; }; class fuse_network : public network { @@ -71,7 +56,7 @@ class fuse_network : public network { fuse_network(const graph&, cholmod_common*, double weight = 1.0); fuse_network(const fuse_network&); - current_info get_current_info() override; + current_info get_current_info() const override; }; class elastic_network : public network { @@ -82,7 +67,7 @@ class elastic_network : public network { elastic_network(const graph&, cholmod_common*, double weight = 0.5); elastic_network(const elastic_network&); - current_info get_current_info() override; + current_info get_current_info() const override; }; class percolation_network : public network { @@ -90,6 +75,6 @@ class percolation_network : public network { percolation_network(const graph&, cholmod_common*); percolation_network(const percolation_network&); - current_info get_current_info() override; + current_info get_current_info() const override; }; 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& c) { +void network::get_cycle_edges_helper(std::set& cycle_edges, + std::set& seen_vertices, unsigned v_prev, + unsigned v_cur) const { + for (unsigned ei : G.dual_vertices[v_cur].ne) { + const std::array& 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 tE; - std::list& bE; +std::set network::get_cycle_edges(unsigned v0) const { + std::set seen_vertices; + std::set cycle_edges; - get_cycle_edges(std::list& bE) : bE(bE) {}; + this->get_cycle_edges_helper(cycle_edges, seen_vertices, v0, v0); - void tree_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - tE.insert(g[e].index); - } + return cycle_edges; +} - void back_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - if (tE.find(g[e].index) == tE.end()) { - bE.push_back(g[e].index); - } - } - }; +std::array network::find_cycle(const std::set& cycle_edges, unsigned v0, unsigned v1) const { + for (unsigned ei : G.dual_vertices[v0].ne) { + const std::array& e = G.dual_edges[ei].v; + unsigned vn = e[0] == v0 ? e[1] : e[0]; + if (vn == v1) { + return {ei}; + } else { + std::list edges = this->get_cycle_edges(vn, v1); + edges.splice(edges.end(), {ei}); + return edges; + } + } + + return {}; +} + +void network::update_backbone(const std::vector& 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& E; - unsigned end; - struct done{}; + public: + const graph& G; + std::array& E; + unsigned end; + struct done {}; - find_cycle(const graph& G, std::array& E, unsigned end) : G(G), E(E), end(end) {} + find_cycle(const graph& G, std::array& E, unsigned end) : G(G), E(E), end(end) {} - void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { - if (v == end) { - throw done{}; - } + void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { + if (v == end) { + throw done{}; } + } - void examine_edge(boost::graph_traits::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::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::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::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& E; + public: + std::set& E; - get_cluster_edges(std::set& E) : E(E) {} + get_cluster_edges(std::set& E) : E(E) {} - void examine_edge(boost::graph_traits::edge_descriptor ed, const Graph& g) { - unsigned e = g[ed].index; - E.insert(e); - } + void examine_edge(boost::graph_traits::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 cedges = {}; - get_cycle_edges vis(cedges); - std::vector 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 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 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> other_cycles = {cycle}; - for (unsigned e : cedges) { - std::array 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 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 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 cedges = this->get_cycle_edges(G.dual_edges[i].v[0], G.dual_edges[i].v[1]); + + std::array 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 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> other_cycles = {cycle}; + for (unsigned e : cedges) { + std::array 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 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 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 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 component_edges = {}; - get_cluster_edges visc(component_edges); - std::vector 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 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 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 component_edges = {}; + get_cluster_edges visc(component_edges); + std::vector 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 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& c) { RE.insert(e); } } - - }; std::set bb_verts; @@ -297,161 +324,189 @@ void network::update_backbone(const std::vector& 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 cedges = {}; - get_cycle_edges vis(cedges); - std::vector 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 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 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> other_cycles = {cycle}; - for (unsigned e : cedges) { - std::array 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 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 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 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 left_edges = {}; - std::set right_edges = {}; - get_tie_edges visc(left_edges, right_edges); - std::vector 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 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 cedges = {}; + get_cycle_edges vis(cedges); + std::vector 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 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 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 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> other_cycles = {cycle}; + for (unsigned e : cedges) { + std::array 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 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 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 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 left_edges = {}; + std::set right_edges = {}; + get_tie_edges visc(left_edges, right_edges); + std::vector 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 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 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& 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; } - -- cgit v1.2.3-54-g00ecf