ran clang-format

3 files changed, 170 insertions, 140 deletions

diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index a5063de..1475c40 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -34,23 +34,21 @@ graph::graph(unsigned Nx, unsigned Ny) {
     vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
     vertices[i].r.y = (double)(i / (Nx / 2));
     signed f = (1 + i / (Nx / 2)) % 2 == 1 ? 1 : -1;
-    vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
-    vertices[i].polygon = {
-      {vertices[i].r.x - 1.0, vertices[i].r.y},
-      {vertices[i].r.x, vertices[i].r.y - 1.0},
-      {vertices[i].r.x + 1.0, vertices[i].r.y},
-      {vertices[i].r.x, vertices[i].r.y + 1.0}
-    };
+    vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)),
+                      (nv + i - Nx / 2) % nv};
+    vertices[i].polygon = {{vertices[i].r.x - 1.0, vertices[i].r.y},
+                           {vertices[i].r.x, vertices[i].r.y - 1.0},
+                           {vertices[i].r.x + 1.0, vertices[i].r.y},
+                           {vertices[i].r.x, vertices[i].r.y + 1.0}};
 
     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 / 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},
-      {dual_vertices[i].r.x + 1.0, vertices[i].r.y},
-      {dual_vertices[i].r.x, vertices[i].r.y + 1.0}
-    };
+    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},
+                                {dual_vertices[i].r.x + 1.0, vertices[i].r.y},
+                                {dual_vertices[i].r.x, vertices[i].r.y + 1.0}};
   }
 
   for (unsigned y = 0; y < Ny; y++) {
@@ -66,7 +64,8 @@ graph::graph(unsigned Nx, unsigned Ny) {
       unsigned dv1 = (Nx * y) / 2 + ((x + (y + 1) % 2) / 2) % (Nx / 2);
       unsigned dv2 = ((Nx * (y + 1)) / 2 + ((x + y % 2) / 2) % (Nx / 2)) % nv;
 
-      dual_edges.push_back({{dv1, dv2}, {0.5 + (double)x, 0.5 + (double)y}, {crossed_x, crossed_y}});
+      dual_edges.push_back(
+          {{dv1, dv2}, {0.5 + (double)x, 0.5 + (double)y}, {crossed_x, crossed_y}});
     }
   }
 
@@ -102,29 +101,23 @@ graph::graph(unsigned Nx, unsigned Ny) {
       dual_vertices[vi].ne.push_back(i);
     }
   }
-
 }
 
-class eulerException: public std::exception
-{
-  virtual const char* what() const throw()
-  {
-    return "The voronoi tessellation generated has the incorrect Euler characteristic for a torus and is malformed.";
+class eulerException : public std::exception {
+  virtual const char* what() const throw() {
+    return "The voronoi tessellation generated has the incorrect Euler characteristic for a torus "
+           "and is malformed.";
   }
 } eulerex;
 
-class clippingException: public std::exception
-{
-  virtual const char* what() const throw()
-  {
+class clippingException : public std::exception {
+  virtual const char* what() const throw() {
     return "An interior site has a clipped edge and the tesselation is malformed.";
   }
 } clippingex;
 
-class triangleException: public std::exception
-{
-  virtual const char* what() const throw()
-  {
+class triangleException : public std::exception {
+  virtual const char* what() const throw() {
     return "A dual-centered triangle has an impossible geometry and the tesselation is malformed.";
   }
 } triex;
@@ -147,29 +140,45 @@ unsigned get_triangle_signature(unsigned j1, unsigned j2, unsigned j3) {
 
   if ((j1 == j2) && (j2 == j3)) {
     return 0;
-  } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) || ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
+  } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) ||
+             ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
     return 1;
-  } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) || ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
+  } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) ||
+             ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
     return 2;
-  } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) || ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
+  } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) ||
+             ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
     return 3;
-  } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) || ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
+  } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) ||
+             ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
     return 4;
-  } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) || ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
+  } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) ||
+             ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
     return 5;
-  } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) || ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
+  } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) ||
+             ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
     return 6;
-  } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) || ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
+  } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) ||
+             ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
     return 7;
-  } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) || ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
+  } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) ||
+             ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
     return 8;
-  } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
+  } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) ||
+             ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) ||
+             ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
     return 9;
