diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index 65088e389a70454a0fcfcb3e7ae806eda8527ca6..f90e5bffdb92a28106aa21f279e98ec9662f5e19 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -6,7 +6,6 @@
#include <arbor/morph/label_dict.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/mprovider.hpp>
-#include <arbor/segment.hpp>
#include <arbor/util/pp_util.hpp>
#include "util/piecewise.hpp"
@@ -32,98 +31,20 @@ struct cable_cell_impl {
using index_type = cable_cell::index_type;
using size_type = cable_cell::size_type;
- cable_cell_impl(const arb::morphology& m,
- const label_dict& dictionary,
- bool compartments_from_discretization):
+ cable_cell_impl(const arb::morphology& m, const label_dict& dictionary):
provider(m, dictionary)
- {
- using point = cable_cell::point_type;
- if (!m.num_branches()) {
- segments.push_back(make_segment<placeholder_segment>());
- parents.push_back(0);
- return;
- }
-
- // Add the soma.
- auto loc = m.samples()[0].loc; // location of soma.
-
- // If there is no spherical root/soma use a zero-radius soma.
- double srad = m.spherical_root()? loc.radius: 0.;
- segments.push_back(make_segment<soma_segment>(srad, point(loc.x, loc.y, loc.z)));
- parents.push_back(-1);
-
- auto& samples = m.samples();
- auto& props = m.sample_props();
- for (auto i: util::make_span(1, m.num_branches())) {
- auto index = util::make_range(m.branch_indexes(i));
-
- // find kind for the branch. Use the tag of the last sample in the branch.
- int tag = samples[index.back()].tag;
- section_kind kind;
- switch (tag) {
- case 1: // soma
- throw cable_cell_error("No support for complex somata (yet)");
- case 2: // axon
- kind = section_kind::axon;
- case 3: // dendrite
- case 4: // apical dendrite
- default: // just take dendrite as default
- kind = section_kind::dendrite;
- }
-
- std::vector<value_type> radii;
- std::vector<cable_cell::point_type> points;
-
- // The current discretization code does not handle collocated points correctly,
- // particularly if they lie at the start of a branch, so we have to skip the first
- // point on a branch if it is collocated with the second point.
- bool skip_first = is_collocated(props[index[1]]);
- for (auto j: util::make_span(skip_first, index.size())) {
- auto& s = samples[index[j]];
- radii.push_back(s.loc.radius);
- points.push_back(cable_cell::point_type(s.loc.x, s.loc.y, s.loc.z));
- }
-
- // Find the id of this branch's parent.
- auto pid = m.branch_parent(i);
- // Adjust pid if a zero-radius soma was used.
- if (!m.spherical_root()) {
- pid = pid==mnpos? 0: pid+1;
- }
- segments.push_back(make_segment<cable_segment>(kind, radii, points));
- parents.push_back(pid);
- if (compartments_from_discretization) {
- int ncolloc = std::count_if(index.begin(), index.end(), [&props](auto i){return is_collocated(props[i]);});
- int ncomp = index.size()-ncolloc-1;
- ncomp -= is_collocated(props[index[0]]);
- segments.back()->as_cable()->set_compartments(ncomp);
- }
- }
- }
+ {}
- cable_cell_impl(): cable_cell_impl({},{},false) {}
+ cable_cell_impl(): cable_cell_impl({},{}) {}
cable_cell_impl(const cable_cell_impl& other):
- parents(other.parents),
provider(other.provider),
region_map(other.region_map),
location_map(other.location_map)
- {
- // unique_ptr's cannot be copy constructed, do a manual assignment
- segments.reserve(other.segments.size());
- for (const auto& s: other.segments) {
- segments.push_back(s->clone());
- }
- }
+ {}
cable_cell_impl(cable_cell_impl&& other) = default;
- // storage for connections
- std::vector<index_type> parents;
-
- // the segments
- std::vector<segment_ptr> segments;
-
// Embedded morphology and labelled region/locset lookup.
mprovider provider;
@@ -158,36 +79,6 @@ struct cable_cell_impl {
return lid_range(first, lid);
}
- void assert_valid_segment(index_type i) const {
- if (i>=segments.size()) {
- throw cable_cell_error("no such segment");
- }
- }
-
- void paint_segment(segment_ptr& s, const mechanism_desc& p) {
- s->add_mechanism(p);
- }
-
- void paint_segment(segment_ptr& s, init_membrane_potential p) {
- s->parameters.init_membrane_potential = p.value;
- }
-
- void paint_segment(segment_ptr& s, axial_resistivity p) {
- s->parameters.axial_resistivity = p.value;
- }
-
- void paint_segment(segment_ptr& s, temperature_K p) {
- s->parameters.temperature_K = p.value;
- }
-
- void paint_segment(segment_ptr& s, membrane_capacitance p) {
- s->parameters.membrane_capacitance = p.value;
- }
-
- void paint_segment(segment_ptr& s, initial_ion_data p) {
- s->parameters.ion_data[p.ion] = p.initial;
- }
-
template <typename T>
mcable_map<T>& get_region_map(const T&) {
return region_map.get<T>();
@@ -210,13 +101,6 @@ struct cable_cell_impl {
if (!mm.insert(c, prop)) {
throw cable_cell_error(util::pprintf("cable {} overpaints", c));
}
-
- if (c.prox_pos!=0 || c.dist_pos!=1) {
- throw cable_cell_error(util::pprintf(
- "cable_cell does not support regions with partial branches: {}", c));
- }
- assert_valid_segment(c.branch);
- paint_segment(segments[c.branch], prop);
}
}
};
@@ -226,53 +110,17 @@ impl_ptr make_impl(cable_cell_impl* c) {
return impl_ptr(c, [](cable_cell_impl* p){delete p;});
}
-cable_cell::cable_cell(const arb::morphology& m,
- const label_dict& dictionary,
- bool compartments_from_discretization):
- impl_(make_impl(new cable_cell_impl(m, dictionary, compartments_from_discretization)))
+cable_cell::cable_cell(const arb::morphology& m, const label_dict& dictionary):
+ impl_(make_impl(new cable_cell_impl(m, dictionary)))
{}
-cable_cell::cable_cell():
- impl_(make_impl(new cable_cell_impl()))
-{}
+cable_cell::cable_cell(): impl_(make_impl(new cable_cell_impl())) {}
cable_cell::cable_cell(const cable_cell& other):
default_parameters(other.default_parameters),
impl_(make_impl(new cable_cell_impl(*other.impl_)))
{}
-size_type cable_cell::num_branches() const {
- return impl_->segments.size();
-}
-
-segment const* cable_cell::parent(index_type index) const {
- impl_->assert_valid_segment(index);
- return impl_->segments[impl_->parents[index]].get();
-}
-
-segment const* cable_cell::segment(index_type index) const {
- impl_->assert_valid_segment(index);
- return impl_->segments[index].get();
-}
-
-const std::vector<segment_ptr>& cable_cell::segments() const {
- return impl_->segments;
-}
-
-const std::vector<index_type>& cable_cell::parents() const {
- return impl_->parents;
-}
-
-value_type cable_cell::segment_length_constant(value_type frequency, index_type segidx,
- const cable_cell_parameter_set& global_defaults) const
-{
- return 0.5/segment_mean_attenuation(frequency, segidx, global_defaults);
-}
-
-bool cable_cell::has_soma() const {
- return !segment(0)->is_placeholder();
-}
-
const concrete_embedding& cable_cell::embedding() const {
return impl_->provider.embedding();
}
@@ -312,82 +160,4 @@ lid_range cable_cell::place(const locset& target, proptype prop) {\
ARB_PP_FOREACH(FWD_PLACE,\
mechanism_desc, i_clamp, gap_junction_site, threshold_detector)
-//
-// TODO: deprectate the following as soon as discretization code catches up with em_morphology
-//
-const soma_segment* cable_cell::soma() const {
- return has_soma()? segment(0)->as_soma(): nullptr;
-}
-
-const cable_segment* cable_cell::cable(index_type index) const {
- impl_->assert_valid_segment(index);
- auto cable = segment(index)->as_cable();
- return cable? cable: throw cable_cell_error("segment is not a cable segment");
-}
-
-std::vector<size_type> cable_cell::compartment_counts() const {
- std::vector<size_type> comp_count;
- comp_count.reserve(num_branches());
- for (const auto& s: segments()) {
- comp_count.push_back(s->num_compartments());
- }
- return comp_count;
-}
-
-size_type cable_cell::num_compartments() const {
- return util::sum_by(impl_->segments,
- [](const segment_ptr& s) { return s->num_compartments(); });
-}
-
-// Approximating wildly by ignoring O(x) effects entirely, the attenuation b
-// over a single cable segment with constant resistivity R and membrane
-// capacitance C is given by:
-//
-// b = 2√(Ï€RCf) · Σ 2L/(√dâ‚€ + √dâ‚)
-//
-// where the sum is taken over each piecewise linear segment of length L
-// with diameters dâ‚€ and dâ‚ at each end.
-
-value_type cable_cell::segment_mean_attenuation(
- value_type frequency, index_type segidx,
- const cable_cell_parameter_set& global_defaults) const
-{
- value_type R = default_parameters.axial_resistivity.value_or(
- global_defaults.axial_resistivity.value());
- value_type C = default_parameters.membrane_capacitance.value_or(
- global_defaults.membrane_capacitance.value());
-
- value_type length_factor = 0; // [1/õm]
-
- if (segidx==0) {
- if (const soma_segment* s = soma()) {
- R = s->parameters.axial_resistivity.value_or(R);
- C = s->parameters.membrane_capacitance.value_or(C);
-
- value_type d = 2*s->radius();
- length_factor = 1/std::sqrt(d);
- }
- }
- else {
- const cable_segment* s = cable(segidx);
- const auto& lengths = s->lengths();
- const auto& radii = s->radii();
-
- value_type total_length = 0;
- R = s->parameters.axial_resistivity.value_or(R);
- C = s->parameters.membrane_capacitance.value_or(C);
-
- for (std::size_t i = 0; i<lengths.size(); ++i) {
- length_factor += 2*lengths[i]/(std::sqrt(radii[i])+std::sqrt(radii[i+1]));
- total_length += lengths[i];
- }
- length_factor /= total_length;
- }
-
- // R*C is in [s·cm/m²]; need to convert to [s/µm]
- value_type tau_per_um = R*C*1e-8;
-
- return 2*std::sqrt(math::pi<double>*tau_per_um*frequency)*length_factor; // [1/µm]
-}
-
} // namespace arb
diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp
index d1e491ece8534bc0fa314b78950ddfa23a2042a8..6ade3d643516b3239333808cc0e4e7b027615578 100644
--- a/arbor/cable_cell_param.cpp
+++ b/arbor/cable_cell_param.cpp
@@ -139,4 +139,35 @@ locset cv_policy_fixed_per_branch::cv_boundary_points(const cable_cell& cell) co
return points;
}
+locset cv_policy_every_sample::cv_boundary_points(const cable_cell& cell) const {
+ const unsigned nbranch = cell.morphology().num_branches();
+ if (!nbranch) return ls::nil();
+
+ // Ignore interior_forks flag, but if single_root_cv is set, take sample indices only from branches 1+.
+ // Always include branch proximal points, so that forks are trivial.
+
+ std::vector<locset> points;
+
+ if (flags_&cv_policy_flag::single_root_cv) {
+ points.push_back(mlocation{0, 0.});
+ points.push_back(mlocation{0, 1.});
+ for (unsigned i = 1; i<nbranch; ++i) {
+ points.push_back(mlocation{i, 0.});
+ for (msize_t sample_idx: util::make_range(cell.morphology().branch_indexes(i))) {
+ points.push_back(ls::sample(sample_idx));
+ }
+ }
+ }
+ else {
+ for (unsigned i = 0; i<nbranch; ++i) {
+ points.push_back(mlocation{i, 0.});
+ for (msize_t sample_idx: util::make_range(cell.morphology().branch_indexes(i))) {
+ points.push_back(ls::sample(sample_idx));
+ }
+ }
+ }
+
+ return std::accumulate(points.begin(), points.end(), ls::nil(), [](auto& l, auto& p) { return join(l, p); });
+}
+
} // namespace arb
diff --git a/arbor/fvm_compartment.hpp b/arbor/fvm_compartment.hpp
deleted file mode 100644
index 9f5713fc9ffb16bb66afaec9a8291c77e8311d7b..0000000000000000000000000000000000000000
--- a/arbor/fvm_compartment.hpp
+++ /dev/null
@@ -1,261 +0,0 @@
-#pragma once
-
-#include <iterator>
-#include <utility>
-
-#include <arbor/common_types.hpp>
-#include <arbor/math.hpp>
-#include <arbor/util/compat.hpp>
-
-#include "util/iterutil.hpp"
-#include "util/partition.hpp"
-#include "util/rangeutil.hpp"
-#include "util/transform.hpp"
-
-namespace arb {
-
-/// Divided compartments for use with (e.g.) fvm control volume setup
-
-struct semi_compartment {
- using value_type = double;
- using value_pair = std::pair<value_type, value_type>;
-
- value_type length;
- value_type area;
- value_type volume;
- value_pair radii;
-
- semi_compartment& operator+=(const semi_compartment& x) {
- length += x.length;
- area += x.area;
- volume += x.volume;
- radii.second = x.radii.second;
- return *this;
- }
-
- static semi_compartment half_sphere(value_type r1) {
- using namespace math;
- return semi_compartment{
- r1,
- area_sphere(r1)/2,
- volume_sphere(r1)/2,
- {r1, r1}
- };
- }
-
- static semi_compartment frustrum(value_type l, value_type r1, value_type r2) {
- using namespace math;
- return semi_compartment{
- l,
- area_frustrum(l, r1, r2),
- volume_frustrum(l, r1, r2),
- {r1, r2}
- };
- }
-};
-
-struct div_compartment {
- using value_type = typename semi_compartment::value_type;
- using value_pair = typename semi_compartment::value_pair;
- using size_type = cell_local_size_type;
-
- size_type index = 0;
- semi_compartment left;
- semi_compartment right;
-
- div_compartment() = default;
- div_compartment(size_type i, semi_compartment l, semi_compartment r):
- index(i), left(std::move(l)), right(std::move(r))
- {}
-
- value_type length() const { return left.length+right.length; }
- value_type area() const { return left.area+right.area; }
- value_type volume() const { return left.volume+right.volume; }
- value_pair radii() const { return {left.radii.first, right.radii.second}; }
-};
-
-/// Divided compartments can be made from cables with sub-segments by
-/// sampling or integrating or by approximating cable as a single frustrum.
-
-class div_compartment_by_ends {
-public:
- using value_type = div_compartment::value_type;
- using size_type = div_compartment::size_type;
-
- // `lengths` must be a sequence of length at least one.
- // `radii` must be a sequence of length `size(lengths)+1`.
- template <typename Radii, typename Lengths>
- div_compartment_by_ends(size_type n, const Radii& radii, const Lengths& lengths):
- oon_(1/value_type(n)),
- length_(util::sum(lengths)),
- ra_(util::front(radii)),
- rb_(util::back(radii))
- {}
-
- div_compartment operator()(size_type i) const {
- value_type r1 = math::lerp(ra_, rb_, i*oon_);
- value_type rc = math::lerp(ra_, rb_, (i+0.5)*oon_);
- value_type r2 = math::lerp(ra_, rb_, (i+1)*oon_);
- value_type semilength = length_*oon_*0.5;
-
- return div_compartment(
- i,
- semi_compartment::frustrum(semilength, r1, rc),
- semi_compartment::frustrum(semilength, rc, r2)
- );
- }
-
-private:
- value_type oon_;
- value_type length_;
- value_type ra_;
- value_type rb_;
-};
-
-class div_compartment_sampler {
-public:
- using value_type = div_compartment::value_type;
- using size_type = div_compartment::size_type;
-
- // `lengths` must be a sequence of length at least one.
- // `radii` must be a sequence of length `size(lengths)+1`.
- template <typename Radii, typename Lengths>
- div_compartment_sampler(size_type n, const Radii& radii, const Lengths& lengths, bool with_soma = false) {
- // set up offset-to-subsegment lookup and interpolation
- using namespace util;
-
- with_soma_ = with_soma;
- segs_ = make_partition(offsets_, lengths);
- compat::compiler_barrier_if_icc_leq(20160415u);
-
- nseg_ = size(segs_);
- scale_ = with_soma_? (segs_.bounds().second - segs_.front().second)/(n-1): segs_.bounds().second/n;
-
- assign(radii_, radii);
- arb_assert(size(radii_)==size(offsets_));
- }
-
- div_compartment operator()(size_type i) const {
- using namespace math;
-
- auto r1 = radius_at(locate(scale_*i));
- auto rc = radius_at(locate(scale_*(i+0.5)));
- auto r2 = radius_at(locate(scale_*(i+1)));
-
- value_type semilength = 0.5*scale_;
- return div_compartment(
- i,
- semi_compartment::frustrum(semilength, r1, rc),
- semi_compartment::frustrum(semilength, rc, r2)
- );
- }
-
-protected:
- struct sub_segment_index {
- size_type i; // index
- value_type p; // proportion [0,1] along sub-segment
-
- sub_segment_index(): i(0), p(0) {}
-
- sub_segment_index(size_type i_, value_type p_): i(i_), p(p_) {}
- bool operator<(sub_segment_index x) const {
- return i<x.i || (i==x.i && p<x.p);
- }
- };
-
- sub_segment_index locate(value_type x) const {
- arb_assert(x>=0);
-
- auto i = segs_.index(x);
- if (i==segs_.npos) {
- i = nseg_-1;
- }
-
- auto seg = segs_[i];
- if (x>=seg.second) {
- return sub_segment_index(i, 1);
- }
- return sub_segment_index(i, (x-seg.first)/(seg.second-seg.first));
- }
-
- value_type radius_at(sub_segment_index s) const {
- return math::lerp(radii_[s.i], radii_[s.i+1], s.p);
- }
-
- bool with_soma_ = false; //segment passed in includes the soma information in radii[0] and lengths[0]
- size_type nseg_ = 0;
- std::vector<value_type> offsets_;
- std::vector<value_type> radii_;
- value_type scale_ = 0;
- decltype(util::partition_view(offsets_)) segs_;
-};
-
-/// Overrides operator() with a more accurate method
-class div_compartment_integrator: public div_compartment_sampler {
-public:
- template <typename Radii, typename Lengths>
- div_compartment_integrator(size_type n, const Radii& radii, const Lengths& lengths, bool with_soma = false):
- div_compartment_sampler(n, radii, lengths, with_soma) {}
-
- div_compartment operator()(size_type i) const {
- using namespace math;
- sub_segment_index sleft, scentre, sright;
- // If segment includes the soma, divisions are specially calculated
- if (with_soma_) {
- auto soma_l = segs_.front().second;
- if (i == 0) {
- sleft = locate(soma_l/2);
- scentre = locate(soma_l);
- sright = locate(soma_l + scale_/4);
- } else if (i == 1) {
- sleft = locate(soma_l + scale_/4);
- scentre = locate(soma_l + scale_/2);
- sright = locate(soma_l + scale_);
- } else {
- sleft = locate(soma_l + scale_ * (i-1));
- scentre = locate(soma_l + scale_ * (i-0.5));
- sright = locate(soma_l + scale_ * i);
- }
- } else {
- sleft = locate(scale_ * i);
- scentre = locate(scale_ * (i + 0.5));
- sright = locate(scale_ * (i + 1));
- }
-
- return div_compartment(i, integrate(sleft, scentre), integrate(scentre, sright));
- }
-
-protected:
- semi_compartment sub_segment_frustrum(sub_segment_index a, sub_segment_index b) const {
- arb_assert(a.i==b.i && a.p<=b.p);
-
- auto seg = segs_[a.i];
- auto l = (b.p-a.p)*(seg.second-seg.first);
-
- // If segment includes the soma it will be at index 0
- // The soma is represented as a sphere
- if (with_soma_ && a.i==0) {
- return semi_compartment::half_sphere(radii_.front());
- }
- return semi_compartment::frustrum(l, radius_at(a), radius_at(b));
- }
-
- semi_compartment integrate(sub_segment_index a, sub_segment_index b) const {
- sub_segment_index x = std::min(b, sub_segment_index(a.i, 1));
-
- auto s = sub_segment_frustrum(a, x);
-
- while (a.i<b.i) {
- ++a.i;
- a.p = 0;
- x = std::min(b, sub_segment_index(a.i, 1));
- s += sub_segment_frustrum(a, x);
- }
-
- return s;
- }
-};
-
-} // namespace arb
-
-
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index deb1d7806c7bfa9f4edcf355cb0df59009b90e54..56068d6510d32144d8d4a4827d69a7035619e5c1 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -5,13 +5,13 @@
#include <arbor/arbexcept.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/math.hpp>
#include <arbor/morph/mcable_map.hpp>
+#include <arbor/morph/mprovider.hpp>
+#include <arbor/morph/morphology.hpp>
#include <arbor/util/optional.hpp>
-#include "algorithms.hpp"
-#include "fvm_compartment.hpp"
#include "fvm_layout.hpp"
-#include "tree.hpp"
#include "util/maputil.hpp"
#include "util/meta.hpp"
#include "util/partition.hpp"
@@ -21,1022 +21,903 @@
namespace arb {
+using util::append;
+using util::assign;
+using util::assign_by;
using util::count_along;
-using util::keys;
-using util::make_span;
-using util::optional;
-using util::subrange_view;
-using util::transform_view;
+using util::pw_elements;
+using util::pw_element;
+using util::sort;
+using util::sort_by;
using util::value_by_key;
-// Convenience routines
-
namespace {
- template <typename ResizableContainer, typename Index>
- void extend_to(ResizableContainer& c, const Index& i) {
- if (util::size(c)<=i) {
- c.resize(i+1);
+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>
+util::optional<V> operator|(const util::optional<V>& a, const util::optional<V>& b) {
+ return a? a: b;
+}
+
+// Given sorted vectors a, b, return sorted vector with unique elements v
+// such that v is present in a or b.
+template <typename V>
+std::vector<V> unique_union(const std::vector<V>& a, const std::vector<V>& b) {
+ std::vector<V> u;
+
+ auto ai = a.begin();
+ auto ae = a.end();
+ auto bi = b.begin();
+ auto be = b.end();
+
+ while (ai!=ae && bi!=be) {
+ const V& elem = *ai<*bi? *ai++: *bi++;
+ if (u.empty() || u.back()!=elem) {
+ u.push_back(elem);
}
}
- template <typename ResizableContainer, typename Index>
- auto& extend_at(ResizableContainer& c, const Index& i) {
- if (util::size(c)<=i) {
- c.resize(i+1);
+ while (ai!=ae) {
+ const V& elem = *ai++;
+ if (u.empty() || u.back()!=elem) {
+ u.push_back(elem);
}
- return c[i];
}
- struct compartment_model {
- arb::tree tree;
- std::vector<tree::int_type> parent_index;
- std::vector<tree::int_type> segment_index;
-
- explicit compartment_model(const cable_cell& cell) {
- tree = arb::tree(cell.parents());
- auto counts = cell.compartment_counts();
- for (unsigned i = 0; i < cell.segments().size(); i++) {
- if (!cell.segment(i)->is_soma() && cell.parent(i)->is_soma()) {
- counts[i]++;
- }
- }
- parent_index = make_parent_index(tree, counts);
- segment_index = algorithms::make_index(counts);
+ while (bi!=be) {
+ const V& elem = *bi++;
+ if (u.empty() || u.back()!=elem) {
+ u.push_back(elem);
}
+ }
+ return u;
+}
+
+} // anonymous namespace
+
+// 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) {}
};
- struct cv_param {
- fvm_size_type cv;
- std::vector<fvm_value_type> params;
- fvm_size_type target;
+ 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 = zip(pw, pw_elements<>({cable.prox_pos, cable.dist_pos}));
+ }
+ return pw;
+}
+
+// Construct cv_geometry for cell from locset describing CV boundary points.
+cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
+ struct mloc_hash {
+ std::size_t operator()(const mlocation& loc) const { return loc.branch ^ std::hash<double>()(loc.pos); }
};
+ using mlocation_map = std::unordered_map<mlocation, fvm_size_type, mloc_hash>;
+ auto pop = [](auto& vec) { auto h = vec.back(); return vec.pop_back(), h; };
- ARB_DEFINE_LEXICOGRAPHIC_ORDERING(cv_param,(a.cv,a.params,a.target),(b.cv,b.params,b.target))
+ cv_geometry geom;
+ const auto& mp = cell.provider();
+ const auto& m = mp.morphology();
- template <typename V>
- optional<V> operator|(const optional<V>& a, const optional<V>& b) {
- return a? a: b;
+ if (mp.morphology().empty()) {
+ geom.cell_cv_divs = {0, 0};
+ return geom;
}
- // For each segment given by the provided sorted sequence of segment
- // indices, call `action` for each CV intersecting the segment, starting
- // from the most proximal.
- //
- // By the ordering of CVs and segments in the discretization, with the
- // exception of the most proximal CV in a segment, each CV will be visited
- // once, and the visited CVs will be in increasing order. The most proximal
- // CV (the 'parent' CV) may be visited multiple times.
- //
- // Action is a functional that takes the following arguments:
- //
- // size_type cv_index The index into the total (sorted) list of
- // CVs that constitute the segments.
- //
- // index_type cv The CV number (within the discretization).
- //
- // value_type cv_area The area of the CV that lies within the
- // current segment.
- //
- // index_type seg_index The index into the provided sequence of
- // the provided segment_indices.
- //
- // index_type seg The segment currently being iterated over.
-
- template <typename Seq, typename Action>
- void for_each_cv_in_segments(const fvm_discretization& D, const Seq& segment_indices, const Action& action) {
- using index_type = fvm_index_type;
- using size_type = fvm_size_type;
-
- std::unordered_map<index_type, size_type> parent_cv_indices;
- size_type cv_index = 0;
-
- index_type seg_index = 0;
- for (const auto& seg: segment_indices) {
- const segment_info& seg_info = D.segments[seg];
-
- if (seg_info.has_parent()) {
- index_type cv = seg_info.parent_cv;
-
- size_type i = parent_cv_indices.insert({cv, cv_index}).first->second;
- if (i==cv_index) {
- ++cv_index;
- }
+ auto canon = [&m](mlocation loc) { return canonical(m, loc); };
+ auto origin = [&canon](mcable cab) { return canon(mlocation{cab.branch, cab.prox_pos}); };
+ auto terminus = [&canon](mcable cab) { return canon(mlocation{cab.branch, cab.dist_pos}); };
+
+ mlocation_list locs = thingify(lset, mp);
+
+ std::vector<msize_t> n_cv_cables;
+ std::vector<mlocation> next_cv_head;
+ next_cv_head.push_back({mnpos, 0});
+
+ mcable_list cables, all_cables;
+ mlocation_list ends;
+ std::vector<msize_t> branches;
+
+ mlocation_map head_count;
+ unsigned extra_cv_count = 0;
+
+ while (!next_cv_head.empty()) {
+ mlocation h = pop(next_cv_head);
- action(i, cv, seg_info.parent_cv_area, seg_index, seg);
+ cables.clear();
+ branches.clear();
+ branches.push_back(h.branch);
+
+ 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});
}
- for (index_type cv = seg_info.proximal_cv; cv < seg_info.distal_cv; ++cv) {
- index_type i = cv_index++;
- action(i, cv, D.cv_area[cv], seg_index, seg);
+ // 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.push_back(*it);
+ }
+ 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);
+ }
}
+ }
- index_type cv = seg_info.distal_cv;
- size_type i = cv_index++;
+ auto empty_cable = [](mcable c) { return c.prox_pos==c.dist_pos; };
+ cables.erase(std::remove_if(cables.begin(), cables.end(), empty_cable), cables.end());
- action(i, seg_info.distal_cv, seg_info.distal_cv_area, seg_index, seg);
- parent_cv_indices.insert({cv, i});
- ++seg_index;
+ if (!cables.empty()) {
+ if (++head_count[origin(cables.front())]==2) {
+ ++extra_cv_count;
+ }
+
+ n_cv_cables.push_back(cables.size());
+ sort(cables);
+ append(all_cables, std::move(cables));
}
}
-} // namespace
-
-// Cable segment discretization
-// ----------------------------
-//
-// Each compartment i straddles the ith control volume on the right
-// and the jth control volume on the left, where j is the parent index
-// of i.
-//
-// Dividing the comparment into two halves, the centre face C
-// corresponds to the shared face between the two control volumes,
-// the surface areas in each half contribute to the surface area of
-// the respective control volumes, and the volumes and lengths of
-// each half are used to calculate the flux coefficients that
-// for the connection between the two control volumes and which
-// is stored in `face_conductance[i]`.
-//
-//
-// +------- cv j --------+------- cv i -------+
-// | | |
-// v v v
-// ____________________________________________
-// | ........ | ........ | | |
-// | ........ L ........ C R |
-// |__________|__________|__________|_________|
-// ^ ^
-// | |
-// +--- compartment i ---+
-//
-// The first control volume of any cell corresponds to the soma
-// and the first half of the first cable compartment of that cell.
-//
-//
-// Face conductance computation
-// ----------------------------
-//
-// The conductance between two adjacent CVs is computed as follows,
-// computed in terms of the two half CVs on either side of the interface,
-// correspond to the regions L–C and C–R in the diagram above.
-//
-// The conductance itself is approximated by the weighted harmonic mean
-// of the mean linear conductivities in each half, corresponding to
-// the two-point flux approximation in 1-D.
-//
-// Mean linear conductivities:
-//
-// g₠= 1/h₠∫₠A(x)/R dx
-// g₂ = 1/h₂ ∫₂ A(x)/R dx
-//
-// where A(x) is the cross-sectional area, R is the bulk resistivity,
-// and h is the width of the region. The integrals are taken over the
-// half CVs as described above.
-//
-// Equivalently, in terms of the semi-compartment volumes Vâ‚ and Vâ‚‚:
-//
-// gâ‚ = 1/R·Vâ‚/hâ‚
-// g₂ = 1/R·V₂/h₂
-//
-// Weighted harmonic mean, with h = hâ‚+hâ‚‚:
-//
-// g = (hâ‚/h·g₯¹+hâ‚‚/h·g₂¯¹)¯¹
-// = 1/R · hVâ‚Vâ‚‚/(h₂²Vâ‚+h₲Vâ‚‚)
-//
-
-fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells, const cable_cell_parameter_set& global_defaults) {
+ geom.cv_cables.reserve(all_cables.size()+extra_cv_count);
+ geom.cv_parent.reserve(n_cv_cables.size()+extra_cv_count);
+ geom.cv_cables_divs.reserve(n_cv_cables.size()+extra_cv_count+1);
+ geom.cv_cables_divs.push_back(0);
+
+ mlocation_map parent_map;
+
+ unsigned all_cables_index = 0;
+ unsigned cv_index = 0;
+
+ // Multiple CVs meeting at (0,0)?
+ mlocation root{0, 0};
+ unsigned n_top_children = value_by_key(head_count, root).value_or(0);
+ if (n_top_children>1) {
+ // Add initial trical CV.
+ geom.cv_parent.push_back(mnpos);
+ geom.cv_cables.push_back(mcable{0, 0, 0});
+ geom.cv_cables_divs.push_back(geom.cv_cables.size());
+ parent_map[root] = cv_index++;
+ }
- using value_type = fvm_value_type;
- using index_type = fvm_index_type;
- using size_type = fvm_size_type;
+ for (auto n_cables: n_cv_cables) {
+ mlocation head = origin(all_cables[all_cables_index]);
+ msize_t parent_cv = value_by_key(parent_map, head).value_or(mnpos);
- fvm_discretization D;
+ auto cables = util::subrange_view(all_cables, all_cables_index, all_cables_index+n_cables);
+ std::copy(cables.begin(), cables.end(), std::back_inserter(geom.cv_cables));
- util::make_partition(D.cell_segment_bounds,
- transform_view(cells, [](const cable_cell& c) { return c.num_branches(); }));
+ geom.cv_parent.push_back(parent_cv);
+ geom.cv_cables_divs.push_back(geom.cv_cables.size());
- std::vector<index_type> cell_cv_bounds;
- auto cell_cv_part = make_partition(cell_cv_bounds,
- transform_view(cells, [](const cable_cell& c) {
- unsigned ncv = 0;
- for (unsigned i = 0; i < c.segments().size(); i++) {
- ncv += c.segment(i)->num_compartments();
- if (!c.segment(i)->is_soma() && c.parent(i)->is_soma()) {
- ncv++;
- }
+ auto this_cv = cv_index++;
+ for (auto cable: cables) {
+ mlocation term = terminus(cable);
+
+ unsigned n_children = value_by_key(head_count, term).value_or(0);
+ if (n_children>1) {
+ // Add trivial CV for lindep.
+ geom.cv_parent.push_back(this_cv);
+ geom.cv_cables.push_back(mcable{cable.branch, 1., 1.});
+ geom.cv_cables_divs.push_back((fvm_index_type)geom.cv_cables.size());
+ parent_map[term] = cv_index++;
+ }
+ else {
+ parent_map[term] = this_cv;
}
- return ncv;
- }));
-
- D.ncell = cells.size();
- D.ncv = cell_cv_part.bounds().second;
-
- D.face_conductance.assign(D.ncv, 0.);
- D.cv_area.assign(D.ncv, 0.);
- D.cv_capacitance.assign(D.ncv, 0.);
- D.init_membrane_potential.assign(D.ncv, 0.);
- D.temperature_K.assign(D.ncv, 0.);
- D.diam_um.assign(D.ncv, 0.);
- D.parent_cv.assign(D.ncv, index_type(-1));
- D.cv_to_cell.resize(D.ncv);
- for (auto i: make_span(0, D.ncell)) {
- util::fill(subrange_view(D.cv_to_cell, cell_cv_part[i]), static_cast<index_type>(i));
+ }
+
+ all_cables_index += n_cables;
}
- std::vector<size_type> seg_cv_bounds;
- for (auto i: make_span(0, D.ncell)) {
- const auto& c = cells[i];
- compartment_model cell_graph(c);
- auto cell_cv_ival = cell_cv_part[i];
+ // 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};
- auto cell_cv_base = cell_cv_ival.first;
- for (auto k: make_span(cell_cv_ival)) {
- D.parent_cv[k] = cell_graph.parent_index[k-cell_cv_base]+cell_cv_base;
- }
+ // Build location query map.
+ geom.branch_cv_map.resize(1);
+ std::vector<pw_elements<fvm_size_type>>& bmap = geom.branch_cv_map.back();
- // Electrical defaults from global defaults, possibly overridden by cell.
- auto cm_default = c.default_parameters.membrane_capacitance | global_defaults.membrane_capacitance;
- auto rL_default = c.default_parameters.axial_resistivity | global_defaults.axial_resistivity;
- auto init_vm_default = c.default_parameters.init_membrane_potential | global_defaults.init_membrane_potential;
- auto temp_default = c.default_parameters.temperature_K | global_defaults.temperature_K;
-
- // Compartment index range for each segment in this cell.
- seg_cv_bounds.clear();
- auto seg_cv_part = make_partition(
- seg_cv_bounds,
- transform_view(make_span(c.num_branches()), [&c](const unsigned s) {
- if (!c.segment(s)->is_soma() && c.parent(s)->is_soma()) {
- return c.segment(s)->num_compartments() + 1;
- }
- return c.segment(s)->num_compartments();
- }),
- cell_cv_base);
+ for (auto cv: util::make_span(geom.size())) {
+ for (auto cable: geom.cables(cv)) {
+ if (cable.branch>=bmap.size()) {
+ bmap.resize(cable.branch+1);
+ }
- const auto nseg = seg_cv_part.size();
- if (nseg==0) {
- throw arbor_internal_error("fvm_layout: cannot discretrize cell with no segments");
+ // Ordering of CV ensures CV cables on any given branch are found sequentially.
+ // Omit empty cables.
+ if (cable.prox_pos<cable.dist_pos) {
+ bmap[cable.branch].push_back(cable.prox_pos, cable.dist_pos, cv);
+ }
}
+ }
- // Handle soma (first segment and root of tree) specifically.
- const auto soma = c.segment(0)->as_soma();
- if (!soma) {
- throw arbor_internal_error("fvm_layout: first segment of cell must be soma");
- }
- else if (soma->num_compartments()!=1) {
- throw arbor_internal_error("fvm_layout: soma must have exactly one compartment");
- }
+ return geom;
+}
- segment_info soma_info;
+namespace impl {
+ using std::begin;
+ using std::end;
+ using std::next;
- size_type soma_cv = cell_cv_base;
- value_type soma_area = math::area_sphere(soma->radius());
+ template <typename Seq>
+ auto tail(Seq& seq) { return util::make_range(next(begin(seq)), end(seq)); };
- soma_info.proximal_cv = soma_cv;
- soma_info.distal_cv = soma_cv;
- soma_info.distal_cv_area = soma_area;
- D.segments.push_back(soma_info);
+ template <typename Container, typename Offset, typename Seq>
+ void append_offset(Container& ctr, Offset offset, const Seq& rhs) {
+ for (const auto& x: rhs) {
+ // Preserve -1 'npos' values.
+ ctr.push_back(x+1==0? x: offset+x);
+ }
+ }
+}
- index_type soma_segment_index = D.segments.size()-1;
- D.parent_segment.push_back(soma_segment_index);
+// Merge CV geometry lists in-place.
+cv_geometry& append(cv_geometry& geom, const cv_geometry& right) {
+ using impl::tail;
+ using impl::append_offset;
- // Other segments must all be cable segments.
- for (size_type j = 1; j<nseg; ++j) {
- const auto& seg_cv_ival = seg_cv_part[j];
- const auto ncv = seg_cv_ival.second-seg_cv_ival.first;
+ if (!right.n_cell()) {
+ return geom;
+ }
- segment_info seg_info;
+ if (!geom.n_cell()) {
+ geom = right;
+ return geom;
+ }
- const auto cable = c.segment(j)->as_cable();
- if (!cable) {
- throw arbor_internal_error("fvm_layout: non-root segments of cell must be cable segments");
- }
+ auto append_divs = [](auto& left, const auto& right) {
+ if (left.empty()) {
+ left = right;
+ }
+ else if (!right.empty()) {
+ append_offset(left, left.back(), tail(right));
+ }
+ };
- const auto& params = cable->parameters;
- auto cm = (params.membrane_capacitance | cm_default).value(); // [F/m²]
- auto rL = (params.axial_resistivity | rL_default).value(); // [Ω·cm]
- auto init_vm = (params.init_membrane_potential | init_vm_default).value(); // [mV]
- auto temp = (params.temperature_K | temp_default).value(); // [mV]
+ auto geom_n_cv = geom.size();
+ auto geom_n_cell = geom.n_cell();
- bool soma_parent = c.parent(j)->as_soma() ? true : false; //segment's parent is a soma
+ append(geom.cv_cables, right.cv_cables);
+ append_divs(geom.cv_cables_divs, right.cv_cables_divs);
- auto radii = cable->radii();
- auto lengths = cable->lengths();
+ append_offset(geom.cv_parent, geom_n_cv, right.cv_parent);
+ append_offset(geom.cv_to_cell, geom_n_cell, right.cv_to_cell);
+ append_divs(geom.cell_cv_divs, right.cell_cv_divs);
- // If segment has soma parent, send soma information to div_compartment_integrator
- if (soma_parent) {
- radii.insert(radii.begin(), soma->radius());
- lengths.insert(lengths.begin(), soma->radius()*2);
- }
+ append(geom.branch_cv_map, right.branch_cv_map);
+ return geom;
+}
- auto divs = div_compartment_integrator(ncv, radii, lengths, soma_parent);
+// Combine two fvm_cv_geometry groups in-place.
+fvm_cv_discretization& append(fvm_cv_discretization& dczn, const fvm_cv_discretization& right) {
+ append(dczn.geometry, right.geometry);
- seg_info.parent_cv = soma_parent ? seg_cv_ival.first : D.parent_cv[seg_cv_ival.first];
- seg_info.parent_cv_area = soma_parent ? divs(0).right.area + divs(1).left.area : divs(0).left.area;
- seg_info.soma_parent = soma_parent;
+ append(dczn.face_conductance, right.face_conductance);
+ append(dczn.cv_area, right.cv_area);
+ append(dczn.cv_capacitance, right.cv_capacitance);
+ append(dczn.init_membrane_potential, right.init_membrane_potential);
+ append(dczn.temperature_K, right.temperature_K);
+ append(dczn.diam_um, right.diam_um);
- seg_info.proximal_cv = soma_parent ? seg_cv_ival.first + 1 : seg_cv_ival.first;
- seg_info.distal_cv = seg_cv_ival.second-1;
- seg_info.distal_cv_area = divs(ncv-1).right.area;
+ return dczn;
+}
- D.segments.push_back(seg_info);
- if (soma_parent) {
- D.parent_segment.push_back(soma_segment_index);
+fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global_dflt) {
+ const auto& dflt = cell.default_parameters;
+ fvm_cv_discretization D;
+
+ cv_policy pol = (dflt.discretization | global_dflt.discretization).value_or(default_cv_policy());
+ locset cv_ends = pol.cv_boundary_points(cell);
+ D.geometry = cv_geometry_from_ends(cell, cv_ends);
+
+ if (D.geometry.empty()) return D;
+
+ auto n_cv = D.geometry.size();
+ D.face_conductance.resize(n_cv);
+ D.cv_area.resize(n_cv);
+ D.cv_capacitance.resize(n_cv);
+ D.init_membrane_potential.resize(n_cv);
+ D.temperature_K.resize(n_cv);
+ D.diam_um.resize(n_cv);
+
+ double dflt_resistivity = *(dflt.axial_resistivity | global_dflt.axial_resistivity);
+ double dflt_capacitance = *(dflt.membrane_capacitance | global_dflt.membrane_capacitance);
+ double dflt_potential = *(dflt.init_membrane_potential | global_dflt.init_membrane_potential);
+ double dflt_temperature = *(dflt.temperature_K | global_dflt.temperature_K);
+
+ const auto& embedding = cell.embedding();
+ for (auto i: count_along(D.geometry.cv_parent)) {
+ auto cv_cables = D.geometry.cables(i);
+
+ // Computing face_conductance:
+ //
+ // Flux between adjacemt CVs is computed as if there were no membrane currents, and with the CV voltage
+ // values taken to be exact at a reference point in each CV:
+ // * If the CV is unbranched, the reference point is taken to be the CV midpoint.
+ // * If the CV is branched, the reference point is taken to be closest branch point to
+ // the interface between the two CVs.
+
+ D.face_conductance[i] = 0;
+
+ fvm_index_type p = D.geometry.cv_parent[i];
+ if (p!=-1) {
+ auto parent_cables = D.geometry.cables(p);
+
+ msize_t bid = cv_cables.front().branch;
+ double parent_refpt = 0;
+ double cv_refpt = 1;
+
+ if (cv_cables.size()==1) {
+ mcable cv_cable = cv_cables.front();
+ cv_refpt = 0.5*(cv_cable.prox_pos+cv_cable.dist_pos);
}
- else {
- auto opt_index = util::binary_search_index(D.segments, seg_info.parent_cv,
- [](const segment_info& seg_info) { return seg_info.distal_cv; });
-
- if (!opt_index) {
- throw arbor_internal_error("fvm_layout: could not find parent segment");
+ if (parent_cables.size()==1) {
+ mcable parent_cable = parent_cables.front();
+ // A trivial parent CV with a zero-length cable might not
+ // be on the same branch.
+ if (parent_cable.branch==bid) {
+ parent_refpt = 0.5*(parent_cable.prox_pos+parent_cable.dist_pos);
}
- D.parent_segment.push_back(*opt_index);
}
+ mcable span{bid, parent_refpt, cv_refpt};
+ double resistance = embedding.integrate_ixa(bid,
+ pw_over_cable(cell.region_assignments().get<axial_resistivity>(), span, dflt_resistivity));
+ D.face_conductance[i] = 100/resistance; // 100 scales to µS.
+ }
- for (auto i: make_span(seg_cv_ival)) {
- const auto& div = divs(i-seg_cv_ival.first);
- auto j = D.parent_cv[i];
+ D.cv_area[i] = 0;
+ D.cv_capacitance[i] = 0;
+ D.init_membrane_potential[i] = 0;
+ D.diam_um[i] = 0;
+ double cv_length = 0;
- auto h1 = div.left.length; // [µm]
- auto V1 = div.left.volume; // [µm³]
- auto h2 = div.right.length; // [µm]
- auto V2 = div.right.volume; // [µm³]
- auto h = h1+h2;
+ for (mcable c: cv_cables) {
+ D.cv_area[i] += embedding.integrate_area(c);
- auto linear_conductivity = 1/rL*h*V1*V2/(h2*h2*V1+h1*h1*V2); // [S·cm¯¹·µm²] ≡ [10²·µS·µm]
- constexpr double unit_scale = 1e2;
- D.face_conductance[i] = unit_scale * linear_conductivity / h; // [µS]
+ D.cv_capacitance[i] += embedding.integrate_area(c.branch,
+ pw_over_cable(cell.region_assignments().get<membrane_capacitance>(), c, dflt_capacitance));
- auto al = div.left.area; // [µm²]
- auto ar = div.right.area; // [µm²]
- auto dr = div.right.radii.second*2; // [µm]
+ D.init_membrane_potential[i] += embedding.integrate_area(c.branch,
+ pw_over_cable(cell.region_assignments().get<init_membrane_potential>(), c, dflt_potential));
- D.cv_area[j] += al; // [µm²]
- D.cv_capacitance[j] += al*cm; // [pF]
- D.init_membrane_potential[j] += al*init_vm; // [mV·µm²]
- D.temperature_K[j] += al*temp; // [K·µm²]
+ D.temperature_K[i] += embedding.integrate_area(c.branch,
+ pw_over_cable(cell.region_assignments().get<temperature_K>(), c, dflt_temperature));
- D.cv_area[i] += ar; // [µm²]
- D.cv_capacitance[i] += ar*cm; // [pF]
- D.init_membrane_potential[i] += ar*init_vm; // [mV·µm²]
- D.temperature_K[i] += ar*temp; // [K·µm²]
- D.diam_um[i] = dr; // [µm]
- }
+ cv_length += embedding.integrate_length(c);
}
- auto soma_cm = (soma->parameters.membrane_capacitance | cm_default).value(); // [F/m²]
- auto soma_init_vm = (soma->parameters.init_membrane_potential | init_vm_default).value(); // [mV]
- auto soma_temp = (soma->parameters.temperature_K | temp_default).value(); // [mV]
-
- D.cv_area[soma_cv] = soma_area; // [µm²]
- D.cv_capacitance[soma_cv] = soma_area*soma_cm; // [pF]
- D.init_membrane_potential[soma_cv] = soma_area*soma_init_vm; // [mV·µm²]
- D.temperature_K[soma_cv] = soma_area*soma_temp; // [K·µm²]
- D.diam_um[soma_cv] = soma->radius()*2; // [µm]
- }
-
- // Rescale CV init_vm and temperature values to get area-weighted means.
- for (auto i: make_span(0, D.ncv)) {
- if (D.cv_area[i]) {
- D.init_membrane_potential[i] /= D.cv_area[i]; // [mV]
- D.temperature_K[i] /= D.cv_area[i]; // [mV]
+ if (D.cv_area[i]>0) {
+ D.init_membrane_potential[i] /= D.cv_area[i];
+ D.temperature_K[i] /= D.cv_area[i];
+ }
+ if (cv_length>0) {
+ D.diam_um[i] = D.cv_area[i]/(cv_length*math::pi<double>);
}
}
- // Number of CVs per cell is exactly number of compartments.
- D.cell_cv_bounds = std::move(cell_cv_bounds);
return D;
}
-// Build up mechanisms.
-//
-// Processing procedes in the following stages:
-//
-// I. Collect segment mechanism info from the cell descriptions into temporary
-// data structures for density mechanism, point mechanisms, and ion channels.
-//
-// II. Build mechanism and ion configuration in `fvm_mechanism_data`:
-// IIa. Ion channel CVs.
-// IIb. Density mechanism CVs, parameter values; ion channel default concentration contributions.
-// IIc. Point mechanism CVs, parameter values, and targets.
-
-fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop, const std::vector<cable_cell>& cells, const fvm_discretization& D) {
- using util::assign;
- using util::sort_by;
- using util::optional;
-
- using value_type = fvm_value_type;
- using index_type = fvm_index_type;
- using size_type = fvm_size_type;
-
- using string_set = std::unordered_set<std::string>;
- using string_index_map = std::unordered_map<std::string, size_type>;
+fvm_cv_discretization fvm_cv_discretize(const std::vector<cable_cell>& cells,
+ const cable_cell_parameter_set& global_defaults)
+{
+ fvm_cv_discretization combined;
- const mechanism_catalogue& catalogue = *gprop.catalogue;
- const cable_cell_parameter_set& gparam = gprop.default_parameters;
+ for (const auto& c: cells) {
+ append(combined, fvm_cv_discretize(c, global_defaults));
+ }
+ return combined;
+}
- fvm_mechanism_data mechdata;
+// CVs are absolute (taken from combined discretization) so do not need to be shifted.
+// Only target numbers need to be shifted.
+fvm_mechanism_data& append(fvm_mechanism_data& left, const fvm_mechanism_data& right) {
+ using impl::append_offset;
- // I. Collect segment mechanism info from cells.
+ fvm_size_type target_offset = left.n_target;
- // Temporary table for density mechanism info, mapping mechanism name to tuple of:
- // 1. Vector of segment indices and mechanism parameter settings where mechanism occurs.
- // 2. Set of the names of parameters that are anywhere modified.
- // 3. Pointer to mechanism metadata from catalogue.
+ for (const auto& kv: right.ions) {
+ fvm_ion_config& L = left.ions[kv.first];
+ const fvm_ion_config& R = kv.second;
- struct density_mech_data {
- std::vector<std::pair<size_type, const mechanism_desc*>> segments;
- string_set paramset;
- std::shared_ptr<mechanism_info> info = nullptr;
- };
- std::unordered_map<std::string, density_mech_data> density_mech_table;
+ append(L.cv, R.cv);
+ append(L.init_iconc, R.init_iconc);
+ append(L.init_econc, R.init_econc);
+ append(L.reset_iconc, R.reset_iconc);
+ append(L.reset_econc, R.reset_econc);
+ append(L.init_revpot, R.init_revpot);
+ }
- // Temporary table for revpot mechanism info, mapping mechanism name to tuple of:
- // 1. Sorted map from cell index to mechanism parameter settings.
- // 2. Set of the names of parameters that are anywhere modified.
- // 3. Pointer to mechanism metadat from catalogue.
+ for (const auto& kv: right.mechanisms) {
+ if (!left.mechanisms.count(kv.first)) {
+ fvm_mechanism_config& L = left.mechanisms[kv.first];
- struct revpot_mech_data {
- std::map<size_type, const mechanism_desc*> cells;
- string_set paramset;
- std::shared_ptr<mechanism_info> info = nullptr;
- };
- std::unordered_map<std::string, revpot_mech_data> revpot_mech_table;
-
- // Temporary table for point mechanism info, mapping mechanism name to tuple:
- // 1. Vector of point info: CV, index into cell group targets, parameter settings.
- // 2. Set of the names of parameters that are anywhere modified.
- // 3. Mechanism parameter settings.
-
- struct point_mech_data {
- struct point_data {
- size_type cv;
- size_type target_index;
- const mechanism_desc* desc;
- };
- std::vector<point_data> points;
- string_set paramset;
- std::shared_ptr<mechanism_info> info = nullptr;
- };
- std::unordered_map<std::string, point_mech_data> point_mech_table;
+ L = kv.second;
+ for (auto& t: L.target) t += target_offset;
+ }
+ else {
+ fvm_mechanism_config& L = left.mechanisms[kv.first];
+ const fvm_mechanism_config& R = kv.second;
+
+ L.kind = R.kind;
+ append(L.cv, R.cv);
+ append(L.multiplicity, R.multiplicity);
+ append(L.norm_area, R.norm_area);
+ append_offset(L.target, target_offset, R.target);
+
+ arb_assert(util::is_sorted_by(L.param_values, util::first));
+ arb_assert(util::is_sorted_by(R.param_values, util::first));
+ arb_assert(L.param_values.size()==R.param_values.size());
+
+ for (auto j: count_along(R.param_values)) {
+ arb_assert(L.param_values[j].first==R.param_values[j].first);
+ append(L.param_values[j].second, R.param_values[j].second);
+ }
+ }
+ }
+ left.n_target += right.n_target;
- // Built-in stimulus mechanism data is dealt with especially below.
- // Record for each stimulus the CV and clamp data.
+ return left;
+}
- std::vector<std::pair<size_type, i_clamp>> stimuli;
+fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop,
+ const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx);
- // Temporary table for presence of ion channels, mapping ion name to a _sorted_
- // collection of per-segment ion data, viz. a map from segment index to
- // initial internal and external concentrations and reversal potentials, plus
- // information about whether there is a mechanism that requires this ion
- // on this segment, and if that ion writes to internal or external concentration.
+fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop,
+ const std::vector<cable_cell>& cells, const fvm_cv_discretization& D)
+{
+ fvm_mechanism_data combined;
+ for (auto cell_idx: count_along(cells)) {
+ append(combined, fvm_build_mechanism_data(gprop, cells[cell_idx], D, cell_idx));
+ }
+ return combined;
+}
- struct ion_segment_data {
- cable_cell_ion_data ion_data;
- bool mech_requires = false;
- bool mech_writes_iconc = false;
- bool mech_writes_econc = false;
- };
+fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop,
+ const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx)
+{
+ using size_type = fvm_size_type;
+ using index_type = fvm_index_type;
+ using value_type = fvm_value_type;
- std::unordered_map<std::string, std::map<index_type, ion_segment_data>> ion_segments;
+ const mechanism_catalogue& catalogue = *gprop.catalogue;
+ const auto& embedding = cell.embedding();
- // Temporary table for presence of mechanisms that read the reversal potential
- // of an ion channel, mapping ion name and cell index to a _sorted_
- // collection of segment indices.
+ const auto& global_dflt = gprop.default_parameters;
+ const auto& dflt = cell.default_parameters;
- std::unordered_map<std::string, std::unordered_map<size_type, std::set<size_type>>> ion_revpot_segments;
+ fvm_mechanism_data M;
- auto verify_ion_usage =
- [&gprop](const std::string& mech_name, const mechanism_info* info)
- {
+ // Verify mechanism ion usage, parameter values.
+ auto verify_mechanism = [&gprop](const mechanism_info& info, const mechanism_desc& desc) {
const auto& global_ions = gprop.ion_species;
- for (const auto& ion: info->ions) {
+ for (const auto& pv: desc.values()) {
+ if (!info.parameters.count(pv.first)) {
+ throw no_such_parameter(desc.name(), pv.first);
+ }
+ if (!info.parameters.at(pv.first).valid(pv.second)) {
+ throw invalid_parameter_value(desc.name(), pv.first, pv.second);
+ }
+ }
+
+ for (const auto& ion: info.ions) {
const auto& ion_name = ion.first;
const auto& ion_dep = ion.second;
if (!global_ions.count(ion_name)) {
throw cable_cell_error(
- "mechanism "+mech_name+" uses ion "+ion_name+ " which is missing in global properties");
+ "mechanism "+desc.name()+" uses ion "+ion_name+ " which is missing in global properties");
}
if (ion_dep.verify_ion_charge) {
if (ion_dep.expected_ion_charge!=global_ions.at(ion_name)) {
throw cable_cell_error(
- "mechanism "+mech_name+" uses ion "+ion_name+ " expecting a different valence");
+ "mechanism "+desc.name()+" uses ion "+ion_name+ " expecting a different valence");
}
}
- }
- };
- auto update_paramset_and_validate =
- [&catalogue,&verify_ion_usage]
- (const mechanism_desc& desc, std::shared_ptr<mechanism_info>& info, string_set& paramset)
- {
- auto& name = desc.name();
- if (!info) {
- info.reset(new mechanism_info(catalogue[name]));
- verify_ion_usage(name, info.get());
- }
- for (const auto& pv: desc.values()) {
- if (!paramset.count(pv.first)) {
- if (!info->parameters.count(pv.first)) {
- throw no_such_parameter(name, pv.first);
- }
- if (!info->parameters.at(pv.first).valid(pv.second)) {
- throw invalid_parameter_value(name, pv.first, pv.second);
- }
- paramset.insert(pv.first);
+ if (ion_dep.write_reversal_potential && (ion_dep.write_concentration_int || ion_dep.write_concentration_ext)) {
+ throw cable_cell_error("mechanism "+desc.name()+" writes both reversal potential and concentration");
}
}
};
- auto cell_segment_part = D.cell_segment_part();
- size_type target_count = 0;
-
- for (auto cell_idx: make_span(0, D.ncell)) {
- auto& cell = cells[cell_idx];
- auto seg_range = cell_segment_part[cell_idx];
- size_type target_offset = target_count;
+ // Track ion usage of mechanisms so that ions are only instantiated where required.
+ std::unordered_map<std::string, std::vector<index_type>> ion_support;
+ auto update_ion_support = [&ion_support](const mechanism_info& info, const std::vector<index_type>& cvs) {
+ arb_assert(util::is_sorted(cvs));
- auto add_ion_segment =
- [&gparam, &cell, &ion_segments]
- (const std::string& ion_name, size_type segment_idx, const cable_cell_parameter_set& seg_param, const ion_dependency* iondep = nullptr)
- {
- const auto& global_ion_data = gparam.ion_data;
- const auto& cell_ion_data = cell.default_parameters.ion_data;
- const auto& seg_ion_data = seg_param.ion_data;
-
- auto& ion_entry = ion_segments[ion_name];
-
- // New entry?
- util::optional<cable_cell_ion_data> opt_ion_data;
- if (!ion_entry.count(segment_idx)) {
- opt_ion_data = value_by_key(seg_ion_data, ion_name);
- if (!opt_ion_data) {
- opt_ion_data = value_by_key(cell_ion_data, ion_name);
- }
- if (!opt_ion_data) {
- opt_ion_data = value_by_key(global_ion_data, ion_name);
- }
- if (!opt_ion_data) {
- throw arbor_internal_error("missing entry for ion "+ion_name+" in cable_cell global defaults");
- }
- }
-
- auto& ion_entry_seg = ion_entry[segment_idx];
- if (opt_ion_data) {
- ion_entry_seg.ion_data = *opt_ion_data;
- }
- if (iondep) {
- ion_entry_seg.mech_requires = true;
- ion_entry_seg.mech_writes_iconc |= iondep->write_concentration_int;
- ion_entry_seg.mech_writes_econc |= iondep->write_concentration_ext;
- }
- };
-
- for (auto segment_idx: make_span(seg_range)) {
- const segment_ptr& seg = cell.segments()[segment_idx-seg_range.first];
-
- for (const mechanism_desc& desc: seg->mechanisms()) {
- const auto& name = desc.name();
-
- density_mech_data& entry = density_mech_table[name];
- update_paramset_and_validate(desc, entry.info, entry.paramset);
- entry.segments.emplace_back(segment_idx, &desc);
-
- for (const auto& ion_entry: entry.info->ions) {
- const std::string& ion_name = ion_entry.first;
- const ion_dependency& iondep = ion_entry.second;
-
- add_ion_segment(ion_name, segment_idx, seg->parameters, &iondep);
-
- if (ion_entry.second.read_reversal_potential) {
- ion_revpot_segments[ion_name][cell_idx].insert(segment_idx);
- }
- }
- }
+ for (const auto& ion: util::keys(info.ions)) {
+ auto& support = ion_support[ion];
+ support = unique_union(support, cvs);
}
+ };
- for (const auto& cellsyn_by_name: cell.synapses()) {
- const auto& name = cellsyn_by_name.first;
- point_mech_data& entry = point_mech_table[name];
+ std::unordered_map<std::string, mcable_map<double>> init_iconc_mask;
+ std::unordered_map<std::string, mcable_map<double>> init_econc_mask;
- for (const auto& placed_synapse: cellsyn_by_name.second) {
- mlocation loc = placed_synapse.loc;
- cell_lid_type target = target_offset + placed_synapse.lid;
- ++target_count;
- const mechanism_desc& desc = placed_synapse.item;
+ // Density mechanisms:
- update_paramset_and_validate(desc, entry.info, entry.paramset);
+ for (const auto& entry: cell.region_assignments().get<mechanism_desc>()) {
+ const std::string& name = entry.first;
+ mechanism_info info = catalogue[name];
- size_type cv = D.branch_location_cv(cell_idx, loc);
- entry.points.push_back({cv, target, &desc});
+ std::vector<double> param_dflt;
+ fvm_mechanism_config config;
+ config.kind = mechanismKind::density;
- const segment_ptr& seg = cell.segments()[loc.branch];
- size_type segment_idx = D.cell_segment_bounds[cell_idx]+loc.branch;
+ std::vector<std::string> param_names;
+ assign(param_names, util::keys(info.parameters));
+ sort(param_names);
- for (const auto& ion_entry: entry.info->ions) {
- const std::string& ion_name = ion_entry.first;
- const ion_dependency& iondep = ion_entry.second;
+ for (auto& p: param_names) {
+ config.param_values.emplace_back(p, std::vector<value_type>{});
+ param_dflt.push_back(info.parameters.at(p).default_value);
+ }
- add_ion_segment(ion_name, segment_idx, seg->parameters, &iondep);
+ mcable_map<double> support;
+ std::vector<mcable_map<double>> param_maps;
- if (ion_entry.second.read_reversal_potential) {
- ion_revpot_segments[ion_name][cell_idx].insert(segment_idx);
- }
+ {
+ std::unordered_map<std::string, mcable_map<double>> keyed_param_maps;
+ for (auto& on_cable: entry.second) {
+ verify_mechanism(info, on_cable.second);
+ mcable cable = on_cable.first;
+
+ support.insert(cable, 1.);
+ for (auto param_assign: on_cable.second.values()) {
+ keyed_param_maps[param_assign.first].insert(cable, param_assign.second);
}
}
+ for (auto& p: param_names) {
+ param_maps.push_back(std::move(keyed_param_maps[p]));
+ }
}
- for (const auto& loc_clamp: cell.stimuli()) {
- size_type cv = D.branch_location_cv(cell_idx, loc_clamp.loc);
- stimuli.push_back({cv, loc_clamp.item});
- }
-
- // Add segments to ion_segments map which intersect with existing segments, so
- // that each CV with an ion value is 100% covered.
+ std::vector<double> param_on_cv(config.param_values.size());
- for (auto& e: ion_segments) {
- const auto& name = e.first;
- const auto& ion_entry = e.second;
+ for (auto cv: D.geometry.cell_cvs(cell_idx)) {
+ double area = 0;
+ util::fill(param_on_cv, 0.);
- for (auto segment_idx: make_span(seg_range)) {
- index_type parent_segment_idx = D.parent_segment[segment_idx];
- const segment_ptr& parent_seg = cell.segments()[parent_segment_idx-seg_range.first];
+ for (mcable c: D.geometry.cables(cv)) {
+ double area_on_cable = embedding.integrate_area(c.branch, pw_over_cable(support, c, 0.));
+ if (!area_on_cable) continue;
- if (ion_entry.count(segment_idx)) {
- add_ion_segment(name, parent_segment_idx, parent_seg->parameters);
+ area += area_on_cable;
+ for (auto i: count_along(param_on_cv)) {
+ param_on_cv[i] += embedding.integrate_area(c.branch, pw_over_cable(param_maps[i], c, param_dflt[i]));
}
}
- for (auto segment_idx: make_span(seg_range)) {
- const segment_ptr& seg = cell.segments()[segment_idx-seg_range.first];
- index_type parent_segment_idx = D.parent_segment[segment_idx];
-
- if (ion_entry.count(parent_segment_idx)) {
- add_ion_segment(name, segment_idx, seg->parameters);
+ if (area>0) {
+ double oo_cv_area = 1./D.cv_area[cv];
+ config.cv.push_back(cv);
+ config.norm_area.push_back(area*oo_cv_area);
+ for (auto i: count_along(param_on_cv)) {
+ config.param_values[i].second.push_back(param_on_cv[i]*oo_cv_area);
}
}
}
-
- // Maintain a map of reversal potential mechanism to written ions for this
- // cell, to ensure consistency with assignments from the cell and global
- // parameters.
-
- std::unordered_multimap<std::string, std::string> revpot_to_ion;
- std::unordered_map<std::string, std::string> ion_to_revpot;
-
- auto add_revpot = [&](const std::string& ion, const mechanism_desc& desc) {
- const auto& name = desc.name();
-
- revpot_mech_data& entry = revpot_mech_table[name];
- update_paramset_and_validate(desc, entry.info, entry.paramset);
- ion_to_revpot[ion] = desc.name();
- for (auto dep: entry.info->ions) {
- if (dep.second.write_reversal_potential) {
- revpot_to_ion.insert({desc.name(), dep.first});
+ for (const auto& iondep: info.ions) {
+ if (iondep.second.write_concentration_int) {
+ for (auto c: support) {
+ bool ok = init_iconc_mask[iondep.first].insert(c.first, 0.);
+ if (!ok) {
+ throw cable_cell_error("overlapping ion concentration writing mechanism "+name);
+ }
}
}
- entry.cells.insert({cell_idx, &desc});
- };
-
- const auto& cellrevpot = cell.default_parameters.reversal_potential_method;
- for (const auto& revpot: cellrevpot) {
- add_revpot(revpot.first, revpot.second);
- }
-
- const auto& globalrevpot = gparam.reversal_potential_method;
- for (const auto& revpot: globalrevpot) {
- if (!cellrevpot.count(revpot.first)) {
- add_revpot(revpot.first, revpot.second);
+ if (iondep.second.write_concentration_ext) {
+ for (auto c: support) {
+ bool ok = init_econc_mask[iondep.first].insert(c.first, 0.);
+ if (!ok) {
+ throw cable_cell_error("overlapping ion concentration writing mechanism "+name);
+ }
+ }
}
}
- // Ensure that if a revpot mechanism writes to multiple ions, that
- // that mechanism is associated with all of them.
- for (auto& entry: revpot_to_ion) {
- auto declared_revpot = value_by_key(ion_to_revpot, entry.second);
- if (!declared_revpot || declared_revpot.value()!=entry.first) {
- throw cable_cell_error("mechanism "+entry.first+" writes reversal potential for "
- +entry.second+", which is not configured to use "+entry.first);
- }
- }
+ update_ion_support(info, config.cv);
+ M.mechanisms[name] = std::move(config);
}
+ // Synapses:
- // II. Build ion and mechanism configs.
+ struct synapse_instance {
+ size_type cv;
+ std::map<std::string, double> param_value; // uses ordering of std::map
+ size_type target_index;
+ };
- // Shared temporary lookup info across mechanism instances, set by build_param_data.
- string_index_map param_index; // maps parameter name to parameter index
- std::vector<std::string> param_name; // indexed by parameter index
- std::vector<value_type> param_default; // indexed by parameter index
+ for (const auto& entry: cell.synapses()) {
+ const std::string& name = entry.first;
+ mechanism_info info = catalogue[name];
+ std::map<std::string, double> default_param_value;
+ std::vector<synapse_instance> sl;
- auto build_param_data =
- [¶m_name, ¶m_index, ¶m_default](const string_set& paramset, const mechanism_info* info)
- {
- assign(param_name, paramset);
- auto nparam = paramset.size();
+ for (const auto& kv: info.parameters) {
+ default_param_value[kv.first] = kv.second.default_value;
+ }
- assign(param_default, transform_view(param_name,
- [info](const std::string& p) { return info->parameters.at(p).default_value; }));
+ for (const placed<mechanism_desc>& pm: entry.second) {
+ verify_mechanism(info, pm.item);
- param_index.clear();
- for (auto i: make_span(0, nparam)) {
- param_index[param_name[i]] = i;
- }
- return nparam;
- };
+ synapse_instance in;
+ in.param_value = default_param_value;
- // IIa. Ion channel CVs.
+ for (const auto& kv: pm.item.values()) {
+ in.param_value.at(kv.first) = kv.second;
+ }
- for (const auto& ionseg: ion_segments) {
- auto& seg_ion_map = ion_segments[ionseg.first];
- fvm_ion_config& ion_config = mechdata.ions[ionseg.first];
- auto& seg_ion_data = ionseg.second;
+ in.target_index = pm.lid;
+ in.cv = D.geometry.location_cv(cell_idx, pm.loc);
+ sl.push_back(std::move(in));
+ }
- for_each_cv_in_segments(D, keys(seg_ion_data),
- [&ion_config, &seg_ion_map](size_type cv_index, index_type cv, value_type area, index_type seg_index, index_type seg) {
- if (seg_ion_map[seg].mech_requires) {
- if (!util::binary_search_index(ion_config.cv, cv)) {
- ion_config.cv.push_back(cv);
- }
- }
- });
+ // Permute synapse instances so that they are in increasing order
+ // (lexicographically) by CV, param_value set, and target, so that
+ // instances in the same CV with the same parameter values are adjacent.
+ // cv_order[i] is the index of the ith instance by this ordering.
- ion_config.init_iconc.resize(ion_config.cv.size());
- ion_config.init_econc.resize(ion_config.cv.size());
- ion_config.reset_iconc.resize(ion_config.cv.size());
- ion_config.reset_econc.resize(ion_config.cv.size());
- ion_config.init_revpot.resize(ion_config.cv.size());
+ std::vector<size_type> cv_order;
+ assign(cv_order, count_along(sl));
+ sort_by(cv_order, [&](size_type i) {
+ return std::tie(sl[i].cv, sl[i].param_value, sl[i].target_index);
+ });
- for_each_cv_in_segments(D, keys(seg_ion_data),
- [&ion_config, &seg_ion_map, &D](size_type cv_index, index_type cv, value_type area, index_type seg_index, index_type seg) {
- auto opt_i = util::binary_search_index(ion_config.cv, cv);
- if (!opt_i) return;
+ bool coalesce = catalogue[name].linear && gprop.coalesce_synapses;
- std::size_t i = *opt_i;
- auto& seg_ion_entry = seg_ion_map[seg];
+ fvm_mechanism_config config;
+ config.kind = mechanismKind::point;
+ for (auto& pentry: default_param_value) {
+ config.param_values.emplace_back(pentry.first, std::vector<value_type>{});
+ }
- value_type weight = area/D.cv_area[cv];
- ion_config.reset_iconc[i] += weight*seg_ion_entry.ion_data.init_int_concentration;
- ion_config.reset_econc[i] += weight*seg_ion_entry.ion_data.init_ext_concentration;
+ const synapse_instance* prev = nullptr;
+ for (auto i: cv_order) {
+ const auto& in = sl[i];
- if (!seg_ion_entry.mech_writes_iconc) {
- ion_config.init_iconc[i] += weight*seg_ion_entry.ion_data.init_int_concentration;
+ if (coalesce && prev && prev->cv==in.cv && prev->param_value==in.param_value) {
+ ++config.multiplicity.back();
+ }
+ else {
+ config.cv.push_back(in.cv);
+ if (coalesce) {
+ config.multiplicity.push_back(1);
}
- if (!seg_ion_entry.mech_writes_econc) {
- ion_config.init_econc[i] += weight*seg_ion_entry.ion_data.init_ext_concentration;
+ unsigned j = 0;
+ for (auto& pentry: in.param_value) {
+ arb_assert(config.param_values[j].first==pentry.first);
+ config.param_values[j++].second.push_back(pentry.second);
}
+ }
+ config.target.push_back(in.target_index);
- // Reversal potentials are not area weighted, and are overridden at the
- // per-cell level by any supplied revpot mechanisms.
- ion_config.init_revpot[i] = seg_ion_entry.ion_data.init_reversal_potential;
- });
- }
-
- // IIb. Reversal potential mechanism CVs and parameters.
-
- for (const auto& entry: revpot_mech_table) {
- const std::string& name = entry.first;
- const mechanism_info& info = *entry.second.info;
-
- fvm_mechanism_config& config = mechdata.mechanisms[name];
- config.kind = mechanismKind::revpot;
-
- auto nparam = build_param_data(entry.second.paramset, &info);
- config.param_values.resize(nparam);
- for (auto pidx: make_span(nparam)) {
- config.param_values[pidx].first = param_name[pidx];
+ prev = ∈
}
- for (auto& cell_entry: entry.second.cells) {
- auto cell_idx = cell_entry.first;
- std::vector<index_type> segment_indices;
-
- for (auto& mech_ion_dep_entry: info.ions) {
- auto& ion_name = mech_ion_dep_entry.first;
- auto& ion_dep = mech_ion_dep_entry.second;
+ // If synapse uses an ion, add to ion support.
+ update_ion_support(info, config.cv);
- if (!ion_dep.write_reversal_potential) continue;
+ M.n_target += config.target.size();
+ M.mechanisms[name] = std::move(config);
+ }
- const auto& segments = value_by_key(ion_revpot_segments[ion_name], cell_idx);
- if (!segments) continue;
+ // Stimuli:
- for (auto seg_idx: segments.value()) {
- segment_indices.push_back(seg_idx);
- }
- }
+ if (!cell.stimuli().empty()) {
+ const auto& stimuli = cell.stimuli();
- const mechanism_desc& desc = *cell_entry.second;
- std::vector<value_type> pval = param_default;
+ std::vector<size_type> stimuli_cv;
+ assign_by(stimuli_cv, stimuli, [&D, cell_idx](auto& p) { return D.geometry.location_cv(cell_idx, p.loc); });
- for (auto pidx: make_span(nparam)) {
- if (auto opt_v = value_by_key(desc.values(), param_name[pidx])) {
- pval[pidx] = opt_v.value();
- }
- }
-
- size_type config_offset = config.cv.size();
- for_each_cv_in_segments(D, segment_indices,
- [&](size_type cv_index, index_type cv, auto, auto, auto) {
- extend_at(config.cv, config_offset+cv_index) = cv;
+ std::vector<size_type> cv_order;
+ assign(cv_order, count_along(stimuli));
+ sort_by(cv_order, [&](size_type i) { return stimuli_cv[i]; });
- for (auto pidx: make_span(nparam)) {
- extend_at(config.param_values[pidx].second, config_offset+cv_index) = pval[pidx];
- }
- });
+ fvm_mechanism_config config;
+ config.kind = mechanismKind::point;
+ // (param_values entries must be ordered by parameter name)
+ config.param_values = {{"amplitude", {}}, {"delay", {}}, {"duration", {}}};
+
+ for (auto i: cv_order) {
+ config.cv.push_back(stimuli_cv[i]);
+ config.param_values[0].second.push_back(stimuli[i].item.amplitude);
+ config.param_values[1].second.push_back(stimuli[i].item.delay);
+ config.param_values[2].second.push_back(stimuli[i].item.duration);
}
- }
- // Remove any reversal potential mechanisms that ultimately have no extent.
- for (auto i = mechdata.mechanisms.begin(); i!=mechdata.mechanisms.end(); ) {
- i = i->second.cv.empty()? mechdata.mechanisms.erase(i): std::next(i);
+ M.mechanisms["_builtin_stimulus"] = std::move(config);
}
- // IIc. Density mechanism CVs, parameters and ionic default concentration contributions.
+ // Ions:
- // Ameliorate area sum rounding errors by clamping normalized area contributions to [0, 1]
- // and rounding values within an epsilon of 0 or 1 to that value.
- auto trim = [](value_type& v) {
- constexpr value_type eps = std::numeric_limits<value_type>::epsilon()*4;
- v = v<eps? 0: v+eps>1? 1: v;
- };
+ auto initial_ion_data_map = cell.region_assignments().get<initial_ion_data>();
- for (const auto& entry: density_mech_table) {
- const std::string& name = entry.first;
- fvm_mechanism_config& config = mechdata.mechanisms[name];
- config.kind = mechanismKind::density;
-
- auto nparam = build_param_data(entry.second.paramset, entry.second.info.get());
-
- // In order to properly account for partially overriden parameters in CVs
- // that are shared between segments, we need to track not only the area-weighted
- // sum of parameter values, but also the total area for each CV for each parameter
- // that has been overriden — the remaining area demands a contribution from the
- // parameter default value.
+ for (const auto& ion_cvs: ion_support) {
+ const std::string& ion = ion_cvs.first;
- std::vector<std::vector<value_type>> param_value(nparam);
- std::vector<std::vector<value_type>> param_area_contrib(nparam);
+ fvm_ion_config config;
+ config.cv = ion_cvs.second;
- for_each_cv_in_segments(D, keys(entry.second.segments),
- [&](size_type cv_index, index_type cv, value_type area, index_type seg_index, index_type seg)
- {
- const mechanism_desc& mech_desc = *entry.second.segments[seg_index].second;
+ auto n_cv = config.cv.size();
+ config.init_iconc.resize(n_cv);
+ config.init_econc.resize(n_cv);
+ config.reset_iconc.resize(n_cv);
+ config.reset_econc.resize(n_cv);
+ config.init_revpot.resize(n_cv);
- extend_at(config.cv, cv_index) = cv;
+ cable_cell_ion_data ion_data = *(value_by_key(dflt.ion_data, ion) | value_by_key(global_dflt.ion_data, ion));
+ double dflt_iconc = ion_data.init_int_concentration;
+ double dflt_econc = ion_data.init_ext_concentration;
+ double dflt_rvpot = ion_data.init_reversal_potential;
- for (auto& kv: mech_desc.values()) {
- int pidx = param_index.at(kv.first);
- value_type v = kv.second;
+ const mcable_map<initial_ion_data>& ion_on_cable = initial_ion_data_map[ion];
- extend_at(param_area_contrib[pidx], cv_index) += area;
- extend_at(param_value[pidx], cv_index) += area*v;
- }
-
- extend_at(config.norm_area, cv_index) += area;
- });
+ auto pw_times = [](const pw_elements<double>& a, const pw_elements<double>& b) {
+ return zip(a, b, [](double left, double right, pw_element<double> a, pw_element<double> b) { return a.second*b.second; });
+ };
- // Complete parameter values with default values.
+ for (auto i: count_along(config.cv)) {
+ auto cv = config.cv[i];
+ if (D.cv_area[cv]==0) continue;
- config.param_values.resize(nparam);
- for (auto pidx: make_span(0, nparam)) {
- value_type default_value = param_default[pidx];
- config.param_values[pidx].first = param_name[pidx];
+ for (mcable c: D.geometry.cables(cv)) {
+ auto iconc = pw_over_cable(ion_on_cable, c, dflt_iconc, [](auto x) { return x.initial.init_int_concentration; });
+ auto econc = pw_over_cable(ion_on_cable, c, dflt_econc, [](auto x) { return x.initial.init_ext_concentration; });
+ auto rvpot = pw_over_cable(ion_on_cable, c, dflt_rvpot, [](auto x) { return x.initial.init_reversal_potential; });
- auto& values = config.param_values[pidx].second;
- values.resize(config.cv.size());
+ config.reset_iconc[i] += embedding.integrate_area(c.branch, iconc);
+ config.reset_econc[i] += embedding.integrate_area(c.branch, econc);
+ config.init_revpot[i] += embedding.integrate_area(c.branch, rvpot);
- for (auto i: count_along(config.cv)) {
- value_type v = param_value[pidx][i];
- value_type cv_area = D.cv_area[config.cv[i]];
- value_type remaining_area = cv_area-param_area_contrib[pidx][i];
+ auto iconc_masked = pw_times(pw_over_cable(init_iconc_mask[ion], c, 1.), iconc);
+ auto econc_masked = pw_times(pw_over_cable(init_econc_mask[ion], c, 1.), econc);
- values[i] = (v+remaining_area*default_value)/cv_area;
+ config.init_iconc[i] += embedding.integrate_area(c.branch, iconc_masked);
+ config.init_econc[i] += embedding.integrate_area(c.branch, econc_masked);
}
- }
- // Normalize norm_area entries.
-
- for (auto i: count_along(config.cv)) {
- config.norm_area[i] /= D.cv_area[config.cv[i]];
- trim(config.norm_area[i]);
+ double oo_cv_area = 1./D.cv_area[cv];
+ config.reset_iconc[i] *= oo_cv_area;
+ config.reset_econc[i] *= oo_cv_area;
+ config.init_revpot[i] *= oo_cv_area;
+ config.init_iconc[i] *= oo_cv_area;
+ config.init_econc[i] *= oo_cv_area;
}
- }
-
- // II.3 Point mechanism CVs, targets, parameters and stimuli.
- for (const auto& entry: point_mech_table) {
- const std::string& name = entry.first;
- const auto& points = entry.second.points;
-
- auto nparam = build_param_data(entry.second.paramset, entry.second.info.get());
- std::vector<std::vector<value_type>> param_value(nparam);
-
- // Permute points in this mechanism so that they are in increasing CV order;
- // cv_order[i] is the index of the ith point by increasing CV.
+ M.ions[ion] = std::move(config);
+ }
- mechdata.ntarget += points.size();
+ std::unordered_map<std::string, mechanism_desc> revpot_tbl;
+ std::unordered_set<std::string> revpot_specified;
- std::vector<size_type> cv_order;
- assign(cv_order, count_along(points));
- sort_by(cv_order, [&](size_type i) { return points[i].cv; });
-
- fvm_mechanism_config& config = mechdata.mechanisms[name];
- config.kind = mechanismKind::point;
+ for (const auto& ion: util::keys(gprop.ion_species)) {
+ if (auto maybe_revpot = value_by_key(dflt.reversal_potential_method, ion)
+ | value_by_key(global_dflt.reversal_potential_method, ion))
+ {
+ const mechanism_desc& revpot = *maybe_revpot;
+ mechanism_info info = catalogue[revpot.name()];
+ verify_mechanism(info, revpot);
+ revpot_specified.insert(ion);
+
+ bool writes_this_revpot = false;
+ for (auto& iondep: info.ions) {
+ if (iondep.second.write_reversal_potential) {
+ if (revpot_tbl.count(iondep.first)) {
+ auto& existing_revpot_desc = revpot_tbl.at(iondep.first);
+ if (existing_revpot_desc.name() != revpot.name() || existing_revpot_desc.values() != revpot.values()) {
+ throw cable_cell_error("inconsistent revpot ion assignment for mechanism "+revpot.name());
+ }
+ }
+ else {
+ revpot_tbl[iondep.first] = revpot;
+ }
- // Generate config.cv: contains cv of group of synapses that can be coalesced into one instance
- // Generate config.param_values: contains parameters of group of synapses that can be coalesced into one instance
- // Generate multiplicity: contains number of synapses in each coalesced group of synapses
- // Generate target: contains the synapse target number
- if (catalogue[name].linear && gprop.coalesce_synapses) {
- // cv_param_vec used to lexicographically sort the cv, parameters and target, which are stored in that order
- std::vector<cv_param> cv_param_vec(cv_order.size());
-
- for (unsigned i = 0; i < cv_order.size(); ++i) {
- auto loc = cv_order[i];
- std::vector<value_type> p(nparam);
- for (auto pidx: make_span(0, nparam)) {
- value_type pdefault = param_default[pidx];
- const std::string& pname = param_name[pidx];
- p[pidx] = value_by_key(points[loc].desc->values(), pname).value_or(pdefault);
+ writes_this_revpot |= iondep.first==ion;
}
- cv_param_vec[i] = cv_param{points[loc].cv, p, points[loc].target_index};
}
- std::sort(cv_param_vec.begin(), cv_param_vec.end());
-
- auto identical_synapse = [](const cv_param& i, const cv_param& j) {
- return i.cv==j.cv && i.params==j.params;
- };
-
- config.param_values.resize(nparam);
- for (auto pidx: make_span(0, nparam)) {
- config.param_values[pidx].first = param_name[pidx];
+ if (!writes_this_revpot) {
+ throw cable_cell_error("revpot mechanism for ion "+ion+" does not write this reversal potential");
}
- config.target.reserve(cv_param_vec.size());
- for (auto i = cv_param_vec.begin(), j = i; i!=cv_param_vec.end(); i = j) {
- ++j;
- while (j!=cv_param_vec.end() && identical_synapse(*i, *j)) ++j;
-
- auto mergeable = util::make_range(i, j);
-
- config.cv.push_back((*i).cv);
- for (auto pidx: make_span(0, nparam)) {
- config.param_values[pidx].second.push_back((*i).params[pidx]);
- }
- config.multiplicity.push_back(mergeable.size());
+ // Only instantiate if the ion is used.
+ if (M.ions.count(ion)) {
+ // Revpot mechanism already configured? Add cvs for this ion too.
+ if (M.mechanisms.count(revpot.name())) {
+ fvm_mechanism_config& config = M.mechanisms[revpot.name()];
+ config.cv = unique_union(config.cv, M.ions[ion].cv);
+ config.norm_area.assign(config.cv.size(), 1.);
- for (auto e: mergeable) {
- config.target.push_back(e.target);
+ for (auto& pv: config.param_values) {
+ pv.second.assign(config.cv.size(), pv.second.front());
+ }
}
- }
- }
- else {
- assign(config.cv, transform_view(cv_order, [&](size_type j) { return points[j].cv; }));
- assign(config.target, transform_view(cv_order, [&](size_type j) { return points[j].target_index; }));
+ else {
+ fvm_mechanism_config config;
+ config.kind = mechanismKind::revpot;
+ config.cv = M.ions[ion].cv;
+ config.norm_area.assign(config.cv.size(), 1.);
+
+ std::map<std::string, double> param_value; // uses ordering of std::map
+ for (const auto& kv: info.parameters) {
+ param_value[kv.first] = kv.second.default_value;
+ }
- config.param_values.resize(nparam);
- for (auto pidx: make_span(0, nparam)) {
- value_type pdefault = param_default[pidx];
- const std::string& pname = param_name[pidx];
+ for (auto& kv: revpot.values()) {
+ param_value[kv.first] = kv.second;
+ }
- config.param_values[pidx].first = pname;
+ for (auto& kv: param_value) {
+ config.param_values.emplace_back(kv.first, std::vector<value_type>(config.cv.size(), kv.second));
+ }
- auto& values = config.param_values[pidx].second;
- assign(values, transform_view(cv_order,
- [&](size_type j) { return value_by_key(points[j].desc->values(), pname).value_or(pdefault); }));
+ M.mechanisms[revpot.name()] = std::move(config);
+ }
}
}
}
- // Sort stimuli by ascending CV and construct parameter vectors.
- if (!stimuli.empty()) {
- fvm_mechanism_config& stim_config = mechdata.mechanisms["_builtin_stimulus"];
- using cv_clamp = const std::pair<size_type, i_clamp>&;
-
- auto stim_cv_field = [](cv_clamp p) { return p.first; };
- sort_by(stimuli, stim_cv_field);
- assign(stim_config.cv, transform_view(stimuli, stim_cv_field));
-
- stim_config.param_values.resize(3);
-
- stim_config.param_values[0].first = "delay";
- assign(stim_config.param_values[0].second,
- transform_view(stimuli, [](cv_clamp p) { return p.second.delay; }));
-
- stim_config.param_values[1].first = "duration";
- assign(stim_config.param_values[1].second,
- transform_view(stimuli, [](cv_clamp p) { return p.second.duration; }));
-
- stim_config.param_values[2].first = "amplitude";
- assign(stim_config.param_values[2].second,
- transform_view(stimuli, [](cv_clamp p) { return p.second.amplitude; }));
+ // Confirm that all ions written to by a revpot have a corresponding entry in a reversal_potential_method table.
+ for (auto& kv: revpot_tbl) {
+ if (!revpot_specified.count(kv.first)) {
+ throw cable_cell_error("revpot mechanism "+kv.second.name()+" also writes to ion "+kv.first);
+ }
}
- return mechdata;
+ return M;
}
} // namespace arb
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index d02475dae0a0854584dbad5b7d3d0ea9352e678b..3e634cd372f107c76562f79be62ae3b96c3655d9 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -1,5 +1,9 @@
#pragma once
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
#include <arbor/fvm_types.hpp>
#include <arbor/cable_cell.hpp>
#include <arbor/mechanism.hpp>
@@ -7,61 +11,91 @@
#include <arbor/mechcat.hpp>
#include <arbor/recipe.hpp>
-#include "fvm_compartment.hpp"
+#include "util/piecewise.hpp"
+#include "util/rangeutil.hpp"
#include "util/span.hpp"
namespace arb {
-// Discretization data for an unbranched segment.
+// CV geometry as determined by per-cell CV boundary points.
-struct segment_info {
- using value_type = fvm_value_type;
+struct cv_geometry {
+ using size_type = fvm_size_type;
using index_type = fvm_index_type;
- bool soma_parent = false; // segments parent is soma
- value_type parent_cv_area = 0;
- value_type distal_cv_area = 0;
+ 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_to_cell; // Maps CV index to cell index.
+ std::vector<index_type> cell_cv_divs; // Partitions CV indices by cell.
- static constexpr index_type npos = -1;
+ // 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;
- index_type parent_cv = npos; // npos => no parent.
- index_type proximal_cv = 0; // First CV in segment, excluding parent.
- index_type distal_cv = 0; // Last CV in segment (may be shared with other segments).
+ auto cables(size_type cv_index) const {
+ auto partn = util::partition_view(cv_cables_divs);
+ return util::subrange_view(cv_cables, partn[cv_index]);
+ }
- bool has_parent() const { return parent_cv!=npos; }
+ std::pair<index_type, index_type> cell_cv_interval(size_type cell_idx) const {
+ auto partn = util::partition_view(cell_cv_divs);
+ return partn[cell_idx];
+ }
- // Range of CV-indices for segment, excluding parent.
- std::pair<index_type, index_type> cv_range() const {
- return {proximal_cv, 1+distal_cv};
+ auto cell_cvs(size_type cell_idx) const {
+ auto partn = util::partition_view(cell_cv_divs);
+ return util::make_span(partn[cell_idx]);
}
- // Position is proportional distal distance along segment, in [0, 1).
- index_type cv_by_position(double pos) const {
- index_type n = distal_cv+1-proximal_cv;
- index_type i = static_cast<index_type>(n*pos+0.5);
- if (i>0) {
- return proximal_cv+(i-1);
- }
- else {
- return parent_cv==npos? proximal_cv: parent_cv;
- }
+ size_type size() const {
+ arb_assert((cv_parent.empty() && cv_cables_divs.empty() &&
+ cv_cables.empty() && cv_to_cell.empty())
+ ||
+ (cv_parent.size()+1 == cv_cables_divs.size() &&
+ cv_parent.size() == cv_to_cell.size() &&
+ (unsigned)cv_to_cell.back()+1 == cell_cv_divs.size()-1));
+
+ return cv_parent.size();
+ }
+
+ bool empty() const {
+ return size()==0;
+ }
+
+ size_type n_cell() const {
+ return cell_cv_divs.empty()? 0: cell_cv_divs.size()-1;
+ }
+
+ size_type location_cv(size_type cell_idx, mlocation loc) const {
+ return cell_cv_divs.at(cell_idx)+branch_cv_map.at(cell_idx).at(loc.branch)(loc.pos).second;
}
};
-// Discretization of morphologies and electrical properties for
-// cells in a cell group.
+// Combine two cv_geometry groups in-place.
+// (Returns reference to first argument.)
+cv_geometry& append(cv_geometry&, const cv_geometry&);
-struct fvm_discretization {
- using value_type = fvm_value_type;
+// 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.
+//
+// diam_um is taken to be the diameter of a CV with constant diameter and same
+// extent which has the same surface area (i.e. cv_area/(Ï€L) where L is the
+// total length of the cables comprising the CV.)
+
+struct fvm_cv_discretization {
using size_type = fvm_size_type;
- using index_type = fvm_index_type; // In particular, used for CV indices.
+ using index_type = fvm_index_type;
+ using value_type = fvm_value_type;
- size_type ncell;
- size_type ncv;
+ cv_geometry geometry;
- // Note: if CV j has no parent, parent_cv[j] = j. TODO: confirm!
- std::vector<index_type> parent_cv;
- std::vector<index_type> cv_to_cell;
+ bool empty() const { return geometry.empty(); }
+ size_type size() const { return geometry.size(); }
+ size_type n_cell() const { return geometry.n_cell(); }
std::vector<value_type> face_conductance; // [µS]
std::vector<value_type> cv_area; // [µm²]
@@ -69,34 +103,16 @@ struct fvm_discretization {
std::vector<value_type> init_membrane_potential; // [mV]
std::vector<value_type> temperature_K; // [K]
std::vector<value_type> diam_um; // [µm]
-
- std::vector<segment_info> segments;
-
- // If segment has no parent segment, parent_segment[j] = j.
- std::vector<index_type> parent_segment;
-
- std::vector<size_type> cell_segment_bounds; // Partitions segment indices by cell.
- std::vector<index_type> cell_cv_bounds; // Partitions CV indices by cell.
-
- auto cell_segment_part() const {
- return util::partition_view(cell_segment_bounds);
- }
-
- auto cell_cv_part() const {
- return util::partition_view(cell_cv_bounds);
- }
-
- size_type branch_location_cv(size_type cell_index, mlocation loc) const {
- auto cell_segs = cell_segment_part()[cell_index];
-
- size_type seg = loc.branch+cell_segs.first;
- arb_assert(seg<cell_segs.second);
- return segments[seg].cv_by_position(loc.pos);
- }
};
-fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells, const cable_cell_parameter_set& params);
+// Combine two fvm_cv_geometry groups in-place.
+// (Returns reference to first argument.)
+fvm_cv_discretization& append(fvm_cv_discretization&, const fvm_cv_discretization&);
+
+// Construct fvm_cv_discretization from one or more cells.
+fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global_dflt);
+fvm_cv_discretization fvm_cv_discretize(const std::vector<cable_cell>& cells, const cable_cell_parameter_set& global_defaults);
// Post-discretization data for point and density mechanism instantiation.
@@ -153,9 +169,9 @@ struct fvm_mechanism_data {
std::unordered_map<std::string, fvm_ion_config> ions;
// Total number of targets (point-mechanism points)
- std::size_t ntarget = 0;
+ std::size_t n_target = 0;
};
-fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop, const std::vector<cable_cell>& cells, const fvm_discretization& D);
+fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop, const std::vector<cable_cell>& cells, const fvm_cv_discretization& D);
} // namespace arb
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 32642764db8f1696963ce028a0d89a6c6f136246..15580a16633350c052ebb972c411a53b5fb95d8c 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -66,7 +66,7 @@ public:
const std::vector<cable_cell>& cells,
const std::vector<cell_gid_type>& gids,
const recipe& rec,
- const fvm_discretization& D);
+ const fvm_cv_discretization& D);
// Generates indom index for every gid, guarantees that gids belonging to the same supercell are in the same intdom
// Fills cell_to_intdom map; returns number of intdoms
@@ -397,14 +397,15 @@ void fvm_lowered_cell_impl<B>::initialize(
// Discretize cells, build matrix.
- fvm_discretization D = fvm_discretize(cells, global_props.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, global_props.default_parameters);
- std::vector<index_type> cv_to_intdom(D.ncv);
- std::transform(D.cv_to_cell.begin(), D.cv_to_cell.end(), cv_to_intdom.begin(),
+ std::vector<index_type> cv_to_intdom(D.size());
+ std::transform(D.geometry.cv_to_cell.begin(), D.geometry.cv_to_cell.end(), cv_to_intdom.begin(),
[&cell_to_intdom](index_type i){ return cell_to_intdom[i]; });
- arb_assert(D.ncell == ncell);
- matrix_ = matrix<backend>(D.parent_cv, D.cell_cv_bounds, D.cv_capacitance, D.face_conductance, D.cv_area, cell_to_intdom);
+ arb_assert(D.n_cell() == ncell);
+ matrix_ = matrix<backend>(D.geometry.cv_parent, D.geometry.cell_cv_divs,
+ D.cv_capacitance, D.face_conductance, D.cv_area, cell_to_intdom);
sample_events_ = sample_event_stream(num_intdoms);
// Discretize mechanism data.
@@ -439,7 +440,7 @@ void fvm_lowered_cell_impl<B>::initialize(
}
}
- target_handles.resize(mech_data.ntarget);
+ target_handles.resize(mech_data.n_target);
unsigned mech_id = 0;
for (auto& m: mech_data.mechanisms) {
@@ -517,7 +518,7 @@ void fvm_lowered_cell_impl<B>::initialize(
cell_gid_type gid = gids[cell_idx];
for (auto entry: cells[cell_idx].detectors()) {
- detector_cv.push_back(D.branch_location_cv(cell_idx, entry.loc));
+ detector_cv.push_back(D.geometry.location_cv(cell_idx, entry.loc));
detector_threshold.push_back(entry.item.threshold);
}
@@ -525,7 +526,7 @@ void fvm_lowered_cell_impl<B>::initialize(
probe_info pi = rec.get_probe({gid, j});
auto where = any_cast<cell_probe_address>(pi.address);
- auto cv = D.branch_location_cv(cell_idx, where.location);
+ auto cv = D.geometry.location_cv(cell_idx, where.location);
probe_handle handle;
switch (where.kind) {
@@ -553,21 +554,20 @@ template <typename B>
std::vector<fvm_gap_junction> fvm_lowered_cell_impl<B>::fvm_gap_junctions(
const std::vector<cable_cell>& cells,
const std::vector<cell_gid_type>& gids,
- const recipe& rec, const fvm_discretization& D) {
+ const recipe& rec, const fvm_cv_discretization& D) {
std::vector<fvm_gap_junction> v;
std::unordered_map<cell_gid_type, std::vector<unsigned>> gid_to_cvs;
- for (auto cell_idx: util::make_span(0, D.ncell)) {
+ for (auto cell_idx: util::make_span(0, D.n_cell())) {
+ if (!rec.num_gap_junction_sites(gids[cell_idx])) continue;
- if (rec.num_gap_junction_sites(gids[cell_idx])) {
- gid_to_cvs[gids[cell_idx]].reserve(rec.num_gap_junction_sites(gids[cell_idx]));
+ gid_to_cvs[gids[cell_idx]].reserve(rec.num_gap_junction_sites(gids[cell_idx]));
+ const auto& cell_gj = cells[cell_idx].gap_junction_sites();
- const auto& cell_gj = cells[cell_idx].gap_junction_sites();
- for (auto gj : cell_gj) {
- auto cv = D.branch_location_cv(cell_idx, gj.loc);
- gid_to_cvs[gids[cell_idx]].push_back(cv);
- }
+ for (auto gj : cell_gj) {
+ auto cv = D.geometry.location_cv(cell_idx, gj.loc);
+ gid_to_cvs[gids[cell_idx]].push_back(cv);
}
}
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index b87d5e6efa9c2211b542a61f0dafd13a68bd456d..dc742b5a2ea7d36b7445afa36ad66208d8a96938 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -15,7 +15,6 @@
#include <arbor/morph/mprovider.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
-#include <arbor/segment.hpp>
#include <arbor/util/typed_map.hpp>
namespace arb {
@@ -79,7 +78,7 @@ using cable_cell_region_map = static_typed_map<region_assignment,
using cable_cell_location_map = static_typed_map<location_assignment,
mechanism_desc, i_clamp, gap_junction_site, threshold_detector>;
-// High-level abstract representation of a cell and its segments
+// High-level abstract representation of a cell.
class cable_cell {
public:
using index_type = cell_lid_type;
@@ -101,18 +100,13 @@ public:
cable_cell(cable_cell&& other) = default;
/// construct from morphology
- cable_cell(const class morphology& m,
- const label_dict& dictionary={},
- bool compartments_from_discretization=false);
+ cable_cell(const class morphology& m, const label_dict& dictionary={});
/// Access to morphology and embedding
const concrete_embedding& embedding() const;
const arb::morphology& morphology() const;
const mprovider& provider() const;
- // the number of branches in the cell
- size_type num_branches() const;
-
// Set cell-wide default physical and ion parameters.
void set_default(init_membrane_potential prop) {
@@ -139,39 +133,6 @@ public:
default_parameters.reversal_potential_method[prop.ion] = prop.method;
}
- // All of the members marked with LEGACY below will be removed once
- // the discretization code has moved from consuming segments to em_morphology.
-
- // LEGACY
- bool has_soma() const;
-
- // LEGACY
- const class segment* parent(index_type index) const;
- // LEGACY
- const class segment* segment(index_type index) const;
-
- // access pointer to the soma
- // returns nullptr if the cell has no soma
- // LEGACY
- const soma_segment* soma() const;
-
- // access pointer to a cable segment
- // will throw an cable_cell_error exception if
- // the cable index is not valid
- // LEGACY
- const cable_segment* cable(index_type index) const;
-
- // LEGACY
- const std::vector<segment_ptr>& segments() const;
-
- // return a vector with the compartment count for each segment in the cell
- // LEGACY
- std::vector<size_type> compartment_counts() const;
-
- // The total number of compartments in the discretised cell.
- // LEGACY
- size_type num_compartments() const;
-
// Painters and placers.
//
// Used to describe regions and locations where density channels, stimuli,
@@ -221,27 +182,6 @@ public:
const cable_cell_region_map& region_assignments() const;
const cable_cell_location_map& location_assignments() const;
- // Checks that two cells have the same
- // - number and type of segments
- // - volume and area properties of each segment
- // - number of compartments in each segment
- // (note: just used for testing: move to test code?)
- friend bool cell_basic_equality(const cable_cell&, const cable_cell&);
-
- // Public view of parent indices vector.
- const std::vector<index_type>& parents() const;
-
- // Approximate per-segment mean attenuation b(f) at given frequency f,
- // ignoring membrane resistance [1/µm].
- value_type segment_mean_attenuation(value_type frequency, index_type segidx,
- const cable_cell_parameter_set& global_defaults) const;
-
- // Estimate of length constant λ(f) = 1/2 · 1/b(f), following
- // Hines and Carnevale (2001), "NEURON: A Tool for Neuroscientists",
- // Neuroscientist 7, pp. 123-135.
- value_type segment_length_constant(value_type frequency, index_type segidx,
- const cable_cell_parameter_set& global_defaults) const;
-
private:
std::unique_ptr<cable_cell_impl, void (*)(cable_cell_impl*)> impl_;
};
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index 3cc08754792c3372cd70767500b1bc5d6bf789d4..2e041e7dfd4e823984fe297cb269d6352a06cd56 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -271,6 +271,20 @@ private:
cv_policy_flag::value flags_;
};
+struct cv_policy_every_sample: cv_policy_base {
+ explicit cv_policy_every_sample(cv_policy_flag::value flags = cv_policy_flag::none):
+ flags_(flags) {}
+
+ cv_policy_base_ptr clone() const override {
+ return cv_policy_base_ptr(new cv_policy_every_sample(*this));
+ }
+
+ locset cv_boundary_points(const cable_cell&) const override;
+
+private:
+ cv_policy_flag::value flags_;
+};
+
inline cv_policy default_cv_policy() {
return cv_policy_fixed_per_branch(1);
}
@@ -294,7 +308,7 @@ struct cable_cell_parameter_set {
std::unordered_map<std::string, cable_cell_ion_data> ion_data;
std::unordered_map<std::string, mechanism_desc> reversal_potential_method;
- cv_policy discretization = default_cv_policy();
+ util::optional<cv_policy> discretization;
};
extern cable_cell_parameter_set neuron_parameter_defaults;
diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp
index da2688ec57e41234d9d0ac08c6515842caa7fa6f..d4785d89a8e9bd53260ff89cadaed3f56083d505 100644
--- a/arbor/include/arbor/morph/locset.hpp
+++ b/arbor/include/arbor/morph/locset.hpp
@@ -138,6 +138,9 @@ locset named(std::string);
// The null (empty) set.
locset nil();
+// Proportional location on every branch.
+locset on_branches(double pos);
+
} // namespace ls
// union of two locsets
diff --git a/arbor/include/arbor/segment.hpp b/arbor/include/arbor/segment.hpp
deleted file mode 100644
index b87bcf36e8a54bb65d0c7fb41895e53c4f7f6031..0000000000000000000000000000000000000000
--- a/arbor/include/arbor/segment.hpp
+++ /dev/null
@@ -1,345 +0,0 @@
-#pragma once
-
-#include <algorithm>
-#include <cmath>
-#include <memory>
-#include <numeric>
-#include <stdexcept>
-#include <string>
-#include <unordered_map>
-#include <vector>
-
-#include <arbor/assert.hpp>
-#include <arbor/cable_cell_param.hpp>
-#include <arbor/common_types.hpp>
-#include <arbor/math.hpp>
-#include <arbor/mechinfo.hpp>
-#include <arbor/point.hpp>
-#include <arbor/util/optional.hpp>
-
-namespace arb {
-
-enum class section_kind {
- soma,
- dendrite,
- axon,
- none
-};
-
-// forward declarations of segment specializations
-class soma_segment;
-class cable_segment;
-
-// abstract base class for a cell segment
-class segment {
-public:
- using value_type = double;
- using size_type = cell_local_size_type;
- using point_type = point<value_type>;
-
- // (Yet more motivation for a separate morphology description class!)
- virtual std::unique_ptr<segment> clone() const = 0;
-
- section_kind kind() const {
- return kind_;
- }
-
- bool is_soma() const
- {
- return kind_==section_kind::soma;
- }
-
- bool is_dendrite() const
- {
- return kind_==section_kind::dendrite;
- }
-
- bool is_axon() const
- {
- return kind_==section_kind::axon;
- }
-
- virtual size_type num_compartments() const = 0;
- virtual void set_compartments(size_type) = 0;
-
- virtual ~segment() = default;
-
- virtual cable_segment* as_cable()
- {
- return nullptr;
- }
-
- virtual soma_segment* as_soma()
- {
- return nullptr;
- }
-
- virtual const cable_segment* as_cable() const
- {
- return nullptr;
- }
-
- virtual const soma_segment* as_soma() const
- {
- return nullptr;
- }
-
- virtual bool is_placeholder() const
- {
- return false;
- }
-
- util::optional<mechanism_desc&> mechanism(const std::string& name) {
- auto it = std::find_if(mechanisms_.begin(), mechanisms_.end(),
- [&](mechanism_desc& m) { return m.name()==name; });
- return it==mechanisms_.end()? util::nullopt: util::just(*it);
- }
-
- void add_mechanism(mechanism_desc mech) {
- auto m = mechanism(mech.name());
- if (m) {
- *m = std::move(mech);
- }
- else {
- mechanisms_.push_back(std::move(mech));
- }
- }
-
- const std::vector<mechanism_desc>& mechanisms() {
- return mechanisms_;
- }
-
- const std::vector<mechanism_desc>& mechanisms() const {
- return mechanisms_;
- }
-
- cable_cell_parameter_set parameters; // override cell and global defaults
-
-protected:
- segment(section_kind kind): kind_(kind) {}
-
- section_kind kind_;
- std::vector<mechanism_desc> mechanisms_;
-};
-
-class placeholder_segment : public segment {
-public:
- placeholder_segment(): segment(section_kind::none) {}
-
- std::unique_ptr<segment> clone() const override {
- // use default copy constructor
- return std::make_unique<placeholder_segment>(*this);
- }
-
- bool is_placeholder() const override
- {
- return true;
- }
-
- size_type num_compartments() const override
- {
- return 0;
- }
-
- virtual void set_compartments(size_type) override {}
-};
-
-class soma_segment : public segment {
-public:
- soma_segment() = delete;
-
- explicit soma_segment(value_type r, point_type c = point_type{}):
- segment(section_kind::soma), radius_{r}, center_(c) {}
-
- std::unique_ptr<segment> clone() const override {
- // use default copy constructor
- return std::make_unique<soma_segment>(*this);
- }
-
- value_type radius() const
- {
- return radius_;
- }
-
- point_type const& center() const
- {
- return center_;
- }
-
- soma_segment* as_soma() override
- {
- return this;
- }
-
- const soma_segment* as_soma() const override
- {
- return this;
- }
-
- /// soma has one and one only compartments
- size_type num_compartments() const override
- {
- return 1;
- }
-
- void set_compartments(size_type n) override {}
-
-private :
- // store the center and radius of the soma
- // note that the center may be undefined
- value_type radius_;
- point_type center_;
-};
-
-class cable_segment : public segment {
-public:
- using base = segment;
- using base::kind_;
- using base::value_type;
- using base::point_type;
-
- // delete the default constructor
- cable_segment() = delete;
-
- // constructors for a cable with no location information
- cable_segment(section_kind k, std::vector<value_type> r, std::vector<value_type> lens):
- segment(k), radii_(std::move(r)), lengths_(std::move(lens))
- {
- arb_assert(kind_==section_kind::dendrite || kind_==section_kind::axon);
- }
-
- cable_segment(section_kind k, value_type r1, value_type r2, value_type len):
- cable_segment{k, {r1, r2}, decltype(lengths_){len}}
- {}
-
- // constructor that lets the user describe the cable as a
- // seriew of radii and locations
- cable_segment(section_kind k, std::vector<value_type> r, std::vector<point_type> p):
- segment(k), radii_(std::move(r)), locations_(std::move(p))
- {
- arb_assert(kind_==section_kind::dendrite || kind_==section_kind::axon);
- update_lengths();
- }
-
- // helper that lets user specify with the radius and location
- // of just the end points of the cable
- // i.e. describing the cable as a single frustrum
- cable_segment(
- section_kind k,
- value_type r1,
- value_type r2,
- point_type const& p1,
- point_type const& p2
- ):
- cable_segment(k, {r1, r2}, {p1, p2})
- {}
-
- std::unique_ptr<segment> clone() const override {
- // use default copy constructor
- return std::make_unique<cable_segment>(*this);
- }
-
- value_type length() const
- {
- return std::accumulate(lengths_.begin(), lengths_.end(), value_type{});
- }
-
- bool has_locations() const
- {
- return locations_.size() > 0;
- }
-
- // the number sub-segments that define the cable segment
- size_type num_sub_segments() const
- {
- return radii_.size()-1;
- }
-
- std::vector<value_type> const& lengths() const
- {
- return lengths_;
- }
-
- std::vector<value_type> const& radii() const
- {
- return radii_;
- }
-
- cable_segment* as_cable() override
- {
- return this;
- }
-
- const cable_segment* as_cable() const override
- {
- return this;
- }
-
- size_type num_compartments() const override
- {
- return num_compartments_;
- }
-
- void set_compartments(size_type n) override
- {
- if(n<1) {
- throw std::out_of_range(
- "number of compartments in a segment must be one or more"
- );
- }
- num_compartments_ = n;
- }
-
- value_type radius(value_type loc) const
- {
- if(loc>=1.) return radii_.back();
- if(loc<=0.) return radii_.front();
-
- auto len = length();
- value_type pos = loc*len;
-
- // This could be cached using a partial sum.
- // In fact a lot of this stuff can be cached if
- // we find ourselves having to do it over and over again.
- // The time to cache it might be when update_lengths() is called.
- auto sum = value_type(0);
- size_type i = 0;
- for (i = 0; i<num_sub_segments(); ++i) {
- if(sum+lengths_[i]>pos) {
- break;
- }
- sum += lengths_[i];
- }
-
- auto rel = (len - sum)/lengths_[i];
-
- return rel*radii_[i] + (1.-rel)*radii_[i+1];
- }
-
-private:
- void update_lengths() {
- if (locations_.size()) {
- lengths_.resize(num_sub_segments());
- for (size_type i=0; i<num_sub_segments(); ++i) {
- lengths_[i] = norm(locations_[i] - locations_[i+1]);
- }
- }
- }
-
- size_type num_compartments_ = 1;
- std::vector<value_type> radii_;
- std::vector<value_type> lengths_;
- std::vector<point_type> locations_;
-};
-
-using segment_ptr = std::unique_ptr<segment>;
-
-/// Helper for constructing segments in a segment_ptr unique pointer wrapper.
-/// Forwards the supplied arguments to construct a segment of type SegmentType.
-/// e.g. auto my_cable = make_segment<cable>(section_kind::dendrite, ... );
-template <typename SegmentType, typename... Args>
-segment_ptr make_segment(Args&&... args) {
- return segment_ptr(new SegmentType(std::forward<Args>(args)...));
-}
-
-} // namespace arb
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index 715ddcc0e3b3dcb866cde30dcd2ddd53e282c3e6..1f7eafaa123195bda775cca63d6a4c638029b1ff 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -118,6 +118,29 @@ std::ostream& operator<<(std::ostream& o, const root_& x) {
return o << "(root)";
}
+// Proportional location on every branch.
+
+struct on_branches_ { double pos; };
+
+locset on_branches(double pos) {
+ return locset{on_branches_{pos}};
+}
+
+mlocation_list thingify_(const on_branches_& ob, const mprovider& p) {
+ msize_t n_branch = p.morphology().num_branches();
+
+ mlocation_list locs;
+ locs.reserve(n_branch);
+ for (msize_t b = 0; b<n_branch; ++b) {
+ locs.push_back({b, ob.pos});
+ }
+ return locs;
+}
+
+std::ostream& operator<<(std::ostream& o, const on_branches_& x) {
+ return o << "(on_branchs " << x.pos << ")";
+}
+
// Named locset.
struct named_: locset_tag {
diff --git a/arbor/util/piecewise.hpp b/arbor/util/piecewise.hpp
index a7a11ab029374f2321e6a5e6ee92762011d85877..4e4d42d1e7cd777a7abf27c8283524e67d7cdefa 100644
--- a/arbor/util/piecewise.hpp
+++ b/arbor/util/piecewise.hpp
@@ -146,6 +146,14 @@ struct pw_elements {
else return partn.index(x);
}
+ value_type operator()(double x) const {
+ size_type i = index_of(x);
+ if (i==npos) {
+ throw std::range_error("position outside support");
+ }
+ return (*this)[i];
+ }
+
// mutating operations:
void reserve(size_type n) {
@@ -411,8 +419,11 @@ auto zip(const pw_elements<A>& a, const pw_elements<B>& b, Combine combine = {})
if (rmin<lmax) return z;
double left = lmax;
- pw_size_type ai = a.intervals().index(left);
- pw_size_type bi = b.intervals().index(left);
+ pw_size_type ai = a.index_of(left);
+ pw_size_type bi = b.index_of(left);
+
+ arb_assert(ai!=(pw_size_type)-1);
+ arb_assert(bi!=(pw_size_type)-1);
if (rmin==left) {
z.push_back(left, left, combine(left, left, a[ai], b[bi]));
@@ -422,15 +433,17 @@ auto zip(const pw_elements<A>& a, const pw_elements<B>& b, Combine combine = {})
double a_right = a.interval(ai).second;
double b_right = b.interval(bi).second;
- while (left<rmin) {
+ for (;;) {
double right = std::min(a_right, b_right);
right = std::min(right, rmin);
z.push_back(left, right, combine(left, right, a[ai], b[bi]));
- if (a_right<=right) {
+ if (right==rmin) break;
+
+ if (a_right==right) {
a_right = a.interval(++ai).second;
}
- if (b_right<=right) {
+ if (b_right==right) {
b_right = b.interval(++bi).second;
}
diff --git a/arbor/util/rangeutil.hpp b/arbor/util/rangeutil.hpp
index 224d228753858837b083736f65a4a0bb43a8ce5e..4e45d8c01e7c306c6cdb0bb497a44502a517b679 100644
--- a/arbor/util/rangeutil.hpp
+++ b/arbor/util/rangeutil.hpp
@@ -305,6 +305,31 @@ bool is_sorted(const Seq& seq) {
return std::is_sorted(canon.begin(), canon.end());
}
+// Range version of std::equal.
+
+template <
+ typename Seq1, typename Seq2,
+ typename Eq = std::equal_to<void>,
+ typename = util::enable_if_sequence_t<const Seq1&>,
+ typename = util::enable_if_sequence_t<const Seq2&>
+>
+bool equal(const Seq1& seq1, const Seq2& seq2, Eq p = Eq{}) {
+ using std::begin;
+ using std::end;
+
+ auto i1 = begin(seq1);
+ auto e1 = end(seq1);
+
+ auto i2 = begin(seq2);
+ auto e2 = end(seq2);
+
+ while (i1!=e1 && i2!=e2) {
+ if (!p(*i1, *i2)) return false;
+ ++i1;
+ ++i2;
+ }
+ return i1==e1 && i2==e2;
+}
// Test if sequence is sorted after apply projection `proj` to elements.
// (TODO: this will perform unnecessary copies if `proj` returns a reference;
diff --git a/example/dryrun/branch_cell.hpp b/example/dryrun/branch_cell.hpp
index bbb818188a3694f4c3df060c6744c5a7246fc697..b4181c868811803391b61d8a9bf650534da68a8e 100644
--- a/example/dryrun/branch_cell.hpp
+++ b/example/dryrun/branch_cell.hpp
@@ -103,7 +103,7 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
d.set("soma", tagged(1));
d.set("dendrites", join(tagged(3), tagged(4)));
- arb::cable_cell cell(arb::morphology(tree, true), d, true);
+ arb::cable_cell cell(arb::morphology(tree, true), d);
cell.paint("soma", "hh");
cell.paint("dendrites", "pas");
@@ -120,5 +120,8 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
cell.place(arb::mlocation{1, 0.5}, "expsyn");
}
+ // Make a CV between every sample in the sample tree.
+ cell.default_parameters.discretization = arb::cv_policy_every_sample();
+
return cell;
}
diff --git a/example/dryrun/dryrun.cpp b/example/dryrun/dryrun.cpp
index ea21ea4ff3a1417672ca1ca78d381355cfc23fc4..d78df681a68538dcfed96aaf07f936e4eb51eff8 100644
--- a/example/dryrun/dryrun.cpp
+++ b/example/dryrun/dryrun.cpp
@@ -148,7 +148,6 @@ struct cell_stats {
size_type ncells = 0;
int nranks = 1;
size_type nsegs = 0;
- size_type ncomp = 0;
cell_stats(arb::recipe& r, run_params params) {
#ifdef ARB_MPI_ENABLED
@@ -161,14 +160,11 @@ struct cell_stats {
size_type b = rank*cells_per_rank;
size_type e = (rank+1)*cells_per_rank;
size_type nsegs_tmp = 0;
- size_type ncomp_tmp = 0;
for (size_type i=b; i<e; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs_tmp += c.num_branches();
- ncomp_tmp += c.num_compartments();
+ nsegs_tmp += c.morphology().num_branches();
}
MPI_Allreduce(&nsegs_tmp, &nsegs, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
- MPI_Allreduce(&ncomp_tmp, &ncomp, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
}
#else
if(!params.dry_run) {
@@ -176,8 +172,7 @@ struct cell_stats {
ncells = r.num_cells();
for (size_type i = 0; i < ncells; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs += c.num_branches();
- ncomp += c.num_compartments();
+ nsegs += c.morphology().num_branches();
}
}
#endif
@@ -187,12 +182,10 @@ struct cell_stats {
for (size_type i = 0; i < params.num_cells_per_rank; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs += c.num_branches();
- ncomp += c.num_compartments();
+ nsegs += c.morphology().num_branches();
}
nsegs *= params.num_ranks;
- ncomp *= params.num_ranks;
}
}
@@ -200,8 +193,7 @@ struct cell_stats {
return o << "cell stats: "
<< s.nranks << " ranks; "
<< s.ncells << " cells; "
- << s.nsegs << " branches; "
- << s.ncomp << " compartments.";
+ << s.nsegs << " branches. ";
}
};
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index 32c789d9950c31ea81f491fdc11808a4110e94ad..f59726a41bfc0c08da9c890fa0f6be06f1b97954 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -146,7 +146,6 @@ struct cell_stats {
using size_type = unsigned;
size_type ncells = 0;
size_type nsegs = 0;
- size_type ncomp = 0;
cell_stats(arb::recipe& r) {
#ifdef ARB_MPI_ENABLED
@@ -158,20 +157,16 @@ struct cell_stats {
size_type b = rank*cells_per_rank;
size_type e = (rank==nranks-1)? ncells: (rank+1)*cells_per_rank;
size_type nsegs_tmp = 0;
- size_type ncomp_tmp = 0;
for (size_type i=b; i<e; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
nsegs_tmp += c.num_branches();
- ncomp_tmp += c.num_compartments();
}
MPI_Allreduce(&nsegs_tmp, &nsegs, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
- MPI_Allreduce(&ncomp_tmp, &ncomp, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
#else
ncells = r.num_cells();
for (size_type i=0; i<ncells; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs += c.num_branches();
- ncomp += c.num_compartments();
+ nsegs += c.morphology().num_branches();
}
#endif
}
@@ -179,8 +174,7 @@ struct cell_stats {
friend std::ostream& operator<<(std::ostream& o, const cell_stats& s) {
return o << "cell stats: "
<< s.ncells << " cells; "
- << s.nsegs << " branchess; "
- << s.ncomp << " compartments.";
+ << s.nsegs << " branchess.";
}
};
diff --git a/example/ring/branch_cell.hpp b/example/ring/branch_cell.hpp
index bbb818188a3694f4c3df060c6744c5a7246fc697..b4181c868811803391b61d8a9bf650534da68a8e 100644
--- a/example/ring/branch_cell.hpp
+++ b/example/ring/branch_cell.hpp
@@ -103,7 +103,7 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
d.set("soma", tagged(1));
d.set("dendrites", join(tagged(3), tagged(4)));
- arb::cable_cell cell(arb::morphology(tree, true), d, true);
+ arb::cable_cell cell(arb::morphology(tree, true), d);
cell.paint("soma", "hh");
cell.paint("dendrites", "pas");
@@ -120,5 +120,8 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
cell.place(arb::mlocation{1, 0.5}, "expsyn");
}
+ // Make a CV between every sample in the sample tree.
+ cell.default_parameters.discretization = arb::cv_policy_every_sample();
+
return cell;
}
diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp
index 2220d6f4eb0ac84d8da60dcf2578b2f27b29bf3e..611ca77c9d5f54034a4ee697f017263896eaef5c 100644
--- a/example/ring/ring.cpp
+++ b/example/ring/ring.cpp
@@ -140,7 +140,6 @@ struct cell_stats {
using size_type = unsigned;
size_type ncells = 0;
size_type nsegs = 0;
- size_type ncomp = 0;
cell_stats(arb::recipe& r) {
#ifdef ARB_MPI_ENABLED
@@ -152,20 +151,16 @@ struct cell_stats {
size_type b = rank*cells_per_rank;
size_type e = (rank==nranks-1)? ncells: (rank+1)*cells_per_rank;
size_type nsegs_tmp = 0;
- size_type ncomp_tmp = 0;
for (size_type i=b; i<e; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
nsegs_tmp += c.num_branches();
- ncomp_tmp += c.num_compartments();
}
MPI_Allreduce(&nsegs_tmp, &nsegs, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
- MPI_Allreduce(&ncomp_tmp, &ncomp, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
#else
ncells = r.num_cells();
for (size_type i=0; i<ncells; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs += c.num_branches();
- ncomp += c.num_compartments();
+ nsegs += c.morphology().num_branches();
}
#endif
}
@@ -173,8 +168,7 @@ struct cell_stats {
friend std::ostream& operator<<(std::ostream& o, const cell_stats& s) {
return o << "cell stats: "
<< s.ncells << " cells; "
- << s.nsegs << " branches; "
- << s.ncomp << " compartments.";
+ << s.nsegs << " branches.";
}
};
diff --git a/example/single/single.cpp b/example/single/single.cpp
index 3e7eb01758e13d17f68d09d7560be4b86f2f97e3..d42f03b86ed0f906db515e3e7a8512eca07bd16a 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -20,6 +20,7 @@ struct options {
double t_end = 20;
double dt = 0.025;
float syn_weight = 0.01;
+ arb::cv_policy policy = arb::default_cv_policy();
};
options parse_options(int argc, char** argv);
@@ -27,8 +28,9 @@ arb::morphology default_morphology();
arb::morphology read_swc(const std::string& path);
struct single_recipe: public arb::recipe {
- explicit single_recipe(arb::morphology m): morpho(std::move(m)) {
+ explicit single_recipe(arb::morphology m, arb::cv_policy pol): morpho(std::move(m)) {
gprop.default_parameters = arb::neuron_parameter_defaults;
+ gprop.default_parameters.discretization = pol;
}
arb::cell_size_type num_cells() const override { return 1; }
@@ -58,26 +60,15 @@ struct single_recipe: public arb::recipe {
using arb::reg::tagged;
dict.set("soma", tagged(1));
dict.set("dend", join(tagged(3), tagged(4), tagged(42)));
- arb::cable_cell c(morpho, dict, false);
+ arb::cable_cell c(morpho, dict);
// Add HH mechanism to soma, passive channels to dendrites.
c.paint("soma", "hh");
c.paint("dend", "pas");
- // Discretize dendrites according to the NEURON d-lambda rule.
- /* skip this during refactoring of the morphology interface
- for (std::size_t i=1; i<c.num_branches(); ++i) {
- arb::cable_segment* branch = c.cable(i);
-
- double dx = c.segment_length_constant(100., i, gprop.default_parameters)*0.3;
- unsigned n = std::ceil(branch->length()/dx);
- branch->set_compartments(n);
- }
- */
-
// Add synapse to last branch.
- arb::cell_lid_type last_branch = c.num_branches()-1;
+ arb::cell_lid_type last_branch = c.morphology().num_branches()-1;
arb::mlocation end_last_branch = { last_branch, 1. };
c.place(end_last_branch, "exp2syn");
@@ -91,7 +82,7 @@ struct single_recipe: public arb::recipe {
int main(int argc, char** argv) {
try {
options opt = parse_options(argc, argv);
- single_recipe R(opt.swc_file.empty()? default_morphology(): read_swc(opt.swc_file));
+ single_recipe R(opt.swc_file.empty()? default_morphology(): read_swc(opt.swc_file), opt.policy);
auto context = arb::make_context();
arb::simulation sim(R, arb::partition_load_balance(R, context), context);
@@ -138,8 +129,11 @@ options parse_options(int argc, char** argv) {
else if (auto swc = parse_opt<std::string>(arg, 'm', "morphology")) {
opt.swc_file = swc.value();
}
+ else if (auto nseg = parse_opt<unsigned>(arg, 'n', "cv-per-branch")) {
+ opt.policy = arb::cv_policy_fixed_per_branch(nseg.value(), arb::cv_policy_flag::single_root_cv);
+ }
else {
- usage(argv[0], "[-m|--morphology SWCFILE] [-d|--dt TIME] [-t|--t-end TIME] [-w|--weight WEIGHT]");
+ usage(argv[0], "[-m|--morphology SWCFILE] [-d|--dt TIME] [-t|--t-end TIME] [-w|--weight WEIGHT] [-n|--cv-per-branch N]");
std::exit(1);
}
}
diff --git a/python/cells.cpp b/python/cells.cpp
index 8150dffe46313ff7f01c2d56ffb39bb81a4512d1..6eaa3cb09f5caf1133c56913b0ff75998360c1f4 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -138,7 +138,7 @@ arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& p
d.set("soma", tagged(1));
d.set("dendrites", join(tagged(3), tagged(4)));
- arb::cable_cell cell(arb::morphology(tree, true), d, true);
+ arb::cable_cell cell(arb::morphology(tree, true), d);
cell.paint("soma", "hh");
cell.paint("dendrites", "pas");
diff --git a/test/common_cells.hpp b/test/common_cells.hpp
index 02821167bc0c550413c8b763bf4c2af1fbc5468f..83d972ae9ad7f953df0a4f14dc662196488ceebe 100644
--- a/test/common_cells.hpp
+++ b/test/common_cells.hpp
@@ -1,7 +1,6 @@
#include <cmath>
#include <arbor/cable_cell.hpp>
-#include <arbor/segment.hpp>
#include <arbor/mechinfo.hpp>
#include <arbor/morph/label_dict.hpp>
#include <arbor/recipe.hpp>
@@ -13,6 +12,7 @@ class soma_cell_builder {
sample_tree tree;
std::vector<msize_t> branch_distal_id;
std::unordered_map<std::string, int> tag_map;
+ locset cv_boundaries = mlocation{0, 1.};
int tag_count = 0;
// Get tag id of region.
@@ -56,7 +56,11 @@ public:
p = tree.append(p, {{0,0,z+len,r2}, tag});
branch_distal_id.push_back(p);
- return branch_distal_id.size()-1;
+ msize_t bid = branch_distal_id.size()-1;
+ for (int i = 0; i<ncomp; ++i) {
+ cv_boundaries = sum(cv_boundaries, mlocation{bid, (2*i+1.)/(2.*ncomp)});
+ }
+ return bid;
}
cable_cell make_cell() const {
@@ -79,9 +83,9 @@ public:
}
// Make cable_cell from sample tree and dictionary.
- // The true flag is used to force the discretization to make compartments
- // at sample points.
- return cable_cell(tree, dict, true);
+ cable_cell c(tree, dict);
+ c.default_parameters.discretization = cv_policy_explicit(cv_boundaries);
+ return c;
}
};
@@ -146,7 +150,7 @@ inline cable_cell make_cell_ball_and_stick(bool with_stim = true) {
}
/*
- * Create cell with a soma and three-segment dendrite with single branch point:
+ * Create cell with a soma and three-branch dendrite with single branch point:
*
* O----======
*
diff --git a/test/ubench/CMakeLists.txt b/test/ubench/CMakeLists.txt
index 21e67ca9f46789d1ebea5843ac796631fcf79db6..bfc4805ad7f831d819d1955163ffe5cadd377615 100644
--- a/test/ubench/CMakeLists.txt
+++ b/test/ubench/CMakeLists.txt
@@ -7,7 +7,7 @@ set(bench_sources
default_construct.cpp
event_setup.cpp
event_binning.cpp
- mech_vec.cpp
+ # mech_vec.cpp --- requires rework for new API.
task_system.cpp
)
diff --git a/test/ubench/mech_vec.cpp b/test/ubench/mech_vec.cpp
index 1691002a72c7a9c4d9c15603100c2d02249e879a..157d044fabe0cf2cb6aa23ee8e7a2d6a4a70c3dc 100644
--- a/test/ubench/mech_vec.cpp
+++ b/test/ubench/mech_vec.cpp
@@ -2,6 +2,9 @@
//
// Start with pas (passive dendrite) mechanism
+// NOTE: This targets an earlier version of the Arbor API and
+// will need to be reworked in order to compile.
+
#include <fstream>
#include <arbor/cable_cell.hpp>
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 07d3d52fcca410e6dd3fc76ce7b9336f2fd9f5b1..f851b7d01190d3f3478e58f96a1a5fc13563ce15 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -80,9 +80,9 @@ set(unit_sources
test_algorithms.cpp
test_any.cpp
test_backend.cpp
- test_cable_cell.cpp
- test_compartments.cpp
test_counter.cpp
+ test_cv_geom.cpp
+ test_cv_layout.cpp
test_cv_policy.cpp
test_cycle.cpp
test_domain_decomposition.cpp
@@ -127,7 +127,6 @@ set(unit_sources
test_probe.cpp
test_range.cpp
test_ratelem.cpp
- test_segment.cpp
test_schedule.cpp
test_spike_source.cpp
test_scope_exit.cpp
diff --git a/test/unit/common.hpp b/test/unit/common.hpp
index 46e3f5cd0fc867a17bace5209e152e4202b79c9f..2436caa3866361183ecdda7116e56b4df45f5b2b 100644
--- a/test/unit/common.hpp
+++ b/test/unit/common.hpp
@@ -11,6 +11,15 @@
#include "../gtest.h"
+// Pair printer.
+
+namespace std {
+ template <typename A, typename B>
+ std::ostream& operator<<(std::ostream& out, const std::pair<A, B>& p) {
+ return out << '(' << p.first << ',' << p.second << ')';
+ }
+}
+
namespace testing {
// Sentinel for C-style strings, for use with range-related tests.
diff --git a/test/unit/common_morphologies.hpp b/test/unit/common_morphologies.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..666fd0d88d43c4be9787e0dc0d5e6f42bec562c8
--- /dev/null
+++ b/test/unit/common_morphologies.hpp
@@ -0,0 +1,47 @@
+#pragma once
+
+// A set of morphologies for testing discretization.
+
+#include <utility>
+#include <vector>
+#include <arbor/morph/morphology.hpp>
+
+namespace common_morphology {
+
+inline std::vector<arb::msample> make_morph_samples(unsigned n) {
+ std::vector<arb::msample> ms;
+ for (unsigned i = 0; i<n; ++i) ms.push_back({{0., 0., (double)i, 0.5}, 5});
+ return ms;
+}
+
+// Test morphologies for CV determination:
+// Sample points have radius 0.5, giving an initial branch length of 1.0
+// for morphologies with spherical roots.
+
+static const arb::morphology m_empty;
+
+// spherical root, one branch
+static const arb::morphology m_sph_b1{arb::sample_tree(make_morph_samples(1), {arb::mnpos}), true};
+
+// regular root, one branch
+static const arb::morphology m_reg_b1{arb::sample_tree(make_morph_samples(2), {arb::mnpos, 0u}), false};
+
+// spherical root, six branches
+static const arb::morphology m_sph_b6{arb::sample_tree(make_morph_samples(8), {arb::mnpos, 0u, 1u, 0u, 3u, 4u, 4u, 4u}), true};
+
+// regular root, six branches
+static const arb::morphology m_reg_b6{arb::sample_tree(make_morph_samples(7), {arb::mnpos, 0u, 1u, 1u, 2u, 2u, 2u}), false};
+
+// regular root, six branches, mutiple top level branches.
+static const arb::morphology m_mlt_b6{arb::sample_tree(make_morph_samples(7), {arb::mnpos, 0u, 1u, 1u, 0u, 4u, 4u}), false};
+
+static std::pair<const char*, arb::morphology> test_morphologies[] = {
+ {"m_empty", m_empty},
+ {"m_sph_b1", m_sph_b1},
+ {"m_reg_b1", m_reg_b1},
+ {"m_sph_b6", m_sph_b6},
+ {"m_reg_b6", m_reg_b6},
+ {"m_mlt_b6", m_mlt_b6}
+};
+
+} // namespace common_morphology
diff --git a/test/unit/test_cable_cell.cpp b/test/unit/test_cable_cell.cpp
deleted file mode 100644
index 62c54706886dee11cf1b5ba829e9c9ede03fb1b1..0000000000000000000000000000000000000000
--- a/test/unit/test_cable_cell.cpp
+++ /dev/null
@@ -1,81 +0,0 @@
-#include "../gtest.h"
-
-#include <arbor/cable_cell.hpp>
-#include <arbor/math.hpp>
-#include <arbor/morph/locset.hpp>
-
-#include "io/sepval.hpp"
-
-#include "tree.hpp"
-
-using namespace arb;
-using ::arb::math::pi;
-
-TEST(cable_cell, soma) {
- // test that insertion of a soma works
- // define with centre point @ (0,0,1)
- double soma_radius = 3.2;
-
- arb::sample_tree samples;
- samples.append({0,0,0,soma_radius,1});
- auto c = cable_cell(arb::morphology(samples));
-
- EXPECT_EQ(c.has_soma(), true);
-
- auto s = c.soma();
- EXPECT_EQ(s->radius(), soma_radius);
-}
-
-TEST(cable_cell, multiple_cables) {
- double soma_radius = std::pow(3./(4.*pi<double>), 1./3.);
-
- // Generate a cylindrical cable segment of length 1/pi and radius 1
- // volume = 1
- // area = 2
- // Returns the distal point of the added cable.
- auto append_branch = [soma_radius](sample_tree& stree, int proximal) {
- constexpr int tag = 2;
- if (!proximal) {
- double z = soma_radius;
- proximal = stree.append(0, {0,0,z, 1/pi<double>, tag});
- }
- return stree.append(proximal, msample{0, 0, stree.samples()[proximal].loc.z+1, 1/pi<double>, tag});
- };
-
- // cell strucure with branch numbers
- //
- // 0
- // / \.
- // 1 2
- // / \.
- // 3 4
-
- arb::sample_tree samples;
- samples.append({0,0,-soma_radius,soma_radius,1});
-
- // hook the dendrite and axons
- append_branch(samples, 0);
- append_branch(samples, 0);
- append_branch(samples, 2);
- append_branch(samples, 2);
-
- auto c = cable_cell(arb::morphology(samples, true));
-
- EXPECT_EQ(c.num_branches(), 5u);
-
- // construct the graph
- tree con(c.parents());
-
- auto no_parent = tree::no_parent;
- EXPECT_EQ(con.num_segments(), 5u);
- EXPECT_EQ(con.parent(0), no_parent);
- EXPECT_EQ(con.parent(1), 0u);
- EXPECT_EQ(con.parent(2), 0u);
- EXPECT_EQ(con.parent(3), 1u);
- EXPECT_EQ(con.parent(4), 1u);
- EXPECT_EQ(con.num_children(0), 2u);
- EXPECT_EQ(con.num_children(1), 2u);
- EXPECT_EQ(con.num_children(2), 0u);
- EXPECT_EQ(con.num_children(3), 0u);
- EXPECT_EQ(con.num_children(4), 0u);
-}
diff --git a/test/unit/test_compartments.cpp b/test/unit/test_compartments.cpp
deleted file mode 100644
index deb5996ee18d806e63a4d26fa9d6552af3fbc136..0000000000000000000000000000000000000000
--- a/test/unit/test_compartments.cpp
+++ /dev/null
@@ -1,259 +0,0 @@
-#include <cmath>
-#include <string>
-
-#include "../gtest.h"
-
-#include <arbor/math.hpp>
-
-#include "fvm_compartment.hpp"
-#include "util/rangeutil.hpp"
-#include "util/span.hpp"
-#include "util/transform.hpp"
-
-using namespace arb;
-using namespace arb::math;
-
-using arb::util::make_span;
-using arb::util::transform_view;
-using arb::util::sum;
-
-// Divided compartments
-// (FVM-friendly compartment data)
-
-template <std::size_t N>
-struct pw_cable_data {
- double radii[N+1];
- double lengths[N];
-
- static constexpr std::size_t nseg() { return N; }
- double r1() const { return radii[0]; }
- double r2() const { return radii[N]; }
- double length() const { return sum(lengths); }
-
- double area() const {
- return sum(transform_view(make_span(N),
- [&](unsigned i) { return area_frustrum(lengths[i], radii[i], radii[i+1]); }));
- }
-
- double volume() const {
- return sum(transform_view(make_span(N),
- [&](unsigned i) { return volume_frustrum(lengths[i], radii[i], radii[i+1]); }));
- }
-};
-
-pw_cable_data<1> cable_one = {
- {2.0, 5.0},
- {10.0}
-};
-
-pw_cable_data<4> cable_linear = {
- {2.0, 3.5, 6.0, 6.5, 6.75},
- {3.0, 5.0, 1.0, 0.5}
-};
-
-pw_cable_data<4> cable_jumble = {
- {2.0, 6.0, 3.5, 6.75, 6.5},
- {3.0, 5.0, 1.0, 0.5}
-};
-
-void expect_equal_divs(const div_compartment& da, const div_compartment& db) {
- EXPECT_EQ(da.index, db.index);
-
- double eps = std::numeric_limits<double>::epsilon();
- double e1 = std::min(da.length(), db.length())*8*eps;
- double e2 = std::min(da.area(), db.area())*8*eps;
- double e3 = std::min(da.volume(), db.volume())*8*eps;
-
- EXPECT_NEAR(da.left.length, db.left.length, e1);
- EXPECT_NEAR(da.left.area, db.left.area, e2);
- EXPECT_NEAR(da.left.volume, db.left.volume, e3);
- EXPECT_NEAR(da.left.radii.first, db.left.radii.first, e1);
- EXPECT_NEAR(da.left.radii.second, db.left.radii.second, e1);
-
- EXPECT_NEAR(da.right.length, db.right.length, e1);
- EXPECT_NEAR(da.right.area, db.right.area, e2);
- EXPECT_NEAR(da.right.volume, db.right.volume, e3);
- EXPECT_NEAR(da.right.radii.first, db.right.radii.first, e1);
- EXPECT_NEAR(da.right.radii.second, db.right.radii.second, e1);
-}
-
-TEST(compartments, div_ends) {
- {
- div_compartment_by_ends divcomps{1, cable_one.radii, cable_one.lengths};
-
- auto d = divcomps(0);
- auto r1 = cable_one.radii[0];
- auto r2 = cable_one.radii[1];
- auto l = cable_one.lengths[0];
-
- EXPECT_DOUBLE_EQ(r1, d.radii().first);
- EXPECT_DOUBLE_EQ(r2, d.radii().second);
- EXPECT_DOUBLE_EQ(l, d.length());
- EXPECT_DOUBLE_EQ(area_frustrum(l, r2, r1), d.area());
- EXPECT_DOUBLE_EQ(volume_frustrum(l, r2, r1), d.volume());
-
- auto sl = l/2.0;
- auto rc = 0.5*(r1+r2);
-
- div_compartment expected{
- 0,
- semi_compartment{sl, area_frustrum(sl, r1, rc), volume_frustrum(sl, r1, rc), {r1, rc}},
- semi_compartment{sl, area_frustrum(sl, rc, r2), volume_frustrum(sl, rc, r2), {rc, r2}}
- };
-
- SCOPED_TRACE("cable_one");
- expect_equal_divs(expected, d);
- }
-
- {
- // for a linear cable, expect this compartment maker to
- // create consistent compartments
-
- constexpr unsigned ncomp = 7;
- div_compartment_by_ends divlin{ncomp, cable_linear.radii, cable_linear.lengths};
-
- auto r1 = cable_linear.r1();
- auto r2 = cable_linear.r2();
- auto l = cable_linear.length();
-
- pw_cable_data<1> one = { {r1, r2}, {l} };
- div_compartment_by_ends divone{ncomp, one.radii, one.lengths};
-
- for (unsigned i=0; i<ncomp; ++i) {
- SCOPED_TRACE("cable_linear compartment "+std::to_string(i));
- auto da = divlin(i);
- auto db = divone(i);
-
- EXPECT_DOUBLE_EQ(l/ncomp, da.length());
- expect_equal_divs(da, db);
- }
- }
-}
-
-TEST(compartments, div_sample) {
- // expect by_ends and sampler to give same results on linear cable
- {
- constexpr unsigned ncomp = 7;
- div_compartment_sampler divsampler{ncomp, cable_linear.radii, cable_linear.lengths};
- div_compartment_by_ends divends{ncomp, cable_linear.radii, cable_linear.lengths};
-
- auto l = cable_linear.length();
- for (unsigned i=0; i<ncomp; ++i) {
- SCOPED_TRACE("cable_linear compartment "+std::to_string(i));
- auto da = divsampler(i);
- auto db = divends(i);
-
- EXPECT_DOUBLE_EQ(l/ncomp, da.length());
- expect_equal_divs(da, db);
- }
- }
-
- // expect (up to rounding) correct total area and volume if compartments
- // align with sub-segments; when they don't align, expect error to decrease
- // with ncomp
- {
- double area_expected = cable_jumble.area();
- double volume_expected = cable_jumble.volume();
- double eps = std::numeric_limits<double>::epsilon();
-
- double common_dx = 0.5; // depends on cable_jumble.lengths;
- // check our common_dx actually is ok
- for (double l: cable_jumble.lengths) {
- ASSERT_DOUBLE_EQ(l/common_dx, std::round(l/common_dx));
- }
-
- double length = cable_jumble.length();
- unsigned nbase = std::round(length/common_dx);
-
- for (unsigned m: {1u, 3u, 7u}) {
- unsigned ncomp = m*nbase;
- div_compartment_sampler divs{ncomp, cable_jumble.radii, cable_jumble.lengths};
-
- double area = sum(transform_view(make_span(ncomp),
- [&](unsigned i) { return divs(i).area(); }));
-
- double volume = sum(transform_view(make_span(ncomp),
- [&](unsigned i) { return divs(i).volume(); }));
-
- double e2 = std::min(area, area_expected)*ncomp*eps;
- double e3 = std::min(volume, volume_expected)*ncomp*eps;
-
- SCOPED_TRACE("cable_jumble ncomp "+std::to_string(ncomp));
-
- EXPECT_NEAR(area_expected, area, e2);
- EXPECT_NEAR(volume_expected, volume, e3);
- }
-
- double coarse_area = area_frustrum(length, cable_jumble.r1(), cable_jumble.r2());
- double coarse_volume = volume_frustrum(length, cable_jumble.r1(), cable_jumble.r2());
- double area_error = std::abs(area_expected-coarse_area);
- double volume_error = std::abs(volume_expected-coarse_volume);
-
- for (unsigned m: {1u, 10u, 100u}) {
- unsigned ncomp = m*nbase+1u;
- div_compartment_sampler divs{ncomp, cable_jumble.radii, cable_jumble.lengths};
-
- double area = sum(transform_view(make_span(ncomp),
- [&](unsigned i) { return divs(i).area(); }));
-
- double volume = sum(transform_view(make_span(ncomp),
- [&](unsigned i) { return divs(i).volume(); }));
-
- SCOPED_TRACE("cable_jumble ncomp "+std::to_string(ncomp));
-
- double err = std::abs(area_expected-area);
- EXPECT_LT(err, area_error);
- area_error = err;
-
- err = std::abs(volume_expected-volume);
- EXPECT_LT(err, volume_error);
- volume_error = err;
- }
- }
-}
-
-TEST(compartments, div_integrator) {
- // expect integrator and sampler to give same results on linear cable
- {
- constexpr unsigned ncomp = 7;
- div_compartment_sampler divintegrator{ncomp, cable_linear.radii, cable_linear.lengths};
- div_compartment_by_ends divends{ncomp, cable_linear.radii, cable_linear.lengths};
-
- auto l = cable_linear.length();
- for (unsigned i=0; i<ncomp; ++i) {
- SCOPED_TRACE("cable_linear compartment "+std::to_string(i));
- auto da = divintegrator(i);
- auto db = divends(i);
-
- EXPECT_DOUBLE_EQ(l/ncomp, da.length());
- expect_equal_divs(da, db);
- }
- }
-
- // expect integrator to give same (up to rounding) total areas and volumes
- // as the cable
- {
- double area_expected = cable_jumble.area();
- double volume_expected = cable_jumble.volume();
- double eps = std::numeric_limits<double>::epsilon();
-
- for (unsigned ncomp: make_span(1u, 23u)) {
- div_compartment_integrator divs{ncomp, cable_jumble.radii, cable_jumble.lengths};
-
- double area = sum(transform_view(make_span(ncomp),
- [&](unsigned i) { return divs(i).area(); }));
-
- double volume = sum(transform_view(make_span(ncomp),
- [&](unsigned i) { return divs(i).volume(); }));
-
- double e2 = std::min(area, area_expected)*ncomp*eps;
- double e3 = std::min(volume, volume_expected)*ncomp*eps;
-
- SCOPED_TRACE("cable_jumble ncomp "+std::to_string(ncomp));
-
- EXPECT_NEAR(area_expected, area, e2);
- EXPECT_NEAR(volume_expected, volume, e3);
- }
- }
-}
-
diff --git a/test/unit/test_cv_geom.cpp b/test/unit/test_cv_geom.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a14b007e55e503663c8cd0cfed5c1eddadab542
--- /dev/null
+++ b/test/unit/test_cv_geom.cpp
@@ -0,0 +1,323 @@
+#include <algorithm>
+#include <utility>
+
+#include <arbor/util/optional.hpp>
+#include <arbor/cable_cell.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/locset.hpp>
+
+#include "fvm_layout.hpp"
+#include "util/rangeutil.hpp"
+
+#include "common.hpp"
+#include "common_morphologies.hpp"
+#include "../common_cells.hpp"
+
+using namespace arb;
+using util::make_span;
+
+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());
+
+ EXPECT_TRUE(geom.cv_parent.empty());
+ EXPECT_TRUE(geom.cv_cables.empty());
+ EXPECT_TRUE(geom.cv_cables_divs.empty());
+
+ EXPECT_EQ(0u, geom.size()); // size()/empty() reflects number of CVs.
+ EXPECT_EQ(1u, geom.n_cell()); // can have no CVs but >0 cells.
+}
+
+TEST(cv_geom, trivial) {
+ using namespace common_morphology;
+
+ for (auto& p: test_morphologies) {
+ if (p.second.empty()) continue;
+
+ SCOPED_TRACE(p.first);
+ cable_cell cell{p.second};
+ 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());
+
+ EXPECT_EQ(1u, geom1.size());
+ 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()));
+
+ EXPECT_EQ(geom3.cv_cables, geom4.cv_cables);
+ if (m.branch_children(mnpos).size()==1) {
+ EXPECT_EQ(geom1.cv_cables, geom4.cv_cables);
+ }
+
+ mcable_list all_cables = thingify(reg::all(), cell.provider());
+ EXPECT_TRUE(testing::seq_eq(all_cables, geom1.cables(0)));
+ }
+}
+
+TEST(cv_geom, one_cv_per_branch) {
+ using namespace common_morphology;
+
+ for (auto& p: test_morphologies) {
+ if (p.second.empty()) continue;
+ SCOPED_TRACE(p.first);
+
+ cable_cell cell{p.second};
+ auto& m = cell.morphology();
+
+ cv_geometry geom = cv_geometry_from_ends(cell, ls::on_branches(0));
+
+ // Expect trivial CVs at every fork point, and single-cable CVs for each branch.
+ std::vector<unsigned> seen_branches(m.num_branches(), 0);
+ auto n_branch_child = [&m](msize_t b) { return m.branch_children(b).size(); };
+ for (auto i: make_span(geom.size())) {
+ auto cables = geom.cables(i);
+
+ ASSERT_EQ(1u, cables.size());
+ auto c = cables.front();
+
+ if (c.prox_pos==c.dist_pos) {
+ if (c.branch==0 && c.prox_pos==0) {
+ EXPECT_TRUE(n_branch_child(mnpos)>1);
+ }
+ else {
+ EXPECT_EQ(1., c.prox_pos);
+ EXPECT_TRUE(n_branch_child(c.branch)>1);
+ }
+ }
+ else {
+ ++seen_branches[c.branch];
+ EXPECT_EQ(1., seen_branches[c.branch]);
+ EXPECT_EQ(0., c.prox_pos);
+ EXPECT_EQ(1., c.dist_pos);
+
+ // Confirm parent CV is fork CV:
+ if (i>0) {
+ mlocation pfork = canonical(m, mlocation{c.branch, 0.});
+
+ auto pcables = geom.cables(geom.cv_parent[i]);
+ ASSERT_EQ(1u, pcables.size());
+
+ mcable p = pcables.front();
+ EXPECT_EQ(pfork.branch, p.branch);
+ EXPECT_EQ(p.prox_pos, p.dist_pos);
+
+ if (p.branch==0) {
+ EXPECT_TRUE(p.prox_pos==0 || p.prox_pos==1);
+ }
+ else {
+ EXPECT_EQ(1., p.prox_pos);
+ }
+ }
+ }
+ }
+
+ EXPECT_TRUE(std::find(seen_branches.begin(), seen_branches.end(), 0)==seen_branches.end());
+ }
+}
+
+TEST(cv_geom, midpoints) {
+ using namespace common_morphology;
+
+ // Place CV boundaries at the midpoints of each branch.
+ for (auto& p: test_morphologies) {
+ if (p.second.empty()) continue;
+ SCOPED_TRACE(p.first);
+
+ cable_cell cell{p.second};
+ auto& m = cell.morphology();
+
+ cv_geometry geom = cv_geometry_from_ends(cell, ls::on_branches(0.5));
+
+ // Expect CVs to be either: covering fork points, with one cable per branch
+ // at the fork (for a multiple-root-branch morphology, this would be treating
+ // (0, 0) as a fork); or the last halves of terminal branches or the first half
+ // of a unique root branch.
+
+ auto n_branch_child = [&m](msize_t b) { return m.branch_children(b).size(); };
+ for (auto i: make_span(geom.size())) {
+ auto cables = geom.cables(i);
+
+ if (i==0) {
+ // Expect inital half of single branch cell, or branched CV around (0,0).
+ if (cables.size()==1) {
+ EXPECT_EQ(1u, n_branch_child(mnpos));
+ auto c = cables.front();
+ EXPECT_EQ(0u, c.branch);
+ EXPECT_EQ(0.0, c.prox_pos);
+ EXPECT_EQ(0.5, c.dist_pos);
+ }
+ else {
+ EXPECT_TRUE(n_branch_child(mnpos)>1);
+ for (auto& c: cables) {
+ auto x = canonical(m, mlocation{c.branch, 0.});
+ EXPECT_EQ(0u, x.branch);
+
+ EXPECT_EQ(0.0, c.prox_pos);
+ EXPECT_EQ(0.5, c.dist_pos);
+ }
+ }
+ }
+ else {
+ // Expect final half of terminal branch or a branched CV around an interior fork.
+ if (cables.size()==1) {
+ // Terminal segment, or initial segment of 1-branch cell.
+ auto c = cables.front();
+ EXPECT_EQ(0.5, c.prox_pos);
+ EXPECT_EQ(1.0, c.dist_pos);
+ EXPECT_EQ(0u, n_branch_child(c.branch));
+ }
+ else {
+ auto prox_cable = cables.front();
+ EXPECT_EQ(0.5, prox_cable.prox_pos);
+ EXPECT_EQ(1.0, prox_cable.dist_pos);
+
+ msize_t prox_branch = prox_cable.branch;
+ EXPECT_EQ(1+n_branch_child(prox_branch), cables.size());
+
+ for (unsigned j = 1; j<cables.size(); ++j) {
+ auto& c = cables[j];
+ EXPECT_EQ(0.0, c.prox_pos);
+ EXPECT_EQ(0.5, c.dist_pos);
+ auto x = canonical(m, mlocation{c.branch, 0.});
+ EXPECT_EQ(prox_branch, x.branch);
+ }
+ }
+ }
+ }
+ }
+}
+
+TEST(cv_geom, weird) {
+ // m_reg_b6 has the following branch structure:
+ //
+ // ---0---+---1---+---3---
+ // | |
+ // | +---4---
+ // 2 |
+ // | +---5---
+ // |
+ //
+ // By placing CV boundary points at (1,0) and (4,0), we
+ // should obtain 3 CVs 'o', '+' and '=' as:
+ //
+ //
+ // oooooooo+++++++++++++++
+ // o +
+ // o +=======
+ // o +
+ // o ++++++++
+ // o
+ //
+ // CV 0 will comprise branches 0 and 2; CV 1 branches 1, 3, 5;
+ // and CV 2 branch 4.
+
+ using C = mcable;
+ 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}});
+
+ ASSERT_EQ(3u, geom.size());
+
+ mcable_list expected0 = {C{0u, 0., 1.}, C{2u, 0., 1.}};
+ EXPECT_TRUE(seq_eq(expected0, geom.cables(0)));
+
+ mcable_list expected1 = {C{1u, 0., 1.}, C{3u, 0., 1.}, C{5u, 0., 1.}};
+ EXPECT_TRUE(seq_eq(expected1, geom.cables(1)));
+
+ mcable_list expected2 = {C{4u, 0., 1.}};
+ EXPECT_TRUE(seq_eq(expected2, geom.cables(2)));
+}
+
+TEST(cv_geom, location_cv) {
+ using namespace common_morphology;
+
+ cable_cell cell{m_reg_b6};
+
+ // Two CVs per branch, plus trivial CV at forks.
+ cv_geometry geom = cv_geometry_from_ends(cell,
+ join(ls::on_branches(0.0), ls::on_branches(0.5), ls::on_branches(1.)));
+
+ // Confirm CVs are all only one cable, and either trivial or half a branch.
+ for (auto cv: geom.cell_cvs(0)) {
+ auto cables = geom.cables(cv);
+ ASSERT_EQ(1u, cables.size());
+
+ mcable cable = cables.front();
+ ASSERT_TRUE(cable.prox_pos==cable.dist_pos ||
+ (cable.prox_pos==0 && cable.dist_pos==0.5 )||
+ (cable.prox_pos==0.5 && cable.dist_pos==1.));
+ }
+
+ // Where ambiguous, CV found should be one with a non-trivial cable containing the location.
+ for (auto bid: util::make_span(cell.morphology().num_branches())) {
+ for (double pos: {0., 0.3, 0.5, 0.7, 1.}) {
+ mlocation loc{bid, pos};
+ SCOPED_TRACE(loc);
+
+ auto cv = geom.location_cv(0, loc);
+
+ mcable cable = geom.cables(cv).front();
+ EXPECT_TRUE(cable.prox_pos<cable.dist_pos);
+ EXPECT_TRUE(cable.prox_pos<=pos);
+ EXPECT_TRUE(cable.dist_pos>=pos);
+ EXPECT_TRUE(cable.branch==loc.branch);
+ }
+ }
+
+ // CV 0 will comprise branches 0 and 2; CV 1 branches 1, 3, 5;
+ // and CV 2 branch 4 (see test cv_geom.weird above).
+ cv_geometry gweird = cv_geometry_from_ends(cell, mlocation_list{{1, 0}, {4,0}});
+
+ // Expect location (4, 0.) to be in CV 2, but colocated locations
+ // (3, 0.), (5, 0.) and (1, 1.) to be in CV 1.
+ EXPECT_EQ(2u, gweird.location_cv(0, mlocation{4, 0.}));
+ EXPECT_EQ(1u, gweird.location_cv(0, mlocation{3, 0.}));
+ EXPECT_EQ(1u, gweird.location_cv(0, mlocation{5, 0.}));
+ EXPECT_EQ(1u, gweird.location_cv(0, mlocation{1, 1.}));
+}
+
+TEST(cv_geom, multicell) {
+ using namespace common_morphology;
+ using index_type = cv_geometry::index_type;
+
+ cable_cell cell = cable_cell{m_reg_b6};
+
+ cv_geometry geom = cv_geometry_from_ends(cell, ls::on_branches(0.5));
+ unsigned n_cv = geom.size();
+
+ cv_geometry geom2 = geom;
+ append(geom2, geom);
+
+ ASSERT_EQ(2*n_cv, geom2.size());
+ for (unsigned i = 0; i<n_cv; ++i) {
+ EXPECT_EQ(geom.cv_parent[i], geom2.cv_parent[i]);
+
+ if (geom2.cv_parent[i]==-1) {
+ EXPECT_EQ(-1, geom2.cv_parent[i+n_cv]);
+ }
+ else {
+ EXPECT_EQ(geom2.cv_parent[i]+(int)n_cv, geom2.cv_parent[i+n_cv]);
+ }
+ EXPECT_EQ(0, geom2.cv_to_cell[i]);
+ EXPECT_EQ(1, geom2.cv_to_cell[i+n_cv]);
+
+ mcable_list cables, cables1, cables2;
+ util::assign(cables, geom.cables(i));
+ util::assign(cables1, geom2.cables(i));
+ util::assign(cables2, geom2.cables(i+n_cv));
+ EXPECT_EQ(cables, cables1);
+ EXPECT_EQ(cables, cables2);
+ }
+
+ 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));
+}
diff --git a/test/unit/test_cv_layout.cpp b/test/unit/test_cv_layout.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4d5f9710c8b51f1d98a1e40f86be6a4042777783
--- /dev/null
+++ b/test/unit/test_cv_layout.cpp
@@ -0,0 +1,201 @@
+#include <algorithm>
+#include <utility>
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/math.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/util/optional.hpp>
+
+#include "fvm_layout.hpp"
+#include "util/span.hpp"
+
+#include "common.hpp"
+#include "common_morphologies.hpp"
+#include "../common_cells.hpp"
+
+using namespace arb;
+using util::make_span;
+
+TEST(cv_layout, empty) {
+ using namespace common_morphology;
+
+ cable_cell empty_cell{m_empty};
+ fvm_cv_discretization D = fvm_cv_discretize(empty_cell, neuron_parameter_defaults);
+
+ EXPECT_TRUE(D.empty());
+ EXPECT_EQ(0u, D.size());
+ EXPECT_EQ(1u, D.n_cell());
+
+ EXPECT_EQ(0u, D.face_conductance.size());
+ EXPECT_EQ(0u, D.cv_area.size());
+ EXPECT_EQ(0u, D.cv_capacitance.size());
+ EXPECT_EQ(0u, D.init_membrane_potential.size());
+ EXPECT_EQ(0u, D.temperature_K.size());
+ EXPECT_EQ(0u, D.diam_um.size());
+}
+
+TEST(cv_layout, trivial) {
+ using namespace common_morphology;
+
+ auto params = neuron_parameter_defaults;
+ params.discretization = cv_policy_explicit(ls::nil());
+
+ // For each cell, check size, confirm area is morphological area from
+ // embedding, and that membrane-properties are equal to defaults.
+
+ std::vector<cable_cell> cells;
+ unsigned n_cv = 0;
+ for (auto& p: test_morphologies) {
+ cells.emplace_back(p.second);
+ n_cv += !p.second.empty(); // one cv per non-empty cell
+ }
+
+ auto n_cells = cells.size();
+ fvm_cv_discretization D = fvm_cv_discretize(cells, params);
+
+ EXPECT_EQ(n_cv, D.size());
+ for (unsigned i = 0; i<n_cells; ++i) {
+ auto cv_indices = util::make_span(D.geometry.cell_cv_interval(i));
+ if (test_morphologies[i].second.empty()) {
+ ASSERT_TRUE(cv_indices.empty());
+ continue;
+ }
+ else {
+ ASSERT_EQ(1u, cv_indices.size());
+ }
+
+ auto cv = cv_indices.front();
+
+ EXPECT_DOUBLE_EQ(params.temperature_K.value(), D.temperature_K[cv]);
+ EXPECT_DOUBLE_EQ(params.init_membrane_potential.value(), D.init_membrane_potential[cv]);
+
+ double total_area = 0;
+ unsigned n_branch = cells[i].morphology().num_branches();
+ const auto& embedding = cells[i].embedding();
+ for (unsigned b = 0; b<n_branch; ++b) {
+ total_area += embedding.integrate_area(mcable{b, 0., 1.});
+ }
+
+ EXPECT_DOUBLE_EQ(total_area, D.cv_area[cv]);
+ EXPECT_DOUBLE_EQ(total_area*params.membrane_capacitance.value(), D.cv_capacitance[cv]);
+ }
+}
+
+TEST(cv_layout, cable) {
+ auto morph = common_morphology::m_reg_b1; // one branch, cable constant radius.
+
+ auto params = neuron_parameter_defaults;
+ params.init_membrane_potential = 0;
+
+ cable_cell c(morph);
+ c.paint(reg::cable(0, 0.0, 0.2), init_membrane_potential{10});
+ c.paint(reg::cable(0, 0.2, 0.7), init_membrane_potential{20});
+ c.paint(reg::cable(0, 0.7, 1.0), init_membrane_potential{30});
+
+ params.discretization = cv_policy_explicit(ls::nil());
+ fvm_cv_discretization D = fvm_cv_discretize(c, params);
+
+ ASSERT_EQ(1u, D.size());
+ EXPECT_DOUBLE_EQ(0.2*10+0.5*20+0.3*30, D.init_membrane_potential[0]);
+
+ params.discretization = cv_policy_explicit(ls::location(0, 0.3));
+ D = fvm_cv_discretize(c, params);
+
+ ASSERT_EQ(2u, D.size());
+ EXPECT_DOUBLE_EQ((0.2*10+0.1*20)/0.3, D.init_membrane_potential[0]);
+ EXPECT_DOUBLE_EQ((0.4*20+0.3*30)/0.7, D.init_membrane_potential[1]);
+}
+
+TEST(cv_layout, cable_conductance) {
+ auto morph = common_morphology::m_reg_b1; // one branch, cable constant radius.
+ const double rho = 5.; // [Ω·cm]
+
+ auto params = neuron_parameter_defaults;
+ params.axial_resistivity = rho;
+
+ cable_cell c(morph);
+ double radius = c.embedding().radius(mlocation{0, 0.5});
+ double length = c.embedding().branch_length(0);
+
+ params.discretization = cv_policy_explicit(ls::location(0, 0.3));
+ fvm_cv_discretization D = fvm_cv_discretize(c, params);
+
+ ASSERT_EQ(2u, D.size());
+
+ // Face conductance should be conductance between (relative) points 0.15 and 0.65.
+ double xa = math::pi<double>*radius*radius; // [µm^2]
+ double l = (0.65-0.15)*length; // [µm]
+ double sigma = 100 * xa/(l*rho); // [µS]
+
+ EXPECT_DOUBLE_EQ(0., D.face_conductance[0]);
+ EXPECT_DOUBLE_EQ(sigma, D.face_conductance[1]);
+}
+
+TEST(cv_layout, zero_size_cv) {
+ // Six branches; branches 0, 1 and 2 meet at (0, 1); branches
+ // 2, 3, 4, and 5 meet at (2, 1). Terminal branches are 1, 3, 4, and 5.
+ auto morph = common_morphology::m_reg_b6;
+ cable_cell cell(morph);
+
+ auto params = neuron_parameter_defaults;
+ const double rho = 5.; // [Ω·cm]
+ const double pi = math::pi<double>;
+ params.axial_resistivity = rho;
+
+ // With one CV per branch, expect reference points for face conductance
+ // to be at (0, 0.5); (0, 1); (1, 0.5); (2, 1); (3, 0.5); (4, 0.5); (5, 0.5).
+ // The first CV should be all of branch 0; the second CV should be the
+ // zero-size CV at the branch point (0, 1).
+ params.discretization = cv_policy_fixed_per_branch(1);
+ fvm_cv_discretization D = fvm_cv_discretize(cell, params);
+
+ unsigned cv_a = 0, cv_x = 1;
+ ASSERT_TRUE(util::equal(mcable_list{mcable{0, 0, 1}}, D.geometry.cables(cv_a)));
+ ASSERT_TRUE(util::equal(mcable_list{mcable{0, 1, 1}}, D.geometry.cables(cv_x)));
+
+ // Find the two CV children of CV x.
+ unsigned cv_b = -1, cv_c = -1;
+ for (unsigned i=2; i<D.size(); ++i) {
+ if ((unsigned)D.geometry.cv_parent[i]==cv_x) {
+ if (cv_b==(unsigned)-1) cv_b = i;
+ else if (cv_c==(unsigned)-1) cv_c = i;
+ else FAIL();
+ }
+ }
+
+ ASSERT_EQ(1u, D.geometry.cables(cv_b).size());
+ ASSERT_EQ(1u, D.geometry.cables(cv_c).size());
+ if (D.geometry.cables(cv_b).front().branch>D.geometry.cables(cv_c).front().branch) {
+ std::swap(cv_b, cv_c);
+ }
+
+ ASSERT_TRUE(util::equal(mcable_list{mcable{1, 0, 1}}, D.geometry.cables(cv_b)));
+ ASSERT_TRUE(util::equal(mcable_list{mcable{2, 0, 1}}, D.geometry.cables(cv_c)));
+
+ // All non-conductance values for zero-size cv_x should be zero.
+ EXPECT_EQ(0., D.cv_area[cv_x]);
+ EXPECT_EQ(0., D.cv_capacitance[cv_x]);
+ EXPECT_EQ(0., D.init_membrane_potential[cv_x]);
+ EXPECT_EQ(0., D.temperature_K[cv_x]);
+ EXPECT_EQ(0., D.diam_um[cv_x]);
+
+ // Face conductance for zero-size cv_x:
+ double l_x = cell.embedding().branch_length(0);
+ double r_x = cell.embedding().radius(mlocation{0, 0.5});
+ double sigma_x = 100 * pi * r_x * r_x / (l_x/2 * rho); // [µS]
+ EXPECT_DOUBLE_EQ(sigma_x, D.face_conductance[cv_x]);
+
+ // Face conductance for child CV cv_b:
+ double l_b = cell.embedding().branch_length(1);
+ double r_b = cell.embedding().radius(mlocation{1, 0.5});
+ double sigma_b = 100 * pi * r_b * r_b / (l_b/2 * rho); // [µS]
+ EXPECT_DOUBLE_EQ(sigma_b, D.face_conductance[cv_b]);
+
+ // Face conductance for child CV cv_c:
+ // (Distal reference point is at end of branch, so l_c not l_c/2 below.)
+ double l_c = cell.embedding().branch_length(1);
+ double r_c = cell.embedding().radius(mlocation{1, 0.5});
+ double sigma_c = 100 * pi * r_c * r_c / (l_c * rho); // [µS]
+ EXPECT_DOUBLE_EQ(sigma_c, D.face_conductance[cv_c]);
+}
diff --git a/test/unit/test_cv_policy.cpp b/test/unit/test_cv_policy.cpp
index 1f42c41f72fa7b17b654f4950405fd61789ad45e..bfb9c43369264e15307c38b798bc57d8cf028b35 100644
--- a/test/unit/test_cv_policy.cpp
+++ b/test/unit/test_cv_policy.cpp
@@ -225,3 +225,51 @@ TEST(cv_policy, max_extent) {
}
}
+TEST(cv_policy, every_sample) {
+ using namespace cv_policy_flag;
+
+ // Cell with root branch and two child branches, with multiple samples per branch.
+ // Fork is at (0., 0., 4.0).
+
+ std::vector<msample> ms;
+
+ ms.push_back({{ 0., 0., 0., 0.5}, 5});
+ for (auto i: make_span(4)) ms.push_back({{ 0., 0., i+1., 0.5}, 5});
+ for (auto i: make_span(4)) ms.push_back({{ 0., i+1., 4.0, 0.5}, 5});
+ for (auto i: make_span(4)) ms.push_back({{i+1., 0, 4.0, 0.5}, 5});
+
+ std::vector<msize_t> parents = {mnpos, 0, 1, 2, 3, 4, 5, 6, 7, 4, 9, 10, 11 };
+ morphology m{sample_tree(ms, parents), false};
+
+ // Including all samples:
+ {
+ cable_cell cell(m);
+ cv_policy pol = cv_policy_every_sample();
+ mlocation_list points = thingify(pol.cv_boundary_points(cell), cell.provider());
+ util::sort(points);
+
+ mlocation_list expected = {
+ {0, 0}, {0, 0.25}, {0, 0.5}, {0, 0.75}, {0, 1.},
+ {1, 0}, {1, 0.25}, {1, 0.5}, {1, 0.75}, {1, 1.},
+ {2, 0}, {2, 0.25}, {2, 0.5}, {2, 0.75}, {2, 1.}
+ };
+
+ EXPECT_EQ(expected, points);
+ }
+
+ // Single root CV:
+ {
+ cable_cell cell(m);
+ cv_policy pol = cv_policy_every_sample(single_root_cv);
+ mlocation_list points = thingify(pol.cv_boundary_points(cell), cell.provider());
+ util::sort(points);
+
+ mlocation_list expected = {
+ {0, 0}, {0, 1.},
+ {1, 0}, {1, 0.25}, {1, 0.5}, {1, 0.75}, {1, 1.},
+ {2, 0}, {2, 0.25}, {2, 0.5}, {2, 0.75}, {2, 1.}
+ };
+
+ EXPECT_EQ(expected, points);
+ }
+}
diff --git a/test/unit/test_event_delivery.cpp b/test/unit/test_event_delivery.cpp
index 148289e3c4c1d9cd7bda9d10f8f975bc290a8c04..dea55b57179b2c92159a9d96190f3e672af19f47 100644
--- a/test/unit/test_event_delivery.cpp
+++ b/test/unit/test_event_delivery.cpp
@@ -33,7 +33,7 @@ struct test_recipe: public n_cable_cell_recipe {
label_dict d;
d.set("soma", arb::reg::tagged(1));
- cable_cell c(st, d, true);
+ cable_cell c(st, d);
c.place(mlocation{0, 0.5}, "expsyn");
c.place(mlocation{0, 0.5}, threshold_detector{-64});
c.place(mlocation{0, 0.5}, gap_junction_site{});
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index a17d2bb1553382a3a73acb85165446379b2ba62e..2b6fda52e5be9b475dd139522ed3f0d5b186f6f3 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -1,10 +1,10 @@
#include <string>
#include <vector>
-#include <arbor/util/optional.hpp>
-#include <arbor/mechcat.hpp>
-#include <arbor/math.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/math.hpp>
+#include <arbor/mechcat.hpp>
+#include <arbor/util/optional.hpp>
#include "fvm_layout.hpp"
#include "util/maputil.hpp"
@@ -24,40 +24,6 @@ using util::count_along;
using util::value_by_key;
namespace {
- double area(const segment* s) {
- if (auto soma = s->as_soma()) {
- return math::area_sphere(soma->radius());
- }
- else if (auto cable = s->as_cable()) {
- unsigned nc = cable->num_sub_segments();
- double a = 0;
- for (unsigned i = 0; i<nc; ++i) {
- a += math::area_frustrum(cable->lengths()[i], cable->radii()[i], cable->radii()[i+1]);
- }
- return a;
- }
- else {
- return 0;
- }
- }
-
- double volume(const segment* s) {
- if (auto soma = s->as_soma()) {
- return math::volume_sphere(soma->radius());
- }
- else if (auto cable = s->as_cable()) {
- unsigned nc = cable->num_sub_segments();
- double v = 0;
- for (unsigned i = 0; i<nc; ++i) {
- v += math::volume_frustrum(cable->lengths()[i], cable->radii()[i], cable->radii()[i+1]);
- }
- return v;
- }
- else {
- return 0;
- }
- }
-
std::vector<cable_cell> two_cell_system() {
std::vector<cable_cell> cells;
@@ -117,188 +83,12 @@ namespace {
}
void check_two_cell_system(std::vector<cable_cell>& cells) {
- ASSERT_EQ(2u, cells[0].num_branches());
- ASSERT_EQ(cells[0].segment(1)->num_compartments(), 4u);
- ASSERT_EQ(cells[1].num_branches(), 4u);
- ASSERT_EQ(cells[1].segment(1)->num_compartments(), 4u);
- ASSERT_EQ(cells[1].segment(2)->num_compartments(), 4u);
- ASSERT_EQ(cells[1].segment(3)->num_compartments(), 4u);
+ ASSERT_EQ(2u, cells.size());
+ ASSERT_EQ(2u, cells[0].morphology().num_branches());
+ ASSERT_EQ(4u, cells[1].morphology().num_branches());
}
} // namespace
-TEST(fvm_layout, topology) {
- std::vector<cable_cell> cells = two_cell_system();
- check_two_cell_system(cells);
-
- fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
-
- // Expected CV layouts for cells, segment indices in paren.
- //
- // Cell 0:
- //
- // CV: | 0 ][1| 2 | 3 | 4 |5|
- // [soma (0)][ segment (1) ]
- //
- // Cell 1:
- //
- // CV: | 6 ][7| 8 | 9 | 10| 11 | 12 | 13 | 14 | 15|
- // [soma (2)][ segment (3) ][ segment (4) ]
- // [ segment (5) ]
- // | 16 | 17 | 18 | 19|
-
- EXPECT_EQ(2u, D.ncell);
- EXPECT_EQ(20u, D.ncv);
-
- unsigned nseg = 6;
- EXPECT_EQ(nseg, D.segments.size());
-
- // General sanity checks:
-
- ASSERT_EQ(D.ncell, D.cell_segment_part().size());
- ASSERT_EQ(D.ncell, D.cell_cv_part().size());
-
- ASSERT_EQ(D.ncv, D.parent_cv.size());
- ASSERT_EQ(D.ncv, D.cv_to_cell.size());
- ASSERT_EQ(D.ncv, D.face_conductance.size());
- ASSERT_EQ(D.ncv, D.cv_area.size());
- ASSERT_EQ(D.ncv, D.cv_capacitance.size());
-
- // Partitions of CVs and segments by cell:
-
- using spair = std::pair<fvm_size_type, fvm_size_type>;
- using ipair = std::pair<fvm_index_type, fvm_index_type>;
-
- EXPECT_EQ(spair(0, 2), D.cell_segment_part()[0]);
- EXPECT_EQ(spair(2, nseg), D.cell_segment_part()[1]);
-
- EXPECT_EQ(ipair(0, 6), D.cell_cv_part()[0]);
- EXPECT_EQ(ipair(6, D.ncv), D.cell_cv_part()[1]);
-
- // Segment and CV parent relationships:
-
- using ivec = std::vector<fvm_index_type>;
-
- EXPECT_EQ(ivec({0,0,1,2,3,4,6,6,7,8,9,10,11,12,13,14,11,16,17,18}), D.parent_cv);
-
- EXPECT_FALSE(D.segments[0].has_parent());
- EXPECT_EQ(1, D.segments[1].parent_cv);
-
- EXPECT_FALSE(D.segments[2].has_parent());
- EXPECT_EQ(7, D.segments[3].parent_cv);
- EXPECT_EQ(11, D.segments[4].parent_cv);
- EXPECT_EQ(11, D.segments[5].parent_cv);
-
- // Segment CV ranges (half-open, exclusing parent):
-
- EXPECT_EQ(ipair(0,1), D.segments[0].cv_range());
- EXPECT_EQ(ipair(2,6), D.segments[1].cv_range());
- EXPECT_EQ(ipair(6,7), D.segments[2].cv_range());
- EXPECT_EQ(ipair(8,12), D.segments[3].cv_range());
- EXPECT_EQ(ipair(12,16), D.segments[4].cv_range());
- EXPECT_EQ(ipair(16,20), D.segments[5].cv_range());
-
- // CV to cell index:
-
- for (auto ci: make_span(D.ncell)) {
- for (auto cv: make_span(D.cell_cv_part()[ci])) {
- EXPECT_EQ(ci, (fvm_size_type)D.cv_to_cell[cv]);
- }
- }
-}
-
-TEST(fvm_layout, diam_and_area) {
- std::vector<cable_cell> cells = two_cell_system();
- check_two_cell_system(cells);
-
- fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
-
- // Note: stick models have constant diameter segments.
- // Refer to comment above for CV vs. segment layout.
-
- EXPECT_FLOAT_EQ(12.6157, D.diam_um[0]);
- EXPECT_FLOAT_EQ(1.0 , D.diam_um[1]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[2]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[3]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[4]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[5]);
-
- EXPECT_FLOAT_EQ(14.0, D.diam_um[6]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[7]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[8]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[9]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[10]);
- EXPECT_FLOAT_EQ(1.0, D.diam_um[11]);
- EXPECT_FLOAT_EQ(0.8, D.diam_um[12]);
- EXPECT_FLOAT_EQ(0.8, D.diam_um[13]);
- EXPECT_FLOAT_EQ(0.8, D.diam_um[14]);
- EXPECT_FLOAT_EQ(0.8, D.diam_um[15]);
- EXPECT_FLOAT_EQ(0.7, D.diam_um[16]);
- EXPECT_FLOAT_EQ(0.7, D.diam_um[17]);
- EXPECT_FLOAT_EQ(0.7, D.diam_um[18]);
- EXPECT_FLOAT_EQ(0.7, D.diam_um[19]);
-
- std::vector<double> A;
- for (auto ci: make_span(D.ncell)) {
- for (auto si: make_span(cells[ci].num_branches())) {
- A.push_back(area(cells[ci].segment(si)));
- }
- }
-
- unsigned n = 4; // compartments per dendritic segment
- EXPECT_FLOAT_EQ(A[0], D.cv_area[0]);
- EXPECT_FLOAT_EQ(A[1]/(2*n), D.cv_area[1]);
- EXPECT_FLOAT_EQ(A[1]/n, D.cv_area[2]);
- EXPECT_FLOAT_EQ(A[1]/n, D.cv_area[3]);
- EXPECT_FLOAT_EQ(A[1]/n, D.cv_area[4]);
- EXPECT_FLOAT_EQ(A[1]/(2*n), D.cv_area[5]);
-
- EXPECT_FLOAT_EQ(A[2], D.cv_area[6]);
- EXPECT_FLOAT_EQ(A[3]/(2*n), D.cv_area[7]);
- EXPECT_FLOAT_EQ(A[3]/n, D.cv_area[8]);
- EXPECT_FLOAT_EQ(A[3]/n, D.cv_area[9]);
- EXPECT_FLOAT_EQ(A[3]/n, D.cv_area[10]);
- EXPECT_FLOAT_EQ((A[3]+A[4]+A[5])/(2*n), D.cv_area[11]);
- EXPECT_FLOAT_EQ(A[4]/n, D.cv_area[12]);
- EXPECT_FLOAT_EQ(A[4]/n, D.cv_area[13]);
- EXPECT_FLOAT_EQ(A[4]/n, D.cv_area[14]);
- EXPECT_FLOAT_EQ(A[4]/(2*n), D.cv_area[15]);
- EXPECT_FLOAT_EQ(A[5]/n, D.cv_area[16]);
- EXPECT_FLOAT_EQ(A[5]/n, D.cv_area[17]);
- EXPECT_FLOAT_EQ(A[5]/n, D.cv_area[18]);
- EXPECT_FLOAT_EQ(A[5]/(2*n), D.cv_area[19]);
-
- // Confirm proportional allocation of surface capacitance:
-
- // CV 9 should have area-weighted sum of the specific
- // capacitance from segments 3, 4 and 5 (cell 1 segments
- // 1, 2 and 3 respectively).
-
- double cm1 = 0.017, cm2 = 0.013, cm3 = 0.018;
-
- double c = A[3]/(2*n)*cm1+A[4]/(2*n)*cm2+A[5]/(2*n)*cm3;
- EXPECT_FLOAT_EQ(c, D.cv_capacitance[11]);
-
- double cm0 = neuron_parameter_defaults.membrane_capacitance.value();
- c = A[2]*cm0;
- EXPECT_FLOAT_EQ(c, D.cv_capacitance[6]);
-
- // Confirm face conductance within a constant diameter
- // equals a/h·1/rL where a is the cross sectional
- // area, and h is the compartment length (given the
- // regular discretization).
-
- auto cable = cells[1].segment(2)->as_cable();
- double a = volume(cable)/cable->length();
- EXPECT_FLOAT_EQ(math::pi<double>*0.8*0.8/4, a);
-
- double rL = 90;
- double h = cable->length()/4;
- double g = a/h/rL; // [µm·S/cm]
- g *= 100; // [µS]
-
- EXPECT_FLOAT_EQ(g, D.face_conductance[13]);
-}
-
TEST(fvm_layout, mech_index) {
std::vector<cable_cell> cells = two_cell_system();
check_two_cell_system(cells);
@@ -312,7 +102,7 @@ TEST(fvm_layout, mech_index) {
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
- fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
auto& hh_config = M.mechanisms.at("hh");
@@ -320,7 +110,6 @@ TEST(fvm_layout, mech_index) {
auto& exp2syn_config = M.mechanisms.at("exp2syn");
using ivec = std::vector<fvm_index_type>;
- using fvec = std::vector<fvm_value_type>;
// HH on somas of two cells, with CVs 0 and 5.
// Proportional area contrib: soma area/CV area.
@@ -328,9 +117,6 @@ TEST(fvm_layout, mech_index) {
EXPECT_EQ(mechanismKind::density, hh_config.kind);
EXPECT_EQ(ivec({0,6}), hh_config.cv);
- fvec norm_area({area(cells[0].soma())/D.cv_area[0], area(cells[1].soma())/D.cv_area[6]});
- EXPECT_TRUE(testing::seq_almost_eq<double>(norm_area, hh_config.norm_area));
-
// Three expsyn synapses, two 0.4 along segment 1, and one 0.4 along segment 5.
// These two synapses can be coalesced into 1 synapse
// 0.4 along => second (non-parent) CV for segment.
@@ -423,13 +209,12 @@ TEST(fvm_layout, coalescing_synapses) {
{
cable_cell cell = make_cell_ball_and_stick();
- // Add synapses of two varieties.
cell.place(mlocation{1, 0.3}, "expsyn");
cell.place(mlocation{1, 0.5}, "expsyn");
cell.place(mlocation{1, 0.7}, "expsyn");
cell.place(mlocation{1, 0.9}, "expsyn");
- fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -445,7 +230,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.place(mlocation{1, 0.7}, "expsyn");
cell.place(mlocation{1, 0.9}, "exp2syn");
- fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -464,7 +249,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.place(mlocation{1, 0.7}, "expsyn");
cell.place(mlocation{1, 0.9}, "expsyn");
- fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -480,7 +265,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.place(mlocation{1, 0.7}, "expsyn");
cell.place(mlocation{1, 0.9}, "exp2syn");
- fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -500,7 +285,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.place(mlocation{1, 0.7}, "expsyn");
cell.place(mlocation{1, 0.7}, "expsyn");
- fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -516,7 +301,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 0.1, 0.2));
cell.place(mlocation{1, 0.7}, syn_desc("expsyn", 0.1, 0.2));
- fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
std::vector<exp_instance> instances{
@@ -542,7 +327,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 0, 2));
cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 1, 2));
- fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
std::vector<exp_instance> instances{
@@ -571,7 +356,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.place(mlocation{1, 0.7}, syn_desc_2("exp2syn", 2, 1));
cell.place(mlocation{1, 0.7}, syn_desc_2("exp2syn", 2, 2));
- fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
for (auto &instance: {exp_instance(2, L{0,2,5}, 1, 2),
@@ -618,7 +403,7 @@ TEST(fvm_layout, synapse_targets) {
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
- fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
ASSERT_EQ(1u, M.mechanisms.count("expsyn"));
@@ -675,9 +460,9 @@ TEST(fvm_layout, density_norm_area) {
// Test area-weighted linear combination of density mechanism parameters.
// Create a cell with 4 segments:
- // - Soma (segment 0) plus three dendrites (1, 2, 3) meeting at a branch point.
+ // - Soma (branch 0) plus three dendrites (1, 2, 3) meeting at a branch point.
// - HH mechanism on all segments.
- // - Dendritic segments are given 3 compartments each.
+ // - Discretize with 3 CVs per non-soma branch, centred on forks.
//
// The CV corresponding to the branch point should comprise the terminal
// 1/6 of segment 1 and the initial 1/6 of segments 2 and 3.
@@ -736,14 +521,12 @@ TEST(fvm_layout, density_norm_area) {
std::vector<double> expected_gkbar(ncv, dflt_gkbar);
std::vector<double> expected_gl(ncv, dflt_gl);
- auto div_by_ends = [](const cable_segment* cable) {
- return div_compartment_by_ends(cable->num_compartments(), cable->radii(), cable->lengths());
- };
- auto& segs = cells[0].segments();
- double soma_area = area(segs[0].get());
- auto seg1_divs = div_by_ends(segs[1]->as_cable());
- auto seg2_divs = div_by_ends(segs[2]->as_cable());
- auto seg3_divs = div_by_ends(segs[3]->as_cable());
+ // Last 1/6 of segment 1
+ double seg1_area_right = cells[0].embedding().integrate_area(mcable{1, 5./6., 1.});
+ // First 1/6 of segment 2
+ double seg2_area_left = cells[0].embedding().integrate_area(mcable{2, 0., 1./6.});
+ // First 1/6 of segment 3
+ double seg3_area_left = cells[0].embedding().integrate_area(mcable{3, 0., 1./6.});
// CV 0: soma
// CV1: left of segment 1
@@ -754,8 +537,8 @@ TEST(fvm_layout, density_norm_area) {
expected_gl[3] = seg1_gl;
// CV 4: mix of right of segment 1 and left of segments 2 and 3.
- expected_gkbar[4] = wmean(seg1_divs(2).right.area, dflt_gkbar, seg2_divs(0).left.area, seg2_gkbar, seg3_divs(0).left.area, seg3_gkbar);
- expected_gl[4] = wmean(seg1_divs(2).right.area, seg1_gl, seg2_divs(0).left.area, dflt_gl, seg3_divs(0).left.area, seg3_gl);
+ expected_gkbar[4] = wmean(seg1_area_right, dflt_gkbar, seg2_area_left, seg2_gkbar, seg3_area_left, seg3_gkbar);
+ expected_gl[4] = wmean(seg1_area_right, seg1_gl, seg2_area_left, dflt_gl, seg3_area_left, seg3_gl);
// CV 5-7: just segment 2
expected_gkbar[5] = seg2_gkbar;
@@ -773,26 +556,9 @@ TEST(fvm_layout, density_norm_area) {
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
- fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
- // Check CV area assumptions.
- // Note: area integrator used here and in `fvm_multicell` may differ, and so areas computed may
- // differ some due to rounding area, even given that we're dealing with simple truncated cones
- // for segments. Check relative error within a tolerance of (say) 10 epsilon.
-
- double area_relerr = 10*std::numeric_limits<double>::epsilon();
- EXPECT_TRUE(testing::near_relative(D.cv_area[0],
- soma_area, area_relerr));
- EXPECT_TRUE(testing::near_relative(D.cv_area[1],
- seg1_divs(0).left.area, area_relerr));
- EXPECT_TRUE(testing::near_relative(D.cv_area[2],
- seg1_divs(0).right.area+seg1_divs(1).left.area, area_relerr));
- EXPECT_TRUE(testing::near_relative(D.cv_area[4],
- seg1_divs(2).right.area+seg2_divs(0).left.area+seg3_divs(0).left.area, area_relerr));
- EXPECT_TRUE(testing::near_relative(D.cv_area[7],
- seg2_divs(2).right.area, area_relerr));
-
// Grab the HH parameters from the mechanism.
EXPECT_EQ(1u, M.mechanisms.size());
@@ -814,7 +580,7 @@ TEST(fvm_layout, valence_verify) {
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
- fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, neuron_parameter_defaults);
mechanism_catalogue testcat = make_unit_test_catalogue();
gprop.catalogue = &testcat;
@@ -903,7 +669,7 @@ TEST(fvm_layout, ion_weights) {
std::vector<cable_cell> cells{std::move(c)};
- fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
ASSERT_EQ(1u, M.ions.count("ca"s));
@@ -957,7 +723,7 @@ TEST(fvm_layout, revpot) {
auto test_gprop = gprop;
test_gprop.default_parameters.reversal_potential_method["b"] = write_eb_ec;
- fvm_discretization D = fvm_discretize(cells, test_gprop.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, test_gprop.default_parameters);
EXPECT_THROW(fvm_build_mechanism_data(test_gprop, cells, D), cable_cell_error);
}
@@ -968,7 +734,7 @@ TEST(fvm_layout, revpot) {
test_gprop.default_parameters.reversal_potential_method["c"] = write_eb_ec;
cells[1].default_parameters.reversal_potential_method["c"] = "write_eX/c";
- fvm_discretization D = fvm_discretize(cells, test_gprop.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, test_gprop.default_parameters);
EXPECT_THROW(fvm_build_mechanism_data(test_gprop, cells, D), cable_cell_error);
}
@@ -977,12 +743,12 @@ TEST(fvm_layout, revpot) {
cell1_prop.reversal_potential_method["b"] = write_eb_ec;
cell1_prop.reversal_potential_method["c"] = write_eb_ec;
- fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
// Only CV which needs write_multiple_eX/x=b,y=c is the soma (first CV)
// of the second cell.
- auto soma1_index = D.cell_cv_bounds[1];
+ auto soma1_index = D.geometry.cell_cv_divs[1];
ASSERT_EQ(1u, M.mechanisms.count(write_eb_ec.name()));
EXPECT_EQ((std::vector<fvm_index_type>(1, soma1_index)), M.mechanisms.at(write_eb_ec.name()).cv);
}
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 434e74245fc896418e9ea673ae5fb4f810493d46..5f7b116ebe0fa9570832583d0ed7a698e4007778 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -9,7 +9,6 @@
#include <arbor/load_balance.hpp>
#include <arbor/math.hpp>
#include <arbor/cable_cell.hpp>
-#include <arbor/segment.hpp>
#include <arbor/recipe.hpp>
#include <arbor/sampling.hpp>
#include <arbor/simulation.hpp>
@@ -26,6 +25,7 @@
#include "util/meta.hpp"
#include "util/maputil.hpp"
#include "util/rangeutil.hpp"
+#include "util/transform.hpp"
#include "common.hpp"
#include "unit_test_catalogue.hpp"
@@ -249,8 +249,8 @@ TEST(fvm_lowered, target_handles) {
make_cell_ball_and_3stick()
};
- EXPECT_EQ(cells[0].num_branches(), 2u);
- EXPECT_EQ(cells[1].num_branches(), 4u);
+ EXPECT_EQ(cells[0].morphology().num_branches(), 2u);
+ EXPECT_EQ(cells[1].morphology().num_branches(), 4u);
// (in increasing target order)
cells[0].place(mlocation{1, 0.4}, "expsyn");
@@ -334,7 +334,7 @@ TEST(fvm_lowered, stimulus) {
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
- fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
const auto& A = D.cv_area;
std::vector<target_handle> targets;
@@ -702,8 +702,8 @@ TEST(fvm_lowered, point_ionic_current) {
// Test area-weighted linear combination of ion species concentrations
TEST(fvm_lowered, weighted_write_ion) {
- // Create a cell with 4 segments (same morphopology as in fvm_layout.ion_weights test):
- // - Soma (segment 0) plus three dendrites (1, 2, 3) meeting at a branch point.
+ // Create a cell with 4 branches (same morphopology as in fvm_layout.ion_weights test):
+ // - Soma (branch 0) plus three dendrites (1, 2, 3) meeting at a branch point.
// - Dendritic segments are given 1 compartments each.
//
// /
@@ -714,7 +714,7 @@ TEST(fvm_lowered, weighted_write_ion) {
// d3
//
// The CV corresponding to the branch point should comprise the terminal
- // 1/2 of segment 1 and the initial 1/2 of segments 2 and 3.
+ // 1/2 of branch 1 and the initial 1/2 of branches 2 and 3.
//
// Geometry:
// soma 0: radius 5 µm
@@ -738,11 +738,11 @@ TEST(fvm_lowered, weighted_write_ion) {
const double con_int = 80;
const double con_ext = 120;
- // Ca ion reader test_kinlva on CV 2 and 3 via segment 2:
- c.segments()[2] ->add_mechanism("test_kinlva");
+ // Ca ion reader test_kinlva on CV 2 and 3 via branch 2:
+ c.paint(reg::branch(2), "test_kinlva");
- // Ca ion writer test_ca on CV 2 and 4 via segment 3:
- c.segments()[3] ->add_mechanism("test_ca");
+ // Ca ion writer test_ca on CV 2 and 4 via branch 3:
+ c.paint(reg::branch(3), "test_ca");
cable1d_recipe rec(c);
rec.add_ion("ca", 2, con_int, con_ext, 0.0);
@@ -764,7 +764,7 @@ TEST(fvm_lowered, weighted_write_ion) {
std::vector<double> ion_init_iconc = util::assign_from(ion.init_Xi_);
std::vector<double> expected_init_iconc = {0.75*con_int, 1.*con_int, 0};
- EXPECT_EQ(expected_init_iconc, ion_init_iconc);
+ EXPECT_TRUE(testing::seq_almost_eq<double>(expected_init_iconc, ion_init_iconc));
auto test_ca = dynamic_cast<multicore::mechanism*>(find_mechanism(fvcell, "test_ca"));
@@ -791,7 +791,7 @@ TEST(fvm_lowered, weighted_write_ion) {
ion.init_concentration();
test_ca->write_ions();
std::vector<double> ion_iconc = util::assign_from(ion.Xi_);
- EXPECT_EQ(expected_iconc, ion_iconc);
+ EXPECT_TRUE(testing::seq_almost_eq<double>(expected_iconc, ion_iconc));
}
TEST(fvm_lowered, gj_coords_simple) {
@@ -837,7 +837,7 @@ TEST(fvm_lowered, gj_coords_simple) {
cells.push_back(std::move(c));
}
- fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, neuron_parameter_defaults);
std::vector<cell_gid_type> gids = {0, 1};
auto GJ = fvcell.fvm_gap_junctions(cells, gids, rec, D);
@@ -854,8 +854,6 @@ TEST(fvm_lowered, gj_coords_simple) {
}
TEST(fvm_lowered, gj_coords_complex) {
- using pair = std::pair<int, int>;
-
class gap_recipe: public recipe {
public:
gap_recipe() {}
@@ -868,23 +866,26 @@ TEST(fvm_lowered, gj_coords_complex) {
std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
std::vector<gap_junction_connection> conns;
switch (gid) {
- case 0 : return {
- gap_junction_connection({2, 0}, {0, 1}, 0.01),
- gap_junction_connection({1, 0}, {0, 0}, 0.03),
- gap_junction_connection({1, 1}, {0, 0}, 0.04)
- };
- case 1 : return {
- gap_junction_connection({0, 0}, {1, 0}, 0.03),
- gap_junction_connection({0, 0}, {1, 1}, 0.04),
- gap_junction_connection({2, 1}, {1, 2}, 0.02),
- gap_junction_connection({2, 2}, {1, 3}, 0.01)
- };
- case 2 : return {
- gap_junction_connection({0, 1}, {2, 0}, 0.01),
- gap_junction_connection({1, 2}, {2, 1}, 0.02),
- gap_junction_connection({1, 3}, {2, 2}, 0.01)
- };
- default : return {};
+ case 0:
+ return {
+ gap_junction_connection({2, 0}, {0, 1}, 0.01),
+ gap_junction_connection({1, 0}, {0, 0}, 0.03),
+ gap_junction_connection({1, 1}, {0, 0}, 0.04)
+ };
+ case 1:
+ return {
+ gap_junction_connection({0, 0}, {1, 0}, 0.03),
+ gap_junction_connection({0, 0}, {1, 1}, 0.04),
+ gap_junction_connection({2, 1}, {1, 2}, 0.02),
+ gap_junction_connection({2, 2}, {1, 3}, 0.01)
+ };
+ case 2:
+ return {
+ gap_junction_connection({0, 1}, {2, 0}, 0.01),
+ gap_junction_connection({1, 2}, {2, 1}, 0.02),
+ gap_junction_connection({1, 3}, {2, 2}, 0.01)
+ };
+ default : return {};
}
return conns;
}
@@ -901,8 +902,10 @@ TEST(fvm_lowered, gj_coords_complex) {
b0.add_branch(0, 8, 0.3, 0.2, 4, "dend");
auto c0 = b0.make_cell();
- c0.place(mlocation{1, 1}, gap_junction_site{});
- c0.place(mlocation{1, 0.5}, gap_junction_site{});
+ mlocation c0_gj[2] = {{1, 1}, {1, 0.5}};
+
+ c0.place(c0_gj[0], gap_junction_site{});
+ c0.place(c0_gj[1], gap_junction_site{});
soma_cell_builder b1(1.4);
b1.add_branch(0, 12, 0.3, 0.5, 6, "dend");
@@ -910,10 +913,13 @@ TEST(fvm_lowered, gj_coords_complex) {
b1.add_branch(1, 5, 0.2, 0.2, 5, "dend");
auto c1 = b1.make_cell();
- c1.place(mlocation{2, 1}, gap_junction_site{});
- c1.place(mlocation{1, 1}, gap_junction_site{});
- c1.place(mlocation{1, 0.45}, gap_junction_site{});
- c1.place(mlocation{1, 0.1}, gap_junction_site{});
+ mlocation c1_gj[4] = {{2, 1}, {1, 1}, {1, 0.45}, {1, 0.1}};
+
+ c1.place(c1_gj[0], gap_junction_site{});
+ c1.place(c1_gj[1], gap_junction_site{});
+ c1.place(c1_gj[2], gap_junction_site{});
+ c1.place(c1_gj[3], gap_junction_site{});
+
soma_cell_builder b2(2.9);
b2.add_branch(0, 4, 0.3, 0.5, 2, "dend");
@@ -923,9 +929,11 @@ TEST(fvm_lowered, gj_coords_complex) {
b2.add_branch(2, 4, 0.2, 0.2, 2, "dend");
auto c2 = b2.make_cell();
- c2.place(mlocation{1, 0.5}, gap_junction_site{});
- c2.place(mlocation{4, 1}, gap_junction_site{});
- c2.place(mlocation{2, 1}, gap_junction_site{});
+ mlocation c2_gj[3] = {{1, 0.5}, {4, 1}, {2, 1}};
+
+ c2.place(c2_gj[0], gap_junction_site{});
+ c2.place(c2_gj[1], gap_junction_site{});
+ c2.place(c2_gj[2], gap_junction_site{});
std::vector<cable_cell> cells{std::move(c0), std::move(c1), std::move(c2)};
@@ -935,32 +943,48 @@ TEST(fvm_lowered, gj_coords_complex) {
gap_recipe rec;
fvcell.fvm_intdom(rec, gids, cell_to_intdom);
- fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
+ fvm_cv_discretization D = fvm_cv_discretize(cells, neuron_parameter_defaults);
- auto GJ = fvcell.fvm_gap_junctions(cells, gids, rec, D);
+ int c0_gj_cv[2];
+ for (int i = 0; i<2; ++i) c0_gj_cv[i] = D.geometry.location_cv(0, c0_gj[i]);
+
+ int c1_gj_cv[4];
+ for (int i = 0; i<4; ++i) c1_gj_cv[i] = D.geometry.location_cv(1, c1_gj[i]);
+
+ int c2_gj_cv[3];
+ for (int i = 0; i<3; ++i) c2_gj_cv[i] = D.geometry.location_cv(2, c2_gj[i]);
+
+ std::vector<fvm_gap_junction> GJ = fvcell.fvm_gap_junctions(cells, gids, rec, D);
EXPECT_EQ(10u, GJ.size());
auto weight = [&](fvm_value_type g, fvm_index_type i){
return g * 1e3 / D.cv_area[i];
};
- std::vector<pair> expected_loc = {{5, 16}, {5,13}, {3,24}, {16, 5}, {13,5} ,{10,31}, {8, 27}, {24,3}, {31,10}, {27, 8}};
- std::vector<double> expected_weight = {
- weight(0.03, 5), weight(0.04, 5), weight(0.01, 3), weight(0.03, 16), weight(0.04, 13),
- weight(0.02, 10), weight(0.01, 8), weight(0.01, 24), weight(0.02, 31), weight(0.01, 27)
+ std::vector<fvm_gap_junction> expected = {
+ {{c0_gj_cv[0], c1_gj_cv[0]}, weight(0.03, c0_gj_cv[0])},
+ {{c0_gj_cv[0], c1_gj_cv[1]}, weight(0.04, c0_gj_cv[0])},
+ {{c0_gj_cv[1], c2_gj_cv[0]}, weight(0.01, c0_gj_cv[1])},
+ {{c1_gj_cv[0], c0_gj_cv[0]}, weight(0.03, c1_gj_cv[0])},
+ {{c1_gj_cv[1], c0_gj_cv[0]}, weight(0.04, c1_gj_cv[1])},
+ {{c1_gj_cv[2], c2_gj_cv[1]}, weight(0.02, c1_gj_cv[2])},
+ {{c1_gj_cv[3], c2_gj_cv[2]}, weight(0.01, c1_gj_cv[3])},
+ {{c2_gj_cv[0], c0_gj_cv[1]}, weight(0.01, c2_gj_cv[0])},
+ {{c2_gj_cv[1], c1_gj_cv[2]}, weight(0.02, c2_gj_cv[1])},
+ {{c2_gj_cv[2], c1_gj_cv[3]}, weight(0.01, c2_gj_cv[2])}
};
- for (unsigned i = 0; i < GJ.size(); i++) {
- bool found = false;
- for (unsigned j = 0; j < expected_loc.size(); j++) {
- if (expected_loc[j].first == GJ[i].loc.first && expected_loc[j].second == GJ[i].loc.second) {
- found = true;
- EXPECT_EQ(expected_weight[j], GJ[i].weight);
- break;
- }
- }
- EXPECT_TRUE(found);
- }
+ using util::sort_by;
+ using util::transform_view;
+
+ auto gj_loc = [](const fvm_gap_junction gj) { return gj.loc; };
+ auto gj_weight = [](const fvm_gap_junction gj) { return gj.weight; };
+
+ sort_by(GJ, [](fvm_gap_junction gj) { return gj.loc; });
+ sort_by(expected, [](fvm_gap_junction gj) { return gj.loc; });
+
+ EXPECT_TRUE(testing::seq_eq(transform_view(expected, gj_loc), transform_view(GJ, gj_loc)));
+ EXPECT_TRUE(testing::seq_almost_eq<double>(transform_view(expected, gj_weight), transform_view(GJ, gj_weight)));
}
TEST(fvm_lowered, cell_group_gj) {
@@ -1019,8 +1043,8 @@ TEST(fvm_lowered, cell_group_gj) {
auto num_dom0 = fvcell.fvm_intdom(rec, gids_cg0, cell_to_intdom0);
auto num_dom1 = fvcell.fvm_intdom(rec, gids_cg1, cell_to_intdom1);
- fvm_discretization D0 = fvm_discretize(cell_group0, neuron_parameter_defaults);
- fvm_discretization D1 = fvm_discretize(cell_group1, neuron_parameter_defaults);
+ fvm_cv_discretization D0 = fvm_cv_discretize(cell_group0, neuron_parameter_defaults);
+ fvm_cv_discretization D1 = fvm_cv_discretize(cell_group1, neuron_parameter_defaults);
auto GJ0 = fvcell.fvm_gap_junctions(cell_group0, gids_cg0, rec, D0);
auto GJ1 = fvcell.fvm_gap_junctions(cell_group1, gids_cg1, rec, D1);
diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp
index d4d51b7e572ada70f2dbf3d73a54d84ed2491c2c..cfe20b754594932d8e6a992185301c3c59a21f42 100644
--- a/test/unit/test_morphology.cpp
+++ b/test/unit/test_morphology.cpp
@@ -578,10 +578,6 @@ TEST(morphology, swc) {
// Test that the morphology contains the expected number of branches.
auto m = arb::morphology(sm);
EXPECT_EQ(31u, m.num_branches());
-
- // Confirm that converting to a cable_cell generates the same number of branches.
- auto c = arb::cable_cell(m, {}, false);
- EXPECT_EQ(31u, c.num_branches());
}
TEST(morphology, minset) {
diff --git a/test/unit/test_range.cpp b/test/unit/test_range.cpp
index 39a80f45ba027f057fffb3c09dc4a1d06797bd31..ed48a7923c35d618bef62f6e4cee2a3c99e0180b 100644
--- a/test/unit/test_range.cpp
+++ b/test/unit/test_range.cpp
@@ -1,6 +1,7 @@
#include "../gtest.h"
#include <algorithm>
+#include <array>
#include <iterator>
#include <functional>
#include <list>
@@ -668,6 +669,63 @@ TEST(range, is_sorted_by) {
EXPECT_TRUE(util::is_sorted_by(seq, [](int x) { return x+2; }, std::greater<int>{}));
}
+template <typename V>
+struct repeat_iterator {
+ typedef std::input_iterator_tag iterator_category;
+ typedef const V& reference;
+ typedef const V* pointer;
+ typedef V value_type;
+ typedef std::ptrdiff_t difference_type;
+
+ V v;
+ repeat_iterator(V v): v(std::move(v)) {}
+
+ bool operator==(const repeat_iterator<V>& i) const { return true; }
+ bool operator!=(const repeat_iterator<V>& i) const { return false; }
+ reference operator*() const { return v; }
+ repeat_iterator& operator++() { return *this; }
+ repeat_iterator& operator++(int) { return *this; }
+ pointer operator->() const { return &v; }
+};
+
+struct never_t {
+ template <typename X> friend bool operator==(never_t, const X&) { return false; }
+ template <typename X> friend bool operator==(const X&, never_t) { return false; }
+
+ template <typename X> friend bool operator!=(never_t, const X&) { return true; }
+ template <typename X> friend bool operator!=(const X&, never_t) { return true; }
+};
+static never_t never;
+
+template <typename V>
+auto repeat(V v) { return util::make_range(repeat_iterator<V>(std::move(v)), never); }
+
+TEST(range, equal) {
+ // Finite containers
+ unsigned a[5] = { 1, 3, 2, 5, 4};
+ std::array<unsigned, 5> b = { 1, 3, 2, 5, 4};
+
+ EXPECT_TRUE(util::equal(a, b));
+
+ a[3] = 10;
+ EXPECT_FALSE(util::equal(a, b));
+
+ unsigned abis[6] = { 1, 3, 2, 5, 4, 6};
+ EXPECT_FALSE(util::equal(abis, b));
+
+ std::vector<std::string> empty1;
+ std::vector<std::string> empty2;
+ EXPECT_TRUE(util::equal(empty1, empty2));
+
+ empty2.push_back("hello");
+ EXPECT_FALSE(util::equal(empty1, empty2));
+
+ // Infinite sequence
+ unsigned c[3] = { 2, 2, 2 };
+ EXPECT_FALSE(util::equal(c, repeat(2u)));
+ EXPECT_FALSE(util::equal(repeat(5u), repeat(2u)));
+}
+
TEST(range, reverse) {
// make a C string into a sentinel-terminated range
auto cstr = [](const char* s) { return util::make_range(s, null_terminated); };
diff --git a/test/unit/test_segment.cpp b/test/unit/test_segment.cpp
deleted file mode 100644
index 509f889c3ac381a0f7ef5ced0fa984ca294b86c9..0000000000000000000000000000000000000000
--- a/test/unit/test_segment.cpp
+++ /dev/null
@@ -1,69 +0,0 @@
-#include <vector>
-
-#include "../gtest.h"
-
-#include <arbor/math.hpp>
-#include <arbor/segment.hpp>
-
-using namespace arb;
-
-TEST(segment, kinfs) {
- using ::arb::math::pi;
-
- {
- auto s = make_segment<soma_segment>(1.0);
- EXPECT_EQ(s->kind(), section_kind::soma);
- }
-
- {
- auto s = make_segment<soma_segment>(1.0, point<double>(0., 1., 2.));
- EXPECT_EQ(s->kind(), section_kind::soma);
- }
-
- double length = 1./pi<double>;
- double radius = 1.;
-
- // single cylindrical frustrum
- {
- auto s = make_segment<cable_segment>(section_kind::dendrite, radius, radius, length);
- EXPECT_EQ(s->kind(), section_kind::dendrite);
- }
-
- // cable made up of three identical cylindrical frustrums
- {
- auto s =
- make_segment<cable_segment>(
- section_kind::axon,
- std::vector<double>{radius, radius, radius, radius},
- std::vector<double>{length, length, length}
- );
-
- EXPECT_EQ(s->kind(), section_kind::axon);
- }
-
- // single frustrum of length 1 and radii 1 and 2
- // the centre of each end are at the origin (0,0,0) and (0,1,0)
- {
- auto s =
- make_segment<cable_segment>(
- section_kind::dendrite,
- 1, 2,
- point<double>(0,0,0), point<double>(0,1,0)
- );
-
- EXPECT_EQ(s->kind(), section_kind::dendrite);
- }
-
- // cable made up of three frustrums
- // that emulate the frustrum from the previous single-frustrum case
- {
- auto s =
- make_segment<cable_segment>(
- section_kind::axon,
- std::vector<double>{1, 1.5, 2},
- std::vector<point<double>>{ {0,0,0}, {0,0.5,0}, {0,1,0} }
- );
-
- EXPECT_EQ(s->kind(), section_kind::axon);
- }
-}
diff --git a/test/unit/test_swcio.cpp b/test/unit/test_swcio.cpp
index c12fd36cc4c4e61c2e04438c354e3d1128867083..3a004c5e155e85d3f85786f244087d9664409672 100644
--- a/test/unit/test_swcio.cpp
+++ b/test/unit/test_swcio.cpp
@@ -406,144 +406,3 @@ TEST(swc_parser, raw)
}
}
-/***
- * If you have proper unit tests for the steps from
- * {sample_tree -> morphology -> cable_cell}
- * then these tests can be simplified to check that {swc -> sample_tree} works
-TEST(swc_io, cell_construction) {
- using namespace arb;
-
- //
- // 0
- // |
- // 1
- // |
- // 2
- // / \.
- // 3 4
- // \.
- // 5
- //
-
- std::stringstream is;
- is << "1 1 0 0 0 2.1 -1\n";
- is << "2 3 0.1 1.2 1.2 1.3 1\n";
- is << "3 3 1.0 2.0 2.2 1.1 2\n";
- is << "4 3 1.5 3.3 1.3 2.2 3\n";
- is << "5 3 2.5 5.3 2.5 0.7 3\n";
- is << "6 3 3.5 2.3 3.7 3.4 5\n";
-
- using point_type = point<double>;
- std::vector<point_type> points = {
- { 0.0, 0.0, 0.0 },
- { 0.1, 1.2, 1.2 },
- { 1.0, 2.0, 2.2 },
- { 1.5, 3.3, 1.3 },
- { 2.5, 5.3, 2.5 },
- { 3.5, 2.3, 3.7 },
- };
-
- // swc -> morphology
- auto morph = swc_as_morphology(parse_swc_file(is));
-
- cable_cell cell = make_cable_cell(morph, true);
- EXPECT_TRUE(cell.has_soma());
- EXPECT_EQ(4u, cell.num_segments());
-
- EXPECT_DOUBLE_EQ(norm(points[1]-points[2]), cell.cable(1)->length());
- EXPECT_DOUBLE_EQ(norm(points[2]-points[3]), cell.cable(2)->length());
- EXPECT_DOUBLE_EQ(norm(points[2]-points[4]) + norm(points[4]-points[5]),
- cell.cable(3)->length());
-
-
- // Check each segment separately
- EXPECT_TRUE(cell.segment(0)->is_soma());
- EXPECT_EQ(2.1, cell.soma()->radius());
- EXPECT_EQ(point_type(0, 0, 0), cell.soma()->center());
-
- for (auto i = 1u; i < cell.num_segments(); ++i) {
- EXPECT_TRUE(cell.segment(i)->is_dendrite());
- }
-
- EXPECT_EQ(1u, cell.cable(1)->num_sub_segments());
- EXPECT_EQ(1u, cell.cable(2)->num_sub_segments());
- EXPECT_EQ(2u, cell.cable(3)->num_sub_segments());
-
- // We asked to use the same discretization as in the SWC, so check number of compartments too.
- EXPECT_EQ(1u, cell.cable(1)->num_compartments());
- EXPECT_EQ(1u, cell.cable(2)->num_compartments());
- EXPECT_EQ(2u, cell.cable(3)->num_compartments());
-
- // Check the radii
- EXPECT_EQ(1.3, cell.cable(1)->radius(0));
- EXPECT_EQ(1.1, cell.cable(1)->radius(1));
-
- EXPECT_EQ(1.1, cell.cable(2)->radius(0));
- EXPECT_EQ(2.2, cell.cable(2)->radius(1));
-
- EXPECT_EQ(1.1, cell.cable(3)->radius(0));
- EXPECT_EQ(3.4, cell.cable(3)->radius(1));
-
- auto len_ratio = norm(points[2]-points[4]) / cell.cable(3)->length();
- EXPECT_NEAR(.7, cell.cable(3)->radius(len_ratio), 1e-15);
-
- // Double-check radii at joins are equal
- EXPECT_EQ(cell.cable(1)->radius(1),
- cell.cable(2)->radius(0));
-
- EXPECT_EQ(cell.cable(1)->radius(1),
- cell.cable(3)->radius(0));
-}
-
-void expect_morph_eq(const morphology& a, const morphology& b) {
- EXPECT_EQ(a.soma.r, b.soma.r);
- EXPECT_EQ(a.sections.size(), b.sections.size());
-
- for (unsigned i = 0; i<a.sections.size(); ++i) {
- const section_geometry& r = a.sections[i];
- const section_geometry& s = b.sections[i];
-
- EXPECT_EQ(r.points.size(), s.points.size());
- unsigned n = r.points.size();
- if (n>1) {
- auto r0 = r.points[0];
- auto s0 = s.points[0];
- EXPECT_EQ(r0.r, s0.r);
-
- // TODO: check lengths for each line segment instead.
- EXPECT_DOUBLE_EQ(r.length, s.length);
-
- for (unsigned i = 1; i<n; ++i) {
- auto r1 = r.points[i];
- auto s1 = s.points[i];
-
- EXPECT_EQ(r1.r, s1.r);
- }
- }
- }
-}
-
-// check that simple ball and stick model with one dendrite attached to a soma
-// which is used in the validation tests can be loaded from file and matches
-// the one generated with the C++ interface
-TEST(swc_parser, from_file_ball_and_stick) {
- std::string datadir{DATADIR};
- auto fname = datadir + "/ball_and_stick.swc";
- std::ifstream fid(fname);
- if (!fid.is_open()) {
- std::cerr << "unable to open file " << fname << "... skipping test\n";
- return;
- }
-
- // read the file as morhpology
- auto bas_morph = swc_as_morphology(parse_swc_file(fid));
-
- // compare with expected morphology
- morphology expected;
-
- expected.soma = {0., 0., 0., 6.30785};
- expected.add_section({{6.30785, 0., 0., 0.5}, {206.30785, 0., 0., 0.5}}, 0u, section_kind::dendrite);
-
- expect_morph_eq(expected, bas_morph);
-}
-*/