diff --git a/arbor/backends/event.hpp b/arbor/backends/event.hpp index d14bfdeaf1b269acbfc8ec7eaff3f37d4b5f9b37..b69f834107940f845dd01c9c6994c4af354f7e20 100644 --- a/arbor/backends/event.hpp +++ b/arbor/backends/event.hpp @@ -49,7 +49,7 @@ inline deliverable_event_data event_data(const deliverable_event& ev) { } -// Sample events (scalar values) +// Sample events (raw values from back-end state). using probe_handle = const fvm_value_type*; diff --git a/arbor/fvm_lowered_cell.hpp b/arbor/fvm_lowered_cell.hpp index 7cd80f5902d8f35491b6dd1ec911eabc1eb5715d..83f934586ae4e7a6a11a16390d9fd7548be220ff 100644 --- a/arbor/fvm_lowered_cell.hpp +++ b/arbor/fvm_lowered_cell.hpp @@ -1,8 +1,11 @@ #pragma once +#include <cstddef> #include <memory> +#include <unordered_map> #include <vector> +#include <arbor/assert.hpp> #include <arbor/common_types.hpp> #include <arbor/cable_cell.hpp> #include <arbor/fvm_types.hpp> @@ -16,6 +19,7 @@ #include "sampler_map.hpp" #include "util/meta.hpp" #include "util/range.hpp" +#include "util/transform.hpp" namespace arb { @@ -28,7 +32,7 @@ struct fvm_integration_result { // A sample for a probe may be derived from multiple 'raw' sampled // values from the backend. -// An `fvm_probe_info` object represents the mapping between +// An `fvm_probe_data` 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. @@ -91,13 +95,13 @@ struct missing_probe_info { std::array<probe_handle, 0> raw_handles; }; -struct fvm_probe_info { - 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)) {} +struct fvm_probe_data { + fvm_probe_data() = default; + fvm_probe_data(fvm_probe_scalar p): info(std::move(p)) {} + fvm_probe_data(fvm_probe_interpolated p): info(std::move(p)) {} + fvm_probe_data(fvm_probe_multi p): info(std::move(p)) {} + fvm_probe_data(fvm_probe_weighted_multi p): info(std::move(p)) {} + fvm_probe_data(fvm_probe_membrane_currents p): info(std::move(p)) {} util::variant< missing_probe_info, @@ -124,6 +128,26 @@ struct fvm_probe_info { explicit operator bool() const { return !util::get_if<missing_probe_info>(info); } }; +// Samplers are tied to probe ids, but one probe id may +// map to multiple probe representations within the mc_cell_group. + +struct probe_association_map { + // Keys are probe id. + + std::unordered_map<cell_member_type, probe_tag> tag; + std::unordered_multimap<cell_member_type, fvm_probe_data> data; + + std::size_t size() const { + arb_assert(tag.size()==data.size()); + return data.size(); + } + + // Return range of fvm_probe_data values associated with probe_id. + auto data_on(cell_member_type probe_id) const { + return util::transform_view(util::make_range(data.equal_range(probe_id)), util::second); + } +}; + // Common base class for FVM implementation on host or gpu back-end. struct fvm_lowered_cell { @@ -134,7 +158,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<fvm_probe_info>& probe_map) = 0; + probe_association_map& 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 d289526e836831569192832082edb30ecae8f9b8..173cad5e744d8b71b3726d5f8693b12eabcc0cf5 100644 --- a/arbor/fvm_lowered_cell_impl.hpp +++ b/arbor/fvm_lowered_cell_impl.hpp @@ -57,7 +57,7 @@ public: const recipe& rec, std::vector<fvm_index_type>& cell_to_intdom, std::vector<target_handle>& target_handles, - probe_association_map<fvm_probe_info>& probe_map) override; + probe_association_map& probe_map) override; fvm_integration_result integrate( value_type tfinal, @@ -141,7 +141,8 @@ private: } // Translate cell probe descriptions into probe handles etc. - fvm_probe_info resolve_probe_address( + void resolve_probe_address( + std::vector<fvm_probe_data>& probe_data, // out parameter const std::vector<cable_cell>& cells, std::size_t cell_idx, const util::any& paddr, @@ -356,7 +357,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<fvm_probe_info>& probe_map) + probe_association_map& probe_map) { using util::any_cast; using util::count_along; @@ -533,6 +534,7 @@ void fvm_lowered_cell_impl<Backend>::initialize( std::vector<index_type> detector_cv; std::vector<value_type> detector_threshold; + std::vector<fvm_probe_data> probe_data; for (auto cell_idx: make_span(ncell)) { cell_gid_type gid = gids[cell_idx]; @@ -542,12 +544,19 @@ void fvm_lowered_cell_impl<Backend>::initialize( detector_threshold.push_back(entry.item.threshold); } - 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(cells, cell_idx, pi.address, D, mech_data, target_handles, mechptr_by_name); + std::vector<probe_info> rec_probes = rec.get_probes(gid); + for (cell_lid_type i: count_along(rec_probes)) { + probe_info& pi = rec_probes[i]; + resolve_probe_address(probe_data, cells, cell_idx, std::move(pi.address), + D, mech_data, target_handles, mechptr_by_name); - if (info) { - probe_map.insert({pi.id, {std::move(info), pi.tag}}); + if (!probe_data.empty()) { + cell_member_type probe_id{gid, i}; + probe_map.tag[probe_id] = pi.tag; + + for (auto& data: probe_data) { + probe_map.data.insert({probe_id, std::move(data)}); + } } } } @@ -653,7 +662,7 @@ 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 +// Resolution of probe addresses into a specific fvm_probe_data draws upon data // from the cable cell, the discretization, the target handle map, and the // back-end shared state. // @@ -663,6 +672,7 @@ fvm_size_type fvm_lowered_cell_impl<Backend>::fvm_intdom( template <typename Backend> struct probe_resolution_data { + std::vector<fvm_probe_data>& result; typename Backend::shared_state* state; const cable_cell& cell; const std::size_t cell_idx; @@ -701,7 +711,8 @@ struct probe_resolution_data { }; template <typename Backend> -fvm_probe_info fvm_lowered_cell_impl<Backend>::resolve_probe_address( +void fvm_lowered_cell_impl<Backend>::resolve_probe_address( + std::vector<fvm_probe_data>& probe_data, const std::vector<cable_cell>& cells, std::size_t cell_idx, const util::any& paddr, @@ -710,8 +721,9 @@ 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) { + probe_data.clear(); probe_resolution_data<Backend> prd{ - state_.get(), cells[cell_idx], cell_idx, D, M, handles, mech_instance_by_name}; + probe_data, state_.get(), cells[cell_idx], cell_idx, D, M, handles, mech_instance_by_name}; using V = util::any_visitor< cable_probe_membrane_voltage, @@ -732,25 +744,28 @@ fvm_probe_info fvm_lowered_cell_impl<Backend>::resolve_probe_address( 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{}; }); + [&prd](auto& probe_addr) { resolve_probe(probe_addr, prd); }, + [] { throw cable_cell_error("unrecognized probe type"), fvm_probe_data{}; }); return V::visit(visitor, paddr); } template <typename B> -fvm_probe_info resolve_probe(const cable_probe_membrane_voltage& p, probe_resolution_data<B>& R) { +void 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); - return fvm_probe_interpolated{ - {data+in.proximal_cv, data+in.distal_cv}, - {in.proximal_coef, in.distal_coef}, - p.location}; + for (mlocation loc: thingify(p.locations, R.cell.provider())) { + fvm_voltage_interpolant in = fvm_interpolate_voltage(R.cell, R.D, R.cell_idx, loc); + + R.result.push_back(fvm_probe_interpolated{ + {data+in.proximal_cv, data+in.distal_cv}, + {in.proximal_coef, in.distal_coef}, + loc}); + } } template <typename B> -fvm_probe_info resolve_probe(const cable_probe_membrane_voltage_cell& p, probe_resolution_data<B>& R) { +void resolve_probe(const cable_probe_membrane_voltage_cell& p, probe_resolution_data<B>& R) { fvm_probe_multi r; mcable_list cables; @@ -763,29 +778,35 @@ fvm_probe_info resolve_probe(const cable_probe_membrane_voltage_cell& p, probe_r } r.metadata = std::move(cables); r.shrink_to_fit(); - return r; + + R.result.push_back(std::move(r)); } template <typename B> -fvm_probe_info resolve_probe(const cable_probe_axial_current& p, probe_resolution_data<B>& R) { +void 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}; + for (mlocation loc: thingify(p.locations, R.cell.provider())) { + fvm_voltage_interpolant in = fvm_axial_current(R.cell, R.D, R.cell_idx, loc); + + R.result.push_back(fvm_probe_interpolated{ + {data+in.proximal_cv, data+in.distal_cv}, + {in.proximal_coef, in.distal_coef}, + loc}); + } } 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}; +void resolve_probe(const cable_probe_total_ion_current_density& p, probe_resolution_data<B>& R) { + for (mlocation loc: thingify(p.locations, R.cell.provider())) { + R.result.push_back(fvm_probe_scalar{ + {R.state->current_density.data() + R.D.geometry.location_cv(R.cell_idx, loc, cv_prefer::cv_nonempty)}, + loc}); + } } template <typename B> -fvm_probe_info resolve_probe(const cable_probe_total_ion_current_cell& p, probe_resolution_data<B>& R) { +void 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)) { @@ -800,11 +821,11 @@ fvm_probe_info resolve_probe(const cable_probe_total_ion_current_cell& p, probe_ } } r.shrink_to_fit(); - return r; + R.result.push_back(std::move(r)); } template <typename B> -fvm_probe_info resolve_probe(const cable_probe_total_current_cell& p, probe_resolution_data<B>& R) { +void 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); @@ -829,28 +850,32 @@ fvm_probe_info resolve_probe(const cable_probe_total_current_cell& p, probe_reso r.cv_cables_divs.push_back(r.metadata.size()); } r.shrink_to_fit(); - return r; + R.result.push_back(std::move(r)); } template <typename B> -fvm_probe_info resolve_probe(const cable_probe_density_state& p, probe_resolution_data<B>& R) { +void 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 {}; + if (!data) 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 {}; + auto support = R.mechanism_support(p.mechanism); + for (mlocation loc: thingify(p.locations, R.cell.provider())) { + if (!support.intersects(loc)) continue; - return fvm_probe_scalar{{data+*opt_i}, p.location}; + fvm_index_type cv = R.D.geometry.location_cv(R.cell_idx, loc, cv_prefer::cv_nonempty); + auto opt_i = util::binary_search_index(R.M.mechanisms.at(p.mechanism).cv, cv); + if (!opt_i) continue; + + R.result.push_back(fvm_probe_scalar{{data+*opt_i}, loc}); + } } template <typename B> -fvm_probe_info resolve_probe(const cable_probe_density_state_cell& p, probe_resolution_data<B>& R) { +void 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 {}; + if (!data) return; mextent support = R.mechanism_support(p.mechanism); auto& mech_cvs = R.M.mechanisms.at(p.mechanism).cv; @@ -869,22 +894,22 @@ fvm_probe_info resolve_probe(const cable_probe_density_state_cell& p, probe_reso } r.metadata = std::move(cables); r.shrink_to_fit(); - return r; + R.result.push_back(std::move(r)); } template <typename B> -fvm_probe_info resolve_probe(const cable_probe_point_state& p, probe_resolution_data<B>& R) { +void 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 {}; + 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 (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 {}; + 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; @@ -896,13 +921,13 @@ fvm_probe_info resolve_probe(const cable_probe_point_state& p, probe_resolution_ cable_probe_point_info metadata{p.target, multiplicity.empty()? 1u: multiplicity.at(mech_index), loc}; - return fvm_probe_scalar{{data+mech_index}, metadata}; + R.result.push_back(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) { +void 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 {}; + if (!data) return; unsigned id = R.mech_instance_by_name.at(p.mechanism)->mechanism_id(); auto& multiplicity = R.M.mechanisms.at(p.mechanism).multiplicity; @@ -932,20 +957,22 @@ fvm_probe_info resolve_probe(const cable_probe_point_state_cell& p, probe_resolu r.metadata = std::move(metadata); r.shrink_to_fit(); - return r; + R.result.push_back(std::move(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 {}; +void resolve_probe(const cable_probe_ion_current_density& p, probe_resolution_data<B>& R) { + for (mlocation loc: thingify(p.locations, R.cell.provider())) { + auto opt_i = R.ion_location_index(p.ion, loc); + if (!opt_i) continue; - return fvm_probe_scalar{{R.state->ion_data.at(p.ion).iX_.data()+*opt_i}, p.location}; + R.result.push_back(fvm_probe_scalar{{R.state->ion_data.at(p.ion).iX_.data()+*opt_i}, loc}); + } } 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 {}; +void 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(); @@ -966,28 +993,32 @@ fvm_probe_info resolve_probe(const cable_probe_ion_current_cell& p, probe_resolu } } r.metadata.shrink_to_fit(); - return r; + R.result.push_back(std::move(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 {}; +void resolve_probe(const cable_probe_ion_int_concentration& p, probe_resolution_data<B>& R) { + for (mlocation loc: thingify(p.locations, R.cell.provider())) { + auto opt_i = R.ion_location_index(p.ion, loc); + if (!opt_i) continue; - return fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xi_.data()+*opt_i}, p.location}; + R.result.push_back(fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xi_.data()+*opt_i}, loc}); + } } 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 {}; +void resolve_probe(const cable_probe_ion_ext_concentration& p, probe_resolution_data<B>& R) { + for (mlocation loc: thingify(p.locations, R.cell.provider())) { + auto opt_i = R.ion_location_index(p.ion, loc); + if (!opt_i) continue; - return fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xo_.data()+*opt_i}, p.location}; + R.result.push_back(fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xo_.data()+*opt_i}, loc}); + } } // 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) { +void 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; @@ -1001,19 +1032,19 @@ fvm_probe_info resolve_ion_conc_common(const std::vector<fvm_index_type>& ion_cv } r.metadata = std::move(cables); r.shrink_to_fit(); - return r; + R.result.push_back(std::move(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); +void resolve_probe(const cable_probe_ion_int_concentration_cell& p, probe_resolution_data<B>& R) { + if (!R.state->ion_data.count(p.ion)) 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); +void resolve_probe(const cable_probe_ion_ext_concentration_cell& p, probe_resolution_data<B>& R) { + if (!R.state->ion_data.count(p.ion)) 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 4d2f5d8b9057a358bbf223c1f78cae70cf5d7393..574750a60a20f14feae01411793b20e904b9a912 100644 --- a/arbor/include/arbor/cable_cell.hpp +++ b/arbor/include/arbor/cable_cell.hpp @@ -15,6 +15,7 @@ #include <arbor/morph/mprovider.hpp> #include <arbor/morph/morphology.hpp> #include <arbor/morph/primitives.hpp> +#include <arbor/util/hash_def.hpp> #include <arbor/util/typed_map.hpp> namespace arb { @@ -52,7 +53,10 @@ using cable_sample_range = std::pair<const double*, const double*>; // * `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. - +// +// Scalar probes which are described by a locset expression will generate multiple +// calls to an attached sampler, one per valid location matched by the expression. +// // Metadata for point process probes. struct cable_probe_point_info { cell_lid_type target; // Target number of point process instance on cell. @@ -64,7 +68,7 @@ struct cable_probe_point_info { // Sample value type: `double` // Sample metadata type: `mlocation` struct cable_probe_membrane_voltage { - mlocation location; + locset locations; }; // Voltage estimate [mV], reported against each cable in each control volume. Not interpolated. @@ -76,14 +80,14 @@ struct cable_probe_membrane_voltage_cell {}; // Sample value type: `double` // Sample metadata type: `mlocation` struct cable_probe_axial_current { - mlocation location; + locset locations; }; // 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; + locset locations; }; // Total ionic current [nA] across membrance _excluding_ capacitive current across components of the cell. @@ -100,7 +104,7 @@ struct cable_probe_total_current_cell {}; // Sample value type: `double` // Sample metadata type: `mlocation` struct cable_probe_density_state { - mlocation location; + locset locations; std::string mechanism; std::string state; }; @@ -135,7 +139,7 @@ struct cable_probe_point_state_cell { // Sample value type: `double` // Sample metadata type: `mlocation` struct cable_probe_ion_current_density { - mlocation location; + locset locations; std::string ion; }; @@ -150,7 +154,7 @@ struct cable_probe_ion_current_cell { // Sample value type: `double` // Sample metadata type: `mlocation` struct cable_probe_ion_int_concentration { - mlocation location; + locset locations; std::string ion; }; @@ -165,7 +169,7 @@ struct cable_probe_ion_int_concentration_cell { // Sample value type: `double` // Sample metadata type: `mlocation` struct cable_probe_ion_ext_concentration { - mlocation location; + locset locations; std::string ion; }; @@ -335,3 +339,5 @@ private: }; } // namespace arb + +ARB_DEFINE_HASH(arb::cable_probe_point_info, a.target, a.multiplicity, a.loc); diff --git a/arbor/include/arbor/common_types.hpp b/arbor/include/arbor/common_types.hpp index 11b11efd850c9bce75ec71026d59d6ee442575b9..fb7a91e8b83bf2daa64fc1453cf8fffec7f8d427 100644 --- a/arbor/include/arbor/common_types.hpp +++ b/arbor/include/arbor/common_types.hpp @@ -13,6 +13,7 @@ #include <type_traits> #include <arbor/util/lexcmp_def.hpp> +#include <arbor/util/hash_def.hpp> namespace arb { @@ -96,29 +97,4 @@ std::ostream& operator<<(std::ostream& o, backend_kind k); } // namespace arb -namespace std { - template <> struct hash<arb::cell_member_type> { - std::size_t operator()(const arb::cell_member_type& m) const { - using namespace arb; - if (sizeof(std::size_t)>sizeof(cell_gid_type)) { - constexpr unsigned shift = 8*sizeof(cell_gid_type); - - std::size_t k = m.gid; - k <<= (shift/2); // dodge gcc shift warning when other branch taken - k <<= (shift/2); - k += m.index; - return std::hash<std::size_t>{}(k); - } - else { - constexpr std::size_t prime1 = 93481; - constexpr std::size_t prime2 = 54517; - - std::size_t k = prime1; - k = k*prime2 + m.gid; - k = k*prime2 + m.index; - return k; - } - } - }; -} - +ARB_DEFINE_HASH(arb::cell_member_type, a.gid, a.index) diff --git a/arbor/include/arbor/load_balance.hpp b/arbor/include/arbor/load_balance.hpp index 8849f45c5a2e075eb16702d45b2fe92c0faf0f9b..7a021f19192e076f7b4af744648ae6b25cf54a6e 100644 --- a/arbor/include/arbor/load_balance.hpp +++ b/arbor/include/arbor/load_balance.hpp @@ -1,5 +1,8 @@ #pragma once +#include <cstddef> +#include <unordered_map> + #include <arbor/context.hpp> #include <arbor/domain_decomposition.hpp> #include <arbor/recipe.hpp> diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp index 8848bdd714b0cfdc6c852a1fe4d288c1529e6cf4..a77864fa643242b08655dc25004ade29852b8149 100644 --- a/arbor/include/arbor/morph/primitives.hpp +++ b/arbor/include/arbor/morph/primitives.hpp @@ -5,6 +5,7 @@ #include <ostream> #include <vector> +#include <arbor/util/hash_def.hpp> #include <arbor/util/lexcmp_def.hpp> // @@ -139,3 +140,6 @@ ARB_PROP(terminal) ARB_PROP(collocated) } // namespace arb + +ARB_DEFINE_HASH(arb::mcable, a.branch, a.prox_pos, a.dist_pos); +ARB_DEFINE_HASH(arb::mlocation, a.branch, a.pos); diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp index 9754408384034c3daa1ab9b9fff2169712c5574a..2d8593f0da85a4d83b398b6434050557fba7c696 100644 --- a/arbor/include/arbor/recipe.hpp +++ b/arbor/include/arbor/recipe.hpp @@ -1,11 +1,8 @@ #pragma once -#include <cstddef> -#include <memory> -#include <unordered_map> -#include <stdexcept> +#include <utility> +#include <vector> -#include <arbor/arbexcept.hpp> #include <arbor/common_types.hpp> #include <arbor/event_generator.hpp> #include <arbor/util/unique_any.hpp> @@ -13,11 +10,19 @@ namespace arb { struct probe_info { - cell_member_type id; probe_tag tag; // Address type will be specific to cell kind of cell `id.gid`. util::any address; + + probe_info(probe_info&) = default; + probe_info(const probe_info&) = default; + probe_info(probe_info&&) = default; + + // Implicit ctor uses tag of zero. + template <typename X> + probe_info(X&& x, probe_tag tag = 0): + tag(tag), address(std::forward<X>(x)) {} }; /* Recipe descriptions are cell-oriented: in order that the building @@ -66,7 +71,6 @@ public: virtual cell_size_type num_sources(cell_gid_type) const { return 0; } virtual cell_size_type num_targets(cell_gid_type) const { return 0; } - virtual cell_size_type num_probes(cell_gid_type) const { return 0; } virtual cell_size_type num_gap_junction_sites(cell_gid_type gid) const { return gap_junctions_on(gid).size(); } @@ -80,8 +84,8 @@ public: return {}; } - virtual probe_info get_probe(cell_member_type probe_id) const { - throw bad_probe_id(probe_id); + virtual std::vector<probe_info> get_probes(cell_gid_type gid) const { + return {}; } // Global property type will be specific to given cell kind. diff --git a/arbor/include/arbor/sampling.hpp b/arbor/include/arbor/sampling.hpp index b02e81ab64bbc199f1677efd660479af78b6a1c4..fa34494f0f308c579165075b9319cc4c52552c0e 100644 --- a/arbor/include/arbor/sampling.hpp +++ b/arbor/include/arbor/sampling.hpp @@ -16,15 +16,20 @@ inline cell_member_predicate one_probe(cell_member_type pid) { return [pid](cell_member_type x) { return pid==x; }; } +struct probe_metadata { + cell_member_type id; // probe id + probe_tag tag; // probe tag associated with id + unsigned index; // index of probe source within those supplied by probe id + util::any_ptr meta; // probe-specific metadata +}; + struct sample_record { time_type time; util::any_ptr data; }; 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 + void (probe_metadata, std::size_t, // number of sample records const sample_record* // pointer to first sample record )>; diff --git a/arbor/include/arbor/simple_sampler.hpp b/arbor/include/arbor/simple_sampler.hpp index a95ae8c8deaba94659505c5d7fa0f9a8b399ba91..3e98d21933647e0f6708f0335f545553908b622c 100644 --- a/arbor/include/arbor/simple_sampler.hpp +++ b/arbor/include/arbor/simple_sampler.hpp @@ -22,16 +22,76 @@ 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. +// `trace_data` wraps a std::vector of trace_entry with +// a copy of the probe-specific metadata associated with a probe. +// +// If `Meta` is void, ignore any metadta. template <typename V, typename Meta = void> -struct trace_data: public std::vector<trace_entry<V>> { - util::optional<Meta> metadata; +struct trace_data: private std::vector<trace_entry<V>> { + Meta meta; + explicit operator bool() const { return !this->empty(); } + + using base = std::vector<trace_entry<V>>; + using base::size; + using base::empty; + using base::at; + using base::begin; + using base::end; + using base::clear; + using base::resize; + using base::push_back; + using base::emplace_back; + using base::operator[]; }; template <typename V> -struct trace_data<V, void>: public std::vector<trace_entry<V>> {}; +struct trace_data<V, void>: private std::vector<trace_entry<V>> { + explicit operator bool() const { return !this->empty(); } + + using base = std::vector<trace_entry<V>>; + using base::size; + using base::empty; + using base::at; + using base::begin; + using base::end; + using base::clear; + using base::resize; + using base::push_back; + using base::emplace_back; + using base::operator[]; +}; + +// `trace_vector` wraps a vector of `trace_data`. +// +// When there are multiple probes associated with a probe id, +// the ith element will correspond to the sample data obtained +// from the probe with index i. +// +// The method `get(i)` returns a const reference to the ith +// element if it exists, or else to an empty trace_data value. + +template <typename V, typename Meta = void> +struct trace_vector: private std::vector<trace_data<V, Meta>> { + const trace_data<V, Meta>& get(std::size_t i) const { + return i<this->size()? (*this)[i]: empty_trace; + } + + using base = std::vector<trace_data<V, Meta>>; + using base::size; + using base::empty; + using base::at; + using base::begin; + using base::end; + using base::clear; + using base::resize; + using base::push_back; + using base::emplace_back; + using base::operator[]; + +private: + trace_data<V, Meta> empty_trace; +}; // Add a bit of smarts for collecting variable-length samples which are // passed back as a pair of pointers describing a range; these can @@ -68,49 +128,57 @@ struct trace_push_back<std::vector<V>> { template <typename V, typename Meta = void> class simple_sampler { public: - explicit simple_sampler(trace_data<V, Meta>& trace): trace_(trace) {} + explicit simple_sampler(trace_vector<V, Meta>& 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) { - // 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"); - } + // TODO: C++17 use if constexpr to test for Meta = void case. + void operator()(probe_metadata pm, std::size_t n, const sample_record* recs) { + const Meta* m = util::any_cast<const Meta*>(pm.meta); + if (!m) { + throw std::runtime_error("unexpected metadata type in simple_sampler"); + } + + if (trace_.size()<=pm.index) { + trace_.resize(pm.index+1); + } + + if (trace_[pm.index].empty()) { + trace_[pm.index].meta = *m; } for (std::size_t i = 0; i<n; ++i) { - if (!trace_push_back<V>::push_back(trace_, recs[i])) { + if (!trace_push_back<V>::push_back(trace_[pm.index], recs[i])) { throw std::runtime_error("unexpected sample type in simple_sampler"); } } } private: - trace_data<V, Meta>& trace_; + trace_vector<V, Meta>& trace_; }; template <typename V> class simple_sampler<V, void> { public: - explicit simple_sampler(trace_data<V, void>& trace): trace_(trace) {} + explicit simple_sampler(trace_vector<V, void>& trace): trace_(trace) {} + + void operator()(probe_metadata pm, std::size_t n, const sample_record* recs) { + if (trace_.size()<=pm.index) { + trace_.resize(pm.index+1); + } - 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])) { + if (!trace_push_back<V>::push_back(trace_[pm.index], recs[i])) { throw std::runtime_error("unexpected sample type in simple_sampler"); } } } private: - trace_data<V, void>& trace_; + trace_vector<V, void>& trace_; }; template <typename V, typename Meta> -inline simple_sampler<V, Meta> make_simple_sampler(trace_data<V, Meta>& trace) { +inline simple_sampler<V, Meta> make_simple_sampler(trace_vector<V, Meta>& trace) { return simple_sampler<V, Meta>(trace); } diff --git a/arbor/include/arbor/symmetric_recipe.hpp b/arbor/include/arbor/symmetric_recipe.hpp index 70f1568315d764e5f315ec494255c97799d14fbc..917ec892c022d3857e8d05d8dc20bcdc66e54abc 100644 --- a/arbor/include/arbor/symmetric_recipe.hpp +++ b/arbor/include/arbor/symmetric_recipe.hpp @@ -47,10 +47,6 @@ public: return tiled_recipe_->num_targets(i % tiled_recipe_->num_cells()); } - cell_size_type num_probes(cell_gid_type i) const override { - return tiled_recipe_->num_probes(i % tiled_recipe_->num_cells()); - } - // Only function that calls the underlying tile's function on the same gid. // This is because applying transformations to event generators is not straightforward. std::vector<event_generator> event_generators(cell_gid_type i) const override { @@ -73,9 +69,9 @@ public: return conns; } - probe_info get_probe(cell_member_type probe_id) const override { - probe_id.gid %= tiled_recipe_->num_cells(); - return tiled_recipe_->get_probe(probe_id); + std::vector<probe_info> get_probes(cell_gid_type i) const override { + i %= tiled_recipe_->num_cells(); + return tiled_recipe_->get_probes(i); } util::any get_global_properties(cell_kind ck) const override { diff --git a/arbor/include/arbor/util/hash_def.hpp b/arbor/include/arbor/util/hash_def.hpp new file mode 100644 index 0000000000000000000000000000000000000000..fc6770cc4125b30aa7f1282547c25f2a0e3a220d --- /dev/null +++ b/arbor/include/arbor/util/hash_def.hpp @@ -0,0 +1,46 @@ +#pragma once + +/* + * Macro definitions for defining hash functions for compound objects. + * + * Use: + * + * To define a std::hash overload for a record type xyzzy + * with fields foo, bar and baz: + * + * ARB_DEFINE_HASH(xyzzy, a.foo, a.bar, a.baz) + * + * The explicit use of 'a' in the macro invocation is just to simplify + * the implementation. + * + * The macro must be used outside of any namespace. + */ + +#include <cstddef> +#include <typeindex> + +// Helpers for forming hash values of compounds objects. + +namespace arb { + +inline std::size_t hash_value_combine(std::size_t n) { + return n; +} + +template <typename... T> +std::size_t hash_value_combine(std::size_t n, std::size_t m, T... tail) { + constexpr std::size_t prime2 = 54517; + return hash_value_combine(prime2*n + m, tail...); +} + +template <typename... T> +std::size_t hash_value(const T&... t) { + constexpr std::size_t prime1 = 93481; + return hash_value_combine(prime1, std::hash<T>{}(t)...); +} +} + +#define ARB_DEFINE_HASH(type,...)\ +namespace std {\ +template <> struct hash<type> { std::size_t operator()(const type& a) const { return ::arb::hash_value(__VA_ARGS__); }};\ +} diff --git a/arbor/include/arbor/util/lexcmp_def.hpp b/arbor/include/arbor/util/lexcmp_def.hpp index a6f0760c0d26b65279dbc79b4c8838688fdd11e1..57f0062ad7524f975f67373ef94eeb78e3f6eac8 100644 --- a/arbor/include/arbor/util/lexcmp_def.hpp +++ b/arbor/include/arbor/util/lexcmp_def.hpp @@ -9,7 +9,7 @@ * To define comparison operations for a record type xyzzy * with fields foo, bar and baz: * - * DEFINE_LEXICOGRAPHIC_ORDERING(xyzzy,(a.foo,a.bar,a.baz),(b.foo,b.bar,b.baz)) + * ARB_DEFINE_LEXICOGRAPHIC_ORDERING(xyzzy,(a.foo,a.bar,a.baz),(b.foo,b.bar,b.baz)) * * The explicit use of 'a' and 'b' in the second and third parameters * is needed only to save a heroic amount of preprocessor macro diff --git a/arbor/mc_cell_group.cpp b/arbor/mc_cell_group.cpp index 454069358894ba9f547ee4f6d50c6cc8259aea01..f5d85ff572bac32d9c3672726828d560b74ffe23 100644 --- a/arbor/mc_cell_group.cpp +++ b/arbor/mc_cell_group.cpp @@ -43,9 +43,7 @@ mc_cell_group::mc_cell_group(const std::vector<cell_gid_type>& gids, const recip util::transform_view(gids_, [&rec](cell_gid_type i) { return rec.num_targets(i); })); std::size_t n_targets = target_handle_divisions_.back(); - // Pre-allocate space to store handles, probe map. - auto n_probes = util::sum_by(gids_, [&rec](cell_gid_type i) { return rec.num_probes(i); }); - probe_map_.reserve(n_probes); + // Pre-allocate space to store handles. target_handles_.reserve(n_targets); // Construct cell implementation, retrieving handles and maps. @@ -86,6 +84,8 @@ struct sampler_call_info { sampler_function sampler; cell_member_type probe_id; probe_tag tag; + unsigned index; + const fvm_probe_data* pdata_ptr; // Offsets are into lowered cell sample time and event arrays. sample_size_type begin_offset; @@ -112,7 +112,7 @@ void run_samples( std::vector<sample_record>&, fvm_probe_scratch&) { - throw arbor_internal_error("invalid fvm_probe_info in sampler map"); + throw arbor_internal_error("invalid fvm_probe_data in sampler map"); } void run_samples( @@ -135,7 +135,7 @@ void run_samples( // 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()); }, + [&](auto& metadata) { sc.sampler({sc.probe_id, sc.tag, sc.index, &metadata}, n_sample, sample_records.data()); }, p.metadata); } @@ -166,7 +166,7 @@ void run_samples( 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()); + sc.sampler({sc.probe_id, sc.tag, sc.index, &p.metadata}, n_sample, sample_records.data()); } void run_samples( @@ -198,7 +198,7 @@ void run_samples( // 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()); }, + [&](auto& metadata) { sc.sampler({sc.probe_id, sc.tag, sc.index, &metadata}, n_sample, sample_records.data()); }, p.metadata); } @@ -243,8 +243,9 @@ void run_samples( } // (Unlike fvm_probe_multi, we only have mcable_list metadata.) - sc.sampler(sc.probe_id, sc.tag, &p.metadata, n_sample, sample_records.data()); + sc.sampler({sc.probe_id, sc.tag, sc.index, &p.metadata}, n_sample, sample_records.data()); } + void run_samples( const fvm_probe_membrane_currents& p, const sampler_call_info& sc, @@ -304,19 +305,18 @@ void run_samples( 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()); + sc.sampler({sc.probe_id, sc.tag, sc.index, &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); + util::visit([&](auto& x) {run_samples(x, sc, raw_times, raw_samples, sample_records, scratch); }, sc.pdata_ptr->info); } void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& event_lanes) { @@ -402,20 +402,22 @@ 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*n_raw}); - for (auto t: sample_times) { + probe_tag tag = probe_map_.tag.at(pid); + unsigned index = 0; + for (const fvm_probe_data& pdata: probe_map_.data_on(pid)) { + call_info.push_back({sa.sampler, pid, tag, index++, &pdata, n_samples, n_samples + n_times*pdata.n_raw()}); auto intdom = cell_to_intdom_[cell_index]; - for (probe_handle h: p.handle.raw_handle_range()) { - sample_event ev{t, (cell_gid_type)intdom, {h, n_samples++}}; - sample_events.push_back(ev); - } - if (sa.policy==sampling_policy::exact) { - target_handle h(-1, 0, intdom); - exact_sampling_events.push_back({t, h, 0.f}); + + for (auto t: sample_times) { + for (probe_handle h: pdata.raw_handle_range()) { + sample_event ev{t, (cell_gid_type)intdom, {h, n_samples++}}; + sample_events.push_back(ev); + } + if (sa.policy==sampling_policy::exact) { + target_handle h(-1, 0, intdom); + exact_sampling_events.push_back({t, h, 0.f}); + } } } } @@ -462,8 +464,7 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e reserve_scratch(scratch, max_samples_per_call); 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); + run_samples(sc, result.sample_time.data(), result.sample_value.data(), sample_records, scratch); } PL(); @@ -481,7 +482,7 @@ void mc_cell_group::add_sampler(sampler_association_handle h, cell_member_predic schedule sched, sampler_function fn, sampling_policy policy) { std::vector<cell_member_type> probeset = - util::assign_from(util::filter(util::keys(probe_map_), probe_ids)); + util::assign_from(util::filter(util::keys(probe_map_.tag), probe_ids)); if (!probeset.empty()) { sampler_map_.add(h, sampler_association{std::move(sched), std::move(fn), std::move(probeset), policy}); diff --git a/arbor/mc_cell_group.hpp b/arbor/mc_cell_group.hpp index d11da380d31abd77ede679d52f726b0c26ba3567..8285a4f0648f1c5239b98cf12ae32e9f9df00821 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<fvm_probe_info> probe_map_; + probe_association_map probe_map_; // Collection of samplers to be run against probes in this group. sampler_association_map sampler_map_; diff --git a/arbor/partition_load_balance.cpp b/arbor/partition_load_balance.cpp index 87e1b4043e139947465a09500f714978f89f3118..5236658e2ba27cb3bd098ff534344af0ddc69b6f 100644 --- a/arbor/partition_load_balance.cpp +++ b/arbor/partition_load_balance.cpp @@ -1,7 +1,9 @@ #include <queue> +#include <unordered_map> #include <unordered_set> #include <vector> +#include <arbor/arbexcept.hpp> #include <arbor/domain_decomposition.hpp> #include <arbor/load_balance.hpp> #include <arbor/recipe.hpp> diff --git a/arbor/sampler_map.hpp b/arbor/sampler_map.hpp index 0366885d478bb2c4685ba9e309b96973502f05db..50d0ff9b4a00df551752a65d2810a206e5fb2659 100644 --- a/arbor/sampler_map.hpp +++ b/arbor/sampler_map.hpp @@ -61,16 +61,4 @@ public: auto end() { return assoc_view().end(); } }; -// Manage associations between probe ids, probe tags, and (lowered cell) probe handles. - -template <typename Handle> -struct probe_association { - using probe_handle_type = Handle; - probe_handle_type handle; - probe_tag tag; -}; - -template <typename Handle> -using probe_association_map = std::unordered_map<cell_member_type, probe_association<Handle>>; - } // namespace arb diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index 2e87f4a45cf4da9f19f52acfd0a158c2d6a4681b..7e32ec774cc7218f03c68849fd714d5e94debcf3 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -2,6 +2,7 @@ #include <set> #include <vector> +#include <arbor/arbexcept.hpp> #include <arbor/context.hpp> #include <arbor/domain_decomposition.hpp> #include <arbor/generic_event.hpp> diff --git a/arbor/spike_source_cell_group.cpp b/arbor/spike_source_cell_group.cpp index 5b309e40348d8d8dae5db5d48d7890f419b91c5f..c5521b306f14ba6cd4f1cbb4f936a74b665d8853 100644 --- a/arbor/spike_source_cell_group.cpp +++ b/arbor/spike_source_cell_group.cpp @@ -1,5 +1,6 @@ #include <exception> +#include <arbor/arbexcept.hpp> #include <arbor/recipe.hpp> #include <arbor/spike_source_cell.hpp> #include <arbor/schedule.hpp> diff --git a/doc/cpp_cable_cell.rst b/doc/cpp_cable_cell.rst index 8418546d4cc9ffc2b77767f29b642c1494ded28d..02f13b3ed3335ae4a2d29ea032630f711ad7b22c 100644 --- a/doc/cpp_cable_cell.rst +++ b/doc/cpp_cable_cell.rst @@ -308,20 +308,26 @@ 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. +Cable cell probe addresses that are described by a ``locset`` may generate more +than one concrete probe: there will be one per location in the locset that is +satisfiable. Sampler callback functions can distinguish between different +probes with the same address and id by examining their index and/or +probe-sepcific metadata found in the ``probe_metadata`` parameter. + Membrane voltage ^^^^^^^^^^^^^^^^ .. code:: struct cable_probe_membrane_voltage { - mlocation location; + locset locations; }; -Queries cell membrane potential at the specified location. +Queries cell membrane potential at each site in ``locations``. * Sample value: ``double``. Membrane potential in millivolts. -* Metadata: ``mlocation``. Location as given in the probe address. +* Metadata: ``mlocation``. Location of probe. .. code:: @@ -343,14 +349,15 @@ Axial current .. code:: struct cable_probe_axial_current { - mlocation location; + locset locations; }; -Estimate intracellular current at given location in the distal direction. +Estimate intracellular current at each site in ``locations``, +in the distal direction. * Sample value: ``double``. Current in nanoamperes. -* Metadata: ``mlocation``. Location as given in the probe address. +* Metadata: ``mlocation``. Location as of probe. Transmembrane current @@ -359,15 +366,16 @@ Transmembrane current .. code:: struct cable_probe_ion_current_density { - mlocation location; + locset locations; std::string ion; }; -Membrance current density attributed to a particular ion at a given location. +Membrance current density attributed to a particular ion at +each site in ``locations``. * Sample value: ``double``. Current density in amperes per square metre. -* Metadata: ``mlocation``. Location as given in the probe address. +* Metadata: ``mlocation``. Location of probe. .. code:: @@ -389,14 +397,14 @@ Membrane current attributed to a particular ion across components of the cell. .. code:: struct cable_probe_total_ion_current_density { - mlocation location; + locset locations; }; -Membrane current density at gvien location _excluding_ capacitive currents. +Membrane current density at given locations _excluding_ capacitive currents. * Sample value: ``double``. Current density in amperes per square metre. -* Metadata: ``mlocation``. Location as given in the probe address. +* Metadata: ``mlocation``. Location of probe. .. code:: @@ -433,15 +441,15 @@ Ion concentration .. code:: struct cable_probe_ion_int_concentration { - mlocation location; + locset locations; std::string ion; }; -Ionic internal concentration of ion at given location. +Ionic internal concentration of ion at each site in ``locations``. * Sample value: ``double``. Ion concentration in millimoles per litre. -* Metadata: ``mlocation``. Location as given in the probe address. +* Metadata: ``mlocation``. Location of probe. .. code:: @@ -467,11 +475,11 @@ Ionic external concentration of ion across components of the cell. std::string ion; }; -Ionic internal concentration of ion at given location. +Ionic external concentration of ion at each site in ``locations``. * Sample value: ``double``. Ion concentration in millimoles per litre. -* Metadata: ``mlocation``. Location as given in the probe address. +* Metadata: ``mlocation``. Location of probe. .. code:: @@ -497,14 +505,14 @@ Mechanism state .. code:: struct cable_probe_density_state { - mlocation location; + locset locations; 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. +Value of state variable in a density mechanism in each site in ``locations``. +If the mechanism is not defined at a particular site, that site is ignored. * Sample value: ``double``. State variable value. diff --git a/doc/sampling_api.rst b/doc/sampling_api.rst index 95e596f4852bc8918492e58959fdb0b68334c654..4deac7fe3d848a9c3fd5a93f3cfd42e9171a9c7c 100644 --- a/doc/sampling_api.rst +++ b/doc/sampling_api.rst @@ -39,32 +39,36 @@ probe will be specific to a particular cell type. using probe_tag = int; struct probe_info { - cell_member_type id; // cell gid, index of probe probe_tag tag; // opaque key, returned in sample record any address; // cell-type specific location info - }; - - probe_info recipe::get_probe(cell_member_type probe_id); + template <typename X> + probe_info(X&& a, probe_tag tag = 0): + tag(tag), address(std::forward<X>(x)) {} + }; -The ``id`` field in the ``probe_info`` struct will be the same value as -the ``probe_id`` used in the ``get_probe`` call. + std::vector<probe_info> recipe::get_probes(cell_gid_type gid); -The ``get_probe()`` method is intended for use by cell group -implementations to set up sampling data structures ahead of time and for -putting in place any structures or information in the concrete cell -implementations to allow monitoring. The ``tag`` field has no semantics for the engine. It is provided merely as a way of passing additional metadata about a probe to any sampler that polls it, with a view to samplers that handle multiple probes, possibly with different value types. -Probe addresses are now decoupled from the cell descriptions themselves — +Probe addresses are decoupled from the cell descriptions themselves — this allows a recipe implementation to construct probes independently of the cells themselves. It is the responsibility of a cell group implementation to parse the probe address objects wrapped in the ``any address`` field. +The _k_th element of the vector returned by ``get_probes(gid)`` is +identified with a probe-id: ``cell_member_type{gid, k}``. + +One probe address may describe more than one concrete probe, depending +upon the interpretation of the probe address by the cell group. In this +instance, each of the concerete probes will be associated with the +same probe-id. Samplers can distinguish between different probes with +the same id by their probe index (see below). + Samplers and sample records --------------------------- @@ -76,18 +80,29 @@ will be passed to a sampler function or function object: .. code-block:: cpp + struct probe_metadata { + cell_member_type id; // probe id + probe_tag tag; // probe tag associated with id + unsigned index; // index of probe source within those supplied by probe id + util::any_ptr meta; // probe-specific metadata + }; + using sampler_function = - std::function<void (cell_member_type, probe_tag, any_ptr, size_t, const sample_record*)>; + std::function<void (probe_metadta, size_t, const sample_record*)>; + +where the parameters are respectively the probe metadata, the number of +samples, and finally 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_id``, identifies the probe by its probe-id (see above). -The ``probe_tag`` is the key given in the ``probe_info`` returned by -the recipe. +The ``probe_tag`` in the metadata 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 +The ``index`` identifies which of the possibly multiple probes associated +with the probe-id is the source of the samples. + +The ``any_ptr`` value iin the metadata 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 @@ -109,7 +124,8 @@ probe address. The data pointed to by ``data``, and the sample records themselves, are only guaranteed to be valid for the duration of the call to the sampler -function. A simple sampler implementation for ``double`` data might be: +function. A simple sampler implementation for ``double`` data, assuming +one probe per probe id, might be as follows: .. container:: example-code @@ -122,13 +138,13 @@ 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, any_ptr, size_t n, const sample_record* records) { + void operator()(probe_metadata pm, size_t n, const sample_record* records) { for (size_t i=0; i<n; ++i) { const auto& rec = records[i]; const double* data = any_cast<const double*>(rec.data); assert(data); - samples[id].emplace_back(rec.time, *data); + samples[pm.id].emplace_back(rec.time, *data); } } }; @@ -288,10 +304,6 @@ and probe metdata. handles and tuples (*schedule*, *sampler*, *probe set*, *policy*), with thread-safe accessors. -``probe_association_map<Handle>`` is a type alias for an unordered map between -probe ids and tuples (*probe handle*, *probe tag*), where the *probe handle* -is a cell-group specific accessor that allows efficient polling. - Batched sampling in ``mc_cell_group`` ------------------------------------- @@ -313,7 +325,7 @@ Given an association tuple (*schedule*, *sampler*, *probe set*, *policy*) where has (non-zero) *n* sample times in the current integration interval, the ``mc_cell_group`` will call the *sampler* callback once for probe in *probe set*, with *n* sample values. -In addition to the ``lax`` sampling polocy, ``mc_cell_group`` supports the ``exact`` +In addition to the ``lax`` sampling policy, ``mc_cell_group`` supports the ``exact`` policy. Integration steps will be shortened such that any sample times associated with an ``exact`` policy can be satisfied precisely. diff --git a/example/brunel/brunel.cpp b/example/brunel/brunel.cpp index 1c37b7b8486fcc2535e97c63b67a08d0ec3abde2..36c74da9cb8a31f805874afa9d69e23057809e75 100644 --- a/example/brunel/brunel.cpp +++ b/example/brunel/brunel.cpp @@ -168,14 +168,6 @@ public: return 1; } - cell_size_type num_probes(cell_gid_type) const override { - return 0; - } - - probe_info get_probe(cell_member_type probe_id) const override { - return {}; - } - private: // Number of excitatory cells. cell_size_type ncells_exc_; diff --git a/example/dryrun/dryrun.cpp b/example/dryrun/dryrun.cpp index e40891c73a98c55624374d94822dc6fc60a3f7d2..53b5ce902cef34d8b87288eab7c00f8432fdc1dd 100644 --- a/example/dryrun/dryrun.cpp +++ b/example/dryrun/dryrun.cpp @@ -17,7 +17,6 @@ #include <arbor/morph/primitives.hpp> #include <arbor/profile/meter_manager.hpp> #include <arbor/profile/profiler.hpp> -#include <arbor/simple_sampler.hpp> #include <arbor/simulation.hpp> #include <arbor/symmetric_recipe.hpp> #include <arbor/recipe.hpp> @@ -120,16 +119,9 @@ public: return gens; } - // There is one probe (for measuring voltage at the soma) on the cell. - cell_size_type num_probes(cell_gid_type gid) const override { - return 1; - } - - 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::cable_probe_membrane_voltage{loc}}; + std::vector<arb::probe_info> get_probes(cell_gid_type gid) const override { + // One probe per cell, sampling membrane voltage at end of soma. + return {arb::cable_probe_membrane_voltage{arb::mlocation{0, 0.0}}}; } private: diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp index 4f3001fcf6dcf1fb4af705b1695bb6dfb5e96b92..a2db630be5aef7d7ccee9c9c1764fff3ca070b49 100644 --- a/example/gap_junctions/gap_junctions.cpp +++ b/example/gap_junctions/gap_junctions.cpp @@ -55,7 +55,7 @@ using arb::cell_kind; using arb::time_type; // Writes voltage trace as a json file. -void write_trace_json(const std::vector<arb::trace_data<double>>& trace, unsigned rank); +void write_trace_json(const std::vector<arb::trace_vector<double>>& trace, unsigned rank); // Generate a cell. arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncells, double stim_duration); @@ -93,15 +93,10 @@ public: return {arb::cell_connection({gid - 1, 0}, {gid, 0}, params_.event_weight, params_.event_min_delay)}; } - // There is one probe (for measuring voltage at the soma) on the cell. - cell_size_type num_probes(cell_gid_type gid) const override { - return 1; - } - - arb::probe_info get_probe(cell_member_type id) const override { - // Measure at the soma. + std::vector<arb::probe_info> get_probes(cell_gid_type gid) const override { + // Measure membrane voltage at end of soma. arb::mlocation loc{0, 1.}; - return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}}; + return {arb::cable_probe_membrane_voltage{loc}}; } arb::util::any get_global_properties(cell_kind k) const override { @@ -193,14 +188,13 @@ int main(int argc, char** argv) { auto sched = arb::regular_schedule(0.025); // This is where the voltage samples will be stored as (time, value) pairs - std::vector<arb::trace_data<double>> voltage(decomp.num_local_cells); + std::vector<arb::trace_vector<double>> voltage_traces(decomp.num_local_cells); // Now attach the sampler at probe_id, with sampling schedule sched, writing to voltage unsigned j=0; for (auto g : decomp.groups) { for (auto i : g.gids) { - auto t = recipe.get_probe({i, 0}); - sim.add_sampler(arb::one_probe(t.id), sched, arb::make_simple_sampler(voltage[j++])); + sim.add_sampler(arb::one_probe({i, 0}), sched, arb::make_simple_sampler(voltage_traces[j++])); } } @@ -244,7 +238,7 @@ int main(int argc, char** argv) { // Write the samples to a json file. if (params.print_all) { - write_trace_json(voltage, arb::rank(context)); + write_trace_json(voltage_traces, arb::rank(context)); } auto report = arb::profile::make_meter_report(meters, context); @@ -258,8 +252,8 @@ int main(int argc, char** argv) { return 0; } -void write_trace_json(const std::vector<arb::trace_data<double>>& trace, unsigned rank) { - for (unsigned i = 0; i < trace.size(); i++) { +void write_trace_json(const std::vector<arb::trace_vector<double>>& traces, unsigned rank) { + for (unsigned i = 0; i < traces.size(); i++) { std::string path = "./voltages_" + std::to_string(rank) + "_" + std::to_string(i) + ".json"; @@ -273,7 +267,7 @@ void write_trace_json(const std::vector<arb::trace_data<double>>& trace, unsigne auto &jt = json["data"]["time"]; auto &jy = json["data"]["voltage"]; - for (const auto &sample: trace[i]) { + for (const auto &sample: traces[i].at(0)) { jt.push_back(sample.t); jy.push_back(sample.v); } diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp index 9b62067f37fc22d2beeff2ea5380da471f1e5b57..5a0aa4d23f1bf40cdd4d9191162b84d01b8f8e79 100644 --- a/example/generators/generators.cpp +++ b/example/generators/generators.cpp @@ -115,19 +115,12 @@ public: return gens; } - // There is one probe (for measuring voltage at the soma) on the cell - cell_size_type num_probes(cell_gid_type gid) const override { - assert(gid==0); // There is only one cell in the model - return 1; - } - - arb::probe_info get_probe(cell_member_type id) const override { - assert(id.gid==0); // There is one cell, - assert(id.index==0); // with one probe. + std::vector<arb::probe_info> get_probes(cell_gid_type gid) const override { + assert(gid==0); // There is only one cell, - // Measure at the soma + // Measure membrane voltage at end of soma. arb::mlocation loc{0, 0.0}; - return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}}; + return {arb::cable_probe_membrane_voltage{loc}}; } }; @@ -153,7 +146,7 @@ int main() { // The schedule for sampling is 10 samples every 1 ms. auto sched = arb::regular_schedule(0.1); // This is where the voltage samples will be stored as (time, value) pairs - arb::trace_data<double> voltage; + arb::trace_vector<double> voltage; // Now attach the sampler at probe_id, with sampling schedule sched, writing to voltage sim.add_sampler(arb::one_probe(probe_id), sched, arb::make_simple_sampler(voltage)); @@ -161,7 +154,7 @@ int main() { sim.run(100, 0.01); // Write the samples to a json file. - write_trace_json(voltage); + write_trace_json(voltage.at(0)); } void write_trace_json(const arb::trace_data<double>& trace) { diff --git a/example/probe-demo/probe-demo.cpp b/example/probe-demo/probe-demo.cpp index 832dc4158a745e578e1da4ab4690ddb47fc681fb..895dc659995a2e35ba2f228cd4c7cc036c848578 100644 --- a/example/probe-demo/probe-demo.cpp +++ b/example/probe-demo/probe-demo.cpp @@ -83,8 +83,8 @@ struct options { 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*); +void vector_sampler(arb::probe_metadata, std::size_t, const arb::sample_record*); +void scalar_sampler(arb::probe_metadata, std::size_t, const arb::sample_record*); struct cable_recipe: public arb::recipe { arb::cable_cell_global_properties gprop; @@ -98,11 +98,10 @@ struct cable_recipe: public arb::recipe { } 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) + std::vector<arb::probe_info> get_probes(arb::cell_gid_type) const override { + return {probe_addr}; // (use default tag value 0) } arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { @@ -160,8 +159,8 @@ int main(int argc, char** argv) { } } -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); +void scalar_sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { + auto* loc = any_cast<const arb::mlocation*>(pm.meta); assert(loc); std::cout << std::fixed << std::setprecision(4); @@ -173,8 +172,8 @@ void scalar_sampler(arb::cell_member_type, arb::probe_tag, any_ptr meta, std::si } } -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); +void vector_sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { + auto* cables_ptr = any_cast<const arb::mcable_list*>(pm.meta); assert(cables_ptr); unsigned n_cable = cables_ptr->size(); @@ -200,19 +199,21 @@ bool parse_options(options& opt, int& argc, char** argv) { auto do_help = [&]() { usage(argv0, help_msg); }; + using L = arb::mlocation; + // 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"}; }}}, + {"v", {"v", true, [](double x) { return arb::cable_probe_membrane_voltage{L{0, x}}; }}}, + {"i_axial", {"i_axial", true, [](double x) { return arb::cable_probe_axial_current{L{0, x}}; }}}, + {"j_ion", {"j_ion", true, [](double x) { return arb::cable_probe_total_ion_current_density{L{0, x}}; }}}, + {"j_na", {"j_na", true, [](double x) { return arb::cable_probe_ion_current_density{L{0, x}, "na"}; }}}, + {"j_k", {"j_k", true, [](double x) { return arb::cable_probe_ion_current_density{L{0, x}, "k"}; }}}, + {"c_na", {"c_na", true, [](double x) { return arb::cable_probe_ion_int_concentration{L{0, x}, "na"}; }}}, + {"c_k", {"c_k", true, [](double x) { return arb::cable_probe_ion_int_concentration{L{0, x}, "k"}; }}}, + {"hh_m", {"hh_m", true, [](double x) { return arb::cable_probe_density_state{L{0, x}, "hh", "m"}; }}}, + {"hh_h", {"hh_h", true, [](double x) { return arb::cable_probe_density_state{L{0, x}, "hh", "h"}; }}}, + {"hh_n", {"hh_n", true, [](double x) { return arb::cable_probe_density_state{L{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{}; }}}, diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp index 8f7acde008655bdf2f4aaca5f68e875fa77ca7d8..00a427f7e260da11af11714fc222f805831c4960 100644 --- a/example/ring/ring.cpp +++ b/example/ring/ring.cpp @@ -3,6 +3,7 @@ * */ +#include <cassert> #include <fstream> #include <iomanip> #include <iostream> @@ -108,15 +109,10 @@ public: return gens; } - // There is one probe (for measuring voltage at the soma) on the cell. - cell_size_type num_probes(cell_gid_type gid) const override { - return 1; - } - - arb::probe_info get_probe(cell_member_type id) const override { - // Measure at the soma. + std::vector<arb::probe_info> get_probes(cell_gid_type gid) const override { + // Measure membrane voltage at end of soma. arb::mlocation loc{0, 0.0}; - return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}}; + return {arb::cable_probe_membrane_voltage{loc}}; } arb::util::any get_global_properties(arb::cell_kind) const override { @@ -186,7 +182,7 @@ int main(int argc, char** argv) { // The schedule for sampling is 10 samples every 1 ms. auto sched = arb::regular_schedule(1); // This is where the voltage samples will be stored as (time, value) pairs - arb::trace_data<double> voltage; + arb::trace_vector<double> voltage; // Now attach the sampler at probe_id, with sampling schedule sched, writing to voltage sim.add_sampler(arb::one_probe(probe_id), sched, arb::make_simple_sampler(voltage)); @@ -229,7 +225,9 @@ int main(int argc, char** argv) { } // Write the samples to a json file. - if (root) write_trace_json(voltage); + if (root) { + write_trace_json(voltage.at(0)); + } auto report = arb::profile::make_meter_report(meters, context); std::cout << report; diff --git a/example/single/single.cpp b/example/single/single.cpp index 07d3cdac7bd74882cc55cff1887f9f2bbc6fdd4c..ac8e67e52aea49a5eed8741136b54c7c7a04ff58 100644 --- a/example/single/single.cpp +++ b/example/single/single.cpp @@ -34,17 +34,16 @@ struct single_recipe: public arb::recipe { } 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 1; } - arb::probe_info get_probe(arb::cell_member_type probe_id) const override { + std::vector<arb::probe_info> get_probes(arb::cell_gid_type) const override { arb::mlocation mid_soma = {0, 0.5}; 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. + // Probe info consists of a probe address and an optional tag, for use + // by the sampler callback. - return {probe_id, 0, probe}; + return {arb::probe_info{probe, 0}}; } arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { @@ -89,8 +88,8 @@ int main(int argc, char** argv) { // Attach a sampler to the probe described in the recipe, sampling every 0.1 ms. - arb::trace_data<double> trace; - sim.add_sampler(arb::all_probes, arb::regular_schedule(0.1), arb::make_simple_sampler(trace)); + arb::trace_vector<double> traces; + sim.add_sampler(arb::all_probes, arb::regular_schedule(0.1), arb::make_simple_sampler(traces)); // Trigger the single synapse (target is gid 0, index 0) at t = 1 ms with // the given weight. @@ -100,8 +99,7 @@ int main(int argc, char** argv) { sim.run(opt.t_end, opt.dt); - std::cout << std::fixed << std::setprecision(4); - for (auto entry: trace) { + for (auto entry: traces.at(0)) { std::cout << entry.t << ", " << entry.v << "\n"; } } diff --git a/python/recipe.cpp b/python/recipe.cpp index 1497307d380c5a08ee44d2b46a0683e6b87de13c..45c5cfa3b2bf34a2dfca9e9da58dba816d0ad306 100644 --- a/python/recipe.cpp +++ b/python/recipe.cpp @@ -29,12 +29,12 @@ arb::util::unique_any py_recipe_shim::get_cell_description(arb::cell_gid_type gi "Python error already thrown"); } -arb::probe_info cable_probe(std::string kind, arb::cell_member_type id, arb::mlocation loc) { +arb::probe_info cable_loc_probe(std::string kind, arb::cell_member_type id, arb::mlocation loc) { if (kind == "voltage") { - return arb::probe_info{id, 0, arb::cable_probe_membrane_voltage{loc}}; + return arb::cable_probe_membrane_voltage{loc}; } else if (kind == "ionic current density") { - return arb::probe_info{id, 0, arb::cable_probe_total_ion_current_density{loc}}; + return arb::cable_probe_total_ion_current_density{loc}; } else throw pyarb_error(util::pprintf("unrecognized probe kind: {}", kind)); }; @@ -161,25 +161,22 @@ void register_recipe(pybind11::module& m) { .def("gap_junctions_on", &py_recipe::gap_junctions_on, "gid"_a, "A list of the gap junctions connected to gid, [] by default.") - .def("num_probes", &py_recipe::num_probes, + .def("get_probes", &py_recipe::get_probes, "gid"_a, - "The number of probes on gid, 0 by default.") - .def("get_probe", &py_recipe::get_probe, - "id"_a, - "The probe(s) to allow monitoring, must be provided if num_probes() returns a non-zero value.") + "The probes to allow monitoring.") // TODO: py_recipe::global_properties .def("__str__", [](const py_recipe&){return "<arbor.recipe>";}) .def("__repr__", [](const py_recipe&){return "<arbor.recipe>";}); // Probes - m.def("cable_probe", &cable_probe, + m.def("cable_probe", &cable_loc_probe, "Description of a probe at a location on a cable cell with id available for monitoring data of kind "\ "where kind is one of 'voltage' or 'ionic current density'.", "kind"_a, "id"_a, "location"_a); pybind11::class_<arb::probe_info> probe(m, "probe"); probe - .def("__repr__", [](const arb::probe_info& p){return util::pprintf("<arbor.probe: cell {}, probe {}>", p.id.gid, p.id.index);}) - .def("__str__", [](const arb::probe_info& p){return util::pprintf("<arbor.probe: cell {}, probe {}>", p.id.gid, p.id.index);}); + .def("__repr__", [](const arb::probe_info& p){return util::pprintf("<arbor.probe: tag {}>", p.tag);}) + .def("__str__", [](const arb::probe_info& p){return util::pprintf("<arbor.probe: tag {}>", p.tag);}); } } // namespace pyarb diff --git a/python/recipe.hpp b/python/recipe.hpp index d2e177688d94dbaeb9d2c76acc577447d4ca4ec1..51aa070ef4cafe07ff74881ed737021acd94c48e 100644 --- a/python/recipe.hpp +++ b/python/recipe.hpp @@ -53,11 +53,8 @@ public: virtual std::vector<arb::gap_junction_connection> gap_junctions_on(arb::cell_gid_type) const { return {}; } - virtual arb::cell_size_type num_probes(arb::cell_gid_type) const { - return 0; - } - virtual arb::probe_info get_probe (arb::cell_member_type id) const { - throw pyarb_error(util::pprintf("bad probe id {}", id)); + virtual std::vector<arb::probe_info> get_probes(arb::cell_gid_type gid) const { + return {}; } //TODO: virtual pybind11::object global_properties(arb::cell_kind kind) const {return pybind11::none();}; }; @@ -100,12 +97,8 @@ public: PYBIND11_OVERLOAD(std::vector<arb::gap_junction_connection>, py_recipe, gap_junctions_on, gid); } - arb::cell_size_type num_probes(arb::cell_gid_type gid) const override { - PYBIND11_OVERLOAD(arb::cell_size_type, py_recipe, num_probes, gid); - } - - arb::probe_info get_probe(arb::cell_member_type id) const override { - PYBIND11_OVERLOAD(arb::probe_info, py_recipe, get_probe, id); + std::vector<arb::probe_info> get_probes(arb::cell_gid_type gid) const override { + PYBIND11_OVERLOAD(std::vector<arb::probe_info>, py_recipe, get_probes, gid); } }; @@ -160,12 +153,8 @@ public: return try_catch_pyexception([&](){ return impl_->gap_junctions_on(gid); }, msg); } - arb::cell_size_type num_probes(arb::cell_gid_type gid) const override { - return try_catch_pyexception([&](){ return impl_->num_probes(gid); }, msg); - } - - arb::probe_info get_probe(arb::cell_member_type id) const override { - return try_catch_pyexception([&](){ return impl_->get_probe(id); }, msg); + std::vector<arb::probe_info> get_probes(arb::cell_gid_type gid) const override { + return try_catch_pyexception([&](){ return impl_->get_probes(gid); }, msg); } // TODO: make thread safe diff --git a/python/sampling.cpp b/python/sampling.cpp index d3a18b56055234648f7163f75b052422b3fde3e1..50d68a42c1585ba20d437a3f6dbdd49831979bd3 100644 --- a/python/sampling.cpp +++ b/python/sampling.cpp @@ -60,8 +60,8 @@ struct sample_callback { sample_store(state) {} - 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); + void operator() (arb::probe_metadata pm, std::size_t n, const arb::sample_record* recs) { + auto& v = sample_store->probe_buffer(pm.id); for (std::size_t i = 0; i<n; ++i) { if (auto p = arb::util::any_cast<const double*>(recs[i].data)) { v.push_back({recs[i].time, *p}); diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp index bde7da44e98d4243dc4206147facde7cf775e2e2..b71987cf3e256d6d944fe1ec342bcccc1effff1b 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, arb::util::any_ptr meta, std::size_t n, const arb::sample_record* recs) { + void operator()(arb::probe_metadata, 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)) { @@ -49,7 +49,7 @@ struct trace_callback { trace_.v.push_back(*p); } else { - throw std::runtime_error("unexpected sample type in simple_sampler"); + throw std::runtime_error("unexpected sample type"); } } } @@ -105,19 +105,13 @@ struct single_cell_recipe: arb::recipe { // probes - virtual arb::cell_size_type num_probes(arb::cell_gid_type) const override { - return probes_.size(); - } - - virtual arb::probe_info get_probe(arb::cell_member_type probe_id) const override { - // Test that a valid probe site is requested. - if (probe_id.gid || probe_id.index>=probes_.size()) { - throw arb::bad_probe_id(probe_id); - } - + virtual std::vector<arb::probe_info> get_probes(arb::cell_gid_type gid) const override { // 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::cable_probe_membrane_voltage{loc}}; + std::vector<arb::probe_info> pinfo; + for (auto& p: probes_) { + pinfo.push_back(arb::cable_probe_membrane_voltage{p.site}); + } + return pinfo; } // gap junctions diff --git a/test/simple_recipes.hpp b/test/simple_recipes.hpp index 331148055ce581197a6e83bd6bfa575fe50f8aa6..24051eb005951838cac5e424f5d377563905ac47 100644 --- a/test/simple_recipes.hpp +++ b/test/simple_recipes.hpp @@ -24,19 +24,13 @@ public: cell_gprop_.default_parameters = neuron_parameter_defaults; } - cell_size_type num_probes(cell_gid_type i) const override { - return probes_.count(i)? probes_.at(i).size(): 0; - } - - virtual probe_info get_probe(cell_member_type probe_id) const override { - return probes_.at(probe_id.gid).at(probe_id.index); + std::vector<probe_info> get_probes(cell_gid_type i) const override { + if (!probes_.count(i)) return {}; + return probes_.at(i); } virtual void add_probe(cell_gid_type gid, probe_tag tag, util::any address) { - auto& pvec_ = probes_[gid]; - - cell_member_type probe_id{gid, cell_lid_type(pvec_.size())}; - pvec_.push_back({probe_id, tag, std::move(address)}); + probes_[gid].emplace_back(std::move(address), tag); } util::any get_global_properties(cell_kind k) const override { diff --git a/test/unit-distributed/test_communicator.cpp b/test/unit-distributed/test_communicator.cpp index bfa3eef27f524165f5568fdbf80bd10445230021..16c3527de7f5c2be2135c454cdf20f4cdbb75281 100644 --- a/test/unit-distributed/test_communicator.cpp +++ b/test/unit-distributed/test_communicator.cpp @@ -203,7 +203,6 @@ namespace { cell_size_type num_sources(cell_gid_type) const override { return 1; } cell_size_type num_targets(cell_gid_type) const override { return 1; } - cell_size_type num_probes(cell_gid_type) const override { return 0; } std::vector<cell_connection> connections_on(cell_gid_type gid) const override { // a single connection from the preceding cell, i.e. a ring @@ -266,7 +265,6 @@ namespace { cell_size_type num_sources(cell_gid_type) const override { return 1; } cell_size_type num_targets(cell_gid_type) const override { return size_; } - cell_size_type num_probes(cell_gid_type) const override { return 0; } std::vector<cell_connection> connections_on(cell_gid_type gid) const override { std::vector<cell_connection> cons; diff --git a/test/unit-distributed/test_domain_decomposition.cpp b/test/unit-distributed/test_domain_decomposition.cpp index e4a903d630d87f8580cd3128bdfbc099373dee1e..3b5c336cc7833c5d8822977b21ee00e64624bf26 100644 --- a/test/unit-distributed/test_domain_decomposition.cpp +++ b/test/unit-distributed/test_domain_decomposition.cpp @@ -53,7 +53,6 @@ namespace { cell_size_type num_sources(cell_gid_type) const override { return 0; } cell_size_type num_targets(cell_gid_type) const override { return 0; } - cell_size_type num_probes(cell_gid_type) const override { return 0; } std::vector<cell_connection> connections_on(cell_gid_type) const override { return {}; diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp index c10c138fa8919a00b8e18a1a2d4ff10237322a84..62294e578d5cfbeb3e4e7a90acb91b609226d56e 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<fvm_probe_info> probe_map; + probe_association_map 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<fvm_probe_info> probe_map; + probe_association_map 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<fvm_probe_info> probe_map; + probe_association_map 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}}); - cable_probe_total_ion_current_density where{{1, 0.3}}; + cable_probe_total_ion_current_density where{mlocation{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<fvm_probe_info> probe_map; + probe_association_map probe_map; arb::execution_context context(resources); fvm_cell fvcell(context); @@ -491,10 +491,10 @@ TEST(fvm_lowered, derived_mechs) { std::vector<double> samples[3]; sampler_function sampler = - [&](cell_member_type pid, probe_tag, util::any_ptr, std::size_t n, const sample_record* records) { + [&](probe_metadata pm, 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); + samples[pm.id.gid].push_back(v); } }; @@ -531,7 +531,7 @@ TEST(fvm_lowered, read_valence) { std::vector<target_handle> targets; std::vector<fvm_index_type> cell_to_intdom; - probe_association_map<fvm_probe_info> probe_map; + probe_association_map probe_map; { std::vector<cable_cell> cells(1); @@ -685,7 +685,7 @@ TEST(fvm_lowered, ionic_currents) { std::vector<target_handle> targets; std::vector<fvm_index_type> cell_to_intdom; - probe_association_map<fvm_probe_info> probe_map; + probe_association_map probe_map; fvm_cell fvcell(context); fvcell.initialize({0}, rec, cell_to_intdom, targets, probe_map); @@ -730,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<fvm_probe_info> probe_map; + probe_association_map probe_map; fvm_cell fvcell(context); fvcell.initialize({0}, rec, cell_to_intdom, targets, probe_map); @@ -810,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<fvm_probe_info> probe_map; + probe_association_map probe_map; fvm_cell fvcell(context); fvcell.initialize({0}, rec, cell_to_intdom, targets, probe_map); diff --git a/test/unit/test_lif_cell_group.cpp b/test/unit/test_lif_cell_group.cpp index 2247ef483158a4bff3954d5c30809f9723ddecdf..1041705f4542bd4f0faf813ef853dd3665ea35ef 100644 --- a/test/unit/test_lif_cell_group.cpp +++ b/test/unit/test_lif_cell_group.cpp @@ -72,15 +72,6 @@ public: cell_size_type num_targets(cell_gid_type) const override { return 1; } - cell_size_type num_probes(cell_gid_type) const override { - return 0; - } - probe_info get_probe(cell_member_type probe_id) const override { - return {}; - } - std::vector<event_generator> event_generators(cell_gid_type) const override { - return {}; - } private: cell_size_type n_lif_cells_; @@ -125,15 +116,6 @@ public: cell_size_type num_targets(cell_gid_type) const override { return 1; } - cell_size_type num_probes(cell_gid_type) const override { - return 0; - } - probe_info get_probe(cell_member_type probe_id) const override { - return {}; - } - std::vector<event_generator> event_generators(cell_gid_type) const override { - return {}; - } private: cell_size_type ncells_; diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp index 1450a4e13e6bc0f24c792bc44ae20f15fe8ba6eb..ee17642978dd83596772b83ff5a8016884742f14 100644 --- a/test/unit/test_probe.cpp +++ b/test/unit/test_probe.cpp @@ -112,35 +112,35 @@ void run_v_i_probe_test(const context& ctx) { std::vector<target_handle> targets; std::vector<fvm_index_type> cell_to_intdom; - probe_association_map<fvm_probe_info> probe_map; + probe_association_map probe_map; fvm_cell lcell(*ctx); lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map); - EXPECT_EQ(3u, rec.num_probes(0)); + EXPECT_EQ(3u, rec.get_probes(0).size()); EXPECT_EQ(3u, probe_map.size()); - EXPECT_EQ(10, probe_map.at({0, 0}).tag); - EXPECT_EQ(20, probe_map.at({0, 1}).tag); - EXPECT_EQ(30, probe_map.at({0, 2}).tag); + EXPECT_EQ(10, probe_map.tag.at({0, 0})); + EXPECT_EQ(20, probe_map.tag.at({0, 1})); + EXPECT_EQ(30, probe_map.tag.at({0, 2})); - auto get_probe_handle = [&](cell_member_type x, unsigned i = 0) { - return probe_map.at(x).handle.raw_handle_range()[i]; + auto get_probe_raw_handle = [&](cell_member_type x, unsigned i = 0) { + return probe_map.data_on(x).front().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)); + ASSERT_TRUE(util::get_if<fvm_probe_interpolated>(probe_map.data_on({0, 0}).front().info)); + ASSERT_TRUE(util::get_if<fvm_probe_interpolated>(probe_map.data_on({0, 0}).front().info)); + ASSERT_TRUE(util::get_if<fvm_probe_scalar>(probe_map.data_on({0, 2}).front().info)); - 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}); + probe_handle p0a = get_probe_raw_handle({0, 0}, 0); + probe_handle p0b = get_probe_raw_handle({0, 0}, 1); + probe_handle p1a = get_probe_raw_handle({0, 1}, 0); + probe_handle p1b = get_probe_raw_handle({0, 1}, 1); + probe_handle p2 = get_probe_raw_handle({0, 2}); // Ball-and-stick cell with default discretization policy should // have three CVs, one for branch 0, one trivial one covering the @@ -208,14 +208,14 @@ void run_v_cell_probe_test(const context& ctx) { std::vector<target_handle> targets; std::vector<fvm_index_type> cell_to_intdom; - probe_association_map<fvm_probe_info> probe_map; + probe_association_map 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); + const fvm_probe_multi* h_ptr = util::get_if<fvm_probe_multi>(probe_map.data_on({0, 0}).front().info); ASSERT_TRUE(h_ptr); auto& h = *h_ptr; @@ -268,25 +268,25 @@ void run_expsyn_g_probe_test(const context& ctx) { std::vector<target_handle> targets; std::vector<fvm_index_type> cell_to_intdom; - probe_association_map<fvm_probe_info> probe_map; + probe_association_map 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, rec.get_probes(0).size()); EXPECT_EQ(2u, probe_map.size()); - ASSERT_EQ(1u, probe_map.count({0, 0})); - ASSERT_EQ(1u, probe_map.count({0, 1})); + ASSERT_EQ(1u, probe_map.data.count({0, 0})); + ASSERT_EQ(1u, probe_map.data.count({0, 1})); - EXPECT_EQ(10, probe_map.at({0, 0}).tag); - EXPECT_EQ(20, probe_map.at({0, 1}).tag); + EXPECT_EQ(10, probe_map.tag.at({0, 0})); + EXPECT_EQ(20, probe_map.tag.at({0, 1})); - auto probe_scalar_handle = [&](cell_member_type x) { - return probe_map.at(x).handle.raw_handle_range()[0]; + auto get_probe_raw_handle = [&](cell_member_type x, unsigned i = 0) { + return probe_map.data_on(x).front().raw_handle_range()[i]; }; - probe_handle p0 = probe_scalar_handle({0, 0}); - probe_handle p1 = probe_scalar_handle({0, 1}); + probe_handle p0 = get_probe_raw_handle({0, 0}); + probe_handle p1 = get_probe_raw_handle({0, 1}); // Expect initial probe values to be intial synapse g == 0. @@ -380,7 +380,7 @@ void run_expsyn_g_cell_probe_test(const context& ctx) { std::vector<target_handle> targets; std::vector<fvm_index_type> cell_to_intdom; - probe_association_map<fvm_probe_info> probe_map; + probe_association_map probe_map; fvm_cell lcell(*ctx); lcell.initialize({0, 1}, rec, cell_to_intdom, targets, probe_map); @@ -407,7 +407,7 @@ void run_expsyn_g_cell_probe_test(const context& ctx) { 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); + const auto* h_ptr = util::get_if<fvm_probe_multi>(probe_map.data_on({i, 0}).front().info); ASSERT_TRUE(h_ptr); const auto* m_ptr = util::get_if<std::vector<cable_probe_point_info>>(h_ptr->metadata); @@ -537,7 +537,7 @@ void run_ion_density_probe_test(const context& ctx) { std::vector<target_handle> targets; std::vector<fvm_index_type> cell_to_intdom; - probe_association_map<fvm_probe_info> probe_map; + probe_association_map probe_map; fvm_cell lcell(*ctx); lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map); @@ -547,24 +547,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(13u, rec.num_probes(0)); + EXPECT_EQ(13u, rec.get_probes(0).size()); 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]; + auto get_probe_raw_handle = [&](cell_member_type x, unsigned i = 0) { + return probe_map.data_on(x).front().raw_handle_range()[i]; }; - 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_scalar_handle({0, 7}); - EXPECT_EQ(0u, probe_map.count({0, 8})); - probe_handle write_ca2_s_cv1 = probe_scalar_handle({0, 9}); - probe_handle write_ca2_s_cv2 = probe_scalar_handle({0, 10}); + probe_handle ca_int_cv0 = get_probe_raw_handle({0, 0}); + probe_handle ca_int_cv1 = get_probe_raw_handle({0, 1}); + probe_handle ca_int_cv2 = get_probe_raw_handle({0, 2}); + probe_handle ca_ext_cv0 = get_probe_raw_handle({0, 3}); + probe_handle ca_ext_cv1 = get_probe_raw_handle({0, 4}); + probe_handle ca_ext_cv2 = get_probe_raw_handle({0, 5}); + EXPECT_EQ(0u, probe_map.data.count({0, 6})); + probe_handle na_int_cv2 = get_probe_raw_handle({0, 7}); + EXPECT_EQ(0u, probe_map.data.count({0, 8})); + probe_handle write_ca2_s_cv1 = get_probe_raw_handle({0, 9}); + probe_handle write_ca2_s_cv2 = get_probe_raw_handle({0, 10}); // Ion concentrations should have been written in initialization. // For CV 1, calcium concentration should be mean of the two values @@ -596,9 +596,9 @@ void run_ion_density_probe_test(const context& ctx) { // 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); + auto* p_ptr = util::get_if<fvm_probe_multi>(probe_map.data_on({0, 11}).front().info); ASSERT_TRUE(p_ptr); - fvm_probe_multi& na_int_all_info = *p_ptr; + const 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); @@ -614,9 +614,9 @@ void run_ion_density_probe_test(const context& ctx) { 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); + p_ptr = util::get_if<fvm_probe_multi>(probe_map.data_on({0, 12}).front().info); ASSERT_TRUE(p_ptr); - fvm_probe_multi& ca_ext_all_info = *p_ptr; + const 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); @@ -713,7 +713,7 @@ void run_partial_density_probe_test(const context& ctx) { std::vector<target_handle> targets; std::vector<fvm_index_type> cell_to_intdom; - probe_association_map<fvm_probe_info> probe_map; + probe_association_map probe_map; fvm_cell lcell(*ctx); lcell.initialize({0, 1}, rec, cell_to_intdom, targets, probe_map); @@ -721,12 +721,12 @@ void run_partial_density_probe_test(const context& ctx) { // 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(10u, rec.num_probes(0)); - EXPECT_EQ(10u, rec.num_probes(1)); + EXPECT_EQ(10u, rec.get_probes(0).size()); + EXPECT_EQ(10u, rec.get_probes(1).size()); EXPECT_EQ(10u, 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]; + auto get_probe_raw_handle = [&](cell_member_type x, unsigned i = 0) { + return probe_map.data_on(x).front().raw_handle_range()[i]; }; // Check probe values against expected values. @@ -736,10 +736,10 @@ void run_partial_density_probe_test(const context& ctx) { 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)); + EXPECT_EQ(0u, probe_map.data.count(probe_id)); } else { - probe_handle h = probe_scalar_handle(probe_id); + probe_handle h = get_probe_raw_handle(probe_id); EXPECT_DOUBLE_EQ(tp.expected[gid], deref(h)); } } @@ -810,14 +810,13 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) { 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) + [&](probe_metadata pm, 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); + if (pm.tag==1) { // (whole cell probe) + const mcable_list* m = util::any_cast<const mcable_list*>(pm.meta); ASSERT_NE(nullptr, m); // Metadata should comprise one cable per CV. ASSERT_EQ(n_cv, m->size()); @@ -832,15 +831,15 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) { } else { // axial current probe // Probe id tells us which axial current this is. - ASSERT_LT(probe_id.index, n_axial_probe); + ASSERT_LT(pm.id.index, n_axial_probe); - const mlocation* m = util::any_cast<const mlocation*>(metadata); + const mlocation* m = util::any_cast<const mlocation*>(pm.meta); 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; + i_axial.at(pm.id.index) = *s; } }); @@ -866,7 +865,8 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) { // Run given cells taking samples from the provied probes on one of the cells. -// Use default mechanism catalogue augmented by unit test specific mechanisms. +// +// Use the default mechanism catalogue augmented by unit test specific mechanisms. // (Timestep fixed at 0.025 ms). template <typename SampleData, typename SampleMeta = void> @@ -891,7 +891,7 @@ auto run_simple_samplers( }; simulation sim(rec, partition_load_balance(rec, ctx, phints), ctx); - std::vector<trace_data<SampleData, SampleMeta>> traces(n_probe); + std::vector<trace_vector<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])); } @@ -912,6 +912,39 @@ auto run_simple_sampler( return run_simple_samplers<SampleData, SampleMeta>(ctx, t_end, cells, probe_cell, {probe_addr}, when).at(0); } +template <typename Backend> +void run_multi_probe_test(const context& ctx) { + // Construct and run thorugh simple sampler a probe defined over + // cell terminal points; check metadata and values. + + // m_mlt_b6 has terminal branches 1, 2, 4, and 5. + cable_cell cell(common_morphology::m_mlt_b6); + + // Paint mechanism on branches 1, 2, and 5, omitting branch 4. + cell.paint(reg::branch(1), mechanism_desc("param_as_state").set("p", 10.)); + cell.paint(reg::branch(2), mechanism_desc("param_as_state").set("p", 20.)); + cell.paint(reg::branch(5), mechanism_desc("param_as_state").set("p", 50.)); + + auto tracev = run_simple_sampler<double, mlocation>(ctx, 0.1, {cell}, 0, cable_probe_density_state{ls::terminal(), "param_as_state", "s"}, {0.}); + + // Expect to have received a sample on each of the terminals of branches 1, 2, and 5. + ASSERT_EQ(3u, tracev.size()); + + std::vector<std::pair<mlocation, double>> vals; + for (auto& trace: tracev) { + ASSERT_EQ(1u, trace.size()); + vals.push_back({trace.meta, trace[0].v}); + } + + util::sort(vals); + EXPECT_EQ((mlocation{1, 1.}), vals[0].first); + EXPECT_EQ((mlocation{2, 1.}), vals[1].first); + EXPECT_EQ((mlocation{5, 1.}), vals[2].first); + EXPECT_EQ(10., vals[0].second); + EXPECT_EQ(20., vals[1].second); + EXPECT_EQ(50., vals[2].second); +} + template <typename Backend> void run_v_sampled_probe_test(const context& ctx) { cable_cell bs = make_cell_ball_and_stick(false); @@ -930,13 +963,17 @@ void run_v_sampled_probe_test(const context& ctx) { 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); + auto trace0 = run_simple_sampler<double, mlocation>(ctx, t_end, cells, 0, cable_probe_membrane_voltage{probe_loc}, when).at(0); + ASSERT_TRUE(trace0); + + EXPECT_EQ(probe_loc, trace0.meta); 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); + auto trace1 = run_simple_sampler<double, mlocation>(ctx, t_end, cells, 1, cable_probe_membrane_voltage{probe_loc}, when).at(0); + ASSERT_TRUE(trace1); + + EXPECT_EQ(probe_loc, trace1.meta); 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); @@ -997,10 +1034,10 @@ void run_total_current_probe_test(const context& ctx) { 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}); + cable_probe_total_current_cell{}, {tau, 20*tau}).at(0); 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}); + cable_probe_total_ion_current_cell{}, {tau, 20*tau}).at(0); ASSERT_EQ(2u, traces[i].size()); ASSERT_EQ(2u, ion_traces[i].size()); @@ -1014,8 +1051,8 @@ void run_total_current_probe_test(const context& ctx) { // 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); + ASSERT_EQ((n_cv_per_branch+(int)interior_forks)*n_branch, traces[i].meta.size()); + EXPECT_EQ(ion_traces[i].meta, traces[i].meta); 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()); @@ -1100,10 +1137,8 @@ void run_exact_sampling_probe_test(const context& ctx) { return cell_kind::cable; } - cell_size_type num_probes(cell_gid_type) const override { return 1; } - - probe_info get_probe(cell_member_type pid) const override { - return {pid, 0, cable_probe_membrane_voltage{mlocation{1, 0.5}}}; + std::vector<probe_info> get_probes(cell_gid_type) const override { + return {cable_probe_membrane_voltage{mlocation{1, 0.5}}}; } cell_size_type num_targets(cell_gid_type) const override { return 1; } @@ -1138,7 +1173,7 @@ void run_exact_sampling_probe_test(const context& ctx) { // 1. Membrane voltage is similar with and without exact sampling. // 2. Sample times are in fact exact with exact sampling. - std::vector<trace_data<double>> lax_traces(4), exact_traces(4); + std::vector<trace_vector<double>> lax_traces(4), exact_traces(4); const double max_dt = 0.001; const double t_end = 1.; @@ -1167,16 +1202,21 @@ void run_exact_sampling_probe_test(const context& ctx) { exact_sim.run(t_end, max_dt); for (unsigned i = 0; i<n_cell; ++i) { - ASSERT_EQ(n_sample_time, lax_traces.at(i).size()); - ASSERT_EQ(n_sample_time, exact_traces.at(i).size()); + ASSERT_EQ(1u, lax_traces.at(i).size()); + ASSERT_EQ(n_sample_time, lax_traces.at(i).at(0).size()); + ASSERT_EQ(1u, exact_traces.at(i).size()); + ASSERT_EQ(n_sample_time, exact_traces.at(i).at(0).size()); } for (unsigned i = 0; i<n_cell; ++i) { + auto& lax_trace = lax_traces[i][0]; + auto& exact_trace = exact_traces[i][0]; + for (unsigned j = 0; j<n_sample_time; ++j) { - EXPECT_NE(sched_times.at(j), lax_traces.at(i).at(j).t); - EXPECT_EQ(sched_times.at(j), exact_traces.at(i).at(j).t); + EXPECT_NE(sched_times.at(j), lax_trace.at(j).t); + EXPECT_EQ(sched_times.at(j), exact_trace.at(j).t); - EXPECT_TRUE(testing::near_relative(lax_traces.at(i).at(j).v, exact_traces.at(i).at(j).v, 0.01)); + EXPECT_TRUE(testing::near_relative(lax_trace.at(j).v, exact_trace.at(j).v, 0.01)); } } } @@ -1187,7 +1227,8 @@ void run_exact_sampling_probe_test(const context& ctx) { #undef PROBE_TESTS #define PROBE_TESTS \ v_i, v_cell, v_sampled, expsyn_g, expsyn_g_cell, \ - ion_density, axial_and_ion_current_sampled, partial_density, total_current, exact_sampling + ion_density, axial_and_ion_current_sampled, partial_density, total_current, exact_sampling, \ + multi #undef RUN_MULTICORE #define RUN_MULTICORE(x) \