summaryrefslogtreecommitdiff
path: root/lib/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src')
-rw-r--r--lib/src/network.cpp130
1 files changed, 63 insertions, 67 deletions
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 <sstream>
+#include <iostream>
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<unsigned>& cycle_edges,
- std::set<unsigned>& seen_vertices, unsigned v_prev,
- unsigned v_cur) const {
+ std::set<unsigned>& 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<unsigned> network::get_cycle_edges(unsigned v0) const {
}
bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<unsigned>& 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<unsigned, 2>& sig, const std::set<uns
}
std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0,
- unsigned v1) const {
+ unsigned v1) const {
std::array<unsigned, 2> sig = {0, 0};
this->find_cycle_helper(sig, cycle_edges, v0, v0, v1);
return sig;
}
-bool network::get_cycle_helper(std::array<unsigned, 2>& sig, std::set<unsigned>& edges, const std::set<unsigned>& cycle_edges,
- unsigned vPrev, unsigned vCur, unsigned vEnd) const {
+bool network::get_cycle_helper(std::array<unsigned, 2>& sig, std::set<unsigned>& edges,
+ const std::set<unsigned>& 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<unsigned, 2>& sig, std::set<unsigned>&
return false;
}
-std::pair<std::array<unsigned, 2>, std::set<unsigned>> network::get_cycle(const std::set<unsigned>& cycle_edges, unsigned v0,
- unsigned v1) const {
+std::pair<std::array<unsigned, 2>, std::set<unsigned>>
+network::get_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, unsigned v1) const {
std::set<unsigned> edges;
std::array<unsigned, 2> sig = {0, 0};
this->get_cycle_helper(sig, edges, cycle_edges, v0, v0, v1);
@@ -202,7 +207,7 @@ std::set<unsigned> network::get_cluster_edges(unsigned v0) const {
}
void network::get_tie_flaps_helper(std::set<unsigned>& 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<std::set<unsigned>> network::get_tie_flaps(unsigned v0) const {
return tie_flaps;
}
-void network::update_backbone(const std::vector<double>& 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<double>& c, const std::array<double, 2>& 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<unsigned, 2> 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<unsigned, 2> c) {
+ return (c[0] % 2 == 0 && c[1] % 2 == 1);
+ };
+ auto cycle_check_wrap_y = [](std::array<unsigned, 2> 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<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]);
+ 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<double>& c) {
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]);
+ 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]);
+ 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]};
+ 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<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));
- })) {
+ 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<double>& c) {
}
}
}
- } else if (std::any_of(all_sigs.begin(), all_sigs.end(),
- [](std::array<unsigned, 2> 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<double>& 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<double>& c) {
all_sigs.push_back(new_sig);
}
- 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);
- })) {
+ 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<double>& c) {
}
}
}
- } else if (std::any_of(all_sigs.begin(), all_sigs.end(),
- [](std::array<unsigned, 2> 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) {}