CV policies with unit tests. (#878)

* Add new `cv_policy` class that wraps classes derived from `cv_policy_base`, which in turn provide a facility for producing CV boundaries on morphology trees from a `cable_cell`.
* Implement fixed-CVs-per-branch and maximum-CV-length policies, with flags that allow for single-CV oots and for CVs that terminate at fork points or which are centred on fork points.
* Expose locset `join` expressions.
* Add `empty()` test method on `morphology`.
* Expose branch length on `em_morphology` via member function.
* Account for the length of a spherical root branch in `em_morphology`.
* Test for validity but do not canonicalize locations when thingifying.
* Rename `root_sample_has_same_tag_as_child` to `root_sample_tag_differs_from_children` to better reflect sense of test.
* Add unit tests for CV policies, branch_length method.

arbor/cable_cell_param.cpp

 #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

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))
     {}
 };
 

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() {}

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;
 

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.) {

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

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) {

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];

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

test/unit/test_cv_policy.cpp

+#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());
+        }
+    }
+}
+

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);
     }