From 3f7f20f21f583ca2de566bea08a87eac4b17ad29 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 23 Sep 2019 23:19:13 -0400 Subject: successfully implemented the backbone trimming with homespun graph algorithms, measurements have mostly been disabled and need to be migrated --- lib/src/graph.cpp | 22 ++- lib/src/network.cpp | 472 +++++++++++++++++++++------------------------------- 2 files changed, 211 insertions(+), 283 deletions(-) (limited to 'lib/src') diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index bfdb952..a5063de 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -44,7 +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].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (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}, @@ -93,6 +93,16 @@ graph::graph(unsigned Nx, unsigned Ny) { } } + for (vertex& v : dual_vertices) { + v.ne.reserve(v.nd.size()); + } + + for (unsigned i = 0; i < dual_edges.size(); i++) { + for (unsigned vi : dual_edges[i].v) { + dual_vertices[vi].ne.push_back(i); + } + } + } class eulerException: public std::exception @@ -411,6 +421,16 @@ void graph::helper(unsigned nv, std::mt19937& rng) { } } + for (vertex& v : dual_vertices) { + v.ne.reserve(v.nd.size()); + } + + for (unsigned i = 0; i < dual_edges.size(); i++) { + for (unsigned vi : dual_edges[i].v) { + dual_vertices[vi].ne.push_back(i); + } + } + jcv_diagram_free(&diagram); } diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 2a187c5..0cf50c7 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -77,6 +77,7 @@ void network::fracture(hooks& m, bool one_axis) { void network::get_cycle_edges_helper(std::set& cycle_edges, std::set& seen_vertices, unsigned v_prev, unsigned v_cur) const { + seen_vertices.insert(v_cur); for (unsigned ei : G.dual_vertices[v_cur].ne) { if (fuses[ei]) { const std::array& e = G.dual_edges[ei].v; @@ -84,9 +85,8 @@ void network::get_cycle_edges_helper(std::set& cycle_edges, 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); + if (it != seen_vertices.end()) { + cycle_edges.insert(ei); } else { this->get_cycle_edges_helper(cycle_edges, seen_vertices, v_cur, vn); } @@ -107,16 +107,21 @@ std::set network::get_cycle_edges(unsigned v0) const { bool network::find_cycle_helper(std::array& sig, const std::set& cycle_edges, unsigned vPrev, unsigned vCur, unsigned vEnd) const { for (unsigned ei : G.dual_vertices[vCur].ne) { if (fuses[ei]) { - const std::array& e = G.dual_edges[ei].v; - unsigned vn = e[0] == vCur ? e[1] : e[0]; - if (vn != vPrev) { - if (vn == vEnd) { - return true; - } else { - if (find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) { + auto it = cycle_edges.find(ei); + if (it == cycle_edges.end()) { + const std::array& e = G.dual_edges[ei].v; + unsigned vn = e[0] == vCur ? e[1] : e[0]; + if (vn != vPrev) { + if (vn == vEnd) { if (G.dual_edges[ei].crossings[0]) sig[0]++; if (G.dual_edges[ei].crossings[1]) sig[1]++; return true; + } else { + if (this->find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) { + if (G.dual_edges[ei].crossings[0]) sig[0]++; + if (G.dual_edges[ei].crossings[1]) sig[1]++; + return true; + } } } } @@ -127,202 +132,169 @@ bool network::find_cycle_helper(std::array& sig, const std::set network::find_cycle(const std::set& cycle_edges, unsigned v0, unsigned v1) const { - std::array sig; - this->find_cycle_helper(sig, cycle_edges, v0, v0, + std::array sig = {0, 0}; + this->find_cycle_helper(sig, cycle_edges, v0, v0, v1); + return sig; } -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 {}; - - find_cycle(const graph& G, std::array& E, unsigned end) : G(G), E(E), end(end) {} +void network::get_cluster_edges_helper(std::set& seen_edges, unsigned v) const { + for (unsigned ei : G.vertices[v].ne) { + if (!backbone[ei]) { + auto it = seen_edges.find(ei); + if (it == seen_edges.end()) { + const std::array& e = G.edges[ei].v; + unsigned vn = e[0] == v ? e[1] : e[0]; - void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { - if (v == end) { - throw done{}; + seen_edges.insert(ei); + this->get_cluster_edges_helper(seen_edges, vn); } } + } +} - 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]++; - } +std::set network::get_cluster_edges(unsigned v0) const { + std::set cluster_edges; + this->get_cluster_edges_helper(cluster_edges, v0); + return cluster_edges; +} - 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 network::get_tie_flaps_helper(std::set& added_edges, unsigned v0, unsigned vCur) const { + for (unsigned ei : G.vertices[vCur].ne) { + if (!backbone[ei]) { + auto it = added_edges.find(ei); + if (it == added_edges.end()) { + const std::array& e = G.edges[ei].v; + unsigned vn = e[0] == vCur ? e[1] : e[0]; + + if (vn != v0) { + added_edges.insert(ei); + this->get_tie_flaps_helper(added_edges, v0, vn); + } + } } - }; + } +} - class get_cluster_edges : public boost::default_dfs_visitor { - public: - std::set& E; +std::list> network::get_tie_flaps(unsigned v0) const { + std::list> tie_flaps; + + for (unsigned ei : G.vertices[v0].ne) { + if (!backbone[ei]) { + bool seen_edge = false; + for (const std::set& flap : tie_flaps) { + auto it = flap.find(ei); + if (it != flap.end()) { + seen_edge = true; + break; + } + } - get_cluster_edges(std::set& E) : E(E) {} + if (!seen_edge) { + tie_flaps.push_back({ei}); - void examine_edge(boost::graph_traits::edge_descriptor ed, const Graph& g) { - unsigned e = g[ed].index; - E.insert(e); + const std::array& e = G.edges[ei].v; + unsigned vn = e[0] == v0 ? e[1] : e[0]; + + this->get_tie_flaps_helper(tie_flaps.back(), v0, vn); + } } - }; + } + + return tie_flaps; +} + +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(); + // First, we will check for lollipops! for (unsigned i = 0; i < G.edges.size(); i++) { 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 = this->get_cycle_edges(G.dual_edges[i].v[0], G.dual_edges[i].v[1]); - - std::array cycle = {0, 0}; + // This is a candidate lollipop stem. First, we will identify any + // cycles in the dual cluster that impinges on the stem and mark them + // by an edge that uniquely severs each. + std::set cedges = this->get_cycle_edges(G.dual_edges[i].v[0]); + + // Now, we create a crossing signature for each cycle. First, we do it + // for the new cycle that would form by removing the stem. + std::array base_sig = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]); if (G.dual_edges[i].crossings[0]) - cycle[0]++; + base_sig[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); - } + base_sig[1]++; - 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]; - } + // Then, we do it for each of the existing cycles we found in the + // previous step. + std::list> all_sigs = {base_sig}; + for (unsigned e : cedges) { + std::array new_sig_1 = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[e].v[0]); + std::array new_sig_2 = this->find_cycle(cedges, G.dual_edges[i].v[1], G.dual_edges[e].v[1]); + std::array new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]}; + + if (G.dual_edges[i].crossings[0]) + new_sig[0]++; + if (G.dual_edges[i].crossings[1]) + new_sig[1]++; + if (G.dual_edges[e].crossings[0]) + new_sig[0]++; + if (G.dual_edges[e].crossings[1]) + new_sig[1]++; + + all_sigs.push_back(new_sig); + } + + // Now, having found all cycles involving the candidate stem, we check + // if any of them spans the torus and therefore can carry current. + if (std::any_of(all_sigs.begin(), all_sigs.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)); + })) { + // If it does, we remove it from the backbone! + seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + backbone[i] = true; + + // We're not done yet: we've severed the stem but the pop remains! We + // check each side of the lollipop and sum up all of the currents + // that span an edge. Any component without a sufficiently large net + // current must be disconnected from the current-carrying cluster and + // can also be removed from the backbone. + for (unsigned j = 0; j < 2; j++) { + std::set cluster_edges = this->get_cluster_edges(G.edges[i].v[j]); + std::array total_current = {0.0, 0.0}; + + for (unsigned e : cluster_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_current[0] += c[e]; + } else { + total_current[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_current[1] += c[e]; + } else { + total_current[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_current[0]) < 1.0 / G.edges.size() && + fabs(total_current[1]) < 1.0 / G.edges.size()) { + for (unsigned e : cluster_edges) { + backbone[e] = true; } } } } - 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& LE; - std::set& RE; - unsigned start_edge; - bool initialized; - bool done_left; - unsigned disc; - unsigned dr; - - get_tie_edges(std::set& LE, std::set& RE) : LE(LE), RE(RE) { - initialized = false; - done_left = false; - disc = 0; - dr = 0; - } - - void discover_vertex(boost::graph_traits::vertex_descriptor vd, const Graph& g) { - if (!initialized && disc > 0) { - start_edge = vd; - initialized = true; - } - disc++; - } - - void finish_vertex(boost::graph_traits::vertex_descriptor vd, const Graph& g) { - if (vd == start_edge) { - done_left = true; - } - } - - void examine_edge(boost::graph_traits::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); - } - } - }; - + // Now, we will check for bow ties! std::set bb_verts; for (unsigned i = 0; i < G.edges.size(); i++) { @@ -385,139 +357,75 @@ void network::update_backbone(const std::vector& c) { 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::set cedges = this->get_cycle_edges(d1); - std::array cycle = {0, 0}; + std::array base_sig = this->find_cycle(cedges, d1, d2); for (unsigned k = p1; k < p2; k++) { if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) - cycle[0]++; + base_sig[0]++; if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) - cycle[1]++; + base_sig[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); + + // Then, we do it for each of the existing cycles we found in the + // previous step. + std::list> all_sigs = {base_sig}; + for (unsigned e : cedges) { + std::array new_sig_1 = this->find_cycle(cedges, d1, G.dual_edges[e].v[0]); + std::array new_sig_2 = this->find_cycle(cedges, d2, G.dual_edges[e].v[1]); + std::array new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]}; + + for (unsigned k = p1; k < p2; k++) { + if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) + new_sig[0]++; + if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) + new_sig[1]++; } + if (G.dual_edges[e].crossings[0]) + new_sig[0]++; + if (G.dual_edges[e].crossings[1]) + new_sig[1]++; - 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]; - } - } - } + all_sigs.push_back(new_sig); + } + + if (std::any_of(all_sigs.begin(), all_sigs.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::array total_cur_right = {0.0, 0.0}; - for (unsigned e : right_edges) { + std::list> flaps = this->get_tie_flaps(v); + + for (const std::set& flap : flaps) { + std::array total_current = {0.0, 0.0}; + for (unsigned e : flap) { 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]; + total_current[0] += c[e]; } else { - total_cur_right[0] -= c[e]; + total_current[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]; + total_current[1] += c[e]; } else { - total_cur_right[1] -= c[e]; + total_current[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) { + if (fabs(total_current[0]) < 1.0 / G.edges.size() && + fabs(total_current[1]) < 1.0 / G.edges.size()) { + for (unsigned e : flap) { 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); - } } } } -- cgit v1.2.3-70-g09d2