diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 42a66fc7693b13594690bbcb0a4c97f2e7b0d257..a6a497a9d21c0bb9ab0ea718f4138e43e2115c11 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -15,6 +15,7 @@ set(arbor_sources
cable_cell_param.cpp
cell_group_factory.cpp
common_types_io.cpp
+ cv_policy.cpp
execution_context.cpp
gpu_context.cpp
event_binner.cpp
diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp
index 417229bb7e296fabe9b93f8ccc05b76fbf3d5508..b20a93a1b9067e3d0ef36106603380d8ff05f9a3 100644
--- a/arbor/cable_cell_param.cpp
+++ b/arbor/cable_cell_param.cpp
@@ -5,11 +5,8 @@
#include <arbor/cable_cell.hpp>
#include <arbor/cable_cell_param.hpp>
-#include <arbor/morph/locset.hpp>
#include "util/maputil.hpp"
-#include "util/rangeutil.hpp"
-#include "util/span.hpp"
namespace arb {
@@ -71,109 +68,4 @@ cable_cell_parameter_set neuron_parameter_defaults = {
},
};
-// Discretization policy implementations:
-
-locset cv_policy_max_extent::cv_boundary_points(const cable_cell& cell) const {
- const unsigned nbranch = cell.morphology().num_branches();
- const auto& embed = cell.embedding();
- if (!nbranch || max_extent_<=0) return ls::nil();
-
- std::vector<mlocation> points;
-
- unsigned bidx = 0;
- if (flags_&cv_policy_flag::single_root_cv) {
- points.push_back({0, 0.});
- points.push_back({0, 1.});
- bidx = 1;
- }
-
- const double oomax_extent = 1./max_extent_;
- while (bidx<nbranch) {
- unsigned ncv = std::ceil(embed.branch_length(bidx)*oomax_extent);
- double ooncv = 1./ncv;
-
- if (flags_&cv_policy_flag::interior_forks) {
- for (unsigned i = 0; i<ncv; ++i) {
- points.push_back({bidx, (1+2*i)*ooncv/2});
- }
- }
- else {
- for (unsigned i = 0; i<ncv; ++i) {
- points.push_back({bidx, i*ooncv});
- }
- points.push_back({bidx, 1.});
- }
- ++bidx;
- }
-
- util::sort(points);
- return points;
-}
-
-locset cv_policy_fixed_per_branch::cv_boundary_points(const cable_cell& cell) const {
- const unsigned nbranch = cell.morphology().num_branches();
- if (!nbranch) return ls::nil();
-
- std::vector<mlocation> points;
-
- unsigned bidx = 0;
- if (flags_&cv_policy_flag::single_root_cv) {
- points.push_back({0, 0.});
- points.push_back({0, 1.});
- bidx = 1;
- }
-
- double ooncv = 1./cv_per_branch_;
- while (bidx<nbranch) {
- if (flags_&cv_policy_flag::interior_forks) {
- for (unsigned i = 0; i<cv_per_branch_; ++i) {
- points.push_back({bidx, (1+2*i)*ooncv/2});
- }
- }
- else {
- for (unsigned i = 0; i<cv_per_branch_; ++i) {
- points.push_back({bidx, i*ooncv});
- }
- points.push_back({bidx, 1.});
- }
- ++bidx;
- }
-
- util::sort(points);
- return points;
-}
-
-locset cv_policy_every_sample::cv_boundary_points(const cable_cell& cell) const {
- const unsigned nbranch = cell.morphology().num_branches();
- if (!nbranch) return ls::nil();
-
- bool single_root = cell.morphology().spherical_root() || (flags_&cv_policy_flag::single_root_cv);
-
- // Ignore interior_forks flag, but if single_root_cv is set, take sample indices only from branches 1+.
- // Always include branch proximal points, so that forks are trivial.
-
-
- if (single_root) {
- std::vector<msize_t> samples;
- for (unsigned i = 1; i<nbranch; ++i) {
- util::append(samples, util::make_range(cell.morphology().branch_indexes(i)));
- }
- util::sort(samples);
- samples.erase(std::unique(samples.begin(), samples.end()), samples.end());
-
- return join(
- ls::on_branches(0.),
- ls::location(0, 1.),
- std::accumulate(samples.begin(), samples.end(), ls::nil(),
- [](auto&& l, auto&& r) { return sum(std::move(l), ls::sample(r)); }));
- }
- else {
- auto samples = util::make_span(cell.morphology().num_samples());
- return join(
- ls::on_branches(0.),
- std::accumulate(samples.begin(), samples.end(), ls::nil(),
- [](auto&& l, auto&& r) { return sum(std::move(l), ls::sample(r)); }));
- }
-}
-
} // namespace arb
diff --git a/arbor/cv_policy.cpp b/arbor/cv_policy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fa39de7d481c78f0eac865f6b763498636eb1b34
--- /dev/null
+++ b/arbor/cv_policy.cpp
@@ -0,0 +1,158 @@
+#include <utility>
+#include <vector>
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/cv_policy.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/region.hpp>
+
+#include "util/rangeutil.hpp"
+#include "util/span.hpp"
+
+// Discretization policy implementations:
+
+namespace arb {
+
+static auto unique_sum = [](auto&&... lss) {
+ return ls::support(sum(std::forward<decltype(lss)>(lss)...));
+};
+
+// Combinators:
+// cv_policy_plus_ represents the result of operator+,
+// cv_policy_bar_ represents the result of operator|.
+
+struct cv_policy_plus_: cv_policy_base {
+ cv_policy_plus_(const cv_policy& lhs, const cv_policy& rhs):
+ lhs_(lhs), rhs_(rhs) {}
+
+ cv_policy_base_ptr clone() const override {
+ return cv_policy_base_ptr(new cv_policy_plus_(*this));
+ }
+
+ locset cv_boundary_points(const cable_cell& c) const override {
+ return unique_sum(lhs_.cv_boundary_points(c), rhs_.cv_boundary_points(c));
+ }
+
+ region domain() const override { return join(lhs_.domain(), rhs_.domain()); }
+
+ cv_policy lhs_, rhs_;
+};
+
+cv_policy operator+(const cv_policy& lhs, const cv_policy& rhs) {
+ return cv_policy_plus_(lhs, rhs);
+}
+
+struct cv_policy_bar_: cv_policy_base {
+ cv_policy_bar_(const cv_policy& lhs, const cv_policy& rhs):
+ lhs_(lhs), rhs_(rhs) {}
+
+ cv_policy_base_ptr clone() const override {
+ return cv_policy_base_ptr(new cv_policy_bar_(*this));
+ }
+
+ locset cv_boundary_points(const cable_cell& c) const override {
+ return unique_sum(ls::restrict(lhs_.cv_boundary_points(c), complement(rhs_.domain())), rhs_.cv_boundary_points(c));
+ }
+
+ region domain() const override { return join(lhs_.domain(), rhs_.domain()); }
+
+ cv_policy lhs_, rhs_;
+};
+
+cv_policy operator|(const cv_policy& lhs, const cv_policy& rhs) {
+ return cv_policy_bar_(lhs, rhs);
+}
+
+// Public policy implementations:
+
+locset cv_policy_explicit::cv_boundary_points(const cable_cell& cell) const {
+ return
+ ls::support(
+ util::foldl(
+ [this](locset l, const auto& comp) {
+ return sum(std::move(l), ls::restrict(locs_, comp));
+ },
+ ls::boundary(domain_),
+ components(cell.morphology(), thingify(domain_, cell.provider()))));
+}
+
+locset cv_policy_max_extent::cv_boundary_points(const cable_cell& cell) const {
+ const unsigned nbranch = cell.morphology().num_branches();
+ const auto& embed = cell.embedding();
+ if (!nbranch || max_extent_<=0) return ls::nil();
+
+ std::vector<mlocation> points;
+ double oomax_extent = 1./max_extent_;
+ auto comps = components(cell.morphology(), thingify(domain_, cell.provider()));
+
+ for (auto& comp: comps) {
+ for (mcable c: comp) {
+ double cable_length = embed.integrate_length(c);
+ unsigned ncv = std::ceil(cable_length*oomax_extent);
+ double scale = (c.dist_pos-c.prox_pos)/ncv;
+
+ if (flags_&cv_policy_flag::interior_forks) {
+ for (unsigned i = 0; i<ncv; ++i) {
+ points.push_back({c.branch, c.prox_pos+(1+2*i)*scale/2});
+ }
+ }
+ else {
+ for (unsigned i = 0; i<ncv; ++i) {
+ points.push_back({c.branch, c.prox_pos+i*scale});
+ }
+ points.push_back({c.branch, c.dist_pos});
+ }
+ }
+ }
+
+ util::sort(points);
+ return unique_sum(locset(std::move(points)), ls::cboundary(domain_));
+}
+
+locset cv_policy_fixed_per_branch::cv_boundary_points(const cable_cell& cell) const {
+ const unsigned nbranch = cell.morphology().num_branches();
+ if (!nbranch) return ls::nil();
+
+ std::vector<mlocation> points;
+ double ooncv = 1./cv_per_branch_;
+ auto comps = components(cell.morphology(), thingify(domain_, cell.provider()));
+
+ for (auto& comp: comps) {
+ for (mcable c: comp) {
+ double scale = (c.dist_pos-c.prox_pos)*ooncv;
+
+ if (flags_&cv_policy_flag::interior_forks) {
+ for (unsigned i = 0; i<cv_per_branch_; ++i) {
+ points.push_back({c.branch, c.prox_pos+(1+2*i)*scale/2});
+ }
+ }
+ else {
+ for (unsigned i = 0; i<cv_per_branch_; ++i) {
+ points.push_back({c.branch, c.prox_pos+i*scale});
+ }
+ points.push_back({c.branch, c.dist_pos});
+ }
+ }
+ }
+
+ util::sort(points);
+ return unique_sum(locset(std::move(points)), ls::cboundary(domain_));
+}
+
+locset cv_policy_every_sample::cv_boundary_points(const cable_cell& cell) const {
+ const unsigned nbranch = cell.morphology().num_branches();
+ if (!nbranch) return ls::nil();
+
+ // Always include branch proximal points, so that forks are trivial.
+
+ return
+ unique_sum(ls::cboundary(domain_),
+ ls::restrict(
+ util::foldl(
+ [](locset l, msize_t sidx) { return sum(std::move(l), ls::sample(sidx)); },
+ ls::on_branches(0),
+ util::make_span(cell.morphology().num_samples())),
+ domain_));
+}
+
+} // namespace arb
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index 2e041e7dfd4e823984fe297cb269d6352a06cd56..804ab131cfe3a2820ca28ddef47394f388955c3f 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -1,6 +1,7 @@
#pragma once
#include <arbor/arbexcept.hpp>
+#include <arbor/cv_policy.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/util/optional.hpp>
@@ -144,151 +145,6 @@ struct ion_reversal_potential_method {
mechanism_desc method;
};
-// FVM discretization policies/hints.
-//
-// CV polices, given a cable cell, provide a locset comprising
-// CV boundary points to be used by the discretization. The
-// discretization need not adopt the boundary points 100% faithfully;
-// for example, it may elide empty CVs or perform other transformations
-// for numeric fidelity or performance reasons.
-//
-// The cv_policy class is a value-like wrapper for actual
-// policies that derive from `cv_policy_base`. At present,
-// there are only three policies implemented, described below.
-// The intent is to provide more sophisticated policies in the
-// near future, specifically one based on the 'd-lambda' rule.
-//
-// cv_policy_explicit:
-// Simply use the provided locset.
-//
-// cv_policy_fixed_per_branch:
-// Use the same number of CVs for each branch.
-//
-// cv_policy_max_extent:
-// Use as many CVs as required to ensure that no CV has
-// a length longer than a given value.
-//
-// Except for the explicit policy, CV policies may also take various
-// flags (implemented as bitwise orable enums) to modify their
-// behaviour. In general, future CV policies may choose to ignore
-// flag values, but should respect them if their semantics are
-// relevant.
-//
-// cv_policy_flag::interior_forks:
-// Position CVs so as to include fork points, as opposed
-// to positioning them so that fork points are at the
-// boundaries of CVs.
-//
-// cv_policy_flag::single_root_cv:
-// Always treat the root branch as a single CV regardless.
-
-class cable_cell;
-
-struct cv_policy_base {
- virtual locset cv_boundary_points(const cable_cell& cell) const = 0;
- virtual std::unique_ptr<cv_policy_base> clone() const = 0;
- virtual ~cv_policy_base() {}
-};
-
-using cv_policy_base_ptr = std::unique_ptr<cv_policy_base>;
-
-struct cv_policy {
- cv_policy(const cv_policy_base& ref) { // implicit
- policy_ptr = ref.clone();
- }
-
- cv_policy(const cv_policy& other):
- policy_ptr(other.policy_ptr->clone()) {}
-
- cv_policy& operator=(const cv_policy& other) {
- policy_ptr = other.policy_ptr->clone();
- return *this;
- }
-
- cv_policy(cv_policy&&) = default;
- cv_policy& operator=(cv_policy&&) = default;
-
- locset cv_boundary_points(const cable_cell& cell) const {
- return policy_ptr->cv_boundary_points(cell);
- }
-
-private:
- cv_policy_base_ptr policy_ptr;
-};
-
-// Common flags for CV policies; bitwise composable.
-namespace cv_policy_flag {
- using value = unsigned;
- enum : unsigned {
- none = 0,
- interior_forks = 1<<0,
- single_root_cv = 1<<1
- };
-}
-
-struct cv_policy_explicit: cv_policy_base {
- explicit cv_policy_explicit(locset locs): locs_(std::move(locs)) {}
-
- cv_policy_base_ptr clone() const override {
- return cv_policy_base_ptr(new cv_policy_explicit(*this));
- }
-
- locset cv_boundary_points(const cable_cell&) const override {
- return locs_;
- }
-
-private:
- locset locs_;
-};
-
-struct cv_policy_max_extent: cv_policy_base {
- explicit cv_policy_max_extent(double max_extent, cv_policy_flag::value flags = cv_policy_flag::none):
- max_extent_(max_extent), flags_(flags) {}
-
- cv_policy_base_ptr clone() const override {
- return cv_policy_base_ptr(new cv_policy_max_extent(*this));
- }
-
- locset cv_boundary_points(const cable_cell&) const override;
-
-private:
- double max_extent_;
- cv_policy_flag::value flags_;
-};
-
-struct cv_policy_fixed_per_branch: cv_policy_base {
- explicit cv_policy_fixed_per_branch(unsigned cv_per_branch, cv_policy_flag::value flags = cv_policy_flag::none):
- cv_per_branch_(cv_per_branch), flags_(flags) {}
-
- cv_policy_base_ptr clone() const override {
- return cv_policy_base_ptr(new cv_policy_fixed_per_branch(*this));
- }
-
- locset cv_boundary_points(const cable_cell&) const override;
-
-private:
- unsigned cv_per_branch_;
- cv_policy_flag::value flags_;
-};
-
-struct cv_policy_every_sample: cv_policy_base {
- explicit cv_policy_every_sample(cv_policy_flag::value flags = cv_policy_flag::none):
- flags_(flags) {}
-
- cv_policy_base_ptr clone() const override {
- return cv_policy_base_ptr(new cv_policy_every_sample(*this));
- }
-
- locset cv_boundary_points(const cable_cell&) const override;
-
-private:
- cv_policy_flag::value flags_;
-};
-
-inline cv_policy default_cv_policy() {
- return cv_policy_fixed_per_branch(1);
-}
-
// Cable cell ion and electrical defaults.
// Parameters can be given as per-cell and global defaults via
diff --git a/arbor/include/arbor/cv_policy.hpp b/arbor/include/arbor/cv_policy.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f59fe418fc6a027daca96dbc3e975a2c4d66c45a
--- /dev/null
+++ b/arbor/include/arbor/cv_policy.hpp
@@ -0,0 +1,206 @@
+#pragma once
+
+#include <memory>
+#include <utility>
+
+#include <arbor/morph/region.hpp>
+#include <arbor/morph/locset.hpp>
+
+// FVM discretization policies/hints.
+//
+// CV polices, given a cable cell, provide a locset comprising CV boundary
+// points to be used by the discretization. The discretization need not adopt
+// the boundary points 100% faithfully; for example, it may elide empty CVs or
+// perform other transformations for numeric fidelity or performance reasons.
+//
+// The cv_policy class is a value-like wrapper for actual policies that derive
+// from `cv_policy_base`. At present, there are only a handful of policies
+// implemented, described below. The intent is to provide more sophisticated
+// policies in the near future, specifically one based on the 'd-lambda' rule.
+//
+// cv_policy_explicit:
+// Simply use the provided locset.
+//
+// cv_policy_single:
+// One CV for whole region.
+//
+// cv_policy_fixed_per_branch:
+// Use the same number of CVs for each branch.
+//
+// cv_policy_max_extent:
+// Use as many CVs as required to ensure that no CV has
+// a length longer than a given value.
+//
+// The policies above can be restricted to apply only to a given region of a
+// cell morphology. If a region is supplied, the CV policy is applied to the
+// completion of each connected component of the morphology within the region,
+// as if the subtree described by each component were itself a full cell
+// morphology. The boundary points of these completed components are always
+// included as boundary points provided by the policy.
+//
+// Except for the single and explicit policies, CV policies may also take
+// various flags (implemented as bitwise orable enums) to modify their
+// behaviour. In general, future CV policies may choose to ignore flag values,
+// but should respect them if their semantics are relevant.
+//
+// cv_policy_flag::interior_forks:
+// Position CVs so as to include fork points, as opposed to positioning
+// them so that fork points are at the boundaries of CVs.
+//
+// CV policy objects can be combined. For two CV polices A & B,
+//
+// A + B:
+// Use the CV boundary points from both A and B.
+//
+// A | B:
+// Use the CV boundary points from A except for on the region where
+// B is defined, where the boundary points from B are used.
+
+namespace arb {
+
+class cable_cell;
+
+struct cv_policy_base {
+ virtual locset cv_boundary_points(const cable_cell& cell) const = 0;
+ virtual region domain() const = 0;
+ virtual std::unique_ptr<cv_policy_base> clone() const = 0;
+ virtual ~cv_policy_base() {}
+};
+
+using cv_policy_base_ptr = std::unique_ptr<cv_policy_base>;
+
+struct cv_policy {
+ cv_policy(const cv_policy_base& ref) { // implicit
+ policy_ptr = ref.clone();
+ }
+
+ cv_policy(const cv_policy& other):
+ policy_ptr(other.policy_ptr->clone()) {}
+
+ cv_policy& operator=(const cv_policy& other) {
+ policy_ptr = other.policy_ptr->clone();
+ return *this;
+ }
+
+ cv_policy(cv_policy&&) = default;
+ cv_policy& operator=(cv_policy&&) = default;
+
+ locset cv_boundary_points(const cable_cell& cell) const {
+ return policy_ptr->cv_boundary_points(cell);
+ }
+
+ region domain() const {
+ return policy_ptr->domain();
+ }
+
+private:
+ cv_policy_base_ptr policy_ptr;
+};
+
+cv_policy operator+(const cv_policy&, const cv_policy&);
+cv_policy operator|(const cv_policy&, const cv_policy&);
+
+
+// Common flags for CV policies; bitwise composable.
+namespace cv_policy_flag {
+ using value = unsigned;
+ enum : unsigned {
+ none = 0,
+ interior_forks = 1<<0
+ };
+}
+
+struct cv_policy_explicit: cv_policy_base {
+ explicit cv_policy_explicit(locset locs, region domain = reg::all()):
+ locs_(std::move(locs)), domain_(std::move(domain)) {}
+
+ cv_policy_base_ptr clone() const override {
+ return cv_policy_base_ptr(new cv_policy_explicit(*this));
+ }
+
+ locset cv_boundary_points(const cable_cell&) const override;
+ region domain() const override { return domain_; }
+
+private:
+ locset locs_;
+ region domain_;
+};
+
+struct cv_policy_single: cv_policy_base {
+ explicit cv_policy_single(region domain = reg::all()):
+ domain_(domain) {}
+
+ cv_policy_base_ptr clone() const override {
+ return cv_policy_base_ptr(new cv_policy_single(*this));
+ }
+
+ locset cv_boundary_points(const cable_cell&) const override {
+ return ls::cboundary(domain_);
+ }
+ region domain() const override { return domain_; }
+
+private:
+ region domain_;
+};
+
+struct cv_policy_max_extent: cv_policy_base {
+ cv_policy_max_extent(double max_extent, region domain, cv_policy_flag::value flags = cv_policy_flag::none):
+ max_extent_(max_extent), domain_(std::move(domain)), flags_(flags) {}
+
+ explicit cv_policy_max_extent(double max_extent, cv_policy_flag::value flags = cv_policy_flag::none):
+ max_extent_(max_extent), domain_(reg::all()), flags_(flags) {}
+
+ cv_policy_base_ptr clone() const override {
+ return cv_policy_base_ptr(new cv_policy_max_extent(*this));
+ }
+
+ locset cv_boundary_points(const cable_cell&) const override;
+ region domain() const override { return domain_; }
+
+private:
+ double max_extent_;
+ region domain_;
+ cv_policy_flag::value flags_;
+};
+
+struct cv_policy_fixed_per_branch: cv_policy_base {
+ cv_policy_fixed_per_branch(unsigned cv_per_branch, region domain, cv_policy_flag::value flags = cv_policy_flag::none):
+ cv_per_branch_(cv_per_branch), domain_(std::move(domain)), flags_(flags) {}
+
+ explicit cv_policy_fixed_per_branch(unsigned cv_per_branch, cv_policy_flag::value flags = cv_policy_flag::none):
+ cv_per_branch_(cv_per_branch), domain_(reg::all()), flags_(flags) {}
+
+ cv_policy_base_ptr clone() const override {
+ return cv_policy_base_ptr(new cv_policy_fixed_per_branch(*this));
+ }
+
+ locset cv_boundary_points(const cable_cell&) const override;
+ region domain() const override { return domain_; }
+
+private:
+ unsigned cv_per_branch_;
+ region domain_;
+ cv_policy_flag::value flags_;
+};
+
+struct cv_policy_every_sample: cv_policy_base {
+ explicit cv_policy_every_sample(region domain = reg::all()):
+ domain_(std::move(domain)) {}
+
+ cv_policy_base_ptr clone() const override {
+ return cv_policy_base_ptr(new cv_policy_every_sample(*this));
+ }
+
+ locset cv_boundary_points(const cable_cell&) const override;
+ region domain() const override { return domain_; }
+
+private:
+ region domain_;
+ cv_policy_flag::value flags_;
+};
+
+inline cv_policy default_cv_policy() {
+ return cv_policy_fixed_per_branch(1);
+}
+
+} // namespace arb
diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp
index a815b67f05c48294eaa7506e765aaa7f31e5f437..d6e3772b5970a83443530d70f2ba4c943eedd326 100644
--- a/arbor/include/arbor/morph/locset.hpp
+++ b/arbor/include/arbor/morph/locset.hpp
@@ -139,9 +139,16 @@ locset nil();
// Most distal points of a region.
locset most_distal(region reg);
-// Most proximal point of a region.
+// Most proximal points of a region.
locset most_proximal(region reg);
+// Boundary points of a region.
+locset boundary(region reg);
+
+// Completed boundary points of a region.
+// (Boundary of completed components.)
+locset cboundary(region reg);
+
// Returns all locations in a locset that are also in the region.
locset restrict(locset ls, region reg);
@@ -152,10 +159,14 @@ locset uniform(region reg, unsigned left, unsigned right, uint64_t seed);
// Proportional location on every branch.
locset on_branches(double pos);
+// Support of a locset (x s.t. x in locset).
+locset support(locset);
+
} // namespace ls
// Union of two locsets.
locset join(locset, locset);
+
// Multiset sum of two locsets.
locset sum(locset, locset);
diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
index 7c00efed00577a27f46ac025d416e62a43adb965..208c9926f6941ddd7d5cc3e1355921b2132b3d34 100644
--- a/arbor/include/arbor/morph/morphology.hpp
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -61,18 +61,6 @@ public:
friend std::ostream& operator<<(std::ostream&, const morphology&);
};
-// Morphology utility functions.
-
-// Find the set of locations in an mlocation_list for which there
-// are no other locations that are more proximal in that list.
-mlocation_list minset(const morphology&, const mlocation_list&);
-
-// Find the set of locations in an mlocation_list for which there
-// are no other locations that are more distal in the list.
-mlocation_list maxset(const morphology&, const mlocation_list&);
-
-mlocation canonical(const morphology&, mlocation);
-
// Represent a (possibly empty or disconnected) region on a morphology.
//
// Wraps an mcable_list, and satisfies the additional constraint that
@@ -135,8 +123,50 @@ private:
mcable_list cables_;
};
+// Morphology utility functions.
+
+mlocation canonical(const morphology&, mlocation);
+
+// Find the set of locations in an mlocation_list for which there
+// are no other locations that are more proximal in that list.
+mlocation_list minset(const morphology&, const mlocation_list&);
+
+// Find the set of locations in an mlocation_list for which there
+// are no other locations that are more distal in the list.
+mlocation_list maxset(const morphology&, const mlocation_list&);
+
// Reduced representation of an extent, excluding zero-length cables
// that are covered by more proximal or non-zero-length cables.
mcable_list canonical(const morphology& m, const mextent& a);
+// Determine the components of an extent.
+//
+// Let T be the topological tree described by a morphology and C be the
+// cover, comprising the disjoint union of unit intervals, one per branch.
+//
+// Let π be the projection from C onto T.
+//
+// Locations in C are ordered by distality: (b1, x1) < (b2, x2) if branches b1
+// and b2 are the same and x1<x2, or else if b1 is a more proximal branch than
+// b2.
+//
+// Locations in T are ordered by distality: given points a and b in C,
+// π(a) < π(b) if a<b and π(a) is not equal to π(b).
+//
+// (NOTE: the notion of the cover may be extended in the future to include
+// a 'most proximal point' (-1, 1) which projects to the root of the tree,
+// and which is strictly more proximal than any other point in the cover.)
+//
+// Let two locations a,b in an extent X of C be directed-path-connected if
+// there is an order-preserving map p: [0, 1] -> C such that π∘p is a
+// path in T, with p(0) = a and p(1) = b.
+//
+// The components E_i of an extent X are subsets such that for all x and y
+// in E_i, there exists a location a with both a, x and a, y
+// directed-path-connected in X, and such that for all x in E_i and all y in
+// E_j, with i not equal to j, x and y are not directed-path-connected in X.
+
+std::vector<mextent> components(const morphology& m, const mextent&);
+
+
} // namespace arb
diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp
index 7fc8e87b256ebd9c0e507d0f3cc979f0c764eab8..f325cc20e6baba6641bbd3f8415f79012ada1fbb 100644
--- a/arbor/include/arbor/morph/primitives.hpp
+++ b/arbor/include/arbor/morph/primitives.hpp
@@ -88,6 +88,7 @@ bool test_invariants(const mlocation_list&);
mlocation_list sum(const mlocation_list&, const mlocation_list&);
mlocation_list join(const mlocation_list&, const mlocation_list&);
mlocation_list intersection(const mlocation_list&, const mlocation_list&);
+mlocation_list support(mlocation_list);
// Describe an unbranched cable in the morphology.
//
@@ -104,8 +105,12 @@ struct mcable {
double prox_pos; // ∈ [0,1]
double dist_pos; // ∈ [0,1]
- friend mlocation prox_loc(const mcable&);
- friend mlocation dist_loc(const mcable&);
+ friend mlocation prox_loc(const mcable& c) {
+ return {c.branch, c.prox_pos};
+ }
+ friend mlocation dist_loc(const mcable& c) {
+ return {c.branch, c.dist_pos};
+ }
// branch ≠npos, and 0 ≤ prox_pos ≤ dist_pos ≤ 1
friend bool test_invariants(const mcable&);
diff --git a/arbor/include/arbor/morph/region.hpp b/arbor/include/arbor/morph/region.hpp
index 19fbcbaaa2ceef56aebb3cff151a62b413e78e98..5e60d945a372d1763ba28feda44709ceb4d68871 100644
--- a/arbor/include/arbor/morph/region.hpp
+++ b/arbor/include/arbor/morph/region.hpp
@@ -158,8 +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 including all covers of included fork points.
+// (Pre-image of projection onto the topological tree.)
+region complete(region);
// Region associated with a name.
region named(std::string);
@@ -172,4 +173,10 @@ region join(region, region);
// Intersection of two regions.
region intersect(region, region);
+// Closed complement of a region.
+region complement(region);
+
+// (Closure of) set difference of two regions.
+region difference(region a, region b);
+
} // namespace arb
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index 0e4df9c32ed6dc358831e22656b541bae7f37101..393f9b025a35939fcd6ae1430b641e806eb5acfb 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -242,6 +242,80 @@ std::ostream& operator<<(std::ostream& o, const most_proximal_& x) {
return o << "(proximal \"" << x.reg << "\")";
}
+// Boundary points of a region.
+//
+// The boundary points of a region R are defined as the most proximal
+// and most distal locations in the components of R.
+
+struct boundary_: locset_tag {
+ explicit boundary_(region reg): reg(std::move(reg)) {}
+ region reg;
+};
+
+locset boundary(region reg) {
+ return locset(boundary_(std::move(reg)));
+};
+
+mlocation_list thingify_(const boundary_& n, const mprovider& p) {
+ std::vector<mextent> comps = components(p.morphology(), thingify(n.reg, p));
+
+ mlocation_list L;
+
+ for (const mextent& comp: comps) {
+ mlocation_list proximal_set, distal_set;
+
+ for (const mcable& c: comp) {
+ proximal_set.push_back({c.branch, c.prox_pos});
+ distal_set.push_back({c.branch, c.dist_pos});
+ }
+
+ L = sum(L, minset(p.morphology(), proximal_set));
+ L = sum(L, maxset(p.morphology(), distal_set));
+ }
+ return support(std::move(L));
+}
+
+std::ostream& operator<<(std::ostream& o, const boundary_& x) {
+ return o << "(boundary " << x.reg << ")";
+}
+
+// Completed boundary points of a region.
+//
+// The completed boundary is the boundary of the completion of
+// each component.
+
+struct cboundary_: locset_tag {
+ explicit cboundary_(region reg): reg(std::move(reg)) {}
+ region reg;
+};
+
+locset cboundary(region reg) {
+ return locset(cboundary_(std::move(reg)));
+};
+
+mlocation_list thingify_(const cboundary_& n, const mprovider& p) {
+ std::vector<mextent> comps = components(p.morphology(), thingify(n.reg, p));
+
+ mlocation_list L;
+
+ for (const mextent& comp: comps) {
+ mlocation_list proximal_set, distal_set;
+
+ mextent ccomp = thingify(reg::complete(comp), p);
+ for (const mcable& c: ccomp.cables()) {
+ proximal_set.push_back({c.branch, c.prox_pos});
+ distal_set.push_back({c.branch, c.dist_pos});
+ }
+
+ L = sum(L, minset(p.morphology(), proximal_set));
+ L = sum(L, maxset(p.morphology(), distal_set));
+ }
+ return support(std::move(L));
+}
+
+std::ostream& operator<<(std::ostream& o, const cboundary_& x) {
+ return o << "(cboundary " << x.reg << ")";
+}
// Uniform locset.
@@ -348,6 +422,25 @@ std::ostream& operator<<(std::ostream& o, const lsum& x) {
return o << "(sum " << x.lhs << " " << x.rhs << ")";
}
+// Support of point set.
+
+struct lsup_: locset_tag {
+ locset arg;
+ lsup_(locset arg): arg(std::move(arg)) {}
+};
+
+locset support(locset arg) {
+ return locset{lsup_{std::move(arg)}};
+}
+
+mlocation_list thingify_(const lsup_& P, const mprovider& p) {
+ return support(thingify(P.arg, p));
+};
+
+std::ostream& operator<<(std::ostream& o, const lsup_& x) {
+ return o << "(support " << x.arg << ")";
+}
+
// Restrict a locset on to a region: returns all locations in the locset that
// are also in the region.
diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp
index 4f131e6be09633d07e38309bcc06ff26dd7919e7..623eb7d9bc0d51c79fde9724796606f075c08dfc 100644
--- a/arbor/morph/morphology.cpp
+++ b/arbor/morph/morphology.cpp
@@ -447,5 +447,32 @@ mextent join(const mextent& a, const mextent& b) {
return m;
}
+std::vector<mextent> components(const morphology& m, const mextent& ex) {
+ std::unordered_map<mlocation, unsigned> component_index;
+ std::vector<mcable_list> component_cables;
+
+ for (mcable c: ex) {
+ mlocation head = canonical(m, prox_loc(c));
+
+ unsigned index;
+ if (component_index.count(head)) {
+ index = component_index.at(head);
+ }
+ else {
+ index = component_cables.size();
+ component_cables.push_back({});
+ }
+
+ component_cables[index].push_back(c);
+ component_index[dist_loc(c)] = index;
+ }
+
+ std::vector<mextent> components;
+ for (auto& cl: component_cables) {
+ components.emplace_back(std::move(cl));
+ }
+ return components;
+}
+
} // namespace arb
diff --git a/arbor/morph/primitives.cpp b/arbor/morph/primitives.cpp
index f98ddfcf5d0423b4f736f8fb64fce390d01372b9..40d9061b4ad430697f8236017c611095a9d265a7 100644
--- a/arbor/morph/primitives.cpp
+++ b/arbor/morph/primitives.cpp
@@ -7,6 +7,7 @@
#include "io/sepval.hpp"
#include "util/span.hpp"
#include "util/rangeutil.hpp"
+#include "util/unique.hpp"
namespace arb {
@@ -125,16 +126,13 @@ mlocation_list intersection(const mlocation_list& lhs, const mlocation_list& rhs
return L;
}
-bool test_invariants(const mcable& c) {
- return (0.<=c.prox_pos && c.prox_pos<=c.dist_pos && c.dist_pos<=1.) && c.branch!=mnpos;
-}
-
-mlocation prox_loc(const mcable& c) {
- return {c.branch, c.prox_pos};
+mlocation_list support(mlocation_list L) {
+ util::unique_in_place(L);
+ return L;
}
-mlocation dist_loc(const mcable& c) {
- return {c.branch, c.dist_pos};
+bool test_invariants(const mcable& c) {
+ return (0.<=c.prox_pos && c.prox_pos<=c.dist_pos && c.dist_pos<=1.) && c.branch!=mnpos;
}
bool test_invariants(const mcable_list& l) {
diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
index 90c30690c75d2f410f8a8961939bbcf71a2cf013..89bfce428cf09eaf90e274edbc60c544f99dafd2 100644
--- a/arbor/morph/region.cpp
+++ b/arbor/morph/region.cpp
@@ -558,7 +558,7 @@ struct super_ {
region reg;
};
-region super(region r) {
+region complete(region r) {
return region(super_{std::move(r)});
}
@@ -617,7 +617,7 @@ mextent thingify_(const super_& r, const mprovider& p) {
}
std::ostream& operator<<(std::ostream& o, const super_& r) {
- return o << "(super " << r.reg << ")";
+ return o << "(complete " << r.reg << ")";
}
@@ -654,11 +654,69 @@ std::ostream& operator<<(std::ostream& o, const reg_or& x) {
return o << "(join " << x.lhs << " " << x.rhs << ")";
}
+
+// Complement of a region.
+
+struct reg_not: region_tag {
+ region r;
+ explicit reg_not(region r): r(std::move(r)) {}
+};
+
+mextent thingify_(const reg_not& P, const mprovider& p) {
+ auto nb = p.morphology().num_branches();
+ mcable_list result;
+
+ mextent rex = thingify(P.r, p);
+ auto rex_i = rex.begin();
+
+ for (auto i: util::make_span(nb)) {
+ while (rex_i!=rex.end() && rex_i->branch<i) ++rex_i;
+
+ double x = 0;
+ while (rex_i!=rex.end() && rex_i->branch==i) {
+ double y = rex_i->prox_pos;
+ if (y>x) {
+ result.push_back(mcable{i, x, y});
+ }
+
+ x = rex_i->dist_pos;
+ ++rex_i;
+ }
+
+ if (x<1) {
+ result.push_back(mcable{i, x, 1});
+ }
+ }
+
+ return mextent(result);
+}
+
+std::ostream& operator<<(std::ostream& o, const reg_not& x) {
+ return o << "(complement " << x.r << ")";
+}
+
+
+// Closed set difference of two regions.
+
+struct reg_minus: region_tag {
+ region lhs;
+ region rhs;
+ reg_minus(region lhs, region rhs): lhs(std::move(lhs)), rhs(std::move(rhs)) {}
+};
+
+mextent thingify_(const reg_minus& P, const mprovider& p) {
+ return thingify(intersect(std::move(P.lhs), complement(std::move(P.rhs))), p);
+}
+
+std::ostream& operator<<(std::ostream& o, const reg_minus& x) {
+ return o << "(difference " << x.lhs << " " << x.rhs << ")";
+}
+
} // namespace reg
-// The intersect and join operations in the arb:: namespace with region so that
-// ADL allows for construction of expressions with regions without having
-// to namespace qualify the intersect/join.
+// The intersect, join, complement and difference operations are in the arb::
+// namespace with region so that ADL allows for construction of expressions
+// with regions without having to namespace qualify these operations.
region intersect(region l, region r) {
return region{reg::reg_and(std::move(l), std::move(r))};
@@ -668,6 +726,14 @@ region join(region l, region r) {
return region{reg::reg_or(std::move(l), std::move(r))};
}
+region complement(region r) {
+ return region{reg::reg_not(std::move(r))};
+}
+
+region difference(region l, region r) {
+ return region{reg::reg_minus(std::move(l), std::move(r))};
+}
+
region::region() {
*this = reg::nil();
}
diff --git a/arbor/util/rangeutil.hpp b/arbor/util/rangeutil.hpp
index 4e45d8c01e7c306c6cdb0bb497a44502a517b679..7e6a5fda30db450992cc4ed81f695fb4dd5dc9fa 100644
--- a/arbor/util/rangeutil.hpp
+++ b/arbor/util/rangeutil.hpp
@@ -399,6 +399,39 @@ range<Rev, Rev> reverse_view(Seq&& seq) {
return range<Rev, Rev>(Rev(strict.right), Rev(strict.left));
}
+// Left fold (accumulate) over sequence.
+//
+// Note that the order of arguments follows the application order;
+// schematically:
+//
+// foldl f a [] = a
+// foldl f a (b:bs) = foldl f (f a b) bs
+//
+// The binary operator f will be invoked once per element in the
+// sequence in turn, with the running accumulator as the first
+// argument. If the iterators for the sequence deference to a
+// mutable lvalue, then mutation of the value in the input sequence
+// is explicitly permitted.
+//
+// std::accumulate(begin, end, init, f) is equivalent to
+// util::foldl(f, init, util::make_range(begin, end)).
+
+template <typename Seq, typename Acc, typename BinOp>
+auto foldl(BinOp f, Acc a, Seq&& seq) {
+ using std::begin;
+ using std::end;
+
+ auto b = begin(seq);
+ auto e = end(seq);
+
+ while (b!=e) {
+ a = f(std::move(a), *b);
+ ++b;
+ }
+ return a;
+}
+
+
} // namespace util
} // namespace arb
diff --git a/doc/cpp_cable_cell.rst b/doc/cpp_cable_cell.rst
index 02f13b3ed3335ae4a2d29ea032630f711ad7b22c..f60e322d0bd1267e8c5be8c1199a292b3b0e1677 100644
--- a/doc/cpp_cable_cell.rst
+++ b/doc/cpp_cable_cell.rst
@@ -141,6 +141,10 @@ value should be taken from the cell or global parameter set.
Local areal capacitance of the cell membrane, in Farads per square metre.
+ .. cpp:member:: util::optional<cv_policy> discretization
+
+ Method by which CV boundaries are determined when the cell is discretized.
+ See :ref:`cv-policies`.
Default parameters for a cell are given by the :cpp:expr:`default_parameters`
field in the :cpp:type:`cable_cell` object. This is a value of type :cpp:type:`cable_cell_parameter_set`,
@@ -567,3 +571,117 @@ with which it is associated.
* Metadata: ``std::vector<cable_probe_point_info>``. Target metadata for each
associated target.
+
+
+.. _cv-policies:
+
+Discretization and CV policies
+------------------------------
+
+For the purpose of simulation, cable cells are decomposed into discrete
+subcomponents called *control volumes* (CVs), following the finite volume method
+terminology. Each control volume comprises a connected subset of the
+morphology. Each fork point in the morphology will be the responsibility of
+a single CV, and as a special case a zero-volume CV can be used to represent
+a single fork point in isolation.
+
+The CVs are uniquely determined by a set of *B* of ``mlocation`` boundary points.
+For each non-terminal point *h* in *B*, there is a CV comprising the points
+{*x*: *h* ≤ *x* and ¬∃ *y* ∈ *B* s.t *h* < *y* < *x*}, where < and ≤ refer to the
+geometrical partial order of locations on the morphology. A fork point is
+owned by a CV if and only if all of its corresponding representative locations
+are in the CV.
+
+The set of boundary points used by the simulator is determined by a *CV policy*.
+These are objects of type ``cv_policy``, which has the following
+public methods:
+
+.. cpp:class:: cv_policy
+
+ .. cpp:function:: locset cv_boundary_points(const cable_cell&) const
+
+ Return a locset describing the boundary points for CVs on the given cell.
+
+ .. cpp:function:: region domain() const
+
+ Give the subset of a cell morphology on which this policy has been declared,
+ as a morphological ``region`` expression.
+
+Specific CV policy objects are created by functions described below (strictly
+speaking, these are class constructors for classes are implicit converted to
+``cv_policy`` objects). These all take a ``region`` parameter that restrict the
+domain of applicability of that policy; this facility is useful for specifying
+differing discretizations on different parts of a cell morphology. When a CV
+policy is constrained in this manner, the boundary of the domain will always
+constitute part of the CV boundary point set.
+
+CV policies can be combined with ``+`` and ``|`` operators. For two policies
+*A* and *B*, *A*Â +Â *B* is a policy which gives boundary points from both *A*
+and *B*, while *A*Â |Â *B* is a policy which gives all the boundary points from
+*B* together with those from *A* which do not within the domain of *B*.
+The domain of *A*Â +Â *B* and *A*Â |Â *B* is the union of the domains of *A* and
+*B*.
+
+``cv_policy_single``
+^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ cv_policy_single(region domain = reg::all())
+
+Use one CV for the whole cell, or one for each connected component of the
+supplied domain.
+
+``cv_policy_explicit``
+^^^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ cv_policy_explicit(locset locs, region domain = reg::all())
+
+Use the points given by ``locs`` for CV boundaries, optionally restricted to the
+supplied domain.
+
+``cv_policy_every_sample``
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ cv_policy_every_sample(region domain = reg::all())
+
+Use every sample point in the morpholgy definition as a CV boundary, optionally
+restricted to the supplied domain. Each fork point in the domain is
+represented by a trivial CV.
+
+``cv_policy_fixed_per_branch``
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ cv_policy_fixed_per_branch(unsigned cv_per_branch, region domain, cv_policy_flag::value flags = cv_policy_flag::none);
+
+ cv_policy_fixed_per_branch(unsigned cv_per_branch, cv_policy_flag::value flags = cv_policy_flag::none):
+
+For each branch in each connected component of the domain (or the whole cell,
+if no domain is given), evenly distribute boundary points along the branch so
+as to produce exactly ``cv_per_branch`` CVs.
+
+By default, CVs will terminate at branch ends. If the flag
+``cv_policy_flag::interior_forks`` is given, fork points will be included in
+non-trivial, branched CVs and CVs covering terminal points in the morphology
+will be half-sized.
+
+
+``cv_policy_max_extent``
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ cv_policy_max_extent(double max_extent, region domain, cv_policy_flag::value flags = cv_policy_flag::none);
+
+ cv_policy_max_extent(double max_extent, cv_policy_flag::value flags = cv_policy_flag::none):
+
+As for ``cv_policy_fixed_per_branch``, save that the number of CVs on any
+given branch will be chosen to be the smallest number that ensures no
+CV will have an extent on the branch longer than ``max_extent`` micrometres.
+
diff --git a/python/morph_parse.cpp b/python/morph_parse.cpp
index 6c261029ced2b92d0715c0db9a5ed1bd37561a15..ba35b9a0c4a3bffc56a76be49745b28a0cdac8ab 100644
--- a/python/morph_parse.cpp
+++ b/python/morph_parse.cpp
@@ -196,7 +196,7 @@ std::unordered_multimap<std::string, evaluator> eval_map {
{"proximal_interval", make_call<arb::locset>(
[](arb::locset ls){return arb::reg::proximal_interval(std::move(ls), std::numeric_limits<double>::max());},
"'proximal_interval' with 1 argument: (start:locset)")},
- {"super", make_call<arb::region>(arb::reg::super,
+ {"complete", make_call<arb::region>(arb::reg::complete,
"'super' with 1 argment: (reg:region)")},
{"radius_lt",make_call<arb::region, double>(arb::reg::radius_lt,
"'radius_lt' with 2 arguments: (reg:region radius:real)")},
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 966ff9fdba5219813029e9d171c1b946a3436bbb..c13aa40bc82116bf3c43122b6cc00068e6acad18 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -120,6 +120,7 @@ set(unit_sources
test_merge_events.cpp
test_merge_view.cpp
test_morphology.cpp
+ test_morph_components.cpp
test_morph_embedding.cpp
test_morph_expr.cpp
test_morph_place.cpp
diff --git a/test/unit/common_morphologies.hpp b/test/unit/common_morphologies.hpp
index a3b63bf4ac720e7b3c9d662faea0867de0f2899a..de00aac7e494db36fb21f6fadf3da04c810a4929 100644
--- a/test/unit/common_morphologies.hpp
+++ b/test/unit/common_morphologies.hpp
@@ -20,29 +20,20 @@ inline std::vector<arb::msample> make_morph_samples(unsigned n) {
static const arb::morphology m_empty;
-// spherical root, one branch
-static const arb::morphology m_sph_b1{arb::sample_tree(make_morph_samples(1), {arb::mnpos}), true};
-
-// regular root, one branch
+// One branch.
static const arb::morphology m_reg_b1{arb::sample_tree(make_morph_samples(2), {arb::mnpos, 0u}), false};
-// spherical root, six branches:
-// branch 0 (spherical root) has child branches 1 and 2; branch 2 has child branches 3, 4 and 5.
-static const arb::morphology m_sph_b6{arb::sample_tree(make_morph_samples(8), {arb::mnpos, 0u, 1u, 0u, 3u, 4u, 4u, 4u}), true};
-
-// regular root, six branches
+// Six branches:
// branch 0 has child branches 1 and 2; branch 2 has child branches 3, 4 and 5.
static const arb::morphology m_reg_b6{arb::sample_tree(make_morph_samples(7), {arb::mnpos, 0u, 1u, 1u, 2u, 2u, 2u}), false};
-// regular root, six branches, mutiple top level branches.
+// Six branches, mutiple top level branches:
// branch 0 has child branches 1 and 2; branch 3 has child branches 4 and 5.
static const arb::morphology m_mlt_b6{arb::sample_tree(make_morph_samples(7), {arb::mnpos, 0u, 1u, 1u, 0u, 4u, 4u}), false};
static std::pair<const char*, arb::morphology> test_morphologies[] = {
{"m_empty", m_empty},
- {"m_sph_b1", m_sph_b1},
{"m_reg_b1", m_reg_b1},
- {"m_sph_b6", m_sph_b6},
{"m_reg_b6", m_reg_b6},
{"m_mlt_b6", m_mlt_b6}
};
diff --git a/test/unit/test_cv_geom.cpp b/test/unit/test_cv_geom.cpp
index f71bb9e3d7ed5e270ec1f73142df25ab304a850d..a007a86b7277905b06a4b9ed8abae7c01f7302eb 100644
--- a/test/unit/test_cv_geom.cpp
+++ b/test/unit/test_cv_geom.cpp
@@ -127,8 +127,8 @@ 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();
+ auto complete = [] (const arb::morphology& m, arb::mcable c) {
+ return thingify(arb::reg::complete(arb::region(c)), arb::mprovider(m)).cables();
};
for (auto& p: test_morphologies) {
@@ -160,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(super(m,c), cables));
+ EXPECT_TRUE(testing::seq_eq(complete(m,c), cables));
}
else {
ASSERT_EQ(1u, cables.size());
@@ -171,7 +171,7 @@ TEST(cv_geom, one_cv_per_branch) {
// Confirm parent CV is fork CV:
if (i>0) {
- auto fork_ext = super(m, {c.branch, 0});
+ auto fork_ext = complete(m, {c.branch, 0});
mcable_list pcables = util::assign_from(geom.cables(geom.cv_parent[i]));
ASSERT_TRUE(testing::cablelist_eq(fork_ext, pcables));
}
@@ -311,8 +311,8 @@ TEST(cv_geom, location_cv) {
return mextent(cl);
};
- auto super = [] (const arb::morphology& m, arb::mcable c) {
- return thingify(arb::reg::super(arb::region(c)), arb::mprovider(m)).cables();
+ auto complete = [] (const arb::morphology& m, arb::mcable c) {
+ return thingify(arb::reg::complete(arb::region(c)), arb::mprovider(m)).cables();
};
// Two CVs per branch, plus trivial CV at forks.
@@ -334,7 +334,7 @@ TEST(cv_geom, location_cv) {
ASSERT_TRUE(cable0.prox_pos==cable0.dist_pos);
mcable_list clist = util::assign_from(cables);
- ASSERT_TRUE(testing::cablelist_eq(super(m, cable0), clist));
+ ASSERT_TRUE(testing::cablelist_eq(complete(m, cable0), clist));
}
}
diff --git a/test/unit/test_cv_layout.cpp b/test/unit/test_cv_layout.cpp
index af731d0f0503dd2541ef6ad8ec886e8ffa67b304..8502318d9e66b6b03a3a8f7255baee08b4c21625 100644
--- a/test/unit/test_cv_layout.cpp
+++ b/test/unit/test_cv_layout.cpp
@@ -47,6 +47,10 @@ TEST(cv_layout, trivial) {
std::vector<cable_cell> cells;
unsigned n_cv = 0;
for (auto& p: test_morphologies) {
+ // Skip morpohologies with more than one root branch, becaue
+ // they are not 'connected', and will generate multiple CVs.
+ if (p.second.branch_children(mnpos).size()>1u) continue;
+
cells.emplace_back(p.second);
n_cv += !p.second.empty(); // one cv per non-empty cell
}
diff --git a/test/unit/test_cv_policy.cpp b/test/unit/test_cv_policy.cpp
index bfb9c43369264e15307c38b798bc57d8cf028b35..bba0fa1531736602163a966dfae4eb70365375ce 100644
--- a/test/unit/test_cv_policy.cpp
+++ b/test/unit/test_cv_policy.cpp
@@ -9,113 +9,107 @@
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/morph/mprovider.hpp>
+#include <arbor/morph/region.hpp>
#include "util/filter.hpp"
#include "util/rangeutil.hpp"
#include "util/span.hpp"
-#include "common.hpp"
-#include "unit_test_catalogue.hpp"
-#include "../common_cells.hpp"
+#include "common_morphologies.hpp"
+#include "morph_pred.hpp"
using namespace arb;
using util::make_span;
+using testing::locset_eq;
+using testing::region_eq;
+using testing::mlocationlist_eq;
+
+using namespace common_morphology;
namespace {
- std::vector<msample> make_samples(unsigned n) {
- std::vector<msample> ms;
- for (auto i: make_span(n)) ms.push_back({{0., 0., (double)i, 0.5}, 5});
- return ms;
+ template <typename... A>
+ locset as_locset(mlocation head, A... tail) {
+ return join(locset(head), locset(tail)...);
}
+}
- // Test morphologies for CV determination:
- // Samples points have radius 0.5, giving an initial branch length of 1.0
- // for morphologies with spherical roots.
+TEST(cv_policy, single) {
+ // In all instances, expect the boundary points to correspond to
+ // the extremal points of the completions of the components of the
+ // supplied region.
- const morphology m_empty;
+ cable_cell cell(m_mlt_b6);
+ for (region reg:
+ {reg::all(), reg::branch(2), reg::cable(3, 0.25, 1.),
+ join(reg::cable(1, 0.75, 1), reg::branch(3), reg::cable(2, 0, 0.5)),
+ join(reg::cable(2, 0, 0.5), reg::branch(3), reg::cable(4, 0, 0.5))})
+ {
+ locset expected = ls::cboundary(reg);
+ EXPECT_TRUE(locset_eq(cell.provider(), ls::cboundary(reg), cv_policy_single(reg).cv_boundary_points(cell)));
+ }
- // spherical root, one branch
- const morphology m_sph_b1{sample_tree(make_samples(1), {mnpos}), true};
+ EXPECT_TRUE(locset_eq(cell.provider(), cv_policy_single().cv_boundary_points(cell), cv_policy_single(reg::all()).cv_boundary_points(cell)));
+}
- // regular root, one branch
- const morphology m_reg_b1{sample_tree(make_samples(2), {mnpos, 0u}), false};
+TEST(cv_policy, explicit_policy) {
+ using L = mlocation;
+ locset lset = as_locset(L{0, 0}, L{0, 0.5}, L{0, 1.}, L{1, 0.5}, L{4, 0.2});
- // spherical root, six branches
- const morphology m_sph_b6{sample_tree(make_samples(8), {mnpos, 0u, 1u, 0u, 3u, 4u, 4u, 4u}), true};
+ cv_policy pol = cv_policy_explicit(lset);
+ for (auto& m: {m_reg_b6, m_mlt_b6}) {
+ cable_cell cell(m);
- // regular root, six branches
- const morphology m_reg_b6{sample_tree(make_samples(7), {mnpos, 0u, 1u, 1u, 2u, 2u, 2u}), false};
+ locset result = pol.cv_boundary_points(cell);
+ locset expected = join(ls::boundary(reg::all()), lset);
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, result));
+ }
- // regular root, six branches, mutiple top level branches.
- const morphology m_mlt_b6{sample_tree(make_samples(7), {mnpos, 0u, 1u, 1u, 0u, 4u, 4u}), false};
+ // With cables 1 and 2, expect to pick up (1, 0.5) from locset,
+ // and cable ends (1, 0), (1, 1), (2, 0), (2, 1), as the two
+ // cables constitute two components.
- template <typename... A>
- locset as_locset(mlocation head, A... tail) {
- return join(locset(head), locset(tail)...);
+ region b12 = join(reg::branch(1), reg::branch(2));
+ pol = cv_policy_explicit(lset, b12);
+ for (auto& m: {m_reg_b6, m_mlt_b6}) {
+ cable_cell cell(m);
+
+ locset result = pol.cv_boundary_points(cell);
+ locset expected = as_locset(L{1, 0}, L{1, 0.5}, L{1, 1}, L{2, 0}, L{2, 1});
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, result));
}
-}
-TEST(cv_policy, explicit_policy) {
- using L = mlocation;
- locset lset = as_locset(L{0, 0}, L{0, 0.5}, L{0, 1.}, L{1, 0.5}, L{4, 0.2});
+ // Taking the completion of the two cables, the boundary of the region
+ // will be (0, 1), (1, 1), (2, 1) for m_mlt_b6.
- cv_policy pol = cv_policy_explicit(lset);
- for (auto& m: {m_sph_b6, m_reg_b6, m_mlt_b6}) {
+ pol = cv_policy_explicit(lset, reg::complete(b12));
+ for (auto& m: {m_mlt_b6}) {
cable_cell cell(m);
locset result = pol.cv_boundary_points(cell);
- EXPECT_EQ(thingify(lset, cell.provider()), thingify(result, cell.provider()));
+ locset expected = as_locset(L{0, 1}, L{1, 0.5}, L{1, 1}, L{2, 1});
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, result));
}
}
TEST(cv_policy, empty_morphology) {
- // Any policy applied to an empty morphology should give an empty locset,
- // with the exception of cv_policy_explicit (this is still being debated).
+ // Any policy applied to an empty morphology should give an empty locset.
using namespace cv_policy_flag;
cv_policy policies[] = {
cv_policy_fixed_per_branch(3),
- cv_policy_fixed_per_branch(3, single_root_cv|interior_forks),
+ cv_policy_fixed_per_branch(3, interior_forks),
cv_policy_max_extent(0.234),
- cv_policy_max_extent(0.234, single_root_cv|interior_forks)
+ cv_policy_max_extent(0.234, interior_forks),
+ cv_policy_single(),
+ cv_policy_single(reg::all()),
+ cv_policy_explicit(ls::location(0, 0))
};
cable_cell cell(m_empty);
- auto empty_loclist = thingify(ls::nil(), cell.provider());
for (auto& pol: policies) {
- EXPECT_EQ(empty_loclist, thingify(pol.cv_boundary_points(cell), cell.provider()));
- }
-}
-
-TEST(cv_policy, single_root_cv) {
- // For policies that respect the single_root_cv flag, the boundary points should
- // be the same as if the flag were not provided, except that the points on branch 0
- // should only ever be (0, 0) and (0, 1) if the morphology is not empty.
-
- using namespace cv_policy_flag;
-
- std::pair<cv_policy, cv_policy> policy_pairs[] = {
- {cv_policy_fixed_per_branch(3), cv_policy_fixed_per_branch(3, single_root_cv)},
- {cv_policy_fixed_per_branch(3, interior_forks), cv_policy_fixed_per_branch(3, single_root_cv|interior_forks)},
- {cv_policy_max_extent(0.234), cv_policy_max_extent(0.234, single_root_cv)},
- {cv_policy_max_extent(0.234, interior_forks), cv_policy_max_extent(0.234, single_root_cv|interior_forks)}
- };
-
- for (auto& morph: {m_sph_b1, m_reg_b1, m_sph_b6, m_reg_b6, m_mlt_b6}) {
- cable_cell cell(morph);
-
- for (auto& polpair: policy_pairs) {
- mlocation_list p1 = thingify(polpair.first.cv_boundary_points(cell), cell.provider());
- mlocation_list p2 = thingify(polpair.second.cv_boundary_points(cell), cell.provider());
-
- auto p1_no_b0 = util::filter(p1, [](mlocation l) { return l.branch>0; });
- mlocation_list expected = {{0,0}, {0,1}};
- expected.insert(expected.end(), p1_no_b0.begin(), p1_no_b0.end());
-
- EXPECT_EQ(expected, p2);
- }
+ EXPECT_TRUE(locset_eq(cell.provider(), ls::nil(), pol.cv_boundary_points(cell)));
}
}
@@ -123,49 +117,83 @@ TEST(cv_policy, fixed_per_branch) {
using namespace cv_policy_flag;
using L = mlocation;
- // root branch only
- for (auto& morph: {m_sph_b1, m_reg_b1}) {
- cable_cell cell(morph);
+ // Root branch only:
+ {
+ cable_cell cell(m_reg_b1);
{
// boundary fork points
cv_policy pol = cv_policy_fixed_per_branch(4);
- locset points = pol.cv_boundary_points(cell);
locset expected = as_locset(L{0, 0}, L{0, 0.25}, L{0, 0.5}, L{0, 0.75}, L{0, 1});
- EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
}
{
// interior fork points
cv_policy pol = cv_policy_fixed_per_branch(4, interior_forks);
locset points = pol.cv_boundary_points(cell);
- locset expected = as_locset(L{0, 0.125}, L{0, 0.375}, L{0, 0.625}, L{0, 0.875});
- EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
+ locset expected = as_locset(L{0, 0}, L{0, 0.125}, L{0, 0.375}, L{0, 0.625}, L{0, 0.875}, L{0, 1});
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
}
}
- // spherical root, six branches and multiple top level branches cases:
- // expected points are the same.
- for (auto& morph: {m_sph_b6, m_mlt_b6}) {
- cable_cell cell(morph);
+ // Multiple top level branches:
+ // top level branches are 0 and 3, terminal branches are 1, 2, 4 and 5.
+ {
+ cable_cell cell(m_mlt_b6);
{
- // boundary fork points
+ // With boundary fork points:
cv_policy pol = cv_policy_fixed_per_branch(2);
- locset points = pol.cv_boundary_points(cell);
locset expected = as_locset(
L{0, 0}, L{0, 0.5}, L{0,1}, L{1, 0}, L{1, 0.5}, L{1,1}, L{2, 0}, L{2, 0.5}, L{2,1},
L{3, 0}, L{3, 0.5}, L{3,1}, L{4, 0}, L{4, 0.5}, L{4,1}, L{5, 0}, L{5, 0.5}, L{5,1}
);
- EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
}
{
- // interior fork points
+ // With interior fork points:
cv_policy pol = cv_policy_fixed_per_branch(2, interior_forks);
- locset points = pol.cv_boundary_points(cell);
locset expected = as_locset(
- L{0, 0.25}, L{0, 0.75}, L{1, 0.25}, L{1, 0.75}, L{2, 0.25}, L{2, 0.75},
- L{3, 0.25}, L{3, 0.75}, L{4, 0.25}, L{4, 0.75}, L{5, 0.25}, L{5, 0.75}
+ L{0, 0}, L{0, 0.25}, L{0, 0.75},
+ L{1, 0.25}, L{1, 0.75}, L{1, 1.0},
+ L{2, 0.25}, L{2, 0.75}, L{2, 1.0},
+ L{3, 0}, L{3, 0.25}, L{3, 0.75},
+ L{4, 0.25}, L{4, 0.75}, L{4, 1.0},
+ L{5, 0.25}, L{5, 0.75}, L{5, 1.0}
);
- EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
+ }
+ }
+
+ // Restrict to an incomplete subtree (distal half of branch 0 and all of branch 2)
+ // in m_mlt_b6 morphology.
+ {
+ cable_cell cell(m_mlt_b6);
+ region reg = mcable_list{{0, 0.5, 1.}, {2, 0., 1.}};
+ {
+ // With two per branch and fork points as boundaries, expect to see:
+ // (0, 0.5), (0, 0.75), (0, 1) on branch 0;
+ // (2, 0), (2, 0.5), (2, 1) on branch 2;
+ // (1, 0) on branch 1.
+ cv_policy pol = cv_policy_fixed_per_branch(2, reg);
+ locset expected = as_locset(
+ L{0, 0.5}, L{0, 0.75}, L{0, 1},
+ L{1, 0},
+ L{2, 0}, L{2, 0.5}, L{2, 1}
+ );
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
+ }
+ {
+ // With two per branch and interior forks, expect to see:
+ // (0, 0.5), (0, 0.625), (0, 0.0875) on branch 0;
+ // (2, 0.25), (2, 0.75), (2, 1) on branch 2;
+ // (1, 0) on branch 1.
+ cv_policy pol = cv_policy_fixed_per_branch(2, reg, interior_forks);
+ locset expected = as_locset(
+ L{0, 0.5}, L{0, 0.625}, L{0, 0.875},
+ L{1, 0},
+ L{2, 0.25}, L{2, 0.75}, L{2, 1}
+ );
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
}
}
}
@@ -174,28 +202,41 @@ TEST(cv_policy, max_extent) {
using namespace cv_policy_flag;
using L = mlocation;
- // root branch only
- for (auto& morph: {m_sph_b1, m_reg_b1}) {
- cable_cell cell(morph);
+ // Root branch only:
+ {
+ cable_cell cell(m_reg_b1);
ASSERT_EQ(1.0, cell.embedding().branch_length(0));
{
- // extent of 0.25 should give exact fp calculation, giving
+ // Extent of 0.25 should give exact fp calculation, giving
// 4 CVs on the root branch.
cv_policy pol = cv_policy_max_extent(0.25);
- locset points = pol.cv_boundary_points(cell);
locset expected = as_locset(L{0, 0}, L{0, 0.25}, L{0, 0.5}, L{0, 0.75}, L{0, 1});
- EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
+ }
+ {
+ // Same, but applied to cable (0, 0.25, 0.75) should give 2 Cvs.
+ cv_policy pol = cv_policy_max_extent(0.25, reg::cable(0, 0.25, 0.75));
+ locset expected = as_locset(L{0, 0.25}, L{0, 0.5}, L{0, 0.75});
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
+
}
{
+ // Interior forks:
cv_policy pol = cv_policy_max_extent(0.25, interior_forks);
- locset points = pol.cv_boundary_points(cell);
- locset expected = as_locset(L{0, 0.125}, L{0, 0.375}, L{0, 0.625}, L{0, 0.875});
- EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
+ locset expected = as_locset(L{0, 0}, L{0, 0.125}, L{0, 0.375}, L{0, 0.625}, L{0, 0.875}, L{0, 1});
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
+ }
+ {
+ // Interior forks but restricted to sub-cable.
+ cv_policy pol = cv_policy_max_extent(0.25, reg::cable(0, 0.25, 0.75), interior_forks);
+ locset expected = as_locset(L{0, 0.25}, L{0, 0.375}, L{0, 0.625}, L{0, 0.75});
+ EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
+
}
}
- // cell with varying branch lengths; extent not exact fraction.
+ // Cell with varying branch lengths; extent not exact fraction:
{
cable_cell cell(m_mlt_b6);
ASSERT_EQ(1.0, cell.embedding().branch_length(0));
@@ -206,7 +247,7 @@ TEST(cv_policy, max_extent) {
ASSERT_EQ(2.0, cell.embedding().branch_length(5));
{
- // max extent of 0.6 should give two CVs on branches of length 1,
+ // Max extent of 0.6 should give two CVs on branches of length 1,
// four CVs on branches of length 2, and seven CVs on the branch of length 4.
cv_policy pol = cv_policy_max_extent(0.6);
mlocation_list points = thingify(pol.cv_boundary_points(cell), cell.provider());
@@ -217,7 +258,7 @@ TEST(cv_policy, max_extent) {
{1, 0}, {1, 0.5}, {1, 1},
{2, 0}, {2, 0.25}, {2, 0.5}, {2, 0.75}, {2, 1}
};
- EXPECT_EQ(expected_b012, points_b012);
+ EXPECT_TRUE(mlocationlist_eq(expected_b012, points_b012));
mlocation_list points_b3 = util::assign_from(util::filter(points, [](mlocation l) { return l.branch==3; }));
EXPECT_EQ(8u, points_b3.size());
@@ -245,31 +286,89 @@ TEST(cv_policy, every_sample) {
{
cable_cell cell(m);
cv_policy pol = cv_policy_every_sample();
- mlocation_list points = thingify(pol.cv_boundary_points(cell), cell.provider());
- util::sort(points);
mlocation_list expected = {
- {0, 0}, {0, 0.25}, {0, 0.5}, {0, 0.75}, {0, 1.},
- {1, 0}, {1, 0.25}, {1, 0.5}, {1, 0.75}, {1, 1.},
- {2, 0}, {2, 0.25}, {2, 0.5}, {2, 0.75}, {2, 1.}
+ {0, 0}, {0, 0.25}, {0, 0.5}, {0, 0.75}, {0, 1},
+ {1, 0}, {1, 0.25}, {1, 0.5}, {1, 0.75}, {1, 1},
+ {2, 0}, {2, 0.25}, {2, 0.5}, {2, 0.75}, {2, 1}
};
- EXPECT_EQ(expected, points);
+ EXPECT_TRUE(locset_eq(cell.provider(), locset(expected), pol.cv_boundary_points(cell)));
}
-
- // Single root CV:
+ // Restricting to the two child branches (disconnected):
{
cable_cell cell(m);
- cv_policy pol = cv_policy_every_sample(single_root_cv);
- mlocation_list points = thingify(pol.cv_boundary_points(cell), cell.provider());
- util::sort(points);
+ region reg = join(reg::branch(1), reg::branch(2));
+ cv_policy pol = cv_policy_every_sample(reg);
+ // Get samples from branches 1 and 2, plus boundary points from completions of each
+ // branch, viz. (0, 1), (2, 0), (1, 1) from branch 1 and (0, 1), (1, 0), (2, 1) from
+ // branch 2.
mlocation_list expected = {
- {0, 0}, {0, 1.},
- {1, 0}, {1, 0.25}, {1, 0.5}, {1, 0.75}, {1, 1.},
- {2, 0}, {2, 0.25}, {2, 0.5}, {2, 0.75}, {2, 1.}
+ {0, 1},
+ {1, 0}, {1, 0.25}, {1, 0.5}, {1, 0.75}, {1, 1},
+ {2, 0}, {2, 0.25}, {2, 0.5}, {2, 0.75}, {2, 1}
};
- EXPECT_EQ(expected, points);
+ EXPECT_TRUE(locset_eq(cell.provider(), locset(expected), pol.cv_boundary_points(cell)));
+ }
+}
+
+TEST(cv_policy, domain) {
+ using namespace cv_policy_flag;
+
+ region reg1 = join(reg::branch(1), reg::cable(2, 0, 0.5));
+ region reg2 = join(reg::branch(1), reg::cable(2, 0.5, 1), reg::cable(4, 0, 1));
+
+ cable_cell cell(m_mlt_b6);
+
+ EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_single(reg1).domain()));
+ EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_fixed_per_branch(3, reg1).domain()));
+ EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_fixed_per_branch(3, reg1, interior_forks).domain()));
+ EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_max_extent(3, reg1).domain()));
+ EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_max_extent(3, reg1, interior_forks).domain()));
+ EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_every_sample(reg1).domain()));
+
+ EXPECT_TRUE(region_eq(cell.provider(), join(reg1, reg2), (cv_policy_single(reg1)+cv_policy_single(reg2)).domain()));
+ EXPECT_TRUE(region_eq(cell.provider(), join(reg1, reg2), (cv_policy_single(reg1)|cv_policy_single(reg2)).domain()));
+}
+
+TEST(cv_policy, combinators) {
+ auto unique_sum = [](auto... a) { return ls::support(sum(locset(a)...)); };
+
+ cable_cell cell(m_reg_b6);
+ auto eval_locset_eq = [&cell](const locset& a, const cv_policy& p) {
+ return locset_eq(cell.provider(), a, p.cv_boundary_points(cell));
+ };
+
+ {
+ mlocation_list locs1{{0, 0.5}, {1, 0.25}, {2, 1}};
+ mlocation_list locs2{{0, 0.75}, {1, 0.25}, {4, 0}};
+ locset all_bdy = ls::boundary(reg::all());
+
+ ASSERT_TRUE(eval_locset_eq(unique_sum(all_bdy, locs1), cv_policy_explicit(locs1)));
+ ASSERT_TRUE(eval_locset_eq(unique_sum(all_bdy, locs2), cv_policy_explicit(locs2)));
+
+ EXPECT_TRUE(eval_locset_eq(unique_sum(all_bdy, locs1, locs2), cv_policy_explicit(locs1)+cv_policy_explicit(locs2)));
+ EXPECT_TRUE(eval_locset_eq(unique_sum(all_bdy, locs2), cv_policy_explicit(locs1)|cv_policy_explicit(locs2)));
+ }
+
+ {
+ region reg12 = join(reg::branch(1), reg::branch(2));
+ region reg23 = join(reg::branch(2), reg::branch(3));
+
+ cv_policy pol12 = cv_policy_explicit(ls::on_branches(0.5), reg12);
+ cv_policy pol23 = cv_policy_explicit(ls::on_branches(0.75), reg23);
+
+ using L = mlocation;
+
+ ASSERT_TRUE(eval_locset_eq(unique_sum(ls::boundary(reg12), L{1, 0.5}, L{2, 0.5}), pol12));
+ ASSERT_TRUE(eval_locset_eq(unique_sum(ls::boundary(reg23), L{2, 0.75}, L{3, 0.75}), pol23));
+
+ EXPECT_TRUE(eval_locset_eq(unique_sum(ls::boundary(reg12), ls::boundary(reg23), L{1, 0.5}, L{2, 0.5}, L{2, 0.75}, L{3, 0.75}),
+ pol12+pol23));
+
+ EXPECT_TRUE(eval_locset_eq(unique_sum(ls::boundary(reg12), ls::boundary(reg23), L{1, 0.5}, L{2, 0.75}, L{3, 0.75}),
+ pol12|pol23));
}
}
diff --git a/test/unit/test_morph_components.cpp b/test/unit/test_morph_components.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..32207ffc1b18de6a2ae32a3a69a627c35633e188
--- /dev/null
+++ b/test/unit/test_morph_components.cpp
@@ -0,0 +1,81 @@
+#include <cmath>
+#include <vector>
+
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/mprovider.hpp>
+
+#include "util/rangeutil.hpp"
+
+#include "../test/gtest.h"
+#include "common_morphologies.hpp"
+
+using namespace arb;
+
+TEST(morph, subbranch_components) {
+ using namespace common_morphology;
+
+ const auto& m = m_reg_b1; // simple cable
+ mextent ex(mcable_list{mcable{0, 0.125, 0.25}, mcable{0, 0.5, 0.625}, mcable{0, 0.625, 0.75}, mcable{0, 0.875, 1.}});
+ auto comps = components(m, ex);
+
+ ASSERT_EQ(3u, comps.size());
+ ASSERT_FALSE(util::any_of(comps, [](const auto& x) { return x.empty(); }));
+
+ util::sort_by(comps, [](const mextent& x) { return prox_loc(x.front()); });
+
+ EXPECT_EQ((mcable_list{mcable{0, 0.125, 0.25}}), comps[0].cables());
+ EXPECT_EQ((mcable_list{mcable{0, 0.5, 0.75}}), comps[1].cables());
+ EXPECT_EQ((mcable_list{mcable{0, 0.875, 1.}}), comps[2].cables());
+}
+
+TEST(morph, subtree_components) {
+ using namespace common_morphology;
+
+ const auto& m = m_mlt_b6; // cell with two 'Y's (branches 0, 1, 2 and branches 3, 4, 5) meeting at root.
+
+ // Component semantics has that initial segments of branches from a common fork are _not_ regarded
+ // as connected; a final segment of a branch and an initial segment of a child branch _are_ regarded
+ // as connected.
+
+ std::pair<mcable_list, std::vector<mcable_list>> test_cases[] =
+ {
+ // Full cell gives two components (one for each 'Y').
+ {
+ {mcable{0, 0, 1}, mcable{1, 0, 1}, mcable{2, 0, 1}, mcable{3, 0, 1}, mcable{4, 0, 1}, mcable{5, 0, 1}},
+ {
+ {mcable{0, 0, 1}, mcable{1, 0, 1}, mcable{2, 0, 1}},
+ {mcable{3, 0, 1}, mcable{4, 0, 1}, mcable{5, 0, 1}}
+ }
+ },
+
+ // Siblings are separated.
+ {
+ {mcable{1, 0, 1}, mcable{2, 0, 1}, mcable{4, 0, 1}, mcable{5, 0, 1}},
+ {
+ {mcable{1, 0, 1}},
+ {mcable{2, 0, 1}},
+ {mcable{4, 0, 1}},
+ {mcable{5, 0, 1}}
+ }
+ },
+
+ // Parent-child branches are connected if they include the fork point.
+ {
+ {mcable{0, 0.5, 1}, mcable{2, 0, 0.5}, mcable{3, 0.5, 1}, mcable{5, 0.5, 1}},
+ {
+ {mcable{0, 0.5, 1}, mcable{2, 0, 0.5}},
+ {mcable{3, 0.5, 1}},
+ {mcable{5, 0.5, 1}}
+ }
+ },
+ };
+
+ for (auto tc: test_cases) {
+ auto comps = components(m, mextent(tc.first));
+ util::sort_by(comps, [](const mextent& x) { return prox_loc(x.front()); });
+
+ std::vector<mcable_list> result = util::assign_from(util::transform_view(comps, [](const auto& x) { return x.cables(); }));
+ EXPECT_EQ(tc.second, result);
+ }
+}
+
diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp
index 08e52d1e919ebc46b3f66ab4c966c8837ce63b76..40307aa0b8630b3266ce70afc84c1edd195e3535 100644
--- a/test/unit/test_morph_expr.cpp
+++ b/test/unit/test_morph_expr.cpp
@@ -34,8 +34,10 @@ 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 s = reg::complete(t3);
auto all = reg::all();
+ auto df = difference(t2, t1);
+ auto cm = complement(t3);
EXPECT_EQ(to_string(c1), "(cable 1 0 1)");
EXPECT_EQ(to_string(c2), "(cable 4 0.125 0.5)");
@@ -43,12 +45,14 @@ 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(s), "(complete (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(all), "(all)");
+ EXPECT_EQ(to_string(df), "(difference (tag 2) (tag 1))");
+ EXPECT_EQ(to_string(cm), "(complement (tag 3))");
}
TEST(region, invalid_mcable) {
@@ -64,6 +68,7 @@ TEST(locset, expr_repn) {
auto term = ls::terminal();
auto samp = ls::sample(42);
auto loc = ls::location(2, 0.5);
+ auto bdy = ls::boundary(reg::tagged(1));
EXPECT_EQ(to_string(root), "(root)");
EXPECT_EQ(to_string(term), "(terminal)");
@@ -72,6 +77,7 @@ TEST(locset, expr_repn) {
EXPECT_EQ(to_string(sum(root, term, samp, loc)), "(sum (sum (sum (root) (terminal)) (sample 42)) (location 2 0.5))");
EXPECT_EQ(to_string(samp), "(sample 42)");
EXPECT_EQ(to_string(loc), "(location 2 0.5)");
+ EXPECT_EQ(to_string(bdy), "(boundary (tag 1))");
}
TEST(locset, invalid_mlocation) {
@@ -348,6 +354,10 @@ TEST(region, thingify_simple_morphologies) {
auto t1 = reg::tagged(1);
auto t2 = reg::tagged(2);
auto all = reg::all();
+ auto q1 = reg::cable(0, 0.0, 0.25);
+ auto q2 = reg::cable(0, 0.25, 0.5);
+ auto q3 = reg::cable(0, 0.5, 0.75);
+ auto q4 = reg::cable(0, 0.75, 1);
// Concrete cable lists
cl h1_{{0, 0.0, 0.5}};
@@ -372,6 +382,11 @@ TEST(region, thingify_simple_morphologies) {
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}}));
+ EXPECT_TRUE(region_eq(mp, complement(join(q1, q3)), join(q2, q4)));
+ EXPECT_TRUE(region_eq(mp, complement(q2), join(q1, q3, q4)));
+ EXPECT_TRUE(region_eq(mp, difference(h1, q1), q2));
+ EXPECT_TRUE(region_eq(mp, difference(h1, q3), h1));
+
// Check round-trip of implicit region conversions.
// (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());
@@ -487,7 +502,7 @@ TEST(region, thingify_moderate_morphologies) {
using reg::branch;
using reg::all;
using reg::cable;
- using reg::super;
+ using reg::complete;
auto soma = tagged(1);
auto axon = tagged(2);
@@ -517,12 +532,12 @@ 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{}));
- 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(complete(apic), complete(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(complete(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}}));
+ EXPECT_TRUE(region_eq(mp, intersect(complete(axon), dend), cl{{0,0,0}, {3,0,0}}));
+ EXPECT_TRUE(region_eq(mp, intersect(complete(dend), axon), cl{{1,0,0}, {1,1,1}}));
// Test distal and proximal interavls
auto start0_ = location(0, 0 );
@@ -771,7 +786,7 @@ TEST(region, thingify_complex_morphologies) {
{
mprovider mp(m);
using reg::cable;
- using reg::super;
+ using reg::complete;
using ls::most_distal;
using ls::most_proximal;
@@ -797,8 +812,8 @@ TEST(region, thingify_complex_morphologies) {
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}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(complete(reg_f_)), mp), mlocation_list{{1,1}}));
+ EXPECT_TRUE(mlocationlist_eq(thingify(most_proximal(complete(reg_g_)), mp), mlocation_list{{4,1}}));
}
}
{
diff --git a/test/unit/test_range.cpp b/test/unit/test_range.cpp
index ed48a7923c35d618bef62f6e4cee2a3c99e0180b..ff1e22765734a64b0b5aaca867cd2ff0dd0259f9 100644
--- a/test/unit/test_range.cpp
+++ b/test/unit/test_range.cpp
@@ -735,3 +735,33 @@ TEST(range, reverse) {
EXPECT_EQ("olleh"s, rev);
}
+
+TEST(range, foldl) {
+ // Check invocation order:
+ int xs[] = {1, 2, 3, 4};
+ int result = util::foldl(
+ [](int l, int r) { return 10*l+r; },
+ 0,
+ xs);
+
+ EXPECT_EQ(1234, result);
+
+ // Check mutability permission:
+ result = util::foldl(
+ [](int l, int& r) { return ++r, 10*l+r; },
+ 0,
+ xs);
+ EXPECT_EQ(2345, result);
+ EXPECT_EQ(5,xs[3]);
+
+ // Check works with move-only values:
+ nocopy<int> ns[] = {1, 2, 3, 4};
+ auto plus = [](nocopy<int> a, nocopy<int> b) { return nocopy<int>(10*a.value + b.value); };
+ auto nc_result = util::foldl(
+ [&plus](auto a, auto& b) { return plus(std::move(a), std::move(b)); },
+ nocopy<int>(0),
+ ns);
+
+ EXPECT_EQ(1234, nc_result.value);
+ EXPECT_EQ(0, ns[3].value);
+}