diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp index 2ec43f27ac6d750d7e531de9282d7ba34b725a40..2d643732fd54e0f1c6bcf72a723627dcea4f4504 100644 --- a/arbor/cable_cell.cpp +++ b/arbor/cable_cell.cpp @@ -94,10 +94,13 @@ struct cable_cell_impl { template <typename Property> void paint(const region& reg, const Property& prop) { - mcable_list cables = thingify(reg, provider); + mextent cables = thingify(reg, provider); auto& mm = get_region_map(prop); for (auto c: cables) { + // Skip zero-length cables in extent: + if (c.prox_pos==c.dist_pos) continue; + if (!mm.insert(c, prop)) { throw cable_cell_error(util::pprintf("cable {} overpaints", c)); } @@ -108,7 +111,7 @@ struct cable_cell_impl { return thingify(l, provider); } - mcable_list concrete_region(const region& r) const { + mextent concrete_region(const region& r) const { return thingify(r, provider); } }; @@ -145,7 +148,7 @@ mlocation_list cable_cell::concrete_locset(const locset& l) const { return impl_->concrete_locset(l); } -mcable_list cable_cell::concrete_region(const region& r) const { +mextent cable_cell::concrete_region(const region& r) const { return impl_->concrete_region(r); } diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp index 1a7522920efc7590cbe774738dbba523206e475b..f8b62eb1e545c89222efa2f2ec1f131889588666 100644 --- a/arbor/include/arbor/cable_cell.hpp +++ b/arbor/include/arbor/cable_cell.hpp @@ -182,7 +182,7 @@ public: mlocation_list concrete_locset(const locset&) const; // Access to a concrete list of cable segments for a region. - mcable_list concrete_region(const region&) const; + mextent concrete_region(const region&) const; // Generic access to painted and placed items. const cable_cell_region_map& region_assignments() const; diff --git a/arbor/include/arbor/morph/morphexcept.hpp b/arbor/include/arbor/morph/morphexcept.hpp index 2c97947daeb420b9cbd01b76221778160cda8e47..9620f6dc02fb1a71973ec4173037a7a5cc0a32d4 100644 --- a/arbor/include/arbor/morph/morphexcept.hpp +++ b/arbor/include/arbor/morph/morphexcept.hpp @@ -26,6 +26,10 @@ struct invalid_mcable: morphology_error { mcable cable; }; +struct invalid_mcable_list: morphology_error { + invalid_mcable_list(); +}; + struct invalid_sample_parent: morphology_error { invalid_sample_parent(msize_t parent, msize_t tree_size); msize_t parent; diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp index bd36d65fb78e7ba3ee718b30b8e1a44b68294840..4285e3aec034b67f22f6876b8fc535a0d6dc4627 100644 --- a/arbor/include/arbor/morph/morphology.hpp +++ b/arbor/include/arbor/morph/morphology.hpp @@ -36,9 +36,11 @@ public: msize_t num_samples() const; // The parent branch of branch b. + // Return mnpos if branch has no parent. msize_t branch_parent(msize_t b) const; // The child branches of branch b. + // If b is mnpos, return root branches. const std::vector<msize_t>& branch_children(msize_t b) const; // Branches with no children. @@ -64,4 +66,67 @@ 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. +// +// 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. +struct mextent { + mextent() = default; + mextent(const mextent&) = default; + mextent(mextent&&) = default; + + mextent& operator=(const mextent&) = default; + mextent& operator=(mextent&&) = default; + + mextent(const morphology&, const mcable_list&); + + bool test_invariants() const; // check invariant (I) above. + bool test_invariants(const morphology&) const; // check invariants (I) and (II) above. + + const mcable_list& cables() const { + return cables_; + } + + bool operator==(const mextent& a) const { return cables_==a.cables_; } + bool operator!=(const mextent& a) const { return cables_!=a.cables_; } + + bool intersects(const mcable_list& a) const; + bool intersects(const mcable& a) const { return intersects(mcable_list{a}); } + + bool intersects(const mextent& a) const { + return intersects(a.cables()); + } + bool intersects(mlocation loc) const { + return intersects(mcable{loc.branch, loc.pos, loc.pos}); + } + + friend mextent intersect(const mextent& a, const mextent& b); + friend mextent join(const mextent& a, const mextent& b); + + // Forward const container operations: + decltype(auto) cbegin() const { return cables_.cbegin(); } + decltype(auto) begin() const { return cables_.begin(); } + decltype(auto) cend() const { return cables_.cend(); } + decltype(auto) end() const { return cables_.end(); } + + bool empty() const { return cables_.empty(); } + std::size_t size() const { return cables_.size(); } + const mcable& front() const { return cables_.front(); } + const mcable& back() const { return cables_.back(); } + +private: + mcable_list cables_; +}; + + } // namespace arb diff --git a/arbor/include/arbor/morph/mprovider.hpp b/arbor/include/arbor/morph/mprovider.hpp index 727550ccc957cce8700ef872e85bcbff2a8fb5be..0169df17d7cc3cb74800d75b4c0227ed45c0c379 100644 --- a/arbor/include/arbor/morph/mprovider.hpp +++ b/arbor/include/arbor/morph/mprovider.hpp @@ -17,7 +17,7 @@ struct mprovider { explicit mprovider(arb::morphology m): mprovider(m, nullptr) {} // Throw exception on missing or recursive definition. - const mcable_list& region(const std::string& name) const; + const mextent& region(const std::string& name) const; const mlocation_list& locset(const std::string& name) const; // Read-only access to morphology and constructed embedding. @@ -34,7 +34,7 @@ private: struct circular_def {}; // Maps are mutated only during initialization phase of mprovider. - mutable std::unordered_map<std::string, util::either<mcable_list, circular_def>> regions_; + mutable std::unordered_map<std::string, util::either<mextent, circular_def>> regions_; mutable std::unordered_map<std::string, util::either<mlocation_list, circular_def>> locsets_; // Non-null only during initialization phase. diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp index fa6cb3eab733136757acfc7a715a0b090a746b0a..8848bdd714b0cfdc6c852a1fe4d288c1529e6cf4 100644 --- a/arbor/include/arbor/morph/primitives.hpp +++ b/arbor/include/arbor/morph/primitives.hpp @@ -83,6 +83,11 @@ mlocation_list join(const mlocation_list&, const mlocation_list&); mlocation_list intersection(const mlocation_list&, const mlocation_list&); // Describe an unbranched cable in the morphology. +// +// Cables are a representation of a closed interval of a branch in a morphology. +// They may be zero-length, and fork points in the morphology may have multiple, +// equivalent zero-length cable representations. + struct mcable { // The id of the branch on which the cable lies. msize_t branch; diff --git a/arbor/include/arbor/morph/region.hpp b/arbor/include/arbor/morph/region.hpp index 40db2f63b5a2f3e469e60d0774fa6786d40d714b..bce98b639088b8eadd32c03c419332d1e8938520 100644 --- a/arbor/include/arbor/morph/region.hpp +++ b/arbor/include/arbor/morph/region.hpp @@ -53,15 +53,16 @@ public: return *this; } - // Implicit conversion from mcable or mcable_list. + // Implicit conversion from mcable, mcable_list, or mextent. region(mcable); - region(const mcable_list&); + region(mextent); + region(mcable_list); // Implicitly convert string to named region expression. region(std::string label); region(const char* label); - friend mcable_list thingify(const region& r, const mprovider& m) { + friend mextent thingify(const region& r, const mprovider& m) { return r.impl_->thingify(m); } @@ -90,7 +91,7 @@ private: virtual ~interface() {} virtual std::unique_ptr<interface> clone() = 0; virtual std::ostream& print(std::ostream&) = 0; - virtual mcable_list thingify(const mprovider&) = 0; + virtual mextent thingify(const mprovider&) = 0; }; std::unique_ptr<interface> impl_; @@ -104,7 +105,7 @@ private: return std::unique_ptr<interface>(new wrap<Impl>(wrapped)); } - virtual mcable_list thingify(const mprovider& m) override { + virtual mextent thingify(const mprovider& m) override { return thingify_(wrapped, m); } diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp index 5accd8818f597fc8f8eb004c1c9104960ca6f1ba..54ad51e090af6acedec4b91e6ddb9a1afd378d15 100644 --- a/arbor/morph/locset.cpp +++ b/arbor/morph/locset.cpp @@ -203,20 +203,40 @@ locset most_distal(region reg) { return locset(most_distal_{std::move(reg)}); } +template <typename X> +void unique_in_place(std::vector<X>& v) { + if (v.empty()) return; + + auto write = v.begin(); + auto read = write; + + while (++read!=v.end()) { + if (*read==*write) continue; + if (++write!=read) *write = std::move(*read); + } + + v.erase(++write, v.end()); +} + mlocation_list thingify_(const most_distal_& n, const mprovider& p) { mlocation_list L; - auto cables = thingify(n.reg, p); + 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; + 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({bid, pos}); + L.push_back(canonical(p.morphology(), mlocation{bid, pos})); while (bid != mnpos) { branches_visited.insert(bid); bid = p.morphology().branch_parent(bid); @@ -224,6 +244,7 @@ mlocation_list thingify_(const most_distal_& n, const mprovider& p) { } util::sort(L); + unique_in_place(L); return L; } @@ -243,10 +264,11 @@ locset most_proximal(region reg) { } mlocation_list thingify_(const most_proximal_& n, const mprovider& p) { - auto cables = thingify(n.reg, p); - arb_assert(test_invariants(cables)); + auto extent = thingify(n.reg, p); + if (extent.empty()) return {}; - auto most_prox = cables.front(); + arb_assert(extent.test_invariants(p.morphology())); + auto most_prox = extent.cables().front(); return {{most_prox.branch, most_prox.prox_pos}}; } @@ -274,7 +296,8 @@ mlocation_list thingify_(const uniform_& u, const mprovider& p) { auto embed = p.embedding(); // Thingify the region and store relevant data - auto reg_cables = thingify(u.reg, p); + mextent reg_extent = thingify(u.reg, p); + const mcable_list& reg_cables = reg_extent.cables(); std::vector<double> lengths_bounds; auto lengths_part = util::make_partition(lengths_bounds, diff --git a/arbor/morph/morphexcept.cpp b/arbor/morph/morphexcept.cpp index 21463f20c185bacee26c211109f16397946dc00e..3364127573b08269f343f197d3dc08fad05c8f0c 100644 --- a/arbor/morph/morphexcept.cpp +++ b/arbor/morph/morphexcept.cpp @@ -25,6 +25,10 @@ invalid_mcable::invalid_mcable(mcable cable): cable(cable) {} +invalid_mcable_list::invalid_mcable_list(): + morphology_error("bad mcable_list") +{} + invalid_sample_parent::invalid_sample_parent(msize_t parent, msize_t tree_size): morphology_error(pprintf("invalid sample parent {} for a sample tree of size {}", parent, tree_size)), parent(parent), diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp index 3b59fef868ee5c0999ba03186427ab52dda21d59..ead854d28155bde26c0422cc34d4e392a0aaedd8 100644 --- a/arbor/morph/morphology.cpp +++ b/arbor/morph/morphology.cpp @@ -8,13 +8,14 @@ #include <arbor/morph/primitives.hpp> #include "morph/mbranch.hpp" +#include "util/mergeview.hpp" #include "util/rangeutil.hpp" #include "util/span.hpp" using arb::util::make_span; namespace arb { -namespace impl{ +namespace impl { std::vector<mbranch> branches_from_parent_index(const std::vector<msize_t>& parents, const std::vector<point_prop>& props, @@ -279,5 +280,206 @@ 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) { + 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); + } + 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; + +} + +mextent::mextent(const morphology& m, const mcable_list& cables): + cables_(build_mextent_cables(m, cables)) {} + +bool mextent::test_invariants() const { + // Checks for sortedness: + if (!arb::test_invariants(cables_)) return false; + + // Check for intersections: + for (auto i: util::count_along(cables_)) { + if (!i) continue; + + const auto& c = cables_[i]; + const auto& p = cables_[i-1]; + if (p.branch==c.branch && p.dist_pos>=c.prox_pos) return false; + } + + return true; +} + +bool mextent::test_invariants(const morphology& m) const { + // Check for sortedness, intersections: + if (!test_invariants()) return false; + + // 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; +} + +bool mextent::intersects(const mcable_list& a) const { + arb_assert(arb::test_invariants(a)); + + // Early exit? + if (empty() || a.empty() || + cables_.front().branch>a.back().branch || + cables_.back().branch<a.front().branch) + { + return false; + } + + auto from = cables_.begin(); + for (auto& c: a) { + auto i = std::lower_bound(from, cables_.end(), c); + + if (i!=cables_.end() && i->branch==c.branch) { + arb_assert(i->prox_pos>=c.prox_pos); + if (i->prox_pos<=c.dist_pos) return true; + } + + if (i!=cables_.begin()) { + auto j = std::prev(i); + if (j->branch==c.branch) { + arb_assert(j->prox_pos<c.prox_pos); + if (j->dist_pos>=c.prox_pos) return true; + } + } + + from = i; + } + + return false; +} + +mextent intersect(const mextent& a, const mextent& b) { + auto precedes = [](mcable x, mcable y) { + return x.branch<y.branch || (x.branch==y.branch && x.dist_pos<y.prox_pos); + }; + + mextent m; + auto ai = a.cables().begin(); + auto ae = a.cables().end(); + auto bi = b.cables().begin(); + auto be = b.cables().end(); + + while (ai!=ae && bi!=be) { + if (precedes(*ai, *bi)) { + ++ai; + } + else if (precedes(*bi, *ai)) { + ++bi; + } + else { + m.cables_.push_back(mcable{ai->branch, + std::max(ai->prox_pos, bi->prox_pos), + std::min(ai->dist_pos, bi->dist_pos)}); + if (ai->dist_pos<bi->dist_pos) { + ++ai; + } + else { + ++bi; + } + } + } + return m; +} + +mextent join(const mextent& a, const mextent& b) { + mextent m; + mcable_list& cs = m.cables_; + + for (auto c: util::merge_view(a.cables(), b.cables())) { + 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 m; +} + } // namespace arb diff --git a/arbor/morph/mprovider.cpp b/arbor/morph/mprovider.cpp index 8b8cb216c91e519e9007b05d62728412581c71c9..38e0226c0f2b6a57235290240432c926def2c3f1 100644 --- a/arbor/morph/mprovider.cpp +++ b/arbor/morph/mprovider.cpp @@ -61,7 +61,7 @@ static const auto& try_lookup(const mprovider& provider, const std::string& name } } -const mcable_list& mprovider::region(const std::string& name) const { +const mextent& mprovider::region(const std::string& name) const { const auto* regions_ptr = label_dict_ptr? &(label_dict_ptr->regions()): nullptr; return try_lookup(*this, name, regions_, regions_ptr, circular_def{}); } diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp index c98313a31e3a0e44e4462174b250b9eaa10d62a1..75663487806e8f0135b2899faaba0c9856d03986 100644 --- a/arbor/morph/region.cpp +++ b/arbor/morph/region.cpp @@ -1,4 +1,3 @@ -#include <set> #include <string> #include <vector> #include <stack> @@ -8,6 +7,7 @@ #include <arbor/morph/morphexcept.hpp> #include <arbor/morph/mprovider.hpp> #include <arbor/morph/region.hpp> +#include <arbor/util/optional.hpp> #include "util/span.hpp" #include "util/strprintf.hpp" @@ -17,160 +17,12 @@ namespace arb { namespace reg { -// Head and tail of an mcable as mlocations. -inline mlocation head(mcable c) { - return mlocation{c.branch, c.prox_pos}; -} - -inline mlocation tail(mcable c) { - return mlocation{c.branch, c.dist_pos}; -} - -// Returns true iff cable sections a and b: -// 1. are on the same branch -// 2. overlap, i.e. their union is not empty. -bool is_disjoint(const mcable& a, const mcable& b) { - if (a.branch!=b.branch) return true; - return a<b? a.dist_pos<b.prox_pos: b.dist_pos<a.prox_pos; -} - -// Take the union of two cable sections that are not disjoint. -mcable make_union(const mcable& a, const mcable& b) { - assert(!is_disjoint(a,b)); - return mcable{a.branch, std::min(a.prox_pos, b.prox_pos), std::max(a.dist_pos, b.dist_pos)}; -} - -// Take the intersection of two cable sections that are not disjoint. -mcable make_intersection(const mcable& a, const mcable& b) { - assert(!is_disjoint(a,b)); - - return mcable{a.branch, std::max(a.prox_pos, b.prox_pos), std::min(a.dist_pos, b.dist_pos)}; -} - -mcable_list merge(const mcable_list& v) { - if (v.size()<2) return v; - std::vector<mcable> L; - L.reserve(v.size()); - auto c = v.front(); - auto it = v.begin()+1; - while (it!=v.end()) { - if (!is_disjoint(c, *it)) { - c = make_union(c, *it); - } - else { - L.push_back(c); - c = *it; - } - ++it; - } - L.push_back(c); - return L; -} - -// List of colocated mlocations, excluding the parameter. -mlocation_list colocated(mlocation loc, const morphology& m) { - mlocation_list L{}; - - // Note: if the location is not at the end of a branch, there are no - // other colocated points. - - if (loc.pos==0) { - // Include head of each branch with same parent, - // and end of parent branch if not mnpos. - - auto p = m.branch_parent(loc.branch); - if (p!=mnpos) L.push_back({p, 1}); - - for (auto b: m.branch_children(p)) { - if (b!=loc.branch) L.push_back({b, 0}); - } - } - else if (loc.pos==1) { - // Include head of each child branch. - - for (auto b: m.branch_children(loc.branch)) { - L.push_back({b, 0}); - } - } - - return L; -} - -// Insert a zero-length region at the start of each child branch for every cable -// that includes the end of a branch. -// Insert a zero-length region at the end of the parent branch for each cable -// that includes the start of the branch. -mcable_list cover(mcable_list cables, const morphology& m) { - mcable_list L = cables; - - for (auto& c: cables) { - for (auto& x: colocated(head(c), m)) { - L.push_back({x.branch, x.pos, x.pos}); - } - for (auto& x: colocated(tail(c), m)) { - L.push_back({x.branch, x.pos, x.pos}); - } - } - - util::sort(L); - return L; -} - -mcable_list remove_cover(mcable_list cables, const morphology& m) { - // Find all zero-length cables at the end of cables, and convert to - // their canonical representation. - for (auto& c: cables) { - if (c.dist_pos==0 || c.prox_pos==1) { - auto cloc = canonical(m, head(c)); - c = {cloc.branch, cloc.pos, cloc.pos}; - } - } - // Some new zero-length cables may be out of order, so sort - // the cables. - util::sort(cables); - - // Remove multiple copies of zero-length cables if present. - return merge(cables); -} - -// Remove zero-length regions from a sorted cable_list that are in the cover of another region in the list. -mcable_list remove_covered_points(mcable_list cables, const morphology& m) { - struct branch_index_pair { - msize_t bid; - unsigned lid; - }; - std::vector<branch_index_pair> end_branches; - std::vector<unsigned> erase_indices; - - for (unsigned i = 0; i < cables.size(); i++) { - auto c = cables[i]; - // Save zero length cables at the distal end of a branch - // Or at the proximal end of the soma - if ((c.prox_pos==1 && c.dist_pos==1) || - (c.prox_pos==0 && c.dist_pos==0)) { - end_branches.push_back({c.branch, i}); - } - // Look for branches that are children of the cables saved in end_branches - else if (c.prox_pos==0) { - auto parent = m.branch_parent(c.branch); - if (parent==mnpos) parent = 0; - - auto it = std::find_if(end_branches.begin(), end_branches.end(), - [&](branch_index_pair& p) { return p.bid==parent;}); - - if (it!=end_branches.end()) { - erase_indices.push_back((*it).lid); - end_branches.erase(it); - } - } - } - - util::sort(erase_indices); - for (auto it = erase_indices.rbegin(); it != erase_indices.rend(); it++) { - cables.erase(cables.begin() + *it); - } +util::optional<mcable> intersect(const mcable& a, const mcable& b) { + if (a.branch!=b.branch) return util::nullopt; - return cables; + double prox = std::max(a.prox_pos, b.prox_pos); + double dist = std::min(a.dist_pos, b.dist_pos); + return prox<=dist? util::just(mcable{a.branch, prox, dist}): util::nullopt; } // Empty region. @@ -181,8 +33,8 @@ region nil() { return region{nil_{}}; } -mcable_list thingify_(const nil_& x, const mprovider&) { - return {}; +mextent thingify_(const nil_& x, const mprovider&) { + return mextent{}; } std::ostream& operator<<(std::ostream& o, const nil_& x) { @@ -209,17 +61,71 @@ region branch(msize_t bid) { return cable(bid, 0, 1); } -mcable_list thingify_(const cable_& reg, const mprovider& p) { +mextent thingify_(const cable_& reg, const mprovider& p) { if (reg.cable.branch>=p.morphology().num_branches()) { throw no_such_branch(reg.cable.branch); } - return {reg.cable}; + return mextent(p.morphology(), mcable_list{{reg.cable}}); } std::ostream& operator<<(std::ostream& o, const cable_& c) { return o << c.cable; } +// Exlicit list of cables. +// (Not part of front-end API: used by region ctor.) + +struct cable_list_: region_tag { + explicit cable_list_(mcable_list cs): cables(std::move(cs)) {} + mcable_list cables; +}; + +region cable_list(mcable_list cs) { + if (!test_invariants(cs)) { throw invalid_mcable_list(); } + return region(cable_list_{std::move(cs)}); +} + +mextent thingify_(const cable_list_& reg, const mprovider& p) { + if (reg.cables.empty()) { + return mextent{}; + } + + auto last_branch = reg.cables.back().branch; + if (last_branch>=p.morphology().num_branches()) { + throw no_such_branch(last_branch); + } + return mextent(p.morphology(), reg.cables); +} + +std::ostream& operator<<(std::ostream& o, const cable_list_& x) { + o << "(cable_list"; + for (auto c: x.cables) { o << ' ' << c; } + return o << ')'; +} + +// Explicit extent. +// (Not part of front-end API: used by region ctor.) + +struct extent_: region_tag { + explicit extent_(mextent x): extent(std::move(x)) {} + mextent extent; +}; + +region extent(mextent x) { + arb_assert(x.test_invariants()); + return region(extent_{std::move(x)}); +} + +mextent thingify_(const extent_& x, const mprovider& p) { + arb_assert(x.extent.test_invariants(p.morphology())); + return x.extent; +} + +std::ostream& operator<<(std::ostream& o, const extent_& x) { + o << "(extent"; + for (auto c: x.extent.cables()) { o << ' ' << c; } + return o << ')'; +} // Region with all segments with the same numeric tag. @@ -232,7 +138,7 @@ region tagged(int id) { return region(tagged_{id}); } -mcable_list thingify_(const tagged_& reg, const mprovider& p) { +mextent thingify_(const tagged_& reg, const mprovider& p) { const auto& m = p.morphology(); const auto& e = p.embedding(); size_t nb = m.num_branches(); @@ -275,10 +181,8 @@ mcable_list thingify_(const tagged_& reg, const mprovider& p) { start = std::find_if(last, end, matches); } } - if (L.size()<L.capacity()/4) { - L.shrink_to_fit(); - } - return L; + + return mextent(m, L); } std::ostream& operator<<(std::ostream& o, const tagged_& t) { @@ -293,21 +197,21 @@ region all() { return region(all_{}); } -mcable_list thingify_(const all_&, const mprovider& p) { +mextent thingify_(const all_&, const mprovider& p) { auto nb = p.morphology().num_branches(); mcable_list branches; branches.reserve(nb); for (auto i: util::make_span(nb)) { branches.push_back({i,0,1}); } - return branches; + return mextent(p.morphology(), branches); } std::ostream& operator<<(std::ostream& o, const all_& t) { return o << "(all)"; } -// Region with all segments distal from another region +// Region comprising points up to `distance` distal to a point in `start`. struct distal_interval_ { locset start; @@ -318,7 +222,7 @@ region distal_interval(locset start, double distance) { return region(distal_interval_{start, distance}); } -mcable_list thingify_(const distal_interval_& reg, const mprovider& p) { +mextent thingify_(const distal_interval_& reg, const mprovider& p) { const auto& m = p.morphology(); const auto& e = p.embedding(); @@ -370,14 +274,16 @@ mcable_list thingify_(const distal_interval_& reg, const mprovider& p) { first_branch = false; } } - return remove_covered_points(remove_cover(L, m), m); + + util::sort(L); + return mextent(m, L); } std::ostream& operator<<(std::ostream& o, const distal_interval_& d) { return o << "(distal_interval " << d.start << " " << d.distance << ")"; } -// Region with all segments proximal from another region +// Region comprising points up to `distance` proximal to a point in `end`. struct proximal_interval_ { locset end; @@ -388,7 +294,7 @@ region proximal_interval(locset end, double distance) { return region(proximal_interval_{end, distance}); } -mcable_list thingify_(const proximal_interval_& reg, const mprovider& p) { +mextent thingify_(const proximal_interval_& reg, const mprovider& p) { const auto& m = p.morphology(); const auto& e = p.embedding(); @@ -422,26 +328,30 @@ mcable_list thingify_(const proximal_interval_& reg, const mprovider& p) { L.push_back({branch, prox_pos, dist_pos}); } } - return remove_cover(L, m); + + util::sort(L); + return mextent(m, L); } std::ostream& operator<<(std::ostream& o, const proximal_interval_& d) { return o << "(distal_interval " << d.end << " " << d.distance << ")"; } -mcable_list radius_cmp(const mprovider& p, region r, double v, comp_op op) { +mextent radius_cmp(const mprovider& p, region r, double val, comp_op op) { const auto& e = p.embedding(); + auto reg_extent = thingify(r, p); - std::vector<mcable> L; - auto reg = thingify(r, p); - auto val = v; - for (auto c: reg) { - for (auto r: e.radius_cmp(c.branch, val, op)) { - if (is_disjoint(c, r)) continue; - L.push_back(make_intersection(c, r)); + msize_t bid = mnpos; + mcable_list cmp_cables; + + for (auto c: reg_extent) { + if (bid != c.branch) { + bid = c.branch; + util::append(cmp_cables, e.radius_cmp(bid, val, op)); } } - return remove_cover(L, p.morphology()); + + return intersect(reg_extent, mextent(p.morphology(), cmp_cables)); } // Region with all segments with radius less than r @@ -454,7 +364,7 @@ region radius_lt(region reg, double val) { return region(radius_lt_{reg, val}); } -mcable_list thingify_(const radius_lt_& r, const mprovider& p) { +mextent thingify_(const radius_lt_& r, const mprovider& p) { return radius_cmp(p, r.reg, r.val, comp_op::lt); } @@ -472,7 +382,7 @@ region radius_le(region reg, double val) { return region(radius_le_{reg, val}); } -mcable_list thingify_(const radius_le_& r, const mprovider& p) { +mextent thingify_(const radius_le_& r, const mprovider& p) { return radius_cmp(p, r.reg, r.val, comp_op::le); } @@ -490,7 +400,7 @@ region radius_gt(region reg, double val) { return region(radius_gt_{reg, val}); } -mcable_list thingify_(const radius_gt_& r, const mprovider& p) { +mextent thingify_(const radius_gt_& r, const mprovider& p) { return radius_cmp(p, r.reg, r.val, comp_op::gt); } @@ -508,7 +418,7 @@ region radius_ge(region reg, double val) { return region(radius_gt_{reg, val}); } -mcable_list thingify_(const radius_ge_& r, const mprovider& p) { +mextent thingify_(const radius_ge_& r, const mprovider& p) { return radius_cmp(p, r.reg, r.val, comp_op::ge); } @@ -516,7 +426,7 @@ std::ostream& operator<<(std::ostream& o, const radius_ge_& r) { return o << "(radius_ge " << r.reg << " " << r.val << ")"; } -mcable_list projection_cmp(const mprovider& p, double v, comp_op op) { +mextent projection_cmp(const mprovider& p, double v, comp_op op) { const auto& m = p.morphology(); const auto& e = p.embedding(); @@ -525,7 +435,7 @@ mcable_list 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 remove_cover(L, p.morphology()); + return mextent(p.morphology(), L); } // Region with all segments with projection less than val @@ -537,7 +447,7 @@ region projection_lt(double val) { return region(projection_lt_{val}); } -mcable_list thingify_(const projection_lt_& r, const mprovider& p) { +mextent thingify_(const projection_lt_& r, const mprovider& p) { return projection_cmp(p, r.val, comp_op::lt); } @@ -554,7 +464,7 @@ region projection_le(double val) { return region(projection_le_{val}); } -mcable_list thingify_(const projection_le_& r, const mprovider& p) { +mextent thingify_(const projection_le_& r, const mprovider& p) { return projection_cmp(p, r.val, comp_op::le); } @@ -571,7 +481,7 @@ region projection_gt(double val) { return region(projection_gt_{val}); } -mcable_list thingify_(const projection_gt_& r, const mprovider& p) { +mextent thingify_(const projection_gt_& r, const mprovider& p) { return projection_cmp(p, r.val, comp_op::gt); } @@ -588,7 +498,7 @@ region projection_ge(double val) { return region(projection_ge_{val}); } -mcable_list thingify_(const projection_ge_& r, const mprovider& p) { +mextent thingify_(const projection_ge_& r, const mprovider& p) { return projection_cmp(p, r.val, comp_op::ge); } @@ -633,7 +543,7 @@ region named(std::string name) { return region(named_{std::move(name)}); } -mcable_list thingify_(const named_& n, const mprovider& p) { +mextent thingify_(const named_& n, const mprovider& p) { return p.region(n.name); } @@ -650,38 +560,8 @@ struct reg_and: region_tag { reg_and(region lhs, region rhs): lhs(std::move(lhs)), rhs(std::move(rhs)) {} }; -mcable_list thingify_(const reg_and& P, const mprovider& p) { - auto& m = p.morphology(); - - using cable_it = std::vector<mcable>::const_iterator; - using cable_it_pair = std::pair<cable_it, cable_it>; - - auto lhs = cover(thingify(P.lhs, p), m); - auto rhs = cover(thingify(P.rhs, p), m); - - // Perform intersection - cable_it_pair it{lhs.begin(), rhs.begin()}; - cable_it_pair end{lhs.end(), rhs.end()}; - std::vector<mcable> L; - - bool at_end = it.first==end.first || it.second==end.second; - while (!at_end) { - bool first_less = *(it.first) < *(it.second); - auto& lhs = first_less? it.first: it.second; - auto& rhs = first_less? it.second: it.first; - if (!is_disjoint(*lhs, *rhs)) { - L.push_back(make_intersection(*lhs, *rhs)); - } - if (dist_loc(*lhs) < dist_loc(*rhs)) { - ++lhs; - } - else { - ++rhs; - } - at_end = it.first==end.first || it.second==end.second; - } - - return remove_covered_points(remove_cover(L, m), m); +mextent thingify_(const reg_and& P, const mprovider& p) { + return intersect(thingify(P.lhs, p), thingify(P.rhs, p)); } std::ostream& operator<<(std::ostream& o, const reg_and& x) { @@ -697,15 +577,8 @@ struct reg_or: region_tag { reg_or(region lhs, region rhs): lhs(std::move(lhs)), rhs(std::move(rhs)) {} }; -mcable_list thingify_(const reg_or& P, const mprovider& p) { - auto lhs = thingify(P.lhs, p); - auto rhs = thingify(P.rhs, p); - mcable_list L; - L.resize(lhs.size() + rhs.size()); - - std::merge(lhs.begin(), lhs.end(), rhs.begin(), rhs.end(), L.begin()); - - return merge(L); +mextent thingify_(const reg_or& P, const mprovider& p) { + return join(thingify(P.lhs, p), thingify(P.rhs, p)); } std::ostream& operator<<(std::ostream& o, const reg_or& x) { @@ -744,9 +617,12 @@ region::region(mcable c) { *this = reg::cable(c.branch, c.prox_pos, c.dist_pos); } -region::region(const mcable_list& cl) { - *this = std::accumulate(cl.begin(), cl.end(), reg::nil(), - [](auto& rg, auto& p) { return join(rg, region(p)); }); +region::region(mcable_list cl) { + *this = reg::cable_list(std::move(cl)); +} + +region::region(mextent x) { + *this = reg::extent(std::move(x)); } } // namespace arb diff --git a/arbor/util/iterutil.hpp b/arbor/util/iterutil.hpp index 90a376ee77b6150cfa75d3a106635dc9a22f593f..1556a50f53832acd55870d0e856c422502b7b8e6 100644 --- a/arbor/util/iterutil.hpp +++ b/arbor/util/iterutil.hpp @@ -89,10 +89,11 @@ decltype(auto) back(Seq& seq) { * present rvalues on dereference. */ template <typename V> -struct pointer_proxy: public V { - pointer_proxy(const V& v): V(v) {} - pointer_proxy(V&& v): V(std::move(v)) {} - const V* operator->() const { return this; } +struct pointer_proxy { + V v; + pointer_proxy(const V& v): v(v) {} + pointer_proxy(V&& v): v(std::move(v)) {} + const V* operator->() const { return &v; } }; /* diff --git a/arbor/util/mergeview.hpp b/arbor/util/mergeview.hpp new file mode 100644 index 0000000000000000000000000000000000000000..86943b54371e237927cf59395c2c287e0969ef5f --- /dev/null +++ b/arbor/util/mergeview.hpp @@ -0,0 +1,163 @@ +#pragma once + +#include <iterator> +#include <stdexcept> +#include <type_traits> +#include <utility> + +#include "util/iterutil.hpp" +#include "util/meta.hpp" +#include "util/range.hpp" + +namespace arb { +namespace util { + +template <typename A, typename B, typename = void> +struct has_common_type: std::false_type {}; + +template <typename A, typename B> +struct has_common_type<A, B, util::void_t<std::common_type_t<A, B>>>: std::true_type {}; + +template <typename A, typename B, typename X, typename = void> +struct common_type_or_else { using type = X; }; + +template <typename A, typename B, typename X> +struct common_type_or_else<A, B, X, util::void_t<std::common_type_t<A, B>>> { + using type = std::common_type_t<A, B>; +}; + +template <typename A, typename B, typename X> +using common_type_or_else_t = typename common_type_or_else<A, B, X>::type; + +template <typename LeftI, typename LeftS, typename RightI, typename RightS> +class merge_iterator { + mutable LeftI left_; + mutable RightI right_; + + LeftS left_sentinel_; + RightS right_sentinel_; + + mutable enum which { + undecided, left_next, right_next, done + } next_ = undecided; + + void resolve_next() const { + if (next_!=undecided) return; + + if (left_==left_sentinel_) { + next_ = right_==right_sentinel_? done: right_next; + } + else if (right_==right_sentinel_) { + next_ = left_next; + } + else { + next_ = *left_<*right_? left_next: right_next; + } + } + +public: + using value_type = std::common_type_t< + typename std::iterator_traits<LeftI>::value_type, + typename std::iterator_traits<RightI>::value_type>; + + using pointer = common_type_or_else_t<LeftI, RightI, pointer_proxy<value_type>>; + + using reference = std::conditional_t< + has_common_type<LeftI, RightI>::value, + typename std::iterator_traits<common_type_or_else_t<LeftI, RightI, char*>>::reference, + value_type>; + + using difference_type = typename std::iterator_traits<LeftI>::difference_type; + using iterator_category = std::forward_iterator_tag; + + merge_iterator(): next_(done) {} + + template <typename LSeq, typename RSeq> + merge_iterator(LSeq&& lseq, RSeq&& rseq) { + using std::begin; + using std::end; + + left_ = begin(lseq); + left_sentinel_ = end(lseq); + right_ = begin(rseq); + right_sentinel_ = end(rseq); + } + + bool operator==(const merge_iterator& x) const { + resolve_next(); + x.resolve_next(); + + if (next_==done && x.next_==done) return true; + return next_==x.next_ && left_==x.left_ && right_==x.right_; + } + + bool operator!=(const merge_iterator& x) const { return !(*this==x); } + + reference operator*() const { + resolve_next(); + + switch (next_) { + case left_next: return *left_; + case right_next: return *right_; + default: + throw std::range_error("derefence past end of sequence"); + } + } + + template <typename L = LeftI, typename R = RightI, + typename = std::enable_if_t<has_common_type<L, R>::value>> + std::common_type_t<L, R> operator->() const { + resolve_next(); + return next_==left_next? left_: right_; + } + + template <typename L = LeftI, typename R = RightI, + typename = std::enable_if_t<!has_common_type<L, R>::value>> + pointer_proxy<value_type> operator->() const { + return pointer_proxy<value_type>(*this); + } + + merge_iterator& operator++() { + resolve_next(); + + switch (next_) { + case undecided: + throw std::logic_error("internal error: unexpected undecided"); + case done: return *this; + case left_next: + ++left_; + break; + case right_next: + ++right_; + break; + } + next_ = undecided; + return *this; + } + + merge_iterator& operator++(int) { + auto me = *this; + ++*this; + return me; + } + + merge_iterator& operator=(const merge_iterator&) = default; +}; + +template <typename LSeq, typename RSeq> +auto merge_view(LSeq&& left, RSeq&& right) { + using std::begin; + using std::end; + using LIter = decltype(begin(left)); + using LSentinel = decltype(end(left)); + using RIter = decltype(begin(right)); + using RSentinel = decltype(end(right)); + using MIter = merge_iterator<LIter, LSentinel, RIter, RSentinel>; + + return make_range(MIter(left, right), MIter()); +} + + +} // namespace util +} // namespace arb + diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index f851b7d01190d3f3478e58f96a1a5fc13563ce15..cafd9b78c5371ec462a5ce0eec998e434ee0bed9 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -112,6 +112,7 @@ set(unit_sources test_mechcat.cpp test_mechinfo.cpp test_merge_events.cpp + test_merge_view.cpp test_morphology.cpp test_morph_embedding.cpp test_morph_expr.cpp diff --git a/test/unit/morph_pred.hpp b/test/unit/morph_pred.hpp new file mode 100644 index 0000000000000000000000000000000000000000..ebe57bd95bd9dbdb4450c51f87286d624f1e32da --- /dev/null +++ b/test/unit/morph_pred.hpp @@ -0,0 +1,84 @@ +#pragma once + +// Predicates for morphology testing. + +#include "../gtest.h" + +#include <arbor/morph/morphology.hpp> +#include <arbor/morph/primitives.hpp> +#include <arbor/morph/region.hpp> + +#include "util/span.hpp" + +namespace testing { + +inline ::testing::AssertionResult mlocation_eq(arb::mlocation a, arb::mlocation b) { + if (a.branch!=b.branch) { + return ::testing::AssertionFailure() + << "cables " << a << " and " << b << " differ"; + } + + using FP = testing::internal::FloatingPoint<double>; + if (FP(a.pos).AlmostEquals(FP(b.pos))) { + return ::testing::AssertionSuccess(); + } + else { + return ::testing::AssertionFailure() + << "mlocations " << a << " and " << b << " differ"; + } +} + +inline ::testing::AssertionResult cable_eq(arb::mcable a, arb::mcable b) { + if (a.branch!=b.branch) { + return ::testing::AssertionFailure() + << "cables " << a << " and " << b << " differ"; + } + + using FP = testing::internal::FloatingPoint<double>; + if (FP(a.prox_pos).AlmostEquals(FP(b.prox_pos)) && FP(a.dist_pos).AlmostEquals(FP(b.dist_pos))) { + return ::testing::AssertionSuccess(); + } + else { + return ::testing::AssertionFailure() + << "cables " << a << " and " << b << " differ"; + } +} + +inline ::testing::AssertionResult cablelist_eq(const arb::mcable_list& as, const arb::mcable_list& bs) { + if (as.size()!=bs.size()) { + return ::testing::AssertionFailure() + << "cablelists " << as << " and " << bs << " differ"; + } + + for (auto i: arb::util::count_along(as)) { + auto result = cable_eq(as[i], bs[i]); + if (!result) return ::testing::AssertionFailure() + << "cablelists " << as << " and " << bs << " differ"; + } + return ::testing::AssertionSuccess(); +} + +inline ::testing::AssertionResult extent_eq(const arb::mextent& xa, const arb::mextent& xb) { + return cablelist_eq(xa.cables(), xb.cables()); +} + +inline ::testing::AssertionResult region_eq(const arb::mprovider& p, arb::region a, arb::region b) { + return extent_eq(thingify(a, p), thingify(b, p)); +} + +inline ::testing::AssertionResult mlocationlist_eq(const arb::mlocation_list& as, const arb::mlocation_list& bs) { + if (as.size()!=bs.size()) { + return ::testing::AssertionFailure() + << "cablelists " << as << " and " << bs << " differ"; + } + + for (auto i: arb::util::count_along(as)) { + auto result = mlocation_eq(as[i], bs[i]); + if (!result) return ::testing::AssertionFailure() + << "mlocation lists " << as << " and " << bs << " differ"; + } + return ::testing::AssertionSuccess(); +} + +} // namespace testing + diff --git a/test/unit/test_cv_geom.cpp b/test/unit/test_cv_geom.cpp index 0a14b007e55e503663c8cd0cfed5c1eddadab542..0a0fdbaefa4f2df2f6ab6a94df9bb00f0ce4c192 100644 --- a/test/unit/test_cv_geom.cpp +++ b/test/unit/test_cv_geom.cpp @@ -30,6 +30,10 @@ TEST(cv_geom, empty) { EXPECT_EQ(1u, geom.n_cell()); // can have no CVs but >0 cells. } +static bool region_eq(const mprovider& p, region a, region b) { + return thingify(a, p)==thingify(b, p); +} + TEST(cv_geom, trivial) { using namespace common_morphology; @@ -57,8 +61,8 @@ TEST(cv_geom, trivial) { EXPECT_EQ(geom1.cv_cables, geom4.cv_cables); } - mcable_list all_cables = thingify(reg::all(), cell.provider()); - EXPECT_TRUE(testing::seq_eq(all_cables, geom1.cables(0))); + mcable_list geom1_cables = util::assign_from(geom1.cables(0)); + EXPECT_TRUE(region_eq(cell.provider(), reg::all(), geom1_cables)); } } diff --git a/test/unit/test_merge_view.cpp b/test/unit/test_merge_view.cpp new file mode 100644 index 0000000000000000000000000000000000000000..427b1ea40d75fb0143f755a82b1493801c29e2dc --- /dev/null +++ b/test/unit/test_merge_view.cpp @@ -0,0 +1,111 @@ +#include "../gtest.h" + +#include <forward_list> +#include <vector> + +#include <util/mergeview.hpp> + +#include "common.hpp" + +using namespace arb; + +static auto cstr_range(const char* b) { return util::make_range(b, testing::null_terminated); } + +TEST(mergeview, ctor_eq_iter) { + using std::begin; + using std::end; + + // Make view from sequences with the same iterator types. + { + int a[3] = {1, 3, 5}; + int b[3] = {2, 4, 7}; + + auto merged = util::merge_view(a, b); + EXPECT_TRUE((std::is_same<decltype(merged.begin())::pointer, int*>::value)); + } + + { + const char c[5] = "fish"; + const char* x = "cakes"; + + auto merged = util::merge_view(c, cstr_range(x)); + EXPECT_TRUE((std::is_same<decltype(merged.begin())::pointer, const char*>::value)); + } +} + +TEST(mergeview, ctor_eq_compat) { + // Make view from sequences with compatible iterator types. + { + int a[3] = {1, 3, 5}; + const int b[3] = {2, 4, 7}; + + auto merged = util::merge_view(a, b); + EXPECT_TRUE((std::is_same<decltype(merged.begin())::pointer, const int*>::value)); + } + + { + const std::vector<int> a = {1, 3, 5}; + std::vector<int> b = {2, 4, 6}; + + auto merged = util::merge_view(a, b); + EXPECT_TRUE((std::is_same<decltype(merged.begin())::pointer, std::vector<int>::const_iterator>::value)); + } +} + +TEST(mergeview, ctor_value_compat) { + // Make view from sequences with compatible value types. + int a[3] = {1, 3, 5}; + double b[3] = {2, 4, 7}; + + auto merged = util::merge_view(a, b); + EXPECT_TRUE((std::is_same<std::decay_t<decltype(*merged.begin())>, double>::value)); +} + +TEST(mergeview, empty) { + std::vector<int> a, b; + auto merged = util::merge_view(a, b); + EXPECT_EQ(0, std::distance(merged.begin(), merged.end())); + + b = {1, 3, 4}; + merged = util::merge_view(a, b); + EXPECT_EQ(3, std::distance(merged.begin(), merged.end())); + EXPECT_TRUE(testing::seq_eq(merged, b)); + + std::swap(a, b); + merged = util::merge_view(a, b); + EXPECT_EQ(3, std::distance(merged.begin(), merged.end())); + EXPECT_TRUE(testing::seq_eq(merged, a)); +} + +TEST(mergeview, merge) { + int a[] = {1, 3, 5}; + double b[] = {1., 2.5, 7.}; + double expected[] = {1., 1., 2.5, 3., 5., 7.}; + + EXPECT_TRUE(testing::seq_eq(expected, util::merge_view(a, b))); + EXPECT_TRUE(testing::seq_eq(expected, util::merge_view(b, a))); +} + +TEST(mergeview, iter_ptr) { + struct X { + X(double v): v(v) {} + double v; + bool written = false; + void write(std::vector<double>& buf) { written = true; buf.push_back(v); } + + bool operator<(X b) const { return v<b.v; } + }; + + std::vector<double> out; + X a[] = {2., 4., 4., 6.}; + X b[] = {5., 9.}; + + auto merged = util::merge_view(a, b); + for (auto i = merged.begin(); i!=merged.end(); ++i) { + i->write(out); + } + + EXPECT_EQ((std::vector<double>{2., 4., 4., 5., 6., 9.}), out); + for (auto& x: a) EXPECT_TRUE(x.written); + for (auto& x: b) EXPECT_TRUE(x.written); +} diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp index 9f2d9449fda9801f796f5ed91809d00a8db232b5..ff5ee103e113dde1c54951e1018d374ea398f621 100644 --- a/test/unit/test_morph_expr.cpp +++ b/test/unit/test_morph_expr.cpp @@ -14,75 +14,14 @@ #include "util/span.hpp" #include "util/strprintf.hpp" +#include "morph_pred.hpp" + using namespace arb; using embedding = embed_pwlin; -namespace arb { - namespace reg { - mcable_list remove_covered_points(mcable_list cables, const morphology& m); - - } -} - -::testing::AssertionResult mlocation_eq(mlocation a, mlocation b) { - if (a.branch!=b.branch) { - return ::testing::AssertionFailure() - << "cables " << a << " and " << b << " differ"; - } - - using FP = testing::internal::FloatingPoint<double>; - if (FP(a.pos).AlmostEquals(FP(b.pos))) { - return ::testing::AssertionSuccess(); - } - else { - return ::testing::AssertionFailure() - << "mlocations " << a << " and " << b << " differ"; - } -} - -::testing::AssertionResult cable_eq(mcable a, mcable b) { - if (a.branch!=b.branch) { - return ::testing::AssertionFailure() - << "cables " << a << " and " << b << " differ"; - } - - using FP = testing::internal::FloatingPoint<double>; - if (FP(a.prox_pos).AlmostEquals(FP(b.prox_pos)) && FP(a.dist_pos).AlmostEquals(FP(b.dist_pos))) { - return ::testing::AssertionSuccess(); - } - else { - return ::testing::AssertionFailure() - << "cables " << a << " and " << b << " differ"; - } -} - -::testing::AssertionResult cablelist_eq(const mcable_list& as, const mcable_list& bs) { - if (as.size()!=bs.size()) { - return ::testing::AssertionFailure() - << "cablelists " << as << " and " << bs << " differ"; - } - - for (auto i: util::count_along(as)) { - auto result = cable_eq(as[i], bs[i]); - if (!result) return ::testing::AssertionFailure() - << "cablelists " << as << " and " << bs << " differ"; - } - return ::testing::AssertionSuccess(); -} - -::testing::AssertionResult mloctionlist_eq(const mlocation_list& as, const mlocation_list& bs) { - if (as.size()!=bs.size()) { - return ::testing::AssertionFailure() - << "cablelists " << as << " and " << bs << " differ"; - } - - for (auto i: util::count_along(as)) { - auto result = mlocation_eq(as[i], bs[i]); - if (!result) return ::testing::AssertionFailure() - << "mlocation lists " << as << " and " << bs << " differ"; - } - return ::testing::AssertionSuccess(); -} +using testing::region_eq; +using testing::cablelist_eq; +using testing::mlocationlist_eq; TEST(region, expr_repn) { using util::to_string; @@ -416,27 +355,25 @@ TEST(region, thingify_simple_morphologies) { cl all_{{0, 0, 1}}; cl empty_{}; - EXPECT_EQ(thingify(h1, mp), h1_); - EXPECT_EQ(thingify(h2, mp), h2_); - EXPECT_EQ(thingify(join(h1, h2), mp), all_); - EXPECT_EQ(thingify(intersect(h1, h2), mp), (cl{{0, 0.5, 0.5}})); - - EXPECT_TRUE(cablelist_eq(thingify(t1, mp), t1_)); - EXPECT_TRUE(cablelist_eq(thingify(t2, mp), t2_)); - EXPECT_TRUE(cablelist_eq(thingify(intersect(h1, h1), mp), h1_)); - EXPECT_TRUE(cablelist_eq(thingify(intersect(t1, t1), mp), t1_)); - EXPECT_TRUE(cablelist_eq(thingify(join(t1, t2), mp), all_)); - EXPECT_TRUE(cablelist_eq(thingify(intersect(all, t1), mp), t1_)); - EXPECT_TRUE(cablelist_eq(thingify(intersect(all, t2), mp), t2_)); - EXPECT_TRUE(cablelist_eq(thingify(join(all, t1), mp), all_)); - EXPECT_TRUE(cablelist_eq(thingify(join(all, t2), mp), all_)); - EXPECT_TRUE(cablelist_eq(thingify(join(h1, t1), mp), (cl{{0, 0, 0.7}}))); - EXPECT_TRUE(cablelist_eq(thingify(join(h1, t2), mp), (cl{{0, 0, 0.5}, {0, 0.7, 1}}))); - EXPECT_TRUE(cablelist_eq(thingify(intersect(h2, t1), mp), (cl{{0, 0.5, 0.7}}))); + EXPECT_TRUE(region_eq(mp, join(h1, h2), all_)); + EXPECT_TRUE(region_eq(mp, intersect(h1, h2), cl{{0, 0.5, 0.5}})); + EXPECT_TRUE(region_eq(mp, t1, t1_)); + EXPECT_TRUE(region_eq(mp, t2, t2_)); + EXPECT_TRUE(region_eq(mp, intersect(h1, h1), h1_)); + EXPECT_TRUE(region_eq(mp, intersect(t1, t1), t1_)); + EXPECT_TRUE(region_eq(mp, join(t1, t2), all_)); + EXPECT_TRUE(region_eq(mp, intersect(all, t1), t1_)); + EXPECT_TRUE(region_eq(mp, intersect(all, t2), t2_)); + EXPECT_TRUE(region_eq(mp, join(all, t1), all_)); + EXPECT_TRUE(region_eq(mp, join(all, t2), all_)); + EXPECT_TRUE(region_eq(mp, join(h1, t1), cl{{0, 0, 0.7}})); + EXPECT_TRUE(region_eq(mp, join(h1, t2), cl{{0, 0, 0.5}, {0, 0.7, 1}})); + EXPECT_TRUE(region_eq(mp, intersect(h2, t1), cl{{0, 0.5, 0.7}})); // Check round-trip of implicit region conversions. - EXPECT_EQ((mcable_list{{0, 0.3, 0.6}}), thingify(region(mcable{0, 0.3, 0.6}), mp)); - EXPECT_TRUE(cablelist_eq(t2_, thingify(region(t2_), mp))); + // (No fork points in cables, so extent should not including anyhing extra). + EXPECT_EQ((mcable_list{{0, 0.3, 0.6}}), thingify(region(mcable{0, 0.3, 0.6}), mp).cables()); + EXPECT_TRUE(cablelist_eq(t2_, thingify(region(t2_), mp).cables())); } @@ -484,18 +421,18 @@ TEST(region, thingify_simple_morphologies) { auto reg5_ = distal_interval(start1_, 0); auto reg6_ = proximal_interval(start1_, 0); - EXPECT_EQ(thingify(tagged(1), mp), (mcable_list{{0,0,1}})); - EXPECT_EQ(thingify(tagged(2), mp), (mcable_list{{2,0,1}})); - EXPECT_EQ(thingify(tagged(3), mp), (mcable_list{{1,0,1}})); - EXPECT_EQ(thingify(join(tagged(1), tagged(2), tagged(3)), mp), (mcable_list{{0,0,1}, {1,0,1}, {2,0,1}})); - EXPECT_EQ(thingify(join(tagged(1), tagged(2), tagged(3)), mp), thingify(all(), mp)); - EXPECT_EQ(thingify(reg0_, mp), (mcable_list{{1,0,0.5}})); - EXPECT_EQ(thingify(reg1_, mp), (mcable_list{{0,0.5,1}, {1,0,0.8}, {2,0,0.8}})); - EXPECT_EQ(thingify(reg2_, mp), (mcable_list{{1,0.5,1}})); - EXPECT_EQ(thingify(reg3_, mp), (mcable_list{{0, 0.75, 1}, {1,0,1}})); - EXPECT_EQ(thingify(reg4_, mp), (mcable_list{{1,1,1}})); - EXPECT_EQ(thingify(reg5_, mp), (mcable_list{{0,1,1}})); - EXPECT_EQ(thingify(reg6_, mp), (mcable_list{{0,1,1}})); + EXPECT_TRUE(region_eq(mp, tagged(1), mcable_list{{0,0,1}})); + EXPECT_TRUE(region_eq(mp, tagged(2), mcable_list{{2,0,1}})); + EXPECT_TRUE(region_eq(mp, tagged(3), mcable_list{{1,0,1}})); + EXPECT_TRUE(region_eq(mp, join(tagged(1), tagged(2), tagged(3)), mcable_list{{0,0,1}, {1,0,1}, {2,0,1}})); + EXPECT_TRUE(region_eq(mp, join(tagged(1), tagged(2), tagged(3)), all())); + EXPECT_TRUE(region_eq(mp, reg0_, mcable_list{{1,0,0.5}})); + EXPECT_TRUE(region_eq(mp, reg1_, mcable_list{{0,0.5,1}, {1,0,0.8}, {2,0,0.8}})); + 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}})); } } @@ -570,18 +507,18 @@ TEST(region, thingify_moderate_morphologies) { mcable c_end1_{1,1,1}; mcable c_root_{0,0,0}; - EXPECT_EQ(thingify(all(), mp), all_); - EXPECT_EQ(thingify(axon, mp), (cl{b1_})); - EXPECT_EQ(thingify(dend, mp), (cl{b0_,b3_})); - EXPECT_EQ(thingify(apic, mp), (cl{b2_})); - EXPECT_EQ(thingify(join(dend, apic), mp), (cl{b0_,b2_,b3_})); - EXPECT_EQ(thingify(join(axon, join(dend, apic)), mp), all_); + EXPECT_TRUE(region_eq(mp, all(), all_)); + EXPECT_TRUE(region_eq(mp, axon, cl{b1_})); + EXPECT_TRUE(region_eq(mp, dend, cl{b0_,b3_})); + EXPECT_TRUE(region_eq(mp, apic, cl{b2_})); + EXPECT_TRUE(region_eq(mp, join(dend, apic), cl{b0_,b2_,b3_})); + EXPECT_TRUE(region_eq(mp, join(axon, join(dend, apic)), all_)); // Test that intersection correctly generates zero-length cables at // parent-child interfaces. - EXPECT_EQ(thingify(intersect(apic, dend), mp), (cl{c_end1_})); - EXPECT_EQ(thingify(intersect(apic, axon), mp), (cl{c_end1_})); - EXPECT_EQ(thingify(intersect(axon, dend), mp), (cl{c_root_, c_end1_})); + 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_})); // Test distal and proximal interavls auto start0_ = location(0, 0 ); @@ -600,39 +537,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(cablelist_eq(thingify(distal_interval(start0_, 1000), mp), (mcable_list{{0,0,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(distal_interval(quar_1_, 150), mp), (mcable_list{{1,0.25,1}, {2,0,0.75}, {3,0,0.375}}))); - EXPECT_TRUE(cablelist_eq(thingify(distal_interval(mid1_, 1000), mp), (mcable_list{{1,0.5,1}, {2,0,1}, {3,0,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(distal_interval(mid1_, 150), mp), (mcable_list{{1,0.5,1}, {2,0,1}, {3,0,0.5}}))); - EXPECT_TRUE(cablelist_eq(thingify(distal_interval(end1_, 100), mp), (mcable_list{{2,0,1},{3,0,0.5}}))); - EXPECT_TRUE(cablelist_eq(thingify(distal_interval(join(quar_1_, mid1_), 150), mp), (mcable_list{{1,0.25,1}, {2,0,1}, {3,0,0.5}}))); - EXPECT_TRUE(cablelist_eq(thingify(distal_interval(join(quar_1_, loc_3_1_), 150), mp), (mcable_list{{1,0.25,1}, {2,0,0.75}, {3,0,0.375}, {3,0.65,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(distal_interval(join(quar_1_, loc_3_1_), 150), mp), (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), 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}})); // Proximal from point and/or interval - EXPECT_TRUE(cablelist_eq(thingify(proximal_interval(mid3_, 100), mp), (mcable_list{{3,0,0.5}}))); - EXPECT_TRUE(cablelist_eq(thingify(proximal_interval(mid3_, 150), mp), (mcable_list{{1,0.5,1}, {3,0,0.5}}))); - EXPECT_TRUE(cablelist_eq(thingify(proximal_interval(end2_, 150), mp), (mcable_list{{1,0.5,1}, {2,0,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(proximal_interval(end2_, 500), mp), (mcable_list{{1,0,1}, {2,0,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(proximal_interval(loc_3_0_, 100), mp), (mcable_list{{1,0.8,1}, {3,0,0.4}}))); - EXPECT_TRUE(cablelist_eq(thingify(proximal_interval(join(loc_3_0_, mid2_), 120), mp), (mcable_list{{1,0.3,1}, {2,0,0.5}, {3, 0, 0.4}}))); + 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}})); // Test radius_lt and radius_gt - EXPECT_TRUE(cablelist_eq(thingify(radius_lt(all(), 2), mp), (mcable_list{{0,0,0.55}, {1,0,0.325}, {3,0.375,0.75}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_lt(all(), 3), mp), (mcable_list{{0,0,1}, {1,0,0.55}, {2,6.0/9.0,1}, {3,0.25,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_gt(all(), 2), mp), (mcable_list{{0,0.55,1}, {1,0.325,1}, {2,0,1}, {3,0,0.375}, {3,0.75,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_gt(all(), 3), mp), (mcable_list{{1,0.55,1}, {2,0,6.0/9.0}, {3,0,0.25}}))); - - EXPECT_TRUE(cablelist_eq(thingify(radius_le(all(), 2), mp), (mcable_list{{0,0,0.55}, {1,0,0.325}, {2,1,1}, {3,0.375,0.75}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_le(all(), 3), mp), (mcable_list{{0,0,1}, {1,0,0.55}, {2,6.0/9.0,1}, {3,0.25,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_ge(all(), 2), mp), (mcable_list{{0,0.55,1}, {1,0.325,1}, {2,0,1}, {3,0,0.375}, {3,0.75,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_ge(all(), 3), mp), (mcable_list{{1,0.55,1}, {2,0,6.0/9.0}, {3,0,0.25}}))); - - EXPECT_TRUE(cablelist_eq(thingify(radius_lt(reg_a_, 2), mp), (mcable_list{{0,0.1,0.4},{3,0.375,0.4}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_gt(reg_a_, 2), mp), (mcable_list{{2,0,1},{3,0.1,0.375}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_lt(reg_b_, 2), mp), (mcable_list{{0,0.1,0.4}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_gt(reg_c_, 2), mp), (mcable_list{{0,0.55,0.7},{2,0,0.5},{3,0.1,0.375},{3,0.9,1}}))); - EXPECT_TRUE(cablelist_eq(thingify(radius_gt(reg_d_, 2), mp), (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), 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}})); // Test some more interesting intersections and unions. @@ -643,15 +580,15 @@ TEST(region, thingify_moderate_morphologies) { // |xxxxxxxxx|xxxxxxxxx| ror auto lhs = b13; auto rhs = join(cable(1,.2,.7), cable(3,.3,.6)); - auto rand = cl{ {1,.2,.7}, {3,.3,.6}}; - auto ror = cl{ {1,.0,1.}, {3,.0,1.}}; - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + auto rand = cl{{1,.2,.7}, {3,.3,.6}}; + auto ror = cl{{1,.0,1.}, {3,.0,1.}}; + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // Assert communtativity std::swap(lhs, rhs); - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // 123456789 123456789 // | ----- | ---- | lhs @@ -662,13 +599,13 @@ TEST(region, thingify_moderate_morphologies) { rhs = join(cable(1,.2,.7), cable(3,.3,.6)); rand = cl{ {1,.3,.7}, {3,.3,.5}}; ror = cl{ {1,.2,.8}, {3,.1,.6}}; - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // Assert communtativity std::swap(lhs, rhs); - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // 123456789 123456789 // | -- - | --- --- | lhs @@ -679,13 +616,13 @@ TEST(region, thingify_moderate_morphologies) { rhs = join(cable(1,.2,.7), cable(3,.3,.6)); rand = cl{ {1,.2,.3}, {1,.4,.5}, {3,.3,.4}, {3,.5,.6}}; ror = cl{ {1,.1,.7}, {3,.1,.9}}; - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // Assert communtativity std::swap(lhs, rhs); - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // b1 // 123456789 @@ -697,8 +634,8 @@ TEST(region, thingify_moderate_morphologies) { rhs = cable(1,0,.5); rand = cl{{1,0,.5}}; ror = cl{{1,0,.5}}; - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // b3 // 123456789 @@ -710,8 +647,8 @@ TEST(region, thingify_moderate_morphologies) { rhs = cable(3,0,.5); rand = cl{{3,0,.5}}; ror = cl{{3,0,.5}}; - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // b0 b1 // 123456789 123456789 @@ -723,13 +660,13 @@ TEST(region, thingify_moderate_morphologies) { rhs = cable(1,0,.5); rand = cl{{0,0,0}}; ror = cl{{0,0,.5},{1,0,.5}}; - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // Assert communtativity std::swap(lhs, rhs); - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // b2 b3 // 123456789 123456789 @@ -741,13 +678,13 @@ TEST(region, thingify_moderate_morphologies) { rhs = cable(3,0,.5); rand = cl{{1,1,1}}; ror = cl{{2,0,.5},{3,0,.5}}; - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // Assert communtativity std::swap(lhs, rhs); - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // b0 b3 // 123456789 123456789 @@ -759,20 +696,19 @@ TEST(region, thingify_moderate_morphologies) { rhs = join(cable(0,0,.7), cable(3,0,.3)); rand = cl{{0,0,.5},{3,0,.3}}; ror = cl{{0,0,.7},{3,0,.5}}; - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); // Assert communtativity std::swap(lhs, rhs); - EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand); - EXPECT_EQ(thingify(join(lhs, rhs), mp), ror); + EXPECT_TRUE(region_eq(mp, intersect(lhs, rhs), rand)); + EXPECT_TRUE(region_eq(mp, join(lhs, rhs), ror)); } } TEST(region, thingify_complex_morphologies) { using pvec = std::vector<msize_t>; using svec = std::vector<msample>; - using cl = mcable_list; { pvec parents = {mnpos, 0, 1, 0, 3, 4, 5, 5, 7, 7, 4, 10}; svec samples = { @@ -791,20 +727,6 @@ TEST(region, thingify_complex_morphologies) { }; sample_tree sm(samples, parents); auto m = morphology(sm, false); - { - auto in = cl{{0,0,0},{1,0,0.5},{1,1,1},{2,0,1},{2,1,1},{3,1,1},{4,0,1},{5,1,1},{7,0,1}}; - auto out = reg::remove_covered_points(in, m); - - auto expected = cl{{1,0,0.5},{2,0,1},{3,1,1},{4,0,1},{5,1,1},{7,0,1}}; - EXPECT_TRUE(cablelist_eq(out, expected)); - } - { - auto in = cl{{0,0,0},{1,0,0.5},{1,1,1},{2,1,1},{3,1,1},{4,0,1},{5,1,1},{7,0,1}}; - auto out = reg::remove_covered_points(in, m); - - auto expected = cl{{1,0,0.5},{3,1,1},{4,0,1},{5,1,1},{7,0,1}}; - EXPECT_TRUE(cablelist_eq(out, expected)); - } { mprovider mp(m); using reg::cable; @@ -818,19 +740,19 @@ TEST(region, thingify_complex_morphologies) { 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)); - EXPECT_TRUE(mloctionlist_eq(thingify(most_distal(reg_a_), mp), mlocation_list{{0,0.9},{1,0.4}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_distal(reg_b_), mp), mlocation_list{{0,0.9},{1,0.5}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_distal(reg_c_), mp), mlocation_list{{0,0.9},{2,0.5}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_distal(reg_d_), mp), mlocation_list{{3,0.1},{4,0.6}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_distal(reg_e_), mp), mlocation_list{{5,0.6}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_distal(reg_f_), mp), mlocation_list{{5,0.6},{7,1}})); - - EXPECT_TRUE(mloctionlist_eq(thingify(most_proximal(reg_a_), mp), mlocation_list{{0,0}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_proximal(reg_b_), mp), mlocation_list{{0,0}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_proximal(reg_c_), mp), mlocation_list{{0,0}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_proximal(reg_d_), mp), mlocation_list{{2,0}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_proximal(reg_e_), mp), mlocation_list{{2,0}})); - EXPECT_TRUE(mloctionlist_eq(thingify(most_proximal(reg_f_), mp), mlocation_list{{2,0}})); + 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}})); + EXPECT_TRUE(mlocationlist_eq(thingify(most_distal(reg_c_), mp), mlocation_list{{0,0.9},{2,0.5}})); + EXPECT_TRUE(mlocationlist_eq(thingify(most_distal(reg_d_), mp), mlocation_list{{3,0.1},{4,0.6}})); + 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}})); } } { @@ -854,6 +776,7 @@ TEST(region, thingify_complex_morphologies) { mprovider mp(morphology(sm, false)); using reg::all; + using reg::nil; using reg::z_dist_from_root_lt; using reg::z_dist_from_root_le; using reg::z_dist_from_root_gt; @@ -861,48 +784,45 @@ TEST(region, thingify_complex_morphologies) { using reg::cable; // Test projection - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_lt(0), mp), (mcable_list{}))); - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_ge(0), mp), thingify(all(), mp))); - - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_le(100), mp), thingify(all(), mp))); - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_gt(100), mp), (mcable_list{}))); - - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_le(90), mp), thingify(all(), mp))); - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_gt(90), mp), (mcable_list{}))); - - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_lt(20), mp), - (mcable_list{{0,0,1}, - {1,0,0.578250901781922829}, - {2,0.61499300915417734997,0.8349970039232188642}, - {3,0,0.179407353580315756} - }))); - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_ge(20), mp), - (mcable_list{{0,1,1}, - {1,0.578250901781922829,1}, - {2,0,0.61499300915417734997}, - {2,0.8349970039232188642,1}, - {3,0.179407353580315756,1}, - {4,0,1}, - {5,0,1} - }))); - EXPECT_TRUE(cablelist_eq(thingify(join(z_dist_from_root_lt(20), z_dist_from_root_ge(20)), mp), thingify(all(), mp))); - - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_le(50), mp), - (mcable_list{{0,0,1}, - {1,0,1}, - {2,0,0.2962417607888518767}, - {2,0.4499900130773962142,1}, - {3,0,0.4485183839507893905}, - {3,0.7691110303704736343,1}, - {4,0,0.0869615364994152821}, - {5,0,0.25} - }))); - EXPECT_TRUE(cablelist_eq(thingify(z_dist_from_root_gt(50), mp), - (mcable_list{{2,0.2962417607888518767,0.4499900130773962142}, - {3,0.4485183839507893905,0.7691110303704736343}, - {4,0.0869615364994152821,1}, - {5,0.25,1}}))); - - EXPECT_TRUE(cablelist_eq(thingify(join(z_dist_from_root_le(50), z_dist_from_root_gt(50)), mp), thingify(all(), mp))); + EXPECT_TRUE(region_eq(mp, z_dist_from_root_lt(0), nil())); + EXPECT_TRUE(region_eq(mp, z_dist_from_root_gt(0), all())); + + EXPECT_TRUE(region_eq(mp, z_dist_from_root_le(100), all())); + EXPECT_TRUE(region_eq(mp, z_dist_from_root_gt(100), nil())); + + EXPECT_TRUE(region_eq(mp, z_dist_from_root_le(90), all())); + EXPECT_TRUE(region_eq(mp, z_dist_from_root_gt(90), nil())); + + EXPECT_TRUE(region_eq(mp, z_dist_from_root_lt(20), + mcable_list{{0,0,1}, + {1,0,0.578250901781922829}, + {2,0.61499300915417734997,0.8349970039232188642}, + {3,0,0.179407353580315756}})); + EXPECT_TRUE(region_eq(mp, z_dist_from_root_ge(20), + mcable_list{{0,1,1}, + {1,0.578250901781922829,1}, + {2,0,0.61499300915417734997}, + {2,0.8349970039232188642,1}, + {3,0.179407353580315756,1}, + {4,0,1}, + {5,0,1}})); + EXPECT_TRUE(region_eq(mp, join(z_dist_from_root_lt(20), z_dist_from_root_ge(20)), all())); + + EXPECT_TRUE(region_eq(mp, z_dist_from_root_le(50), + mcable_list{{0,0,1}, + {1,0,1}, + {2,0,0.2962417607888518767}, + {2,0.4499900130773962142,1}, + {3,0,0.4485183839507893905}, + {3,0.7691110303704736343,1}, + {4,0,0.0869615364994152821}, + {5,0,0.25}})); + EXPECT_TRUE(region_eq(mp, z_dist_from_root_gt(50), + mcable_list{{2,0.2962417607888518767,0.4499900130773962142}, + {3,0.4485183839507893905,0.7691110303704736343}, + {4,0.0869615364994152821,1}, + {5,0.25,1}})); + + EXPECT_TRUE(region_eq(mp, join(z_dist_from_root_le(50), z_dist_from_root_gt(50)), all())); } } diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp index cfe20b754594932d8e6a992185301c3c59a21f42..5bf310274a867acbb928d1d3ae32135cfb0392f8 100644 --- a/test/unit/test_morphology.cpp +++ b/test/unit/test_morphology.cpp @@ -13,6 +13,8 @@ #include "morph/mbranch.hpp" #include "util/span.hpp" +#include "morph_pred.hpp" + // Test basic functions on properties of mpoint TEST(morphology, mpoint) { using mp = arb::mpoint; @@ -631,3 +633,104 @@ TEST(morphology, minset) { } } +// Testing mextent; intersection and join operations are +// exercised by region/locset thingifies in test_morph_expr.cpp. + +TEST(mextent, invariants) { + 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{{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}}); + ASSERT_TRUE(x2.test_invariants(m)); + EXPECT_TRUE(cablelist_eq(cl{{0, 0, 0}, {1, 0., 0.75}}, x2.cables())); + + mextent x3(m, 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())); +} + +TEST(mextent, intersects) { + 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{{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})); + EXPECT_FALSE(x1.intersects(mlocation{1, 0.})); + EXPECT_FALSE(x1.intersects(mlocation{1, 1.})); + + EXPECT_FALSE(x1.intersects(mcable{1, 0., 0.2})); + EXPECT_TRUE(x1.intersects(mcable{1, 0., 0.25})); + EXPECT_TRUE(x1.intersects(mcable{1, 0.4, 0.6})); + EXPECT_TRUE(x1.intersects(mcable{1, 0.2, 0.8})); + 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}}); + EXPECT_TRUE(x2.intersects(mlocation{1, 0.})); + EXPECT_TRUE(x2.intersects(mlocation{0, 0.})); + EXPECT_FALSE(x2.intersects(mlocation{0, 1.})); + + mextent x3(m, 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{3, 1.})); +}