Add canonical repn for mextent cables. (#989)

Provide canonical function for mextent that returns a reduced
mcable_list representation which excludes redundant zero-length
cables.

arbor/include/arbor/morph/morphology.hpp

@@ -128,5 +128,8 @@ private:
     mcable_list cables_;
 };
 
+// Reduced representation of an extent, excluding zero-length cables
+// that are covered by more proximal or non-zero-length cables.
+mcable_list canonical(const morphology& m, const mextent& a);
 
 } // namespace arb

arbor/morph/morphology.cpp

@@ -481,5 +481,26 @@ mextent join(const mextent& a, const mextent& b) {
     return m;
 }
 
+mcable_list canonical(const morphology& m, const mextent& a) {
+    // For zero-length cables representing isolated points, keep
+    // only the most proximal. All other zero-length cables should be
+    // elided.
+
+    mcable_list result;
+    std::unordered_set<msize_t> remove_set;
+
+    for (auto& c: a.cables()) {
+        if (c.prox_pos==0 && c.dist_pos>0) {
+            remove_set.insert(m.branch_parent(c.branch));
+        }
+    }
+    for (auto& c: a.cables()) {
+        if (c.prox_pos==1 && remove_set.count(c.branch)) continue;
+        if (c.dist_pos==0 && m.branch_parent(c.branch)!=mnpos) continue;
+        result.push_back(c);
+    }
+    return result;
+}
+
 } // namespace arb
 

test/unit/test_morphology.cpp

@@ -734,3 +734,56 @@ TEST(mextent, intersects) {
     EXPECT_TRUE(x3.intersects(mlocation{3, 0.}));
     EXPECT_FALSE(x3.intersects(mlocation{3, 1.}));
 }
+
+TEST(mextent, canonical) {
+    using namespace arb;
+    using testing::cablelist_eq;
+
+    using pvec = std::vector<msize_t>;
+    using svec = std::vector<msample>;
+    using cl = mcable_list;
+
+    // Eight samples
+    //                  no-sphere
+    //          sample   branch
+    //            0         0
+    //           1 3       0 1
+    //          2   4     0   1
+    //             5 6       2 3
+    //                7         3
+    pvec parents = {mnpos, 0, 1, 0, 3, 4, 4, 6};
+    svec samples = {
+        {{  0,  0,  0,  2}, 3},
+        {{ 10,  0,  0,  2}, 3},
+        {{100,  0,  0,  2}, 3},
+        {{  0, 10,  0,  2}, 3},
+        {{  0,100,  0,  2}, 3},
+        {{100,100,  0,  2}, 3},
+        {{  0,200,  0,  2}, 3},
+        {{  0,300,  0,  2}, 3},
+    };
+
+    morphology m(sample_tree(samples, parents));
+
+    mextent x1(m, cl{{0, 0.3, 0.7}, {1, 0.5, 0.7}, {3, 0.4, 1}});
+    EXPECT_TRUE(cablelist_eq(canonical(m, x1), x1.cables()));
+
+    mcable_list cl2{{1, 0.5, 1}, {2, 0, 0.5}}; // a reduced representation
+    mextent x2(m, cl2);
+    EXPECT_FALSE(cablelist_eq(cl2, x2.cables()));
+    EXPECT_TRUE(cablelist_eq(cl2, canonical(m, x2)));
+
+    mcable_list cl3{{3, 0, 0}}; // extent is cover of fork point
+    mextent x3(m, cl3);
+    EXPECT_TRUE(cablelist_eq(cl{{1, 1, 1}}, canonical(m, x3)));
+
+    mcable_list cl4{{1, 1, 1}, {3, 0, 0.5}}; // canonical repn excludes zero-length cable
+    mextent x4(m, cl4);
+    EXPECT_TRUE(cablelist_eq(cl{{3, 0, 0.5}}, canonical(m, x4)));
+
+    // Heads of top-level branches should be preserved.
+    mcable_list cl5{{1, 0, 0}};
+    mextent x5(m, cl5);
+    EXPECT_TRUE(cablelist_eq(cl{{0, 0, 0}, {1, 0, 0}}, x5.cables()));
+    EXPECT_TRUE(cablelist_eq(x5.cables(), canonical(m, x5)));
+}