partially fixed problems that arise in small systems by reworking the way that boundary edges are identified

1 files changed, 21 insertions, 12 deletions

diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 4d4ed2d..3812f43 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -7,7 +7,7 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
     double v0y = G.vertices[G.edges[i].v[0]].r.y;
     double v1y = G.vertices[G.edges[i].v[1]].r.y;
 
-    if (fabs(v0y - v1y) > G.L.y / 2) {
+    if (G.edges[i].crossings[1]) {
       bool ind = v0y < v1y ? 0 : 1;
 
       ((double *)b->x)[G.edges[i].v[ind]] += 1.0;
@@ -19,8 +19,6 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
 
   cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz,  1, CHOLMOD_REAL, c);
 
-  t->nnz = nnz;
-
   for (unsigned int i = 0; i < G.vertices.size(); i++) {
     ((CHOL_INT *)t->i)[i] = i;
     ((CHOL_INT *)t->j)[i] = i;
@@ -29,26 +27,37 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
 
   unsigned int terms = G.vertices.size();
 
+  std::unordered_map<std::array<unsigned int, 2>, unsigned int> known_edges;
+
   for (unsigned int i = 0; i < G.edges.size(); i++) {
     unsigned int v0 = G.edges[i].v[0];
     unsigned int v1 = G.edges[i].v[1];
 
+    ((double *)t->x)[v0]++;
+    ((double *)t->x)[v1]++;
+
     unsigned int s0 = v0 < v1 ? v0 : v1;
     unsigned int s1 = v0 < v1 ? v1 : v0;
 
-    ((CHOL_INT *)t->i)[terms] = s1;
-    ((CHOL_INT *)t->j)[terms] = s0;
-    ((double *)t->x)[terms] = -1.0;
+    auto it = known_edges.find({s0, s1});
 
-    ((double *)t->x)[v0]++;
-    ((double *)t->x)[v1]++;
+    if (it == known_edges.end()) {
+      ((CHOL_INT *)t->i)[terms] = s1;
+      ((CHOL_INT *)t->j)[terms] = s0;
+      ((double *)t->x)[terms] = -1.0;
 
-    terms++;
+      known_edges[{s0, s1}] = terms;
+      terms++;
+    } else {
+      ((double *)t->x)[it->second] -= 1.0;
+    }
   }
 
   ((double *)t->x)[0]++;
 
-  cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, nnz, c);
+  t->nnz = terms;
+
+  cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c);
   CHOL_F(free_triplet)(&t, c);
   factor = CHOL_F(analyze)(laplacian, c);
   CHOL_F(factorize)(laplacian, factor, c);
@@ -138,7 +147,7 @@ void network::break_edge(unsigned int e, bool unbreak) {
   double v0y = G.vertices[v0].r.y;
   double v1y = G.vertices[v1].r.y;
 
-  if (fabs(v0y - v1y) > G.L.y / 2) {
+  if (G.edges[e].crossings[1]) {
     bool ind = v0y < v1y ? unbreak : !unbreak;
 
     ((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
@@ -172,7 +181,7 @@ current_info network::get_current_info() {
       double v0y = G.vertices[G.edges[i].v[0]].r.y;
       double v1y = G.vertices[G.edges[i].v[1]].r.y;
 
-      if (fabs(v0y - v1y) > G.L.y / 2) {
+      if (G.edges[i].crossings[1]) {
         if (v0y > v1y) {
           currents[i] += 1.0;
           total_current += currents[i];