summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/include/network.hpp8
-rw-r--r--lib/src/network.cpp78
2 files changed, 71 insertions, 15 deletions
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 29bf55a..5b474bb 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -35,9 +35,10 @@ private:
void update_backbone(const std::vector<double>& c);
void break_edge(unsigned, bool unbreak = false);
void get_cycle_edges_helper(std::set<unsigned>&, std::set<unsigned>&, unsigned, unsigned) const;
- std::set<unsigned> get_cycle_edges(unsigned) const;
bool find_cycle_helper(std::array<unsigned, 2>&, const std::set<unsigned>&, unsigned, unsigned,
unsigned) const;
+ bool get_cycle_helper(std::array<unsigned, 2>&, std::set<unsigned>&, const std::set<unsigned>&, unsigned, unsigned,
+ unsigned) const;
std::array<unsigned, 2> find_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
void get_cluster_edges_helper(std::set<unsigned>&, unsigned) const;
std::set<unsigned> get_cluster_edges(unsigned) const;
@@ -56,7 +57,10 @@ public:
network(const graph&, cholmod_common*);
void set_thresholds(double, std::mt19937&);
- void fracture(hooks&, bool one_axis = true);
+ void fracture(hooks&);
+
+ std::set<unsigned> get_cycle_edges(unsigned) const;
+ std::pair<std::array<unsigned, 2>, std::set<unsigned>> get_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
virtual current_info get_current_info() {
current_info empty;
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index eef7e9a..3e460a9 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -6,8 +6,8 @@ class nofuseException : public std::exception {
} nofuseex;
network::network(const graph& G, cholmod_common* c)
- : G(G), C(G.dual_vertices.size()), fuses(G.edges.size()), backbone(G.edges.size()),
- thresholds(G.edges.size()), px(G, 0, c), py(G, 1, 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()) {}
void network::set_thresholds(double beta, std::mt19937& rng) {
if (beta == 0.0) {
@@ -29,25 +29,22 @@ void network::set_thresholds(double beta, std::mt19937& rng) {
}
}
-void network::fracture(hooks& m, bool one_axis) {
+void network::fracture(hooks& m) {
m.pre_fracture(*this);
while (true) {
- double min_cond = 1.0 / G.edges.size();
+ const double min_cond = 1.0 / G.edges.size();
current_info c = this->get_current_info();
- if (px.sol.conductivity[0] > min_cond) {
- this->update_backbone(px.sol.currents);
- } else {
- this->update_backbone(py.sol.currents);
- }
- if (c.conductivity[0] < min_cond && c.conductivity[1] < min_cond) {
+ if (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) {
break;
}
- if (one_axis && (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond)) {
- break;
+ if (px.sol.conductivity[0] > min_cond) {
+ this->update_backbone(px.sol.currents);
+ } else {
+ this->update_backbone(py.sol.currents);
}
unsigned max_pos = UINT_MAX;
@@ -144,6 +141,48 @@ std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edge
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 {
+ for (unsigned ei : G.dual_vertices[vCur].ne) {
+ if (fuses[ei]) {
+ auto it = cycle_edges.find(ei);
+ if (it == cycle_edges.end()) {
+ const std::array<unsigned, 2>& 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::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);
+ return {sig, edges};
+}
+
void network::get_cluster_edges_helper(std::set<unsigned>& seen_edges, unsigned v) const {
for (unsigned ei : G.vertices[v].ne) {
if (!backbone[ei]) {
@@ -264,7 +303,7 @@ void network::update_backbone(const std::vector<double>& 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!
+ // If none 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;
@@ -301,6 +340,14 @@ 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));
+ })) {
+ // 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;
}
}
}
@@ -450,6 +497,11 @@ 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;
}
}
}