#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, 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) { /* zero beta doesn't make any sense computationally, we interpret it as * "infinity" */ for (long double& threshold : thresholds) { threshold = 1.0; } } else { std::uniform_real_distribution d(0.0, 1.0); for (long double& threshold : thresholds) { threshold = std::numeric_limits::lowest(); while (threshold == std::numeric_limits::lowest()) { threshold = logl(d(rng)) / (long double)beta; } } } } void network::fracture(hooks& m) { m.pre_fracture(*this); while (true) { const double min_cond = 1.0 / G.edges.size(); current_info c = this->get_current_info(); 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, c.conductivity); } else { this->update_backbone(py.sol.currents, c.conductivity); } unsigned max_pos = UINT_MAX; long double max_val = std::numeric_limits::lowest(); for (unsigned i = 0; i < G.edges.size(); i++) { if (!backbone[i]) { long double val = logl(c.currents[i]) - thresholds[i]; if (val > max_val) { max_val = val; max_pos = i; } } } if (max_pos == UINT_MAX) { throw nofuseex; } m.bond_broken(*this, c, max_pos); this->break_edge(max_pos); } m.post_fracture(*this); } 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; unsigned vn = e[0] == v_cur ? e[1] : e[0]; if (vn != v_prev) { if (seen_vertices.contains(vn)) { cycle_edges.insert(ei); } else { this->get_cycle_edges_helper(cycle_edges, seen_vertices, v_cur, vn); } } } } } std::set network::get_cycle_edges(unsigned v0) const { std::set seen_vertices; std::set cycle_edges; this->get_cycle_edges_helper(cycle_edges, seen_vertices, v0, v0); return cycle_edges; } 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]) { if (!cycle_edges.contains(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) { 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; } } } } } } return false; } std::array network::find_cycle(const std::set& cycle_edges, unsigned v0, 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 { for (unsigned ei : G.dual_vertices[vCur].ne) { if (fuses[ei]) { if (!cycle_edges.contains(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) { edges.insert(ei); if (G.dual_edges[ei].crossings[0]) sig[0]++; if (G.dual_edges[ei].crossings[1]) sig[1]++; return true; } else { if (this->get_cycle_helper(sig, edges, cycle_edges, vCur, vn, vEnd)) { edges.insert(ei); if (G.dual_edges[ei].crossings[0]) sig[0]++; if (G.dual_edges[ei].crossings[1]) sig[1]++; return true; } } } } } } return false; } 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); return {sig, edges}; } void network::get_cluster_edges_helper(std::set& seen_edges, unsigned v) const { for (unsigned ei : G.vertices[v].ne) { if (!backbone[ei]) { if (!seen_edges.contains(ei)) { const std::array& e = G.edges[ei].v; unsigned vn = e[0] == v ? e[1] : e[0]; seen_edges.insert(ei); this->get_cluster_edges_helper(seen_edges, vn); } } } } 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 network::get_tie_flaps_helper(std::set& added_edges, unsigned v0, unsigned vCur) const { for (unsigned ei : G.vertices[vCur].ne) { if (!backbone[ei]) { if (!added_edges.contains(ei)) { 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); } } } } } 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) { if (flap.contains(ei)) { seen_edge = true; break; } } if (!seen_edge) { tie_flaps.push_back({ei}); 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, 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_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. 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]) base_sig[0]++; if (G.dual_edges[i].crossings[1]) base_sig[1]++; // 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(), cycle_check_current)) { // If none does, we remove it from the backbone! 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 // 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_current[1] += c[e]; } else { total_current[1] -= c[e]; } } } 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; } } } } 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_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; } } } } // Now, we will check for bow ties! std::set bb_verts; // First we get a list of all vertices that remain in the backbone. for (unsigned i = 0; i < G.edges.size(); i++) { if (!backbone[i]) { bb_verts.insert(G.edges[i].v[0]); bb_verts.insert(G.edges[i].v[1]); } } for (unsigned v : bb_verts) { bool found_pair = false; unsigned d1, d2, p1, p2; // For each vertex, we check to see if two of the faces adjacent to the // vertex are in the same component of the dual lattice and if so, we make // sure they do not belong to a contiguously connected series of such // faces. 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 (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_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; unsigned ie2 = 0; // yuck, think of something more elegant? 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; } } 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_intervening2 && !any_intervening1) { for (unsigned k = il + 1; k <= ig; k++) { if (!backbone[G.vertices[v].ne[k]]) { ie1++; } } } 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 all these conditions are true, we process the pair. The same // cycle analysis as in the lollipop case is done. 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::set cedges = this->get_cycle_edges(d1); 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]) base_sig[0]++; if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) base_sig[1]++; } 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]++; all_sigs.push_back(new_sig); } if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) { // one of our pairs qualifies for trimming! 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. std::list> flaps = this->get_tie_flaps(v); // All the bonds in each flap without current are removed from // the backbone. 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_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_current[1] += c[e]; } else { total_current[1] -= c[e]; } } } 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; } } } } 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; } } } } } } } } void network::break_edge(unsigned e, bool unbreak) { fuses[e] = !unbreak; backbone[e] = !unbreak; C.add_bond(G.dual_edges[e]); px.break_edge(e, unbreak); py.break_edge(e, unbreak); } std::string network::write() { std::string output; current_info c = this->get_current_info(); output += "\"fuses\"->{"; for (unsigned i = 0; i < G.edges.size(); i++) { if (!fuses[i]) { output += std::to_string(i) + ","; } } output.pop_back(); output += "},\"backbone\"->{"; for (unsigned i = 0; i < G.edges.size(); i++) { if (!backbone[i]) { output += std::to_string(i) + ","; } } output.pop_back(); output += "},\"thresholds\"->{"; for (const long double& t : thresholds) { output += std::to_string(t) + ","; } output.pop_back(); 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) + ","; } output.pop_back(); output += "}," + G.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 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(); current_info ctot; ctot.currents.resize(G.edges.size()); ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]}; if (!done_x) { for (unsigned i = 0; i < G.edges.size(); i++) { ctot.currents[i] = fabs(px.sol.currents[i] / px.sol.conductivity[0]); } } return ctot; } elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight) : network(G, true, c), weight(weight) {} elastic_network::elastic_network(const elastic_network& n) : network(n), weight(n.weight) {} current_info elastic_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[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++) { ctot.currents[i] = weight * fabs(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] = (1 - weight) * fabs(px.sol.currents[i]) / px.sol.conductivity[0]; } } else if (!done_x && !done_y) { for (unsigned i = 0; i < G.edges.size(); i++) { ctot.currents[i] = sqrt(pow((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0], 2) + pow(weight * py.sol.currents[i] / py.sol.conductivity[1], 2)); } } return ctot; } percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, true, c) {} percolation_network::percolation_network(const percolation_network& n) : network(n) {} current_info percolation_network::get_current_info() { current_info ctot; ctot.currents.resize(G.edges.size(), 0.0); px.solve(fuses); py.solve(fuses); ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]}; for (unsigned i = 0; i < G.edges.size(); i++) { if (fabs(px.sol.currents[i]) > CURRENT_CUTOFF || fabs(py.sol.currents[i]) > CURRENT_CUTOFF) { ctot.currents[i] = 1.0; } } return ctot; }