From 2faf0e4598c7c046d58107d23145f95db334200c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 16 Mar 2020 18:19:09 -0400 Subject: Implemented new behavior when w=0 or w=1. --- lib/src/network.cpp | 130 +++++++++++++++++++++++++--------------------------- 1 file changed, 63 insertions(+), 67 deletions(-) (limited to 'lib/src') diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 561796b..8944c15 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -1,14 +1,16 @@ #include "network.hpp" #include +#include 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) - : px(G, 0, c), py(G, 1, c), G(G), C(G.dual_vertices.size()), fuses(G.edges.size()), - backbone(G.edges.size()), thresholds(G.edges.size()) {} +network::network(const graph& G, bool two_sides, cholmod_common* c) + : px(G, 0, c), py(G, 1, c), two_sides(two_sides), G(G), + C(G.dual_vertices.size()), fuses(G.edges.size()), backbone(G.edges.size()), + thresholds(G.edges.size()) {} void network::set_thresholds(double beta, std::mt19937& rng) { if (beta == 0.0) { @@ -38,14 +40,16 @@ void network::fracture(hooks& m) { current_info c = this->get_current_info(); - if (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) { + if ((c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) && two_sides) { + break; + } else if (c.conductivity[0] < min_cond && !two_sides) { break; } if (px.sol.conductivity[0] > min_cond) { - this->update_backbone(px.sol.currents); + this->update_backbone(px.sol.currents, c.conductivity); } else { - this->update_backbone(py.sol.currents); + this->update_backbone(py.sol.currents, c.conductivity); } unsigned max_pos = UINT_MAX; @@ -74,8 +78,8 @@ void network::fracture(hooks& m) { } void network::get_cycle_edges_helper(std::set& cycle_edges, - std::set& seen_vertices, unsigned v_prev, - unsigned v_cur) const { + 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]) { @@ -103,7 +107,7 @@ 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 { + unsigned vPrev, unsigned vCur, unsigned vEnd) const { for (unsigned ei : G.dual_vertices[vCur].ne) { if (fuses[ei]) { if (!cycle_edges.contains(ei)) { @@ -134,14 +138,15 @@ bool network::find_cycle_helper(std::array& sig, const std::set network::find_cycle(const std::set& cycle_edges, unsigned v0, - unsigned v1) const { + unsigned v1) const { std::array sig = {0, 0}; this->find_cycle_helper(sig, cycle_edges, v0, v0, v1); return sig; } -bool network::get_cycle_helper(std::array& sig, std::set& edges, const std::set& cycle_edges, - unsigned vPrev, unsigned vCur, unsigned vEnd) const { +bool network::get_cycle_helper(std::array& sig, std::set& edges, + const std::set& cycle_edges, unsigned vPrev, unsigned vCur, + unsigned vEnd) const { for (unsigned ei : G.dual_vertices[vCur].ne) { if (fuses[ei]) { if (!cycle_edges.contains(ei)) { @@ -173,8 +178,8 @@ bool network::get_cycle_helper(std::array& sig, std::set& return false; } -std::pair, std::set> network::get_cycle(const std::set& cycle_edges, unsigned v0, - unsigned v1) const { +std::pair, std::set> +network::get_cycle(const std::set& cycle_edges, unsigned v0, unsigned v1) const { std::set edges; std::array sig = {0, 0}; this->get_cycle_helper(sig, edges, cycle_edges, v0, v0, v1); @@ -202,7 +207,7 @@ std::set network::get_cluster_edges(unsigned v0) const { } void network::get_tie_flaps_helper(std::set& added_edges, unsigned v0, - unsigned vCur) const { + unsigned vCur) const { for (unsigned ei : G.vertices[vCur].ne) { if (!backbone[ei]) { if (!added_edges.contains(ei)) { @@ -245,14 +250,25 @@ std::list> network::get_tie_flaps(unsigned v0) const { 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(); +void network::update_backbone(const std::vector& c, const std::array& cc) { + bool done_x = cc[0] < 1.0 / G.edges.size(); + bool done_y = cc[1] < 1.0 / G.edges.size(); + + auto cycle_check_current = [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)); + }; + auto cycle_check_wrap_x = [](std::array c) { + return (c[0] % 2 == 0 && c[1] % 2 == 1); + }; + auto cycle_check_wrap_y = [](std::array c) { + return (c[0] % 2 == 1 && c[1] % 2 == 0); + }; // 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])) { - if (!seen_pairs.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]})) { + if ((!seen_pairs_x.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}) && !done_x) || (!seen_pairs_y.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}) && !done_y)) { // 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. @@ -261,7 +277,7 @@ void network::update_backbone(const std::vector& c) { // 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]); + this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]); if (G.dual_edges[i].crossings[0]) base_sig[0]++; if (G.dual_edges[i].crossings[1]) @@ -272,11 +288,11 @@ void network::update_backbone(const std::vector& c) { 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]); + 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]); + 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]}; + new_sig_1[1] + new_sig_2[1]}; if (G.dual_edges[i].crossings[0]) new_sig[0]++; @@ -292,13 +308,10 @@ void network::update_backbone(const std::vector& c) { // 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 (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) { // If none does, we remove it from the backbone! - seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + seen_pairs_x[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + seen_pairs_y[{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 @@ -334,14 +347,13 @@ void network::update_backbone(const std::vector& c) { } } } - } else if (std::any_of(all_sigs.begin(), all_sigs.end(), - [](std::array c) { - return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0)); - })) { + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_x)) { // If the bond separates a wrapping path, breaking it would end the // fracture and therefore it will never be removed from the backbone // in this manner. We can skip it in the future. - seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + seen_pairs_x[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_y)) { + seen_pairs_y[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; } } } @@ -371,7 +383,7 @@ void network::update_backbone(const std::vector& c) { 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; - if (!seen_pairs.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]})) { + if ((!seen_pairs_x.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]}) && !done_x) || (!seen_pairs_y.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]}) && !done_y)) { bool any_intervening1 = false; bool any_intervening2 = false; unsigned ie1 = 0; @@ -449,14 +461,10 @@ void network::update_backbone(const std::vector& c) { 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); - })) { + if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) { // one of our pairs qualifies for trimming! - seen_pairs[{d1, d2}] = true; + seen_pairs_x[{d1, d2}] = true; + seen_pairs_y[{d1, d2}] = true; // We separate the flaps of the bowtie (there may be more than // two!) into distinct sets. @@ -490,11 +498,10 @@ void network::update_backbone(const std::vector& c) { } } } - } else if (std::any_of(all_sigs.begin(), all_sigs.end(), - [](std::array c) { - return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0)); - })) { - seen_pairs[{d1, d2}] = true; + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_x)) { + seen_pairs_x[{d1, d2}] = true; + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_y)) { + seen_pairs_y[{d1, d2}] = true; } } } @@ -520,7 +527,7 @@ std::string network::write() { output += "\"fuses\"->{"; for (unsigned i = 0; i < G.edges.size(); i++) { if (!fuses[i]) { - output += std::to_string(i) +","; + output += std::to_string(i) + ","; } } output.pop_back(); @@ -536,7 +543,8 @@ std::string network::write() { output += std::to_string(t) + ","; } output.pop_back(); - output += "},\"conductivity\"->{" + std::to_string(c.conductivity[0]) + "," + std::to_string(c.conductivity[1]); + output += "},\"conductivity\"->{" + std::to_string(c.conductivity[0]) + "," + + std::to_string(c.conductivity[1]); output += "},\"currents\"->{"; for (const double& t : c.currents) { output += std::to_string(t) + ","; @@ -547,35 +555,23 @@ std::string network::write() { return output; }; +fuse_network::fuse_network(const graph& G, cholmod_common* c) : network(G, false, c), weight(1.0) {} -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) {} +fuse_network::fuse_network(const fuse_network& n) : network(n), weight(1.0) {} current_info fuse_network::get_current_info() { px.solve(fuses); py.solve(fuses); bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size(); - bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size(); current_info ctot; ctot.currents.resize(G.edges.size()); ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]}; - if (done_x && !done_y) { - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = fabs(weight * py.sol.currents[i] / py.sol.conductivity[1]); - } - } else if (done_y && !done_x) { - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0]); - } - } else if (!done_x && !done_y) { + if (!done_x) { 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]); + ctot.currents[i] = fabs(px.sol.currents[i] / px.sol.conductivity[0]); } } @@ -583,7 +579,7 @@ current_info fuse_network::get_current_info() { } elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight) - : network(G, c), weight(weight) {} + : network(G, true, c), weight(weight) {} elastic_network::elastic_network(const elastic_network& n) : network(n), weight(n.weight) {} @@ -617,7 +613,7 @@ 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 graph& G, cholmod_common* c) : network(G, true, c) {} percolation_network::percolation_network(const percolation_network& n) : network(n) {} -- cgit v1.2.3-70-g09d2