-  } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
+  } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) ||
+             ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) ||
+             ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
     return 10;
-  } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
+  } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) ||
+             ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) ||
+             ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
     return 11;
-  } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
+  } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) ||
+             ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) ||
+             ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
     return 12;
   } else {
     throw triex;
@@ -200,7 +209,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
   // the coordinates of the lattice, from which a delaunay triangulation
   // and corresponding voronoi tessellation will be built. Everyone is in the
   // rectangle (0, 0) (Lx, Ly)
-  for (vertex &v : vertices) {
+  for (vertex& v : vertices) {
     v = {{L.x * d(rng), L.y * d(rng)}};
   }
 
@@ -243,7 +252,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
         ep = ep->next;
       }
       // for each edge on the site's cell
-      while(e) {
+      while (e) {
         // assess whether the edge needs to be added. only neighboring
         // sites whose index is larger than ours are given edges, so we
         // don't double up later
@@ -275,7 +284,8 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
 
         if (it1 == known_vertices.end()) {
           vi1 = dual_vertices.size();
-          dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{i1},{},{{site->p.x, site->p.y}}});
+          dual_vertices.push_back(
+              {{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)}, {i1}, {}, {{site->p.x, site->p.y}}});
           known_vertices[t1] = vi1;
         } else {
           vi1 = it1->second;
@@ -284,7 +294,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
         }
 
         vertices[i1].nd.push_back(vi1);
-        
+
         if (i1 < i2 || (i1 == i2 && !self_bonded)) {
           if (i1 == i2) {
             self_bonded = true;
@@ -292,12 +302,11 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
           bool crossed_x = x2 != 1;
           bool crossed_y = y2 != 1;
           edges.push_back({{i1, i2},
-              {mod((site->p.x + neighbor->p.x) / 2, L.x),
-               mod((site->p.y + neighbor->p.y) / 2, L.y)},
-               {crossed_x, crossed_y}
-               });
+                           {mod((site->p.x + neighbor->p.x) / 2, L.x),
+                            mod((site->p.y + neighbor->p.y) / 2, L.y)},
+                           {crossed_x, crossed_y}});
 
-          jcv_graphedge *en;
+          jcv_graphedge* en;
           if (e->next == NULL) {
             en = site->edges;
           } else {
@@ -322,20 +331,21 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
 
           if (it2 == known_vertices.end()) {
             vi2 = dual_vertices.size();
-            dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)},{}});
+            dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)}, {}});
             known_vertices[t2] = vi2;
           } else {
             vi2 = it2->second;
           }
 
-          bool dcrossed_x = (unsigned)floor(e->pos[0].x / L.x) != (unsigned)floor(e->pos[1].x / L.x);
-          bool dcrossed_y = (unsigned)floor(e->pos[0].y / L.y) != (unsigned)floor(e->pos[1].y / L.y);
+          bool dcrossed_x =
+              (unsigned)floor(e->pos[0].x / L.x) != (unsigned)floor(e->pos[1].x / L.x);
+          bool dcrossed_y =
+              (unsigned)floor(e->pos[0].y / L.y) != (unsigned)floor(e->pos[1].y / L.y);
 
           dual_edges.push_back({{vi1, vi2},
-              {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
-               mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
-               {dcrossed_x, dcrossed_y}
-               });
+                                {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
+                                 mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
+                                {dcrossed_x, dcrossed_y}});
         }
 
         ep = e;
@@ -348,7 +358,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
     throw eulerex;
   }
 
