From b2e5d1f0f1cd181a34dc62a860589091f52ed760 Mon Sep 17 00:00:00 2001
From: Sam Yates <sam@quux.dropbear.id.au>
Date: Fri, 22 May 2020 15:06:10 +0200
Subject: [PATCH] Extend cable cell probe variety. (#1034)
Fixes #589 and #730, providing cell-wide probes and a correct total trans-membrane current.
Fixes #822, providing spatial interpolation of membrane voltages as governed by cable resistivity and diameter.
* Add partition of point mechanism target vector by cell index to `fvm_mechanism_config`.
* Fix sign error in axial current interpolation.
* Replace `fvm_probe_info` with class that wraps a variant type capturing one of a set of backend-cellgroup probe translation representations: `fvm_probe_scalar`, `fvm_probe_interpolated`, `fvm_probe_multi`, `fvm_probe_weighted_multi`, `fvm_probe_membrance_currents`.
* Refactor `fvm_lowered_cell_impl::resolve_probe_address` so that it avoids the long chain of if-else tests on wrapped `any` type with an invocation of `util::any_visitor`. Each specific cable cell probe address representation is translated into one of the intermediate `fvm_probe_info` variants via an overload of `resolve_probe`.
* Add new non-scalar cell probe address types for all-of-cell values, including point mechanism state, density mechanism state, ion currents and concentrations, membrane voltage by CV cable, and total trans-mebrane currents. The latter are ultimately calculated by computing the net axial current fluxes in each CV.
* Add `util::any_ptr` field to the sampling function interface, pointing to constant-per-simulation-run metadata associated with a probe address. Underlying metadata type varies between probe types.
* Extend `simple_sampler` so that it optionally captures sampler metadata.
* Add special case to `simple_sampler` for vector-valued sample data, allowing it to translate from a pointer-pair range to a std::vector.
* Fix minor bugs in `util::variant`.
* Add to `mc_cell_group` the logic for translating each intermediate probe representation to a sampler callback invocation, reusing temporary scratch vectors across calls.
* Add unit comments to `embed_pwlin_data`.
* Add probe variety demonstration `probe-demo` to examples, covering most of the new cable cell probe classes.
* Add `util::any_visitor` runtime-type dependent dispatcher class.
* Add `util::overload` helper for generating overloaded functional objects from individual functions, lambdas or functionals.
* Extend probe unit tests in `test_probe.cpp` to cover new probe address types.
* Add unit tests for `any_visitor` and `overload`.
---
arbor/backends/gpu/fvm.hpp | 1 +
arbor/backends/multicore/fvm.hpp | 1 +
arbor/fvm_layout.cpp | 46 +-
arbor/fvm_layout.hpp | 6 +-
arbor/fvm_lowered_cell.hpp | 105 ++-
arbor/fvm_lowered_cell_impl.hpp | 443 ++++++++---
arbor/include/arbor/cable_cell.hpp | 129 ++-
arbor/include/arbor/sampling.hpp | 8 +-
arbor/include/arbor/simple_sampler.hpp | 85 +-
arbor/include/arbor/util/any_visitor.hpp | 155 ++++
arbor/include/arbor/util/variant.hpp | 7 +-
arbor/mc_cell_group.cpp | 275 ++++++-
arbor/mc_cell_group.hpp | 2 +-
arbor/morph/embed_pwlin.cpp | 10 +-
doc/cpp_cable_cell.rst | 303 +++++++
doc/sampling_api.rst | 13 +-
example/CMakeLists.txt | 1 +
example/dryrun/dryrun.cpp | 6 +-
example/gap_junctions/gap_junctions.cpp | 2 +-
example/generators/generators.cpp | 15 +-
example/generators/readme.md | 12 +-
example/probe-demo/CMakeLists.txt | 3 +
example/probe-demo/probe-demo.cpp | 247 ++++++
example/ring/ring.cpp | 2 +-
example/single/single.cpp | 2 +-
python/recipe.cpp | 4 +-
python/sampling.cpp | 2 +-
python/single_cell_model.cpp | 4 +-
test/unit/CMakeLists.txt | 3 +
test/unit/common_morphologies.hpp | 5 +-
test/unit/mod/ca_linear.mod | 19 +
test/unit/mod/param_as_state.mod | 24 +
test/unit/test_any_visitor.cpp | 190 +++++
test/unit/test_cv_layout.cpp | 6 +-
test/unit/test_fvm_layout.cpp | 12 +-
test/unit/test_fvm_lowered.cpp | 31 +-
test/unit/test_probe.cpp | 958 +++++++++++++++++++----
test/unit/test_schedule.cpp | 2 +-
test/unit/unit_test_catalogue.cpp | 8 +-
test/unit/unit_test_catalogue.hpp | 2 +-
40 files changed, 2776 insertions(+), 373 deletions(-)
create mode 100644 arbor/include/arbor/util/any_visitor.hpp
create mode 100644 example/probe-demo/CMakeLists.txt
create mode 100644 example/probe-demo/probe-demo.cpp
create mode 100644 test/unit/mod/ca_linear.mod
create mode 100644 test/unit/mod/param_as_state.mod
create mode 100644 test/unit/test_any_visitor.cpp
diff --git a/arbor/backends/gpu/fvm.hpp b/arbor/backends/gpu/fvm.hpp
index 0a56bf3c..0fa52b63 100644
--- a/arbor/backends/gpu/fvm.hpp
+++ b/arbor/backends/gpu/fvm.hpp
@@ -47,6 +47,7 @@ struct backend {
using sample_event_stream = arb::gpu::sample_event_stream;
using shared_state = arb::gpu::shared_state;
+ using ion_state = arb::gpu::ion_state;
static threshold_watcher voltage_watcher(
const shared_state& state,
diff --git a/arbor/backends/multicore/fvm.hpp b/arbor/backends/multicore/fvm.hpp
index 46767b34..d7b233de 100644
--- a/arbor/backends/multicore/fvm.hpp
+++ b/arbor/backends/multicore/fvm.hpp
@@ -45,6 +45,7 @@ struct backend {
using sample_event_stream = arb::multicore::sample_event_stream;
using shared_state = arb::multicore::shared_state;
+ using ion_state = arb::multicore::ion_state;
static threshold_watcher voltage_watcher(
const shared_state& state,
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index 7eb69f63..a950e5d0 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -272,6 +272,17 @@ namespace impl {
ctr.push_back(x+1==0? x: offset+x);
}
}
+
+ template <typename Container>
+ void append_divs(Container& left, const Container& right) {
+ if (left.empty()) {
+ left = right;
+ }
+ else if (!right.empty()) {
+ append_offset(left, left.back(), tail(right));
+ }
+ };
+
}
// Merge CV geometry lists in-place.
@@ -279,6 +290,7 @@ namespace impl {
cv_geometry& append(cv_geometry& geom, const cv_geometry& right) {
using impl::tail;
using impl::append_offset;
+ using impl::append_divs;
if (!right.n_cell()) {
return geom;
@@ -289,15 +301,6 @@ cv_geometry& append(cv_geometry& geom, const cv_geometry& right) {
return geom;
}
- auto append_divs = [](auto& left, const auto& right) {
- if (left.empty()) {
- left = right;
- }
- else if (!right.empty()) {
- append_offset(left, left.back(), tail(right));
- }
- };
-
auto geom_n_cv = geom.size();
auto geom_n_cell = geom.n_cell();
@@ -432,7 +435,23 @@ fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell
if (D.cv_area[i]>0) {
D.init_membrane_potential[i] /= D.cv_area[i];
D.temperature_K[i] /= D.cv_area[i];
+
+ // If parent is trivial, and there is no grandparent, then we can use values from this CV
+ // to get initial values for the parent. (The other case, when there is a grandparent, is
+ // caught below.)
+
+ if (p!=-1 && D.geometry.cv_parent[p]==-1 && D.cv_area[p]==0) {
+ D.init_membrane_potential[p] = D.init_membrane_potential[i];
+ D.temperature_K[p] = D.temperature_K[i];
+ }
+ }
+ else if (p!=-1) {
+ // Use parent CV to get a sensible initial value for voltage and temp on zero-size CVs.
+
+ D.init_membrane_potential[i] = D.init_membrane_potential[p];
+ D.temperature_K[i] = D.temperature_K[p];
}
+
if (cv_length>0) {
D.diam_um[i] = D.cv_area[i]/(cv_length*math::pi<double>);
}
@@ -671,8 +690,8 @@ fvm_voltage_interpolant fvm_axial_current(const cable_cell& cell, const fvm_cv_d
mcable rr_span = mcable{bid, vrefs.proximal.loc.pos , vrefs.distal.loc.pos};
double rr_conductance = 100/embedding.integrate_ixa(rr_span, D.axial_resistivity[cell_idx].at(bid)); // [µS]
- vi.proximal_coef = -rr_conductance;
- vi.distal_coef = rr_conductance;
+ vi.proximal_coef = rr_conductance;
+ vi.distal_coef = -rr_conductance;
}
return vi;
@@ -686,6 +705,7 @@ fvm_voltage_interpolant fvm_axial_current(const cable_cell& cell, const fvm_cv_d
fvm_mechanism_data& append(fvm_mechanism_data& left, const fvm_mechanism_data& right) {
using impl::append_offset;
+ using impl::append_divs;
fvm_size_type target_offset = left.n_target;
@@ -728,7 +748,10 @@ fvm_mechanism_data& append(fvm_mechanism_data& left, const fvm_mechanism_data& r
}
}
}
+
left.n_target += right.n_target;
+ append_divs(left.target_divs, right.target_divs);
+ arb_assert(left.n_target==left.target_divs.back());
return left;
}
@@ -1166,6 +1189,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
}
}
+ M.target_divs = {0u, M.n_target};
return M;
}
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index e3ade279..72221e6b 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -240,6 +240,7 @@ struct fvm_mechanism_config {
std::vector<value_type> norm_area;
// Synapse target number (point mechanisms only).
+ // For each instance index i, there are multiplicity[i] consecutive entries.
std::vector<index_type> target;
// (Non-global) parameters and parameter values across the mechanism instance.
@@ -275,8 +276,11 @@ struct fvm_mechanism_data {
// Ion config, indexed by ion name.
std::unordered_map<std::string, fvm_ion_config> ions;
- // Total number of targets (point-mechanism points)
+ // Total number of targets (point-mechanism points).
std::size_t n_target = 0;
+
+ // Partitions target numbers by cell.
+ std::vector<std::size_t> target_divs;
};
fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop, const std::vector<cable_cell>& cells, const fvm_cv_discretization& D, const arb::execution_context& ctx={});
diff --git a/arbor/fvm_lowered_cell.hpp b/arbor/fvm_lowered_cell.hpp
index 99d3000e..7cd80f59 100644
--- a/arbor/fvm_lowered_cell.hpp
+++ b/arbor/fvm_lowered_cell.hpp
@@ -4,13 +4,17 @@
#include <vector>
#include <arbor/common_types.hpp>
+#include <arbor/cable_cell.hpp>
#include <arbor/fvm_types.hpp>
+#include <arbor/morph/primitives.hpp>
#include <arbor/recipe.hpp>
+#include <arbor/util/variant.hpp>
#include "backends/event.hpp"
#include "backends/threshold_crossing.hpp"
#include "execution_context.hpp"
#include "sampler_map.hpp"
+#include "util/meta.hpp"
#include "util/range.hpp"
namespace arb {
@@ -23,14 +27,101 @@ struct fvm_integration_result {
// A sample for a probe may be derived from multiple 'raw' sampled
// values from the backend.
-//
-// While supported probes are at this point all simple scalar values,
-// fvm_probe_info will be the class that represents the mapping
-// between a single sample result and the back-end raw probe handles.
+
+// An `fvm_probe_info` object represents the mapping between
+// a sample value (possibly non-scalar) presented to a sampler
+// function, and one or more probe handles that reference data
+// in the FVM back-end.
+
+struct fvm_probe_scalar {
+ probe_handle raw_handles[1] = {nullptr};
+ util::variant<mlocation, cable_probe_point_info> metadata;
+};
+
+struct fvm_probe_interpolated {
+ probe_handle raw_handles[2] = {nullptr, nullptr};
+ double coef[2];
+ mlocation metadata;
+};
+
+struct fvm_probe_multi {
+ std::vector<probe_handle> raw_handles;
+ util::variant<mcable_list, std::vector<cable_probe_point_info>> metadata;
+
+ void shrink_to_fit() {
+ raw_handles.shrink_to_fit();
+ util::visit([](auto& v) { v.shrink_to_fit(); }, metadata);
+ }
+};
+
+struct fvm_probe_weighted_multi {
+ std::vector<probe_handle> raw_handles;
+ std::vector<double> weight;
+ mcable_list metadata;
+
+ void shrink_to_fit() {
+ raw_handles.shrink_to_fit();
+ weight.shrink_to_fit();
+ metadata.shrink_to_fit();
+ }
+};
+
+// Trans-membrane currents require special handling!
+struct fvm_probe_membrane_currents {
+ std::vector<probe_handle> raw_handles; // Voltage per CV.
+ std::vector<mcable> metadata; // Cables from each CV, in CV order.
+
+ std::vector<unsigned> cv_parent; // Parent CV index for each CV.
+ std::vector<double> cv_parent_cond; // Face conductance between CV and parent.
+ std::vector<double> weight; // Area of cable : area of CV.
+ std::vector<unsigned> cv_cables_divs; // Partitions metadata by CV index.
+
+ void shrink_to_fit() {
+ raw_handles.shrink_to_fit();
+ metadata.shrink_to_fit();
+ cv_parent.shrink_to_fit();
+ cv_parent_cond.shrink_to_fit();
+ weight.shrink_to_fit();
+ cv_cables_divs.shrink_to_fit();
+ }
+};
+
+struct missing_probe_info {
+ // dummy data...
+ std::array<probe_handle, 0> raw_handles;
+};
struct fvm_probe_info {
- // nullptr => nothing to probe
- probe_handle raw_handle = nullptr;
+ fvm_probe_info() = default;
+ fvm_probe_info(fvm_probe_scalar p): info(std::move(p)) {}
+ fvm_probe_info(fvm_probe_interpolated p): info(std::move(p)) {}
+ fvm_probe_info(fvm_probe_multi p): info(std::move(p)) {}
+ fvm_probe_info(fvm_probe_weighted_multi p): info(std::move(p)) {}
+ fvm_probe_info(fvm_probe_membrane_currents p): info(std::move(p)) {}
+
+ util::variant<
+ missing_probe_info,
+ fvm_probe_scalar,
+ fvm_probe_interpolated,
+ fvm_probe_multi,
+ fvm_probe_weighted_multi,
+ fvm_probe_membrane_currents
+ > info = missing_probe_info{};
+
+ auto raw_handle_range() const {
+ return util::make_range(
+ util::visit(
+ [](auto& i) -> std::pair<const probe_handle*, const probe_handle*> {
+ using util::data;
+ using util::size;
+ return {data(i.raw_handles), data(i.raw_handles)+size(i.raw_handles)};
+ },
+ info));
+ }
+
+ sample_size_type n_raw() const { return raw_handle_range().size(); }
+
+ explicit operator bool() const { return !util::get_if<missing_probe_info>(info); }
};
// Common base class for FVM implementation on host or gpu back-end.
@@ -43,7 +134,7 @@ struct fvm_lowered_cell {
const recipe& rec,
std::vector<fvm_index_type>& cell_to_intdom,
std::vector<target_handle>& target_handles,
- probe_association_map<probe_handle>& probe_map) = 0;
+ probe_association_map<fvm_probe_info>& probe_map) = 0;
virtual fvm_integration_result integrate(
fvm_value_type tfinal,
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index f35f446e..d289526e 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -20,6 +20,9 @@
#include <arbor/common_types.hpp>
#include <arbor/cable_cell_param.hpp>
#include <arbor/recipe.hpp>
+#include <arbor/util/any.hpp>
+#include <arbor/util/any_visitor.hpp>
+#include <arbor/util/optional.hpp>
#include "builtin_mechanisms.hpp"
#include "execution_context.hpp"
@@ -54,7 +57,7 @@ public:
const recipe& rec,
std::vector<fvm_index_type>& cell_to_intdom,
std::vector<target_handle>& target_handles,
- probe_association_map<probe_handle>& probe_map) override;
+ probe_association_map<fvm_probe_info>& probe_map) override;
fvm_integration_result integrate(
value_type tfinal,
@@ -139,6 +142,7 @@ private:
// Translate cell probe descriptions into probe handles etc.
fvm_probe_info resolve_probe_address(
+ const std::vector<cable_cell>& cells,
std::size_t cell_idx,
const util::any& paddr,
const fvm_cv_discretization& D,
@@ -352,7 +356,7 @@ void fvm_lowered_cell_impl<Backend>::initialize(
const recipe& rec,
std::vector<fvm_index_type>& cell_to_intdom,
std::vector<target_handle>& target_handles,
- probe_association_map<probe_handle>& probe_map)
+ probe_association_map<fvm_probe_info>& probe_map)
{
using util::any_cast;
using util::count_along;
@@ -483,15 +487,16 @@ void fvm_lowered_cell_impl<Backend>::initialize(
layout.weight[i] = 1000/D.cv_area[cv];
// (builtin stimulus, for example, has no targets)
+ if (config.target.empty()) continue;
- if (!config.target.empty()) {
- if(!config.multiplicity.empty()) {
- for (auto j: make_span(multiplicity_part[i])) {
- target_handles[config.target[j]] = target_handle(mech_id, i, cv_to_intdom[cv]);
- }
- } else {
- target_handles[config.target[i]] = target_handle(mech_id, i, cv_to_intdom[cv]);
- };
+ target_handle handle(mech_id, i, cv_to_intdom[cv]);
+ if (config.multiplicity.empty()) {
+ target_handles[config.target[i]] = handle;
+ }
+ else {
+ for (auto j: make_span(multiplicity_part[i])) {
+ target_handles[config.target[j]] = handle;
+ }
}
}
break;
@@ -539,10 +544,10 @@ void fvm_lowered_cell_impl<Backend>::initialize(
for (cell_lid_type j: make_span(rec.num_probes(gid))) {
probe_info pi = rec.get_probe({gid, j});
- fvm_probe_info info = resolve_probe_address(cell_idx, pi.address, D, mech_data, target_handles, mechptr_by_name);
+ fvm_probe_info info = resolve_probe_address(cells, cell_idx, pi.address, D, mech_data, target_handles, mechptr_by_name);
- if (info.raw_handle) {
- probe_map.insert({pi.id, {info.raw_handle, pi.tag}});
+ if (info) {
+ probe_map.insert({pi.id, {std::move(info), pi.tag}});
}
}
}
@@ -648,8 +653,56 @@ fvm_size_type fvm_lowered_cell_impl<Backend>::fvm_intdom(
return intdom_id;
}
+// Resolution of probe addresses into a specific fvm_probe_info draws upon data
+// from the cable cell, the discretization, the target handle map, and the
+// back-end shared state.
+//
+// `resolve_probe_address` collates this data into a `probe_resolution_data`
+// struct which is then passed on to the specific resolution procedure
+// determined by the type of the user-supplied probe address.
+
+template <typename Backend>
+struct probe_resolution_data {
+ typename Backend::shared_state* state;
+ const cable_cell& cell;
+ const std::size_t cell_idx;
+ const fvm_cv_discretization& D;
+ const fvm_mechanism_data& M;
+ const std::vector<target_handle>& handles;
+ const std::unordered_map<std::string, mechanism*>& mech_instance_by_name;
+
+ // Backend state data for a given mechanism and state variable.
+ const fvm_value_type* mechanism_state(const std::string& name, const std::string& state_var) const {
+ mechanism* m = util::value_by_key(mech_instance_by_name, name).value_or(nullptr);
+ if (!m) return nullptr;
+
+ const fvm_value_type* data = Backend::mechanism_field_data(m, state_var);
+ if (!data) throw cable_cell_error("no state variable '"+state_var+"' in mechanism '"+name+"'");
+
+ return data;
+ }
+
+ // Extent of density mechanism on cell.
+ mextent mechanism_support(const std::string& name) const {
+ auto& mech_map = cell.region_assignments().template get<mechanism_desc>();
+ auto opt_mm = util::value_by_key(mech_map, name);
+
+ return opt_mm? opt_mm->support(): mextent{};
+ };
+
+ // Index into ion data from location.
+ util::optional<fvm_index_type> ion_location_index(const std::string& ion, mlocation loc) const {
+ if (state->ion_data.count(ion)) {
+ return util::binary_search_index(M.ions.at(ion).cv,
+ fvm_index_type(D.geometry.location_cv(cell_idx, loc, cv_prefer::cv_nonempty)));
+ }
+ return util::nullopt;
+ }
+};
+
template <typename Backend>
fvm_probe_info fvm_lowered_cell_impl<Backend>::resolve_probe_address(
+ const std::vector<cable_cell>& cells,
std::size_t cell_idx,
const util::any& paddr,
const fvm_cv_discretization& D,
@@ -657,102 +710,310 @@ fvm_probe_info fvm_lowered_cell_impl<Backend>::resolve_probe_address(
const std::vector<target_handle>& handles,
const std::unordered_map<std::string, mechanism*>& mech_instance_by_name)
{
- using util::any_cast;
- using util::value_by_key;
- using util::binary_search_index;
- using namespace cv_prefer;
-
- // Probe address can be one of a number of cable cell-specific
- // probe types; dispatch on type and compute probe info accordingly.
- //
- // Mechanisms or ions not being instantiated at a particular address
- // is not treated as an error; a request for a state variable that
- // does not exist on a mechanism is, however.
-
- auto location_cv = [&](auto* p, cv_prefer::type prefer) -> fvm_index_type {
- return D.geometry.location_cv(cell_idx, p->location, prefer);
- };
+ probe_resolution_data<Backend> prd{
+ state_.get(), cells[cell_idx], cell_idx, D, M, handles, mech_instance_by_name};
+
+ using V = util::any_visitor<
+ cable_probe_membrane_voltage,
+ cable_probe_membrane_voltage_cell,
+ cable_probe_axial_current,
+ cable_probe_total_ion_current_density,
+ cable_probe_total_ion_current_cell,
+ cable_probe_total_current_cell,
+ cable_probe_density_state,
+ cable_probe_density_state_cell,
+ cable_probe_point_state,
+ cable_probe_point_state_cell,
+ cable_probe_ion_current_density,
+ cable_probe_ion_current_cell,
+ cable_probe_ion_int_concentration,
+ cable_probe_ion_int_concentration_cell,
+ cable_probe_ion_ext_concentration,
+ cable_probe_ion_ext_concentration_cell>;
+
+ auto visitor = util::overload(
+ [&prd](auto& probe_addr) { return resolve_probe(probe_addr, prd); },
+ [] { return throw cable_cell_error("unrecognized probe type"), fvm_probe_info{}; });
+
+ return V::visit(visitor, paddr);
+}
- struct mechanism_data {
- const fvm_value_type* state_data = nullptr;
- unsigned id;
- explicit operator bool() const { return state_data; }
- };
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_membrane_voltage& p, probe_resolution_data<B>& R) {
+ const fvm_value_type* data = R.state->voltage.data();
+ fvm_voltage_interpolant in = fvm_interpolate_voltage(R.cell, R.D, R.cell_idx, p.location);
- auto lookup_mechanism_data = [&](auto* p) -> mechanism_data {
- if (mechanism* m = value_by_key(mech_instance_by_name, p->mechanism).value_or(nullptr)) {
- const fvm_value_type* data = Backend::mechanism_field_data(m, p->state);
- if (!data) {
- throw cable_cell_error("no state variable '"+p->state+"' in mechanism '"+p->mechanism+"'");
- }
- return {data, m->mechanism_id()};
- }
- return {};
- };
+ return fvm_probe_interpolated{
+ {data+in.proximal_cv, data+in.distal_cv},
+ {in.proximal_coef, in.distal_coef},
+ p.location};
+}
- auto ion_location_index = [&](auto* p) -> util::optional<fvm_index_type> {
- if (state_->ion_data.count(p->ion)) {
- auto cv = location_cv(p, cv_nonempty);
- return util::binary_search_index(M.ions.at(p->ion).cv, cv);
- }
- return util::nullopt;
- };
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_membrane_voltage_cell& p, probe_resolution_data<B>& R) {
+ fvm_probe_multi r;
+ mcable_list cables;
- if (auto *p = any_cast<cell_probe_membrane_voltage>(&paddr)) {
- const fvm_value_type* src = state_->voltage.data() + location_cv(p, cv_empty);
- return fvm_probe_info{src};
- }
- else if (auto* p = any_cast<cell_probe_total_ionic_current_density>(&paddr)) {
- const fvm_value_type* src = state_->current_density.data() + location_cv(p, cv_nonempty);
- return fvm_probe_info{src};
+ for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) {
+ const double* ptr = R.state->voltage.data()+cv;
+ for (auto cable: R.D.geometry.cables(cv)) {
+ r.raw_handles.push_back(ptr);
+ cables.push_back(cable);
+ }
}
- else if (auto* p = any_cast<cell_probe_density_state>(&paddr)) {
- if (auto data = lookup_mechanism_data(p)) {
- auto cv = location_cv(p, cv_nonempty);
- if (auto opt_i = binary_search_index(M.mechanisms.at(p->mechanism).cv, cv)) {
- return fvm_probe_info{data.state_data+opt_i.value()};
+ r.metadata = std::move(cables);
+ r.shrink_to_fit();
+ return r;
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_axial_current& p, probe_resolution_data<B>& R) {
+ const fvm_value_type* data = R.state->voltage.data();
+ fvm_voltage_interpolant in = fvm_axial_current(R.cell, R.D, R.cell_idx, p.location);
+
+ return fvm_probe_interpolated{
+ {data+in.proximal_cv, data+in.distal_cv},
+ {in.proximal_coef, in.distal_coef},
+ p.location};
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_total_ion_current_density& p, probe_resolution_data<B>& R) {
+ return fvm_probe_scalar{
+ {R.state->current_density.data() + R.D.geometry.location_cv(R.cell_idx, p.location, cv_prefer::cv_nonempty)},
+ p.location};
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_total_ion_current_cell& p, probe_resolution_data<B>& R) {
+ fvm_probe_weighted_multi r;
+
+ for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) {
+ const double* ptr = R.state->current_density.data()+cv;
+ for (auto cable: R.D.geometry.cables(cv)) {
+ double area = R.cell.embedding().integrate_area(cable); // [µm²]
+ if (area>0) {
+ r.raw_handles.push_back(ptr);
+ r.weight.push_back(0.001*area); // Scale from [µm²·A/m²] to [nA].
+ r.metadata.push_back(cable);
}
}
- return fvm_probe_info{};
}
- else if (auto* p = any_cast<cell_probe_point_state>(&paddr)) {
- if (p->target>=handles.size()) {
- return fvm_probe_info{};
- }
+ r.shrink_to_fit();
+ return r;
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_total_current_cell& p, probe_resolution_data<B>& R) {
+ fvm_probe_membrane_currents r;
+
+ auto cell_cv_ival = R.D.geometry.cell_cv_interval(R.cell_idx);
+ auto cv0 = cell_cv_ival.first;
+
+ util::assign(r.cv_parent, util::transform_view(util::subrange_view(R.D.geometry.cv_parent, cell_cv_ival),
+ [cv0](auto cv) { return cv+1==0? cv: cv-cv0; }));
+ util::assign(r.cv_parent_cond, util::subrange_view(R.D.face_conductance, cell_cv_ival));
+
+ r.cv_cables_divs = {0};
+ for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) {
+ r.raw_handles.push_back(R.state->voltage.data()+cv);
+ double oo_cv_area = R.D.cv_area[cv]>0? 1./R.D.cv_area[cv]: 0;
- if (auto data = lookup_mechanism_data(p)) {
- // Confirm mechanism is actually at this target.
- const auto& th = handles[p->target];
- if (th.mech_id == data.id) {
- return fvm_probe_info{data.state_data+th.mech_index};
+ for (auto cable: R.D.geometry.cables(cv)) {
+ double area = R.cell.embedding().integrate_area(cable); // [µm²]
+ if (area>0) {
+ r.weight.push_back(area*oo_cv_area);
+ r.metadata.push_back(cable);
}
}
- return fvm_probe_info{};
+ r.cv_cables_divs.push_back(r.metadata.size());
}
- else if (auto* p = any_cast<cell_probe_ion_current_density>(&paddr)) {
- if (auto opt_i = ion_location_index(p)) {
- const fvm_value_type* data = state_->ion_data.at(p->ion).iX_.data();
- return fvm_probe_info{data+opt_i.value()};
+ r.shrink_to_fit();
+ return r;
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_density_state& p, probe_resolution_data<B>& R) {
+ const fvm_value_type* data = R.mechanism_state(p.mechanism, p.state);
+ if (!data) return {};
+ if (!R.mechanism_support(p.mechanism).intersects(p.location)) return {};
+
+ fvm_index_type cv = R.D.geometry.location_cv(R.cell_idx, p.location, cv_prefer::cv_nonempty);
+ auto opt_i = util::binary_search_index(R.M.mechanisms.at(p.mechanism).cv, cv);
+ if (!opt_i) return {};
+
+ return fvm_probe_scalar{{data+*opt_i}, p.location};
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_density_state_cell& p, probe_resolution_data<B>& R) {
+ fvm_probe_multi r;
+
+ const fvm_value_type* data = R.mechanism_state(p.mechanism, p.state);
+ if (!data) return {};
+
+ mextent support = R.mechanism_support(p.mechanism);
+ auto& mech_cvs = R.M.mechanisms.at(p.mechanism).cv;
+ mcable_list cables;
+
+ for (auto i: util::count_along(mech_cvs)) {
+ auto cv = mech_cvs[i];
+ auto cv_cables = R.D.geometry.cables(cv);
+ mextent cv_extent = mcable_list(cv_cables.begin(), cv_cables.end());
+ for (auto cable: intersect(cv_extent, support)) {
+ if (cable.prox_pos==cable.dist_pos) continue;
+
+ r.raw_handles.push_back(data+i);
+ cables.push_back(cable);
}
- return fvm_probe_info{};
}
- else if (auto* p = any_cast<cell_probe_ion_int_concentration>(&paddr)) {
- if (auto opt_i = ion_location_index(p)) {
- const fvm_value_type* data = state_->ion_data.at(p->ion).Xi_.data();
- return fvm_probe_info{data+opt_i.value()};
- }
- return fvm_probe_info{};
+ r.metadata = std::move(cables);
+ r.shrink_to_fit();
+ return r;
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_point_state& p, probe_resolution_data<B>& R) {
+ arb_assert(R.handles.size()==R.M.target_divs.back());
+ arb_assert(R.handles.size()==R.M.n_target);
+
+ const fvm_value_type* data = R.mechanism_state(p.mechanism, p.state);
+ if (!data) return {};
+
+ // Convert cell-local target number to cellgroup target number.
+ auto cg_target = p.target + R.M.target_divs[R.cell_idx];
+ if (cg_target>=R.M.target_divs.at(R.cell_idx+1)) return {};
+
+ if (R.handles[cg_target].mech_id!=R.mech_instance_by_name.at(p.mechanism)->mechanism_id()) return {};
+ auto mech_index = R.handles[cg_target].mech_index;
+
+ auto& multiplicity = R.M.mechanisms.at(p.mechanism).multiplicity;
+ auto& placed_instances = R.cell.synapses().at(p.mechanism);
+
+ auto opt_i = util::binary_search_index(placed_instances, p.target, [](auto& item) { return item.lid; });
+ if (!opt_i) throw arbor_internal_error("inconsistent mechanism state");
+ mlocation loc = placed_instances[*opt_i].loc;
+
+ cable_probe_point_info metadata{p.target, multiplicity.empty()? 1u: multiplicity.at(mech_index), loc};
+
+ return fvm_probe_scalar{{data+mech_index}, metadata};
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_point_state_cell& p, probe_resolution_data<B>& R) {
+ const fvm_value_type* data = R.mechanism_state(p.mechanism, p.state);
+ if (!data) return {};
+
+ unsigned id = R.mech_instance_by_name.at(p.mechanism)->mechanism_id();
+ auto& multiplicity = R.M.mechanisms.at(p.mechanism).multiplicity;
+ auto& placed_instances = R.cell.synapses().at(p.mechanism);
+
+ std::size_t cell_targets_base = R.M.target_divs[R.cell_idx];
+ std::size_t cell_targets_end = R.M.target_divs[R.cell_idx+1];
+
+ fvm_probe_multi r;
+ std::vector<cable_probe_point_info> metadata;
+
+ for (auto target: util::make_span(cell_targets_base, cell_targets_end)) {
+ if (R.handles[target].mech_id!=id) continue;
+
+ auto mech_index = R.handles[target].mech_index;
+ r.raw_handles.push_back(data+mech_index);
+
+ auto cell_target = target-cell_targets_base; // Convert to cell-local target index.
+
+ auto opt_i = util::binary_search_index(placed_instances, cell_target, [](auto& item) { return item.lid; });
+ if (!opt_i) throw arbor_internal_error("inconsistent mechanism state");
+ mlocation loc = placed_instances[*opt_i].loc;
+
+ metadata.push_back(cable_probe_point_info{
+ cell_lid_type(cell_target), multiplicity.empty()? 1u: multiplicity.at(mech_index), loc});
}
- else if (auto* p = any_cast<cell_probe_ion_ext_concentration>(&paddr)) {
- if (auto opt_i = ion_location_index(p)) {
- const fvm_value_type* data = state_->ion_data.at(p->ion).Xo_.data();
- return fvm_probe_info{data+opt_i.value()};
+
+ r.metadata = std::move(metadata);
+ r.shrink_to_fit();
+ return r;
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_ion_current_density& p, probe_resolution_data<B>& R) {
+ auto opt_i = R.ion_location_index(p.ion, p.location);
+ if (!opt_i) return {};
+
+ return fvm_probe_scalar{{R.state->ion_data.at(p.ion).iX_.data()+*opt_i}, p.location};
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_ion_current_cell& p, probe_resolution_data<B>& R) {
+ if (!R.state->ion_data.count(p.ion)) return {};
+
+ auto& ion_cvs = R.M.ions.at(p.ion).cv;
+ const fvm_value_type* src = R.state->ion_data.at(p.ion).iX_.data();
+
+ fvm_probe_weighted_multi r;
+ for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) {
+ auto opt_i = util::binary_search_index(ion_cvs, cv);
+ if (!opt_i) continue;
+
+ const double* ptr = src+*opt_i;
+ for (auto cable: R.D.geometry.cables(cv)) {
+ double area = R.cell.embedding().integrate_area(cable); // [µm²]
+ if (area>0) {
+ r.raw_handles.push_back(ptr);
+ r.weight.push_back(0.001*area); // Scale from [µm²·A/m²] to [nA].
+ r.metadata.push_back(cable);
+ }
}
- return fvm_probe_info{};
}
- else {
- throw cable_cell_error("unrecognized probe address type");
+ r.metadata.shrink_to_fit();
+ return r;
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_ion_int_concentration& p, probe_resolution_data<B>& R) {
+ auto opt_i = R.ion_location_index(p.ion, p.location);
+ if (!opt_i) return {};
+
+ return fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xi_.data()+*opt_i}, p.location};
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_ion_ext_concentration& p, probe_resolution_data<B>& R) {
+ auto opt_i = R.ion_location_index(p.ion, p.location);
+ if (!opt_i) return {};
+
+ return fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xo_.data()+*opt_i}, p.location};
+}
+
+// Common implementation for int and ext concentrations across whole cell:
+template <typename B>
+fvm_probe_info resolve_ion_conc_common(const std::vector<fvm_index_type>& ion_cvs, const fvm_value_type* src, probe_resolution_data<B>& R) {
+ fvm_probe_multi r;
+ mcable_list cables;
+
+ for (auto i: util::count_along(ion_cvs)) {
+ for (auto cable: R.D.geometry.cables(ion_cvs[i])) {
+ if (cable.prox_pos!=cable.dist_pos) {
+ r.raw_handles.push_back(src+i);
+ cables.push_back(cable);
+ }
+ }
}
+ r.metadata = std::move(cables);
+ r.shrink_to_fit();
+ return r;
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_ion_int_concentration_cell& p, probe_resolution_data<B>& R) {
+ if (!R.state->ion_data.count(p.ion)) return {};
+ return resolve_ion_conc_common<B>(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xi_.data(), R);
+}
+
+template <typename B>
+fvm_probe_info resolve_probe(const cable_probe_ion_ext_concentration_cell& p, probe_resolution_data<B>& R) {
+ if (!R.state->ion_data.count(p.ion)) return {};
+ return resolve_ion_conc_common<B>(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xo_.data(), R);
}
} // namespace arb
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index ed9d3a92..d942bb47 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -29,51 +29,153 @@ struct lid_range {
begin(b), end(e) {}
};
-// Each kind of probe has its own type for representing its address,
-// as below:
+// `cable_sample_range` is simply a pair of `const double*` pointers describing the sequence
+// of double values associated with the cell-wide sample.
-// Voltage estimate [mV] at `location`, possibly interpolated.
-struct cell_probe_membrane_voltage {
+using cable_sample_range = std::pair<const double*, const double*>;
+
+
+// Each kind of probe has its own type for representing its address, as below.
+//
+// Probe address specifications can be for _scalar_ data, associated with a fixed
+// location or synapse on a cell, or _vector_ data, associated with multiple
+// sites or sub-sections of a cell.
+//
+// Sampler functions receive an `any_ptr` to sampled data. The underlying pointer
+// type is a const pointer to:
+// * `double` for scalar data;
+// * `cable_sample_rang*` for vector data (see definition above).
+//
+// The metadata associated with a probe is also passed to a sampler via an `any_ptr`;
+// the underlying pointer will be a const pointer to one of the following metadata types:
+// * `mlocation` for most scalar queries;
+// * `cable_probe_point_info` for point mechanism state queries;
+// * `mcable_list` for most vector queries;
+// * `std::vector<cable_probe_point_info>` for cell-wide point mechanism state queries.
+
+// Metadata for point process probes.
+struct cable_probe_point_info {
+ cell_lid_type target; // Target number of point process instance on cell.
+ unsigned multiplicity; // Number of combined instances at this site.
+ mlocation loc; // Point on cell morphology where instance is placed.
+};
+
+// Voltage estimate [mV] at `location`, interpolated.
+// Sample value type: `double`
+// Sample metadata type: `mlocation`
+struct cable_probe_membrane_voltage {
mlocation location;
};
-// Total current density [A/m²] across membrane excluding capacitive current at `location`.
-struct cell_probe_total_ionic_current_density {
+// Voltage estimate [mV], reported against each cable in each control volume. Not interpolated.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mcable_list`
+struct cable_probe_membrane_voltage_cell {};
+
+// Axial current estimate [nA] at `location`, interpolated.
+// Sample value type: `double`
+// Sample metadata type: `mlocation`
+struct cable_probe_axial_current {
mlocation location;
};
+// Total current density [A/m²] across membrane _excluding_ capacitive current at `location`.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mlocation`
+struct cable_probe_total_ion_current_density {
+ mlocation location;
+};
+
+// Total ionic current [nA] across membrance _excluding_ capacitive current across components of the cell.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mcable_list`
+struct cable_probe_total_ion_current_cell {};
+
+// Total membrance current [nA] across components of the cell.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mcable_list`
+struct cable_probe_total_current_cell {};
+
// Value of state variable `state` in density mechanism `mechanism` in CV at `location`.
-struct cell_probe_density_state {
+// Sample value type: `double`
+// Sample metadata type: `mlocation`
+struct cable_probe_density_state {
mlocation location;
std::string mechanism;
std::string state;
};
+// Value of state variable `state` in density mechanism `mechanism` across components of the cell.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mcable_list`
+struct cable_probe_density_state_cell {
+ std::string mechanism;
+ std::string state;
+};
+
// Value of state variable `key` in point mechanism `source` at target `target`.
-struct cell_probe_point_state {
+// Sample value type: `double`
+// Sample metadata type: `cable_probe_point_info`
+struct cable_probe_point_state {
cell_lid_type target;
std::string mechanism;
std::string state;
};
+// Value of state variable `key` in point mechanism `source` at every target with this mechanism.
+// Metadata has one entry of type cable_probe_point_info for each matched (possibly coalesced) instance.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `std::vector<cable_probe_point_info>`
+struct cable_probe_point_state_cell {
+ std::string mechanism;
+ std::string state;
+};
+
// Current density [A/m²] across membrane attributed to the ion `source` at `location`.
-struct cell_probe_ion_current_density {
+// Sample value type: `double`
+// Sample metadata type: `mlocation`
+struct cable_probe_ion_current_density {
mlocation location;
std::string ion;
};
-// Ionic internal concentration [mmol/L] of ion `source` at `location`, possibly interpolated.
-struct cell_probe_ion_int_concentration {
+// Total ionic current [nA] attributed to the ion `source` across components of the cell.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mcable_list`
+struct cable_probe_ion_current_cell {
+ std::string ion;
+};
+
+// Ionic internal concentration [mmol/L] of ion `source` at `location`.
+// Sample value type: `double`
+// Sample metadata type: `mlocation`
+struct cable_probe_ion_int_concentration {
mlocation location;
std::string ion;
};
-// Ionic external concentration [mmol/L] of ion `source` at `location`, possibly interpolated.
-struct cell_probe_ion_ext_concentration {
+// Ionic internal concentration [mmol/L] of ion `source` across components of the cell.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mcable_list`
+struct cable_probe_ion_int_concentration_cell {
+ std::string ion;
+};
+
+// Ionic external concentration [mmol/L] of ion `source` at `location`.
+// Sample value type: `double`
+// Sample metadata type: `mlocation`
+struct cable_probe_ion_ext_concentration {
mlocation location;
std::string ion;
};
+// Ionic external concentration [mmol/L] of ion `source` across components of the cell.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mcable_list`
+struct cable_probe_ion_ext_concentration_cell {
+ std::string ion;
+};
+
// Forward declare the implementation, for PIMPL.
struct cable_cell_impl;
@@ -96,6 +198,7 @@ struct placed {
T item;
};
+// Note: lid fields of elements of mlocation_map used in cable_cell are strictly increasing.
template <typename T>
using mlocation_map = std::vector<placed<T>>;
diff --git a/arbor/include/arbor/sampling.hpp b/arbor/include/arbor/sampling.hpp
index fc29c99d..f917e0c6 100644
--- a/arbor/include/arbor/sampling.hpp
+++ b/arbor/include/arbor/sampling.hpp
@@ -21,7 +21,13 @@ struct sample_record {
util::any_ptr data;
};
-using sampler_function = std::function<void (cell_member_type, probe_tag, std::size_t, const sample_record*)>;
+using sampler_function = std::function<
+ void (cell_member_type, // probe id
+ probe_tag, // probe tag associated with probe id
+ util::any_ptr, // pointer to constant metadata
+ std::size_t, // number of sample records
+ const sample_record* // pointer to first sample record
+ )>;
using sampler_association_handle = std::size_t;
diff --git a/arbor/include/arbor/simple_sampler.hpp b/arbor/include/arbor/simple_sampler.hpp
index 2814eb25..a95ae8c8 100644
--- a/arbor/include/arbor/simple_sampler.hpp
+++ b/arbor/include/arbor/simple_sampler.hpp
@@ -12,6 +12,7 @@
#include <arbor/common_types.hpp>
#include <arbor/sampling.hpp>
#include <arbor/util/any_ptr.hpp>
+#include <arbor/util/optional.hpp>
namespace arb {
@@ -21,32 +22,96 @@ struct trace_entry {
V v;
};
+// Trace data wraps a std::vector of trace_entry, optionally with
+// a copy of the metadata associated with a probe.
+
+template <typename V, typename Meta = void>
+struct trace_data: public std::vector<trace_entry<V>> {
+ util::optional<Meta> metadata;
+};
+
+template <typename V>
+struct trace_data<V, void>: public std::vector<trace_entry<V>> {};
+
+// Add a bit of smarts for collecting variable-length samples which are
+// passed back as a pair of pointers describing a range; these can
+// be used to populate a trace of vectors.
+
+template <typename V>
+struct trace_push_back {
+ template <typename Meta>
+ static bool push_back(trace_data<V, Meta>& trace, const sample_record& rec) {
+ if (auto p = util::any_cast<const V*>(rec.data)) {
+ trace.push_back({rec.time, *p});
+ return true;
+ }
+ return false;
+ }
+};
+
template <typename V>
-using trace_data = std::vector<trace_entry<V>>;
+struct trace_push_back<std::vector<V>> {
+ template <typename Meta>
+ static bool push_back(trace_data<std::vector<V>, Meta>& trace, const sample_record& rec) {
+ if (auto p = util::any_cast<const std::vector<V>*>(rec.data)) {
+ trace.push_back({rec.time, *p});
+ return true;
+ }
+ else if (auto p = util::any_cast<const std::pair<const V*, const V*>*>(rec.data)) {
+ trace.push_back({rec.time, std::vector<V>(p->first, p->second)});
+ return true;
+ }
+ return false;
+ }
+};
-template <typename V, typename = std::enable_if_t<std::is_trivially_copyable<V>::value>>
+template <typename V, typename Meta = void>
class simple_sampler {
public:
- explicit simple_sampler(trace_data<V>& trace): trace_(trace) {}
+ explicit simple_sampler(trace_data<V, Meta>& trace): trace_(trace) {}
- void operator()(cell_member_type probe_id, probe_tag tag, std::size_t n, const sample_record* recs) {
- for (std::size_t i = 0; i<n; ++i) {
- if (auto p = util::any_cast<const V*>(recs[i].data)) {
- trace_.push_back({recs[i].time, *p});
+ void operator()(cell_member_type probe_id, probe_tag tag, util::any_ptr meta, std::size_t n, const sample_record* recs) {
+ // TODO: C++17 use if constexpr to test for Meta = void case.
+ if (meta && !trace_.metadata) {
+ if (auto m = util::any_cast<const Meta*>(meta)) {
+ trace_.metadata = *m;
}
else {
+ throw std::runtime_error("unexpected metadata type in simple_sampler");
+ }
+ }
+
+ for (std::size_t i = 0; i<n; ++i) {
+ if (!trace_push_back<V>::push_back(trace_, recs[i])) {
throw std::runtime_error("unexpected sample type in simple_sampler");
}
}
}
private:
- trace_data<V>& trace_;
+ trace_data<V, Meta>& trace_;
};
template <typename V>
-inline simple_sampler<V> make_simple_sampler(trace_data<V>& trace) {
- return simple_sampler<V>(trace);
+class simple_sampler<V, void> {
+public:
+ explicit simple_sampler(trace_data<V, void>& trace): trace_(trace) {}
+
+ void operator()(cell_member_type probe_id, probe_tag tag, util::any_ptr meta, std::size_t n, const sample_record* recs) {
+ for (std::size_t i = 0; i<n; ++i) {
+ if (!trace_push_back<V>::push_back(trace_, recs[i])) {
+ throw std::runtime_error("unexpected sample type in simple_sampler");
+ }
+ }
+ }
+
+private:
+ trace_data<V, void>& trace_;
+};
+
+template <typename V, typename Meta>
+inline simple_sampler<V, Meta> make_simple_sampler(trace_data<V, Meta>& trace) {
+ return simple_sampler<V, Meta>(trace);
}
} // namespace arb
diff --git a/arbor/include/arbor/util/any_visitor.hpp b/arbor/include/arbor/util/any_visitor.hpp
new file mode 100644
index 00000000..348bdaa8
--- /dev/null
+++ b/arbor/include/arbor/util/any_visitor.hpp
@@ -0,0 +1,155 @@
+#pragma once
+
+// Provides `overload` function for constructing a functor with overloaded
+// `operator()` from a number of functions or function objects.
+//
+// Provides the `any_visitor` class, which will call a provided functional
+// with the contained value in a `util::any` if it is one of a fixed set
+// of types.
+
+#include <utility>
+#include <type_traits>
+
+#include <arbor/util/any.hpp>
+
+namespace arb {
+namespace util {
+
+namespace impl {
+
+template <typename> using void_t = void; // TODO: C++17 use std::void_t.
+
+template <typename X, typename Y>
+struct propagate_qualifier { using type = Y; };
+
+template <typename X, typename Y>
+struct propagate_qualifier<const X, Y> { using type = const Y; };
+
+template <typename X, typename Y>
+struct propagate_qualifier<X&, Y> { using type = Y&; };
+
+template <typename X, typename Y>
+struct propagate_qualifier<const X&, Y> { using type = const Y&; };
+
+template <typename X, typename Y>
+using propagate_qualifier_t = typename propagate_qualifier<X, Y>::type;
+
+} // namespace impl
+
+// A type `any_visitor<A, B, ...>` has one public static method
+// `visit(f, a)` where `a` is a possibly const lvalue or rvalue util::any,
+// and `f` is a functional object or function pointer.
+//
+// If `a` contains a value of any of the types `A, B, ...`, `f` will
+// be called with that value. If `a` is an lvalue, it will be passed by
+// lvaue reference; otherwise it will be moved.
+//
+// If `a` contains no value, or a value not in the type list `A, B, ...`,
+// then it will evaluate `f()` if it is defined, or else throw a
+// `bad_any_cast` exception.
+
+template <typename, typename...> struct any_visitor;
+
+template <typename T>
+struct any_visitor<T> {
+ template <typename F, typename = void>
+ struct invoke_or_throw {
+ template <typename A>
+ static auto visit(F&& f, A&& a) {
+ using Q = impl::propagate_qualifier_t<A, T>;
+ return util::any_cast<T>(&a)?
+ std::forward<F>(f)(util::any_cast<Q&&>(std::forward<A>(a))):
+ throw ::arb::util::bad_any_cast();
+ }
+ };
+
+ template <typename F>
+ struct invoke_or_throw<F, impl::void_t<decltype(std::declval<F>()())>> {
+ template <typename A>
+ static auto visit(F&& f, A&& a) {
+ using Q = impl::propagate_qualifier_t<A, T>;
+ return util::any_cast<T>(&a)?
+ std::forward<F>(f)(util::any_cast<Q&&>(std::forward<A>(a))):
+ std::forward<F>(f)();
+ }
+ };
+
+ template <typename F, typename A,
+ typename = std::enable_if_t<std::is_same<util::any, std::decay_t<A>>::value>
+ >
+ static auto visit(F&& f, A&& a) {
+ return invoke_or_throw<F>::visit(std::forward<F>(f), std::forward<A>(a));
+ }
+};
+
+template <typename T, typename U, typename... Rest>
+struct any_visitor<T, U, Rest...> {
+ template <typename F, typename A,
+ typename = std::enable_if_t<std::is_same<util::any, std::decay_t<A>>::value>
+ >
+ static auto visit(F&& f, A&& a) {
+ using Q = impl::propagate_qualifier_t<A, T>;
+ return util::any_cast<T>(&a)?
+ std::forward<F>(f)(util::any_cast<Q&&>(std::forward<A>(a))):
+ any_visitor<U, Rest...>::visit(std::forward<F>(f), std::forward<A>(a));
+ }
+};
+
+namespace impl {
+
+template <typename F, typename... A>
+struct invocable_impl {
+ template <typename G, typename = void>
+ struct test: std::false_type {};
+
+ template <typename G>
+ struct test<G, void_t<decltype(std::declval<G>()(std::declval<A>()...))>>: std::true_type {};
+
+ using type = typename test<F>::type;
+};
+
+template <typename F, typename... A>
+using invocable = typename invocable_impl<F, A...>::type;
+
+template <typename, typename...> struct overload_impl {};
+
+template <typename F1>
+struct overload_impl<F1> {
+ F1 f_;
+
+ overload_impl(F1&& f1): f_(std::forward<F1>(f1)) {}
+
+ template <typename... A, std::enable_if_t<invocable<F1, A...>::value, int> = 0>
+ decltype(auto) operator()(A&&... a) { return f_(std::forward<A>(a)...); }
+};
+
+template <typename F1, typename F2, typename... Fn>
+struct overload_impl<F1, F2, Fn...>: overload_impl<F2, Fn...> {
+ F1 f_;
+
+ overload_impl(F1&& f1, F2&& f2, Fn&&... fn):
+ overload_impl<F2, Fn...>(std::forward<F2>(f2), std::forward<Fn>(fn)...),
+ f_(std::forward<F1>(f1)) {}
+
+ template <typename... A, std::enable_if_t<invocable<F1, A...>::value, int> = 0>
+ decltype(auto) operator()(A&&... a) { return f_(std::forward<A>(a)...); }
+
+ template <typename... A, std::enable_if_t<!invocable<F1, A...>::value, int> = 0>
+ decltype(auto) operator()(A&&... a) {
+ return overload_impl<F2, Fn...>::operator()(std::forward<A>(a)...);
+ }
+};
+
+} // namespace impl
+
+
+// `overload(f, g, h, ...)` returns a functional object whose `operator()` is overloaded
+// with those of `f`, `g`, `h`, ... in decreasing order of precedence.
+
+template <typename... Fn>
+auto overload(Fn&&... fn) {
+ return impl::overload_impl<Fn...>(std::forward<Fn>(fn)...);
+}
+
+} // namespace util
+} // namespace arb
diff --git a/arbor/include/arbor/util/variant.hpp b/arbor/include/arbor/util/variant.hpp
index 09700a1c..5b2715b8 100644
--- a/arbor/include/arbor/util/variant.hpp
+++ b/arbor/include/arbor/util/variant.hpp
@@ -431,7 +431,7 @@ struct variant {
auto get_if() noexcept { return get_if<I>(); }
template <std::size_t I, typename = std::enable_if_t<(I<sizeof...(T))>, typename X = type_select_t<I, T...>>
- const X* get_if() const noexcept { return which_==I? data_ptr<>(): nullptr; }
+ const X* get_if() const noexcept { return which_==I? data_ptr<X>(): nullptr; }
template <typename X, std::size_t I = type_index<X, T...>::value>
auto get_if() const noexcept { return get_if<I>(); }
@@ -502,9 +502,6 @@ struct variant_alternative<I, variant<T...>> { using type = type_select_t<I, T..
template <std::size_t I, typename... T>
struct variant_alternative<I, const variant<T...>> { using type = std::add_const_t<type_select_t<I, T...>>; };
-template <typename Visitor, typename... Variant>
-using visit_return_t = decltype(std::declval<Visitor>()(std::declval<typename variant_alternative<0, std::remove_volatile_t<std::remove_reference_t<Variant>>>::type>()...));
-
} // namespace detail
template <typename... T>
@@ -560,7 +557,7 @@ decltype(auto) get_if(const variant<T...>& v) noexcept { return v.template get_i
template <typename Visitor, typename Variant>
decltype(auto) visit(Visitor&& f, Variant&& v) {
- using R = detail::visit_return_t<Visitor&&, Variant&&>;
+ using R = decltype(f(get<0>(std::forward<Variant>(v))));
if (v.valueless_by_exception()) throw bad_variant_access{};
std::size_t i = v.index();
diff --git a/arbor/mc_cell_group.cpp b/arbor/mc_cell_group.cpp
index ce142817..0032d62a 100644
--- a/arbor/mc_cell_group.cpp
+++ b/arbor/mc_cell_group.cpp
@@ -4,6 +4,7 @@
#include <arbor/assert.hpp>
#include <arbor/common_types.hpp>
+#include <arbor/cable_cell.hpp>
#include <arbor/sampling.hpp>
#include <arbor/recipe.hpp>
#include <arbor/spike.hpp>
@@ -79,6 +80,245 @@ void mc_cell_group::set_binning_policy(binning_kind policy, time_type bin_interv
binners_.resize(gids_.size(), event_binner(policy, bin_interval));
}
+// Probe-type specific sample data marshalling.
+
+struct sampler_call_info {
+ sampler_function sampler;
+ cell_member_type probe_id;
+ probe_tag tag;
+
+ // Offsets are into lowered cell sample time and event arrays.
+ sample_size_type begin_offset;
+ sample_size_type end_offset;
+};
+
+// Working space for computing and collating data for samplers.
+using fvm_probe_scratch = std::tuple<std::vector<double>, std::vector<cable_sample_range>>;
+
+template <typename VoidFn, typename... A>
+void tuple_foreach(VoidFn&& f, std::tuple<A...>& t) {
+ (void)(int []){(f(std::get<A>(t)), 0)...};
+}
+
+void reserve_scratch(fvm_probe_scratch& scratch, std::size_t n) {
+ tuple_foreach([n](auto& v) { v.reserve(n); }, scratch);
+}
+
+void run_samples(
+ const missing_probe_info&,
+ const sampler_call_info&,
+ const fvm_value_type*,
+ const fvm_value_type*,
+ std::vector<sample_record>&,
+ fvm_probe_scratch&)
+{
+ throw arbor_internal_error("invalid fvm_probe_info in sampler map");
+}
+
+void run_samples(
+ const fvm_probe_scalar& p,
+ const sampler_call_info& sc,
+ const fvm_value_type* raw_times,
+ const fvm_value_type* raw_samples,
+ std::vector<sample_record>& sample_records,
+ fvm_probe_scratch&)
+{
+ // Scalar probes do not need scratch space — provided that the user-presented
+ // sample type (double) matches the raw type (fvm_value_type).
+ static_assert(std::is_same<double, fvm_value_type>::value, "require sample value translation");
+
+ sample_size_type n_sample = sc.end_offset-sc.begin_offset;
+ sample_records.clear();
+ for (auto i = sc.begin_offset; i!=sc.end_offset; ++i) {
+ sample_records.push_back(sample_record{time_type(raw_times[i]), &raw_samples[i]});
+ }
+
+ // Metadata may be an mlocation or cell_lid_type (for point mechanism state probes).
+ util::visit(
+ [&](auto& metadata) { sc.sampler(sc.probe_id, sc.tag, &metadata, n_sample, sample_records.data()); },
+ p.metadata);
+}
+
+void run_samples(
+ const fvm_probe_interpolated& p,
+ const sampler_call_info& sc,
+ const fvm_value_type* raw_times,
+ const fvm_value_type* raw_samples,
+ std::vector<sample_record>& sample_records,
+ fvm_probe_scratch& scratch)
+{
+ constexpr sample_size_type n_raw_per_sample = 2;
+ sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample;
+ arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample);
+
+ auto& tmp = std::get<std::vector<double>>(scratch);
+ tmp.clear();
+ sample_records.clear();
+
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ tmp.push_back(p.coef[0]*raw_samples[offset] + p.coef[1]*raw_samples[offset+1]);
+ }
+
+ const auto& ctmp = tmp;
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ sample_records.push_back(sample_record{time_type(raw_times[offset]), &ctmp[j]});
+ }
+
+ sc.sampler(sc.probe_id, sc.tag, &p.metadata, n_sample, sample_records.data());
+}
+
+void run_samples(
+ const fvm_probe_multi& p,
+ const sampler_call_info& sc,
+ const fvm_value_type* raw_times,
+ const fvm_value_type* raw_samples,
+ std::vector<sample_record>& sample_records,
+ fvm_probe_scratch& scratch)
+{
+ const sample_size_type n_raw_per_sample = p.raw_handles.size();
+ sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample;
+ arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample);
+
+ auto& sample_ranges = std::get<std::vector<cable_sample_range>>(scratch);
+ sample_ranges.clear();
+ sample_records.clear();
+
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ sample_ranges.push_back({raw_samples+offset, raw_samples+offset+n_raw_per_sample});
+ }
+
+ const auto& csample_ranges = sample_ranges;
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ sample_records.push_back(sample_record{time_type(raw_times[offset]), &csample_ranges[j]});
+ }
+
+ // Metadata may be an mlocation_list or mcable_list.
+ util::visit(
+ [&](auto& metadata) { sc.sampler(sc.probe_id, sc.tag, &metadata, n_sample, sample_records.data()); },
+ p.metadata);
+}
+
+void run_samples(
+ const fvm_probe_weighted_multi& p,
+ const sampler_call_info& sc,
+ const fvm_value_type* raw_times,
+ const fvm_value_type* raw_samples,
+ std::vector<sample_record>& sample_records,
+ fvm_probe_scratch& scratch)
+{
+ const sample_size_type n_raw_per_sample = p.raw_handles.size();
+ sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample;
+ arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample);
+ arb_assert((unsigned)n_raw_per_sample==p.weight.size());
+
+ auto& sample_ranges = std::get<std::vector<cable_sample_range>>(scratch);
+ sample_ranges.clear();
+ sample_records.clear();
+
+ auto& tmp = std::get<std::vector<double>>(scratch);
+ tmp.clear();
+ tmp.reserve(n_raw_per_sample*n_sample);
+
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ for (sample_size_type i = 0; i<n_raw_per_sample; ++i) {
+ tmp.push_back(raw_samples[offset+i]*p.weight[i]);
+ }
+ }
+
+ const double* tmp_ptr = tmp.data();
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ sample_ranges.push_back({tmp_ptr, tmp_ptr+n_raw_per_sample});
+ tmp_ptr += n_raw_per_sample;
+ }
+
+ const auto& csample_ranges = sample_ranges;
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ sample_records.push_back(sample_record{time_type(raw_times[offset]), &csample_ranges[j]});
+ }
+
+ // (Unlike fvm_probe_multi, we only have mcable_list metadata.)
+ sc.sampler(sc.probe_id, sc.tag, &p.metadata, n_sample, sample_records.data());
+}
+void run_samples(
+ const fvm_probe_membrane_currents& p,
+ const sampler_call_info& sc,
+ const fvm_value_type* raw_times,
+ const fvm_value_type* raw_samples,
+ std::vector<sample_record>& sample_records,
+ fvm_probe_scratch& scratch)
+{
+ const sample_size_type n_raw_per_sample = p.raw_handles.size();
+ sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample;
+ arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample);
+
+ const auto n_cable = p.metadata.size();
+ const auto n_cv = p.cv_parent_cond.size();
+ const auto cables_by_cv = util::partition_view(p.cv_cables_divs);
+
+ auto& sample_ranges = std::get<std::vector<cable_sample_range>>(scratch);
+ sample_ranges.clear();
+
+ auto& tmp = std::get<std::vector<double>>(scratch);
+ tmp.assign(n_cable*n_sample, 0.);
+
+ sample_records.clear();
+
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ auto tmp_base = tmp.data()+j*n_cable;
+
+ // Each CV voltage contributes to the current sum of its parent's cables
+ // and its own cables.
+
+ const double* v = raw_samples+offset;
+ for (auto cv: util::make_span(n_cv)) {
+ fvm_index_type parent_cv = p.cv_parent[cv];
+ if (parent_cv+1==0) continue;
+
+ double cond = p.cv_parent_cond[cv];
+
+ double cv_I = v[cv]*cond;
+ double parent_cv_I = v[parent_cv]*cond;
+
+ for (auto cable_i: util::make_span(cables_by_cv[cv])) {
+ tmp_base[cable_i] -= (cv_I-parent_cv_I)*p.weight[cable_i];
+ }
+
+ for (auto cable_i: util::make_span(cables_by_cv[parent_cv])) {
+ tmp_base[cable_i] += (cv_I-parent_cv_I)*p.weight[cable_i];
+ }
+ }
+
+ sample_ranges.push_back({tmp_base, tmp_base+n_cable});
+ }
+
+ const auto& csample_ranges = sample_ranges;
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ sample_records.push_back(sample_record{time_type(raw_times[offset]), &csample_ranges[j]});
+ }
+
+ sc.sampler(sc.probe_id, sc.tag, &p.metadata, n_sample, sample_records.data());
+}
+
+// Generic run_samples dispatches on probe info variant type.
+void run_samples(
+ const fvm_probe_info& p,
+ const sampler_call_info& sc,
+ const fvm_value_type* raw_times,
+ const fvm_value_type* raw_samples,
+ std::vector<sample_record>& sample_records,
+ fvm_probe_scratch& scratch)
+{
+ util::visit([&](auto& x) {run_samples(x, sc, raw_times, raw_samples, sample_records, scratch); }, p.info);
+}
+
void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& event_lanes) {
time_type tstart = lowered_->time();
@@ -133,8 +373,8 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
//
// For each (schedule, sampler, probe set) in the sampler association
// map that will be triggered in this integration interval, create
- // sample events for the lowered cell, one for each scheduled sample
- // time and probe in the probe set.
+ // sample events for the lowered cell, one or more for each scheduled
+ // sample time and probe in the probe set.
//
// Each event is associated with an offset into the sample data and
// time buffers; these are assigned contiguously such that one call to
@@ -142,16 +382,6 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
// value as defined below, grouping together all the samples of the
// same probe for this callback in this association.
- struct sampler_call_info {
- sampler_function sampler;
- cell_member_type probe_id;
- probe_tag tag;
-
- // Offsets are into lowered cell sample time and event arrays.
- sample_size_type begin_offset;
- sample_size_type end_offset;
- };
-
PE(advance_samplesetup);
std::vector<sampler_call_info> call_info;
@@ -171,14 +401,18 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
for (cell_member_type pid: sa.probe_ids) {
auto cell_index = gid_index_map_.at(pid.gid);
auto p = probe_map_[pid];
+ sample_size_type n_raw = p.handle.n_raw();
- call_info.push_back({sa.sampler, pid, p.tag, n_samples, n_samples+n_times});
+ call_info.push_back({sa.sampler, pid, p.tag, n_samples, n_samples+n_times*n_raw});
for (auto t: sample_times) {
- sample_event ev{t, (cell_gid_type)cell_to_intdom_[cell_index], {p.handle, n_samples++}};
- sample_events.push_back(ev);
+ for (probe_handle h: p.handle.raw_handle_range()) {
+ sample_event ev{t, (cell_gid_type)cell_to_intdom_[cell_index], {h, n_samples++}};
+ sample_events.push_back(ev);
+ }
}
}
+ arb_assert(n_samples==call_info.back().end_offset);
}
// Sample events must be ordered by time for the lowered cell.
@@ -197,13 +431,12 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
std::vector<sample_record> sample_records;
sample_records.reserve(max_samples_per_call);
- for (auto& sc: call_info) {
- sample_records.clear();
- for (auto i = sc.begin_offset; i!=sc.end_offset; ++i) {
- sample_records.push_back(sample_record{time_type(result.sample_time[i]), &result.sample_value[i]});
- }
+ fvm_probe_scratch scratch;
+ reserve_scratch(scratch, max_samples_per_call);
- sc.sampler(sc.probe_id, sc.tag, sc.end_offset-sc.begin_offset, sample_records.data());
+ for (auto& sc: call_info) {
+ run_samples(probe_map_[sc.probe_id].handle, sc, result.sample_time.data(), result.sample_value.data(),
+ sample_records, scratch);
}
PL();
diff --git a/arbor/mc_cell_group.hpp b/arbor/mc_cell_group.hpp
index 70f255dc..d11da380 100644
--- a/arbor/mc_cell_group.hpp
+++ b/arbor/mc_cell_group.hpp
@@ -89,7 +89,7 @@ private:
std::vector<target_handle> target_handles_;
// Maps probe ids to probe handles (from lowered cell) and tags (from probe descriptions).
- probe_association_map<probe_handle> probe_map_;
+ probe_association_map<fvm_probe_info> probe_map_;
// Collection of samplers to be run against probes in this group.
sampler_association_map sampler_map_;
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
index ae3d0f21..f20bfb27 100644
--- a/arbor/morph/embed_pwlin.cpp
+++ b/arbor/morph/embed_pwlin.cpp
@@ -116,11 +116,11 @@ mcable_list data_cmp(const branch_pw_ratpoly<1, 0>& f, unsigned bid, double val,
}
struct embed_pwlin_data {
- branch_pw_ratpoly<1, 0> length;
- branch_pw_ratpoly<1, 0> directed_projection;
- branch_pw_ratpoly<1, 0> radius;
- branch_pw_ratpoly<2, 0> area;
- branch_pw_ratpoly<1, 1> ixa;
+ branch_pw_ratpoly<1, 0> length; // [µm]
+ branch_pw_ratpoly<1, 0> directed_projection; // [µm]
+ branch_pw_ratpoly<1, 0> radius; // [µm]
+ branch_pw_ratpoly<2, 0> area; // [µm²]
+ branch_pw_ratpoly<1, 1> ixa; // [1/µm]
explicit embed_pwlin_data(msize_t n_branch):
length(n_branch),
diff --git a/doc/cpp_cable_cell.rst b/doc/cpp_cable_cell.rst
index 56c42a42..8418546d 100644
--- a/doc/cpp_cable_cell.rst
+++ b/doc/cpp_cable_cell.rst
@@ -256,3 +256,306 @@ constants.
gprop.catalogue = &mycat;
gprop.default_parameters.reversal_potential_method["ca"] = "nernst1998/ca";
+
+Cable cell probes
+-----------------
+
+Various properties of a a cable cell can be sampled via one of the cable cell
+specific probe address described below. They fall into two classes: scalar
+probes are associated with a single real value, such as a membrane voltage
+or mechanism state value at a particular location; vector probes return
+multiple values corresponding to a quantity sampled over a whole cell.
+
+The sample data associated with a cable cell probe will either be a ``double``
+for scalar probes, or a ``cable_sample_range`` describing a half-open range
+of ``double`` values:
+
+.. code::
+
+ using cable_sample_range = std::pair<const double*, const double*>
+
+The probe metadata passed to the sampler will be a const pointer to:
+
+* ``mlocation`` for most scalar probes;
+
+* ``cable_probe_point_info`` for point mechanism state queries;
+
+* ``mcable_list`` for most vector queries;
+
+* ``std::vector<cable_probe_point_info>`` for cell-wide point mechanism state queries.
+
+The type ``cable_probe_point_info`` holds metadata for a single target on a cell:
+
+.. code::
+
+ struct cable_probe_point_info {
+ // Target number of point process instance on cell.
+ cell_lid_type target;
+
+ // Number of combined instances at this site.
+ unsigned multiplicity;
+
+ // Point on cell morphology where instance is placed.
+ mlocation loc;
+ };
+
+Note that the ``multiplicity`` will always be 1 if synapse coalescing is
+disabled.
+
+Cable cell probes that contingently do not correspond to a valid measurable
+quantity are ignored: samplers attached to them will receive no values.
+Mechanism state queries however will throw a ``cable_cell_error`` exception
+at simulation initialization if the requested state variable does not exist
+on the mechanism.
+
+Membrane voltage
+^^^^^^^^^^^^^^^^
+
+.. code::
+
+ struct cable_probe_membrane_voltage {
+ mlocation location;
+ };
+
+Queries cell membrane potential at the specified location.
+
+* Sample value: ``double``. Membrane potential in millivolts.
+
+* Metadata: ``mlocation``. Location as given in the probe address.
+
+
+.. code::
+
+ struct cable_probe_membrane_voltage_cell {};
+
+Queries cell membrane potential across whole cell.
+
+* Sample value: ``cable_sample_range``. Each value is the
+ average membrane potential in millivolts across an unbranched
+ component of the cell, as determined by the discretization.
+
+* Metadata: ``mcable_list``. Each cable in the cable list describes
+ the unbranched component for the corresponding sample value.
+
+Axial current
+^^^^^^^^^^^^^
+
+.. code::
+
+ struct cable_probe_axial_current {
+ mlocation location;
+ };
+
+Estimate intracellular current at given location in the distal direction.
+
+* Sample value: ``double``. Current in nanoamperes.
+
+* Metadata: ``mlocation``. Location as given in the probe address.
+
+
+Transmembrane current
+^^^^^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ struct cable_probe_ion_current_density {
+ mlocation location;
+ std::string ion;
+ };
+
+Membrance current density attributed to a particular ion at a given location.
+
+* Sample value: ``double``. Current density in amperes per square metre.
+
+* Metadata: ``mlocation``. Location as given in the probe address.
+
+
+.. code::
+
+ struct cable_probe_ion_current_cell {
+ std::string ion;
+ };
+
+Membrane current attributed to a particular ion across components of the cell.
+
+* Sample value: ``cable_sample_range``. Each value is the current in
+ nanoamperes across an unbranched component of the cell, as determined
+ by the discretization.
+
+* Metadata: ``mcable_list``. Each cable in the cable list describes
+ the unbranched component for the corresponding sample value.
+
+
+.. code::
+
+ struct cable_probe_total_ion_current_density {
+ mlocation location;
+ };
+
+Membrane current density at gvien location _excluding_ capacitive currents.
+
+* Sample value: ``double``. Current density in amperes per square metre.
+
+* Metadata: ``mlocation``. Location as given in the probe address.
+
+
+.. code::
+
+ struct cable_probe_total_ion_current_cell {};
+
+Membrane current _excluding_ capacitive currents across components of the cell.
+
+* Sample value: ``cable_sample_range``. Each value is the current in
+ nanoamperes across an unbranched component of the cell, as determined
+ by the discretization.
+
+* Metadata: ``mcable_list``. Each cable in the cable list describes
+ the unbranched component for the corresponding sample value.
+
+
+.. code::
+
+ struct cable_probe_total_current_cell {};
+
+Total membrance current across components of the cell.
+
+* Sample value: ``cable_sample_range``. Each value is the current in
+ nanoamperes across an unbranched component of the cell, as determined
+ by the discretization.
+
+* Metadata: ``mcable_list``. Each cable in the cable list describes
+ the unbranched component for the corresponding sample value.
+
+
+Ion concentration
+^^^^^^^^^^^^^^^^^
+
+.. code::
+
+ struct cable_probe_ion_int_concentration {
+ mlocation location;
+ std::string ion;
+ };
+
+Ionic internal concentration of ion at given location.
+
+* Sample value: ``double``. Ion concentration in millimoles per litre.
+
+* Metadata: ``mlocation``. Location as given in the probe address.
+
+
+.. code::
+
+ struct cable_probe_ion_int_concentration_cell {
+ std::string ion;
+ };
+
+Ionic external concentration of ion across components of the cell.
+
+* Sample value: ``cable_sample_range``. Each value is the concentration in
+ millimoles per lire across an unbranched component of the cell, as determined
+ by the discretization.
+
+* Metadata: ``mcable_list``. Each cable in the cable list describes
+ the unbranched component for the corresponding sample value.
+
+
+.. code::
+
+ struct cable_probe_ion_ext_concentration {
+ mlocation location;
+ std::string ion;
+ };
+
+Ionic internal concentration of ion at given location.
+
+* Sample value: ``double``. Ion concentration in millimoles per litre.
+
+* Metadata: ``mlocation``. Location as given in the probe address.
+
+
+.. code::
+
+ struct cable_probe_ion_ext_concentration_cell {
+ std::string ion;
+ };
+
+Ionic external concentration of ion across components of the cell.
+
+* Sample value: ``cable_sample_range``. Each value is the concentration in
+ millimoles per lire across an unbranched component of the cell, as determined
+ by the discretization.
+
+* Metadata: ``mcable_list``. Each cable in the cable list describes
+ the unbranched component for the corresponding sample value.
+
+
+
+Mechanism state
+^^^^^^^^^^^^^^^
+
+.. code::
+
+ struct cable_probe_density_state {
+ mlocation location;
+ std::string mechanism;
+ std::string state;
+ };
+
+
+Value of state variable in a density mechanism in at given location.
+If the mechanism is not defined at the location, the probe is ignored.
+
+* Sample value: ``double``. State variable value.
+
+* Metadata: ``mlocation``. Location as given in the probe address.
+
+
+.. code::
+
+ struct cable_probe_density_state_cell {
+ std::string mechanism;
+ std::string state;
+ };
+
+Value of state variable in adensity mechanism across components of the cell.
+
+* Sample value: ``cable_sample_range``. State variable values from the
+ mechanism across unbranched components of the cell, as determined
+ by the discretization and mechanism extent.
+
+* Metadata: ``mcable_list``. Each cable in the cable list describes
+ the unbranched component for the corresponding sample value.
+
+
+.. code::
+
+ struct cable_probe_point_state {
+ cell_lid_type target;
+ std::string mechanism;
+ std::string state;
+ };
+
+Value of state variable in a point mechanism associated with the given target.
+If the mechanism is not associated with this target, the probe is ignored.
+
+* Sample value: ``double``. State variable value.
+
+* Metadata: ``cable_probe_point_info``. Target number, multiplicity and location.
+
+
+.. code::
+
+ struct cable_probe_point_state_cell {
+ std::string mechanism;
+ std::string state;
+ };
+
+Value of state variable in a point mechanism for each of the targets in the cell
+with which it is associated.
+
+* Sample value: ``cable_sample_range``. State variable values at each associated
+ target.
+
+* Metadata: ``std::vector<cable_probe_point_info>``. Target metadata for each
+ associated target.
diff --git a/doc/sampling_api.rst b/doc/sampling_api.rst
index 1c025511..da4a80b6 100644
--- a/doc/sampling_api.rst
+++ b/doc/sampling_api.rst
@@ -77,14 +77,19 @@ will be passed to a sampler function or function object:
.. code-block:: cpp
using sampler_function =
- std::function<void (cell_member_type, probe_tag, size_t, const sample_record*)>;
+ std::function<void (cell_member_type, probe_tag, any_ptr, size_t, const sample_record*)>;
-where the parameters are respectively the probe id, the tag, the number
-of samples and a pointer to the sequence of sample records.
+where the parameters are respectively the probe id, the tag, a typed
+pointer to probe metadata, the number of samples, and finally a pointer
+to the sequence of sample records.
The ``probe_tag`` is the key given in the ``probe_info`` returned by
the recipe.
+The ``any_ptr`` value points to const probe-specific metadata; the type
+of the metadata will depend upon the probe address specified in the
+``probe_info`` provided by the recipe.
+
One ``sample_record`` struct contains one sample of the probe data at a
given simulation time point:
@@ -117,7 +122,7 @@ function. A simple sampler implementation for ``double`` data might be:
explicit scalar_sample(sample_data& samples): samples(samples) {}
- void operator()(cell_member_type id, probe_tag, size_t n, const sample_record* records) {
+ void operator()(cell_member_type id, probe_tag, any_ptr, size_t n, const sample_record* records) {
for (size_t i=0; i<n; ++i) {
const auto& rec = records[i];
diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index aeed60d2..0c33de29 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -9,3 +9,4 @@ add_subdirectory(bench)
add_subdirectory(ring)
add_subdirectory(gap_junctions)
add_subdirectory(single)
+add_subdirectory(probe-demo)
diff --git a/example/dryrun/dryrun.cpp b/example/dryrun/dryrun.cpp
index c8354c26..e40891c7 100644
--- a/example/dryrun/dryrun.cpp
+++ b/example/dryrun/dryrun.cpp
@@ -3,13 +3,13 @@
*
*/
+#include <cassert>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <nlohmann/json.hpp>
-#include <arbor/assert_macro.hpp>
#include <arbor/common_types.hpp>
#include <arbor/context.hpp>
#include <arbor/cable_cell.hpp>
@@ -129,7 +129,7 @@ public:
// Measure at the soma.
arb::mlocation loc{0, 0.0};
- return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
+ return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}};
}
private:
@@ -164,7 +164,7 @@ int main(int argc, char** argv) {
}
}
#endif
- arb_assert(arb::num_ranks(ctx)==params.num_ranks);
+ assert(arb::num_ranks(ctx)==params.num_ranks);
#ifdef ARB_PROFILE_ENABLED
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index 4168aa76..4f3001fc 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -101,7 +101,7 @@ public:
arb::probe_info get_probe(cell_member_type id) const override {
// Measure at the soma.
arb::mlocation loc{0, 1.};
- return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
+ return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}};
}
arb::util::any get_global_properties(cell_kind k) const override {
diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp
index bc392b43..9b62067f 100644
--- a/example/generators/generators.cpp
+++ b/example/generators/generators.cpp
@@ -6,6 +6,7 @@
* event generators, one inhibitory, and one excitatory, are attached.
*/
+#include <cassert>
#include <fstream>
#include <iomanip>
#include <iostream>
@@ -65,7 +66,7 @@ public:
}
cell_kind get_cell_kind(cell_gid_type gid) const override {
- arb_assert(gid==0); // There is only one cell in the model
+ assert(gid==0); // There is only one cell in the model
return cell_kind::cable;
}
@@ -77,13 +78,13 @@ public:
// The cell has one target synapse, which receives both inhibitory and exchitatory inputs.
cell_size_type num_targets(cell_gid_type gid) const override {
- arb_assert(gid==0); // There is only one cell in the model
+ assert(gid==0); // There is only one cell in the model
return 1;
}
// Return two generators attached to the one cell.
std::vector<arb::event_generator> event_generators(cell_gid_type gid) const override {
- arb_assert(gid==0); // There is only one cell in the model
+ assert(gid==0); // There is only one cell in the model
using RNG = std::mt19937_64;
@@ -116,17 +117,17 @@ public:
// There is one probe (for measuring voltage at the soma) on the cell
cell_size_type num_probes(cell_gid_type gid) const override {
- arb_assert(gid==0); // There is only one cell in the model
+ assert(gid==0); // There is only one cell in the model
return 1;
}
arb::probe_info get_probe(cell_member_type id) const override {
- arb_assert(id.gid==0); // There is one cell,
- arb_assert(id.index==0); // with one probe.
+ assert(id.gid==0); // There is one cell,
+ assert(id.index==0); // with one probe.
// Measure at the soma
arb::mlocation loc{0, 0.0};
- return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
+ return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}};
}
};
diff --git a/example/generators/readme.md b/example/generators/readme.md
index d8858740..fd32e77d 100644
--- a/example/generators/readme.md
+++ b/example/generators/readme.md
@@ -153,15 +153,15 @@ In our case, the cell has one probe, which refers to the voltage at the soma.
}
arb::probe_info get_probe(cell_member_type id) const override {
- // Get the appropriate kind for measuring voltage
- cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
-
// The location at which we measure: position 0 on branch 0.
// The cell has only one branch, branch 0, which is the soma.
- arb::mlocation loc(0, 0.0);
+ arb::mlocation loc{0, 0.0};
+
+ // The thing we are measuring is the membrane potential.
+ arb::cable_probe_membrane_voltage address{loc};
- // Put this together into a `probe_info`
- return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
+ // Put this together into a `probe_info`; the tag value 0 is not used.
+ return arb::probe_info{id, 0, address};
}
```
diff --git a/example/probe-demo/CMakeLists.txt b/example/probe-demo/CMakeLists.txt
new file mode 100644
index 00000000..e7cf7ad7
--- /dev/null
+++ b/example/probe-demo/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_executable(probe-demo EXCLUDE_FROM_ALL probe-demo.cpp)
+add_dependencies(examples probe-demo)
+target_link_libraries(probe-demo PRIVATE arbor ext-tinyopt)
diff --git a/example/probe-demo/probe-demo.cpp b/example/probe-demo/probe-demo.cpp
new file mode 100644
index 00000000..00bafdc8
--- /dev/null
+++ b/example/probe-demo/probe-demo.cpp
@@ -0,0 +1,247 @@
+#include <functional>
+#include <iomanip>
+#include <iostream>
+#include <stdexcept>
+#include <string>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include <arbor/load_balance.hpp>
+#include <arbor/cable_cell.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/region.hpp>
+#include <arbor/simulation.hpp>
+#include <arbor/sampling.hpp>
+#include <arbor/util/any.hpp>
+#include <arbor/util/any_ptr.hpp>
+#include <tinyopt/smolopt.h>
+
+// Simulate a cell modelled as a simple cable with HH dynamics,
+// emitting the results of a user specified probe over time.
+
+using arb::util::any;
+using arb::util::any_cast;
+using arb::util::any_ptr;
+
+const char* help_msg =
+ "[OPTION]... PROBE\n"
+ "\n"
+ " --dt=TIME set simulation dt to TIME [ms]\n"
+ " --until=TIME simulate until TIME [ms]\n"
+ " -n, --n-cv=N discretize with N CVs\n"
+ " -t, --sample=TIME take a sample every TIME [ms]\n"
+ " -x, --at=X take sample at relative position X along cable\n"
+ " -h, --help print extended usage information and exit\n"
+ "\n"
+ "Simulate a simple 1 mm cable with HH dynamics, taking samples according\n"
+ "to PROBE (see below). Unless otherwise specified, the simulation will\n"
+ "use 10 CVs and run for 100 ms with a time step of 0.025 ms.\n"
+ "\n"
+ "Samples are by default taken every 1 ms; for probes that test a specific\n"
+ "point on the cell, that point is 0.5 along the cable (i.e. 500 µm) unless\n"
+ "specified with the --at option.\n"
+ "\n"
+ "PROBE is one of:\n"
+ "\n"
+ " v membrane potential [mV] at X\n"
+ " i_axial axial (distal) current [nA] at X\n"
+ " j_ion total ionic membrane current density [A/m²] at X\n"
+ " j_na sodium ion membrane current density [A/m²] at X\n"
+ " j_k potassium ion membrane current density [A/m²] at X\n"
+ " c_na internal sodium concentration [mmol/L] at X\n"
+ " c_k internal potassium concentration [mmol/L] at X\n"
+ " hh_m HH state variable m at X\n"
+ " hh_h HH state variable h at X\n"
+ " hh_n HH state variable n at X\n"
+ "\n"
+ "where X is the relative position along the cable as described above, or else:\n"
+ "\n"
+ " all_v membrane potential [mV] in each CV\n"
+ " all_i_ion total ionic membrane current [nA] in each CV\n"
+ " all_i_na sodium ion membrane current [nA] in each CV\n"
+ " all_i_k potassium ion membrane current [nA] in each CV\n"
+ " all_i total membrane current [nA] in each CV\n"
+ " all_c_na internal sodium concentration [mmol/L] in each CV\n"
+ " all_c_k internal potassium concentration [mmol/L] in each CV\n"
+ " all_hh_m HH state variable m in each CV\n"
+ " all_hh_h HH state variable h in each CV\n"
+ " all_hh_n HH state variable n in each CV\n";
+
+struct options {
+ double sim_end = 100.0; // [ms]
+ double sim_dt = 0.025; // [ms]
+ double sample_dt = 1.0; // [ms]
+ unsigned n_cv = 10;
+
+ bool scalar_probe = true;
+ any probe_addr;
+ std::string value_name;
+};
+
+const char* argv0 = "";
+bool parse_options(options&, int& argc, char** argv);
+void vector_sampler(arb::cell_member_type, arb::probe_tag, any_ptr, std::size_t, const arb::sample_record*);
+void scalar_sampler(arb::cell_member_type, arb::probe_tag, any_ptr, std::size_t, const arb::sample_record*);
+
+struct cable_recipe: public arb::recipe {
+ arb::cable_cell_global_properties gprop;
+ any probe_addr;
+
+ explicit cable_recipe(any probe_addr, unsigned n_cv):
+ probe_addr(std::move(probe_addr))
+ {
+ gprop.default_parameters = arb::neuron_parameter_defaults;
+ gprop.default_parameters.discretization = arb::cv_policy_fixed_per_branch(n_cv);
+ }
+
+ arb::cell_size_type num_cells() const override { return 1; }
+ arb::cell_size_type num_probes(arb::cell_gid_type) const override { return 1; }
+ arb::cell_size_type num_targets(arb::cell_gid_type) const override { return 0; }
+
+ arb::probe_info get_probe(arb::cell_member_type probe_id) const override {
+ return {probe_id, 0, probe_addr}; // (tag value 0 is not used)
+ }
+
+ arb::cell_kind get_cell_kind(arb::cell_gid_type) const override {
+ return arb::cell_kind::cable;
+ }
+
+ arb::util::any get_global_properties(arb::cell_kind) const override {
+ return gprop;
+ }
+
+ arb::util::unique_any get_cell_description(arb::cell_gid_type) const override {
+ using namespace arb;
+ const double length = 1000; // [µm]
+ const double diam = 1; // [µm]
+
+ sample_tree samples({msample{0, 0, 0, 0.5*diam}, msample{length, 0, 0, 0.5*diam}}, {mnpos, 0u});
+ cable_cell c(samples);
+
+ c.paint(reg::all(), "hh"); // HH mechanism over whole cell.
+ c.place(mlocation{0, 0.}, i_clamp{0., INFINITY, 1.}); // Inject a 1 nA current indefinitely.
+ return c;
+ }
+};
+
+int main(int argc, char** argv) {
+ argv0 = argv[0];
+ try {
+ options opt;
+ if (!parse_options(opt, argc, argv)) {
+ return 0;
+ }
+
+ cable_recipe R(opt.probe_addr, opt.n_cv);
+
+ auto context = arb::make_context();
+ arb::simulation sim(R, arb::partition_load_balance(R, context), context);
+
+ sim.add_sampler(arb::all_probes, arb::regular_schedule(opt.sample_dt),
+ opt.scalar_probe? scalar_sampler: vector_sampler);
+
+ // CSV header for sample output:
+ std::cout << "t, " << (opt.scalar_probe? "x, ": "x0, x1, ") << opt.value_name << '\n';
+
+ sim.run(opt.sim_end, opt.sim_dt);
+ }
+ catch (to::option_error& e) {
+ to::usage(argv0, "[OPTIONS]... PROBE\nTry '--help' for more information.", e.what());
+ return 1;
+ }
+ catch (std::exception& e) {
+ std::cerr << "caught exception: " << e.what() << "\n";
+ return 2;
+ }
+}
+
+void scalar_sampler(arb::cell_member_type, arb::probe_tag, any_ptr meta, std::size_t n, const arb::sample_record* samples) {
+ auto* loc = any_cast<const arb::mlocation*>(meta);
+ assert(loc);
+
+ std::cout << std::fixed << std::setprecision(4);
+ for (std::size_t i = 0; i<n; ++i) {
+ auto* value = any_cast<const double*>(samples[i].data);
+ assert(value);
+
+ std::cout << samples[i].time << ", " << loc->pos << ", " << *value << '\n';
+ }
+}
+
+void vector_sampler(arb::cell_member_type, arb::probe_tag, any_ptr meta, std::size_t n, const arb::sample_record* samples) {
+ auto* cables_ptr = any_cast<const arb::mcable_list*>(meta);
+ assert(cables_ptr);
+ unsigned n_cable = cables_ptr->size();
+
+ std::cout << std::fixed << std::setprecision(4);
+ for (std::size_t i = 0; i<n; ++i) {
+ auto* value_range = any_cast<const arb::cable_sample_range*>(samples[i].data);
+ assert(value_range);
+ assert(n_cable==value_range->second-value_range->first);
+
+ for (unsigned j = 0; j<n_cable; ++j) {
+ arb::mcable where = (*cables_ptr)[j];
+ std::cout << samples[i].time << ", "
+ << where.prox_pos << ", "
+ << where.dist_pos << ", "
+ << value_range->first[j] << '\n';
+ }
+ }
+}
+
+bool parse_options(options& opt, int& argc, char** argv) {
+ using std::get;
+ using namespace to;
+
+ auto do_help = [&]() { usage(argv0, help_msg); };
+
+ // Map probe argument to output variable name, scalarity, and a lambda that makes specific probe address from a location.
+ std::pair<const char*, std::tuple<const char*, bool, std::function<any (double)>>> probe_tbl[] {
+ // located probes
+ {"v", {"v", true, [](double x) { return arb::cable_probe_membrane_voltage{{0, x}}; }}},
+ {"i_axial", {"i_axial", true, [](double x) { return arb::cable_probe_axial_current{{0, x}}; }}},
+ {"j_ion", {"j_ion", true, [](double x) { return arb::cable_probe_total_ion_current_density{{0, x}}; }}},
+ {"j_na", {"j_na", true, [](double x) { return arb::cable_probe_ion_current_density{{0, x}, "na"}; }}},
+ {"j_k", {"j_k", true, [](double x) { return arb::cable_probe_ion_current_density{{0, x}, "k"}; }}},
+ {"c_na", {"c_na", true, [](double x) { return arb::cable_probe_ion_int_concentration{{0, x}, "na"}; }}},
+ {"c_k", {"c_k", true, [](double x) { return arb::cable_probe_ion_int_concentration{{0, x}, "k"}; }}},
+ {"hh_m", {"hh_m", true, [](double x) { return arb::cable_probe_density_state{{0, x}, "hh", "m"}; }}},
+ {"hh_h", {"hh_h", true, [](double x) { return arb::cable_probe_density_state{{0, x}, "hh", "h"}; }}},
+ {"hh_n", {"hh_n", true, [](double x) { return arb::cable_probe_density_state{{0, x}, "hh", "n"}; }}},
+ // all-of-cell probes
+ {"all_v", {"v", false, [](double) { return arb::cable_probe_membrane_voltage_cell{}; }}},
+ {"all_i_ion", {"i_ion", false, [](double) { return arb::cable_probe_total_ion_current_cell{}; }}},
+ {"all_i_na", {"i_na", false, [](double) { return arb::cable_probe_ion_current_cell{"na"}; }}},
+ {"all_i_k", {"i_k", false, [](double) { return arb::cable_probe_ion_current_cell{"k"}; }}},
+ {"all_i", {"i", false, [](double) { return arb::cable_probe_total_current_cell{}; }}},
+ {"all_c_na", {"c_na", false, [](double) { return arb::cable_probe_ion_int_concentration_cell{"na"}; }}},
+ {"all_c_k", {"c_k", false, [](double) { return arb::cable_probe_ion_int_concentration_cell{"k"}; }}},
+ {"all_hh_m", {"hh_m", false, [](double) { return arb::cable_probe_density_state_cell{"hh", "m"}; }}},
+ {"all_hh_h", {"hh_h", false, [](double) { return arb::cable_probe_density_state_cell{"hh", "h"}; }}},
+ {"all_hh_n", {"hh_n", false, [](double) { return arb::cable_probe_density_state_cell{"hh", "n"}; }}}
+ };
+
+ std::tuple<const char*, bool, std::function<any (double)>> probe_spec;
+ double probe_pos = 0.5;
+
+ to::option cli_opts[] = {
+ { to::action(do_help), to::flag, to::exit, "-h", "--help" },
+ { {probe_spec, to::keywords(probe_tbl)}, to::single },
+ { opt.sim_dt, "--dt" },
+ { opt.sim_end, "--until" },
+ { opt.sample_dt, "-t", "--sample" },
+ { probe_pos, "-x", "--at" },
+ { opt.n_cv, "-n", "--n-cv" }
+ };
+
+ if (!to::run(cli_opts, argc, argv+1)) return false;
+ if (!get<2>(probe_spec)) throw to::user_option_error("missing PROBE");
+ if (argv[1]) throw to::user_option_error("unrecognized option");
+
+ opt.value_name = get<0>(probe_spec);
+ opt.scalar_probe = get<1>(probe_spec);
+ opt.probe_addr = get<2>(probe_spec)(probe_pos);
+ return true;
+}
+
diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp
index 561b09df..8f7acde0 100644
--- a/example/ring/ring.cpp
+++ b/example/ring/ring.cpp
@@ -116,7 +116,7 @@ public:
arb::probe_info get_probe(cell_member_type id) const override {
// Measure at the soma.
arb::mlocation loc{0, 0.0};
- return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
+ return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}};
}
arb::util::any get_global_properties(arb::cell_kind) const override {
diff --git a/example/single/single.cpp b/example/single/single.cpp
index 4e62819a..07d3cdac 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -39,7 +39,7 @@ struct single_recipe: public arb::recipe {
arb::probe_info get_probe(arb::cell_member_type probe_id) const override {
arb::mlocation mid_soma = {0, 0.5};
- arb::cell_probe_membrane_voltage probe = {mid_soma};
+ arb::cable_probe_membrane_voltage probe = {mid_soma};
// Probe info consists of: the probe id, a tag value to distinguish this probe
// from others for any attached sampler (unused), and the cell probe address.
diff --git a/python/recipe.cpp b/python/recipe.cpp
index 0ce022d9..1497307d 100644
--- a/python/recipe.cpp
+++ b/python/recipe.cpp
@@ -31,10 +31,10 @@ arb::util::unique_any py_recipe_shim::get_cell_description(arb::cell_gid_type gi
arb::probe_info cable_probe(std::string kind, arb::cell_member_type id, arb::mlocation loc) {
if (kind == "voltage") {
- return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
+ return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}};
}
else if (kind == "ionic current density") {
- return arb::probe_info{id, 0, arb::cell_probe_total_ionic_current_density{loc}};
+ return arb::probe_info{id, 0, arb::cable_probe_total_ion_current_density{loc}};
}
else throw pyarb_error(util::pprintf("unrecognized probe kind: {}", kind));
};
diff --git a/python/sampling.cpp b/python/sampling.cpp
index e281c2c2..d3a18b56 100644
--- a/python/sampling.cpp
+++ b/python/sampling.cpp
@@ -60,7 +60,7 @@ struct sample_callback {
sample_store(state)
{}
- void operator() (arb::cell_member_type probe_id, arb::probe_tag tag, std::size_t n, const arb::sample_record* recs) {
+ void operator() (arb::cell_member_type probe_id, arb::probe_tag tag, arb::util::any_ptr meta, std::size_t n, const arb::sample_record* recs) {
auto& v = sample_store->probe_buffer(probe_id);
for (std::size_t i = 0; i<n; ++i) {
if (auto p = arb::util::any_cast<const double*>(recs[i].data)) {
diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp
index ae69e3e2..bde7da44 100644
--- a/python/single_cell_model.cpp
+++ b/python/single_cell_model.cpp
@@ -41,7 +41,7 @@ struct trace_callback {
trace_callback(trace& t): trace_(t) {}
- void operator()(arb::cell_member_type probe_id, arb::probe_tag tag, std::size_t n, const arb::sample_record* recs) {
+ void operator()(arb::cell_member_type probe_id, arb::probe_tag tag, arb::util::any_ptr meta, std::size_t n, const arb::sample_record* recs) {
// Push each (time, value) pair from the last epoch into trace_.
for (std::size_t i=0; i<n; ++i) {
if (auto p = arb::util::any_cast<const double*>(recs[i].data)) {
@@ -117,7 +117,7 @@ struct single_cell_recipe: arb::recipe {
// For now only voltage can be selected for measurement.
const auto& loc = probes_[probe_id.index].site;
- return arb::probe_info{probe_id, 0, arb::cell_probe_membrane_voltage{loc}};
+ return arb::probe_info{probe_id, 0, arb::cable_probe_membrane_voltage{loc}};
}
// gap junctions
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 398c981d..bd48463e 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -1,8 +1,10 @@
# Build mechanisms used solely in unit tests.
set(test_mechanisms
+ ca_linear
celsius_test
diam_test
+ param_as_state
test_linear_state
test_linear_init
test_linear_init_shuffle
@@ -79,6 +81,7 @@ endforeach()
set(unit_sources
test_algorithms.cpp
test_any.cpp
+ test_any_visitor.cpp
test_backend.cpp
test_counter.cpp
test_cv_geom.cpp
diff --git a/test/unit/common_morphologies.hpp b/test/unit/common_morphologies.hpp
index 666fd0d8..a3b63bf4 100644
--- a/test/unit/common_morphologies.hpp
+++ b/test/unit/common_morphologies.hpp
@@ -26,13 +26,16 @@ static const arb::morphology m_sph_b1{arb::sample_tree(make_morph_samples(1), {a
// 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
+// spherical root, six branches:
+// branch 0 (spherical root) has child branches 1 and 2; branch 2 has child branches 3, 4 and 5.
static const arb::morphology m_sph_b6{arb::sample_tree(make_morph_samples(8), {arb::mnpos, 0u, 1u, 0u, 3u, 4u, 4u, 4u}), true};
// regular root, six branches
+// branch 0 has child branches 1 and 2; branch 2 has child branches 3, 4 and 5.
static const arb::morphology m_reg_b6{arb::sample_tree(make_morph_samples(7), {arb::mnpos, 0u, 1u, 1u, 2u, 2u, 2u}), false};
// regular root, six branches, mutiple top level branches.
+// branch 0 has child branches 1 and 2; branch 3 has child branches 4 and 5.
static const arb::morphology m_mlt_b6{arb::sample_tree(make_morph_samples(7), {arb::mnpos, 0u, 1u, 1u, 0u, 4u, 4u}), false};
static std::pair<const char*, arb::morphology> test_morphologies[] = {
diff --git a/test/unit/mod/ca_linear.mod b/test/unit/mod/ca_linear.mod
new file mode 100644
index 00000000..1cd854cd
--- /dev/null
+++ b/test/unit/mod/ca_linear.mod
@@ -0,0 +1,19 @@
+: Test density mechanism for passive calcium current.
+
+NEURON {
+ SUFFIX ca_linear
+ USEION ca WRITE ica VALENCE 2
+ RANGE g
+}
+
+PARAMETER {
+ g = .001 (S/cm2)
+}
+
+ASSIGNED {}
+
+INITIAL {}
+
+BREAKPOINT {
+ ica = g*v
+}
diff --git a/test/unit/mod/param_as_state.mod b/test/unit/mod/param_as_state.mod
new file mode 100644
index 00000000..463d6be0
--- /dev/null
+++ b/test/unit/mod/param_as_state.mod
@@ -0,0 +1,24 @@
+: Copy range parameter to state variable exactly.
+
+NEURON {
+ SUFFIX param_as_state
+ RANGE p
+}
+
+PARAMETER {
+ p = 1
+}
+
+ASSIGNED {}
+
+STATE {
+ s
+}
+
+INITIAL {
+ s = p
+}
+
+BREAKPOINT {
+}
+
diff --git a/test/unit/test_any_visitor.cpp b/test/unit/test_any_visitor.cpp
new file mode 100644
index 00000000..6c22ec07
--- /dev/null
+++ b/test/unit/test_any_visitor.cpp
@@ -0,0 +1,190 @@
+#include <string>
+#include <type_traits>
+#include <typeinfo>
+
+#include <arbor/util/any.hpp>
+#include <arbor/util/any_visitor.hpp>
+
+#include "../gtest.h"
+#include "common.hpp"
+
+using namespace std::string_literals;
+using arb::util::any;
+using arb::util::any_cast;
+using arb::util::any_visitor;
+using arb::util::bad_any_cast;
+using arb::util::overload;
+
+TEST(any_visitor, simple) {
+ enum { A0, B0, C0 };
+ struct A { int value = A0; };
+ struct B { int value = B0; };
+ struct C { int value = C0; };
+
+ using V = any_visitor<A, B, C>;
+
+ auto get_value = [](auto y) { return y.value; };
+
+ EXPECT_EQ(A0, V::visit(get_value, any(A{})));
+ EXPECT_EQ(B0, V::visit(get_value, any(B{})));
+ EXPECT_EQ(C0, V::visit(get_value, any(C{})));
+}
+
+TEST(any_visitor, heterogeneous) {
+ struct r_base { int value; };
+ struct r_derived: r_base { int other; };
+
+ struct D { r_base item; };
+ struct E { r_derived item; };
+
+ using DE_visitor = any_visitor<D, E>;
+ using E_visitor = any_visitor<E>;
+
+ auto get_item = [](auto y) { return y.item; };
+
+ // Return type is common type across possible returns of visiting functor.
+
+ {
+ auto result = E_visitor::visit(get_item, any(E{}));
+ EXPECT_TRUE((std::is_same<decltype(result), r_derived>::value));
+ }
+ {
+ auto result = DE_visitor::visit(get_item, any(E{}));
+ EXPECT_TRUE((std::is_same<decltype(result), r_base>::value));
+ }
+}
+
+TEST(any_visitor, unmatched) {
+ enum { A0, B0, C0, D0 };
+ struct A { int value = A0; };
+ struct B { int value = B0; };
+ struct C { int value = C0; };
+ struct D { int value = D0; };
+
+ using V = any_visitor<A, B, C>;
+
+ struct catch_unmatched {
+ bool operator()(A) const { return true; }
+ bool operator()(B) const { return true; }
+ bool operator()(C) const { return true; }
+ bool operator()() const { return false; }
+ } f;
+
+ EXPECT_EQ(true, V::visit(f, any(A{})));
+ EXPECT_EQ(true, V::visit(f, any(B{})));
+ EXPECT_EQ(true, V::visit(f, any(C{})));
+ EXPECT_EQ(false, V::visit(f, any(D{})));
+
+ auto get_value = [](auto y) { return y.value; };
+
+ EXPECT_NO_THROW(V::visit(get_value, any(A{})));
+ EXPECT_NO_THROW(V::visit(get_value, any(A{})));
+ EXPECT_NO_THROW(V::visit(get_value, any(A{})));
+ EXPECT_THROW(V::visit(get_value, any(D{})), bad_any_cast);
+}
+
+TEST(any_visitor, preserve_visiting_ref) {
+ struct check_ref {
+ check_ref(int& result): result(result) {}
+ int& result;
+
+ void operator()(...) & { result = 1; };
+ void operator()(...) && { result = 2; };
+ };
+
+ struct A {};
+ struct B {};
+
+ int result = 0;
+ check_ref lref(result);
+
+ any_visitor<A, B>::visit(lref, any(B{}));
+ EXPECT_EQ(1, result);
+
+ result = 0;
+ any_visitor<A, B>::visit(check_ref(result), any(B{}));
+ EXPECT_EQ(2, result);
+}
+
+TEST(any_visitor, preserve_any_qual) {
+ struct A { int x = 0; };
+ struct B { int x = 0; };
+
+ struct check_qual {
+ int operator()(const A&) { return 0; }
+ int operator()(A&) { return 1; }
+ int operator()(A&&) { return 2; }
+
+ int operator()(const B&) { return 3; }
+ int operator()(B&) { return 4; }
+ int operator()(B&&) { return 5; }
+ };
+
+ using V = any_visitor<A, B>;
+ check_qual q;
+
+ any a(A{});
+ const any const_a(A{});
+ EXPECT_EQ(0, V::visit(q, const_a));
+ EXPECT_EQ(1, V::visit(q, a));
+ EXPECT_EQ(2, V::visit(q, any(A{})));
+
+ any b(B{});
+ const any const_b(B{});
+ EXPECT_EQ(3, V::visit(q, const_b));
+ EXPECT_EQ(4, V::visit(q, b));
+ EXPECT_EQ(5, V::visit(q, any(B{})));
+
+ any v(B{10});
+ V::visit([](auto& u) { ++u.x; }, v);
+ EXPECT_EQ(11, any_cast<B>(v).x);
+}
+
+TEST(overload, simple) {
+ struct A { int value; };
+ struct B { double value; };
+
+ auto f = overload([](A a) { return a.value+10; },
+ [](B b) { return b.value+20; });
+
+ EXPECT_TRUE((std::is_same<decltype(f(A{})), int>::value));
+ EXPECT_TRUE((std::is_same<decltype(f(B{})), double>::value));
+
+ EXPECT_EQ(12, f(A{2}));
+ EXPECT_EQ(22.5, f(B{2.5}));
+}
+
+namespace {
+int f1(int) { return 1; }
+int f2(int) { return 2; }
+}
+
+TEST(overload, precedence) {
+ // Overloaded functional object will match the first invocable operator(), not the best match.
+
+ EXPECT_EQ(1, overload(f1, f2)(0));
+ EXPECT_EQ(2, overload(f2, f1)(0));
+
+ EXPECT_EQ(1, overload([](int) { return 1; }, [](double) { return 2; })(0));
+ EXPECT_EQ(1, overload([](int) { return 1; }, [](double) { return 2; })(2.3));
+}
+
+TEST(overload, qualified_match) {
+ struct A {};
+
+ auto f = overload(
+ [](A&&) { return 1; },
+ [](const A&&) { return 2; },
+ [](A&) { return 3; },
+ [](const A&) { return 4; });
+
+ A a;
+ const A const_a;
+
+ EXPECT_EQ(1, f(std::move(a)));
+ EXPECT_EQ(2, f(std::move(const_a)));
+ EXPECT_EQ(3, f(a));
+ EXPECT_EQ(4, f(const_a));
+}
+
+
diff --git a/test/unit/test_cv_layout.cpp b/test/unit/test_cv_layout.cpp
index 62989fd8..af731d0f 100644
--- a/test/unit/test_cv_layout.cpp
+++ b/test/unit/test_cv_layout.cpp
@@ -174,12 +174,12 @@ TEST(cv_layout, zero_size_cv) {
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.
+ // All non-conductance values for zero-size cv_x should be zero, or value of parent.
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]);
+ EXPECT_EQ(D.init_membrane_potential[cv_a], D.init_membrane_potential[cv_x]);
+ EXPECT_EQ(D.temperature_K[cv_a], D.temperature_K[cv_x]);
// Face conductance for zero-size cv_x:
double l_x = cell.embedding().branch_length(0);
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index 54403866..68314cc3 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -990,8 +990,8 @@ TEST(fvm_layout, iinterp) {
EXPECT_TRUE(I.proximal_cv>=D.geometry.cell_cv_interval(cell_idx).first);
double fc = D.face_conductance.at(I.distal_cv);
- EXPECT_DOUBLE_EQ(-fc, I.proximal_coef);
- EXPECT_DOUBLE_EQ(+fc, I.distal_coef);
+ EXPECT_DOUBLE_EQ(+fc, I.proximal_coef);
+ EXPECT_DOUBLE_EQ(-fc, I.distal_coef);
}
}
}
@@ -1039,8 +1039,8 @@ TEST(fvm_layout, iinterp) {
EXPECT_EQ(0, I.proximal_cv);
EXPECT_EQ(1, I.distal_cv);
- EXPECT_EQ(-fc1, I.proximal_coef);
- EXPECT_EQ(+fc1, I.distal_coef);
+ EXPECT_EQ(+fc1, I.proximal_coef);
+ EXPECT_EQ(-fc1, I.distal_coef);
}
// Branch 2:
@@ -1053,7 +1053,7 @@ TEST(fvm_layout, iinterp) {
EXPECT_EQ(0, I.proximal_cv);
EXPECT_EQ(2, I.distal_cv);
- EXPECT_EQ(-fc2, I.proximal_coef);
- EXPECT_EQ(+fc2, I.distal_coef);
+ EXPECT_EQ(+fc2, I.proximal_coef);
+ EXPECT_EQ(-fc2, I.distal_coef);
}
}
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 69a8e30d..c10c138f 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -225,7 +225,7 @@ TEST(fvm_lowered, matrix_init)
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell fvcell(context);
fvcell.initialize({0}, cable1d_recipe(cell), cell_to_intdom, targets, probe_map);
@@ -278,7 +278,7 @@ TEST(fvm_lowered, target_handles) {
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
auto test_target_handles = [&](fvm_cell& cell) {
mechanism *expsyn = find_mechanism(cell, "expsyn");
@@ -361,7 +361,7 @@ TEST(fvm_lowered, stimulus) {
const auto& A = D.cv_area;
std::vector<target_handle> targets;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell fvcell(context);
fvcell.initialize({0}, cable1d_recipe(cells), cell_to_intdom, targets, probe_map);
@@ -449,7 +449,7 @@ TEST(fvm_lowered, derived_mechs) {
cable1d_recipe rec(cells);
rec.catalogue().derive("custom_kin1", "test_kin1", {{"tau", 20.0}});
- cell_probe_total_ionic_current_density where{{1, 0.3}};
+ cable_probe_total_ion_current_density where{{1, 0.3}};
rec.add_probe(0, 0, where);
rec.add_probe(1, 0, where);
rec.add_probe(2, 0, where);
@@ -459,7 +459,7 @@ TEST(fvm_lowered, derived_mechs) {
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
arb::execution_context context(resources);
fvm_cell fvcell(context);
@@ -490,12 +490,13 @@ TEST(fvm_lowered, derived_mechs) {
std::vector<double> samples[3];
- sampler_function sampler = [&](cell_member_type pid, probe_tag, std::size_t n, const sample_record* records) {
- for (std::size_t i = 0; i<n; ++i) {
- double v = *util::any_cast<const double*>(records[i].data);
- samples[pid.gid].push_back(v);
- }
- };
+ sampler_function sampler =
+ [&](cell_member_type pid, probe_tag, util::any_ptr, std::size_t n, const sample_record* records) {
+ for (std::size_t i = 0; i<n; ++i) {
+ double v = *util::any_cast<const double*>(records[i].data);
+ samples[pid.gid].push_back(v);
+ }
+ };
float times[] = {10.f, 20.f};
@@ -530,7 +531,7 @@ TEST(fvm_lowered, read_valence) {
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
{
std::vector<cable_cell> cells(1);
@@ -684,7 +685,7 @@ TEST(fvm_lowered, ionic_currents) {
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell fvcell(context);
fvcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
@@ -729,7 +730,7 @@ TEST(fvm_lowered, point_ionic_current) {
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell fvcell(context);
fvcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
@@ -809,7 +810,7 @@ TEST(fvm_lowered, weighted_write_ion) {
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell fvcell(context);
fvcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index de0acb73..85731963 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -1,22 +1,35 @@
#include "../gtest.h"
+#include <cmath>
+#include <vector>
+
#include <arbor/cable_cell.hpp>
#include <arbor/common_types.hpp>
+#include <arbor/load_balance.hpp>
+#include <arbor/mechanism.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/mechinfo.hpp>
+#include <arbor/sampling.hpp>
+#include <arbor/schedule.hpp>
+#include <arbor/simple_sampler.hpp>
+#include <arbor/simulation.hpp>
#include <arbor/version.hpp>
#include <arborenv/gpu_env.hpp>
#include "backends/event.hpp"
#include "backends/multicore/fvm.hpp"
+#include "backends/multicore/mechanism.hpp"
#ifdef ARB_GPU_ENABLED
#include "backends/gpu/fvm.hpp"
+#include "backends/gpu/mechanism.hpp"
#endif
#include "fvm_lowered_cell_impl.hpp"
#include "memory/gpu_wrappers.hpp"
+#include "util/filter.hpp"
#include "util/rangeutil.hpp"
#include "common.hpp"
+#include "common_morphologies.hpp"
#include "unit_test_catalogue.hpp"
#include "../common_cells.hpp"
#include "../simple_recipes.hpp"
@@ -27,7 +40,6 @@ using multicore_fvm_cell = fvm_lowered_cell_impl<multicore::backend>;
using multicore_shared_state = multicore::backend::shared_state;
ACCESS_BIND(std::unique_ptr<multicore_shared_state> multicore_fvm_cell::*, multicore_fvm_state_ptr, &multicore_fvm_cell::state_);
-
template <typename Backend>
struct backend_access {
using fvm_cell = multicore_fvm_cell;
@@ -36,7 +48,8 @@ struct backend_access {
return *(cell.*multicore_fvm_state_ptr).get();
}
- static fvm_value_type deref(const fvm_value_type* p) { return *p; }
+ template <typename fvm_type>
+ static fvm_type deref(const fvm_type* p) { return *p; }
};
#ifdef ARB_GPU_ENABLED
@@ -53,8 +66,9 @@ struct backend_access<gpu::backend> {
return *(cell.*gpu_fvm_state_ptr).get();
}
- static fvm_value_type deref(const fvm_value_type* p) {
- fvm_value_type r;
+ template <typename fvm_type>
+ static fvm_type deref(const fvm_type* p) {
+ fvm_type r;
memory::gpu_memcpy_d2h(&r, p, sizeof(r));
return r;
}
@@ -62,12 +76,25 @@ struct backend_access<gpu::backend> {
#endif
+// Used in a number of tests below, produce a simple Y-shaped 3-branch cell morphology with
+// linearly tapered branches.
+
+static morphology make_y_morphology() {
+ return morphology(sample_tree(
+ {msample{{0., 0., 0., 1.}, 0},
+ msample{{100., 0., 0., 0.8}, 0},
+ msample{{100., 100., 0., 0.5}, 0},
+ msample{{100., 0, 100., 0.4}, 0}},
+ {mnpos, 0u, 1u, 1u}));
+}
+
template <typename Backend>
void run_v_i_probe_test(const context& ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const fvm_value_type* p) { return backend_access<Backend>::deref(p); };
cable_cell bs = make_cell_ball_and_stick(false);
+ bs.default_parameters.discretization = cv_policy_fixed_per_branch(1);
i_clamp stim(0, 100, 0.3);
bs.place(mlocation{1, 1}, stim);
@@ -78,13 +105,13 @@ void run_v_i_probe_test(const context& ctx) {
mlocation loc1{1, 1};
mlocation loc2{1, 0.3};
- rec.add_probe(0, 10, cell_probe_membrane_voltage{loc0});
- rec.add_probe(0, 20, cell_probe_membrane_voltage{loc1});
- rec.add_probe(0, 30, cell_probe_total_ionic_current_density{loc2});
+ rec.add_probe(0, 10, cable_probe_membrane_voltage{loc0});
+ rec.add_probe(0, 20, cable_probe_membrane_voltage{loc1});
+ rec.add_probe(0, 30, cable_probe_total_ion_current_density{loc2});
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell lcell(*ctx);
lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
@@ -96,45 +123,130 @@ void run_v_i_probe_test(const context& ctx) {
EXPECT_EQ(20, probe_map.at({0, 1}).tag);
EXPECT_EQ(30, probe_map.at({0, 2}).tag);
- probe_handle p0 = probe_map.at({0, 0}).handle;
- probe_handle p1 = probe_map.at({0, 1}).handle;
- probe_handle p2 = probe_map.at({0, 2}).handle;
+ auto get_probe_handle = [&](cell_member_type x, unsigned i = 0) {
+ return probe_map.at(x).handle.raw_handle_range()[i];
+ };
+
+ // Voltage probes are interpolated, so expect fvm_probe_info
+ // to wrap an fvm_probe_interpolated; ion current density is
+ // a scalar, so should wrap fvm_probe_scalar.
+
+ ASSERT_TRUE(util::get_if<fvm_probe_interpolated>(probe_map.at({0, 0}).handle.info));
+ ASSERT_TRUE(util::get_if<fvm_probe_interpolated>(probe_map.at({0, 1}).handle.info));
+ ASSERT_TRUE(util::get_if<fvm_probe_scalar>(probe_map.at({0, 2}).handle.info));
- // Expect initial probe values to be the resting potential
- // for the voltage probes (cell membrane potential should
- // be constant), and zero for the current probe.
+ probe_handle p0a = get_probe_handle({0, 0}, 0);
+ probe_handle p0b = get_probe_handle({0, 0}, 1);
+ probe_handle p1a = get_probe_handle({0, 1}, 0);
+ probe_handle p1b = get_probe_handle({0, 1}, 1);
+ probe_handle p2 = get_probe_handle({0, 2});
+
+ // Ball-and-stick cell with default discretization policy should
+ // have three CVs, one for branch 0, one trivial one covering the
+ // branch point, and one for branch 1.
+ //
+ // Consequently, expect the interpolated voltage probe handles
+ // to be on CVs 0 and 1 for probe 0,0 on branch 0,
+ // and on CVs 1 and 2 for probe 0,1 on branch 1.
auto& state = backend_access<Backend>::state(lcell);
auto& voltage = state.voltage;
+ EXPECT_EQ(voltage.data(), p0a);
+ EXPECT_EQ(voltage.data()+1, p0b);
+ EXPECT_EQ(voltage.data()+1, p1a);
+ EXPECT_EQ(voltage.data()+2, p1b);
+
+ // Expect initial raw probe handle values to be the resting potential for
+ // the voltage probes (cell membrane potential should be constant), and
+ // zero for the current probe.
+
fvm_value_type resting = voltage[0];
EXPECT_NE(0.0, resting);
- // (Probe handles are just pointers in this implementation).
- EXPECT_EQ(resting, deref(p0));
- EXPECT_EQ(resting, deref(p1));
+ EXPECT_EQ(resting, deref(p0a));
+ EXPECT_EQ(resting, deref(p0b));
+ EXPECT_EQ(resting, deref(p1a));
+ EXPECT_EQ(resting, deref(p1a));
EXPECT_EQ(0.0, deref(p2));
// After an integration step, expect voltage probe values
- // to differ from resting, and between each other, and
- // for there to be a non-zero current.
- //
- // First probe, at (0,0), should match voltage in first
- // compartment.
+ // to differ from resting, and for there to be a non-zero current.
lcell.integrate(0.01, 0.0025, {}, {});
- EXPECT_NE(resting, deref(p0));
- EXPECT_NE(resting, deref(p1));
- EXPECT_NE(deref(p0), deref(p1));
+ EXPECT_NE(resting, deref(p0a));
+ EXPECT_NE(resting, deref(p0b));
+ EXPECT_NE(resting, deref(p1a));
+ EXPECT_NE(resting, deref(p1b));
EXPECT_NE(0.0, deref(p2));
+}
+
+template <typename Backend>
+void run_v_cable_probe_test(const context& ctx) {
+ using fvm_cell = typename backend_access<Backend>::fvm_cell;
+
+ // Take the per-cable voltage over a Y-shaped cell with and without
+ // interior forks in the discretization. The metadata will be used
+ // to determine the corresponding CVs for each cable, and the raw
+ // pointer to backend data checked against the expected CV offset.
+
+ cable_cell cell(make_y_morphology());
+
+ std::pair<const char*, cv_policy> test_policies[] = {
+ {"trivial fork", cv_policy_fixed_per_branch(3, cv_policy_flag::none)},
+ {"interior fork", cv_policy_fixed_per_branch(3, cv_policy_flag::interior_forks)},
+ };
+
+ for (auto& testcase: test_policies) {
+ SCOPED_TRACE(testcase.first);
+ cell.default_parameters.discretization = testcase.second;
+
+ cable1d_recipe rec(cell, false);
+ rec.add_probe(0, 0, cable_probe_membrane_voltage_cell{});
+
+ std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
+ probe_association_map<fvm_probe_info> probe_map;
+
+ fvm_cell lcell(*ctx);
+ lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
+
+ ASSERT_EQ(1u, probe_map.size());
+
+ const fvm_probe_multi* h_ptr = util::get_if<fvm_probe_multi>(probe_map.at({0, 0}).handle.info);
+ ASSERT_TRUE(h_ptr);
+ auto& h = *h_ptr;
+
+ const mcable_list* cl_ptr = util::get_if<mcable_list>(h_ptr->metadata);
+ ASSERT_TRUE(cl_ptr);
+ auto& cl = *cl_ptr;
+
+ ASSERT_EQ(h.raw_handles.size(), cl.size());
- fvm_value_type v = voltage[0];
- EXPECT_EQ(v, deref(p0));
+ // Independetly discretize the cell so we can follow cable–CV relationship.
+
+ cv_geometry geom = cv_geometry_from_ends(cell, testcase.second.cv_boundary_points(cell));
+
+ // For each cable in metadata, get CV from geom and confirm raw handle is
+ // state voltage + CV.
+
+ auto& state = backend_access<Backend>::state(lcell);
+ auto& voltage = state.voltage;
+
+ for (auto i: util::count_along(*cl_ptr)) {
+ mlocation cable_mid{cl[i].branch, 0.5*(cl[i].prox_pos+cl[i].dist_pos)};
+ auto cv = geom.location_cv(0, cable_mid, cv_prefer::cv_empty);
+
+ EXPECT_EQ(voltage.data()+cv, h.raw_handles[i]);
+ }
+ }
}
template <typename Backend>
void run_expsyn_g_probe_test(const context& ctx, bool coalesce_synapses = false) {
+ SCOPED_TRACE(coalesce_synapses? "coalesced": "uncoalesced");
+
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const fvm_value_type* p) { return backend_access<Backend>::deref(p); };
@@ -151,24 +263,30 @@ void run_expsyn_g_probe_test(const context& ctx, bool coalesce_synapses = false)
bs.default_parameters.discretization = cv_policy_fixed_per_branch(2);
cable1d_recipe rec(bs, coalesce_synapses);
- rec.add_probe(0, 10, cell_probe_point_state{0u, "expsyn", "g"});
- rec.add_probe(0, 20, cell_probe_point_state{1u, "expsyn", "g"});
+ rec.add_probe(0, 10, cable_probe_point_state{0u, "expsyn", "g"});
+ rec.add_probe(0, 20, cable_probe_point_state{1u, "expsyn", "g"});
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell lcell(*ctx);
lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
EXPECT_EQ(2u, rec.num_probes(0));
EXPECT_EQ(2u, probe_map.size());
+ ASSERT_EQ(1u, probe_map.count({0, 0}));
+ ASSERT_EQ(1u, probe_map.count({0, 1}));
EXPECT_EQ(10, probe_map.at({0, 0}).tag);
EXPECT_EQ(20, probe_map.at({0, 1}).tag);
- probe_handle p0 = probe_map.at({0, 0}).handle;
- probe_handle p1 = probe_map.at({0, 1}).handle;
+ auto probe_scalar_handle = [&](cell_member_type x) {
+ return probe_map.at(x).handle.raw_handle_range()[0];
+ };
+
+ probe_handle p0 = probe_scalar_handle({0, 0});
+ probe_handle p1 = probe_scalar_handle({0, 1});
// Expect initial probe values to be intial synapse g == 0.
@@ -209,6 +327,127 @@ void run_expsyn_g_probe_test(const context& ctx, bool coalesce_synapses = false)
}
}
+template <typename Backend>
+void run_expsyn_g_cable_probe_test(const context& ctx, bool coalesce_synapses = false) {
+ SCOPED_TRACE(coalesce_synapses? "coalesced": "uncoalesced");
+
+ using fvm_cell = typename backend_access<Backend>::fvm_cell;
+ auto deref = [](const auto* p) { return backend_access<Backend>::deref(p); };
+
+ // Place a mixture of expsyn and exp2syn synapses over two cells.
+ // Confirm that a whole-cell expsyn g state probe gives:
+ // * A metadata element for each placed expsyn probe.
+ // * A multiplicity that matches the number of probes sharing a CV if synapses
+ // were coalesced.
+ // * A raw handle that sees an event sent to the corresponding target,
+
+ cv_policy policy = cv_policy_fixed_per_branch(3);
+
+ cable_cell cell(make_y_morphology());
+ cell.default_parameters.discretization = policy;
+
+ std::unordered_map<cell_lid_type, mlocation> expsyn_target_loc_map;
+
+ for (unsigned bid = 0; bid<3u; ++bid) {
+ for (unsigned j = 0; j<10; ++j) {
+ mlocation expsyn_loc{bid, 0.1*j};
+ lid_range target_lids = cell.place(expsyn_loc, "expsyn");
+
+ ASSERT_EQ(1u, target_lids.end-target_lids.begin);
+ expsyn_target_loc_map[target_lids.begin] = expsyn_loc;
+
+ cell.place(mlocation{bid, 0.1*j+0.05}, "exp2syn");
+ }
+ }
+ const unsigned n_expsyn = 30;
+
+ std::vector<cable_cell> cells(2, cell);
+ cable1d_recipe rec(cells, coalesce_synapses);
+
+ rec.add_probe(0, 0, cable_probe_point_state_cell{"expsyn", "g"});
+ rec.add_probe(1, 0, cable_probe_point_state_cell{"expsyn", "g"});
+
+ std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
+ probe_association_map<fvm_probe_info> probe_map;
+
+ fvm_cell lcell(*ctx);
+ lcell.initialize({0, 1}, rec, cell_to_intdom, targets, probe_map);
+
+ // Send an event to each expsyn synapse with a weight = target+100*cell_gid, and
+ // integrate for a tiny time step.
+
+ std::vector<deliverable_event> events;
+ for (unsigned i: {0u, 1u}) {
+ // Cells have the same number of targets, so the offset for cell 1 is exactly...
+ cell_local_size_type cell_offset = i==0? 0: targets.size()/2;
+
+ for (auto target_id: util::keys(expsyn_target_loc_map)) {
+ deliverable_event ev{0., targets.at(target_id+cell_offset), float(target_id+100*i)};
+ events.push_back(ev);
+ }
+ }
+ (void)lcell.integrate(1e-5, 1e-5, events, {});
+
+ // Independently get cv geometry to compute CV indices.
+
+ cv_geometry geom = cv_geometry_from_ends(cells[0], policy.cv_boundary_points(cells[0]));
+ append(geom, cv_geometry_from_ends(cells[1], policy.cv_boundary_points(cells[1])));
+
+ ASSERT_EQ(2u, probe_map.size());
+ for (unsigned i: {0u, 1u}) {
+ const auto* h_ptr = util::get_if<fvm_probe_multi>(probe_map.at({i, 0}).handle.info);
+ ASSERT_TRUE(h_ptr);
+
+ const auto* m_ptr = util::get_if<std::vector<cable_probe_point_info>>(h_ptr->metadata);
+ ASSERT_TRUE(m_ptr);
+
+ const fvm_probe_multi& h = *h_ptr;
+ const std::vector<cable_probe_point_info> m = *m_ptr;
+
+ ASSERT_EQ(h.raw_handles.size(), m.size());
+ ASSERT_EQ(n_expsyn, m.size());
+
+ std::vector<double> expected_coalesced_cv_value(geom.size());
+ std::vector<double> expected_uncoalesced_value(targets.size());
+
+ std::vector<double> target_cv(targets.size(), (unsigned)-1);
+ std::unordered_map<fvm_size_type, unsigned> cv_expsyn_count;
+
+ for (unsigned j = 0; j<n_expsyn; ++j) {
+ ASSERT_EQ(1u, expsyn_target_loc_map.count(m[j].target));
+ EXPECT_EQ(expsyn_target_loc_map.at(m[j].target), m[j].loc);
+
+ auto cv = geom.location_cv(i, m[j].loc, cv_prefer::cv_nonempty);
+ target_cv[j] = cv;
+ ++cv_expsyn_count[cv];
+
+ double event_weight = m[j].target+100*i;
+ expected_uncoalesced_value[j] = event_weight;
+ expected_coalesced_cv_value[cv] += event_weight;
+ }
+
+ for (unsigned j = 0; j<n_expsyn; ++j) {
+ if (coalesce_synapses) {
+ EXPECT_EQ(cv_expsyn_count.at(target_cv[j]), m[j].multiplicity);
+ }
+ else {
+ EXPECT_EQ(1u, m[j].multiplicity);
+ }
+ }
+
+ for (unsigned j = 0; j<n_expsyn; ++j) {
+ double expected_value = coalesce_synapses?
+ expected_coalesced_cv_value[target_cv[j]]:
+ expected_uncoalesced_value[j];
+
+ double value = deref(h.raw_handles[j]);
+
+ EXPECT_NEAR(expected_value, value, 0.01); // g values will have decayed a little.
+ }
+ }
+}
+
template <typename Backend>
void run_ion_density_probe_test(const context& ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
@@ -244,34 +483,39 @@ void run_ion_density_probe_test(const context& ctx) {
rec.catalogue() = cat;
// Probe (0, 0): ca internal on CV 0.
- rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc0, "ca"});
+ rec.add_probe(0, 0, cable_probe_ion_int_concentration{loc0, "ca"});
// Probe (0, 1): ca internal on CV 1.
- rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc1, "ca"});
+ rec.add_probe(0, 0, cable_probe_ion_int_concentration{loc1, "ca"});
// Probe (0, 2): ca internal on CV 2.
- rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc2, "ca"});
+ rec.add_probe(0, 0, cable_probe_ion_int_concentration{loc2, "ca"});
// Probe (0, 3): ca external on CV 0.
- rec.add_probe(0, 0, cell_probe_ion_ext_concentration{loc0, "ca"});
+ rec.add_probe(0, 0, cable_probe_ion_ext_concentration{loc0, "ca"});
// Probe (0, 4): ca external on CV 1.
- rec.add_probe(0, 0, cell_probe_ion_ext_concentration{loc1, "ca"});
+ rec.add_probe(0, 0, cable_probe_ion_ext_concentration{loc1, "ca"});
// Probe (0, 5): ca external on CV 2.
- rec.add_probe(0, 0, cell_probe_ion_ext_concentration{loc2, "ca"});
-
+ rec.add_probe(0, 0, cable_probe_ion_ext_concentration{loc2, "ca"});
+
// Probe (0, 6): na internal on CV 0.
- rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc0, "na"});
+ rec.add_probe(0, 0, cable_probe_ion_int_concentration{loc0, "na"});
// Probe (0, 7): na internal on CV 2.
- rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc2, "na"});
+ rec.add_probe(0, 0, cable_probe_ion_int_concentration{loc2, "na"});
// Probe (0, 8): write_ca2 state 's' in CV 0.
- rec.add_probe(0, 0, cell_probe_density_state{loc0, "write_ca2", "s"});
+ rec.add_probe(0, 0, cable_probe_density_state{loc0, "write_ca2", "s"});
// Probe (0, 9): write_ca2 state 's' in CV 1.
- rec.add_probe(0, 0, cell_probe_density_state{loc1, "write_ca2", "s"});
+ rec.add_probe(0, 0, cable_probe_density_state{loc1, "write_ca2", "s"});
// Probe (0, 10): write_ca2 state 's' in CV 2.
- rec.add_probe(0, 0, cell_probe_density_state{loc2, "write_ca2", "s"});
+ rec.add_probe(0, 0, cable_probe_density_state{loc2, "write_ca2", "s"});
+
+ // Probe (0, 11): na internal on whole cell.
+ rec.add_probe(0, 0, cable_probe_ion_int_concentration_cell{"na"});
+ // Probe (0, 12): ca external on whole cell.
+ rec.add_probe(0, 0, cable_probe_ion_ext_concentration_cell{"ca"});
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell lcell(*ctx);
lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
@@ -281,20 +525,24 @@ void run_ion_density_probe_test(const context& ctx) {
// CV 0, and so probe (0, 8) should have been discarded. All other probes
// should be in the map.
- EXPECT_EQ(11u, rec.num_probes(0));
- EXPECT_EQ(9u, probe_map.size());
+ EXPECT_EQ(13u, rec.num_probes(0));
+ EXPECT_EQ(11u, probe_map.size());
+
+ auto probe_scalar_handle = [&](cell_member_type x) {
+ return util::get<fvm_probe_scalar>(probe_map.at(x).handle.info).raw_handles[0];
+ };
- probe_handle ca_int_cv0 = probe_map.at({0, 0}).handle;
- probe_handle ca_int_cv1 = probe_map.at({0, 1}).handle;
- probe_handle ca_int_cv2 = probe_map.at({0, 2}).handle;
- probe_handle ca_ext_cv0 = probe_map.at({0, 3}).handle;
- probe_handle ca_ext_cv1 = probe_map.at({0, 4}).handle;
- probe_handle ca_ext_cv2 = probe_map.at({0, 5}).handle;
+ probe_handle ca_int_cv0 = probe_scalar_handle({0, 0});
+ probe_handle ca_int_cv1 = probe_scalar_handle({0, 1});
+ probe_handle ca_int_cv2 = probe_scalar_handle({0, 2});
+ probe_handle ca_ext_cv0 = probe_scalar_handle({0, 3});
+ probe_handle ca_ext_cv1 = probe_scalar_handle({0, 4});
+ probe_handle ca_ext_cv2 = probe_scalar_handle({0, 5});
EXPECT_EQ(0u, probe_map.count({0, 6}));
- probe_handle na_int_cv2 = probe_map.at({0, 7}).handle;
+ probe_handle na_int_cv2 = probe_scalar_handle({0, 7});
EXPECT_EQ(0u, probe_map.count({0, 8}));
- probe_handle write_ca2_s_cv1 = probe_map.at({0, 9}).handle;
- probe_handle write_ca2_s_cv2 = probe_map.at({0, 10}).handle;
+ probe_handle write_ca2_s_cv1 = probe_scalar_handle({0, 9});
+ probe_handle write_ca2_s_cv2 = probe_scalar_handle({0, 10});
// Ion concentrations should have been written in initialization.
// For CV 1, calcium concentration should be mean of the two values
@@ -316,129 +564,484 @@ void run_ion_density_probe_test(const context& ctx) {
EXPECT_EQ(2.75, deref(write_ca2_s_cv1));
EXPECT_EQ(2.75, deref(write_ca2_s_cv2));
EXPECT_NE(write_ca2_s_cv1, write_ca2_s_cv2);
+
+ // For the all cell sodium concentration probe, check that metadata
+ // cables comprise a cable for all of CV 1 and a cable for all of CV 2:
+ // half coverage of CV 1 by sodium ion-requiring mechanism implies
+ // ion values are valid across whole CV.
+ //
+ // Implementation detail has that the cables (and corresponding raw handles) are
+ // sorted by CV in the fvm_probe_weighted_multi object; this is assumed
+ // below.
+
+ auto* p_ptr = util::get_if<fvm_probe_multi>(probe_map.at({0, 11}).handle.info);
+ ASSERT_TRUE(p_ptr);
+ fvm_probe_multi& na_int_all_info = *p_ptr;
+
+ auto* m_ptr = util::get_if<mcable_list>(na_int_all_info.metadata);
+ ASSERT_TRUE(m_ptr);
+ mcable_list na_int_all_metadata = *m_ptr;
+
+ ASSERT_EQ(2u, na_int_all_metadata.size());
+ ASSERT_EQ(2u, na_int_all_info.raw_handles.size());
+
+ EXPECT_DOUBLE_EQ(1./3., na_int_all_metadata[0].prox_pos);
+ EXPECT_DOUBLE_EQ(2./3., na_int_all_metadata[0].dist_pos);
+ EXPECT_DOUBLE_EQ(2./3., na_int_all_metadata[1].prox_pos);
+ EXPECT_DOUBLE_EQ(1., na_int_all_metadata[1].dist_pos);
+ EXPECT_EQ(na_int_cv2, na_int_all_info.raw_handles[1]);
+ EXPECT_EQ(na_int_cv2-1, na_int_all_info.raw_handles[0]);
+
+ p_ptr = util::get_if<fvm_probe_multi>(probe_map.at({0, 12}).handle.info);
+ ASSERT_TRUE(p_ptr);
+ fvm_probe_multi& ca_ext_all_info = *p_ptr;
+
+ m_ptr = util::get_if<mcable_list>(ca_ext_all_info.metadata);
+ ASSERT_TRUE(m_ptr);
+ mcable_list ca_ext_all_metadata = *m_ptr;
+
+ ASSERT_EQ(3u, ca_ext_all_metadata.size());
+ ASSERT_EQ(3u, ca_ext_all_info.raw_handles.size());
+
+ EXPECT_DOUBLE_EQ(0., ca_ext_all_metadata[0].prox_pos);
+ EXPECT_DOUBLE_EQ(1./3., ca_ext_all_metadata[1].prox_pos);
+ EXPECT_DOUBLE_EQ(2./3., ca_ext_all_metadata[2].prox_pos);
+ EXPECT_EQ(ca_ext_cv0, ca_ext_all_info.raw_handles[0]);
+ EXPECT_EQ(ca_ext_cv1, ca_ext_all_info.raw_handles[1]);
+ EXPECT_EQ(ca_ext_cv2, ca_ext_all_info.raw_handles[2]);
}
template <typename Backend>
-void run_ion_current_probe_test(const context& ctx) {
+void run_partial_density_probe_test(const context& ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const fvm_value_type* p) { return backend_access<Backend>::deref(p); };
- // Use test mechanism fixed_ica_current, and a derived mechanism for sodium, to
- // write to specific ion currents.
+ // Use test mechanism param_as_state to query averaged state values in CVs with
+ // partial coverage by the mechanism.
auto cat = make_unit_test_catalogue();
- cat.derive("fixed_ina_current", "fixed_ica_current", {}, {{"ca", "na"}});
cable_cell cells[2];
- // Simple constant diameter cable, 3 CVs.
+ // Each cell is a simple constant diameter cable, with 3 CVs each.
- cells[0] = cable_cell(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 1.}, 0}}, {mnpos, 0u}));
+ morphology m(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 1.}, 0}}, {mnpos, 0u}));
+ cells[0] = cable_cell(m);
cells[0].default_parameters.discretization = cv_policy_fixed_per_branch(3);
- // Calcium ions everywhere, half with current density jca0, half with jca1.
- // Sodium ions only on distal half, with current densitry jna1.
-
- const double jca0 = 1.5; // [A/m²]
- const double jca1 = 2.0;
- const double jna1 = 2.5;
-
- // Scaling factor 0.1 is to convert our current densities in [A/m²] to NMODL units [mA/cm²].
-
- cells[0].paint(mcable{0, 0., 0.5}, mechanism_desc("fixed_ica_current").set("current_density", 0.1*jca0));
- cells[0].paint(mcable{0, 0.5, 1.}, mechanism_desc("fixed_ica_current").set("current_density", 0.1*jca1));
- cells[0].paint(mcable{0, 0.5, 1.}, mechanism_desc("fixed_ina_current").set("current_density", 0.1*jna1));
-
- // Make a second cable cell, with same layout but 3 times the current.
-
- cells[1] = cable_cell(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 1.}, 0}}, {mnpos, 0u}));
+ cells[1] = cable_cell(m);
cells[1].default_parameters.discretization = cv_policy_fixed_per_branch(3);
- cells[1].paint(mcable{0, 0., 0.5}, mechanism_desc("fixed_ica_current").set("current_density", 0.3*jca0));
- cells[1].paint(mcable{0, 0.5, 1.}, mechanism_desc("fixed_ica_current").set("current_density", 0.3*jca1));
- cells[1].paint(mcable{0, 0.5, 1.}, mechanism_desc("fixed_ina_current").set("current_density", 0.3*jna1));
-
- // Place probes in each CV on cell 0, plus one in the last CV on cell 1.
-
- mlocation loc0{0, 0.1};
- mlocation loc1{0, 0.5};
- mlocation loc2{0, 0.9};
+ // Paint the mechanism on every second 10% interval of each cell.
+ // Expected values on a CV are the weighted mean of the parameter values
+ // over the intersections of the support and the CV.
+ //
+ // Cell 0: [0.0, 0.1], [0.2, 0.3] [0.4, 0.5] [0.6, 0.7] [0.8, 0.9].
+ // Param values: 2 3 4 5 6
+ // Expected values:
+ // CV 0: 2.5
+ // CV 1: 4.4
+ // CV 2: 5.75
+ //
+ // Cell 1: [0.1, 0.2], [0.3, 0.4] [0.5, 0.6] [0.7, 0.8] [0.9, 1.0].
+ // Param values: 7 8 9 10 11
+ // Expected values:
+ // CV 3: 7.25
+ // CV 4: 8.6
+ // CV 5: 10.5
+
+ auto mk_mech = [](double param) { return mechanism_desc("param_as_state").set("p", param); };
+
+ cells[0].paint(mcable{0, 0.0, 0.1}, mk_mech(2));
+ cells[0].paint(mcable{0, 0.2, 0.3}, mk_mech(3));
+ cells[0].paint(mcable{0, 0.4, 0.5}, mk_mech(4));
+ cells[0].paint(mcable{0, 0.6, 0.7}, mk_mech(5));
+ cells[0].paint(mcable{0, 0.8, 0.9}, mk_mech(6));
+
+ cells[1].paint(mcable{0, 0.1, 0.2}, mk_mech(7));
+ cells[1].paint(mcable{0, 0.3, 0.4}, mk_mech(8));
+ cells[1].paint(mcable{0, 0.5, 0.6}, mk_mech(9));
+ cells[1].paint(mcable{0, 0.7, 0.8}, mk_mech(10));
+ cells[1].paint(mcable{0, 0.9, 1.0}, mk_mech(11));
+
+ // Place probes in the middle of each 10% interval, i.e. at 0.05, 0.15, etc.
+ struct test_probe {
+ double pos;
+ double expected[2]; // Expected value in each cell, NAN => n/a
+ } test_probes[] = {
+ { 0.05, { 2.5, NAN }}, // CV 0, 3
+ { 0.15, { NAN, 7.25 }}, // CV 0, 3
+ { 0.25, { 2.5, NAN }}, // CV 0, 3
+ { 0.35, { NAN, 8.6 }}, // CV 1, 4
+ { 0.45, { 4.4, NAN }}, // CV 1, 4
+ { 0.55, { NAN, 8.6 }}, // CV 1, 4
+ { 0.65, { 4.4, NAN }}, // CV 1, 4
+ { 0.75, { NAN, 10.5 }}, // CV 2, 5
+ { 0.85, { 5.75, NAN }}, // CV 2, 5
+ { 0.95, { NAN, 10.5 }} // CV 2, 5
+ };
cable1d_recipe rec(cells);
rec.catalogue() = cat;
- // Probe (0, 0): ica on CV 0.
- rec.add_probe(0, 0, cell_probe_ion_current_density{loc0, "ca"});
- // Probe (0, 1): ica on CV 1.
- rec.add_probe(0, 0, cell_probe_ion_current_density{loc1, "ca"});
- // Probe (0, 2): ica on CV 2.
- rec.add_probe(0, 0, cell_probe_ion_current_density{loc2, "ca"});
-
- // Probe (0, 3): ina on CV 0.
- rec.add_probe(0, 0, cell_probe_ion_current_density{loc0, "na"});
- // Probe (0, 4): ina on CV 1.
- rec.add_probe(0, 0, cell_probe_ion_current_density{loc1, "na"});
- // Probe (0, 5): ina on CV 2.
- rec.add_probe(0, 0, cell_probe_ion_current_density{loc2, "na"});
-
- // Probe (0, 6): total ion current density on CV 0.
- rec.add_probe(0, 0, cell_probe_total_ionic_current_density{loc0});
- // Probe (0, 7): total ion current density on CV 1.
- rec.add_probe(0, 0, cell_probe_total_ionic_current_density{loc1});
- // Probe (0, 8): total ion current density on CV 2.
- rec.add_probe(0, 0, cell_probe_total_ionic_current_density{loc2});
-
- // Probe (1, 0): ica on CV 5 (CV 2 of cell 1).
- rec.add_probe(1, 0, cell_probe_ion_current_density{loc2, "ca"});
- // Probe (1, 1): total ion current density on CV 5 (CV 2 of cell 1).
- rec.add_probe(1, 0, cell_probe_total_ionic_current_density{loc2});
+ for (auto tp: test_probes) {
+ rec.add_probe(0, 0, cable_probe_density_state{mlocation{0, tp.pos}, "param_as_state", "s"});
+ rec.add_probe(1, 0, cable_probe_density_state{mlocation{0, tp.pos}, "param_as_state", "s"});
+ }
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
- probe_association_map<probe_handle> probe_map;
+ probe_association_map<fvm_probe_info> probe_map;
fvm_cell lcell(*ctx);
lcell.initialize({0, 1}, rec, cell_to_intdom, targets, probe_map);
- // Should be no sodium ion instantiated on CV 0, so probe (0, 3) should
- // have been silently discared.
+ // There should be 10 probes on each cell, but only 10 in total in the probe map,
+ // as only those probes that are in the mechanism support should have an entry.
- EXPECT_EQ(9u, rec.num_probes(0));
- EXPECT_EQ(2u, rec.num_probes(1));
+ EXPECT_EQ(10u, rec.num_probes(0));
+ EXPECT_EQ(10u, rec.num_probes(1));
EXPECT_EQ(10u, probe_map.size());
- probe_handle ica_cv0 = probe_map.at({0, 0}).handle;
- probe_handle ica_cv1 = probe_map.at({0, 1}).handle;
- probe_handle ica_cv2 = probe_map.at({0, 2}).handle;
- EXPECT_EQ(0u, probe_map.count({0, 3}));
- probe_handle ina_cv1 = probe_map.at({0, 4}).handle;
- probe_handle ina_cv2 = probe_map.at({0, 5}).handle;
- probe_handle i_cv0 = probe_map.at({0, 6}).handle;
- probe_handle i_cv1 = probe_map.at({0, 7}).handle;
- probe_handle i_cv2 = probe_map.at({0, 8}).handle;
+ auto probe_scalar_handle = [&](cell_member_type x) {
+ return util::get<fvm_probe_scalar>(probe_map.at(x).handle.info).raw_handles[0];
+ };
- probe_handle ica_cv5 = probe_map.at({1, 0}).handle;
- probe_handle i_cv5 = probe_map.at({1, 1}).handle;
+ // Check probe values against expected values.
+
+ cell_lid_type probe_lid = 0;
+ for (auto tp: test_probes) {
+ for (cell_gid_type gid: {0, 1}) {
+ cell_member_type probe_id{gid, probe_lid};
+ if (std::isnan(tp.expected[gid])) {
+ EXPECT_EQ(0u, probe_map.count(probe_id));
+ }
+ else {
+ probe_handle h = probe_scalar_handle(probe_id);
+ EXPECT_DOUBLE_EQ(tp.expected[gid], deref(h));
+ }
+ }
+ ++probe_lid;
+ }
+}
- // Integrate cell for a bit, and check that currents add up as we expect.
+template <typename Backend>
+void run_axial_and_ion_current_probe_sampled_test(const context& ctx) {
+ // On a passive cable in steady-state, the capacitive membrane current will be zero,
+ // and the axial currents should balance the ionic membrane currents in any CV.
+ //
+ // Membrane currents not associated with an ion are not visible, so we use a custom
+ // mechanism 'ca_linear' that presents a passive, constant conductance membrane current
+ // as a calcium ion current.
- lcell.integrate(0.01, 0.0025, {}, {});
+ auto cat = make_unit_test_catalogue();
+
+ // Cell is a tapered cable with 3 CVs.
+
+ morphology m(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 0.8}, 0}}, {mnpos, 0u}));
+ cable_cell cell(m);
+
+ const unsigned n_cv = 3;
+ cv_policy policy = cv_policy_fixed_per_branch(n_cv);
+ cell.default_parameters.discretization = policy;
+
+ cell.place(mlocation{0, 0}, i_clamp(0, INFINITY, 0.3));
+
+ // The time constant will be membrane capacitance / membrane conductance.
+ // For τ = 0.1 ms, set conductance to 0.01 S/cm² and membrance capacitance
+ // to 0.01 F/m².
+
+ cell.paint(reg::all(), mechanism_desc("ca_linear").set("g", 0.01)); // [S/cm²]
+ cell.default_parameters.membrane_capacitance = 0.01; // [F/m²]
+ const double tau = 0.1; // [ms]
+
+ cable1d_recipe rec(cell);
+ rec.catalogue() = cat;
+
+ // Place axial current probes at CV boundaries and make a cell-wide probe for
+ // total ionic membrane current.
- EXPECT_DOUBLE_EQ(jca0, deref(ica_cv0));
- EXPECT_DOUBLE_EQ((jca0+jca1)/2, deref(ica_cv1));
- EXPECT_DOUBLE_EQ(jca1, deref(ica_cv2));
+ mlocation_list cv_boundaries;
+ util::assign(cv_boundaries,
+ util::filter(thingify(policy.cv_boundary_points(cell), cell.provider()),
+ [](mlocation loc) { return loc.pos!=0 && loc.pos!=1; }));
- EXPECT_DOUBLE_EQ(jna1/2, deref(ina_cv1));
- EXPECT_DOUBLE_EQ(jna1, deref(ina_cv2));
+ ASSERT_EQ(n_cv-1, cv_boundaries.size());
+ const unsigned n_axial_probe = n_cv-1;
- EXPECT_DOUBLE_EQ(jca0, deref(i_cv0));
- EXPECT_DOUBLE_EQ(jna1/2+jca0/2+jca1/2, deref(i_cv1));
- EXPECT_DOUBLE_EQ(jna1+jca1, deref(i_cv2));
+ for (mlocation loc: cv_boundaries) {
+ // Probe ids will be be (0, 0), (0, 1), etc.
+ rec.add_probe(0, 0, cable_probe_axial_current{loc});
+ }
+
+ // Use tag 1 to indicate it's a whole-cell probe.
+ rec.add_probe(0, 1, cable_probe_total_ion_current_cell{});
+
+ partition_hint_map phints = {
+ {cell_kind::cable, {partition_hint::max_size, partition_hint::max_size, true}}
+ };
+ simulation sim(rec, partition_load_balance(rec, ctx, phints), ctx);
+
+ // Take a sample at 20 tau, and run sim for just a bit longer.
+
+ std::vector<double> i_axial(n_axial_probe);
+ std::vector<double> i_memb;
+
+ sim.add_sampler(all_probes, explicit_schedule({20*tau}),
+ [&](cell_member_type probe_id, probe_tag tag, util::any_ptr metadata,
+ std::size_t n_sample, const sample_record* samples)
+ {
+ // Expect exactly one sample.
+ ASSERT_EQ(1u, n_sample);
+
+ if (tag==1) { // (whole cell probe)
+ const mcable_list* m = util::any_cast<const mcable_list*>(metadata);
+ ASSERT_NE(nullptr, m);
+ // Metadata should comprise one cable per CV.
+ ASSERT_EQ(n_cv, m->size());
+
+ const cable_sample_range* s = util::any_cast<const cable_sample_range*>(samples[0].data);
+ ASSERT_NE(nullptr, s);
+ ASSERT_EQ(s->first+n_cv, s->second);
+
+ for (const double* p = s->first; p!=s->second; ++p) {
+ i_memb.push_back(*p);
+ }
+ }
+ else { // axial current probe
+ // Probe id tells us which axial current this is.
+ ASSERT_LT(probe_id.index, n_axial_probe);
+
+ const mlocation* m = util::any_cast<const mlocation*>(metadata);
+ ASSERT_NE(nullptr, m);
+
+ const double* s = util::any_cast<const double*>(samples[0].data);
+ ASSERT_NE(nullptr, s);
+
+ i_axial.at(probe_id.index) = *s;
+ }
+ });
+
+ const double dt = 0.025; // [ms]
+ sim.run(20*tau+dt, dt);
+
+ ASSERT_EQ(n_cv, i_memb.size());
+
+ for (unsigned i = 0; i<n_cv; ++i) {
+ // Axial currents are in the distal (increasing CV index) direction,
+ // while membrane currents are from intra- to extra-cellular medium.
+ //
+ // Net outward flux from CV is axial current on distal side, minus
+ // axial current on proximal side, plus membrane current, and should
+ // sum to zero.
+
+ double net_axial_flux = i<n_axial_probe? i_axial[i]: 0;
+ net_axial_flux -= i>0? i_axial[i-1]: 0;
+
+ EXPECT_TRUE(testing::near_relative(net_axial_flux, -i_memb[i], 1e-6));
+ }
+}
+
+
+// Run given cells taking samples from the provied probes on one of the cells.
+// Use default mechanism catalogue augmented by unit test specific mechanisms.
+// (Timestep fixed at 0.025 ms).
+
+template <typename SampleData, typename SampleMeta = void>
+auto run_simple_samplers(
+ const arb::context& ctx,
+ double t_end,
+ const std::vector<cable_cell>& cells,
+ cell_gid_type probe_cell,
+ const std::vector<util::any>& probe_addrs,
+ const std::vector<double>& when)
+{
+ cable1d_recipe rec(cells, false);
+ rec.catalogue() = make_unit_test_catalogue(global_default_catalogue());
+ unsigned n_probe = probe_addrs.size();
+
+ for (auto& addr: probe_addrs) {
+ rec.add_probe(probe_cell, 0, addr);
+ }
+
+ partition_hint_map phints = {
+ {cell_kind::cable, {partition_hint::max_size, partition_hint::max_size, true}}
+ };
+ simulation sim(rec, partition_load_balance(rec, ctx, phints), ctx);
+
+ std::vector<trace_data<SampleData, SampleMeta>> traces(n_probe);
+ for (unsigned i = 0; i<n_probe; ++i) {
+ sim.add_sampler(one_probe({probe_cell, i}), explicit_schedule(when), make_simple_sampler(traces[i]));
+ }
+
+ sim.run(t_end, 0.025);
+ return traces;
+}
+
+template <typename SampleData, typename SampleMeta = void>
+auto run_simple_sampler(
+ const arb::context& ctx,
+ double t_end,
+ const std::vector<cable_cell>& cells,
+ cell_gid_type probe_cell,
+ const util::any& probe_addr,
+ const std::vector<double>& when)
+{
+ return run_simple_samplers<SampleData, SampleMeta>(ctx, t_end, cells, probe_cell, {probe_addr}, when).at(0);
+}
+
+template <typename Backend>
+void run_v_sampled_test(const context& ctx) {
+ cable_cell bs = make_cell_ball_and_stick(false);
+ bs.default_parameters.discretization = cv_policy_fixed_per_branch(1);
+
+ std::vector<cable_cell> cells = {bs, bs};
+
+ // Add stims, up to 0.5 ms on cell 0, up to 1.0 ms on cell 1, so that
+ // samples at the same point on each cell will give the same value at
+ // 0.3 ms, but different at 0.6 ms.
+
+ cells[0].place(mlocation{1, 1}, i_clamp(0, 0.5, 1.));
+ cells[1].place(mlocation{1, 1}, i_clamp(0, 1.0, 1.));
+ mlocation probe_loc{1, 0.2};
+
+ const double t_end = 1.; // [ms]
+ std::vector<double> when = {0.3, 0.6}; // Sample at 0.3 and 0.6 ms.
+
+ auto trace0 = run_simple_sampler<double, mlocation>(ctx, t_end, cells, 0, cable_probe_membrane_voltage{probe_loc}, when);
+ EXPECT_EQ(2u, trace0.size());
+ EXPECT_EQ(probe_loc, trace0.metadata.value());
+
+ auto trace1 = run_simple_sampler<double, mlocation>(ctx, t_end, cells, 1, cable_probe_membrane_voltage{probe_loc}, when);
+ EXPECT_EQ(2u, trace1.size());
+ EXPECT_EQ(probe_loc, trace1.metadata.value());
+
+ EXPECT_EQ(trace0[0].t, trace1[0].t);
+ EXPECT_EQ(trace0[0].v, trace1[0].v);
+
+ EXPECT_EQ(trace0[1].t, trace1[1].t);
+ EXPECT_NE(trace0[1].v, trace1[1].v);
+}
+
+template <typename Backend>
+void run_total_current_test(const context& ctx) {
+ // Model two passive Y-shaped cells with a similar but not identical
+ // time constant Ï„.
+ //
+ // Sample each cell's total membrane currents at the same time,
+ // approximately equal to the time constants Ï„.
+ //
+ // Net current flux in each cell should be zero, but currents should
+ // differ between the cells.
+
+ cable_cell cell(make_y_morphology());
+
+ const unsigned n_cv_per_branch = 3;
+ const unsigned n_branch = 3;
+
+ // The time constant will be membrane capacitance / membrane conductance.
+ // For τ = 0.1 ms, set conductance to 0.01 S/cm² and membrance capacitance
+ // to 0.01 F/m².
+
+ const double tau = 0.1; // [ms]
+ cell.place(mlocation{0, 0}, i_clamp(0, INFINITY, 0.3));
+
+ cell.paint(reg::all(), mechanism_desc("ca_linear").set("g", 0.01)); // [S/cm²]
+ cell.default_parameters.membrane_capacitance = 0.01; // [F/m²]
+
+ std::vector<cable_cell> cells = {cell, cell};
+
+ // Tweak membrane capacitance on cells[1] so as to change dynamics a bit.
+ cells[1].default_parameters.membrane_capacitance = 0.009; // [F/m²]
+
+ // We'll run each set of tests twice: once with a trivial (zero-volume) CV
+ // at the fork points, and once with a non-trivial CV centred on the fork
+ // point.
+
+ trace_data<std::vector<double>, mcable_list> traces[2];
+ trace_data<std::vector<double>, mcable_list> ion_traces[2];
+
+ // Run the cells sampling at Ï„ and 20Ï„ for both total membrane
+ // current and total membrane ionic current.
+
+ auto run_cells = [&](bool interior_forks) {
+ auto flags = interior_forks? cv_policy_flag::interior_forks: cv_policy_flag::none;
+ cv_policy policy = cv_policy_fixed_per_branch(n_cv_per_branch, flags);
+ for (auto& c: cells) { c.default_parameters.discretization = policy; }
- // Currents on cell 1 should be 3 times those on cell 0.
+ for (unsigned i = 0; i<2; ++i) {
+ SCOPED_TRACE(i);
- EXPECT_DOUBLE_EQ(jca1*3, deref(ica_cv5));
- EXPECT_DOUBLE_EQ((jna1+jca1)*3, deref(i_cv5));
+ const double t_end = 21*tau; // [ms]
+
+ traces[i] = run_simple_sampler<std::vector<double>, mcable_list>(ctx, t_end, cells, i,
+ cable_probe_total_current_cell{}, {tau, 20*tau});
+
+ ion_traces[i] = run_simple_sampler<std::vector<double>, mcable_list>(ctx, t_end, cells, i,
+ cable_probe_total_ion_current_cell{}, {tau, 20*tau});
+
+ ASSERT_EQ(2u, traces[i].size());
+ ASSERT_EQ(2u, ion_traces[i].size());
+
+ // Check metadata size:
+ // * With trivial forks, should have n_cv_per_branch*n_branch cables; zero-length cables
+ // associated with the trivial CVs at forks should not be included.
+ // * With nontrivial forks, we'll have an extra cable at the head of each branch, which
+ // for all but the root branch will be a component cable of the CV on the fork.
+ //
+ // Total membrane current and total ionic mebrane current should have the
+ // same support and same metadata.
+
+ ASSERT_EQ((n_cv_per_branch+(int)interior_forks)*n_branch, traces[i].metadata.value().size());
+ EXPECT_EQ(ion_traces[i].metadata, traces[i].metadata);
+ EXPECT_EQ(ion_traces[i][0].v.size(), traces[i][0].v.size());
+ EXPECT_EQ(ion_traces[i][1].v.size(), traces[i][1].v.size());
+
+ // Check total membrane currents are individually non-zero, but sum is, both
+ // at t=Ï„ (j=0) and t=20Ï„ (j=1).
+
+ for (unsigned j: {0u, 1u}) {
+ double max_abs_i_memb = 0;
+ double sum_i_memb = 0;
+ for (auto i_memb: traces[i][j].v) {
+ EXPECT_NE(0.0, i_memb);
+ max_abs_i_memb = std::max(max_abs_i_memb, std::abs(i_memb));
+ sum_i_memb += i_memb;
+ }
+
+ EXPECT_NEAR(0.0, sum_i_memb, 1e-6*max_abs_i_memb);
+ }
+
+ // Confirm that total and ion currents differ at Ï„ but are close at 20Ï„.
+
+ for (unsigned k = 0; k<traces[i].size(); ++k) {
+ const double rtol_large = 1e-3;
+ EXPECT_FALSE(testing::near_relative(traces[i][0].v.at(k), ion_traces[i][0].v.at(k), rtol_large));
+ }
+
+ for (unsigned k = 0; k<traces[i].size(); ++k) {
+ const double rtol_small = 1e-6;
+ EXPECT_TRUE(testing::near_relative(traces[i][1].v.at(k), ion_traces[i][1].v.at(k), rtol_small));
+ }
+
+ }
+
+ // Total membrane currents should differ between the two cells at t=Ï„.
+
+ for (unsigned k = 0; k<traces[0][0].v.size(); ++k) {
+ EXPECT_NE(traces[0][0].v.at(k), traces[1][0].v.at(k));
+ }
+ };
+
+ {
+ SCOPED_TRACE("trivial fork CV");
+ run_cells(false);
+ }
+
+ {
+ SCOPED_TRACE("non-trival fork CV");
+ run_cells(true);
+ }
}
TEST(probe, multicore_v_i) {
@@ -446,22 +1049,46 @@ TEST(probe, multicore_v_i) {
run_v_i_probe_test<multicore::backend>(ctx);
}
+TEST(probe, multicore_v_cell) {
+ context ctx = make_context();
+ run_v_cable_probe_test<multicore::backend>(ctx);
+}
+
+TEST(probe, multicore_v_sampled) {
+ context ctx = make_context();
+ run_v_sampled_test<multicore::backend>(ctx);
+}
+
TEST(probe, multicore_expsyn_g) {
context ctx = make_context();
- SCOPED_TRACE("uncoalesced synapses");
run_expsyn_g_probe_test<multicore::backend>(ctx, false);
- SCOPED_TRACE("coalesced synapses");
run_expsyn_g_probe_test<multicore::backend>(ctx, true);
}
+TEST(probe, multicore_expsyn_g_cell) {
+ context ctx = make_context();
+ run_expsyn_g_cable_probe_test<multicore::backend>(ctx, false);
+ run_expsyn_g_cable_probe_test<multicore::backend>(ctx, true);
+}
+
TEST(probe, multicore_ion_conc) {
context ctx = make_context();
run_ion_density_probe_test<multicore::backend>(ctx);
}
-TEST(probe, multicore_ion_currents) {
+TEST(probe, multicore_axial_and_ion_current) {
+ context ctx = make_context();
+ run_axial_and_ion_current_probe_sampled_test<multicore::backend>(ctx);
+}
+
+TEST(probe, multicore_partial_density) {
+ context ctx = make_context();
+ run_partial_density_probe_test<multicore::backend>(ctx);
+}
+
+TEST(probe, multicore_total_current) {
context ctx = make_context();
- run_ion_current_probe_test<multicore::backend>(ctx);
+ run_total_current_test<multicore::backend>(ctx);
}
#ifdef ARB_GPU_ENABLED
@@ -472,16 +1099,34 @@ TEST(probe, gpu_v_i) {
}
}
+TEST(probe, gpu_v_cell) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ if (has_gpu(ctx)) {
+ run_v_cable_probe_test<gpu::backend>(ctx);
+ }
+}
+
+TEST(probe, gpu_v_sampled) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ run_v_sampled_test<multicore::backend>(ctx);
+}
+
TEST(probe, gpu_expsyn_g) {
context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
if (has_gpu(ctx)) {
- SCOPED_TRACE("uncoalesced synapses");
run_expsyn_g_probe_test<gpu::backend>(ctx, false);
- SCOPED_TRACE("coalesced synapses");
run_expsyn_g_probe_test<gpu::backend>(ctx, true);
}
}
+TEST(probe, gpu_expsyn_g_cell) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ if (has_gpu(ctx)) {
+ run_expsyn_g_cable_probe_test<gpu::backend>(ctx, false);
+ run_expsyn_g_cable_probe_test<gpu::backend>(ctx, true);
+ }
+}
+
TEST(probe, gpu_ion_conc) {
context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
if (has_gpu(ctx)) {
@@ -489,11 +1134,24 @@ TEST(probe, gpu_ion_conc) {
}
}
-TEST(probe, gpu_ion_currents) {
+TEST(probe, gpu_axial_and_ion_current) {
context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
if (has_gpu(ctx)) {
- run_ion_current_probe_test<gpu::backend>(ctx);
+ run_axial_and_ion_current_probe_sampled_test<gpu::backend>(ctx);
+ }
+}
+
+TEST(probe, gpu_partial_density) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ if (has_gpu(ctx)) {
+ run_partial_density_probe_test<multicore::backend>(ctx);
}
}
-#endif
+TEST(probe, gpu_total_current) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ if (has_gpu(ctx)) {
+ run_total_current_test<gpu::backend>(ctx);
+ }
+}
+#endif
diff --git a/test/unit/test_schedule.cpp b/test/unit/test_schedule.cpp
index c49a1a0e..e03707fe 100644
--- a/test/unit/test_schedule.cpp
+++ b/test/unit/test_schedule.cpp
@@ -186,7 +186,7 @@ struct skew_adaptor {
explicit skew_adaptor(double power): power_(power) {}
result_type operator()() {
- constexpr double scale = max()-min();
+ constexpr double scale = (double)(max()-min());
constexpr double ooscale = 1./scale;
double x = ooscale*(G_()-min());
diff --git a/test/unit/unit_test_catalogue.cpp b/test/unit/unit_test_catalogue.cpp
index 8d01b374..aeb7f39d 100644
--- a/test/unit/unit_test_catalogue.cpp
+++ b/test/unit/unit_test_catalogue.cpp
@@ -7,8 +7,10 @@
#include "backends/multicore/fvm.hpp"
#include "unit_test_catalogue.hpp"
+#include "mechanisms/ca_linear.hpp"
#include "mechanisms/celsius_test.hpp"
#include "mechanisms/diam_test.hpp"
+#include "mechanisms/param_as_state.hpp"
#include "mechanisms/test0_kin_diff.hpp"
#include "mechanisms/test_linear_state.hpp"
#include "mechanisms/test_linear_init.hpp"
@@ -50,11 +52,13 @@ c.register_implementation(#x, testing::make_mechanism_##x<gpu::backend>());
using namespace arb;
-mechanism_catalogue make_unit_test_catalogue() {
- mechanism_catalogue cat;
+mechanism_catalogue make_unit_test_catalogue(const mechanism_catalogue& from) {
+ mechanism_catalogue cat(from);
+ ADD_MECH(cat, ca_linear)
ADD_MECH(cat, celsius_test)
ADD_MECH(cat, diam_test)
+ ADD_MECH(cat, param_as_state)
ADD_MECH(cat, test_linear_state)
ADD_MECH(cat, test_linear_init)
ADD_MECH(cat, test_linear_init_shuffle)
diff --git a/test/unit/unit_test_catalogue.hpp b/test/unit/unit_test_catalogue.hpp
index 6516ffa2..68435023 100644
--- a/test/unit/unit_test_catalogue.hpp
+++ b/test/unit/unit_test_catalogue.hpp
@@ -2,4 +2,4 @@
#include <arbor/mechcat.hpp>
-arb::mechanism_catalogue make_unit_test_catalogue();
+arb::mechanism_catalogue make_unit_test_catalogue(const arb::mechanism_catalogue& from = {});
--
GitLab