From 21c400b2fa970eb357211af27089b48858b5183a Mon Sep 17 00:00:00 2001 From: Ben Cumming <bcumming@cscs.ch> Date: Thu, 9 Apr 2020 10:10:32 +0200 Subject: [PATCH] Extents don't include the cover of fork points by default. (#998) This PR removes the requirement that all cover points are included in a region. The motivation is to allow more flexible definion of regions, particularly the proximal and distal sets thereof. The other motivation is that the author finds it much simpler to reason about, however others find the existing approach easier to reason about. The changes: * `mextent` does not always include points on the cover at fork points. * it still enforces that its constituent cables have no intersections and are ordered. * a `super` region expression adds the cover points to a region. * update `most_proximal` to return the `minset` of the proximal points in a region's cables. * fix some cut and paste errors in comments and printing of locset expressions. --- arbor/include/arbor/morph/morphology.hpp | 28 +++--- arbor/include/arbor/morph/region.hpp | 3 + arbor/morph/locset.cpp | 20 +++-- arbor/morph/morphology.cpp | 103 ++------------------- arbor/morph/region.cpp | 85 ++++++++++++++++-- python/cells.cpp | 2 +- test/unit/test_cv_geom.cpp | 25 ++++-- test/unit/test_morph_expr.cpp | 110 ++++++++++++----------- test/unit/test_morphology.cpp | 89 +++++------------- 9 files changed, 216 insertions(+), 249 deletions(-) diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp index aa9ad61b..de45148a 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 bce98b63..19fbcbaa 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 ad6a9a4e..ca3fbdb3 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 ad93e472..877334c6 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 75663487..2d057f0d 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 852b51a7..a74856bd 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 63dca963..f71bb9e3 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 ff5ee103..35b0db1c 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 fd3bece9..b4423c76 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))); -} -- GitLab