diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
index aa9ad61b88e0969acea52976327e7982c97a1548..de45148a4ee0ea9dc344b39b66d1ad10c88f0bee 100644
--- a/arbor/include/arbor/morph/morphology.hpp
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -67,19 +67,16 @@ mlocation_list minset(const morphology&, const mlocation_list&);
mlocation canonical(const morphology&, mlocation);
// Represent a (possibly empty or disconnected) region on a morphology.
-// Wraps an mcable_list, and satisfies the additional constraints:
//
-// I. Any two cables on the same branch are strictly disjoint, i.e.
-// for cables p and q on the same branch, either p.prox_pos > q.dist_pos
-// or p.dist_pos < q.prox_pos.
+// Wraps an mcable_list, and satisfies the additional constraint that
+// any two cables on the same branch are strictly disjoint, i.e.
+// for cables p and q on the same branch, either p.prox_pos > q.dist_pos
+// or p.dist_pos < q.prox_pos.
//
-// II. For any branch b on the morphology tree that intersects the subset
-// described by the extent, there exists a cable in the extent defined
-// on this branch b.
-//
-// While an mextent can only be built from an mcable_list in conjunction with
-// a morphology, union, intersection, and location membership operations can
-// be performed without one.
+// Union, intersection, and location membership operations can be performed
+// without a morphology.
+// A morphology is required to assert the invariant that an mextent does
+// not contain branches not in the morphology.
struct mextent {
mextent() = default;
mextent(const mextent&) = default;
@@ -88,10 +85,13 @@ struct mextent {
mextent& operator=(const mextent&) = default;
mextent& operator=(mextent&&) = default;
- mextent(const morphology&, const mcable_list&);
+ mextent(const mcable_list&);
- bool test_invariants() const; // check invariant (I) above.
- bool test_invariants(const morphology&) const; // check invariants (I) and (II) above.
+ // Check that the cable segments are valid, and that cables are strictly disjoint.
+ bool test_invariants() const;
+ // Checks the above, along with asserting that only cables in the morphology
+ // are present.
+ bool test_invariants(const morphology&) const;
const mcable_list& cables() const {
return cables_;
diff --git a/arbor/include/arbor/morph/region.hpp b/arbor/include/arbor/morph/region.hpp
index bce98b639088b8eadd32c03c419332d1e8938520..19fbcbaaa2ceef56aebb3cff151a62b413e78e98 100644
--- a/arbor/include/arbor/morph/region.hpp
+++ b/arbor/include/arbor/morph/region.hpp
@@ -158,6 +158,9 @@ region z_dist_from_root_ge(double r);
// Region with all segments in a cell.
region all();
+// The extent of a region, i.e. including all fork cover points.
+region super(region);
+
// Region associated with a name.
region named(std::string);
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index ad6a9a4e6a4c33ede2c0806e7b646acb11c706db..ca3fbdb3800a23aaead3001c4d5bcdb1333aacdb 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -171,7 +171,7 @@ mlocation_list thingify_(const on_branches_& ob, const mprovider& p) {
}
std::ostream& operator<<(std::ostream& o, const on_branches_& x) {
- return o << "(on_branchs " << x.pos << ")";
+ return o << "(on_branches " << x.pos << ")";
}
// Named locset.
@@ -235,10 +235,10 @@ mlocation_list thingify_(const most_distal_& n, const mprovider& p) {
}
std::ostream& operator<<(std::ostream& o, const most_distal_& x) {
- return o << "(locset \"" << x.reg << "\")";
+ return o << "(distal \"" << x.reg << "\")";
}
-// Most distal points of a region
+// Most proximal points of a region
struct most_proximal_: locset_tag {
explicit most_proximal_(region reg): reg(std::move(reg)) {}
@@ -251,15 +251,19 @@ locset most_proximal(region reg) {
mlocation_list thingify_(const most_proximal_& n, const mprovider& p) {
auto extent = thingify(n.reg, p);
- if (extent.empty()) return {};
-
arb_assert(extent.test_invariants(p.morphology()));
- auto most_prox = extent.cables().front();
- return {{most_prox.branch, most_prox.prox_pos}};
+
+ // 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});
+ }
+
+ return minset(p.morphology(), P);
}
std::ostream& operator<<(std::ostream& o, const most_proximal_& x) {
- return o << "(locset \"" << x.reg << "\")";
+ return o << "(proximal \"" << x.reg << "\")";
}
diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp
index ad93e472873d9bd3efe6fa6c9e60aab56a67151a..877334c6bf832164aeb1e32d739c17028ff86251 100644
--- a/arbor/morph/morphology.cpp
+++ b/arbor/morph/morphology.cpp
@@ -280,26 +280,15 @@ mlocation canonical(const morphology& m, mlocation loc) {
return loc;
}
-// Constructing an mextent from an mcable_list consists of taking the union
-// of any intersecting cables, and adding any required zero-length cables
-// around fork-points.
-
-mcable_list build_mextent_cables(const morphology& m, const mcable_list& cables) {
+// Merge overlapping cables so that none of the cables in the output overlap.
+// Used by the mextent constructor.
+mcable_list build_mextent_cables(const mcable_list& cables) {
arb_assert(arb::test_invariants(cables));
- std::unordered_set<msize_t> branch_tails;
-
mcable_list cs;
for (auto& c: cables) {
mcable* prev = cs.empty()? nullptr: &cs.back();
- if (c.prox_pos==0) {
- branch_tails.insert(m.branch_parent(c.branch));
- }
- if (c.dist_pos==1) {
- branch_tails.insert(c.branch);
- }
-
if (prev && prev->branch==c.branch && prev->dist_pos>=c.prox_pos) {
prev->dist_pos = std::max(prev->dist_pos, c.dist_pos);
}
@@ -308,42 +297,16 @@ mcable_list build_mextent_cables(const morphology& m, const mcable_list& cables)
}
}
- if (!branch_tails.empty()) {
- std::vector<mcable> fork_covers;
-
- for (auto b: branch_tails) {
- if (b!=mnpos) fork_covers.push_back(mcable{b, 1., 1.});
- for (auto b_child: m.branch_children(b)) {
- fork_covers.push_back(mcable{b_child, 0., 0.});
- }
- }
- util::sort(fork_covers);
-
- // Merge cables in cs with 0-length cables corresponding to fork covers.
- mcable_list a;
- a.swap(cs);
-
- for (auto c: util::merge_view(a, fork_covers)) {
- mcable* prev = cs.empty()? nullptr: &cs.back();
-
- if (prev && prev->branch==c.branch && prev->dist_pos>=c.prox_pos) {
- prev->dist_pos = std::max(prev->dist_pos, c.dist_pos);
- }
- else {
- cs.push_back(c);
- }
- }
- }
-
return cs;
-
}
-mextent::mextent(const morphology& m, const mcable_list& cables):
- cables_(build_mextent_cables(m, cables)) {}
+mextent::mextent(const mcable_list& cables):
+ cables_(build_mextent_cables(cables)) {}
bool mextent::test_invariants() const {
- // Checks for sortedness:
+ // Checks for:
+ // * validity of each cables.
+ // * sortedness of cables in list.
if (!arb::test_invariants(cables_)) return false;
// Check for intersections:
@@ -365,35 +328,6 @@ bool mextent::test_invariants(const morphology& m) const {
// Too many branches?
if (!empty() && cables_.back().branch>=m.num_branches()) return false;
- // Gather branches which are covered at the proximal or distal end:
- std::unordered_set<msize_t> branch_heads, branch_tails;
- for (auto& c: cables_) {
- if (c.prox_pos==0) branch_heads.insert(c.branch);
- if (c.dist_pos==1) branch_tails.insert(c.branch);
- }
-
- // There should be an entry in branch_tails for parent(j) for all j
- // in branch_heads, and an entry j in branch_heads for every j
- // with parent(j) in branch_tails.
-
- for (auto b: branch_heads) {
- msize_t parent = m.branch_parent(b);
- if (parent==mnpos) {
- branch_tails.insert(mnpos);
- }
- else if (!branch_tails.count(parent)) {
- return false;
- }
- }
-
- for (auto b: branch_tails) {
- for (auto child: m.branch_children(b)) {
- if (!branch_heads.count(child)) {
- return false;
- }
- }
- }
-
return true;
}
@@ -481,26 +415,5 @@ 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
diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
index 75663487806e8f0135b2899faaba0c9856d03986..2d057f0d2b1cc548f4c8218f0b9edf157b73242f 100644
--- a/arbor/morph/region.cpp
+++ b/arbor/morph/region.cpp
@@ -9,6 +9,7 @@
#include <arbor/morph/region.hpp>
#include <arbor/util/optional.hpp>
+#include "util/mergeview.hpp"
#include "util/span.hpp"
#include "util/strprintf.hpp"
#include "util/range.hpp"
@@ -65,7 +66,7 @@ mextent thingify_(const cable_& reg, const mprovider& p) {
if (reg.cable.branch>=p.morphology().num_branches()) {
throw no_such_branch(reg.cable.branch);
}
- return mextent(p.morphology(), mcable_list{{reg.cable}});
+ return mextent(mcable_list{{reg.cable}});
}
std::ostream& operator<<(std::ostream& o, const cable_& c) {
@@ -94,7 +95,7 @@ mextent thingify_(const cable_list_& reg, const mprovider& p) {
if (last_branch>=p.morphology().num_branches()) {
throw no_such_branch(last_branch);
}
- return mextent(p.morphology(), reg.cables);
+ return mextent(reg.cables);
}
std::ostream& operator<<(std::ostream& o, const cable_list_& x) {
@@ -182,7 +183,7 @@ mextent thingify_(const tagged_& reg, const mprovider& p) {
}
}
- return mextent(m, L);
+ return mextent(L);
}
std::ostream& operator<<(std::ostream& o, const tagged_& t) {
@@ -204,7 +205,7 @@ mextent thingify_(const all_&, const mprovider& p) {
for (auto i: util::make_span(nb)) {
branches.push_back({i,0,1});
}
- return mextent(p.morphology(), branches);
+ return mextent(branches);
}
std::ostream& operator<<(std::ostream& o, const all_& t) {
@@ -276,7 +277,7 @@ mextent thingify_(const distal_interval_& reg, const mprovider& p) {
}
util::sort(L);
- return mextent(m, L);
+ return mextent(L);
}
std::ostream& operator<<(std::ostream& o, const distal_interval_& d) {
@@ -330,7 +331,7 @@ mextent thingify_(const proximal_interval_& reg, const mprovider& p) {
}
util::sort(L);
- return mextent(m, L);
+ return mextent(L);
}
std::ostream& operator<<(std::ostream& o, const proximal_interval_& d) {
@@ -351,7 +352,7 @@ mextent radius_cmp(const mprovider& p, region r, double val, comp_op op) {
}
}
- return intersect(reg_extent, mextent(p.morphology(), cmp_cables));
+ return intersect(reg_extent, mextent(cmp_cables));
}
// Region with all segments with radius less than r
@@ -435,7 +436,7 @@ mextent projection_cmp(const mprovider& p, double v, comp_op op) {
for (auto i: util::make_span(m.num_branches())) {
util::append(L, e.projection_cmp(i, val, op));
}
- return mextent(p.morphology(), L);
+ return mextent(L);
}
// Region with all segments with projection less than val
@@ -551,6 +552,74 @@ std::ostream& operator<<(std::ostream& o, const named_& x) {
return o << "(region \"" << x.name << "\")";
}
+// Adds all cover points to a region.
+// Ensures that all valid representations of all fork points in the region are included.
+struct super_ {
+ region reg;
+};
+
+region super(region r) {
+ return region(super_{std::move(r)});
+}
+
+mextent thingify_(const super_& r, const mprovider& p) {
+ const auto& m = p.morphology();
+ auto cables = thingify(r.reg, p).cables();
+ std::unordered_set<msize_t> branch_tails;
+
+ mcable_list cs;
+ for (auto& c: cables) {
+ mcable* prev = cs.empty()? nullptr: &cs.back();
+
+ if (c.prox_pos==0) {
+ branch_tails.insert(m.branch_parent(c.branch));
+ }
+ if (c.dist_pos==1) {
+ branch_tails.insert(c.branch);
+ }
+
+ if (prev && prev->branch==c.branch && prev->dist_pos>=c.prox_pos) {
+ prev->dist_pos = std::max(prev->dist_pos, c.dist_pos);
+ }
+ else {
+ cs.push_back(c);
+ }
+ }
+
+ if (!branch_tails.empty()) {
+ std::vector<mcable> fork_covers;
+
+ for (auto b: branch_tails) {
+ if (b!=mnpos) fork_covers.push_back(mcable{b, 1., 1.});
+ for (auto b_child: m.branch_children(b)) {
+ fork_covers.push_back(mcable{b_child, 0., 0.});
+ }
+ }
+ util::sort(fork_covers);
+
+ // Merge cables in cs with 0-length cables corresponding to fork covers.
+ mcable_list a;
+ a.swap(cs);
+
+ for (auto c: util::merge_view(a, fork_covers)) {
+ mcable* prev = cs.empty()? nullptr: &cs.back();
+
+ if (prev && prev->branch==c.branch && prev->dist_pos>=c.prox_pos) {
+ prev->dist_pos = std::max(prev->dist_pos, c.dist_pos);
+ }
+ else {
+ cs.push_back(c);
+ }
+ }
+ }
+
+ return {cs};
+}
+
+std::ostream& operator<<(std::ostream& o, const super_& r) {
+ return o << "(super " << r.reg << ")";
+}
+
// Intersection of two regions.
diff --git a/python/cells.cpp b/python/cells.cpp
index 852b51a7138e118362c4c74953dbc3df10379230..a74856bd58d46cd8ec0b70f0d19abe1c1cfc5741 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -530,7 +530,7 @@ void register_cells(pybind11::module& m) {
[](arb::cable_cell& c, const char* label) {return c.concrete_locset(label);},
"label"_a, "The locations of the cell morphology for a locset label.")
.def("region",
- [](arb::cable_cell& c, const char* label) {return c.concrete_region(label);},
+ [](arb::cable_cell& c, const char* label) {return c.concrete_region(label).cables();},
"label"_a, "The cable segments of the cell morphology for a region label.")
// Discretization control.
.def("compartments_on_samples",
diff --git a/test/unit/test_cv_geom.cpp b/test/unit/test_cv_geom.cpp
index 63dca9639ff22219ac872fe958a3e787f12b19d7..f71bb9e3d7ed5e270ec1f73142df25ab304a850d 100644
--- a/test/unit/test_cv_geom.cpp
+++ b/test/unit/test_cv_geom.cpp
@@ -5,6 +5,7 @@
#include <arbor/cable_cell.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/locset.hpp>
+#include <arbor/morph/region.hpp>
#include "fvm_layout.hpp"
#include "util/rangeutil.hpp"
@@ -126,6 +127,10 @@ TEST(cv_geom, trivial) {
TEST(cv_geom, one_cv_per_branch) {
using namespace common_morphology;
+ auto super = [] (const arb::morphology& m, arb::mcable c) {
+ return thingify(arb::reg::super(arb::region(c)), arb::mprovider(m)).cables();
+ };
+
for (auto& p: test_morphologies) {
if (p.second.empty()) continue;
SCOPED_TRACE(p.first);
@@ -133,7 +138,8 @@ TEST(cv_geom, one_cv_per_branch) {
cable_cell cell{p.second};
auto& m = cell.morphology();
- cv_geometry geom = cv_geometry_from_ends(cell, sum(ls::on_branches(0), ls::on_branches(1)));
+ cv_geometry geom =
+ cv_geometry_from_ends(cell, sum(ls::on_branches(0), ls::on_branches(1)));
EXPECT_TRUE(verify_cv_children(geom));
// Expect trivial CVs at every fork point, and single-cable CVs for each branch.
@@ -154,7 +160,7 @@ TEST(cv_geom, one_cv_per_branch) {
EXPECT_TRUE(n_branch_child(c.branch)>1);
}
// Cables in trivial CV should be the same as those in the extent over the point.
- EXPECT_TRUE(testing::seq_eq(mextent{m, mcable_list{c}}.cables(), cables));
+ EXPECT_TRUE(testing::seq_eq(super(m,c), cables));
}
else {
ASSERT_EQ(1u, cables.size());
@@ -165,9 +171,9 @@ TEST(cv_geom, one_cv_per_branch) {
// Confirm parent CV is fork CV:
if (i>0) {
- mextent fork_ext{m, mcable_list{{c.branch, 0}}};
+ auto fork_ext = super(m, {c.branch, 0});
mcable_list pcables = util::assign_from(geom.cables(geom.cv_parent[i]));
- ASSERT_TRUE(testing::cablelist_eq(fork_ext.cables(), pcables));
+ ASSERT_TRUE(testing::cablelist_eq(fork_ext, pcables));
}
}
}
@@ -299,10 +305,14 @@ TEST(cv_geom, location_cv) {
cable_cell cell{m_reg_b6};
auto& m = cell.morphology();
- auto cv_extent = [&m](const cv_geometry& geom, auto cv) {
+ auto cv_extent = [](const cv_geometry& geom, auto cv) {
mcable_list cl;
util::assign(cl, geom.cables(cv));
- return mextent{m, cl};
+ return mextent(cl);
+ };
+
+ auto super = [] (const arb::morphology& m, arb::mcable c) {
+ return thingify(arb::reg::super(arb::region(c)), arb::mprovider(m)).cables();
};
// Two CVs per branch, plus trivial CV at forks.
@@ -323,9 +333,8 @@ TEST(cv_geom, location_cv) {
mcable cable0 = cables.front();
ASSERT_TRUE(cable0.prox_pos==cable0.dist_pos);
- mextent fork_ext{m, mcable_list{cable0}};
mcable_list clist = util::assign_from(cables);
- ASSERT_TRUE(testing::cablelist_eq(fork_ext.cables(), clist));
+ ASSERT_TRUE(testing::cablelist_eq(super(m, cable0), clist));
}
}
diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp
index ff5ee103e113dde1c54951e1018d374ea398f621..35b0db1c25d2cc14b36388d7047a9baca9f5d896 100644
--- a/test/unit/test_morph_expr.cpp
+++ b/test/unit/test_morph_expr.cpp
@@ -33,6 +33,7 @@ TEST(region, expr_repn) {
auto t1 = reg::tagged(1);
auto t2 = reg::tagged(2);
auto t3 = reg::tagged(3);
+ auto s = reg::super(t3);
auto all = reg::all();
EXPECT_EQ(to_string(c1), "(cable 1 0 1)");
@@ -41,11 +42,11 @@ TEST(region, expr_repn) {
EXPECT_EQ(to_string(b1), "(cable 1 0 1)");
EXPECT_EQ(to_string(t1), "(tag 1)");
EXPECT_EQ(to_string(t2), "(tag 2)");
+ EXPECT_EQ(to_string(s), "(super (tag 3))");
EXPECT_EQ(to_string(intersect(c1, t2)), "(intersect (cable 1 0 1) (tag 2))");
EXPECT_EQ(to_string(join(c1, t2)), "(join (cable 1 0 1) (tag 2))");
EXPECT_EQ(to_string(join(t1, t2, t3)), "(join (join (tag 1) (tag 2)) (tag 3))");
EXPECT_EQ(to_string(intersect(t1, t2, t3)), "(intersect (intersect (tag 1) (tag 2)) (tag 3))");
- EXPECT_EQ(to_string(intersect(join(c1, t2), c2)), "(intersect (join (cable 1 0 1) (tag 2)) (cable 4 0.125 0.5))");
EXPECT_EQ(to_string(all), "(all)");
}
@@ -397,7 +398,7 @@ TEST(region, thingify_simple_morphologies) {
{{ 0,100, 0, 2}, 2},
};
- // with a spherical root
+ // Build morphology with a spherical root
sample_tree sm(samples, parents);
mprovider mp(morphology(sm, true));
@@ -417,7 +418,7 @@ TEST(region, thingify_simple_morphologies) {
auto reg1_ = distal_interval(mid0_, 74);
auto reg2_ = proximal_interval(end1_, 45);
auto reg3_ = proximal_interval(end1_, 91);
- auto reg4_ = distal_interval(end1_, 0);
+ auto reg4_ = distal_interval(end1_, 100);
auto reg5_ = distal_interval(start1_, 0);
auto reg6_ = proximal_interval(start1_, 0);
@@ -431,8 +432,8 @@ TEST(region, thingify_simple_morphologies) {
EXPECT_TRUE(region_eq(mp, reg2_, mcable_list{{1,0.5,1}}));
EXPECT_TRUE(region_eq(mp, reg3_, mcable_list{{0, 0.75, 1}, {1,0,1}}));
EXPECT_TRUE(region_eq(mp, reg4_, mcable_list{{1,1,1}}));
- EXPECT_TRUE(region_eq(mp, reg5_, mcable_list{{0,1,1}}));
- EXPECT_TRUE(region_eq(mp, reg6_, mcable_list{{0,1,1}}));
+ EXPECT_TRUE(region_eq(mp, reg5_, mcable_list{{1,0,0}}));
+ EXPECT_TRUE(region_eq(mp, reg6_, mcable_list{{1,0,0}}));
}
}
@@ -485,6 +486,7 @@ TEST(region, thingify_moderate_morphologies) {
using reg::branch;
using reg::all;
using reg::cable;
+ using reg::super;
auto soma = tagged(1);
auto axon = tagged(2);
@@ -504,9 +506,6 @@ TEST(region, thingify_moderate_morphologies) {
mcable b3_{3,0,1};
cl all_ = {b0_,b1_,b2_,b3_};
- mcable c_end1_{1,1,1};
- mcable c_root_{0,0,0};
-
EXPECT_TRUE(region_eq(mp, all(), all_));
EXPECT_TRUE(region_eq(mp, axon, cl{b1_}));
EXPECT_TRUE(region_eq(mp, dend, cl{b0_,b3_}));
@@ -516,9 +515,13 @@ TEST(region, thingify_moderate_morphologies) {
// Test that intersection correctly generates zero-length cables at
// parent-child interfaces.
- EXPECT_TRUE(region_eq(mp, intersect(apic, dend), cl{c_end1_}));
- EXPECT_TRUE(region_eq(mp, intersect(apic, axon), cl{c_end1_}));
- EXPECT_TRUE(region_eq(mp, intersect(axon, dend), cl{c_root_, c_end1_}));
+ EXPECT_TRUE(region_eq(mp, intersect(apic, dend), cl{}));
+ EXPECT_TRUE(region_eq(mp, intersect(super(apic), super(dend)), cl{{1,1,1}, {2,0,0}, {3,0,0}}));
+ EXPECT_TRUE(region_eq(mp, intersect(apic, axon), cl{}));
+ EXPECT_TRUE(region_eq(mp, intersect(super(apic), axon), cl{{1,1,1}}));
+ EXPECT_TRUE(region_eq(mp, intersect(axon, dend), cl{}));
+ EXPECT_TRUE(region_eq(mp, intersect(super(axon), dend), cl{{0,0,0}, {3,0,0}}));
+ EXPECT_TRUE(region_eq(mp, intersect(super(dend), axon), cl{{1,0,0}, {1,1,1}}));
// Test distal and proximal interavls
auto start0_ = location(0, 0 );
@@ -537,39 +540,39 @@ TEST(region, thingify_moderate_morphologies) {
auto reg_d_ = join(cable(0,0,0.7), cable(2,0,0.5), cable(3,0.1,0.9));
// Distal from point and/or interval
- EXPECT_TRUE(region_eq(mp, distal_interval(start0_, 1000), mcable_list{{0,0,1}}));
- EXPECT_TRUE(region_eq(mp, distal_interval(quar_1_, 150), mcable_list{{1,0.25,1}, {2,0,0.75}, {3,0,0.375}}));
- EXPECT_TRUE(region_eq(mp, distal_interval(mid1_, 1000), mcable_list{{1,0.5,1}, {2,0,1}, {3,0,1}}));
- EXPECT_TRUE(region_eq(mp, distal_interval(mid1_, 150), mcable_list{{1,0.5,1}, {2,0,1}, {3,0,0.5}}));
- EXPECT_TRUE(region_eq(mp, distal_interval(end1_, 100), mcable_list{{2,0,1},{3,0,0.5}}));
- EXPECT_TRUE(region_eq(mp, distal_interval(join(quar_1_, mid1_), 150), mcable_list{{1,0.25,1}, {2,0,1}, {3,0,0.5}}));
- EXPECT_TRUE(region_eq(mp, distal_interval(join(quar_1_, loc_3_1_), 150), mcable_list{{1,0.25,1}, {2,0,0.75}, {3,0,0.375}, {3,0.65,1}}));
- EXPECT_TRUE(region_eq(mp, distal_interval(join(quar_1_, loc_3_1_), 150), mcable_list{{1,0.25,1}, {2,0,0.75}, {3,0,0.375}, {3,0.65,1}}));
+ EXPECT_TRUE(region_eq(mp, distal_interval(start0_, 1000), cl{{0,0,1}}));
+ EXPECT_TRUE(region_eq(mp, distal_interval(quar_1_, 150), cl{{1,0.25,1}, {2,0,0.75}, {3,0,0.375}}));
+ EXPECT_TRUE(region_eq(mp, distal_interval(mid1_, 1000), cl{{1,0.5,1}, {2,0,1}, {3,0,1}}));
+ EXPECT_TRUE(region_eq(mp, distal_interval(mid1_, 150), cl{{1,0.5,1}, {2,0,1}, {3,0,0.5}}));
+ EXPECT_TRUE(region_eq(mp, distal_interval(end1_, 100), cl{{1,1,1}, {2,0,1}, {3,0,0.5}}));
+ EXPECT_TRUE(region_eq(mp, distal_interval(join(quar_1_, mid1_), 150), cl{{1,0.25,1}, {2,0,1}, {3,0,0.5}}));
+ EXPECT_TRUE(region_eq(mp, distal_interval(join(quar_1_, loc_3_1_), 150), cl{{1,0.25,1}, {2,0,0.75}, {3,0,0.375}, {3,0.65,1}}));
+ EXPECT_TRUE(region_eq(mp, distal_interval(join(quar_1_, loc_3_1_), 150), cl{{1,0.25,1}, {2,0,0.75}, {3,0,0.375}, {3,0.65,1}}));
// Proximal from point and/or interval
- EXPECT_TRUE(region_eq(mp, proximal_interval(mid3_, 100), mcable_list{{3,0,0.5}}));
- EXPECT_TRUE(region_eq(mp, proximal_interval(mid3_, 150), mcable_list{{1,0.5,1}, {3,0,0.5}}));
- EXPECT_TRUE(region_eq(mp, proximal_interval(end2_, 150), mcable_list{{1,0.5,1}, {2,0,1}}));
- EXPECT_TRUE(region_eq(mp, proximal_interval(end2_, 500), mcable_list{{1,0,1}, {2,0,1}}));
- EXPECT_TRUE(region_eq(mp, proximal_interval(loc_3_0_, 100), mcable_list{{1,0.8,1}, {3,0,0.4}}));
- EXPECT_TRUE(region_eq(mp, proximal_interval(join(loc_3_0_, mid2_), 120), mcable_list{{1,0.3,1}, {2,0,0.5}, {3, 0, 0.4}}));
+ EXPECT_TRUE(region_eq(mp, proximal_interval(mid3_, 100), cl{{3,0,0.5}}));
+ EXPECT_TRUE(region_eq(mp, proximal_interval(mid3_, 150), cl{{1,0.5,1}, {3,0,0.5}}));
+ EXPECT_TRUE(region_eq(mp, proximal_interval(end2_, 150), cl{{1,0.5,1}, {2,0,1}}));
+ EXPECT_TRUE(region_eq(mp, proximal_interval(end2_, 500), cl{{1,0,1}, {2,0,1}}));
+ EXPECT_TRUE(region_eq(mp, proximal_interval(loc_3_0_, 100), cl{{1,0.8,1}, {3,0,0.4}}));
+ EXPECT_TRUE(region_eq(mp, proximal_interval(join(loc_3_0_, mid2_), 120), cl{{1,0.3,1}, {2,0,0.5}, {3, 0, 0.4}}));
// Test radius_lt and radius_gt
- EXPECT_TRUE(region_eq(mp, radius_lt(all(), 2), mcable_list{{0,0,0.55}, {1,0,0.325}, {3,0.375,0.75}}));
- EXPECT_TRUE(region_eq(mp, radius_lt(all(), 3), mcable_list{{0,0,1}, {1,0,0.55}, {2,6.0/9.0,1}, {3,0.25,1}}));
- EXPECT_TRUE(region_eq(mp, radius_gt(all(), 2), mcable_list{{0,0.55,1}, {1,0.325,1}, {2,0,1}, {3,0,0.375}, {3,0.75,1}}));
- EXPECT_TRUE(region_eq(mp, radius_gt(all(), 3), mcable_list{{1,0.55,1}, {2,0,6.0/9.0}, {3,0,0.25}}));
-
- EXPECT_TRUE(region_eq(mp, radius_le(all(), 2), mcable_list{{0,0,0.55}, {1,0,0.325}, {2,1,1}, {3,0.375,0.75}}));
- EXPECT_TRUE(region_eq(mp, radius_le(all(), 3), mcable_list{{0,0,1}, {1,0,0.55}, {2,6.0/9.0,1}, {3,0.25,1}}));
- EXPECT_TRUE(region_eq(mp, radius_ge(all(), 2), mcable_list{{0,0.55,1}, {1,0.325,1}, {2,0,1}, {3,0,0.375}, {3,0.75,1}}));
- EXPECT_TRUE(region_eq(mp, radius_ge(all(), 3), mcable_list{{1,0.55,1}, {2,0,6.0/9.0}, {3,0,0.25}}));
-
- EXPECT_TRUE(region_eq(mp, radius_lt(reg_a_, 2), mcable_list{{0,0.1,0.4},{3,0.375,0.4}}));
- EXPECT_TRUE(region_eq(mp, radius_gt(reg_a_, 2), mcable_list{{2,0,1},{3,0.1,0.375}}));
- EXPECT_TRUE(region_eq(mp, radius_lt(reg_b_, 2), mcable_list{{0,0.1,0.4}}));
- EXPECT_TRUE(region_eq(mp, radius_gt(reg_c_, 2), mcable_list{{0,0.55,0.7},{2,0,0.5},{3,0.1,0.375},{3,0.9,1}}));
- EXPECT_TRUE(region_eq(mp, radius_gt(reg_d_, 2), mcable_list{{0,0.55,0.7},{2,0,0.5},{3,0.1,0.375},{3,0.75,0.9}}));
+ EXPECT_TRUE(region_eq(mp, radius_lt(all(), 2), cl{{0,0,0.55}, {1,0,0.325}, {3,0.375,0.75}}));
+ EXPECT_TRUE(region_eq(mp, radius_lt(all(), 3), cl{{0,0,1}, {1,0,0.55}, {2,6.0/9.0,1}, {3,0.25,1}}));
+ EXPECT_TRUE(region_eq(mp, radius_gt(all(), 2), cl{{0,0.55,1}, {1,0.325,1}, {2,0,1}, {3,0,0.375}, {3,0.75,1}}));
+ EXPECT_TRUE(region_eq(mp, radius_gt(all(), 3), cl{{1,0.55,1}, {2,0,6.0/9.0}, {3,0,0.25}}));
+
+ EXPECT_TRUE(region_eq(mp, radius_le(all(), 2), cl{{0,0,0.55}, {1,0,0.325}, {2,1,1}, {3,0.375,0.75}}));
+ EXPECT_TRUE(region_eq(mp, radius_le(all(), 3), cl{{0,0,1}, {1,0,0.55}, {2,6.0/9.0,1}, {3,0.25,1}}));
+ EXPECT_TRUE(region_eq(mp, radius_ge(all(), 2), cl{{0,0.55,1}, {1,0.325,1}, {2,0,1}, {3,0,0.375}, {3,0.75,1}}));
+ EXPECT_TRUE(region_eq(mp, radius_ge(all(), 3), cl{{1,0.55,1}, {2,0,6.0/9.0}, {3,0,0.25}}));
+
+ EXPECT_TRUE(region_eq(mp, radius_lt(reg_a_, 2), cl{{0,0.1,0.4},{3,0.375,0.4}}));
+ EXPECT_TRUE(region_eq(mp, radius_gt(reg_a_, 2), cl{{2,0,1},{3,0.1,0.375}}));
+ EXPECT_TRUE(region_eq(mp, radius_lt(reg_b_, 2), cl{{0,0.1,0.4}}));
+ EXPECT_TRUE(region_eq(mp, radius_gt(reg_c_, 2), cl{{0,0.55,0.7},{2,0,0.5},{3,0.1,0.375},{3,0.9,1}}));
+ EXPECT_TRUE(region_eq(mp, radius_gt(reg_d_, 2), cl{{0,0.55,0.7},{2,0,0.5},{3,0.1,0.375},{3,0.75,0.9}}));
// Test some more interesting intersections and unions.
@@ -654,11 +657,11 @@ TEST(region, thingify_moderate_morphologies) {
// 123456789 123456789
// |xxxxx | | lhs
// | |xxxxx | rhs
- // x | | rand
+ // | | | rand
// |xxxxx |xxxxx | ror
lhs = cable(0,0,.5);
rhs = cable(1,0,.5);
- rand = cl{{0,0,0}};
+ rand = cl{};
ror = cl{{0,0,.5},{1,0,.5}};
EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand));
EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror));
@@ -672,11 +675,11 @@ TEST(region, thingify_moderate_morphologies) {
// 123456789 123456789
// |xxxxx | | lhs
// | |xxxxx | rhs
- // x | | rand
+ // | | | rand
// |xxxxx |xxxxx | ror
lhs = cable(2,0,.5);
rhs = cable(3,0,.5);
- rand = cl{{1,1,1}};
+ rand = cl{};
ror = cl{{2,0,.5},{3,0,.5}};
EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand));
EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror));
@@ -706,6 +709,7 @@ TEST(region, thingify_moderate_morphologies) {
}
}
+
TEST(region, thingify_complex_morphologies) {
using pvec = std::vector<msize_t>;
using svec = std::vector<msample>;
@@ -730,6 +734,7 @@ TEST(region, thingify_complex_morphologies) {
{
mprovider mp(m);
using reg::cable;
+ using reg::super;
using ls::most_distal;
using ls::most_proximal;
@@ -739,6 +744,7 @@ TEST(region, thingify_complex_morphologies) {
auto reg_d_ = join(cable(2,0,0.9), cable(3,0.1,0.1), cable(4,0.1,0.6));
auto reg_e_ = join(cable(2,0,0.9), cable(4,0.1,0.1), cable(5,0.1,0.6));
auto reg_f_ = join(cable(7,0,1), cable(2,0,0.9), cable(4,0.1,0.1), cable(5,0.1,0.6));
+ auto reg_g_ = join(cable(5,0,1), cable(6,0,1));
EXPECT_TRUE(mlocationlist_eq(thingify(most_distal(reg_a_), mp), mlocation_list{{0,0.9},{1,0.4}}));
EXPECT_TRUE(mlocationlist_eq(thingify(most_distal(reg_b_), mp), mlocation_list{{0,0.9},{1,0.5}}));
@@ -747,12 +753,15 @@ TEST(region, thingify_complex_morphologies) {
EXPECT_TRUE(mlocationlist_eq(thingify(most_distal(reg_e_), mp), mlocation_list{{5,0.6}}));
EXPECT_TRUE(mlocationlist_eq(thingify(most_distal(reg_f_), mp), mlocation_list{{5,0.6},{7,1}}));
- EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_a_), mp), mlocation_list{{0,0}}));
- EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_b_), mp), mlocation_list{{0,0}}));
- EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_c_), mp), mlocation_list{{0,0}}));
- EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_d_), mp), mlocation_list{{1,1}}));
- EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_e_), mp), mlocation_list{{1,1}}));
- EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_f_), mp), mlocation_list{{1,1}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_a_), mp), mlocation_list{{0,0}, {1,0.1}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_b_), mp), mlocation_list{{0,0}, {1,0.1}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_c_), mp), mlocation_list{{0,0}, {1,0.1}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_d_), mp), mlocation_list{{2,0}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_e_), mp), mlocation_list{{2,0}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_f_), mp), mlocation_list{{2,0}, {7,0}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(reg_g_), mp), mlocation_list{{5,0}, {6,0}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(super(reg_f_)), mp), mlocation_list{{1,1}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(super(reg_g_)), mp), mlocation_list{{4,1}}));
}
}
{
@@ -800,6 +809,7 @@ TEST(region, thingify_complex_morphologies) {
{3,0,0.179407353580315756}}));
EXPECT_TRUE(region_eq(mp, z_dist_from_root_ge(20),
mcable_list{{0,1,1},
+ {1,0,0},
{1,0.578250901781922829,1},
{2,0,0.61499300915417734997},
{2,0.8349970039232188642,1},
diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp
index fd3bece9d4e26dbf4701f53bbe51cea0660116fd..b4423c7606f6c95b9798a7596194fea44a35c67a 100644
--- a/test/unit/test_morphology.cpp
+++ b/test/unit/test_morphology.cpp
@@ -664,19 +664,31 @@ TEST(mextent, invariants) {
{{ 0,300, 0, 2}, 3},
};
- morphology m(sample_tree(samples, parents));
+ // Instantiate morphology with spherical_root=false.
+ morphology m(sample_tree(samples, parents), false);
- mextent x1(m, cl{{1, 0.25, 0.75}});
+ mextent x1(cl{{1, 0.25, 0.75}});
ASSERT_TRUE(x1.test_invariants(m));
EXPECT_TRUE(cablelist_eq(cl{{1, 0.25, 0.75}}, x1.cables()));
- mextent x2(m, cl{{1, 0., 0.75}});
+ mextent x2(cl{{1, 0., 0.75}});
ASSERT_TRUE(x2.test_invariants(m));
- EXPECT_TRUE(cablelist_eq(cl{{0, 0, 0}, {1, 0., 0.75}}, x2.cables()));
+ EXPECT_TRUE(cablelist_eq(cl{{1, 0., 0.75}}, x2.cables()));
- mextent x3(m, cl{{2, 0., 1.}});
+ mextent x3(cl{{2, 0., 1.}});
ASSERT_TRUE(x3.test_invariants(m));
- EXPECT_TRUE(cablelist_eq(cl{{1, 1, 1}, {2, 0., 1.}, {3, 0., 0.}}, x3.cables()));
+ EXPECT_TRUE(cablelist_eq(cl{{2, 0., 1.}}, x3.cables()));
+
+ // Test that overlapping cables are merged on construction.
+ mextent x4(cl{{0, 0.0, 1.0}, {0, 0.5, 0.7},
+ {1, 0.2, 0.8}, {1, 0.4, 0.6},
+ {2, 0.2, 0.5}, {2, 0.5, 0.6},
+ {3, 0.2, 0.5}, {3, 0.6, 0.7}});
+ ASSERT_TRUE(x4.test_invariants(m));
+ EXPECT_TRUE(cablelist_eq(cl{{0, 0.0, 1.0},
+ {1, 0.2, 0.8},
+ {2, 0.2, 0.6},
+ {3, 0.2, 0.5}, {3, 0.6, 0.7}}, x4.cables()));
}
TEST(mextent, intersects) {
@@ -709,7 +721,7 @@ TEST(mextent, intersects) {
morphology m(sample_tree(samples, parents));
- mextent x1(m, cl{{1, 0.25, 0.75}});
+ mextent x1(cl{{1, 0.25, 0.75}});
EXPECT_TRUE(x1.intersects(mlocation{1, 0.25}));
EXPECT_TRUE(x1.intersects(mlocation{1, 0.5}));
EXPECT_TRUE(x1.intersects(mlocation{1, 0.75}));
@@ -723,67 +735,14 @@ TEST(mextent, intersects) {
EXPECT_TRUE(x1.intersects(mcable{1, 0.75, 1.0}));
EXPECT_FALSE(x1.intersects(mcable{1, 0.8, 1.0}));
- mextent x2(m, cl{{1, 0., 0.75}});
+ mextent x2(cl{{1, 0., 0.75}});
EXPECT_TRUE(x2.intersects(mlocation{1, 0.}));
- EXPECT_TRUE(x2.intersects(mlocation{0, 0.}));
+ EXPECT_FALSE(x2.intersects(mlocation{0, 0.}));
EXPECT_FALSE(x2.intersects(mlocation{0, 1.}));
- mextent x3(m, cl{{2, 0., 1.}});
+ mextent x3(cl{{2, 0., 1.}});
EXPECT_FALSE(x3.intersects(mlocation{1, 0.}));
- EXPECT_TRUE(x3.intersects(mlocation{1, 1.}));
- EXPECT_TRUE(x3.intersects(mlocation{3, 0.}));
+ EXPECT_FALSE(x3.intersects(mlocation{1, 1.}));
+ EXPECT_FALSE(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)));
-}