summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2019-09-23 23:19:13 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2019-09-23 23:19:13 -0400
commit3f7f20f21f583ca2de566bea08a87eac4b17ad29 (patch)
treeea55b8f42f797653ed2bac21ad0bf57bf95068f5 /lib
parent9e3726c6c7762f6292333f85117546acc01f5568 (diff)
downloadfuse_networks-3f7f20f21f583ca2de566bea08a87eac4b17ad29.tar.gz
fuse_networks-3f7f20f21f583ca2de566bea08a87eac4b17ad29.tar.bz2
fuse_networks-3f7f20f21f583ca2de566bea08a87eac4b17ad29.zip
successfully implemented the backbone trimming with homespun graph algorithms, measurements have mostly been disabled and need to be migrated
Diffstat (limited to 'lib')
-rw-r--r--lib/include/network.hpp56
-rw-r--r--lib/src/graph.cpp22
-rw-r--r--lib/src/network.cpp472
3 files changed, 245 insertions, 305 deletions
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index bf9015c..8748ea9 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -18,34 +18,46 @@
#define CURRENT_CUTOFF 1e-10
+class fuse_network;
+class elastic_network;
+class percolation_network;
+
class network {
- private:
- problem px;
- problem py;
- std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+ friend class fuse_network;
+ friend class elastic_network;
+ friend class percolation_network;
- 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;
- std::array<unsigned, 2> find_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
+ private:
+ problem px;
+ problem py;
+ std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+
+ 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;
+ 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;
+ void get_tie_flaps_helper(std::set<unsigned>&, unsigned, unsigned) const;
+ std::list<std::set<unsigned>> get_tie_flaps(unsigned) const;
public:
- const graph& G;
+ const graph& G;
- ClusterTree C;
+ ClusterTree C;
- std::vector<bool> fuses;
- std::vector<bool> backbone;
- std::vector<long double> thresholds;
+ std::vector<bool> fuses;
+ std::vector<bool> backbone;
+ std::vector<long double> thresholds;
- network(const graph&, cholmod_common*);
- network(const network&);
+ network(const graph&, cholmod_common*);
- void set_thresholds(double, std::mt19937&);
- void fracture(hooks&, bool one_axis = true);
+ void set_thresholds(double, std::mt19937&);
+ void fracture(hooks&, bool one_axis = true);
- virtual current_info get_current_info() const {current_info empty; return empty;};
+ virtual current_info get_current_info() {current_info empty; return empty;};
};
class fuse_network : public network {
@@ -56,7 +68,7 @@ class fuse_network : public network {
fuse_network(const graph&, cholmod_common*, double weight = 1.0);
fuse_network(const fuse_network&);
- current_info get_current_info() const override;
+ current_info get_current_info() override;
};
class elastic_network : public network {
@@ -67,7 +79,7 @@ class elastic_network : public network {
elastic_network(const graph&, cholmod_common*, double weight = 0.5);
elastic_network(const elastic_network&);
- current_info get_current_info() const override;
+ current_info get_current_info() override;
};
class percolation_network : public network {
@@ -75,6 +87,6 @@ class percolation_network : public network {
percolation_network(const graph&, cholmod_common*);
percolation_network(const percolation_network&);
- current_info get_current_info() const override;
+ current_info get_current_info() override;
};
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index bfdb952..a5063de 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -44,7 +44,7 @@ graph::graph(unsigned Nx, unsigned Ny) {
dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
dual_vertices[i].r.y = (double)(i / (Nx / 2));
- dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx * (i / (Nx / 2)) + (i % (Nx / 2)), (nv + i - Nx / 2) % nv};
+ dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
dual_vertices[i].polygon = {
{dual_vertices[i].r.x - 1.0, vertices[i].r.y},
{dual_vertices[i].r.x, vertices[i].r.y - 1.0},
@@ -93,6 +93,16 @@ graph::graph(unsigned Nx, unsigned Ny) {
}
}
+ for (vertex& v : dual_vertices) {
+ v.ne.reserve(v.nd.size());
+ }
+
+ for (unsigned i = 0; i < dual_edges.size(); i++) {
+ for (unsigned vi : dual_edges[i].v) {
+ dual_vertices[vi].ne.push_back(i);
+ }
+ }
+
}
class eulerException: public std::exception
@@ -411,6 +421,16 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
}
}
+ for (vertex& v : dual_vertices) {
+ v.ne.reserve(v.nd.size());
+ }
+
+ for (unsigned i = 0; i < dual_edges.size(); i++) {
+ for (unsigned vi : dual_edges[i].v) {
+ dual_vertices[vi].ne.push_back(i);
+ }
+ }
+
jcv_diagram_free(&diagram);
}
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 2a187c5..0cf50c7 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -77,6 +77,7 @@ void network::fracture(hooks& m, bool one_axis) {
void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
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]) {
const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
@@ -84,9 +85,8 @@ void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
if (vn != v_prev) {
auto it = seen_vertices.find(vn);
- if (it == seen_vertices.end()) {
- // if (seen_vertices.contains()) { // need to wait for better c++20 support...
- cycle_edges.insert(vn);
+ if (it != seen_vertices.end()) {
+ cycle_edges.insert(ei);
} else {
this->get_cycle_edges_helper(cycle_edges, seen_vertices, v_cur, vn);
}
@@ -107,16 +107,21 @@ 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 {
for (unsigned ei : G.dual_vertices[vCur].ne) {
if (fuses[ei]) {
- 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) {
- return true;
- } else {
- if (find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) {
+ 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) {
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;
+ }
}
}
}
@@ -127,202 +132,169 @@ 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 {
- std::array<unsigned, 2> sig;
- this->find_cycle_helper(sig, cycle_edges, v0, v0,
+ std::array<unsigned, 2> sig = {0, 0};
+ this->find_cycle_helper(sig, cycle_edges, v0, v0, v1);
+ return sig;
}
-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();
-
- class find_cycle : public boost::default_dfs_visitor {
- public:
- const graph& G;
- std::array<unsigned, 2>& E;
- unsigned end;
- struct done {};
-
- find_cycle(const graph& G, std::array<unsigned, 2>& E, unsigned end) : G(G), E(E), end(end) {}
+void network::get_cluster_edges_helper(std::set<unsigned>& seen_edges, unsigned v) const {
+ for (unsigned ei : G.vertices[v].ne) {
+ if (!backbone[ei]) {
+ auto it = seen_edges.find(ei);
+ if (it == seen_edges.end()) {
+ const std::array<unsigned, 2>& e = G.edges[ei].v;
+ unsigned vn = e[0] == v ? e[1] : e[0];
- void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
- if (v == end) {
- throw done{};
+ seen_edges.insert(ei);
+ this->get_cluster_edges_helper(seen_edges, vn);
}
}
+ }
+}
- void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
- unsigned e = g[ed].index;
- if (G.dual_edges[e].crossings[0])
- E[0]++;
- if (G.dual_edges[e].crossings[1])
- E[1]++;
- }
+std::set<unsigned> network::get_cluster_edges(unsigned v0) const {
+ std::set<unsigned> cluster_edges;
+ this->get_cluster_edges_helper(cluster_edges, v0);
+ return cluster_edges;
+}
- void finish_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
- unsigned e = g[ed].index;
- if (G.dual_edges[e].crossings[0])
- E[0]--;
- if (G.dual_edges[e].crossings[1])
- E[1]--;
+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);
+ if (it == added_edges.end()) {
+ const std::array<unsigned, 2>& 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);
+ }
+ }
}
- };
+ }
+}
- class get_cluster_edges : public boost::default_dfs_visitor {
- public:
- std::set<unsigned>& E;
+std::list<std::set<unsigned>> network::get_tie_flaps(unsigned v0) const {
+ std::list<std::set<unsigned>> tie_flaps;
+
+ for (unsigned ei : G.vertices[v0].ne) {
+ if (!backbone[ei]) {
+ bool seen_edge = false;
+ for (const std::set<unsigned>& flap : tie_flaps) {
+ auto it = flap.find(ei);
+ if (it != flap.end()) {
+ seen_edge = true;
+ break;
+ }
+ }
- get_cluster_edges(std::set<unsigned>& E) : E(E) {}
+ if (!seen_edge) {
+ tie_flaps.push_back({ei});
- void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
- unsigned e = g[ed].index;
- E.insert(e);
+ const std::array<unsigned, 2>& 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<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();
+ // 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])) {
auto it_search = seen_pairs.find({G.dual_edges[i].v[0], G.dual_edges[i].v[1]});
if (it_search == seen_pairs.end()) {
- std::list<unsigned> cedges = this->get_cycle_edges(G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
-
- std::array<unsigned, 2> cycle = {0, 0};
+ // 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<unsigned> 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<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])
- cycle[0]++;
+ base_sig[0]++;
if (G.dual_edges[i].crossings[1])
- cycle[1]++;
- find_cycle vis2(G, cycle, G.dual_edges[i].v[1]);
- std::vector<boost::default_color_type> cm2(G.dual_vertices.size());
- try {
- boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis2,
- boost::make_iterator_property_map(
- cm2.begin(), boost::get(boost::vertex_index, dG), cm2[0]));
- } catch (find_cycle::done const&) {
- std::list<std::array<unsigned, 2>> other_cycles = {cycle};
- for (unsigned e : cedges) {
- std::array<unsigned, 2> tcycle = {0, 0};
- if (G.dual_edges[i].crossings[0])
- tcycle[0]++;
- if (G.dual_edges[i].crossings[1])
- tcycle[1]++;
- if (G.dual_edges[e].crossings[0])
- tcycle[0]++;
- if (G.dual_edges[e].crossings[1])
- tcycle[1]++;
- find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]);
- std::vector<boost::default_color_type> cm3(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, G.dual_edges[i].v[0], vis3,
- boost::make_iterator_property_map(cm3.begin(),
- boost::get(boost::vertex_index, dG), cm3[0]));
- } catch (find_cycle::done const&) {
- }
- find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]);
- std::vector<boost::default_color_type> cm4(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, G.dual_edges[i].v[1], vis4,
- boost::make_iterator_property_map(cm4.begin(),
- boost::get(boost::vertex_index, dG), cm4[0]));
- } catch (find_cycle::done const&) {
- }
- other_cycles.push_back(tcycle);
- }
+ base_sig[1]++;
- if (std::any_of(other_cycles.begin(), other_cycles.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));
- })) {
- seen_pairs[{G.edges[i].v[0], G.edges[i].v[1]}] = true;
- backbone[i] = true;
- boost::remove_edge(G.edges[i].v[0], G.edges[i].v[1], bG);
-
- for (unsigned j = 0; j < 2; j++) {
- std::set<unsigned> component_edges = {};
- get_cluster_edges visc(component_edges);
- std::vector<boost::default_color_type> cm2(G.vertices.size());
- boost::depth_first_visit(
- bG, G.edges[i].v[j], visc,
- boost::make_iterator_property_map(cm2.begin(),
- boost::get(boost::vertex_index, bG), cm2[0]));
-
- std::array<double, 2> total_cur = {0.0, 0.0};
- for (unsigned e : component_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_cur[0] += c[e];
- } else {
- total_cur[0] -= c[e];
- }
+ // 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, 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]++;
+ 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(),
+ [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;
+
+ // 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<unsigned> cluster_edges = this->get_cluster_edges(G.edges[i].v[j]);
+ std::array<double, 2> 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_cur[1] += c[e];
- } else {
- total_cur[1] -= 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_cur[0]) < 1.0 / G.edges.size() &&
- fabs(total_cur[1]) < 1.0 / G.edges.size()) {
- for (unsigned e : component_edges) {
- backbone[e] = true;
- boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
- }
+ 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;
}
}
}
}
- for (unsigned e : cedges) {
- boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
- }
}
}
}
- class get_tie_edges : public boost::default_dfs_visitor {
- public:
- std::set<unsigned>& LE;
- std::set<unsigned>& RE;
- unsigned start_edge;
- bool initialized;
- bool done_left;
- unsigned disc;
- unsigned dr;
-
- get_tie_edges(std::set<unsigned>& LE, std::set<unsigned>& RE) : LE(LE), RE(RE) {
- initialized = false;
- done_left = false;
- disc = 0;
- dr = 0;
- }
-
- void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
- if (!initialized && disc > 0) {
- start_edge = vd;
- initialized = true;
- }
- disc++;
- }
-
- void finish_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
- if (vd == start_edge) {
- done_left = true;
- }
- }
-
- void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
- unsigned e = g[ed].index;
- if (!done_left) {
- LE.insert(e);
- } else if (LE.find(e) == LE.end()) {
- RE.insert(e);
- }
- }
- };
-
+ // Now, we will check for bow ties!
std::set<unsigned> bb_verts;
for (unsigned i = 0; i < G.edges.size(); i++) {
@@ -385,139 +357,75 @@ void network::update_backbone(const std::vector<double>& c) {
d1 = G.vertices[v].nd[il];
d2 = G.vertices[v].nd[ig];
- std::list<unsigned> cedges = {};
- get_cycle_edges vis(cedges);
- std::vector<boost::default_color_type> cm1(G.dual_vertices.size());
- boost::depth_first_visit(
- dG, d1, vis,
- boost::make_iterator_property_map(cm1.begin(),
- boost::get(boost::vertex_index, dG), cm1[0]));
-
- for (unsigned e : cedges) {
- boost::remove_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], dG);
- }
+ std::set<unsigned> cedges = this->get_cycle_edges(d1);
- std::array<unsigned, 2> cycle = {0, 0};
+ std::array<unsigned, 2> 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])
- cycle[0]++;
+ base_sig[0]++;
if (G.dual_edges[G.vertices[v].ne[k]].crossings[1])
- cycle[1]++;
+ base_sig[1]++;
}
- find_cycle vis2(G, cycle, d2);
- std::vector<boost::default_color_type> cm2(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, d1, vis2,
- boost::make_iterator_property_map(cm2.begin(),
- boost::get(boost::vertex_index, dG), cm2[0]));
- } catch (find_cycle::done const&) {
- std::list<std::array<unsigned, 2>> other_cycles = {cycle};
- for (unsigned e : cedges) {
- std::array<unsigned, 2> tcycle = {0, 0};
- for (unsigned k = p1; k < p2; k++) {
- if (G.dual_edges[G.vertices[v].ne[k]].crossings[0])
- tcycle[0]++;
- if (G.dual_edges[G.vertices[v].ne[k]].crossings[1])
- tcycle[1]++;
- }
- if (G.dual_edges[e].crossings[0])
- tcycle[0]++;
- if (G.dual_edges[e].crossings[1])
- tcycle[1]++;
- find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]);
- std::vector<boost::default_color_type> cm3(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, d1, vis3,
- boost::make_iterator_property_map(
- cm3.begin(), boost::get(boost::vertex_index, dG), cm3[0]));
- } catch (find_cycle::done const&) {
- }
- find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]);
- std::vector<boost::default_color_type> cm4(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, d2, vis4,
- boost::make_iterator_property_map(
- cm4.begin(), boost::get(boost::vertex_index, dG), cm4[0]));
- } catch (find_cycle::done const&) {
- }
- other_cycles.push_back(tcycle);
+
+ // 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]};
+
+ 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]++;
- if (std::any_of(other_cycles.begin(), other_cycles.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);
- })) {
- seen_pairs[{d1, d2}] = true;
-
- std::set<unsigned> left_edges = {};
- std::set<unsigned> right_edges = {};
- get_tie_edges visc(left_edges, right_edges);
- std::vector<boost::default_color_type> cm2(G.vertices.size());
- boost::depth_first_visit(
- bG, v, visc,
- boost::make_iterator_property_map(
- cm2.begin(), boost::get(boost::vertex_index, bG), cm2[0]));
-
- std::array<double, 2> total_cur_left = {0.0, 0.0};
- for (unsigned e : left_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_cur_left[0] += c[e];
- } else {
- total_cur_left[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_cur_left[1] += c[e];
- } else {
- total_cur_left[1] -= c[e];
- }
- }
- }
+ 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);
+ })) {
+ seen_pairs[{d1, d2}] = true;
- std::array<double, 2> total_cur_right = {0.0, 0.0};
- for (unsigned e : right_edges) {
+ std::list<std::set<unsigned>> flaps = this->get_tie_flaps(v);
+
+ for (const std::set<unsigned>& flap : flaps) {
+ std::array<double, 2> 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_cur_right[0] += c[e];
+ total_current[0] += c[e];
} else {
- total_cur_right[0] -= c[e];
+ 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_cur_right[1] += c[e];
+ total_current[1] += c[e];
} else {
- total_cur_right[1] -= c[e];
+ total_current[1] -= c[e];
}
}
}
- if (fabs(total_cur_left[0]) < 1.0 / G.edges.size() &&
- fabs(total_cur_left[1]) < 1.0 / G.edges.size()) {
- for (unsigned e : left_edges) {
+ 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;
- boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
- }
- }
- if (fabs(total_cur_right[0]) < 1.0 / G.edges.size() &&
- fabs(total_cur_right[1]) < 1.0 / G.edges.size()) {
- for (unsigned e : right_edges) {
- backbone[e] = true;
- boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
}
}
}
}
- for (unsigned e : cedges) {
- boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
- }
}
}
}