diff options
Diffstat (limited to 'lib/src/network.cpp')
-rw-r--r-- | lib/src/network.cpp | 83 |
1 files changed, 54 insertions, 29 deletions
diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 0cf50c7..eef7e9a 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -11,7 +11,8 @@ network::network(const graph& G, cholmod_common* c) 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" */ + /* zero beta doesn't make any sense computationally, we interpret it as + * "infinity" */ for (long double& threshold : thresholds) { threshold = 1.0; } @@ -104,7 +105,8 @@ std::set<unsigned> network::get_cycle_edges(unsigned v0) const { return cycle_edges; } -bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<unsigned>& cycle_edges, unsigned vPrev, unsigned vCur, unsigned vEnd) const { +bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<unsigned>& cycle_edges, + unsigned vPrev, unsigned vCur, unsigned vEnd) const { for (unsigned ei : G.dual_vertices[vCur].ne) { if (fuses[ei]) { auto it = cycle_edges.find(ei); @@ -113,13 +115,17 @@ bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<uns 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]++; + 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]++; + if (G.dual_edges[ei].crossings[0]) + sig[0]++; + if (G.dual_edges[ei].crossings[1]) + sig[1]++; return true; } } @@ -131,7 +137,8 @@ bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<uns return false; } -std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, unsigned v1) const { +std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, + unsigned v1) const { std::array<unsigned, 2> sig = {0, 0}; this->find_cycle_helper(sig, cycle_edges, v0, v0, v1); return sig; @@ -158,7 +165,8 @@ std::set<unsigned> network::get_cluster_edges(unsigned v0) const { return cluster_edges; } -void network::get_tie_flaps_helper(std::set<unsigned>& added_edges, unsigned v0, unsigned vCur) const { +void network::get_tie_flaps_helper(std::set<unsigned>& added_edges, unsigned v0, + unsigned vCur) const { for (unsigned ei : G.vertices[vCur].ne) { if (!backbone[ei]) { auto it = added_edges.find(ei); @@ -219,7 +227,8 @@ void network::update_backbone(const std::vector<double>& 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<unsigned, 2> base_sig = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]); + std::array<unsigned, 2> 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]) @@ -229,9 +238,12 @@ void network::update_backbone(const std::vector<double>& c) { // previous step. std::list<std::array<unsigned, 2>> all_sigs = {base_sig}; for (unsigned e : cedges) { - std::array<unsigned, 2> new_sig_1 = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[e].v[0]); - std::array<unsigned, 2> new_sig_2 = this->find_cycle(cedges, G.dual_edges[i].v[1], G.dual_edges[e].v[1]); - std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]}; + std::array<unsigned, 2> new_sig_1 = + this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[e].v[0]); + std::array<unsigned, 2> new_sig_2 = + this->find_cycle(cedges, G.dual_edges[i].v[1], G.dual_edges[e].v[1]); + std::array<unsigned, 2> 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]++; @@ -248,10 +260,10 @@ void network::update_backbone(const std::vector<double>& 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<unsigned, 2> c) { - return (c[0] % 2 == 0 && c[1] % 2 == 0) || - ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y)); - })) { + [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)); + })) { // 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; @@ -297,6 +309,7 @@ void network::update_backbone(const std::vector<double>& c) { // Now, we will check for bow ties! std::set<unsigned> 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]); @@ -307,6 +320,10 @@ void network::update_backbone(const std::vector<double>& c) { 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(); @@ -316,12 +333,11 @@ void network::update_backbone(const std::vector<double>& c) { 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; + // 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; @@ -330,7 +346,7 @@ void network::update_backbone(const std::vector<double>& c) { } 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()])) { + G.vertices[v].nd[k % G.vertices[v].nd.size()])) { any_intervening2 = true; break; } @@ -349,6 +365,9 @@ void network::update_backbone(const std::vector<double>& c) { } } } + + // 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; @@ -367,13 +386,14 @@ void network::update_backbone(const std::vector<double>& c) { base_sig[1]++; } - // Then, we do it for each of the existing cycles we found in the - // previous step. std::list<std::array<unsigned, 2>> all_sigs = {base_sig}; for (unsigned e : cedges) { - std::array<unsigned, 2> new_sig_1 = this->find_cycle(cedges, d1, G.dual_edges[e].v[0]); - std::array<unsigned, 2> new_sig_2 = this->find_cycle(cedges, d2, G.dual_edges[e].v[1]); - std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]}; + std::array<unsigned, 2> new_sig_1 = + this->find_cycle(cedges, d1, G.dual_edges[e].v[0]); + std::array<unsigned, 2> new_sig_2 = + this->find_cycle(cedges, d2, G.dual_edges[e].v[1]); + std::array<unsigned, 2> 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]) @@ -390,15 +410,20 @@ void network::update_backbone(const std::vector<double>& c) { } if (std::any_of(all_sigs.begin(), all_sigs.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); - })) { + [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); + })) { + // one of our pairs qualifies for trimming! seen_pairs[{d1, d2}] = true; + // We separate the flaps of the bowtie (there may be more than + // two!) into distinct sets. std::list<std::set<unsigned>> flaps = this->get_tie_flaps(v); + // All the bonds in each flap without current are removed from + // the backbone. for (const std::set<unsigned>& flap : flaps) { std::array<double, 2> total_current = {0.0, 0.0}; for (unsigned e : flap) { |