diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index b4c1b99479d971ce6f9c91f04aa702d8bf287ee3..462ce9f6b48a2f5e2c7fc3a596199ad27d757ad2 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -28,6 +28,7 @@ set(arbor_sources
mechinfo.cpp
memory/gpu_wrappers.cpp
memory/util.cpp
+ morph/cv_data.cpp
morph/embed_pwlin.cpp
morph/label_dict.cpp
morph/locset.cpp
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index d926192cc5bed96a5f9c8c37db74891d530b2c9c..db52f29b0b1bb39517883a0069b720322cef48b5 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -18,6 +18,7 @@
#include "util/meta.hpp"
#include "util/partition.hpp"
#include "util/piecewise.hpp"
+#include "util/pw_over_cable.hpp"
#include "util/rangeutil.hpp"
#include "util/transform.hpp"
#include "util/unique.hpp"
@@ -30,19 +31,13 @@ using util::count_along;
using util::make_span;
using util::pw_elements;
using util::pw_element;
+using util::pw_over_cable;
using util::sort;
using util::sort_by;
using util::stable_sort_by;
using util::value_by_key;
namespace {
-struct get_value {
- template <typename X>
- double operator()(const X& x) const { return x.value; }
-
- double operator()(double x) const { return x; }
-};
-
template <typename V>
std::optional<V> operator|(const std::optional<V>& a, const std::optional<V>& b) {
return a? a: b;
@@ -81,170 +76,23 @@ std::vector<V> unique_union(const std::vector<V>& a, const std::vector<V>& b) {
}
return u;
}
-
-// Convert mcable_map values to a piecewise function over an mcable.
-// The projection gives the map from the values in the mcable_map to the values in the piecewise function.
-template <typename T, typename U, typename Proj = get_value>
-pw_elements<U> pw_over_cable(const mcable_map<T>& mm, mcable cable, U dflt_value, Proj projection = Proj{}) {
- using value_type = typename mcable_map<T>::value_type;
- msize_t bid = cable.branch;
-
- struct as_branch {
- msize_t value;
- as_branch(const value_type& x): value(x.first.branch) {}
- as_branch(const msize_t& x): value(x) {}
- };
-
- auto map_on_branch = util::make_range(
- std::equal_range(mm.begin(), mm.end(), bid,
- [](as_branch a, as_branch b) { return a.value<b.value; }));
-
- if (map_on_branch.empty()) {
- return pw_elements<U>({cable.prox_pos, cable.dist_pos}, {dflt_value});
- }
-
- pw_elements<U> pw;
- for (const auto& el: map_on_branch) {
- double pw_right = pw.empty()? 0: pw.bounds().second;
- if (el.first.prox_pos>pw_right) {
- pw.push_back(pw_right, el.first.prox_pos, dflt_value);
- }
- pw.push_back(el.first.prox_pos, el.first.dist_pos, projection(el.second));
- }
-
- double pw_right = pw.empty()? 0: pw.bounds().second;
- if (pw_right<1.) {
- pw.push_back(pw_right, 1., dflt_value);
- }
-
- if (cable.prox_pos!=0 || cable.dist_pos!=1) {
- pw = pw_zip_with(pw, pw_elements<>({cable.prox_pos, cable.dist_pos}));
- }
- return pw;
-}
} // anonymous namespace
// Building CV geometry
// --------------------
-// Construct cv_geometry for cell from locset describing CV boundary points.
-
-cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
- auto pop = [](auto& vec) { auto h = vec.back(); return vec.pop_back(), h; };
-
- cv_geometry geom;
- const auto& mp = cell.provider();
- const auto& m = mp.morphology();
-
- if (m.empty()) {
- geom.cell_cv_divs = {0, 0};
- return geom;
- }
-
- mlocation_list locs = thingify(lset, mp);
-
- // Filter out root, terminal locations and repeated locations so as to
- // avoid trivial CVs outside of fork points. (This is not necessary for
- // correctness, but is for the convenience of specification by lset.)
-
- auto neither_root_nor_terminal = [&m](mlocation x) {
- return !(x.pos==0 && x.branch==(m.branch_children(mnpos).size()>1u? mnpos: 0)) // root?
- && !(x.pos==1 && m.branch_children(x.branch).empty()); // terminal?
- };
- locs.erase(std::partition(locs.begin(), locs.end(), neither_root_nor_terminal), locs.end());
- util::sort(locs);
- util::unique_in_place(locs);
-
- // Collect cables constituting each CV, maintaining a stack of CV
- // proximal 'head' points, and recursing down branches in the morphology
- // within each CV.
-
- constexpr fvm_index_type no_parent = -1;
- std::vector<std::pair<mlocation, fvm_index_type>> next_cv_head; // head loc, parent cv index
- next_cv_head.emplace_back(mlocation{mnpos, 0}, no_parent);
-
- mcable_list cables;
- std::vector<msize_t> branches;
- geom.cv_cables_divs.push_back(0);
- fvm_index_type cv_index = 0;
-
- while (!next_cv_head.empty()) {
- auto next = pop(next_cv_head);
- mlocation h = next.first;
-
- cables.clear();
- branches.clear();
- branches.push_back(h.branch);
- geom.cv_parent.push_back(next.second);
-
- while (!branches.empty()) {
- msize_t b = pop(branches);
-
- // Find most proximal point in locs on this branch, strictly more distal than h.
- auto it = locs.end();
- if (b!=mnpos && b==h.branch) {
- it = std::upper_bound(locs.begin(), locs.end(), h);
- }
- else if (b!=mnpos) {
- it = std::lower_bound(locs.begin(), locs.end(), mlocation{b, 0});
- }
-
- // If found, use as an end point, and stop descent.
- // Otherwise, recurse over child branches.
- if (it!=locs.end() && it->branch==b) {
- cables.push_back({b, b==h.branch? h.pos: 0, it->pos});
- next_cv_head.emplace_back(*it, cv_index);
- }
- else {
- if (b!=mnpos) {
- cables.push_back({b, b==h.branch? h.pos: 0, 1});
- }
- for (auto& c: m.branch_children(b)) {
- branches.push_back(c);
- }
- }
- }
-
- sort(cables);
- util::append(geom.cv_cables, std::move(cables));
- geom.cv_cables_divs.push_back(geom.cv_cables.size());
- ++cv_index;
- }
-
- auto n_cv = cv_index;
- arb_assert(n_cv>0);
- arb_assert(geom.cv_parent.front()==-1);
- arb_assert(util::all_of(util::subrange_view(geom.cv_parent, 1, n_cv),
- [](auto v) { return v!=no_parent; }));
-
- // Construct CV children mapping by sorting CV indices by parent.
- assign(geom.cv_children, make_span(1, n_cv));
- stable_sort_by(geom.cv_children, [&geom](auto cv) { return geom.cv_parent[cv]; });
-
- geom.cv_children_divs.reserve(n_cv+1);
- geom.cv_children_divs.push_back(0);
-
- auto b = geom.cv_children.begin();
- auto e = geom.cv_children.end();
- auto from = b;
-
- for (fvm_index_type cv = 0; cv<n_cv; ++cv) {
- from = std::partition_point(from, e,
- [cv, &geom](auto i) { return geom.cv_parent[i]<=cv; });
- geom.cv_children_divs.push_back(from-b);
- }
-
- // Fill cv/cell mapping for single cell (index 0).
- geom.cv_to_cell.assign(cv_index, 0);
- geom.cell_cv_divs = {0, (fvm_index_type)cv_index};
+// CV geometry
+cv_geometry::cv_geometry(const cable_cell& cell, const locset& ls):
+ base(cell, ls)
+{
// Build location query map.
- geom.branch_cv_map.resize(1);
- std::vector<pw_elements<fvm_size_type>>& bmap = geom.branch_cv_map.back();
-
- for (auto cv: make_span(n_cv)) {
- for (auto cable: geom.cables(cv)) {
+ auto n_cv = cv_parent.size();
+ branch_cv_map.resize(1);
+ std::vector<util::pw_elements<fvm_size_type>>& bmap = branch_cv_map.back();
+ for (auto cv: util::make_span(n_cv)) {
+ for (auto cable: cables(cv)) {
if (cable.branch>=bmap.size()) {
bmap.resize(cable.branch+1);
}
@@ -253,8 +101,43 @@ cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
bmap[cable.branch].push_back(cable.prox_pos, cable.dist_pos, cv);
}
}
+ cv_to_cell.assign(n_cv, 0);
+ cell_cv_divs = {0, (fvm_index_type)n_cv};
+}
- return geom;
+fvm_size_type cv_geometry::location_cv(size_type cell_idx, mlocation loc, cv_prefer::type prefer) const {
+ auto& pw_cv_offset = branch_cv_map.at(cell_idx).at(loc.branch);
+ auto zero_extent = [&pw_cv_offset](auto j) {
+ return pw_cv_offset.extent(j).first==pw_cv_offset.extent(j).second;
+ };
+
+ auto i = pw_cv_offset.index_of(loc.pos);
+ auto i_max = pw_cv_offset.size()-1;
+ auto cv_prox = pw_cv_offset.extent(i).first;
+
+ // index_of() should have returned right-most matching interval.
+ arb_assert(i==i_max || loc.pos<pw_cv_offset.extent(i+1).first);
+
+ using namespace cv_prefer;
+ switch (prefer) {
+ case cv_distal:
+ break;
+ case cv_proximal:
+ if (loc.pos==cv_prox && i>0) --i;
+ break;
+ case cv_nonempty:
+ if (zero_extent(i)) {
+ if (i>0 && !zero_extent(i-1)) --i;
+ else if (i<i_max && !zero_extent(i+1)) ++i;
+ }
+ break;
+ case cv_empty:
+ if (loc.pos==cv_prox && i>0 && zero_extent(i-1)) --i;
+ break;
+ }
+
+ index_type cv_base = cell_cv_divs.at(cell_idx);
+ return cv_base+pw_cv_offset.value(i);
}
namespace impl {
@@ -338,7 +221,6 @@ fvm_cv_discretization& append(fvm_cv_discretization& dczn, const fvm_cv_discreti
return dczn;
}
-
// FVM discretization
// ------------------
@@ -346,7 +228,7 @@ fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell
const auto& dflt = cell.default_parameters();
fvm_cv_discretization D;
- D.geometry = cv_geometry_from_ends(cell,
+ D.geometry = cv_geometry(cell,
dflt.discretization? dflt.discretization->cv_boundary_points(cell):
global_dflt.discretization? global_dflt.discretization->cv_boundary_points(cell):
default_cv_policy().cv_boundary_points(cell));
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index 8d1356fb9148d02f93672060d47f67b0585e12bd..7d178dc6112fc04df372fbc4669ea791200c37c3 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -9,8 +9,10 @@
#include <arbor/mechinfo.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/recipe.hpp>
+#include <arbor/morph/cv_data.hpp>
#include "execution_context.hpp"
+#include "morph/cv_data.hpp"
#include "util/piecewise.hpp"
#include "util/rangeutil.hpp"
#include "util/span.hpp"
@@ -53,23 +55,20 @@ namespace cv_prefer {
};
}
-struct cv_geometry {
+struct cv_geometry: public cell_cv_data_impl {
+ using base = cell_cv_data_impl;
+
using size_type = fvm_size_type;
using index_type = fvm_index_type;
- std::vector<mcable> cv_cables; // CV unbranched sections, partitioned by CV.
- std::vector<index_type> cv_cables_divs; // Partitions cv_cables by CV index.
- std::vector<index_type> cv_parent; // Index of CV parent or size_type(-1) for a cell root CV.
-
- std::vector<index_type> cv_children; // CV child indices, partitioned by CV, and then in order.
- std::vector<index_type> cv_children_divs; // Paritions cv_children by CV index.
-
std::vector<index_type> cv_to_cell; // Maps CV index to cell index.
std::vector<index_type> cell_cv_divs; // Partitions CV indices by cell.
// CV offset map by cell index then branch. Used for location_cv query.
std::vector<std::vector<util::pw_elements<size_type>>> branch_cv_map;
+ cv_geometry() = default;
+
auto cables(size_type cv_index) const {
auto partn = util::partition_view(cv_cables_divs);
return util::subrange_view(cv_cables, partn[cv_index]);
@@ -113,49 +112,15 @@ struct cv_geometry {
return branch_cv_map.at(cell_idx).size();
}
- size_type location_cv(size_type cell_idx, mlocation loc, cv_prefer::type prefer) const {
- auto& pw_cv_offset = branch_cv_map.at(cell_idx).at(loc.branch);
- auto zero_extent = [&pw_cv_offset](auto j) {
- return pw_cv_offset.extent(j).first==pw_cv_offset.extent(j).second;
- };
-
- auto i = pw_cv_offset.index_of(loc.pos);
- auto i_max = pw_cv_offset.size()-1;
- auto cv_prox = pw_cv_offset.extent(i).first;
-
- // index_of() should have returned right-most matching interval.
- arb_assert(i==i_max || loc.pos<pw_cv_offset.extent(i+1).first);
-
- using namespace cv_prefer;
- switch (prefer) {
- case cv_distal:
- break;
- case cv_proximal:
- if (loc.pos==cv_prox && i>0) --i;
- break;
- case cv_nonempty:
- if (zero_extent(i)) {
- if (i>0 && !zero_extent(i-1)) --i;
- else if (i<i_max && !zero_extent(i+1)) ++i;
- }
- break;
- case cv_empty:
- if (loc.pos==cv_prox && i>0 && zero_extent(i-1)) --i;
- break;
- }
-
- index_type cv_base = cell_cv_divs.at(cell_idx);
- return cv_base+pw_cv_offset.value(i);
- }
+ size_type location_cv(size_type cell_idx, mlocation loc, cv_prefer::type prefer) const;
+
+ cv_geometry(const cable_cell& cell, const locset& ls);
};
// Combine two cv_geometry groups in-place.
// (Returns reference to first argument.)
cv_geometry& append(cv_geometry&, const cv_geometry&);
-// Construct cv_geometry from locset describing boundaries.
-cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset);
-
// Discretization of morphologies and physical properties. Contains cv_geometry
// as above.
//
diff --git a/arbor/include/arbor/morph/cv_data.hpp b/arbor/include/arbor/morph/cv_data.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..58dc1e7c139d04abe9ef3856a79e15bc5a12ef9a
--- /dev/null
+++ b/arbor/include/arbor/morph/cv_data.hpp
@@ -0,0 +1,54 @@
+#pragma once
+
+#include <vector>
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/fvm_types.hpp>
+#include <arbor/morph/embed_pwlin.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/region.hpp>
+
+namespace arb {
+
+struct cell_cv_data_impl;
+
+// Stores info about the CV geometry of a discretized cable-cell
+class cell_cv_data {
+public:
+ // Returns mcables comprising the CV at a given index.
+ mcable_list cables(fvm_size_type index) const;
+
+ // Returns the CV indices of the children of a given CV index.
+ std::vector<fvm_index_type> children(fvm_size_type index) const;
+
+ // Returns the CV index of the parent of a given CV index.
+ fvm_index_type parent(fvm_size_type index) const;
+
+ // Returns total number of CVs.
+ fvm_size_type size() const;
+
+ cell_cv_data(const cable_cell& cell, const locset& lset);
+
+ const mprovider& provider() const {
+ return provider_;
+ }
+
+private:
+ std::unique_ptr<cell_cv_data_impl, void (*)(cell_cv_data_impl*)> impl_;
+
+ // Embedded morphology and labelled region/locset lookup.
+ mprovider provider_;
+};
+
+struct cv_proportion {
+ fvm_size_type idx;
+ fvm_value_type proportion;
+};
+
+// Construct cell_cv_geometry for cell from default cell discretization if it exists.
+std::optional<cell_cv_data> cv_data(const cable_cell& cell);
+
+std::vector<cv_proportion> intersect_region(const region& reg, const cell_cv_data& cvs, bool intergrate_by_length = false);
+
+} //namespace arb
diff --git a/arbor/morph/cv_data.cpp b/arbor/morph/cv_data.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..410bac0cc4bc1d1db0f2b6936f4ea3d9563b1925
--- /dev/null
+++ b/arbor/morph/cv_data.cpp
@@ -0,0 +1,190 @@
+#include <vector>
+
+#include <arbor/common_types.hpp>
+#include <arbor/morph/cv_data.hpp>
+#include <arbor/morph/embed_pwlin.hpp>
+
+#include "morph/cv_data.hpp"
+#include "util/partition.hpp"
+#include "util/pw_over_cable.hpp"
+#include "util/rangeutil.hpp"
+#include "util/span.hpp"
+#include "util/unique.hpp"
+
+namespace arb {
+
+cell_cv_data_impl::cell_cv_data_impl(const cable_cell& cell, const locset& lset) {
+ auto pop = [](auto& vec) { auto h = vec.back(); return vec.pop_back(), h; };
+
+ const auto& mp = cell.provider();
+ const auto& m = mp.morphology();
+
+ if (m.empty()) {
+ return;
+ }
+
+ mlocation_list locs = thingify(lset, mp);
+
+ // Filter out root, terminal locations and repeated locations so as to
+ // avoid trivial CVs outside of fork points. (This is not necessary for
+ // correctness, but is for the convenience of specification by lset.)
+
+ auto neither_root_nor_terminal = [&m](mlocation x) {
+ return !(x.pos==0 && x.branch==(m.branch_children(mnpos).size()>1u? mnpos: 0)) // root?
+ && !(x.pos==1 && m.branch_children(x.branch).empty()); // terminal?
+ };
+ locs.erase(std::partition(locs.begin(), locs.end(), neither_root_nor_terminal), locs.end());
+ util::sort(locs);
+ util::unique_in_place(locs);
+
+ // Collect cables constituting each CV, maintaining a stack of CV
+ // proximal 'head' points, and recursing down branches in the morphology
+ // within each CV.
+
+ constexpr fvm_index_type no_parent = -1;
+ std::vector<std::pair<mlocation, fvm_index_type>> next_cv_head; // head loc, parent cv index
+ next_cv_head.emplace_back(mlocation{mnpos, 0}, no_parent);
+
+ mcable_list cables;
+ std::vector<msize_t> branches;
+ cv_cables_divs.push_back(0);
+ fvm_index_type cv_index = 0;
+
+ while (!next_cv_head.empty()) {
+ auto next = pop(next_cv_head);
+ mlocation h = next.first;
+
+ cables.clear();
+ branches.clear();
+ branches.push_back(h.branch);
+ cv_parent.push_back(next.second);
+
+ while (!branches.empty()) {
+ msize_t b = pop(branches);
+
+ // Find most proximal point in locs on this branch, strictly more distal than h.
+ auto it = locs.end();
+ if (b!=mnpos && b==h.branch) {
+ it = std::upper_bound(locs.begin(), locs.end(), h);
+ }
+ else if (b!=mnpos) {
+ it = std::lower_bound(locs.begin(), locs.end(), mlocation{b, 0});
+ }
+
+ // If found, use as an end point, and stop descent.
+ // Otherwise, recurse over child branches.
+ if (it!=locs.end() && it->branch==b) {
+ cables.push_back({b, b==h.branch? h.pos: 0, it->pos});
+ next_cv_head.emplace_back(*it, cv_index);
+ }
+ else {
+ if (b!=mnpos) {
+ cables.push_back({b, b==h.branch? h.pos: 0, 1});
+ }
+ for (auto& c: m.branch_children(b)) {
+ branches.push_back(c);
+ }
+ }
+ }
+
+ util::sort(cables);
+ util::append(cv_cables, std::move(cables));
+ cv_cables_divs.push_back(cv_cables.size());
+ ++cv_index;
+ }
+
+ auto n_cv = cv_index;
+ arb_assert(n_cv>0);
+ arb_assert(cv_parent.front()==-1);
+ arb_assert(util::all_of(util::subrange_view(cv_parent, 1, n_cv),
+ [](auto v) { return v!=no_parent; }));
+
+ // Construct CV children mapping by sorting CV indices by parent.
+ assign(cv_children, util::make_span(1, n_cv));
+ util::stable_sort_by(cv_children, [this](auto cv) { return cv_parent[cv]; });
+
+ cv_children_divs.reserve(n_cv+1);
+ cv_children_divs.push_back(0);
+
+ auto b = cv_children.begin();
+ auto e = cv_children.end();
+ auto from = b;
+
+ for (fvm_index_type cv = 0; cv<n_cv; ++cv) {
+ from = std::partition_point(from, e,
+ [cv, this](auto i) { return cv_parent[i]<=cv; });
+ cv_children_divs.push_back(from-b);
+ }
+}
+
+mcable_list cell_cv_data::cables(fvm_size_type cv_index) const {
+ auto partn = util::partition_view(impl_->cv_cables_divs);
+ auto view = util::subrange_view(impl_->cv_cables, partn[cv_index]);
+ return mcable_list{view.begin(), view.end()};
+}
+
+std::vector<fvm_index_type> cell_cv_data::children(fvm_size_type cv_index) const {
+ auto partn = util::partition_view(impl_->cv_children_divs);
+ auto view = util::subrange_view(impl_->cv_children, partn[cv_index]);
+ return std::vector<fvm_index_type>{view.begin(), view.end()};
+}
+
+fvm_index_type cell_cv_data::parent(fvm_size_type cv_index) const {
+ return impl_->cv_parent[cv_index];
+}
+
+fvm_size_type cell_cv_data::size() const {
+ return impl_->cv_parent.size();
+}
+
+std::optional<cell_cv_data> cv_data(const cable_cell& cell) {
+ if (auto policy = cell.decorations().defaults().discretization) {
+ return cell_cv_data(cell, policy->cv_boundary_points(cell));
+ }
+ return {};
+}
+
+using impl_ptr = std::unique_ptr<cell_cv_data_impl, void (*)(cell_cv_data_impl*)>;
+impl_ptr make_impl(cell_cv_data_impl* c) {
+ return impl_ptr(c, [](cell_cv_data_impl* p){delete p;});
+}
+
+cell_cv_data::cell_cv_data(const cable_cell& cell, const locset& lset):
+ impl_(make_impl(new cell_cv_data_impl(cell, lset))),
+ provider_(cell.provider())
+{}
+
+std::vector<cv_proportion> intersect_region(const region& reg, const cell_cv_data& geom, bool by_length) {
+ const auto& mp = geom.provider();
+ const auto& embedding = mp.embedding();
+
+ std::vector<cv_proportion> intersection;
+ auto extent = thingify(reg, mp);
+ mcable_map<double> support;
+ for (auto& cable: extent) {
+ support.insert(cable, 1.);
+ }
+ if(support.empty()) {
+ return {};
+ }
+
+ for (auto cv: util::make_span(geom.size())) {
+ double cv_total = 0, reg_on_cv = 0;
+ for (mcable c: geom.cables(cv)) {
+ if (by_length) {
+ cv_total += embedding.integrate_length(c);
+ reg_on_cv += embedding.integrate_length(c.branch, util::pw_over_cable(support, c, 0.));
+ }
+ else {
+ cv_total += embedding.integrate_area(c);
+ reg_on_cv += embedding.integrate_area(c.branch, util::pw_over_cable(support, c, 0.));
+ }
+ }
+ if (reg_on_cv>0) {
+ intersection.push_back({cv, reg_on_cv/cv_total});
+ }
+ }
+ return intersection;
+}
+
+} //namespace arb
diff --git a/arbor/morph/cv_data.hpp b/arbor/morph/cv_data.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1e3e74713e804a93dc5528bc93b3cb0bf8845eb6
--- /dev/null
+++ b/arbor/morph/cv_data.hpp
@@ -0,0 +1,24 @@
+#pragma once
+
+#include <vector>
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/fvm_types.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/morphology.hpp>
+
+namespace arb {
+
+struct cell_cv_data_impl {
+ mcable_list cv_cables; // CV unbranched sections, partitioned by CV.
+ std::vector<fvm_index_type> cv_cables_divs; // Partitions cv_cables by CV index.
+
+ std::vector<fvm_index_type> cv_parent; // Index of CV parent or size_type(-1) for a cell root CV.
+ std::vector<fvm_index_type> cv_children; // CV child indices, partitioned by CV, and then in order.
+ std::vector<fvm_index_type> cv_children_divs; // Paritions cv_children by CV index.
+
+ cell_cv_data_impl() = default;
+ cell_cv_data_impl(const cable_cell&, const locset&);
+};
+
+} // namespace arb
diff --git a/arbor/util/pw_over_cable.hpp b/arbor/util/pw_over_cable.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9345b86b036b83470b75c6551031037128f2d44d
--- /dev/null
+++ b/arbor/util/pw_over_cable.hpp
@@ -0,0 +1,56 @@
+#include <arbor/morph/mcable_map.hpp>
+
+#include "util/piecewise.hpp"
+
+namespace arb {
+namespace util {
+
+namespace impl {
+struct get_value {
+ template <typename X>
+ double operator()(const X& x) const { return x.value; }
+ double operator()(double x) const { return x; }
+};
+} // namespace impl
+
+// Convert mcable_map values to a piecewise function over an mcable.
+// The projection gives the map from the values in the mcable_map to the values in the piecewise function.
+template <typename T, typename U, typename Proj = impl::get_value>
+util::pw_elements<U> pw_over_cable(const mcable_map<T>& mm, mcable cable, U dflt_value, Proj projection = Proj{}) {
+ using value_type = typename mcable_map<T>::value_type;
+ msize_t bid = cable.branch;
+
+ struct as_branch {
+ msize_t value;
+ as_branch(const value_type& x): value(x.first.branch) {}
+ as_branch(const msize_t& x): value(x) {}
+ };
+
+ auto map_on_branch = util::make_range(
+ std::equal_range(mm.begin(), mm.end(), bid, [](as_branch a, as_branch b) { return a.value<b.value; }));
+
+ if (map_on_branch.empty()) {
+ return util::pw_elements<U>({cable.prox_pos, cable.dist_pos}, {dflt_value});
+ }
+
+ util::pw_elements<U> pw;
+ for (const auto& el: map_on_branch) {
+ double pw_right = pw.empty()? 0: pw.bounds().second;
+ if (el.first.prox_pos>pw_right) {
+ pw.push_back(pw_right, el.first.prox_pos, dflt_value);
+ }
+ pw.push_back(el.first.prox_pos, el.first.dist_pos, projection(el.second));
+ }
+
+ double pw_right = pw.empty()? 0: pw.bounds().second;
+ if (pw_right<1.) {
+ pw.push_back(pw_right, 1., dflt_value);
+ }
+
+ if (cable.prox_pos!=0 || cable.dist_pos!=1) {
+ pw = util::pw_zip_with(pw, util::pw_elements<>({cable.prox_pos, cable.dist_pos}));
+ }
+ return pw;
+}
+} // namespace util
+} // namespace arb
\ No newline at end of file
diff --git a/doc/cpp/morphology.rst b/doc/cpp/morphology.rst
index dc432b625b7c9b0d4998318aad2ba326dc9a72ff..f2ea79e56042342d31e1d96c04049747f68f3951 100644
--- a/doc/cpp/morphology.rst
+++ b/doc/cpp/morphology.rst
@@ -336,6 +336,58 @@ 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.
+CV discretization as mcables
+----------------------------
+
+It is often useful for the user to have a detailed view of the CVs generated for a given morphology
+and :ref:`cv-policy <cv-policies>`. For example, while debugging and fine-tuning model implementations,
+it can be helpful to know how many CVs a cable-cell is comprised of, or how many CVs lie on a certain
+region of the cell.
+
+The following classes and functions allow the user to inspect the CVs of a cell or region.
+
+.. cpp:class:: cell_cv_data
+
+ Stores the discretisation data of a cable-cell in terms of CVs and the :term:`mcables <mcable>` comprising each of
+ these CVs.
+
+ .. cpp:function:: mcable_list cables(unsigned idx) const
+
+ Returns an vector of :term:`mcable` representing the CV at a given index.
+
+ .. cpp:function:: std::vector<unsigned> children(unsigned idx) const
+
+ Returns a vector of the indices of the CVs representing the children of the CV at index ``idx``.
+
+ .. cpp:function:: unsigned parent(unsigned idx) const
+
+ Returns the index of the CV representing the parent of the CV at index ``idx``.
+
+ .. cpp:function:: unsigned size() const
+
+ Returns the total number of CVs on the cell.
+
+.. cpp:function:: std::optional<cell_cv_data> cv_data(const cable_cell& cell)
+
+ Constructs a :cpp:class:`cell_cv_data` object representing the CVs comprising the cable-cell according
+ to the :cpp:class:`cv_policy` provided in the :cpp:class:`decor` of the cell. Returns ``std::nullopt_t``
+ if no :cpp:class:`cv_policy` was provided in the decor.
+
+.. cpp:class:: cv_proportion
+
+ .. cpp:member:: unsigned idx
+
+ Index of the CV.
+
+ .. cpp:member:: double proportion
+
+ Proportion of the CV by area.
+
+.. cpp:function:: std::vector<cv_proportion> intersect_region(const region& reg, const cell_cv_data& cvs, bool integrate_by_length=false)
+
+ Returns a vector of :cpp:class:`cv_proportion` identifying the indices of the CVs from the :cpp:class:`cell_cv_data`
+ argument that lie in the provided region, and how much of each CV belongs to that region. The proportion of CV lying
+ in the region is the area proportion if ``integrate_by_length`` is false, otherwise, it is the length proportion.
Supported morphology formats
============================
diff --git a/doc/python/morphology.rst b/doc/python/morphology.rst
index 6e233e183835029186b7255006789d8fd23112df..caa3f3ccedbe32c775e0d0d0f88e127fb2fed54a 100644
--- a/doc/python/morphology.rst
+++ b/doc/python/morphology.rst
@@ -498,6 +498,61 @@ constitute part of the CV boundary point set.
:param float max_etent: The maximum length for generated CVs.
:param str domain: The region on which the policy is applied.
+CV discretization as mcables
+----------------------------
+
+It is often useful for the user to have a detailed view of the CVs generated for a
+given morphology and :ref:`cv-policy <pymorph-cv-policies>`. For example, while
+debugging and fine-tuning model implementations, it can be helpful to know how many CVs
+a cable-cell is comprised of, or how many CVs lie on a certain region of the cell.
+
+The following classes and functions allow the user to inspect the CVs of a cell or
+region.
+
+.. py:class:: cell_cv_data
+
+ Stores the discretisation data of a cable-cell in terms of CVs and the :py:class:`cables <cable>`
+ comprising each of these CVs.
+
+ .. py:method:: cables(idx) -> list[cable]
+
+ Returns a list of :py:class:`cable` representing the CV at a given index ``idx``.
+
+ .. py:method:: children(idx) -> list[int]
+
+ Returns a list of the indices of the CVs representing the children of the CV at index ``idx``.
+
+ .. py:method:: parent(idx) -> int
+
+ Returns the index of the CV representing the parent of the CV at index ``idx``.
+
+ .. py:attribute:: int num_cv
+
+ Returns the total number of CVs on the cell.
+
+.. py:function:: cv_data(cell) -> optional<cell_cv_data>
+
+ Constructs a :py:class:`cell_cv_data` object representing the CVs comprising the cable-cell according
+ to the :py:class:`cv_policy` provided in the :py:class:`decor` of the cell. Returns ``None`` if no
+ :py:class:`cv_policy` was provided in the decor.
+
+ :param cable_cell cell: The cable-cell.
+ :rtype: optional<:py:class:`cell_cv_data`>
+
+.. py:function:: intersect_region(reg, cv_data, integrate_along) -> list[idx, proportion]
+
+ Returns a list of tuples ``[idx, proportion]`` identifying the indices (``idx``) of the CVs from the
+ ``cv_data`` argument that lie in the provided region ``reg``, and how much of each CV belongs to that
+ region (``proportion``). The ``proportion`` is either the area proportion or the length proportion,
+ chosen according to the ``integrate_along`` argument.
+
+ :param str reg: The region on the cable-cell represented as s-expression or a label from the
+ label-dictionary of the cell.
+ :param cell_cv_data cv_data: The cv_data of a cell.
+ :param string integrate_along: Either "area" or "length". Decides whether the proportion of a
+ CV is measured according to the area or length of the CV.
+ :rtype: list[idx, proportion]
+
.. _pyswc:
SWC
diff --git a/python/cells.cpp b/python/cells.cpp
index 1c26f78348b5230669293a96767d42d30b450e8a..da9060d0e95eba820d67de8e3aaaa25685a04a06 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -18,6 +18,7 @@
#include <arbor/cable_cell.hpp>
#include <arbor/lif_cell.hpp>
#include <arbor/cv_policy.hpp>
+#include <arbor/morph/cv_data.hpp>
#include <arbor/morph/label_dict.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/morph/region.hpp>
@@ -345,6 +346,64 @@ void register_cells(pybind11::module& m) {
"domain"_a="(all)", "the domain to which the policy is to be applied",
"Policy to use the same number of CVs for each branch.");
+ // arb::cell_cv_data
+ pybind11::class_<arb::cell_cv_data> cell_cv_data(m, "cell_cv_data",
+ "Provides information on the CVs representing the discretization of a cable-cell.");
+ cell_cv_data
+ .def_property_readonly("num_cv", [](const arb::cell_cv_data& data){return data.size();},
+ "Return the number of CVs in the cell.")
+ .def("cables",
+ [](const arb::cell_cv_data& d, unsigned index) {
+ if (index >= d.size()) throw pybind11::index_error("index out of range");
+ return d.cables(index);
+ },
+ "index"_a, "Return a list of cables representing the CV at the given index.")
+ .def("children",
+ [](const arb::cell_cv_data& d, unsigned index) {
+ if (index >= d.size()) throw pybind11::index_error("index out of range");
+ return d.children(index);
+ },
+ "index"_a,
+ "Return a list of indices of the CVs representing the children of the CV at the given index.")
+ .def("parent",
+ [](const arb::cell_cv_data& d, unsigned index) {
+ if (index >= d.size()) throw pybind11::index_error("index out of range");
+ return d.parent(index);
+ },
+ "index"_a,
+ "Return the index of the CV representing the parent of the CV at the given index.")
+ .def("__str__", [](const arb::cell_cv_data& p){return "<arbor.cell_cv_data>";})
+ .def("__repr__", [](const arb::cell_cv_data& p){return "<arbor.cell_cv_data>";});
+
+ m.def("cv_data", [](const arb::cable_cell& cell) { return arb::cv_data(cell);},
+ "cell"_a, "the cable cell",
+ "Returns a cell_cv_data object representing the CVs comprising the cable-cell according to the "
+ "discretization policy provided in the decor of the cell. Returns None if no CV-policy was provided "
+ "in the decor."
+ );
+
+ m.def("intersect_region",
+ [](const char* reg, const arb::cell_cv_data& cvs, const std::string& integrate_along) {
+ bool integrate_area;
+ if (integrate_along == "area") integrate_area = true;
+ else if (integrate_along == "length") integrate_area = false;
+ else throw pyarb_error(util::pprintf("{} does not name a valid integration axis. "
+ "Only 'area' and 'length' are supported)", integrate_along));
+
+ auto object_vec = arb::intersect_region(arborio::parse_region_expression(reg).unwrap(), cvs, integrate_area);
+ auto tuple_vec = std::vector<pybind11::tuple>(object_vec.size());
+ std::transform(object_vec.begin(), object_vec.end(), tuple_vec.begin(),
+ [](const auto& t) { return pybind11::make_tuple(t.idx, t.proportion); });
+ return tuple_vec;
+ },
+ "reg"_a, "A region on a cell",
+ "data"_a, "The data representing the CVs of the cell.",
+ "integrate_along"_a, "the axis of integration used to determine the proportion of the CV belonging to the region.",
+ "Returns a list of [index, proportion] tuples identifying the CVs present in the region.\n"
+ "`index` is the index of the CV in the cell_cv_data object provided as an argument.\n"
+ "`proportion` is the proportion of the CV (itegrated by area or length) included in the region."
+ );
+
// arb::density
pybind11::class_<arb::density> density(m, "density", "For painting a density mechanism on a region.");
@@ -511,6 +570,8 @@ void register_cells(pybind11::module& m) {
},
"default NEURON cable_global_properties");
+ // arb::decor
+
pybind11::class_<arb::decor> decor(m, "decor",
"Description of the decorations to be applied to a cable cell, that is the painted,\n"
"placed and defaulted properties, mecahanisms, ion species etc.");
diff --git a/test/unit/test_cv_geom.cpp b/test/unit/test_cv_geom.cpp
index 13410e7f861659a5db90a353bb3fc0364bdb4515..2e7d58b476e0680d033ae3fd651b7319057ec2b7 100644
--- a/test/unit/test_cv_geom.cpp
+++ b/test/unit/test_cv_geom.cpp
@@ -2,10 +2,13 @@
#include <utility>
#include <arbor/cable_cell.hpp>
+#include <arbor/morph/cv_data.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/morph/region.hpp>
+#include <arborio/label_parse.hpp>
+
#include "fvm_layout.hpp"
#include "util/rangeutil.hpp"
@@ -71,7 +74,7 @@ TEST(cv_geom, empty) {
using namespace common_morphology;
cable_cell empty_cell{m_empty};
- cv_geometry geom = cv_geometry_from_ends(empty_cell, ls::nil());
+ cv_geometry geom(empty_cell, ls::nil());
EXPECT_TRUE(verify_cv_children(geom));
EXPECT_TRUE(geom.cv_parent.empty());
@@ -97,8 +100,8 @@ TEST(cv_geom, trivial) {
auto& m = cell.morphology();
// Equivalent ways of specifying one CV comprising whole cell:
- cv_geometry geom1 = cv_geometry_from_ends(cell, ls::nil());
- cv_geometry geom2 = cv_geometry_from_ends(cell, ls::terminal());
+ cv_geometry geom1(cell, ls::nil());
+ cv_geometry geom2(cell, ls::terminal());
EXPECT_TRUE(verify_cv_children(geom1));
EXPECT_TRUE(verify_cv_children(geom2));
@@ -107,8 +110,8 @@ TEST(cv_geom, trivial) {
EXPECT_EQ(geom1.cv_cables, geom2.cv_cables);
// These are equivalent too, if there is a single root branch.
- cv_geometry geom3 = cv_geometry_from_ends(cell, ls::root());
- cv_geometry geom4 = cv_geometry_from_ends(cell, join(ls::root(), ls::terminal()));
+ cv_geometry geom3(cell, ls::root());
+ cv_geometry geom4(cell, join(ls::root(), ls::terminal()));
EXPECT_TRUE(verify_cv_children(geom3));
EXPECT_TRUE(verify_cv_children(geom4));
@@ -137,8 +140,8 @@ TEST(cv_geom, one_cv_per_branch) {
cable_cell cell{p.second};
auto& m = cell.morphology();
- cv_geometry geom =
- cv_geometry_from_ends(cell, sum(ls::on_branches(0), ls::on_branches(1)));
+ auto cell_cv_geom = cv_geometry(cell, sum(ls::on_branches(0), ls::on_branches(1)));
+ auto geom = cv_geometry(cell_cv_geom);
EXPECT_TRUE(verify_cv_children(geom));
// Expect trivial CVs at every fork point, and single-cable CVs for each branch.
@@ -192,7 +195,7 @@ TEST(cv_geom, midpoints) {
cable_cell cell{p.second};
auto& m = cell.morphology();
- cv_geometry geom = cv_geometry_from_ends(cell, ls::on_branches(0.5));
+ cv_geometry geom(cell, ls::on_branches(0.5));
EXPECT_TRUE(verify_cv_children(geom));
// Expect CVs to be either: covering fork points, with one cable per branch
@@ -283,7 +286,7 @@ TEST(cv_geom, weird) {
using testing::seq_eq;
cable_cell cell{common_morphology::m_reg_b6};
- cv_geometry geom = cv_geometry_from_ends(cell, mlocation_list{{1, 0}, {4,0}});
+ cv_geometry geom(cell, mlocation_list{{1, 0}, {4,0}});
EXPECT_TRUE(verify_cv_children(geom));
ASSERT_EQ(3u, geom.size());
@@ -315,7 +318,7 @@ TEST(cv_geom, location_cv) {
};
// Two CVs per branch, plus trivial CV at forks.
- cv_geometry geom = cv_geometry_from_ends(cell,
+ cv_geometry geom(cell,
join(ls::on_branches(0.), ls::on_branches(0.5), ls::on_branches(1.)));
// Confirm CVs are either trivial or a single cable covering half a branch.
@@ -450,7 +453,7 @@ TEST(cv_geom, multicell) {
cable_cell cell = cable_cell(m_reg_b6);
- cv_geometry geom = cv_geometry_from_ends(cell, ls::on_branches(0.5));
+ cv_geometry geom(cell, ls::on_branches(0.5));
unsigned n_cv = geom.size();
cv_geometry geom2 = geom;
@@ -482,3 +485,222 @@ TEST(cv_geom, multicell) {
EXPECT_EQ((std::pair<index_type, index_type>(0, n_cv)), geom2.cell_cv_interval(0));
EXPECT_EQ((std::pair<index_type, index_type>(n_cv, 2*n_cv)), geom2.cell_cv_interval(1));
}
+
+TEST(region_cv, empty) {
+ using namespace common_morphology;
+
+ cable_cell empty_cell{m_empty};
+ cell_cv_data cv_data(empty_cell, ls::nil());
+
+ auto all_cv = intersect_region(reg::all(), cv_data);
+ EXPECT_EQ(0u, all_cv.size());
+
+ auto tag1_cv = intersect_region(reg::tagged(1), cv_data);
+ EXPECT_EQ(0u, tag1_cv.size());
+}
+
+TEST(region_cv, trivial) {
+ using namespace common_morphology;
+
+ for (auto& p: test_morphologies) {
+ if (p.second.empty()) continue;
+
+ SCOPED_TRACE(p.first);
+ cable_cell cell{p.second};
+
+ // One CV comprising whole cell:
+ cell_cv_data cell_geom1(cell, ls::nil());
+
+ auto all_cv = intersect_region(reg::all(), cell_geom1);
+ auto tag1_cv = intersect_region(reg::tagged(1), cell_geom1);
+ auto tag2_cv = intersect_region(reg::tagged(2), cell_geom1);
+
+ EXPECT_EQ(1u, all_cv.size());
+ EXPECT_EQ(0u, all_cv.front().idx);
+ EXPECT_EQ(1., all_cv.front().proportion);
+
+ EXPECT_EQ(1u, tag1_cv.size());
+ EXPECT_EQ(0u, tag1_cv.front().idx);
+ EXPECT_EQ(1., tag1_cv.front().proportion);
+
+ EXPECT_EQ(0u, tag2_cv.size());
+ }
+}
+
+TEST(region_cv, custom_geometry) {
+ auto almost_eq = [](const cv_proportion& a, const cv_proportion& b) {
+ return a.idx==b.idx && (fabs(a.proportion-b.proportion)<1e-6);
+ };
+
+ using namespace arborio::literals;
+ segment_tree tree;
+ // the test morphology:
+ //
+ // _ _ y=40
+ // \ /
+ // seg3 \ / seg2
+ // branch2 \ / branch1
+ // \ /
+ // | y=25
+ // |
+ // | seg1
+ // | branch0
+ // |
+ // |
+ // - y=5
+ // | seg0
+ // _ y=-5
+
+ // Root branch.
+ mpoint psoma_p{0, -5, 0, 5};
+ mpoint psoma_d{0, 5, 0, 5};
+
+ msize_t ssoma = tree.append(mnpos, psoma_p, psoma_d, 1);
+
+ mpoint py1_p{0, 5, 0, 1};
+ mpoint py1_d{0, 25, 0, 1};
+
+ msize_t sy1 = tree.append(ssoma, py1_p, py1_d, 3);
+
+ // branch 1 of y: translation (9,15) in one segment
+ mpoint py2_d{ 9, 40, 0, 1};
+ tree.append(sy1, py2_d, 3);
+
+ // branch 2 of y: translation (-9,15) in 2 segments
+ mpoint py3_d{-9, 40, 0, 1};
+ tree.append(sy1, py3_d, 3);
+
+ morphology m(tree);
+
+ label_dict l;
+ l.set("all", reg::all());
+ l.set("soma", reg::tagged(1));
+ l.set("dend", reg::tagged(3));
+ l.set("custom0", reg::cable(0, 0, 0.5));
+ l.set("custom1", reg::cable(1, 0.5, 1));
+ l.set("custom2", join(reg::named("custom0"), reg::named("custom1")));
+ l.set("custom3", reg::cable(0, 1./6., 3./5.));
+ l.set("custom4", reg::cable(0, 3./5., 1));
+ {
+ decor d;
+ // Discretize by segment
+ d.set_default(cv_policy_every_segment());
+ auto cell = cable_cell(m, l, d);
+ auto geom = cv_data(cell);
+ EXPECT_TRUE(geom);
+
+ auto cv0 = cv_proportion{0, 1};
+ auto cv1 = cv_proportion{1, 1};
+ auto cv1_quarter = cv_proportion{1, 0.25};
+ auto cv2 = cv_proportion{3, 1};
+ auto cv3 = cv_proportion{4, 1};
+ auto cv2_half = cv_proportion{3, 0.5};
+
+ auto all_cv = intersect_region("all"_lab, geom.value());
+ EXPECT_EQ(4u, all_cv.size());
+ EXPECT_TRUE(almost_eq(cv0, all_cv[0]));
+ EXPECT_TRUE(almost_eq(cv1, all_cv[1]));
+ EXPECT_TRUE(almost_eq(cv2, all_cv[2]));
+ EXPECT_TRUE(almost_eq(cv3, all_cv[3]));
+
+ auto soma_cv = intersect_region("soma"_lab, geom.value());
+ EXPECT_EQ(1u, soma_cv.size());
+ EXPECT_TRUE(almost_eq(cv0, soma_cv[0]));
+
+ auto dend_cv = intersect_region("dend"_lab, geom.value());
+ EXPECT_EQ(3u, dend_cv.size());
+ EXPECT_TRUE(almost_eq(cv1, dend_cv[0]));
+ EXPECT_TRUE(almost_eq(cv2, dend_cv[1]));
+ EXPECT_TRUE(almost_eq(cv3, dend_cv[2]));
+
+ auto c0_cv = intersect_region("custom0"_lab, geom.value());
+ EXPECT_EQ(2u, c0_cv.size());
+ EXPECT_TRUE(almost_eq(cv0, c0_cv[0]));
+ EXPECT_TRUE(almost_eq(cv1_quarter, c0_cv[1]));
+
+ auto c1_cv = intersect_region("custom1"_lab, geom.value());
+ EXPECT_EQ(1u, c1_cv.size());
+ EXPECT_TRUE(almost_eq(cv2_half, c1_cv[0]));
+
+ auto c2_cv = intersect_region("custom2"_lab, geom.value());
+ EXPECT_EQ(3u, c2_cv.size());
+ EXPECT_TRUE(almost_eq(cv0, c2_cv[0]));
+ EXPECT_TRUE(almost_eq(cv1_quarter, c2_cv[1]));
+ EXPECT_TRUE(almost_eq(cv2_half, c2_cv[2]));
+ }
+ {
+ decor d;
+ // Discretize using explicit locset
+ auto ls = locset({
+ {0, 0},
+ {0, 0.1},
+ {0, 0.5},
+ {0, 0.7},
+ {0, 0.9},
+ {1, 0.1},
+ {1, 1},
+ {2, 0.2},
+ {2, 1}
+ });
+ d.set_default(cv_policy_explicit(ls));
+ auto cell = cable_cell(m, l, d);
+ auto geom = cv_data(cell);
+ EXPECT_TRUE(geom);
+
+ auto cv0 = cv_proportion{0, 1};
+ auto cv1 = cv_proportion{1, 1};
+ auto cv1_part0 = cv_proportion{1, 7./8.};
+ auto cv1_part1 = cv_proportion{1, 1./8.};
+ auto cv1_part2 = cv_proportion{1, 3./4.};
+ auto cv2 = cv_proportion{2, 1};
+ auto cv2_part = cv_proportion{2, 1./2.};
+ auto cv3 = cv_proportion{3, 1};
+ auto cv4 = cv_proportion{4, 1};
+ auto cv4_part = cv_proportion{4, 10/(sqrt(306)+10)};
+ auto cv5 = cv_proportion{5, 1};
+ auto cv5_part = cv_proportion{5, 15./27.};
+ auto cv6 = cv_proportion{6, 1};
+
+ auto all_cv = intersect_region("all"_lab, geom.value());
+ EXPECT_EQ(7u, all_cv.size());
+ EXPECT_TRUE(almost_eq(cv0, all_cv[0]));
+ EXPECT_TRUE(almost_eq(cv1, all_cv[1]));
+ EXPECT_TRUE(almost_eq(cv2, all_cv[2]));
+ EXPECT_TRUE(almost_eq(cv3, all_cv[3]));
+ EXPECT_TRUE(almost_eq(cv4, all_cv[4]));
+ EXPECT_TRUE(almost_eq(cv5, all_cv[5]));
+ EXPECT_TRUE(almost_eq(cv6, all_cv[6]));
+
+ auto soma_cv = intersect_region("soma"_lab, geom.value());
+ EXPECT_EQ(2u, soma_cv.size());
+ EXPECT_TRUE(almost_eq(cv0, soma_cv[0]));
+ EXPECT_TRUE(almost_eq(cv1_part0, soma_cv[1]));
+
+ auto dend_cv = intersect_region("dend"_lab, geom.value());
+ EXPECT_EQ(6u, dend_cv.size());
+ EXPECT_TRUE(almost_eq(cv1_part1, dend_cv[0]));
+ EXPECT_TRUE(almost_eq(cv2, dend_cv[1]));
+ EXPECT_TRUE(almost_eq(cv3, dend_cv[2]));
+ EXPECT_TRUE(almost_eq(cv4, dend_cv[3]));
+ EXPECT_TRUE(almost_eq(cv5, dend_cv[4]));
+ EXPECT_TRUE(almost_eq(cv6, dend_cv[5]));
+ EXPECT_TRUE(almost_eq(cv6, dend_cv[5]));
+
+ auto c2_cv = intersect_region("custom2"_lab, geom.value());
+ EXPECT_EQ(3u, c2_cv.size());
+ EXPECT_TRUE(almost_eq(cv0, c2_cv[0]));
+ EXPECT_TRUE(almost_eq(cv1, c2_cv[1]));
+ EXPECT_TRUE(almost_eq(cv5_part, c2_cv[2]));
+
+ auto c3_cv = intersect_region("custom3"_lab, geom.value());
+ EXPECT_EQ(2u, c3_cv.size());
+ EXPECT_TRUE(almost_eq(cv1_part2, c3_cv[0]));
+ EXPECT_TRUE(almost_eq(cv2_part, c3_cv[1]));
+
+ auto c4_cv = intersect_region("custom4"_lab, geom.value());
+ EXPECT_EQ(3u, c2_cv.size());
+ EXPECT_TRUE(almost_eq(cv2_part, c4_cv[0]));
+ EXPECT_TRUE(almost_eq(cv3, c4_cv[1]));
+ EXPECT_TRUE(almost_eq(cv4_part, c4_cv[2]));
+ }
+}
diff --git a/test/unit/test_morph_place.cpp b/test/unit/test_morph_place.cpp
index fff37b7461c7da3a956cdaba9d8409f903e6256e..bd353fd23005d4850915b922246a9758543ff78a 100644
--- a/test/unit/test_morph_place.cpp
+++ b/test/unit/test_morph_place.cpp
@@ -581,7 +581,7 @@ TEST(place_pwlin, nearest) {
// _ _ y=40
// \ /
// seg4 \ / seg2
- // branch 2 / / branch 1
+ // branch 2 \ / branch 1
// seg3 \ /
// | y=25
// |
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index 02b570b1219ab12a987cb2bbcc73f3b83c96bf44..4608bb1d4c938e634dcf47517e212f5ddba7b162 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -240,7 +240,7 @@ void run_v_cell_probe_test(const context& ctx) {
// Independetly discretize the cell so we can follow cable–CV relationship.
- cv_geometry geom = cv_geometry_from_ends(cell, testcase.second.cv_boundary_points(cell));
+ cv_geometry geom(cell, testcase.second.cv_boundary_points(cell));
// For each cable in metadata, get CV from geom and confirm raw handle is
// state voltage + CV.
@@ -414,8 +414,8 @@ void run_expsyn_g_cell_probe_test(const context& ctx) {
// Independently get cv geometry to compute CV indices.
- cv_geometry geom = cv_geometry_from_ends(cells[0], policy.cv_boundary_points(cells[0]));
- append(geom, cv_geometry_from_ends(cells[1], policy.cv_boundary_points(cells[1])));
+ cv_geometry geom(cells[0], policy.cv_boundary_points(cells[0]));
+ append(geom, cv_geometry(cells[1], policy.cv_boundary_points(cells[1])));
ASSERT_EQ(2u, probe_map.size());
for (unsigned i: {0u, 1u}) {