-  for (vertex &v : dual_vertices) {
+  for (vertex& v : dual_vertices) {
     if (fabs(v.polygon[0].x - v.polygon[1].x) > L.x / 2) {
       if (v.polygon[0].x < L.x / 2) {
         v.polygon[0].x += L.x;
@@ -434,33 +444,31 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
   jcv_diagram_free(&diagram);
 }
 
-graph::coordinate reverse(const graph::coordinate &x) {
-  return {x.y, x.x};
-}
+graph::coordinate reverse(const graph::coordinate& x) { return {x.y, x.x}; }
 
 graph const graph::rotate() {
   graph g2(*this);
 
-  for (graph::vertex &v : g2.vertices) {
+  for (graph::vertex& v : g2.vertices) {
     v.r = reverse(v.r);
-    for (graph::coordinate &r : v.polygon) {
+    for (graph::coordinate& r : v.polygon) {
       r = reverse(r);
     }
   }
 
-  for (graph::edge &e : g2.edges) {
+  for (graph::edge& e : g2.edges) {
     e.r = reverse(e.r);
     e.crossings = {e.crossings[1], e.crossings[0]};
   }
 
-  for (graph::vertex &v : g2.dual_vertices) {
+  for (graph::vertex& v : g2.dual_vertices) {
     v.r = reverse(v.r);
-    for (graph::coordinate &r : v.polygon) {
+    for (graph::coordinate& r : v.polygon) {
       r = reverse(r);
     }
   }
 
-  for (graph::edge &e : g2.dual_edges) {
+  for (graph::edge& e : g2.dual_edges) {
     e.r = reverse(e.r);
     e.crossings = {e.crossings[1], e.crossings[0]};
   }
@@ -469,4 +477,3 @@ graph const graph::rotate() {
 
   return g2;
 }
-
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) {
diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp
index 0645244..ab87b0c 100644
--- a/lib/src/problem.cpp
+++ b/lib/src/problem.cpp
@@ -1,15 +1,12 @@
 
 #include "problem.hpp"
 
-class nanException: public std::exception
-{
-  virtual const char* what() const throw()
-  {
-    return "The linear problem returned NaN.";
-  }
+class nanException : public std::exception {
+  virtual const char* what() const throw() { return "The linear problem returned NaN."; }
 } nanex;
 
-problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), axis(axis), voltcurmat(vcmat), c(c) {
+problem::problem(const graph& G, unsigned axis, cholmod_sparse* vcmat, cholmod_common* c)
+    : G(G), axis(axis), voltcurmat(vcmat), c(c) {
   b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
   for (unsigned i = 0; i < G.edges.size(); i++) {
     graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
@@ -23,19 +20,20 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
         ind = v0.x < v1.x ? 0 : 1;
       }
 
-      ((double *)b->x)[G.edges[i].v[ind]]  += 1.0;
-      ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0;
+      ((double*)b->x)[G.edges[i].v[ind]] += 1.0;
+      ((double*)b->x)[G.edges[i].v[!ind]] -= 1.0;
     }
   }
 
   unsigned nnz = G.vertices.size() + G.edges.size();
 
-  cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz,  1, CHOLMOD_REAL, c);
+  cholmod_triplet* t =
+      CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c);
 
   for (unsigned i = 0; i < G.vertices.size(); i++) {
-    ((CHOL_INT *)t->i)[i] = i;
-    ((CHOL_INT *)t->j)[i] = i;
-    ((double *)t->x)[i] = 0.0;
+    ((CHOL_INT*)t->i)[i] = i;
+    ((CHOL_INT*)t->j)[i] = i;
+    ((double*)t->x)[i] = 0.0;
   }
 
   unsigned terms = G.vertices.size();
@@ -46,8 +44,8 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
     unsigned v0 = G.edges[i].v[0];
     unsigned v1 = G.edges[i].v[1];
 
-    ((double *)t->x)[v0]++;
-    ((double *)t->x)[v1]++;
+    ((double*)t->x)[v0]++;
+    ((double*)t->x)[v1]++;
 
     unsigned s0 = v0 < v1 ? v0 : v1;
     unsigned s1 = v0 < v1 ? v1 : v0;
@@ -55,22 +53,22 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
     auto it = known_edges.find({s0, s1});
 
     if (it == known_edges.end()) {
-      ((CHOL_INT *)t->i)[terms] = s1;
-      ((CHOL_INT *)t->j)[terms] = s0;
-      ((double *)t->x)[terms] = -1.0;
+      ((CHOL_INT*)t->i)[terms] = s1;
+      ((CHOL_INT*)t->j)[terms] = s0;
+      ((double*)t->x)[terms] = -1.0;
 
       known_edges[{s0, s1}] = terms;
       terms++;
     } else {
-      ((double *)t->x)[it->second] -= 1.0;
+      ((double*)t->x)[it->second] -= 1.0;
     }
   }
 
-  ((double *)t->x)[0]++;
+  ((double*)t->x)[0]++;
 
   t->nnz = terms;
 
-  cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c);
+  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);
@@ -79,19 +77,20 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
   sol.currents.resize(G.edges.size());
 }
 
-problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) {
-  cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
+problem::problem(const graph& G, unsigned axis, cholmod_common* c) : problem(G, axis, NULL, c) {
+  cholmod_triplet* t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(),
+                                                2 * G.edges.size(), 0, CHOLMOD_REAL, c);
 
   t->nnz = 2 * G.edges.size();
 
   for (unsigned i = 0; i < G.edges.size(); i++) {
-    ((CHOL_INT *)t->i)[2 * i] = i;
-    ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0];
-    ((double *)t->x)[2 * i] = 1.0;
+    ((CHOL_INT*)t->i)[2 * i] = i;
+    ((CHOL_INT*)t->j)[2 * i] = G.edges[i].v[0];
+    ((double*)t->x)[2 * i] = 1.0;
 
-    ((CHOL_INT *)t->i)[2 * i + 1] = i;
-    ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1];
-    ((double *)t->x)[2 * i + 1] = -1.0;
+    ((CHOL_INT*)t->i)[2 * i + 1] = i;
+    ((CHOL_INT*)t->j)[2 * i + 1] = G.edges[i].v[1];
+    ((double*)t->x)[2 * i + 1] = -1.0;
   }
 
   voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c);
@@ -103,7 +102,7 @@ problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c
   b = CHOL_F(copy_dense)(other.b, c);
   factor = CHOL_F(copy_factor)(other.factor, c);
   voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
-} 
+}
 
 problem::~problem() {
   CHOL_F(free_dense)(&b, c);
@@ -112,13 +111,13 @@ problem::~problem() {
 }
 
 void problem::solve(std::vector<bool>& fuses) {
-  cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+  cholmod_dense* x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
 
-  if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+  if (((double*)x->x)[0] != ((double*)x->x)[0]) {
     throw nanex;
   }
 
-  cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+  cholmod_dense* y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
 
   double alpha[2] = {1, 0};
   double beta[2] = {0, 0};
@@ -130,7 +129,7 @@ void problem::solve(std::vector<bool>& fuses) {
     if (fuses[i]) {
       sol.currents[i] = 0;
     } else {
-      sol.currents[i] = ((double *)y->x)[i];
+      sol.currents[i] = ((double*)y->x)[i];
 
       graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
       graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
@@ -166,15 +165,15 @@ void problem::break_edge(unsigned e, bool unbreak) {
 
   unsigned n = factor->n;
 
-  cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
+  cholmod_sparse* update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
 
   unsigned s1, s2;
   s1 = v0 < v1 ? v0 : v1;
   s2 = v0 < v1 ? v1 : v0;
 
-  CHOL_INT *pp = (CHOL_INT *)update_mat->p;
-  CHOL_INT *ii = (CHOL_INT *)update_mat->i;
-  double *xx = (double *)update_mat->x;
+  CHOL_INT* pp = (CHOL_INT*)update_mat->p;
+  CHOL_INT* ii = (CHOL_INT*)update_mat->i;
+  double* xx = (double*)update_mat->x;
 
   for (unsigned i = 0; i <= s1; i++) {
     pp[i] = 0;
@@ -189,8 +188,8 @@ void problem::break_edge(unsigned e, bool unbreak) {
   xx[0] = 1;
   xx[1] = -1;
 
-  cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
-      update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c);
+  cholmod_sparse* perm_update_mat =
+      CHOL_F(submatrix)(update_mat, (CHOL_INT*)factor->Perm, factor->n, NULL, -1, true, true, c);
 
   CHOL_F(updown)(unbreak, perm_update_mat, factor, c);
 
@@ -208,8 +207,7 @@ void problem::break_edge(unsigned e, bool unbreak) {
       ind = r0.x < r1.x ? unbreak : !unbreak;
     }
 
-    ((double *)b->x)[G.edges[e].v[ind]]  -= 1.0;
-    ((double *)b->x)[G.edges[e].v[!ind]] += 1.0;
+    ((double*)b->x)[G.edges[e].v[ind]] -= 1.0;
+    ((double*)b->x)[G.edges[e].v[!ind]] += 1.0;
   }
 }
-