diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp
index a60bcfb2a2facd822240752a7cb242d0629dd62b..9a13e1a93ac17a32c69b9387a834653698f9e904 100644
--- a/arbor/cable_cell_param.cpp
+++ b/arbor/cable_cell_param.cpp
@@ -1,8 +1,13 @@
#include <cfloat>
#include <cmath>
+#include <numeric>
+#include <vector>
+#include <arbor/cable_cell.hpp>
#include <arbor/cable_cell_param.hpp>
+#include <arbor/morph/locset.hpp>
+#include "morph/em_morphology.hpp"
#include "util/maputil.hpp"
namespace arb {
@@ -64,4 +69,74 @@ cable_cell_local_parameter_set neuron_parameter_defaults = {
0.01
};
+// Discretization policy implementations:
+
+locset cv_policy_max_extent::cv_boundary_points(const cable_cell& cell) const {
+ const auto& emorph = *cell.morphology();
+ const unsigned nbranch = emorph.num_branches();
+ 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(emorph.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;
+ }
+
+ return std::accumulate(points.begin(), points.end(), ls::nil(), [](auto& l, auto& p) { return join(l, ls::location(p)); });
+}
+
+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;
+ }
+
+ return std::accumulate(points.begin(), points.end(), ls::nil(), [](auto& l, auto& p) { return join(l, ls::location(p)); });
+}
+
} // namespace arb
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index 8ef752abae28868b00d2feecb11ebd065d61e77e..0de0dff56e55a77b8136ae36c4b6ba15835dddcb 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -2,8 +2,10 @@
#include <arbor/arbexcept.hpp>
#include <arbor/mechcat.hpp>
+#include <arbor/morph/locset.hpp>
#include <arbor/util/optional.hpp>
+#include <memory>
#include <unordered_map>
#include <string>
@@ -39,7 +41,6 @@ struct threshold_detector {
// Tag type for dispatching cable_cell::place() calls that add gap junction sites.
struct gap_junction_site {};
-
// Mechanism description, viz. mechanism name and
// (non-global) parameter settings. Used to assign
// density and point mechanisms to segments and
@@ -103,6 +104,134 @@ private:
std::unordered_map<std::string, double> param_;
};
+// 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(cv_policy&&) = default;
+
+ cv_policy(const cv_policy& other):
+ policy_ptr(other.policy_ptr->clone()) {}
+
+ cv_policy& operator=(const 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_;
+};
+
+inline cv_policy default_cv_policy() {
+ return cv_policy_fixed_per_branch(1);
+}
+
// Cable cell ion and electrical defaults.
//
// Parameters can be overridden with `cable_cell_local_parameter_set`
@@ -127,11 +256,18 @@ struct cable_cell_local_parameter_set {
struct cable_cell_parameter_set: public cable_cell_local_parameter_set {
std::unordered_map<std::string, mechanism_desc> reversal_potential_method;
+ cv_policy discretization = default_cv_policy();
// We'll need something like this until C++17, for sane initialization syntax.
cable_cell_parameter_set() = default;
- cable_cell_parameter_set(cable_cell_local_parameter_set p, std::unordered_map<std::string, mechanism_desc> m = {}):
- cable_cell_local_parameter_set(std::move(p)), reversal_potential_method(std::move(m))
+ cable_cell_parameter_set(
+ cable_cell_local_parameter_set p,
+ std::unordered_map<std::string, mechanism_desc> m = {},
+ cv_policy d = default_cv_policy()
+ ):
+ cable_cell_local_parameter_set(std::move(p)),
+ reversal_potential_method(std::move(m)),
+ discretization(std::move(d))
{}
};
diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp
index 63e0e6923e1fc4ecf772cb6c7de0c53f43e1eb81..a9a64cbb663122577d07fe19ff076eb115efac26 100644
--- a/arbor/include/arbor/morph/locset.hpp
+++ b/arbor/include/arbor/morph/locset.hpp
@@ -75,6 +75,14 @@ public:
return sum(sum(std::move(l), std::move(r)), std::move(args)...);
}
+ // The union of two location sets.
+ friend locset join(locset, locset);
+
+ template <typename ...Args>
+ friend locset join(locset l, locset r, Args... args) {
+ return join(join(std::move(l), std::move(r)), std::move(args)...);
+ }
+
private:
struct interface {
virtual ~interface() {}
diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
index 16d41ec763479797434a39f533e6b29ec8afcc68..4e90d30a483cd03ebf53da76450622f02403cb64 100644
--- a/arbor/include/arbor/morph/morphology.hpp
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -23,6 +23,9 @@ public:
morphology(sample_tree m);
morphology();
+ // Empty/default-constructed morphology?
+ bool empty() const;
+
// Whether the root of the morphology is spherical.
bool spherical_root() const;
diff --git a/arbor/morph/em_morphology.cpp b/arbor/morph/em_morphology.cpp
index f00865e3c59259730e020fff7370d715451b0b5e..8e49c3a9d42401d95038d5fbbeea28edeb23237a 100644
--- a/arbor/morph/em_morphology.cpp
+++ b/arbor/morph/em_morphology.cpp
@@ -40,6 +40,9 @@ em_morphology::em_morphology(const morphology& m):
auto idx = util::make_range(morph_.branch_indexes(i));
branch_lengths_.push_back(dist2root_[idx.back()]- dist2root_[idx.front()]);
}
+ if (morph_.spherical_root()) {
+ branch_lengths_[0] = samples[0].loc.radius*2;
+ }
// Cache the sample locations.
// Iterate backwards over branches distal to root, so that the parent branch at
@@ -132,13 +135,17 @@ mlocation_list em_morphology::cover(mlocation loc, bool include_loc) const {
return L;
}
-mlocation em_morphology::canonicalize(mlocation loc) const {
+void em_morphology::assert_valid_location(mlocation loc) const {
if (!test_invariants(loc)) {
throw morphology_error(util::pprintf("Invalid location {}", loc));
}
if (loc.branch>=morph_.num_branches()) {
throw morphology_error(util::pprintf("Location {} does not exist in morpology", loc));
}
+}
+
+mlocation em_morphology::canonicalize(mlocation loc) const {
+ assert_valid_location(loc);
// Test if location is at the start of a branch.
if (loc.pos==0.) {
diff --git a/arbor/morph/em_morphology.hpp b/arbor/morph/em_morphology.hpp
index 1e47dd59b5dd4c32ec33fd1e5bfc7f11ba315063..1feb88cda374c17985a269856e2857bce91591d4 100644
--- a/arbor/morph/em_morphology.hpp
+++ b/arbor/morph/em_morphology.hpp
@@ -24,11 +24,24 @@ public:
const morphology& morph() const;
+ // Convenience methods for morphology access
+ // that are forwarded directly to the morphology object:
+
+ auto empty() const { return morph_.empty(); }
+ auto spherical_root() const { return morph_.spherical_root(); }
+ auto num_branches() const { return morph_.num_branches(); }
+ auto num_samples() const { return morph_.num_samples(); }
+ auto branch_parent(msize_t b) const { return morph_.branch_parent(b); }
+ auto branch_children(msize_t b) const { return morph_.branch_children(b); }
+
+ // Access to computed and cached data:
+
mlocation_list terminals() const;
mlocation root() const;
mlocation sample2loc(msize_t sid) const;
+ void assert_valid_location(mlocation) const;
mlocation canonicalize(mlocation) const;
// Find all locations on the morphology that share the same canonoical
@@ -37,6 +50,8 @@ public:
mlocation_list cover(mlocation, bool include_loc=true) const;
mlocation_list minset(const mlocation_list&) const;
+
+ double branch_length(msize_t bid) const { return branch_lengths_.at(bid); }
};
} // namespace arb
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index 55c0752ccf366594db1c6fb211e91a4d85fe72f9..c21bb3d69ded5d2b146f525d1217ebcb99c3a50f 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -125,8 +125,8 @@ locset location(mlocation loc) {
}
mlocation_list thingify_(const location_& x, const em_morphology& m) {
- // canonicalize will throw if the location is not present.
- return {m.canonicalize(x.loc)};
+ m.assert_valid_location(x.loc);
+ return {x.loc};
}
std::ostream& operator<<(std::ostream& o, const location_& x) {
diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp
index ee78c8532ef3a2121aeed365ee0fe1b9cfb53999..34253041750bba913d66d0acf25bd4adc4401b94 100644
--- a/arbor/morph/morphology.cpp
+++ b/arbor/morph/morphology.cpp
@@ -79,7 +79,7 @@ std::vector<mbranch> branches_from_parent_index(const std::vector<msize_t>& pare
}
// Returns false if one of the root's children has the same tag as the root.
-bool root_sample_has_same_tag_as_child(const sample_tree& st) {
+bool root_sample_tag_differs_from_children(const sample_tree& st) {
if (st.empty()) return false;
auto& P = st.parents();
auto& S = st.samples();
@@ -128,7 +128,7 @@ morphology_impl::morphology_impl(sample_tree m, bool use_spherical_root):
morphology_impl::morphology_impl(sample_tree m):
samples_(std::move(m)),
- spherical_root_(impl::root_sample_has_same_tag_as_child(samples_))
+ spherical_root_(impl::root_sample_tag_differs_from_children(samples_))
{
init();
}
@@ -186,6 +186,10 @@ morphology::morphology():
morphology(sample_tree())
{}
+bool morphology::empty() const {
+ return impl_->branches_.empty();
+}
+
// The parent branch of branch b.
msize_t morphology::branch_parent(msize_t b) const {
return impl_->branch_parents_[b];
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 0cd9afad10ba24dd970ea26fb36726f9164a8bc2..910764d41e8830a88b1a3b6792dee921f518e82e 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -77,6 +77,7 @@ set(unit_sources
test_cable_cell.cpp
test_compartments.cpp
test_counter.cpp
+ test_cv_policy.cpp
test_cycle.cpp
test_domain_decomposition.cpp
test_dry_run_context.cpp
diff --git a/test/unit/test_cv_policy.cpp b/test/unit/test_cv_policy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a7a099682dcf54b38f57f660a5d87c343c590b43
--- /dev/null
+++ b/test/unit/test_cv_policy.cpp
@@ -0,0 +1,235 @@
+#include <iterator>
+#include <numeric>
+#include <utility>
+#include <vector>
+
+#include <arbor/util/optional.hpp>
+#include <arbor/cable_cell.hpp>
+#include <arbor/cable_cell_param.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/locset.hpp>
+
+#include "morph/em_morphology.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"
+
+using namespace arb;
+using util::make_span;
+
+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;
+ }
+
+ // Test morphologies for CV determination:
+ // Samples points have radius 0.5, giving an initial branch length of 1.0
+ // for morphologies with spherical roots.
+
+ const morphology m_empty;
+
+ // spherical root, one branch
+ const morphology m_sph_b1{sample_tree(make_samples(1), {mnpos}), true};
+
+ // regular root, one branch
+ const morphology m_reg_b1{sample_tree(make_samples(2), {mnpos, 0u}), false};
+
+ // spherical root, six branches
+ const morphology m_sph_b6{sample_tree(make_samples(8), {mnpos, 0u, 1u, 0u, 3u, 4u, 4u, 4u}), true};
+
+ // regular root, six branches
+ const morphology m_reg_b6{sample_tree(make_samples(7), {mnpos, 0u, 1u, 1u, 2u, 2u, 2u}), false};
+
+ // 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};
+
+ template <typename... A>
+ locset as_locset(mlocation head, A... tail) {
+ return join(ls::location(head), ls::location(tail)...);
+ }
+
+ template <typename Seq>
+ locset as_locset(const Seq& seq) {
+ using std::begin;
+ using std::end;
+ return std::accumulate(begin(seq), end(seq), ls::nil(),
+ [](locset ls, const mlocation& p) { return join(std::move(ls), ls::location(p)); });
+ }
+}
+
+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});
+
+ cv_policy pol = cv_policy_explicit(lset);
+ for (auto& m: {m_sph_b6, m_reg_b6, m_mlt_b6}) {
+ cable_cell cell = make_cable_cell(m);
+
+ locset result = pol.cv_boundary_points(cell);
+ EXPECT_EQ(thingify(lset, *cell.morphology()), thingify(result, *cell.morphology()));
+ }
+}
+
+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).
+
+ 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_max_extent(0.234),
+ cv_policy_max_extent(0.234, single_root_cv|interior_forks)
+ };
+
+ cable_cell cell = make_cable_cell(m_empty);
+ auto empty_loclist = thingify(ls::nil(), *cell.morphology());
+
+ for (auto& pol: policies) {
+ EXPECT_EQ(empty_loclist, thingify(pol.cv_boundary_points(cell), *cell.morphology()));
+ }
+}
+
+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 = make_cable_cell(morph);
+
+ for (auto& polpair: policy_pairs) {
+ mlocation_list p1 = thingify(polpair.first.cv_boundary_points(cell), *cell.morphology());
+ mlocation_list p2 = thingify(polpair.second.cv_boundary_points(cell), *cell.morphology());
+
+ 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);
+ }
+ }
+}
+
+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 = make_cable_cell(morph);
+ {
+ // 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.morphology()), thingify(points, *cell.morphology()));
+ }
+ {
+ // 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.morphology()), thingify(points, *cell.morphology()));
+ }
+ }
+
+ // 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 = make_cable_cell(morph);
+
+ {
+ // 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.morphology()), thingify(points, *cell.morphology()));
+ }
+ {
+ // 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}
+ );
+ EXPECT_EQ(thingify(expected, *cell.morphology()), thingify(points, *cell.morphology()));
+ }
+ }
+}
+
+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 = make_cable_cell(morph);
+ ASSERT_EQ(1.0, cell.morphology()->branch_length(0));
+
+ {
+ // 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.morphology()), thingify(points, *cell.morphology()));
+ }
+ {
+ 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.morphology()), thingify(points, *cell.morphology()));
+ }
+ }
+
+ // cell with varying branch lengths; extent not exact fraction.
+ {
+ cable_cell cell = make_cable_cell(m_mlt_b6);
+ ASSERT_EQ(1.0, cell.morphology()->branch_length(0));
+ ASSERT_EQ(1.0, cell.morphology()->branch_length(1));
+ ASSERT_EQ(2.0, cell.morphology()->branch_length(2));
+ ASSERT_EQ(4.0, cell.morphology()->branch_length(3));
+ ASSERT_EQ(1.0, cell.morphology()->branch_length(4));
+ ASSERT_EQ(2.0, cell.morphology()->branch_length(5));
+
+ {
+ // 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.morphology());
+
+ mlocation_list points_b012 = util::assign_from(util::filter(points, [](mlocation l) { return l.branch<3; }));
+ mlocation_list expected_b012 = {
+ {0, 0}, {0, 0.5}, {0, 1},
+ {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);
+
+ mlocation_list points_b3 = util::assign_from(util::filter(points, [](mlocation l) { return l.branch==3; }));
+ EXPECT_EQ(8u, points_b3.size());
+ }
+ }
+}
+
diff --git a/test/unit/test_em_morphology.cpp b/test/unit/test_em_morphology.cpp
index 39528837f92c891ecc1f141e6ff71de20dcac907..36cbf05f5b1e0a686445c70db71cd47d3097f70a 100644
--- a/test/unit/test_em_morphology.cpp
+++ b/test/unit/test_em_morphology.cpp
@@ -54,6 +54,8 @@ TEST(em_morphology, cache) {
EXPECT_EQ(em.root(), (loc{0,0}));
EXPECT_EQ(em.terminals(), (arb::mlocation_list{{0,1}}));
+
+ EXPECT_EQ(10., em.branch_length(0));
}
// Eight samples
@@ -91,6 +93,13 @@ TEST(em_morphology, cache) {
EXPECT_EQ(em.root(), (loc{0,0}));
EXPECT_EQ(em.terminals(), (arb::mlocation_list{{1,1}, {3,1}, {4,1}}));
+
+ ASSERT_EQ(5u, em.morph().num_branches());
+ EXPECT_EQ(20., em.branch_length(0));
+ EXPECT_EQ(90., em.branch_length(1));
+ EXPECT_EQ(90., em.branch_length(2));
+ EXPECT_EQ(100., em.branch_length(3));
+ EXPECT_EQ(200., em.branch_length(4));
}
{ // No Spherical root
pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
@@ -119,6 +128,12 @@ TEST(em_morphology, cache) {
EXPECT_EQ(em.root(), (loc{0,0}));
EXPECT_EQ(em.terminals(), (arb::mlocation_list{{0,1}, {2,1}, {3,1}}));
+
+ ASSERT_EQ(4u, em.morph().num_branches());
+ EXPECT_EQ(100., em.branch_length(0));
+ EXPECT_EQ(100., em.branch_length(1));
+ EXPECT_EQ(100., em.branch_length(2));
+ EXPECT_EQ(200., em.branch_length(3));
}
}
@@ -348,10 +363,10 @@ TEST(locset, thingify) {
EXPECT_EQ(thingify(midb2, em), (ll{{2,0.5}}));
EXPECT_EQ(thingify(midb1, em), (ll{{1,0.5}}));
EXPECT_EQ(thingify(begb0, em), (ll{{0,0}}));
- EXPECT_EQ(thingify(begb1, em), (ll{{0,1}}));
- EXPECT_EQ(thingify(begb2, em), (ll{{0,1}}));
- EXPECT_EQ(thingify(begb3, em), (ll{{2,1}}));
- EXPECT_EQ(thingify(begb4, em), (ll{{2,1}}));
+ EXPECT_EQ(thingify(begb1, em), (ll{{1,0}}));
+ EXPECT_EQ(thingify(begb2, em), (ll{{2,0}}));
+ EXPECT_EQ(thingify(begb3, em), (ll{{3,0}}));
+ EXPECT_EQ(thingify(begb4, em), (ll{{4,0}}));
}
{
arb::em_morphology em(arb::morphology(sm, false));
@@ -362,9 +377,9 @@ TEST(locset, thingify) {
EXPECT_EQ(thingify(midb2, em), (ll{{2,0.5}}));
EXPECT_EQ(thingify(midb1, em), (ll{{1,0.5}}));
EXPECT_EQ(thingify(begb0, em), (ll{{0,0}}));
- EXPECT_EQ(thingify(begb1, em), (ll{{0,0}}));
- EXPECT_EQ(thingify(begb2, em), (ll{{1,1}}));
- EXPECT_EQ(thingify(begb3, em), (ll{{1,1}}));
+ EXPECT_EQ(thingify(begb1, em), (ll{{1,0}}));
+ EXPECT_EQ(thingify(begb2, em), (ll{{2,0}}));
+ EXPECT_EQ(thingify(begb3, em), (ll{{3,0}}));
// In the absence of a spherical root, there is no branch 4.
EXPECT_THROW(thingify(begb4, em), arb::morphology_error);
}