diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
index de45148a4ee0ea9dc344b39b66d1ad10c88f0bee..7c00efed00577a27f46ac025d416e62a43adb965 100644
--- a/arbor/include/arbor/morph/morphology.hpp
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -62,8 +62,15 @@ public:
};
// Morphology utility functions.
+
+// Find the set of locations in an mlocation_list for which there
+// are no other locations that are more proximal in that list.
mlocation_list minset(const morphology&, const mlocation_list&);
+// Find the set of locations in an mlocation_list for which there
+// are no other locations that are more distal in the list.
+mlocation_list maxset(const morphology&, const mlocation_list&);
+
mlocation canonical(const morphology&, mlocation);
// Represent a (possibly empty or disconnected) region on a morphology.
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index ca3fbdb3800a23aaead3001c4d5bcdb1333aacdb..dc9b7617a00ff65134abc2d1f86933c01933ad54 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -205,33 +205,12 @@ locset most_distal(region reg) {
}
mlocation_list thingify_(const most_distal_& n, const mprovider& p) {
+ // Make a list of the distal ends of each cable segment.
mlocation_list L;
-
- auto cables = thingify(n.reg, p).cables();
- util::sort(cables, [](const auto& l, const auto& r){return (l.branch < r.branch) && (l.dist_pos < r.dist_pos);});
-
- std::unordered_set<msize_t> branches_visited;
- for (auto it = cables.rbegin(); it!=cables.rend(); ++it) {
- auto bid = it->branch;
- auto pos = it->dist_pos;
-
- // Exclude any initial branch points unless they are top-level.
- if (pos==0 && p.morphology().branch_parent(bid)!=mnpos) {
- continue;
- }
-
- // Check if any other points on the branch or any of its children has been added as a distal point
- if (branches_visited.count(bid)) continue;
- L.push_back(canonical(p.morphology(), mlocation{bid, pos}));
- while (bid != mnpos) {
- branches_visited.insert(bid);
- bid = p.morphology().branch_parent(bid);
- }
+ for (auto& c: thingify(n.reg, p)) {
+ L.push_back({c.branch, c.dist_pos});
}
-
- util::sort(L);
- util::unique_in_place(L);
- return L;
+ return maxset(p.morphology(), L);
}
std::ostream& operator<<(std::ostream& o, const most_distal_& x) {
@@ -250,16 +229,13 @@ locset most_proximal(region reg) {
}
mlocation_list thingify_(const most_proximal_& n, const mprovider& p) {
- auto extent = thingify(n.reg, p);
- arb_assert(extent.test_invariants(p.morphology()));
-
// Make a list of the proximal ends of each cable segment.
- mlocation_list P;
- for (const auto& c: extent.cables()) {
- P.push_back({c.branch, c.prox_pos});
+ mlocation_list L;
+ for (auto& c: thingify(n.reg, p)) {
+ L.push_back({c.branch, c.prox_pos});
}
- return minset(p.morphology(), P);
+ return minset(p.morphology(), L);
}
std::ostream& operator<<(std::ostream& o, const most_proximal_& x) {
diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp
index 877334c6bf832164aeb1e32d739c17028ff86251..4f131e6be09633d07e38309bcc06ff26dd7919e7 100644
--- a/arbor/morph/morphology.cpp
+++ b/arbor/morph/morphology.cpp
@@ -1,5 +1,6 @@
#include <iostream>
#include <stack>
+#include <unordered_map>
#include <utility>
#include <arbor/morph/morphexcept.hpp>
@@ -272,6 +273,37 @@ mlocation_list minset(const morphology& m, const mlocation_list& in) {
return L;
}
+mlocation_list maxset(const morphology& m, const mlocation_list& in_) {
+ mlocation_list L;
+
+ // Sort the input in reverse order, so that more distal locations
+ // come first.
+ mlocation_list in = in_;
+ util::sort(in, [](const auto& l, const auto& r) {return r<l;});
+
+ // List of branches that have had a more distal location found.
+ std::unordered_set<msize_t> br;
+ for (auto loc: in) {
+ auto b = loc.branch;
+
+ // A child of this branch has already been visited: a more distal
+ // location has already been found, so we can skip.
+ if (br.count(b)) continue;
+
+ // Add the location to the maxset.
+ L.push_back(loc);
+
+ // Mark the branch and its parents.
+ while (b!=mnpos) {
+ br.insert(b);
+ b = m.branch_parent(b);
+ }
+ }
+
+ std::reverse(L.begin(), L.end());
+ return L;
+}
+
mlocation canonical(const morphology& m, mlocation loc) {
if (loc.pos==0) {
msize_t parent = m.branch_parent(loc.branch);
diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
index 2d057f0d2b1cc548f4c8218f0b9edf157b73242f..3893ab11871faa946d5f45a77be9e8ec92600fc1 100644
--- a/arbor/morph/region.cpp
+++ b/arbor/morph/region.cpp
@@ -335,7 +335,7 @@ mextent thingify_(const proximal_interval_& reg, const mprovider& p) {
}
std::ostream& operator<<(std::ostream& o, const proximal_interval_& d) {
- return o << "(distal_interval " << d.end << " " << d.distance << ")";
+ return o << "(proximal_interval " << d.end << " " << d.distance << ")";
}
mextent radius_cmp(const mprovider& p, region r, double val, comp_op op) {
diff --git a/python/morph_parse.cpp b/python/morph_parse.cpp
index 64288f0f35cde5f1d2b47af21719cbfe6c584f1a..38111045a0ee2cb11a8ff17e5653cd95e3927525 100644
--- a/python/morph_parse.cpp
+++ b/python/morph_parse.cpp
@@ -1,6 +1,7 @@
#include <arbor/util/any.hpp>
#include <arbor/morph/region.hpp>
#include <arbor/morph/locset.hpp>
+#include <limits>
#include "error.hpp"
#include "s_expr.hpp"
@@ -187,8 +188,16 @@ std::unordered_multimap<std::string, evaluator> eval_map {
"'region' with 1 argument: (name:string)")},
{"distal_interval", make_call<arb::locset, double>(arb::reg::distal_interval,
"'distal_interval' with 2 arguments: (start:locset extent:real)")},
+ {"distal_interval", make_call<arb::locset>(
+ [](arb::locset ls){return arb::reg::distal_interval(std::move(ls), std::numeric_limits<double>::max());},
+ "'distal_interval' with 1 argument: (start:locset)")},
{"proximal_interval",make_call<arb::locset, double>(arb::reg::proximal_interval,
"'proximal_interval' with 2 arguments: (start:locset extent:real)")},
+ {"proximal_interval", make_call<arb::locset>(
+ [](arb::locset ls){return arb::reg::proximal_interval(std::move(ls), std::numeric_limits<double>::max());},
+ "'proximal_interval' with 1 argument: (start:locset)")},
+ {"super", make_call<arb::region>(arb::reg::super,
+ "'super' with 1 argment: (reg:region)")},
{"radius_lt",make_call<arb::region, double>(arb::reg::radius_lt,
"'radius_lt' with 2 arguments: (reg:region radius:real)")},
{"radius_le",make_call<arb::region, double>(arb::reg::radius_le,
diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp
index 8ab96c184bac8a6195aa63e673c2a0517d1e8972..48d987488bac96473b15ec5fb95f09c3b47d7fa1 100644
--- a/test/unit/test_morphology.cpp
+++ b/test/unit/test_morphology.cpp
@@ -636,6 +636,57 @@ TEST(morphology, minset) {
}
}
+TEST(morphology, maxset) {
+ 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::morphology m(sm, false);
+
+ EXPECT_EQ((ll{}), maxset(m, ll{}));
+ EXPECT_EQ((ll{{2,0.1}}), maxset(m, ll{{2,0.1}}));
+ EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), maxset(m, ll{{0,0.5}, {1,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}}), maxset(m, ll{{0,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}, {2,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {1,0.5}, {2,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}, {2,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {2,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}, {2,0.5}, {3,1}}), maxset(m, ll{{0,0}, {0,0.5}, {2,0.5}, {3,0}, {3,1}}));
+ EXPECT_EQ((ll{{2,0.5}, {3,0.5}}), maxset(m, ll{{1,0.5}, {2,0.5}, {3,0}, {3,0.2}, {3,0.5}}));
+ }
+ {
+ arb::morphology m(sm, true);
+
+ EXPECT_EQ((ll{}), maxset(m, ll{}));
+ EXPECT_EQ((ll{{2,0.1}}), maxset(m, ll{{2,0.1}}));
+ EXPECT_EQ((ll{{1,0.5}}), maxset(m, ll{{0,0.5}, {1,0.5}}));
+ EXPECT_EQ((ll{{1,0}, {2,0}}), maxset(m, ll{{2,0}, {1,0}, {0,0}}));
+ EXPECT_EQ((ll{{1,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
+ EXPECT_EQ((ll{{1,1}, {3,0.1}, {4,0.7}}), maxset(m, ll{{1,0.5}, {1,1}, {3,0.1}, {4,0.5}, {4,0.7}}));
+ }
+}
+
// Testing mextent; intersection and join operations are
// exercised by region/locset thingifies in test_morph_expr.cpp.