Add minset operation to embedded morphology (#864)

* Add `minset` method for location lists that computes the minimal set by the partial order imposed by the morphology tree.

arbor/morph/em_morphology.cpp

 #include <mutex>
+#include <stack>
 #include <vector>
 
 #include <arbor/morph/error.hpp>
@@ -141,4 +142,40 @@ mlocation em_morphology::canonicalize(mlocation loc) const {
     return loc;
 }
 
+mlocation_list em_morphology::minset(const mlocation_list& in) const {
+    mlocation_list L;
+
+    std::stack<msize_t> stack;
+
+    // All root branches must be searched.
+    for (auto c: morph_.branch_children(mnpos)) {
+        stack.push(c);
+    }
+
+    // Depth-first traversal of the branch tree.
+    while (!stack.empty()) {
+        auto branch = stack.top();
+        stack.pop();
+
+        // Search for a location on branch.
+        auto it = std::lower_bound(in.begin(), in.end(), mlocation{branch, 0});
+
+        // If found, insert to the minset and skip the rest of this sub-tree.
+        if (it!=in.end() && it->branch==branch) {
+            L.push_back(*it);
+            continue;
+        }
+
+        // No location on this branch, so continue searching in this sub-tree.
+        for (auto c: morph_.branch_children(branch)) {
+            stack.push(c);
+        }
+    }
+
+    util::sort(L);
+
+    return L;
+}
+
 } // namespace arb
+

arbor/morph/em_morphology.hpp

@@ -34,6 +34,8 @@ public:
     // representation of loc.
     // If include_loc is false, the input location is excluded from the result.
     mlocation_list cover(mlocation, bool include_loc=true) const;
+
+    mlocation_list minset(const mlocation_list&) const;
 };
 
 } // namespace arb

test/unit/test_em_morphology.cpp

@@ -187,6 +187,57 @@ TEST(em_morphology, cover) {
     }
 }
 
+TEST(em_morphology, minset) {
+    using pvec = std::vector<arb::msize_t>;
+    using svec = std::vector<arb::msample>;
+    using ll = arb::mlocation_list;
+    constexpr auto npos = arb::mnpos;
+
+    // Eight samples
+    //                  no-sphere  sphere
+    //          sample   branch    branch
+    //            0         0         0
+    //           1 3       0 1       1 2
+    //          2   4     0   1     1   2
+    //             5 6       2 3       3 4
+    //                7         3         4
+    pvec parents = {npos, 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},
+    };
+    arb::sample_tree sm(samples, parents);
+    {
+        arb::em_morphology em(arb::morphology(sm, false));
+
+        EXPECT_EQ((ll{}), em.minset(ll{}));
+        EXPECT_EQ((ll{{2,0.1}}), em.minset(ll{{2,0.1}}));
+        EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), em.minset(ll{{0,0.5}, {1,0.5}}));
+        EXPECT_EQ((ll{{0,0.5}}), em.minset(ll{{0,0.5}}));
+        EXPECT_EQ((ll{{0,0}, {1,0}}), em.minset(ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
+        EXPECT_EQ((ll{{0,0}, {1,0.5}}), em.minset(ll{{0,0}, {0,0.5}, {1,0.5}, {2,0.5}}));
+        EXPECT_EQ((ll{{0,0}, {2,0.5}}), em.minset(ll{{0,0}, {0,0.5}, {2,0.5}}));
+        EXPECT_EQ((ll{{0,0}, {2,0.5}, {3,0}}), em.minset(ll{{0,0}, {0,0.5}, {2,0.5}, {3,0}, {3,1}}));
+        EXPECT_EQ((ll{{0,0}, {1,0}}), em.minset(ll{{0,0}, {0,0.5}, {1,0}, {2,0.5}, {3,0}, {3,1}}));
+    }
+    {
+        arb::em_morphology em(arb::morphology(sm, true));
+
+        EXPECT_EQ((ll{}), em.minset(ll{}));
+        EXPECT_EQ((ll{{2,0.1}}), em.minset(ll{{2,0.1}}));
+        EXPECT_EQ((ll{{0,0.5}}), em.minset(ll{{0,0.5}, {1,0.5}}));
+        EXPECT_EQ((ll{{0,0.5}}), em.minset(ll{{0,0.5}}));
+        EXPECT_EQ((ll{{0,0}}), em.minset(ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
+        EXPECT_EQ((ll{{1,0.5}, {3,0.1}, {4,0.5}}), em.minset(ll{{1,0.5}, {1,1}, {3,0.1}, {4,0.5}, {4,0.7}}));
+    }
+}
+
 // Test expressions that describe location sets (locset).
 TEST(locset, expressions) {
     auto root = arb::ls::root();
@@ -625,3 +676,4 @@ TEST(region, thingify) {
         EXPECT_EQ(thingify(join(lhs, rhs), em), ror);
     }
 }
+