diff --git a/.github/workflows/benchmarks.yml b/.github/workflows/benchmarks.yml index be1014b850..0326c18367 100644 --- a/.github/workflows/benchmarks.yml +++ b/.github/workflows/benchmarks.yml @@ -8,12 +8,12 @@ on: jobs: build: name: "Benchmarks" - runs-on: ubuntu-22.04 + runs-on: ubuntu-24.04 strategy: fail-fast: false env: - CC: gcc-12 - CXX: g++-12 + CC: gcc-14 + CXX: g++-14 steps: - name: Get build dependencies run: | diff --git a/.gitignore b/.gitignore index dcca1025f7..e1c109c982 100644 --- a/.gitignore +++ b/.gitignore @@ -106,3 +106,4 @@ _skbuild # generated by test scripts results /_deps/ +/install-* diff --git a/arbor/cable_cell_group.cpp b/arbor/cable_cell_group.cpp index f04ac941d6..e9c41adbfc 100644 --- a/arbor/cable_cell_group.cpp +++ b/arbor/cable_cell_group.cpp @@ -26,9 +26,7 @@ cable_cell_group::cable_cell_group(const std::vector& gids, cell_label_range& cg_sources, cell_label_range& cg_targets, fvm_lowered_cell_ptr lowered): - gids_(gids), lowered_(std::move(lowered)) -{ - + gids_(gids), lowered_(std::move(lowered)) { // Construct cell implementation, retrieving handles and maps. auto fvm_info = lowered_->initialize(gids_, rec); @@ -49,11 +47,9 @@ cable_cell_group::cable_cell_group(const std::vector& gids, void cable_cell_group::reset() { spikes_.clear(); - for (auto &entry: sampler_map_) { entry.second.sched.reset(); } - lowered_->reset(); } @@ -76,211 +72,123 @@ void cable_cell_group::t_deserialize(serializer& ser, const std::string& k) { deserialize(ser, k, *this); } // Working space for computing and collating data for samplers. -using fvm_probe_scratch = std::tuple, std::vector>; - -template -void tuple_foreach(VoidFn&& f, std::tuple& t) { - // executes functions in order (pack expansion) - // uses comma operator (unary left fold) - // avoids potentially overloaded comma operator (cast to void) - std::apply( - [g = std::forward(f)](auto&& ...x){ - (..., static_cast(g(std::forward(x))));}, - t); -} +struct fvm_probe_scratch { + std::vector times; + std::vector values; -void reserve_scratch(fvm_probe_scratch& scratch, std::size_t n) { - tuple_foreach([n](auto& v) { v.reserve(n); }, scratch); -} + void reserve(size_t n) { + values.reserve(n); + times.reserve(n); + } +}; -void run_samples( - const missing_probe_info&, - const sampler_call_info&, - const arb_value_type*, - const arb_value_type*, - std::vector&, - fvm_probe_scratch&) -{ +void run_samples(const missing_probe_info&, + const sampler_call_info&, + const arb_value_type*, + const arb_value_type*, + const fvm_probe_scratch&) { throw arbor_internal_error("invalid fvm_probe_data in sampler map"); } -void run_samples( - const fvm_probe_scalar& p, - const sampler_call_info& sc, - const arb_value_type* raw_times, - const arb_value_type* raw_samples, - std::vector& sample_records, - fvm_probe_scratch&) -{ - // Scalar probes do not need scratch space — provided that the user-presented - // sample type (double) matches the raw type (arb_value_type). - static_assert(std::is_same::value, "require sample value translation"); - - sample_size_type n_sample = sc.end_offset-sc.begin_offset; - sample_records.clear(); - for (auto i = sc.begin_offset; i!=sc.end_offset; ++i) { - sample_records.push_back(sample_record{time_type(raw_times[i]), &raw_samples[i]}); - } - - sc.sampler({sc.probeset_id, sc.index, p.get_metadata_ptr()}, n_sample, sample_records.data()); +template +void do_run_sampler(const sampler_call_info& sc, + std::size_t n_sample, + std::size_t width, + const P& probe, + const fvm_probe_scratch& scratch) { + arb_assert(scratch.values.size() == n_sample * width); + arb_assert(scratch.times.size() == n_sample); + sc.sampler(probe_metadata { .id=sc.probeset_id, + .index=sc.index, + .meta=probe.get_metadata_ptr() }, + sample_records { .n_sample=n_sample, + .width=width, + .time=scratch.times.data(), + .values=const_cast(scratch.values.data()) }); } - -void run_samples( - const fvm_probe_interpolated& p, - const sampler_call_info& sc, - const arb_value_type* raw_times, - const arb_value_type* raw_samples, - std::vector& sample_records, - fvm_probe_scratch& scratch) -{ - constexpr sample_size_type n_raw_per_sample = 2; - sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample; - arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample); - - auto& tmp = std::get>(scratch); - tmp.clear(); - sample_records.clear(); - - for (sample_size_type j = 0; j& sample_records, - fvm_probe_scratch& scratch) -{ + +void run_samples(const fvm_probe_multi& p, + const sampler_call_info& sc, + const arb_value_type* raw_times, + const arb_value_type* raw_samples, + fvm_probe_scratch& scratch) { const sample_size_type n_raw_per_sample = p.raw_handles.size(); - sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample; - arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample); + sample_size_type n_sample = (sc.end_offset - sc.begin_offset)/n_raw_per_sample; + arb_assert((sc.end_offset - sc.begin_offset) == n_sample*n_raw_per_sample); - auto& sample_ranges = std::get>(scratch); - sample_ranges.clear(); - sample_records.clear(); + scratch.times.clear(); + scratch.values.clear(); + scratch.values.reserve(n_raw_per_sample*n_sample); - for (sample_size_type j = 0; j& sample_records, - fvm_probe_scratch& scratch) -{ +void run_samples(const fvm_probe_weighted_multi& p, + const sampler_call_info& sc, + const arb_value_type* raw_times, + const arb_value_type* raw_samples, + fvm_probe_scratch& scratch) { const sample_size_type n_raw_per_sample = p.raw_handles.size(); sample_size_type n_sample = (sc.end_offset - sc.begin_offset)/n_raw_per_sample; arb_assert((sc.end_offset - sc.begin_offset)==n_sample*n_raw_per_sample); arb_assert((unsigned)n_raw_per_sample == p.weight.size()); - auto& sample_ranges = std::get>(scratch); - sample_ranges.clear(); - sample_records.clear(); - - auto& tmp = std::get>(scratch); - tmp.clear(); - tmp.reserve(n_raw_per_sample*n_sample); + scratch.times.clear(); + scratch.values.clear(); + scratch.values.reserve(n_raw_per_sample*n_sample); - for (sample_size_type j = 0; j& sample_records, fvm_probe_scratch& scratch) { const sample_size_type n_raw_per_sample = p.raw_handles.size(); const sample_size_type n_interp_per_sample = n_raw_per_sample/2; - sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample; - arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample); - arb_assert((unsigned)n_interp_per_sample==p.coef[0].size()); - arb_assert((unsigned)n_interp_per_sample==p.coef[1].size()); - - auto& sample_ranges = std::get>(scratch); - sample_ranges.clear(); - sample_records.clear(); + sample_size_type n_sample = (sc.end_offset - sc.begin_offset)/n_raw_per_sample; + arb_assert((sc.end_offset - sc.begin_offset) == n_sample*n_raw_per_sample); + arb_assert((unsigned)n_interp_per_sample == p.coef[0].size()); + arb_assert((unsigned)n_interp_per_sample == p.coef[1].size()); - auto& tmp = std::get>(scratch); - tmp.clear(); - tmp.reserve(n_interp_per_sample*n_sample); + scratch.times.clear(); + scratch.values.clear(); + scratch.values.reserve(n_interp_per_sample*n_sample); - for (sample_size_type j = 0; j& sample_records, - fvm_probe_scratch& scratch) -{ +void run_samples(const fvm_probe_membrane_currents& p, + const sampler_call_info& sc, + const arb_value_type* raw_times, + const arb_value_type* raw_samples, + fvm_probe_scratch& scratch) { const sample_size_type n_raw_per_sample = p.raw_handles.size(); sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample; arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample); @@ -291,21 +199,17 @@ void run_samples( const auto n_stim = p.stim_scale.size(); arb_assert(n_stim+n_cv==(unsigned)n_raw_per_sample); - auto& sample_ranges = std::get>(scratch); - sample_ranges.clear(); - - auto& tmp = std::get>(scratch); - tmp.assign(n_cable*n_sample, 0.); + scratch.times.clear(); + scratch.values.clear(); + scratch.values.assign(n_cable*n_sample, 0.); - sample_records.clear(); - - for (sample_size_type j = 0; j& sample_records, - fvm_probe_scratch& scratch) -{ - std::visit([&](auto& x) {run_samples(x, sc, raw_times, raw_samples, sample_records, scratch); }, sc.pdata_ptr->info); +void run_samples(const sampler_call_info& sc, + const arb_value_type* raw_times, + const arb_value_type* raw_samples, + fvm_probe_scratch& scratch) { + std::visit([&](auto& x) { run_samples(x, sc, raw_times, raw_samples, scratch); }, + sc.pdata_ptr->info); } void cable_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& event_lanes) { @@ -432,14 +326,12 @@ void cable_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange // and then call the callback. PE(advance:sampledeliver); - std::vector sample_records; - sample_records.reserve(max_samples_per_call); fvm_probe_scratch scratch; - reserve_scratch(scratch, max_samples_per_call); + scratch.reserve(max_samples_per_call); for (auto& sc: call_info) { - run_samples(sc, result.sample_time.data(), result.sample_value.data(), sample_records, scratch); + run_samples(sc, result.sample_time.data(), result.sample_value.data(), scratch); } PL(); @@ -447,7 +339,6 @@ void cable_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange // generate spikes with global spike source ids. The threshold crossings // record the local spike source index, which must be converted to a // global index for spike communication. - for (auto c: result.crossings) { spikes_.emplace_back(spike_sources_[c.index], time_type(c.time)); } @@ -481,13 +372,13 @@ void cable_cell_group::remove_all_samplers() { std::vector cable_cell_group::get_probe_metadata(const cell_address_type& probeset_id) const { // SAFETY: Probe associations are fixed after construction, so we do not // need to grab the mutex. - auto data = probe_map_.data_on(probeset_id); + const auto& data = probe_map_.data_on(probeset_id); std::vector result; result.reserve(data.size()); unsigned index = 0; for (const auto& info: data) { - result.push_back({probeset_id, index++, info->get_metadata_ptr()}); + result.push_back({probeset_id, index++, info->get_width(), info->get_metadata_ptr()}); } return result; } diff --git a/arbor/execution_context.hpp b/arbor/execution_context.hpp index b14fe01b34..db72bcc030 100644 --- a/arbor/execution_context.hpp +++ b/arbor/execution_context.hpp @@ -1,7 +1,5 @@ #pragma once -#include - #include #include diff --git a/arbor/fvm_lowered_cell.hpp b/arbor/fvm_lowered_cell.hpp index 896380c019..4612784fee 100644 --- a/arbor/fvm_lowered_cell.hpp +++ b/arbor/fvm_lowered_cell.hpp @@ -34,27 +34,9 @@ namespace arb { // a sample value (possibly non-scalar) presented to a sampler // function, and one or more probe handles that reference data // in the FVM back-end. - -struct fvm_probe_scalar { - probe_handle raw_handles[1] = {nullptr}; - std::variant metadata; - - util::any_ptr get_metadata_ptr() const { - return std::visit([](const auto& x) -> util::any_ptr { return &x; }, metadata); - } -}; - -struct fvm_probe_interpolated { - probe_handle raw_handles[2] = {nullptr, nullptr}; - double coef[2] = {}; - mlocation metadata; - - util::any_ptr get_metadata_ptr() const { return &metadata; } -}; - struct fvm_probe_multi { std::vector raw_handles; - std::variant> metadata; + std::variant> metadata; void shrink_to_fit() { raw_handles.shrink_to_fit(); @@ -62,8 +44,13 @@ struct fvm_probe_multi { } util::any_ptr get_metadata_ptr() const { - return std::visit([](const auto& x) -> util::any_ptr { return &x; }, metadata); + return std::visit([](const auto& x) -> util::any_ptr { return x.data(); }, metadata); } + + std::size_t get_width() const { + return std::visit([](const auto& x) -> std::size_t { return x.size(); }, metadata); + } + }; struct fvm_probe_weighted_multi { @@ -76,29 +63,36 @@ struct fvm_probe_weighted_multi { weight.shrink_to_fit(); metadata.shrink_to_fit(); } - - util::any_ptr get_metadata_ptr() const { return &metadata; } + std::size_t get_width() const { return metadata.size(); } + util::any_ptr get_metadata_ptr() const { return metadata.data(); } }; struct fvm_probe_interpolated_multi { std::vector raw_handles; // First half take coef[0], second half coef[1]. std::vector coef[2]; - mcable_list metadata; + std::variant metadata; void shrink_to_fit() { raw_handles.shrink_to_fit(); coef[0].shrink_to_fit(); coef[1].shrink_to_fit(); - metadata.shrink_to_fit(); + std::visit([](auto& v) { v.shrink_to_fit(); }, metadata); + } + + util::any_ptr get_metadata_ptr() const { + return std::visit([](const auto& x) -> util::any_ptr { return x.data(); }, metadata); + } + + std::size_t get_width() const { + return std::visit([](const auto& x) -> std::size_t { return x.size(); }, metadata); } - util::any_ptr get_metadata_ptr() const { return &metadata; } }; // Trans-membrane currents require special handling! struct fvm_probe_membrane_currents { std::vector raw_handles; // Voltage per CV, followed by stim current densities. - std::vector metadata; // Cables from each CV, in CV order. + mcable_list metadata; // Cables from each CV, in CV order. std::vector cv_parent; // Parent CV index for each CV. std::vector cv_parent_cond; // Face conductance between CV and parent. @@ -119,7 +113,8 @@ struct fvm_probe_membrane_currents { stim_cv.shrink_to_fit(); } - util::any_ptr get_metadata_ptr() const { return &metadata; } + std::size_t get_width() const { return metadata.size(); } + util::any_ptr get_metadata_ptr() const { return metadata.data(); } }; struct missing_probe_info { @@ -127,27 +122,22 @@ struct missing_probe_info { std::array raw_handles; void* metadata = nullptr; - util::any_ptr get_metadata_ptr() const { return util::any_ptr{}; } + std::size_t get_width() const { return 0; } + util::any_ptr get_metadata_ptr() const { return nullptr; } }; 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_interpolated_multi p): info(std::move(p)) {} fvm_probe_data(fvm_probe_membrane_currents p): info(std::move(p)) {} - std::variant< - missing_probe_info, - fvm_probe_scalar, - fvm_probe_interpolated, - fvm_probe_multi, - fvm_probe_weighted_multi, - fvm_probe_interpolated_multi, - fvm_probe_membrane_currents - > info = missing_probe_info{}; + std::variant info = missing_probe_info{}; auto raw_handle_range() const { return util::make_range( @@ -164,6 +154,10 @@ struct fvm_probe_data { return std::visit([](const auto& i) -> util::any_ptr { return i.get_metadata_ptr(); }, info); } + std::size_t get_width() const { + return std::visit([](const auto& i) -> std::size_t { return i.get_width(); }, info); + } + sample_size_type n_raw() const { return raw_handle_range().size(); } explicit operator bool() const { return !std::get_if(&info); } diff --git a/arbor/fvm_lowered_cell_impl.cpp b/arbor/fvm_lowered_cell_impl.cpp index 2a0d15d3b4..971d1dc3f6 100644 --- a/arbor/fvm_lowered_cell_impl.cpp +++ b/arbor/fvm_lowered_cell_impl.cpp @@ -1,6 +1,3 @@ -#include -#include - #include #include #include diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp index adcbab1a97..bad9f2703c 100644 --- a/arbor/fvm_lowered_cell_impl.hpp +++ b/arbor/fvm_lowered_cell_impl.hpp @@ -106,25 +106,24 @@ struct fvm_lowered_cell_impl: public fvm_lowered_cell { } // Translate cell probe descriptions into probe handles etc. - void resolve_probe_address( - std::vector& probe_data, // out parameter - const std::vector& cells, - std::size_t cell_idx, - const std::any& paddr, - const fvm_cv_discretization& D, - const fvm_mechanism_data& M, - const std::vector& handles, - const std::unordered_map& mech_instance_by_name); - - // Add probes to fvm_info::probe_map - void add_probes(const std::vector& gids, - const std::vector& cells, - const recipe& rec, - const fvm_cv_discretization& D, - const std::unordered_map& mechptr_by_name, - const fvm_mechanism_data& mech_data, - const std::vector& target_handles, - probe_association_map& probe_map); + void resolve_probe_address(std::vector& probe_data, // out parameter + const std::vector& cells, + std::size_t cell_idx, + const std::any& paddr, + const fvm_cv_discretization& D, + const fvm_mechanism_data& M, + const std::vector& handles, + const std::unordered_map& mech_instance_by_name); + + // Add probes to fvm_info::probe_map + void add_probes(const std::vector& gids, + const std::vector& cells, + const recipe& rec, + const fvm_cv_discretization& D, + const std::unordered_map& mechptr_by_name, + const fvm_mechanism_data& mech_data, + const std::vector& target_handles, + probe_association_map& probe_map); }; template @@ -336,7 +335,6 @@ fvm_lowered_cell_impl::add_probes(const std::vector& gid const std::vector& target_handles, probe_association_map& probe_map) { auto ncell = gids.size(); - std::vector probe_data; for (auto cell_idx: util::make_span(ncell)) { cell_gid_type gid = gids[cell_idx]; @@ -346,10 +344,9 @@ fvm_lowered_cell_impl::add_probes(const std::vector& gid if (!probe_data.empty()) { cell_address_type addr{gid, pi.tag}; if (probe_map.count(addr)) throw dup_cell_probe(cell_kind::cable, gid, pi.tag); - for (auto& data: probe_data) { - probe_map.insert(addr, std::move(data)); - } + for (auto& data: probe_data) probe_map.insert(addr, std::move(data)); } + probe_data.clear(); } } } @@ -646,57 +643,75 @@ void fvm_lowered_cell_impl::resolve_probe_address(std::vector& handles, const std::unordered_map& mech_instance_by_name) { - probe_data.clear(); probe_resolution_data prd{ - 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, - cable_probe_membrane_voltage_cell, - cable_probe_axial_current, - cable_probe_total_ion_current_density, - cable_probe_total_ion_current_cell, - cable_probe_total_current_cell, - cable_probe_stimulus_current_cell, - cable_probe_density_state, - cable_probe_density_state_cell, - cable_probe_point_state, - cable_probe_point_state_cell, - cable_probe_ion_current_density, - cable_probe_ion_current_cell, - cable_probe_ion_int_concentration, - cable_probe_ion_int_concentration_cell, - cable_probe_ion_diff_concentration, - cable_probe_ion_diff_concentration_cell, - cable_probe_ion_ext_concentration, - cable_probe_ion_ext_concentration_cell>; - - auto visitor = util::overload( - [&prd](auto& probe_addr) { resolve_probe(probe_addr, prd); }, - [] { throw cable_cell_error("unrecognized probe type"), fvm_probe_data{}; }); + .result=probe_data, + .state=state_.get(), + .cell=cells[cell_idx], + .cell_idx=cell_idx, + .D=D, + .M=M, + .handles=handles, + .mech_instance_by_name=mech_instance_by_name + }; + + using V = util::any_visitor; + + auto visitor = util::overload([&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 -void resolve_probe(const cable_probe_membrane_voltage& p, probe_resolution_data& R) { - const arb_value_type* data = R.state->voltage.data(); - - 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}); +void resolve_probe(const cable_probe_membrane_voltage& p, probe_resolution_data& res) { + const arb_value_type* data = res.state->voltage.data(); + + std::vector handles_p, handles_d; + std::vector coef_p; + std::vector coef_d; + mlocation_list meta; + for (const mlocation& loc: thingify(p.locations, res.cell.provider())) { + const auto& in = fvm_interpolate_voltage(res.cell, res.D, res.cell_idx, loc); + handles_p.push_back(data + in.proximal_cv); + handles_d.push_back(data + in.distal_cv); + coef_p.push_back(in.proximal_coef); + coef_d.push_back(in.distal_coef); + meta.push_back(loc); } + util::append(handles_p, handles_d); + handles_p.shrink_to_fit(); + meta.shrink_to_fit(); + coef_p.shrink_to_fit(); + coef_d.shrink_to_fit(); + res.result.push_back(fvm_probe_interpolated_multi{ + .raw_handles=std::move(handles_p), + .coef={std::move(coef_p), std::move(coef_d)}, + .metadata=std::move(meta), + }); } template void resolve_probe(const cable_probe_membrane_voltage_cell& p, probe_resolution_data& R) { fvm_probe_multi r; mcable_list cables; - for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) { const double* ptr = R.state->voltage.data()+cv; for (auto cable: R.D.geometry.cables(cv)) { @@ -706,46 +721,72 @@ void resolve_probe(const cable_probe_membrane_voltage_cell& p, probe_resolution_ } r.metadata = std::move(cables); r.shrink_to_fit(); - R.result.push_back(std::move(r)); } template -void resolve_probe(const cable_probe_axial_current& p, probe_resolution_data& R) { - const arb_value_type* data = R.state->voltage.data(); - - 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}); +void resolve_probe(const cable_probe_axial_current& p, probe_resolution_data& res) { + const arb_value_type* data = res.state->voltage.data(); + + std::vector handles_p, handles_d; + std::vector coef_p; + std::vector coef_d; + mlocation_list meta; + for (const mlocation& loc: thingify(p.locations, res.cell.provider())) { + const auto& in = fvm_axial_current(res.cell, res.D, res.cell_idx, loc); + handles_p.push_back(data + in.proximal_cv); + handles_d.push_back(data + in.distal_cv); + coef_p.push_back(in.proximal_coef); + coef_d.push_back(in.distal_coef); + meta.push_back(loc); } + util::append(handles_p, handles_d); + handles_p.shrink_to_fit(); + meta.shrink_to_fit(); + coef_p.shrink_to_fit(); + coef_d.shrink_to_fit(); + if (meta.empty()) return; + res.result.push_back(fvm_probe_interpolated_multi{ + .raw_handles=std::move(handles_p), + .coef={std::move(coef_p), std::move(coef_d)}, + .metadata=std::move(meta), + }); } template -void resolve_probe(const cable_probe_total_ion_current_density& p, probe_resolution_data& R) { +void resolve_probe(const cable_probe_total_ion_current_density& p, probe_resolution_data& res) { // Use interpolated probe with coeffs 1, -1 to represent difference between accumulated current density and stimulus. - for (mlocation loc: thingify(p.locations, R.cell.provider())) { - arb_index_type cv = R.D.geometry.location_cv(R.cell_idx, loc, cv_prefer::cv_nonempty); - const double* current_cv_ptr = R.state->current_density.data() + cv; - - auto opt_i = util::binary_search_index(R.M.stimuli.cv_unique, cv); - const double* stim_cv_ptr = opt_i? R.state->stim_data.accu_stim_.data()+*opt_i: nullptr; - - R.result.push_back(fvm_probe_interpolated{ - {current_cv_ptr, stim_cv_ptr}, - {1., -1.}, - loc}); + std::vector handles_p, handles_d; + std::vector coef_p; + std::vector coef_d; + mlocation_list meta; + for (const mlocation& loc: thingify(p.locations, res.cell.provider())) { + arb_index_type cv = res.D.geometry.location_cv(res.cell_idx, loc, cv_prefer::cv_nonempty); + auto opt_i = util::binary_search_index(res.M.stimuli.cv_unique, cv); + handles_p.push_back(res.state->current_density.data() + cv); + handles_d.push_back(opt_i ? res.state->stim_data.accu_stim_.data()+ *opt_i: nullptr); + coef_p.push_back( 1.0); + coef_d.push_back(-1.0); + meta.push_back(loc); } + util::append(handles_p, handles_d); + handles_p.shrink_to_fit(); + meta.shrink_to_fit(); + coef_p.shrink_to_fit(); + coef_d.shrink_to_fit(); + if (meta.empty()) return; + res.result.push_back(fvm_probe_interpolated_multi{ + .raw_handles=std::move(handles_p), + .coef={std::move(coef_p), std::move(coef_d)}, + .metadata=std::move(meta), + }); } template void resolve_probe(const cable_probe_total_ion_current_cell& p, probe_resolution_data& R) { fvm_probe_interpolated_multi r; std::vector stim_handles; - + mcable_list meta; for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) { const double* current_cv_ptr = R.state->current_density.data()+cv; auto opt_i = util::binary_search_index(R.M.stimuli.cv_unique, cv); @@ -758,11 +799,12 @@ void resolve_probe(const cable_probe_total_ion_current_cell& p, probe_resolution stim_handles.push_back(stim_cv_ptr); r.coef[0].push_back(0.001*area); // Scale from [µm²·A/m²] to [nA]. r.coef[1].push_back(-r.coef[0].back()); - r.metadata.push_back(cable); + meta.push_back(cable); } } } - + if (meta.empty()) return; + r.metadata = std::move(meta); util::append(r.raw_handles, stim_handles); r.shrink_to_fit(); R.result.push_back(std::move(r)); @@ -805,6 +847,7 @@ void resolve_probe(const cable_probe_total_current_cell& p, probe_resolution_dat r.stim_scale.push_back(0.001*R.D.cv_area[cv]); // Scale from [µm²·A/m²] to [nA]. } r.shrink_to_fit(); + if (r.metadata.empty()) return; R.result.push_back(std::move(r)); } @@ -828,8 +871,8 @@ void resolve_probe(const cable_probe_stimulus_current_cell& p, probe_resolution_ } } } - r.shrink_to_fit(); + if (r.metadata.empty()) return; R.result.push_back(std::move(r)); } @@ -837,19 +880,27 @@ template void resolve_probe(const cable_probe_density_state& p, probe_resolution_data& R) { const auto& mech = p.mechanism; if (!R.mech_instance_by_name.count(mech)) return; - const arb_value_type* data = R.mechanism_state(mech, p.state); + + const auto* data = R.mechanism_state(mech, p.state); if (!data) return; - auto support = R.mechanism_support(mech); - for (mlocation loc: thingify(p.locations, R.cell.provider())) { + const auto& support = R.mechanism_support(mech); + mlocation_list meta; + std::vector handles; + for (const mlocation& loc: thingify(p.locations, R.cell.provider())) { if (!support.intersects(loc)) continue; arb_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(mech).cv, cv); if (!opt_i) continue; - R.result.push_back(fvm_probe_scalar{{data+*opt_i}, loc}); + handles.push_back(data + *opt_i); + meta.push_back(loc); } + handles.shrink_to_fit(); + meta.shrink_to_fit(); + if (meta.empty()) return; + R.result.push_back(fvm_probe_multi{.raw_handles=std::move(handles), .metadata=std::move(meta)}); } template @@ -918,6 +969,8 @@ void resolve_probe(const cable_probe_point_state& p, probe_resolution_data& R const auto& [lr_beg, lr_end] = R.cell .synapse_ranges() .equal_range(t_hash); + std::vector meta; + std::vector handles; for (auto lr = lr_beg; lr != lr_end; ++lr) { const auto& [lid_beg, lid_end] = lr->second; for (auto lid = lid_beg; lid != lid_end; ++lid) { @@ -926,14 +979,22 @@ void resolve_probe(const cable_probe_point_state& p, probe_resolution_data& R const auto& handle = R.handles.at(cg); if (handle.mech_id != mech_id) return; auto mech_index = handle.mech_index; - R.result.push_back(fvm_probe_scalar{{data + mech_index}, - point_info_of(target, - lid, - mech_index, - synapses.at(mech), - R.M.mechanisms.at(mech).multiplicity)}); + meta.push_back(point_info_of(target, + lid, + mech_index, + synapses.at(mech), + R.M.mechanisms.at(mech).multiplicity)); + handles.push_back(data + mech_index); } } + if (meta.empty()) return; + meta.shrink_to_fit(); + handles.shrink_to_fit(); + R.result.push_back( + fvm_probe_multi{ + .raw_handles=handles, + .metadata=meta, + }); } template @@ -987,12 +1048,24 @@ void resolve_probe(const cable_probe_point_state_cell& p, probe_resolution_data< template void resolve_probe(const cable_probe_ion_current_density& p, probe_resolution_data& R) { - for (mlocation loc: thingify(p.locations, R.cell.provider())) { + if(!R.state->ion_data.count(p.ion)) return; + mlocation_list meta; + std::vector handles; + const auto data = R.state->ion_data.at(p.ion).iX_.data(); + for (const auto& loc: thingify(p.locations, R.cell.provider())) { auto opt_i = R.ion_location_index(p.ion, loc); if (!opt_i) continue; - - R.result.push_back(fvm_probe_scalar{{R.state->ion_data.at(p.ion).iX_.data()+*opt_i}, loc}); + handles.push_back(data + *opt_i); + meta.push_back(loc); } + if (meta.empty()) return; + meta.shrink_to_fit(); + handles.shrink_to_fit(); + R.result.push_back( + fvm_probe_multi{ + .raw_handles=handles, + .metadata=meta, + }); } template @@ -1021,42 +1094,55 @@ void resolve_probe(const cable_probe_ion_current_cell& p, probe_resolution_data< } template -void resolve_probe(const cable_probe_ion_int_concentration& p, probe_resolution_data& R) { - const auto& ion = p.ion; - const auto& xi = R.state->ion_data.at(ion).Xi_; - if (xi.empty()) return; - for (mlocation loc: thingify(p.locations, R.cell.provider())) { +void resolve_ion_conc_common(const locset& ls, + const std::string& ion, + const typename B::array& values, + probe_resolution_data& R) { + if (values.empty()) return; + mlocation_list meta; + std::vector handles; + for (const auto& loc: thingify(ls, R.cell.provider())) { auto opt_i = R.ion_location_index(ion, loc); if (!opt_i) continue; - R.result.push_back(fvm_probe_scalar{{xi.data() + *opt_i}, loc}); + handles.push_back(values.data() + *opt_i); + meta.push_back(loc); } + if (meta.empty()) return; + meta.shrink_to_fit(); + handles.shrink_to_fit(); + R.result.push_back( + fvm_probe_multi{ + .raw_handles=handles, + .metadata=meta, + }); +} + +template +void resolve_probe(const cable_probe_ion_int_concentration& p, probe_resolution_data& R) { + const auto& ion = p.ion; + if (!R.state->ion_data.count(ion)) return; + resolve_ion_conc_common(p.locations, ion, R.state->ion_data.at(ion).Xi_, R); } template void resolve_probe(const cable_probe_ion_ext_concentration& p, probe_resolution_data& R) { const auto& ion = p.ion; - const auto& xo = R.state->ion_data.at(ion).Xo_; - for (mlocation loc: thingify(p.locations, R.cell.provider())) { - auto opt_i = R.ion_location_index(ion, loc); - if (!opt_i) continue; - R.result.push_back(fvm_probe_scalar{{xo.data() + *opt_i}, loc}); - } + if (!R.state->ion_data.count(ion)) return; + resolve_ion_conc_common(p.locations, ion, R.state->ion_data.at(ion).Xo_, R); } template void resolve_probe(const cable_probe_ion_diff_concentration& p, probe_resolution_data& R) { const auto& ion = p.ion; - const auto& xd = R.state->ion_data.at(ion).Xd_; - for (mlocation loc: thingify(p.locations, R.cell.provider())) { - auto opt_i = R.ion_location_index(ion, loc); - if (!opt_i) continue; - R.result.push_back(fvm_probe_scalar{{xd.data() + *opt_i}, loc}); - } + if (!R.state->ion_data.count(ion)) return; + resolve_ion_conc_common(p.locations, ion, R.state->ion_data.at(ion).Xd_, R); } // Common implementation for int and ext concentrations across whole cell: template -void resolve_ion_conc_common(const std::vector& ion_cvs, const arb_value_type* src, probe_resolution_data& R) { +void resolve_ion_conc_cell_common(const std::vector& ion_cvs, + const arb_value_type* src, + probe_resolution_data& R) { fvm_probe_multi r; mcable_list cables; for (auto i: util::count_along(ion_cvs)) { @@ -1076,21 +1162,21 @@ template void resolve_probe(const cable_probe_ion_int_concentration_cell& p, probe_resolution_data& R) { if (!R.state->ion_data.count(p.ion)) return; if (R.state->ion_data.at(p.ion).Xi_.empty()) return; - resolve_ion_conc_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xi_.data(), R); + resolve_ion_conc_cell_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xi_.data(), R); } template void resolve_probe(const cable_probe_ion_ext_concentration_cell& p, probe_resolution_data& R) { if (!R.state->ion_data.count(p.ion)) return; if (R.state->ion_data.at(p.ion).Xo_.empty()) return; - resolve_ion_conc_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xo_.data(), R); + resolve_ion_conc_cell_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xo_.data(), R); } template void resolve_probe(const cable_probe_ion_diff_concentration_cell& p, probe_resolution_data& R) { if (!R.state->ion_data.count(p.ion)) return; if (R.state->ion_data.at(p.ion).Xd_.empty()) return; - resolve_ion_conc_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xd_.data(), R); + resolve_ion_conc_cell_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xd_.data(), R); } } // namespace arb diff --git a/arbor/include/arbor/arbexcept.hpp b/arbor/include/arbor/arbexcept.hpp index cf162c12b5..d00edf64a9 100644 --- a/arbor/include/arbor/arbexcept.hpp +++ b/arbor/include/arbor/arbexcept.hpp @@ -1,5 +1,6 @@ #pragma once +#include "arbor/util/any_ptr.hpp" #include #include #include @@ -106,6 +107,32 @@ struct ARB_SYMBOL_VISIBLE zero_thread_requested_error: arbor_exception { unsigned nbt; }; +// Sampling errors +template +struct ARB_SYMBOL_VISIBLE sample_reader_metadata_error: arbor_exception { + sample_reader_metadata_error(util::any_ptr have_): + arbor_exception{"Sample reader: could not cast to metadata type; expected " + expect + " got " + have}, + have{have_.type().name()}, + expect{typeid((E*)nullptr).name()} + {} + + std::string have; + std::string expect; +}; + +template +struct ARB_SYMBOL_VISIBLE sample_reader_value_error: arbor_exception { + sample_reader_value_error(std::any have_): + arbor_exception{"Sample reader: could not cast to value type; expected " + expect + " got " + have}, + have{have_.type().name()}, + expect{typeid((E*)nullptr).name()} + {} + + std::string have; + std::string expect; +}; + + // Domain decomposition errors: struct ARB_SYMBOL_VISIBLE gj_unsupported_domain_decomposition: arbor_exception { diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp index 3ef19cbcf3..6d695f71f1 100644 --- a/arbor/include/arbor/cable_cell.hpp +++ b/arbor/include/arbor/cable_cell.hpp @@ -21,36 +21,10 @@ #include #include #include +#include namespace arb { -// `cable_sample_range` is simply a pair of `const double*` pointers describing the sequence -// of double values associated with the cell-wide sample. - -using cable_sample_range = std::pair; - - -// Each kind of probe has its own type for representing its address, as below. -// -// Probe address specifications can be for _scalar_ data, associated with a fixed -// location or synapse on a cell, or _vector_ data, associated with multiple -// sites or sub-sections of a cell. -// -// Sampler functions receive an `any_ptr` to sampled data. The underlying pointer -// type is a const pointer to: -// * `double` for scalar data; -// * `cable_sample_range*` for vector data (see definition above). -// -// The metadata associated with a probe is also passed to a sampler via an `any_ptr`; -// the underlying pointer will be a const pointer to one of the following metadata types: -// * `mlocation` for most scalar queries; -// * `cable_probe_point_info` for point mechanism state queries; -// * `mcable_list` for most vector queries; -// * `std::vector` 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 ARB_SYMBOL_VISIBLE cable_probe_point_info { cell_tag_type target; // Target tag of point process instance on cell. @@ -59,70 +33,105 @@ struct ARB_SYMBOL_VISIBLE cable_probe_point_info { mlocation loc; // Point on cell morphology where instance is placed. }; -// Voltage estimate [mV] at `location`, interpolated. -// Sample value type: `double` -// Sample metadata type: `mlocation` +// Cable cell type definitions +using cable_sample_type = const double; + +using cable_state_meta_type = const mlocation; +using cable_state_cell_meta_type = const mcable; +using cable_point_meta_type = const cable_probe_point_info; + +template <> +struct probe_value_type_of { + using type = cable_sample_type; +}; + +template <> +struct probe_value_type_of { + using type = cable_sample_type; +}; + +template <> +struct probe_value_type_of { + using type = cable_sample_type; +}; + +// Each kind of probe has its own type for representing its address, as below. +// The metadata associated with a probe is also passed to a sampler via an `any_ptr`; +// the underlying pointer will be a const pointer to the associated metadata type. + +// Voltage estimate [mV] at `location`, +// interpolated. struct ARB_SYMBOL_VISIBLE cable_probe_membrane_voltage { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; locset locations; }; -// Voltage estimate [mV], reported against each cable in each control volume. Not interpolated. -// Sample value type: `cable_sample_range` -// Sample metadata type: `mcable_list` -struct ARB_SYMBOL_VISIBLE cable_probe_membrane_voltage_cell {}; +// Voltage estimate [mV], reported against each cable in each control volume. +// Not interpolated. +struct ARB_SYMBOL_VISIBLE cable_probe_membrane_voltage_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; +}; -// Axial current estimate [nA] at `location`, interpolated. -// Sample value type: `double` -// Sample metadata type: `mlocation` +// Axial current estimate [nA] at `location`, +// interpolated. struct ARB_SYMBOL_VISIBLE cable_probe_axial_current { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; locset locations; }; -// Total current density [A/m²] across membrane _excluding_ capacitive and stimulus current at `location`. -// Sample value type: `cable_sample_range` -// Sample metadata type: `mlocation` +// Total current density [A/m²] across membrane _excluding_ capacitive and +// stimulus current at `location`. struct ARB_SYMBOL_VISIBLE cable_probe_total_ion_current_density { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; locset locations; }; // Total ionic current [nA] across membrane _excluding_ capacitive current across components of the cell. // Sample value type: `cable_sample_range` // Sample metadata type: `mcable_list` -struct ARB_SYMBOL_VISIBLE cable_probe_total_ion_current_cell {}; +struct ARB_SYMBOL_VISIBLE cable_probe_total_ion_current_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; +}; // Total membrane current [nA] across components of the cell _excluding_ stimulus currents. -// Sample value type: `cable_sample_range` -// Sample metadata type: `mcable_list` -struct ARB_SYMBOL_VISIBLE cable_probe_total_current_cell {}; +struct ARB_SYMBOL_VISIBLE cable_probe_total_current_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; +}; // Stimulus currents [nA] across components of the cell. -// Sample value type: `cable_sample_range` -// Sample metadata type: `mcable_list` -struct ARB_SYMBOL_VISIBLE cable_probe_stimulus_current_cell {}; +struct ARB_SYMBOL_VISIBLE cable_probe_stimulus_current_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; +}; // Value of state variable `state` in density mechanism `mechanism` in CV at `location`. -// Sample value type: `double` -// Sample metadata type: `mlocation` struct ARB_SYMBOL_VISIBLE cable_probe_density_state { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; locset locations; std::string mechanism; std::string state; }; // Value of state variable `state` in density mechanism `mechanism` across components of the cell. -// -// Sample value type: `cable_sample_range` -// Sample metadata type: `mcable_list` struct ARB_SYMBOL_VISIBLE cable_probe_density_state_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; std::string mechanism; std::string state; }; // Value of state variable `key` in point mechanism `source` at target `target`. -// -// Sample value type: `double` -// Sample metadata type: `cable_probe_point_info` struct ARB_SYMBOL_VISIBLE cable_probe_point_state { + using value_type = cable_sample_type; + using meta_type = cable_point_meta_type; + cell_tag_type target; std::string mechanism; std::string state; @@ -142,71 +151,70 @@ struct ARB_SYMBOL_VISIBLE cable_probe_point_state { // Value of state variable `key` in point mechanism `source` at every target // with this mechanism. Metadata has one entry of type cable_probe_point_info // for each matched (possibly coalesced) instance. -// -// Sample value type: `cable_sample_range` -// Sample metadata type: `std::vector` struct ARB_SYMBOL_VISIBLE cable_probe_point_state_cell { + using value_type = cable_sample_type; + using meta_type = cable_point_meta_type; std::string mechanism; std::string state; }; // Current density [A/m²] across membrane attributed to the ion `source` at `location`. -// Sample value type: `double` -// Sample metadata type: `mlocation` struct ARB_SYMBOL_VISIBLE cable_probe_ion_current_density { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; locset locations; std::string ion; }; // Total ionic current [nA] attributed to the ion `source` across components of the cell. -// Sample value type: `cable_sample_range` -// Sample metadata type: `mcable_list` struct ARB_SYMBOL_VISIBLE cable_probe_ion_current_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; std::string ion; }; // Ionic internal concentration [mmol/L] of ion `source` at `location`. -// Sample value type: `double` -// Sample metadata type: `mlocation` struct ARB_SYMBOL_VISIBLE cable_probe_ion_int_concentration { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; locset locations; std::string ion; }; // Ionic internal concentration [mmol/L] of ion `source` across components of the cell. -// Sample value type: `cable_sample_range` -// Sample metadata type: `mcable_list` struct ARB_SYMBOL_VISIBLE cable_probe_ion_int_concentration_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; std::string ion; }; // Ionic diffusive concentration [mmol/L] of ion `source` at `location`. -// Sample value type: `double` -// Sample metadata type: `mlocation` struct ARB_SYMBOL_VISIBLE cable_probe_ion_diff_concentration { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; locset locations; std::string ion; }; -// Ionic diffusiev concentration [mmol/L] of ion `source` across components of the cell. -// Sample value type: `cable_sample_range` -// Sample metadata type: `mcable_list` +// Ionic diffusive concentration [mmol/L] of ion `source` across components of the cell. struct ARB_SYMBOL_VISIBLE cable_probe_ion_diff_concentration_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; std::string ion; }; // Ionic external concentration [mmol/L] of ion `source` at `location`. -// Sample value type: `double` -// Sample metadata type: `mlocation` struct ARB_SYMBOL_VISIBLE cable_probe_ion_ext_concentration { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; locset locations; std::string ion; }; // Ionic external concentration [mmol/L] of ion `source` across components of the cell. -// Sample value type: `cable_sample_range` -// Sample metadata type: `mcable_list` struct ARB_SYMBOL_VISIBLE cable_probe_ion_ext_concentration_cell { + using value_type = cable_sample_type; + using meta_type = cable_state_cell_meta_type; std::string ion; }; diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp index 621d36b6fc..f295bdff41 100644 --- a/arbor/include/arbor/cable_cell_param.hpp +++ b/arbor/include/arbor/cable_cell_param.hpp @@ -197,7 +197,7 @@ struct ARB_SYMBOL_VISIBLE ion_diffusivity { ion_diffusivity() = default; ion_diffusivity(const std::string& ion, const U::quantity& m, iexpr scale=1): ion{ion}, value(m.value_as(U::m2/U::s)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [m²/s]."}; + if (std::isnan(value)) throw std::domain_error{"Diffusivity " + ion + ": Value must be finite and in [m²/s]."}; } }; diff --git a/arbor/include/arbor/lif_cell.hpp b/arbor/include/arbor/lif_cell.hpp index 2af3c2bf7c..dbff193e11 100644 --- a/arbor/include/arbor/lif_cell.hpp +++ b/arbor/include/arbor/lif_cell.hpp @@ -5,6 +5,8 @@ #include #include +#include + namespace arb { namespace U = arb::units; @@ -28,8 +30,20 @@ struct ARB_SYMBOL_VISIBLE lif_cell { // LIF probe metadata, to be passed to sampler callbacks. Intentionally left blank. struct ARB_SYMBOL_VISIBLE lif_probe_metadata {}; +using lif_sample_type = const double; +using lif_meta_type = const lif_probe_metadata; + // Voltage estimate [mV]. // Sample value type: `double` -struct ARB_SYMBOL_VISIBLE lif_probe_voltage {}; +struct ARB_SYMBOL_VISIBLE lif_probe_voltage { + using value_type = lif_sample_type; + using meta_type = const lif_probe_metadata; +}; + +template <> +struct probe_value_type_of { + using type = lif_sample_type; +}; + } // namespace arb diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp index ebf952323b..2813172dd8 100644 --- a/arbor/include/arbor/morph/locset.hpp +++ b/arbor/include/arbor/morph/locset.hpp @@ -2,11 +2,8 @@ #include #include -#include #include -#include #include -#include #include #include @@ -118,6 +115,8 @@ namespace ls { // Explicit location on morphology. ARB_ARBOR_API locset location(msize_t branch, double pos); +inline +ARB_ARBOR_API locset location(mlocation loc) { return location(loc.branch, loc.pos); } // Set of terminal nodes on a morphology. ARB_ARBOR_API locset terminal(); diff --git a/arbor/include/arbor/sampling.hpp b/arbor/include/arbor/sampling.hpp index e5501e9fe0..f018fed84d 100644 --- a/arbor/include/arbor/sampling.hpp +++ b/arbor/include/arbor/sampling.hpp @@ -3,8 +3,11 @@ #include #include +#include #include #include +#include +#include namespace arb { @@ -35,23 +38,74 @@ struct one_tag { // User code is responsible for correctly determining the metadata type, // but the value of that metadata must be sufficient to determine the // correct interpretation of sample data provided to sampler callbacks. - struct probe_metadata { - cell_address_type id; // probe id - unsigned index; // index of probe source within those supplied by probe id - util::any_ptr meta; // probe-specific metadata + cell_address_type id; // probe id + unsigned index; // index of probe source within those supplied by probe id + std::size_t width; // width ie count of meta items + util::any_ptr meta; // probe-specific metadata +}; + +struct sample_records { + std::size_t n_sample = 0; // count of sample _rows_ + std::size_t width = 0; // count of sample _columns_ + const time_type* time = nullptr; // pointer to time data + std::any values; // resolves to pointer of probe-specific payload data D of layout D[n_sample][width] }; -struct sample_record { - time_type time; - util::any_ptr data; +template +struct sample_reader { + using value_type = probe_value_type_of_t; + using meta_type = M; + + std::size_t n_row() const { return n_sample_; } + std::size_t n_column() const { return width_; } + + // Retrieve sample value corresponding to + // - time=get_time(i) + // - location=get_metadata(j) + value_type value(std::size_t i, std::size_t j = 0) const { + arb_assert(i < n_sample_); + arb_assert(j < width_); + return values_[i*width_ + j]; + } + + time_type time(std::size_t i) const { + arb_assert(i < n_sample_); + return time_[i]; + } + + meta_type metadata(std::size_t j) const { + arb_assert(j < width_); + return metadata_[j]; + } + + sample_reader(util::any_ptr apm, + const sample_records& sr): + width_(sr.width), + n_sample_(sr.n_sample), + time_(sr.time) + { + using util::any_cast; + metadata_ = any_cast(apm); + if (!metadata_) throw sample_reader_metadata_error{apm}; + using V = sample_reader::value_type; + try { + values_ = any_cast(sr.values); + } + catch(const std::bad_any_cast& e) { + throw sample_reader_value_error{sr.values}; + } + } + + private: + std::size_t width_ = 0; + std::size_t n_sample_ = 0; + const time_type* time_ = nullptr; + value_type* values_ = nullptr; + meta_type* metadata_ = nullptr; }; -using sampler_function = std::function< - void (probe_metadata, - std::size_t, // number of sample records - const sample_record* // pointer to first sample record - )>; +using sampler_function = std::function; using sampler_association_handle = std::size_t; diff --git a/arbor/include/arbor/schedule.hpp b/arbor/include/arbor/schedule.hpp index 8856cceac2..e38f9e1100 100644 --- a/arbor/include/arbor/schedule.hpp +++ b/arbor/include/arbor/schedule.hpp @@ -9,7 +9,6 @@ #include #include -#include #include #include #include @@ -36,11 +35,11 @@ inline time_event_span as_time_event_span(const std::vector& v) { struct ARB_ARBOR_API schedule { schedule(); - template , schedule>>> + template , schedule>>> explicit schedule(const Impl& impl): impl_(new wrap(impl)) {} - template , schedule>>> + template , schedule>>> explicit schedule(Impl&& impl): impl_(new wrap(std::move(impl))) {} diff --git a/arbor/include/arbor/simple_sampler.hpp b/arbor/include/arbor/simple_sampler.hpp index 196f1f1110..ba87b80c69 100644 --- a/arbor/include/arbor/simple_sampler.hpp +++ b/arbor/include/arbor/simple_sampler.hpp @@ -1,175 +1,50 @@ #pragma once -/* - * Simple(st?) implementation of a recorder of scalar - * trace data from a cell probe, with some metadata. - */ - -#include #include #include #include +#include #include #include namespace arb { -template -struct trace_entry { - time_type t; - V v; -}; - -// `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 metadata. - -template -struct trace_data: private std::vector> { - Meta meta; - explicit operator bool() const { return !this->empty(); } - - using base = std::vector>; - 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 -struct trace_data: private std::vector> { - explicit operator bool() const { return !this->empty(); } - - using base = std::vector>; - 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 -struct trace_vector: private std::vector> { - const trace_data& get(std::size_t i) const { - return isize()? (*this)[i]: empty_trace; - } - - using base = std::vector>; - 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 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 -// be used to populate a trace of vectors. - -template -struct trace_push_back { - template - static bool push_back(trace_data& trace, const sample_record& rec) { - if (auto p = util::any_cast(rec.data)) { - trace.push_back({rec.time, *p}); - return true; - } - return false; - } -}; - -template -struct trace_push_back> { - template - static bool push_back(trace_data, Meta>& trace, const sample_record& rec) { - if (auto p = util::any_cast*>(rec.data)) { - trace.push_back({rec.time, *p}); - return true; - } - else if (auto p = util::any_cast*>(rec.data)) { - trace.push_back({rec.time, std::vector(p->first, p->second)}); - return true; - } - return false; - } -}; - -template -class simple_sampler { -public: - explicit simple_sampler(trace_vector& trace): trace_(trace) {} - - void operator()(probe_metadata pm, std::size_t n, const sample_record* recs) { - if constexpr (std::is_void_v) { - if (trace_.size()<=pm.index) { - trace_.resize(pm.index+1); - } - - for (std::size_t i = 0; i::push_back(trace_[pm.index], recs[i])) { - throw std::runtime_error("unexpected sample type in simple_sampler"); - } +// Simple(st?) implementation of a recorder of scalar trace data from a cell +// probe, with some metadata. + +template +struct simple_sampler_result { + using value_type = probe_value_type_of_t; + std::size_t n_sample = 0; + std::size_t width = 0; + std::vector time; + std::vector>> values; + std::vector> metadata; + + void from_reader(const sample_reader& reader) { + n_sample = reader.n_row(); + width = reader.n_column(); + values.resize(width); + for (std::size_t ix = 0ul; ix < reader.n_row(); ++ix) { + time.push_back(reader.time(ix)); + for (std::size_t iy = 0ul; iy < reader.n_column(); ++iy) { + auto v = reader.value(ix, iy); + values[iy].push_back(v); } } - else { - const Meta* m = util::any_cast(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::push_back(trace_[pm.index], recs[i])) { - throw std::runtime_error("unexpected sample type in simple_sampler"); - } - } + for (std::size_t iy = 0ul; iy < reader.n_column(); ++iy) { + metadata.push_back(reader.metadata(iy)); } } - -private: - trace_vector& trace_; }; -template -inline simple_sampler make_simple_sampler(trace_vector& trace) { - return simple_sampler(trace); +template +auto make_simple_sampler(simple_sampler_result& trace) { + return [&trace](const probe_metadata& pm, const sample_records& recs) { + auto reader = sample_reader(pm.meta, recs); + trace.from_reader(reader); + }; } } // namespace arb diff --git a/arbor/include/arbor/util/expected.hpp b/arbor/include/arbor/util/expected.hpp index 1af88eb42f..bb9907a020 100644 --- a/arbor/include/arbor/util/expected.hpp +++ b/arbor/include/arbor/util/expected.hpp @@ -14,8 +14,6 @@ #include #include -#include - namespace arb { namespace util { @@ -87,8 +85,8 @@ struct unexpected { template >, - typename = std::enable_if_t>>, - typename = std::enable_if_t>> + typename = std::enable_if_t>>, + typename = std::enable_if_t>> > explicit unexpected(F&& f): value_(std::forward(f)) {} @@ -205,9 +203,9 @@ struct expected { template < typename S, - typename = std::enable_if_t>>, - typename = std::enable_if_t>>, - typename = std::enable_if_t, remove_cvref_t>>, + typename = std::enable_if_t>>, + typename = std::enable_if_t>>, + typename = std::enable_if_t, std::remove_cvref_t>>, typename = std::enable_if_t> > expected(S&& x): data_(std::in_place_index<0>, std::forward(x)) {} @@ -232,7 +230,7 @@ struct expected { template < typename S, - typename = std::enable_if_t>>, + typename = std::enable_if_t>>, typename = std::enable_if_t>, typename = std::enable_if_t> > diff --git a/arbor/include/arbor/util/extra_traits.hpp b/arbor/include/arbor/util/extra_traits.hpp index 62653158f5..dbbafb5542 100644 --- a/arbor/include/arbor/util/extra_traits.hpp +++ b/arbor/include/arbor/util/extra_traits.hpp @@ -1,17 +1,15 @@ #pragma once namespace arb { -namespace util { -// TODO: C++20 replace with std::remove_cvref, std::remove_cvref_t - -template -struct remove_cvref { - typedef std::remove_cv_t> type; +// Helper class, to be specialized in each cell header, mapping from Metadata to Value types +template +struct probe_value_type_of { + using meta_type = M; + using type = void; }; -template -using remove_cvref_t = typename remove_cvref::type; +template +using probe_value_type_of_t = probe_value_type_of::type; -} // namespace util } // namespace arb diff --git a/arbor/include/arbor/util/unique_any.hpp b/arbor/include/arbor/util/unique_any.hpp index 9cbcf67d14..9c2b48b96c 100644 --- a/arbor/include/arbor/util/unique_any.hpp +++ b/arbor/include/arbor/util/unique_any.hpp @@ -7,7 +7,6 @@ #include #include -#include // A non copyable variant of std::any. The two main use cases for such a // container are: @@ -55,8 +54,8 @@ class ARB_SYMBOL_VISIBLE unique_any { template < typename T, - typename = std::enable_if_t, unique_any>>, - typename = std::enable_if_t, std::any>> + typename = std::enable_if_t, unique_any>>, + typename = std::enable_if_t, std::any>> > unique_any(T&& other) { state_.reset(new model>(std::forward(other))); @@ -69,8 +68,8 @@ class ARB_SYMBOL_VISIBLE unique_any { template < typename T, - typename = std::enable_if_t, unique_any>>, - typename = std::enable_if_t, std::any>> + typename = std::enable_if_t, unique_any>>, + typename = std::enable_if_t, std::any>> > unique_any& operator=(T&& other) { state_.reset(new model>(std::forward(other))); @@ -161,7 +160,7 @@ T* any_cast(unique_any* operand) { template T any_cast(const unique_any& operand) { - using U = remove_cvref_t; + using U = std::remove_cvref_t; static_assert(std::is_constructible::value, "any_cast type can't construct copy of contained object"); @@ -174,7 +173,7 @@ T any_cast(const unique_any& operand) { template T any_cast(unique_any& operand) { - using U = remove_cvref_t; + using U = std::remove_cvref_t; static_assert(std::is_constructible::value, "any_cast type can't construct copy of contained object"); @@ -187,7 +186,7 @@ T any_cast(unique_any& operand) { template T any_cast(unique_any&& operand) { - using U = remove_cvref_t; + using U = std::remove_cvref_t; static_assert(std::is_constructible::value, "any_cast type can't construct copy of contained object"); diff --git a/arbor/lif_cell_group.cpp b/arbor/lif_cell_group.cpp index 5af84c8874..f15b7ef780 100644 --- a/arbor/lif_cell_group.cpp +++ b/arbor/lif_cell_group.cpp @@ -99,6 +99,11 @@ lif_decay(const lif_lowered_cell& cell, double t0, double t1) { return (cell.V_m - cell.E_L)*exp((t0 - t1)/cell.tau_m) + cell.E_L; } +struct lif_samples { + std::vector times; + std::vector values; +}; + // Advances a single cell (lid) with the exact solution (jumps can be arbitrary). // Parameter dt is ignored, since we make jumps between two consecutive spikes. void lif_cell_group::advance_cell(time_type tfinal, @@ -115,7 +120,7 @@ void lif_cell_group::advance_cell(time_type tfinal, // collected sampling data std::unordered_map>> sampled; + lif_samples>> sampled; // samples to process std::size_t n_values = 0; std::vector> samples; @@ -135,7 +140,8 @@ void lif_cell_group::advance_cell(time_type tfinal, for (const auto& pid: assoc.probeset_ids) { if (pid.gid != gid) continue; arb_assert (0 == sampled[hdl].count(pid)); - sampled[hdl][pid].reserve(n_times); + sampled[hdl][pid].times.reserve(n_times); + sampled[hdl][pid].values.reserve(n_times); delta += n_times; } if (delta == 0) continue; @@ -146,10 +152,7 @@ void lif_cell_group::advance_cell(time_type tfinal, std::sort(samples.begin(), samples.end()); int n_samples = samples.size(); int sample_idx = 0; - // Now allocate some scratch space for the probed values, if we don't, - // re-alloc might move our data - std::vector sampled_voltages; - sampled_voltages.reserve(n_values); + // integrate until tfinal using the exact solution of membrane voltage differential equation. for (;;) { const auto event_time = event_idx < n_events ? event_lanes[lid][event_idx].time : tfinal; @@ -208,10 +211,9 @@ void lif_cell_group::advance_cell(time_type tfinal, // we are not in the refractory period, apply decay U = lif_decay(cell, t, time); } - // Store U for later use. - sampled_voltages.push_back(U); // Set up reference to sampled value - sampled[hdl][key].push_back(sample_record{time, {&sampled_voltages.back()}}); + sampled[hdl][key].times.push_back(time); + sampled[hdl][key].values.push_back(U); break; } default: @@ -226,7 +228,6 @@ void lif_cell_group::advance_cell(time_type tfinal, } last_time_updated_[lid] = t; } - arb_assert (sampled_voltages.size() <= n_values); // Now we need to call all sampler callbacks with the data we have collected { std::lock_guard guard(sampler_mex_); @@ -234,8 +235,10 @@ void lif_cell_group::advance_cell(time_type tfinal, const auto& fun = samplers_[k].sampler; for (auto& [id, us]: vs) { auto meta = get_probe_metadata(id)[0]; - fun(meta, us.size(), us.data()); - us.clear(); + fun(meta, sample_records{.n_sample=us.times.size(), + .width=1, + .time=us.times.data(), + .values=const_cast(us.values.data())}); } } } @@ -253,7 +256,7 @@ std::vector lif_cell_group::get_probe_metadata(const cell_addres // SAFETY: Probe associations are fixed after construction, so we do not // need to grab the mutex. if (probes_.count(key)) { - return {probe_metadata{key, 0, &probes_.at(key).metadata}}; + return {probe_metadata{key, 0, 1, &probes_.at(key).metadata}}; } else { return {}; } diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index 9945e25082..eb8ab95296 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -557,9 +557,7 @@ std::vector simulation_state::get_probe_metadata(const cell_addr if (auto linfo = util::value_by_key(gid_to_local_, probeset_id.gid)) { return cell_groups_.at(linfo->group_index)->get_probe_metadata(probeset_id); } - else { - return {}; - } + return {}; } // Simulation class implementations forward to implementation class. @@ -589,11 +587,9 @@ time_type simulation::run(const units::quantity& tfinal, const units::quantity& return impl_->run(tfinal_ms, dt_ms); } -sampler_association_handle simulation::add_sampler( - cell_member_predicate probeset_ids, - schedule sched, - sampler_function f) -{ +sampler_association_handle simulation::add_sampler(cell_member_predicate probeset_ids, + schedule sched, + sampler_function f) { return impl_->add_sampler(std::move(probeset_ids), std::move(sched), std::move(f)); } diff --git a/doc/cpp/probe_sample.rst b/doc/cpp/probe_sample.rst index 912ab99b42..ff6a8a6372 100644 --- a/doc/cpp/probe_sample.rst +++ b/doc/cpp/probe_sample.rst @@ -8,41 +8,22 @@ Cable cell probing and sampling Cable cell probes ----------------- -Various properties of a cable cell can be sampled. They fall into two classes: -scalar probes are associated with a single real value, such as a membrane -voltage or mechanism state value at a particular location; vector probes return -multiple values corresponding to a quantity sampled at a set of points on the -cell. - -The sample data associated with a cable cell probe will either be a ``double`` -for scalar probes, or a ``cable_sample_range`` describing a half-open range of -``double`` values: - -.. code:: - - using cable_sample_range = std::pair - The probe metadata passed to the sampler will be a const pointer to: -* ``mlocation`` for most scalar probes; - -* ``cable_probe_point_info`` for point mechanism state queries; +* ``cable_probe_point_info`` for point mechanism state queries, +* ``mcable`` for cell-wide probes, +* ``mlocation`` for probes on locsets; -* ``mcable_list`` for most vector queries; - -* ``std::vector`` for cell-wide point mechanism state queries. - -The type ``cable_probe_point_info`` holds metadata for a single target on a cell: +where the type ``cable_probe_point_info`` holds metadata for a point process +state variable .. code:: struct cable_probe_point_info { // Target number of point process instance on cell. cell_lid_type target; - // Number of combined instances at this site. unsigned multiplicity; - // Point on cell morphology where instance is placed. mlocation loc; }; @@ -50,18 +31,14 @@ The type ``cable_probe_point_info`` holds metadata for a single target on a cell Note that the ``multiplicity`` will always be 1 if synapse coalescing is disabled. +Value types will always be conforming to ``const double*``. + Cable cell probes that contingently do not correspond to a valid measurable quantity are ignored: samplers attached to them will receive no values. Mechanism state queries, however will throw a ``cable_cell_error`` exception at simulation initialization if the requested state variable does not exist on the mechanism. -Cable cell probeset 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-specific metadata found in the ``probe_metadata`` parameter. - Membrane voltage ^^^^^^^^^^^^^^^^ @@ -74,22 +51,18 @@ Membrane voltage Queries cell membrane potential at each site in ``locations``. * Sample value: ``double``. Membrane potential in millivolts. - * Metadata: ``mlocation``. Location of probe. - .. code:: struct cable_probe_membrane_voltage_cell {}; Queries cell membrane potential across the whole cell. -* Sample value: ``cable_sample_range``. Each value is the - average membrane potential in millivolts across an unbranched - component of the cell, as determined by the discretisation. - -* Metadata: ``mcable_list``. Each cable in the cable list describes - the unbranched component for the corresponding sample value. +* Sample value: the average membrane potential in millivolts across an + unbranched component of the cell, as determined by the discretisation. +* Metadata: ``mcable``. Each cable in the cable list describes the unbranched + component for the corresponding sample value. Axial current ^^^^^^^^^^^^^ @@ -103,11 +76,9 @@ Axial current Estimate intracellular current at each site in ``locations``, in the distal direction. -* Sample value: ``double``. Current in nanoamperes. - +* Sample value: Current in nanoamperes. * Metadata: ``mlocation``. Location as of probe. - Transmembrane current ^^^^^^^^^^^^^^^^^^^^^ @@ -121,11 +92,9 @@ Transmembrane current Membrane current density attributed to a particular ion at each site in ``locations``. -* Sample value: ``double``. Current density in amperes per square metre. - +* Sample value: Current density in amperes per square metre. * Metadata: ``mlocation``. Location of probe. - .. code:: struct cable_probe_ion_current_cell { @@ -134,13 +103,10 @@ each site in ``locations``. Membrane current attributed to a particular ion across components of the cell. -* Sample value: ``cable_sample_range``. Each value is the current in - nanoamperes across an unbranched component of the cell, as determined - by the discretisation. - -* Metadata: ``mcable_list``. Each cable in the cable list describes - the unbranched component for the corresponding sample value. - +* Sample value: the current in nanoamperes across an unbranched component of the + cell, as determined by the discretisation. +* Metadata: ``mcable``. Each cable in the cable list describes the unbranched + component for the corresponding sample value. .. code:: @@ -150,37 +116,31 @@ Membrane current attributed to a particular ion across components of the cell. Membrane current density at given locations _excluding_ capacitive currents. -* Sample value: ``double``. Current density in amperes per square metre. - +* Sample value: Current density in amperes per square metre. * Metadata: ``mlocation``. Location of probe. - .. code:: struct cable_probe_total_ion_current_cell {}; Membrane current _excluding_ capacitive currents and stimuli across components of the cell. -* Sample value: ``cable_sample_range``. Each value is the current in +* Sample value: the current in nanoamperes across an unbranched component of the cell, as determined by the discretisation. - -* Metadata: ``mcable_list``. Each cable in the cable list describes +* Metadata: ``mcable``. Each cable in the cable list describes the unbranched component for the corresponding sample value. - .. code:: struct cable_probe_total_current_cell {}; Total membrane current excluding current stimuli across components of the cell. -* Sample value: ``cable_sample_range``. Each value is the current in - nanoamperes across an unbranched component of the cell, as determined - by the discretisation. - -* Metadata: ``mcable_list``. Each cable in the cable list describes - the unbranched component for the corresponding sample value. +* Sample value: Each value is the current in nanoamperes across an unbranched + component of the cell, as determined by the discretisation. +* Metadata: ``mcable``. Each cable in the cable list describes the unbranched + component for the corresponding sample value. .. code:: @@ -188,13 +148,12 @@ Total membrane current excluding current stimuli across components of the cell. Total stimulus currents applied across components of the cell. -* Sample value: ``cable_sample_range``. Each value is the current in - nanoamperes across an unbranched component of the cell, as determined - by the discretisation. Components of CVs where no stimulus is present - will report a corresponding stimulus value of zero. - -* Metadata: ``mcable_list``. Each cable in the cable list describes - the unbranched component for the corresponding sample value. +* Sample value: Each value is the current in nanoamperes across an unbranched + component of the cell, as determined by the discretisation. Components of CVs + where no stimulus is present will report a corresponding stimulus value of + zero. +* Metadata: ``mcable``. Each cable in the cable list describes the unbranched + component for the corresponding sample value. Ion concentration ^^^^^^^^^^^^^^^^^ @@ -208,10 +167,8 @@ Ion concentration Ionic internal concentration of ion at each site in ``locations``. -* Sample value: ``double``. Ion concentration in millimoles per litre. - -* Metadata: ``mlocation``. Location of probe. - +* Sample value: Ion concentration in millimoles per litre. +* Metadata: ``mlocation``. Location of probe. .. code:: @@ -221,28 +178,23 @@ Ionic internal concentration of ion at each site in ``locations``. Ionic external concentration of ion across components of the cell. -* Sample value: ``cable_sample_range``. Each value is the concentration in - millimoles per lire across an unbranched component of the cell, as determined - by the discretisation. - -* Metadata: ``mcable_list``. Each cable in the cable list describes - the unbranched component for the corresponding sample value. - +* Sample value: the concentration in millimoles per lire across an unbranched + component of the cell, as determined by the discretisation. +* Metadata: ``mcable``. Each cable in the cable list describes the unbranched + component for the corresponding sample value. .. code:: struct cable_probe_ion_ext_concentration { - mlocation location; + locset location; std::string ion; }; Ionic external concentration of ion at each site in ``locations``. -* Sample value: ``double``. Ion concentration in millimoles per litre. - +* Sample value: Ion concentration in millimoles per litre. * Metadata: ``mlocation``. Location of probe. - .. code:: struct cable_probe_ion_ext_concentration_cell { @@ -251,14 +203,11 @@ Ionic external concentration of ion at each site in ``locations``. Ionic external concentration of ion across components of the cell. -* Sample value: ``cable_sample_range``. Each value is the concentration in - millimoles per lire across an unbranched component of the cell, as determined - by the discretisation. - -* Metadata: ``mcable_list``. Each cable in the cable list describes +* Sample value: the concentration in millimoles per litre across an unbranched + component of the cell, as determined by the discretisation. +* Metadata: ``mcable``. Each cable in the cable list describes the unbranched component for the corresponding sample value. - Ionic diffusion concrentration. .. code:: @@ -268,8 +217,13 @@ Ionic diffusion concrentration. std::string ion; }; -Diffusive ionic concentration of the given ``ion`` at the -sites specified by ``locations``. +Diffusive ionic concentration of the given ``ion`` at the sites specified by +``locations``. + +* Sample value: the concentration in millimoles per litre across an unbranched + component of the cell, as determined by the discretisation. +* Metadata: ``mlocation``. Each cable in the cable list describes + the unbranched component for the corresponding sample value. .. code:: @@ -277,8 +231,13 @@ sites specified by ``locations``. std::string ion; }; -Ionic diffusion concrentration attributed to a particular ``ion`` across CVs of the cell. +Ionic diffusion concrentration attributed to a particular ``ion`` across CVs of +the cell. +* Sample value: the concentration in millimoles per litre across an unbranched + component of the cell, as determined by the discretisation. +* Metadata: ``mcable``. Each cable in the cable list describes + the unbranched component for the corresponding sample value. Mechanism state ^^^^^^^^^^^^^^^ @@ -291,14 +250,11 @@ Mechanism state std::string state; }; +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. -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. - -* Metadata: ``mlocation``. Location as given in the probeset address. - +* Sample value: State variable value. +* Metadata: ``mlocation``. Location as given in the probeset address. .. code:: @@ -309,14 +265,12 @@ If the mechanism is not defined at a particular site, that site is ignored. Value of state variable in a density mechanism across components of the cell. -* Sample value: ``cable_sample_range``. State variable values from the - mechanism across unbranched components of the cell, as determined - by the discretisation and mechanism extent. - -* Metadata: ``mcable_list``. Each cable in the cable list describes +* Sample value: State variable values from the mechanism across unbranched + components of the cell, as determined by the discretisation and mechanism + extent. +* Metadata: ``mcable``. Each cable in the cable list describes the unbranched component for the corresponding sample value. - .. code:: struct cable_probe_point_state { @@ -328,11 +282,9 @@ Value of state variable in a density mechanism across components of the cell. Value of state variable in a point mechanism associated with the given target. If the mechanism is not associated with this target, the probe is ignored. -* Sample value: ``double``. State variable value. - +* Sample value: State variable value. * Metadata: ``cable_probe_point_info``. Target number, multiplicity and location. - .. code:: struct cable_probe_point_state_cell { @@ -343,13 +295,10 @@ If the mechanism is not associated with this target, the probe is ignored. Value of state variable in a point mechanism for each of the targets in the cell with which it is associated. -* Sample value: ``cable_sample_range``. State variable values at each associated +* Sample value: State variable values at each associated target. +* Metadata: ``cable_probe_point_info``. Target metadata for each associated target. -* Metadata: ``std::vector``. Target metadata for each - associated target. - - .. _sampling_api: Sampling API @@ -359,7 +308,6 @@ The new API replaces the flexible but irreducibly inefficient scheme where the next sample time for a sampling was determined by the return value of the sampler callback. - Definitions ^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -375,10 +323,8 @@ sampler schedule A function or function object that, given a time interval, returns a list of sample times within that interval. - - Probes -^^^^^^^^^^^^^^^^^^^^^^^^^^^ +^^^^^^ Probes are specified in the recipe objects that are used to initialize a simulation; the specification of the item or value that is subjected to a @@ -429,13 +375,20 @@ will be passed to a sampler function or function object: .. code-block:: cpp struct probe_metadata { - cell_address_type id; // probeset id - unsigned index; // index of probe source within those supplied by probeset id - util::any_ptr meta; // probe-specific metadata + cell_address_type id; // probeset id + unsigned index; // index of probe source within those supplied by probeset id + std::size_t width = 0; // count of sample _columns_ + util::any_ptr meta; // probe-specific metadata }; - using sampler_function = - std::function; + struct sample_records { + std::size_t n_sample; // count of sample _rows_ + std::size_t width; // count of sample _columns_ + const time_type* time; // pointer to time data + std::any values; // resolves to pointer of probe-specific data D[n_sample][width] + }; + + using sampler_function = std::function; where the parameters are respectively the probe metadata, the number of samples, and finally a pointer to the sequence of sample records. @@ -449,53 +402,78 @@ The ``any_ptr`` value in the metadata points to const probe-specific metadata; the type of the metadata will depend upon the probeset address specified in the ``probe_info`` provided by the recipe. -One ``sample_record`` struct contains one sample of the probe data at a -given simulation time point: +The raw data in ``values`` can --- given knowledge of the correct type +information --- be cast to the correct type ``const T*`` and read traversing in +order ``T[n_sample][width]``. Likewise, ``meta`` can be cast to the metadata +type ``const M*`` and yields an array ``M[width]``. -.. container:: api-code +Each probe type has type definitions for the associated value and metadata +types, e.g. + +.. container:: example-code .. code-block:: cpp - struct sample_record { - time_type time; // simulation time of sample - any_ptr data; // sample data - }; + struct cable_probe_membrane_voltage { + using value_type = cable_sample_type; + using meta_type = cable_state_meta_type; + locset locations; + }; -The ``data`` field points to the sample data, wrapped in ``any_ptr`` for -type-checked access. The exact representation will depend on the nature of -the object that is being probed, but it should depend only on the cell type and -probeset address. +Access is made much more convenient through ``sample_reader``, see next section. -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, assuming -one probe per probeset id, might be as follows: +Sample Data Access +^^^^^^^^^^^^^^^^^^ + +The ``sample_reader`` provides a convenient way of accessing data retrieved in a +sampler callback, taking care of casting and the data layout. It can be used as +follows, provided the probe is known .. container:: example-code .. code-block:: cpp - using sample_data = std::map>>; + // This is the probe type we will attach to + using probe_type = cable_probe_membrane_voltage_cell; + + // This is the callback to attach + void callback(const probe_metadata& pm, const sample_records& recs) { + auto reader = sample_reader(pm.meta, recs); - struct scalar_sampler { - sample_data& samples; + for (std::size_t ix = 0ul; ix < reader.n_row(); ++ix) { + auto time = reader.time(ix); + for (std::size_t iy = 0ul; iy < reader.n_column(); ++iy) { + auto value = reader.value(ix, iy); + auto cable = reader.metadata(iy); + // ... use time, cable, value ... + } + } - explicit scalar_sample(sample_data& samples): samples(samples) {} +In general, it provides safe access to the raw samples, time, and metadata and allows +treating ``sample_records`` like tabular data with ``width`` columns containing the +``metadata`` and ``n_sample`` rows containing ``time`` and ``values``. - void operator()(probe_metadata pm, size_t n, const sample_record* records) { - for (size_t i=0; i(rec.data); - assert(data); - samples[pm.id].emplace_back(rec.time, *data); - } - } - }; + .. code-block:: cpp + + template + struct sample_reader { + using meta_type = M; + using value_type = probe_value_type_of_t; -The use of ``any_ptr`` allows type-checked access to the sample data, which -may differ in type from probe to probe. + std::size_t n_row() const { return n_sample_; } + std::size_t n_column() const { return width_; } + // Retrieve sample value corresponding to + // - time=time(i) + // - location=metadata(j) + value_type value(std::size_t i, std::size_t j = 0) const; + // Retrieve i'th time + time_type time(std::size_t i) const; + // Retrieve metadata at j + meta_type metadata(std::size_t j) const; + }; Model and cell group interface ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -510,10 +488,9 @@ Polling rates and sampler functions are set through the using sampler_association_handle = std::size_t; using cell_member_predicate = std::function; - sampler_association_handle simulation::add_sampler( - cell_member_predicate probeset_ids, - schedule sched, - sampler_function fn) + sampler_association_handle simulation::add_sampler(cell_member_predicate probeset_ids, + schedule sched, + sampler_function fn) void simulation::remove_sampler(sampler_association_handle); @@ -572,13 +549,12 @@ corresponding to the current state of the cell group. It should be expected that this difference in time should be no greater the the duration of the integration period (i.e. ``mindelay/2``). - Schedules ^^^^^^^^^ -Schedules represent a non-negative, monotonically increasing sequence -of time points, and are used to specify the sampling schedule in any -given association of a sampler function to a set of probes. +Schedules represent a non-negative, monotonically increasing sequence of time +points, and are used to specify the sampling schedule in any given association +of a sampler function to a set of probes. A ``schedule`` object has two methods: @@ -591,17 +567,17 @@ A ``schedule`` object has two methods: time_event_span events(time_type t0, time_type t1) A ``time_event_span`` is a ``std::pair`` of pointers `const time_type*`, -representing a view into an internally maintained collection of generated -time values. +representing a view into an internally maintained collection of generated time +values. -The ``events(t0, t1)`` method returns a view of monotonically -increasing time values in the half-open interval ``[t0, t1)``. -Successive calls to ``events`` — without an intervening call to ``reset()`` -— must request strictly subsequent intervals. +The ``events(t0, t1)`` method returns a view of monotonically increasing time +values in the half-open interval ``[t0, t1)``. Successive calls to ``events`` — +without an intervening call to ``reset()`` — must request strictly subsequent +intervals. -The data represented by the returned ``time_event_span`` view is valid -for the lifetime of the ``schedule`` object, and is invalidated by any -subsequent call to ``reset()`` or ``events()``. +The data represented by the returned ``time_event_span`` view is valid for the +lifetime of the ``schedule`` object, and is invalidated by any subsequent call +to ``reset()`` or ``events()``. The ``reset()`` method resets the state such that events can be retrieved from again from time zero. A schedule that is reset must then produce @@ -609,9 +585,9 @@ the same sequence of time points, that is, it must exhibit repeatable and deterministic behaviour. The ``schedule`` object itself uses type-erasure to wrap any schedule -implementation class, which can be any copy--constructible class that -provides the methods ``reset()`` and ``events(t0, t1)`` above. Three -schedule implementations are provided by the engine: +implementation class, which can be any copy--constructible class that provides +the methods ``reset()`` and ``events(t0, t1)`` above. Three schedule +implementations are provided by the engine: .. container:: api-code @@ -638,37 +614,39 @@ The ``simulation`` and ``cable_cell_group`` classes use classes defined in ``scheduler_map.hpp`` to simplify the management of sampler--probe associations and probe metadata. -``sampler_association_map`` wraps an ``unordered_map`` between sampler association -handles and tuples (*schedule*, *sampler*, *probe set*), with thread-safe -accessors. - +``sampler_association_map`` wraps an ``unordered_map`` between sampler +association handles and tuples (*schedule*, *sampler*, *probe set*), with +thread-safe accessors. Batched sampling in ``cable_cell_group`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The ``fvm_multicell`` implementations for CPU and GPU simulation of multi-compartment -cable neurons perform sampling in a batched manner: when their integration is -initialized, they take a sequence of ``sample_event`` objects which are used to -populate an implementation-specific ``event_stream`` that describes for each -cell the sample times and what to sample over the integration interval. - -When an integration step for a cell covers a sample event on that cell, the sample -is satisfied with the value from the cell state at the beginning of the time step, -after any postsynaptic spike events have been delivered. - -It is the responsibility of the ``cable_cell_group::advance()`` method to create the sample -events from the entries of its ``sampler_association_map``, and to dispatch the -sampled values to the sampler callbacks after the integration is complete. -Given an association tuple (*schedule*, *sampler*, *probe set*) where the *schedule* -has (non-zero) *n* sample times in the current integration interval, the ``cable_cell_group`` will -call the *sampler* callback once for probe in *probe set*, with *n* sample values. +The ``fvm_multicell`` implementations for CPU and GPU simulation of +multi-compartment cable neurons perform sampling in a batched manner: when their +integration is initialized, they take a sequence of ``sample_event`` objects +which are used to populate an implementation-specific ``event_stream`` that +describes for each cell the sample times and what to sample over the integration +interval. + +When an integration step for a cell covers a sample event on that cell, the +sample is satisfied with the value from the cell state at the beginning of the +time step, after any postsynaptic spike events have been delivered. + +It is the responsibility of the ``cable_cell_group::advance()`` method to create +the sample events from the entries of its ``sampler_association_map``, and to +dispatch the sampled values to the sampler callbacks after the integration is +complete. Given an association tuple (*schedule*, *sampler*, *probe set*) where +the *schedule* has (non-zero) *n* sample times in the current integration +interval, the ``cable_cell_group`` will call the *sampler* callback once for +probe in *probe set*, with *n* sample values. .. note:: - When the time values returned by a call to a schedule's ``events(t0, t1)`` method do not - perfectly coincide with the boundaries of the numerical time step grid, :math:`[t_0, t_0 + dt, - t_0 + 2\, dt, \, \cdots \, , t_1)`, the samples will be taken at the closest possible point in - time. In particular, any sample times :math:`t_s \in \left( t_i - dt/2,~ t_i + dt/2\right]` are + When the time values returned by a call to a schedule's ``events(t0, t1)`` + method do not perfectly coincide with the boundaries of the numerical time + step grid, :math:`[t_0, t_0 + dt, t_0 + 2\, dt, \, \cdots \, , t_1)`, the + samples will be taken at the closest possible point in time. In particular, + any sample times :math:`t_s \in \left( t_i - dt/2,~ t_i + dt/2\right]` are attributed to simulation time step :math:`t_i = t_0 + i\,dt`. @@ -685,5 +663,4 @@ Membrane voltage Queries cell membrane potential. * Sample value: ``double``. Membrane potential (mV). - * Metadata: none diff --git a/example/bench/CMakeLists.txt b/example/bench/CMakeLists.txt index b9a7e1205a..5a77e72a79 100644 --- a/example/bench/CMakeLists.txt +++ b/example/bench/CMakeLists.txt @@ -1,4 +1,4 @@ add_executable(bench EXCLUDE_FROM_ALL bench.cpp) add_dependencies(examples bench) -target_link_libraries(bench PRIVATE arbor arborenv arbor-sup ${json_library_name}) +target_link_libraries(bench PRIVATE arbor arborenv arbor-sup ${json_library_name} fmt::fmt-header-only) diff --git a/example/brunel/CMakeLists.txt b/example/brunel/CMakeLists.txt index 7a20c7a1f8..e7fe8b7760 100644 --- a/example/brunel/CMakeLists.txt +++ b/example/brunel/CMakeLists.txt @@ -1,4 +1,4 @@ add_executable(brunel EXCLUDE_FROM_ALL brunel.cpp) add_dependencies(examples brunel) -target_link_libraries(brunel PRIVATE arbor arborenv arbor-sup ext-tinyopt) +target_link_libraries(brunel PRIVATE arbor arborenv arbor-sup ext-tinyopt fmt::fmt-header-only) diff --git a/example/busyring/ring.cpp b/example/busyring/ring.cpp index b4903f84c9..91d6522b32 100644 --- a/example/busyring/ring.cpp +++ b/example/busyring/ring.cpp @@ -40,21 +40,22 @@ using arb::cable_probe_membrane_voltage; using namespace arborio::literals; namespace U = arb::units; +// result of simple sampler for probe type +using sample_result = arb::simple_sampler_result; + // Writes voltage trace as a json file. -void write_trace_json(std::string fname, const arb::trace_data& trace); +void write_trace_json(const std::string& path, const sample_result&); // Generate a cell. arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params); arb::cable_cell complex_cell(arb::cell_gid_type gid, const cell_parameters& params); -class ring_recipe: public arb::recipe { -public: +struct ring_recipe: public arb::recipe { ring_recipe(ring_params params): num_cells_(params.num_cells), min_delay_(params.min_delay), event_weight_(params.event_weight), - params_(params) - { + params_(params) { gprop.default_parameters = arb::neuron_parameter_defaults; gprop.catalogue.extend(arb::global_allen_catalogue()); @@ -70,24 +71,22 @@ class ring_recipe: public arb::recipe { cell_size_type num_cells() const override { return num_cells_; } cell_kind get_cell_kind(cell_gid_type gid) const override { return cell_kind::cable; } arb::util::unique_any get_cell_description(cell_gid_type gid) const override { - if (params_.cell.complex_cell) { - return complex_cell(gid, params_.cell); - } + if (params_.cell.complex_cell) return complex_cell(gid, params_.cell); return branch_cell(gid, params_.cell); } // Each cell has one incoming connection, from cell with gid-1, // and fan_in-1 random connections with very low weight. std::vector connections_on(cell_gid_type gid) const override { - std::vector cons; const auto ncons = params_.cell.synapses; - cons.reserve(ncons); - const auto s = params_.ring_size; const auto group = gid/s; const auto group_start = s*group; const auto group_end = std::min(group_start+s, num_cells_); - cell_gid_type src = gid==group_start? group_end-1: gid-1; + const auto src = gid == group_start ? group_end-1: gid-1; + + std::vector cons; + cons.reserve(ncons); cons.push_back(arb::cell_connection({src, "d"}, {"p"}, event_weight_, min_delay_*U::ms)); // Used to pick source cell for a connection. @@ -95,14 +94,13 @@ class ring_recipe: public arb::recipe { // Used to pick delay for a connection. std::uniform_real_distribution delay_dist(0, 2*min_delay_); auto src_gen = std::mt19937(gid); - for (unsigned i=1; i event_generators(cell_gid_type gid) const override { - if (gid%params_.ring_size == 0) { - return {arb::explicit_generator_from_milliseconds({"p"}, event_weight_, std::vector{1.0})}; - } else { - return {}; - } + if (gid % params_.ring_size == 0) return {arb::explicit_generator_from_milliseconds({"p"}, event_weight_, std::vector{1.0})}; + return {}; } - std::vector get_probes(cell_gid_type gid) const override { - // Measure at the soma. - arb::mlocation loc{0, 0.0}; - return {{cable_probe_membrane_voltage{loc}, "Um"}}; - } + // Measure at the soma. + std::vector get_probes(cell_gid_type gid) const override { return {{cable_probe_membrane_voltage{arb::mlocation {0, 0.0}}, "Um"}}; } private: cell_size_type num_cells_; @@ -222,8 +214,8 @@ int main(int argc, char** argv) { // Set up the probe that will measure voltage in the cell. - // This is where the voltage samples will be stored as (time, value) pairs - arb::trace_vector voltage; + // This is where the voltage samples will be stored + sample_result voltage; if (params.record_voltage) { // The id of the only probe on the cell: // the cell_member type points to (cell 0, probe 0) @@ -276,10 +268,7 @@ int main(int argc, char** argv) { } // Write the samples to a json file samples were stored on this rank. - if (voltage.size()>0u) { - std::string fname = params.odir + "/" + params.name + "_voltages.json"; - write_trace_json(fname, voltage.at(0)); - } + write_trace_json(params.odir + "/" + params.name + "_voltages.json", voltage); auto report = arb::profile::make_meter_report(meters, context); if (root) { @@ -295,22 +284,19 @@ int main(int argc, char** argv) { return 0; } -void write_trace_json(std::string fname, const arb::trace_data& trace) { +void write_trace_json(const std::string& path, const sample_result& result) { nlohmann::json json; - json["name"] = "ring demo"; + json["name"] = "busyring_demo"; json["units"] = "mV"; - json["cell"] = "0.0"; - json["probe"] = "0"; - - auto& jt = json["data"]["time"]; - auto& jy = json["data"]["voltage"]; - - for (const auto& sample: trace) { - jt.push_back(sample.t); - jy.push_back(sample.v); - } - - std::ofstream file(fname); + json["cell"] = "0"; + json["probe"] = "Um"; + std::stringstream loc; + loc << result.metadata.at(0); + json["location"] = loc.str(); + json["data"]["time"] = result.time; + json["data"]["voltage"] = result.values.at(0); + + std::ofstream file(path); file << std::setw(1) << json << "\n"; } diff --git a/example/diffusion/diffusion.cpp b/example/diffusion/diffusion.cpp index 80cb3fe34d..39831439a6 100644 --- a/example/diffusion/diffusion.cpp +++ b/example/diffusion/diffusion.cpp @@ -44,8 +44,8 @@ struct linear: public recipe { // Setup decor decor; decor.set_default(init_int_concentration{"na", i*U::mM}); - decor.set_default(ion_diffusivity{"na", b*U::mM}); - decor.paint("(tag 1)"_reg, ion_diffusivity{"na", b*U::mM}); + decor.set_default(ion_diffusivity{"na", b*U::m2/U::s}); + decor.paint("(tag 1)"_reg, ion_diffusivity{"na", b*U::m2/U::s}); decor.place("(location 0 0.5)"_ls, synapse("inject/x=bla", {{"alpha", 200.0*l}}), "Zap"); decor.paint("(all)"_reg, density("decay/x=bla")); return cable_cell({tree}, decor); @@ -57,17 +57,15 @@ struct linear: public recipe { std::ofstream out; -void sampler(probe_metadata pm, std::size_t n, const sample_record* samples) { - auto ptr = util::any_cast(pm.meta); - assert(ptr); - auto n_cable = ptr->size(); - out << "time,prox,dist,Xd\n" - << std::fixed << std::setprecision(4); - for (std::size_t i = 0; i(samples[i].data); - for (unsigned j = 0; j(pm.meta, samples); + for (std::size_t ix= 0; ix < reader.n_row(); ++ix) { + auto time = reader.time(ix); + for (std::size_t iy = 0; iy < reader.n_column(); ++iy) { + auto loc = reader.metadata(iy); + auto val = reader.value(ix, iy); + out << time << ',' << loc.prox_pos << ',' << loc.dist_pos << ',' << val << '\n'; } } out << '\n'; diff --git a/example/drybench/CMakeLists.txt b/example/drybench/CMakeLists.txt index f9cecd75c5..68f8f50c0e 100644 --- a/example/drybench/CMakeLists.txt +++ b/example/drybench/CMakeLists.txt @@ -1,4 +1,4 @@ add_executable(drybench EXCLUDE_FROM_ALL drybench.cpp) add_dependencies(examples drybench) -target_link_libraries(drybench PRIVATE arbor arborenv arbor-sup ${json_library_name}) +target_link_libraries(drybench PRIVATE arbor arborenv arbor-sup ${json_library_name} fmt::fmt-header-only) diff --git a/example/dryrun/CMakeLists.txt b/example/dryrun/CMakeLists.txt index e36fe972b9..3b8b932275 100644 --- a/example/dryrun/CMakeLists.txt +++ b/example/dryrun/CMakeLists.txt @@ -1,4 +1,4 @@ add_executable(dryrun EXCLUDE_FROM_ALL dryrun.cpp) add_dependencies(examples dryrun) -target_link_libraries(dryrun PRIVATE arbor arborio arborenv arbor-sup ${json_library_name}) +target_link_libraries(dryrun PRIVATE arbor arborio arborenv arbor-sup ${json_library_name} fmt::fmt-header-only) diff --git a/example/dryrun/dryrun.cpp b/example/dryrun/dryrun.cpp index b49b651f07..30d1b281b8 100644 --- a/example/dryrun/dryrun.cpp +++ b/example/dryrun/dryrun.cpp @@ -45,7 +45,13 @@ struct run_params { bool defaulted = true; }; -void write_trace_json(const arb::trace_data& trace); + +// result of simple sampler for probe type +using sample_result = arb::simple_sampler_result; + +// Writes voltage trace as a json file. +void write_trace_json(const sample_result&); + run_params read_options(int argc, char** argv); using arb::cell_gid_type; @@ -57,8 +63,7 @@ using arb::time_type; // Generate a cell. arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params); -class tile_desc: public arb::tile { -public: +struct tile_desc: public arb::tile { tile_desc(unsigned num_cells, unsigned num_tiles, cell_parameters params, unsigned min_delay): num_cells_(num_cells), num_tiles_(num_tiles), @@ -66,21 +71,13 @@ class tile_desc: public arb::tile { min_delay_(min_delay) {} - cell_size_type num_cells() const override { - return num_cells_; - } + cell_size_type num_cells() const override { return num_cells_; } - cell_size_type num_tiles() const override { - return num_tiles_; - } + cell_size_type num_tiles() const override { return num_tiles_; } - arb::util::unique_any get_cell_description(cell_gid_type gid) const override { - return branch_cell(gid, cell_params_); - } + arb::util::unique_any get_cell_description(cell_gid_type gid) const override { return branch_cell(gid, cell_params_); } - cell_kind get_cell_kind(cell_gid_type gid) const override { - return cell_kind::cable; - } + cell_kind get_cell_kind(cell_gid_type gid) const override { return cell_kind::cable; } std::any get_global_properties(arb::cell_kind) const override { arb::cable_cell_global_properties gprop; @@ -92,11 +89,9 @@ class tile_desc: public arb::tile { // src gid in {0, ..., num_cells_*num_tiles_ - 1}. std::vector connections_on(cell_gid_type gid) const override { std::uniform_int_distribution source_distribution(0, num_cells_*num_tiles_ - 2); - auto src_gen = std::mt19937(gid); auto src = source_distribution(src_gen); if (src>=gid) ++src; - return {arb::cell_connection({src, "detector"}, {"synapse"}, event_weight_, min_delay_*U::ms)}; } @@ -104,17 +99,12 @@ class tile_desc: public arb::tile { // for ALL cells on ALL ranks. This is because the symmetric recipe can not easily // translate the src gid of an event generator std::vector event_generators(cell_gid_type gid) const override { - std::vector gens; - if (gid%20 == 0) { - gens.push_back(arb::explicit_generator_from_milliseconds({"synapse"}, event_weight_, std::vector{1.0})); - } - return gens; + if (gid%20 == 0) return { arb::explicit_generator_from_milliseconds({"synapse"}, event_weight_, std::vector{1.0}) }; + return {}; } - std::vector 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}}, "Um"}}; - } + // One probe per cell, sampling membrane voltage at end of soma. + std::vector get_probes(cell_gid_type gid) const override { return {{arb::cable_probe_membrane_voltage{arb::mlocation{0, 0.0}}, "Um"}}; } private: cell_size_type num_cells_; @@ -152,11 +142,11 @@ int main(int argc, char** argv) { std::cout << sup::mask_stream(root); // Print a banner with information about hardware configuration - std::cout << "gpu: " << (has_gpu(ctx)? "yes": "no") << "\n"; - std::cout << "threads: " << num_threads(ctx) << "\n"; - std::cout << "mpi: " << (has_mpi(ctx)? "yes": "no") << "\n"; - std::cout << "ranks: " << num_ranks(ctx) << "(" << params.num_ranks << ")\n" << std::endl; - std::cout << "run mode: " << distribution_type(ctx) << "\n"; + std::cout << "gpu: " << (has_gpu(ctx)? "yes": "no") << "\n" + << "threads: " << num_threads(ctx) << "\n" + << "mpi: " << (has_mpi(ctx)? "yes": "no") << "\n" + << "ranks: " << num_ranks(ctx) << "(" << params.num_ranks << ")\n" << std::endl + << "run mode: " << distribution_type(ctx) << "\n"; assert(arb::num_ranks(ctx)==params.num_ranks); @@ -176,7 +166,7 @@ int main(int argc, char** argv) { // The schedule for sampling every 1 ms. auto sched = arb::regular_schedule(1*arb::units::ms); // This is where the voltage samples will be stored as (time, value) pairs - arb::trace_vector voltage; + sample_result voltage; // Now attach the sampler at probeset_id, with sampling schedule sched, writing to voltage sim.add_sampler(arb::one_probe(probeset_id), sched, arb::make_simple_sampler(voltage)); @@ -216,7 +206,7 @@ int main(int argc, char** argv) { } } // Write the samples to a json file. - write_trace_json(voltage.at(0)); + write_trace_json(voltage); } auto profile = arb::profile::profiler_summary(); @@ -233,27 +223,26 @@ int main(int argc, char** argv) { return 0; } -void write_trace_json(const arb::trace_data& trace) { +void write_trace_json(const sample_result& result) { std::string path = "./voltages.json"; nlohmann::json json; - json["name"] = "ring demo"; + json["name"] = "ring_demo"; json["units"] = "mV"; - json["cell"] = "0.0"; - json["probe"] = "0"; - - auto& jt = json["data"]["time"]; - auto& jy = json["data"]["voltage"]; - - for (const auto& sample: trace) { - jt.push_back(sample.t); - jy.push_back(sample.v); - } + json["cell"] = "0"; + json["probe"] = "Um"; + std::stringstream loc; + loc << result.metadata.at(0); + json["location"] = loc.str(); + json["data"]["time"] = result.time; + json["data"]["voltage"] = result.values.at(0); std::ofstream file(path); file << std::setw(1) << json << "\n"; } + + run_params read_options(int argc, char** argv) { using sup::param_from_json; diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp index bf20ad6fac..96d71ed8ac 100644 --- a/example/gap_junctions/gap_junctions.cpp +++ b/example/gap_junctions/gap_junctions.cpp @@ -1,13 +1,11 @@ -/* - * A miniapp that demonstrates how to make a model with gap junctions - * - */ +// A miniapp that demonstrates how to make a model with gap junctions #include #include #include #include +#include #include #include @@ -60,8 +58,11 @@ using arb::cell_size_type; using arb::cell_kind; using arb::time_type; +using probe_t = arb::cable_probe_membrane_voltage; +using sample_results = std::vector>; + // Writes voltage trace as a json file. -void write_trace_json(const std::vector>& trace, unsigned rank); +void write_trace_json(const sample_results& traces, unsigned rank); // Generate a cell. arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncells, double stim_duration); @@ -141,11 +142,7 @@ int main(int argc, char** argv) { unsigned nt = arbenv::default_concurrency(); int gpu_id = arbenv::find_private_gpu(MPI_COMM_WORLD); auto context = arb::make_context(arb::proc_allocation{nt, gpu_id}, MPI_COMM_WORLD); - { - int rank; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - root = rank==0; - } + root = arb::rank(context) == 0; #else auto context = arb::make_context(arbenv::default_allocation()); #endif @@ -156,11 +153,14 @@ int main(int argc, char** argv) { std::cout << sup::mask_stream(root); + auto b2a = [](bool c) { return c ; }; + // Print a banner with information about hardware configuration - std::cout << "gpu: " << (has_gpu(context)? "yes": "no") << "\n"; - std::cout << "threads: " << num_threads(context) << "\n"; - std::cout << "mpi: " << (has_mpi(context)? "yes": "no") << "\n"; - std::cout << "ranks: " << num_ranks(context) << "\n" << std::endl; + std::cout << fmt::format("gpu: {}\n" + "threads: {}\n" + "mpi: {}\n" + "ranks: {}\n", + b2a(has_gpu(context)), num_threads(context), b2a(has_mpi(context)), num_ranks(context)); auto params = read_options(argc, argv); @@ -179,12 +179,14 @@ int main(int argc, char** argv) { auto sched = arb::regular_schedule(0.025*U::ms); // This is where the voltage samples will be stored as (time, value) pairs - std::vector> voltage_traces(decomp.num_local_cells()); + + sample_results voltage_traces(decomp.num_local_cells()); // Now attach the sampler at probeset_id, with sampling schedule sched, writing to voltage - unsigned j=0; - for (auto g : decomp.groups()) { - for (auto i : g.gids) { + unsigned j = 0; + + for (const auto& g: decomp.groups()) { + for (auto i: g.gids) { sim.add_sampler(arb::one_probe({i, "Um"}), sched, arb::make_simple_sampler(voltage_traces[j++])); } } @@ -243,24 +245,22 @@ int main(int argc, char** argv) { return 0; } -void write_trace_json(const std::vector>& traces, unsigned rank) { +void write_trace_json(const sample_results& traces, unsigned rank) { for (unsigned i = 0; i < traces.size(); i++) { - std::string path = "./voltages_" + std::to_string(rank) + - "_" + std::to_string(i) + ".json"; + std::string path = fmt::format("./voltages_{}_{}.json", rank, i); nlohmann::json json; - json["name"] = "gj demo: cell " + std::to_string(i); + json["name"] = fmt::format("gj demo: cell {}", i); json["units"] = "mV"; - json["cell"] = std::to_string(i); - json["group"] = std::to_string(rank); + json["cell"] = i; + json["group"] = rank; json["probe"] = "0"; - - auto &jt = json["data"]["time"]; - auto &jy = json["data"]["voltage"]; - - for (const auto &sample: traces[i].at(0)) { - jt.push_back(sample.t); - jy.push_back(sample.v); + json["data"]["time"] = traces[i].time; + for (std::size_t ix = 0; ix < traces[i].width; ++ix) { + std::stringstream ss; + ss << traces[i].metadata[ix]; + json["data"]["meta"][ix] = ss.str(); + json["data"]["voltage"][ix] = traces[i].values[ix]; } std::ofstream file(path); @@ -318,9 +318,7 @@ gap_params read_options(int argc, char** argv) { std::cout << "Loading parameters from file: " << fname << "\n"; std::ifstream f(fname); - if (!f.good()) { - throw std::runtime_error("Unable to open input parameter file: "+fname); - } + if (!f.good()) throw std::runtime_error("Unable to open input parameter file: "+fname); nlohmann::json json; f >> json; diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp index b0600175b9..313a394eeb 100644 --- a/example/generators/generators.cpp +++ b/example/generators/generators.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -33,8 +34,11 @@ using arb::time_type; using namespace arborio::literals; +// result of simple sampler for probe type +using sample_result = arb::simple_sampler_result; + // Writes voltage trace as a json file. -void write_trace_json(const arb::trace_data& trace); +void write_trace_json(const sample_result&); class generator_recipe: public arb::recipe { public: @@ -62,7 +66,7 @@ class generator_recipe: public arb::recipe { // Add one synapse at the soma. // This synapse will be the target for all events, from both // event_generators. - .place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "syn"); + .place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "syn"); return arb::cable_cell(tree, decor, labels); } @@ -136,7 +140,7 @@ int main() { // The schedule for sampling is 10 samples every 1 ms. auto sched = arb::regular_schedule(0.1*arb::units::ms); // This is where the voltage samples will be stored as (time, value) pairs - arb::trace_vector voltage; + sample_result voltage; // Now attach the sampler at probeset_id, with sampling schedule sched, writing to voltage sim.add_sampler(arb::one_probe(probeset_id), sched, arb::make_simple_sampler(voltage)); @@ -144,25 +148,22 @@ int main() { sim.run(100*arb::units::ms, 0.01*arb::units::ms); // Write the samples to a json file. - write_trace_json(voltage.at(0)); + write_trace_json(voltage); } -void write_trace_json(const arb::trace_data& trace) { +void write_trace_json(const sample_result& result) { std::string path = "./voltages.json"; nlohmann::json json; json["name"] = "event_gen_demo"; json["units"] = "mV"; - json["cell"] = "0.0"; - json["probe"] = "0"; - - auto& jt = json["data"]["time"]; - auto& jy = json["data"]["voltage"]; - - for (const auto& sample: trace) { - jt.push_back(sample.t); - jy.push_back(sample.v); - } + json["cell"] = "0"; + json["probe"] = "Um"; + std::stringstream loc; + loc << result.metadata.at(0); + json["location"] = loc.str(); + json["data"]["time"] = result.time; + json["data"]["voltage"] = result.values.at(0); std::ofstream file(path); file << std::setw(1) << json << "\n"; diff --git a/example/lfp/CMakeLists.txt b/example/lfp/CMakeLists.txt index a6dbabfd6e..7f3945d799 100644 --- a/example/lfp/CMakeLists.txt +++ b/example/lfp/CMakeLists.txt @@ -1,4 +1,4 @@ add_executable(lfp EXCLUDE_FROM_ALL lfp.cpp) add_dependencies(examples lfp) -target_link_libraries(lfp PRIVATE arbor arborio ext-tinyopt) +target_link_libraries(lfp PRIVATE arbor arborio ext-tinyopt fmt::fmt-header-only) file(COPY plot-lfp.py DESTINATION "${CMAKE_RUNTIME_OUTPUT_DIRECTORY}") diff --git a/example/lfp/lfp.cpp b/example/lfp/lfp.cpp index e32a1a3b15..b033d8be4f 100644 --- a/example/lfp/lfp.cpp +++ b/example/lfp/lfp.cpp @@ -3,6 +3,8 @@ #include #include +#include + #include #include #include @@ -117,7 +119,7 @@ struct lfp_sampler { p.x -= e.x; p.y -= e.y; p.z -= e.z; - double r = std::sqrt(p.x*p.x+p.y*p.y+p.z*p.z); // [μm] + double r = std::sqrt(p.x*p.x + p.y*p.y + p.z*p.z); // [μm] return coef/r; // [MΩ] }); } @@ -125,18 +127,16 @@ struct lfp_sampler { // On receipt of a sequence of cell-wide current samples, apply response matrix and save results to lfp_voltage. arb::sampler_function callback() { - return [this](arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { + return [this](arb::probe_metadata pm, const arb::sample_records& samples) { std::vector currents; + auto data = std::any_cast(samples.values); lfp_voltage.resize(response.size()); - - for (std::size_t i = 0; i(samples[i].data); - assert(data_ptr); - - for (unsigned j = 0; jfirst, data_ptr->second, response[j].begin(), 0.)); + for (std::size_t ix = 0; ix < samples.n_sample; ++ix) { + lfp_time.push_back(samples.time[ix]); + for (std::size_t j = 0; j < response.size(); ++j) { + lfp_voltage[j].push_back(std::inner_product(data + ix*samples.width, data + (ix + 1)*samples.width, + response[j].begin(), + 0.)); } } }; @@ -149,35 +149,7 @@ struct lfp_sampler { std::vector> response; // [MΩ] }; -// JSON output helpers: - -template -struct as_json_array_wrap { - const T& data; - F fn; - as_json_array_wrap(const T& data, const F& fn): data(data), fn(fn) {} - - friend std::ostream& operator<<(std::ostream& out, const as_json_array_wrap& a) { - out << '['; - bool first = true; - for (auto& x: a.data) out << (!first? ", ": (first=false, "")) << a.fn(x); - return out << ']'; - } -}; - -struct { - template - auto operator()(const F& fn) const { - return [&fn](const auto& data) { return as_json_array_wrap(data, fn); }; - } - - auto operator()() const { - return this->operator()([](const auto& x) { return x; }); - } -} as_json_array; - -// Run simulation. - +// Run simulation.1 int main(int argc, char** argv) { // Configuration const double t_stop = 100; // [ms] @@ -197,41 +169,36 @@ int main(int argc, char** argv) { arb::morphology cell_morphology = any_cast(recipe.get_cell_description(0)).morphology(); arb::place_pwlin placed_cell(cell_morphology); - auto probe0_metadata = sim.get_probe_metadata({0, "Itotal"}); - assert(probe0_metadata.size()==1); // Should only be one probe associated with this id. - arb::mcable_list current_cables = *any_cast(probe0_metadata.at(0).meta); + auto probe_metadata = sim.get_probe_metadata({0, "Itotal"}); + assert(probe_metadata.size() == 1); // Should only be one probe associated with this id. + auto probe0_metadata = probe_metadata.at(0); + auto ptr = any_cast(probe0_metadata.meta); + arb::mcable_list current_cables(ptr, ptr + probe0_metadata.width); lfp_sampler lfp(placed_cell, current_cables, electrodes, 3.0); auto sample_schedule = arb::regular_schedule(sample_dt*U::ms); sim.add_sampler(arb::one_probe({0, "Itotal"}), sample_schedule, lfp.callback()); - arb::trace_vector membrane_voltage; - sim.add_sampler(arb::one_probe({0, "Um"}), sample_schedule, make_simple_sampler(membrane_voltage)); + arb::simple_sampler_result membrane_voltage; + sim.add_sampler(arb::one_probe({0, "Um"}), sample_schedule, arb::make_simple_sampler(membrane_voltage)); - arb::trace_vector ionic_current_density; - sim.add_sampler(arb::one_probe({0, "Iion"}), sample_schedule, make_simple_sampler(ionic_current_density)); + arb::simple_sampler_result ionic_current_density; + sim.add_sampler(arb::one_probe({0, "Iion"}), sample_schedule, arb::make_simple_sampler(ionic_current_density)); - arb::trace_vector synapse_g; - sim.add_sampler(arb::one_probe({0, "expsyn-g"}), sample_schedule, make_simple_sampler(synapse_g)); + arb::simple_sampler_result synapse_g; + sim.add_sampler(arb::one_probe({0, "expsyn-g"}), sample_schedule, arb::make_simple_sampler(synapse_g)); sim.run(t_stop*U::ms, dt*U::ms); - // Output results in JSON format suitable for plotting by plot-lfp.py script. - - auto get_t = [](const auto& x) { return x.t; }; - auto get_v = [](const auto& x) { return x.v; }; - auto scale = [](double s) { return [s](const auto& x) { return x*s; }; }; - auto to_xz = [](const auto& p) { return std::array{p.x, p.z}; }; - // Compute synaptic current from synapse conductance and membrane potential. + // TODO(Sampling, TH) Check this... std::vector syn_i; - assert(synapse_g.get(0).size()==membrane_voltage.get(0).size()); - std::transform(synapse_g.get(0).begin(), synapse_g.get(0).end(), membrane_voltage.get(0).begin(), std::back_inserter(syn_i), - [](arb::trace_entry g, arb::trace_entry v) { - assert(g.t==v.t); - return g.v*v.v; - }); + assert(synapse_g.n_sample == membrane_voltage.n_sample); + assert(synapse_g.width == membrane_voltage.width); + std::transform(synapse_g.values.at(0).begin(), synapse_g.values.at(0).end(), + membrane_voltage.values.at(0).begin(), std::back_inserter(syn_i), + [](double g, double v) { return g*v; }); // Collect points from 2-d morphology in vectors of arrays (x, z, radius), one per branch. // (This process will be simplified with improvements to the place_pwlin API.) @@ -243,38 +210,43 @@ int main(int argc, char** argv) { samples.back().push_back(std::array{seg.dist.x, seg.dist.z, seg.dist.radius}); } } - - auto probe_xz = to_xz(placed_cell.at(membrane_voltage.get(0).meta)); + auto mloc = membrane_voltage.metadata.at(0); + auto to_xz = [](const auto& p) { return std::array{p.x, p.z}; }; + auto probe_xz = to_xz(placed_cell.at(mloc)); std::vector> electrodes_xz; - std::transform(electrodes.begin(), electrodes.end(), std::back_inserter(electrodes_xz), to_xz); - - std::cout << - "{\n" - "\"morphology\": {\n" - "\"unit\": \"μm\",\n" - "\"samples\": " << as_json_array(as_json_array(as_json_array()))(samples) << ",\n" - "\"probe\": " << as_json_array()(probe_xz) << ",\n" - "\"electrodes\": " << as_json_array(as_json_array())(electrodes_xz) << "\n" - "},\n" - "\"extracellular potential\": {\n" - "\"unit\": \"μV\",\n" - "\"time\": " << as_json_array()(lfp.lfp_time) << ",\n" - "\"values\": " << as_json_array(as_json_array(scale(1e3)))(lfp.lfp_voltage) << "\n" - "},\n" - "\"synaptic current\": {\n" - "\"unit\": \"nA\",\n" - "\"time\": " << as_json_array(get_t)(synapse_g.get(0)) << ",\n" - "\"value\": " << as_json_array()(syn_i) << "\n" - "},\n" - "\"membrane potential\": {\n" - "\"unit\": \"mV\",\n" - "\"time\": " << as_json_array(get_t)(membrane_voltage.get(0)) << ",\n" - "\"value\": " << as_json_array(get_v)(membrane_voltage.get(0)) << "\n" - "},\n" - "\"ionic current density\": {\n" - "\"unit\": \"A/m²\",\n" - "\"time\": " << as_json_array(get_t)(ionic_current_density.get(0)) << ",\n" - "\"value\": " << as_json_array(get_v)(ionic_current_density.get(0)) << "\n" - "}\n" - "}\n"; + std::transform(electrodes.begin(), electrodes.end(), + std::back_inserter(electrodes_xz), to_xz); + + auto lfp_uV = lfp.lfp_voltage; + std::for_each(lfp_uV.begin(), lfp_uV.end(), + [](auto& xs) { + std::transform(xs.begin(), xs.end(), + xs.begin(), + [](auto x) { return x*1e3; }); + }); + + nlohmann::json json; + json["morphology"]["unit"] = "μm"; + json["morphology"]["samples"] = samples; + json["morphology"]["probe"] = probe_xz; + json["morphology"]["electrodes"] = electrodes_xz; + + json["ionic current density"]["unit"] = "A/m²"; + json["ionic current density"]["time"] = ionic_current_density.time; + json["ionic current density"]["value"] = ionic_current_density.values.at(0); + + json["membrane potential"]["unit"] = "mV"; + json["membrane potential"]["time"] = membrane_voltage.time; + json["membrane potential"]["value"] = membrane_voltage.values.at(0); + + json["extracellular potential"]["unit"] = "μV"; + json["extracellular potential"]["time"] = lfp.lfp_time; + json["extracellular potential"]["values"] = lfp_uV; + + json["synaptic current"]["unit"] = "nA"; + json["synaptic current"]["time"] = synapse_g.time; + json["synaptic current"]["value"] = syn_i; + + std::cout << json << "\n"; + } diff --git a/example/network_description/network_description.cpp b/example/network_description/network_description.cpp index 8f80814be7..5fe3321a1e 100644 --- a/example/network_description/network_description.cpp +++ b/example/network_description/network_description.cpp @@ -5,11 +5,9 @@ #include #include -#include #include #include #include -#include #include #include @@ -60,18 +58,19 @@ ring_params read_options(int argc, char** argv); using arb::cell_gid_type; using arb::cell_kind; using arb::cell_lid_type; -using arb::cell_member_type; using arb::cell_size_type; using arb::time_type; +// result of simple sampler for probe type +using sample_result = arb::simple_sampler_result; + // Writes voltage trace as a json file. -void write_trace_json(const arb::trace_data& trace); +void write_trace_json(const sample_result&); // Generate a cell. arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params); -class ring_recipe: public arb::recipe { -public: +struct ring_recipe: public arb::recipe { ring_recipe(unsigned num_cells, cell_parameters params, unsigned min_delay): num_cells_(num_cells), cell_params_(params), @@ -187,10 +186,10 @@ int main(int argc, char** argv) { std::cout << sup::mask_stream(root); // Print a banner with information about hardware configuration - std::cout << "gpu: " << (has_gpu(context) ? "yes" : "no") << "\n"; - std::cout << "threads: " << num_threads(context) << "\n"; - std::cout << "mpi: " << (has_mpi(context) ? "yes" : "no") << "\n"; - std::cout << "ranks: " << num_ranks(context) << "\n" << std::endl; + std::cout << "gpu: " << (has_gpu(context) ? "yes" : "no") << "\n" + << "threads: " << num_threads(context) << "\n" + << "mpi: " << (has_mpi(context) ? "yes" : "no") << "\n" + << "ranks: " << num_ranks(context) << "\n" << std::endl; auto params = read_options(argc, argv); @@ -211,7 +210,7 @@ int main(int argc, char** argv) { // The schedule for sampling is 10 samples every 1 ms. auto sched = arb::regular_schedule(1*arb::units::ms); // This is where the voltage samples will be stored as (time, value) pairs - arb::trace_vector voltage; + sample_result voltage; // Now attach the sampler at probeset_id, with sampling schedule sched, writing to voltage sim.add_sampler(arb::one_probe(probeset_id), sched, arb::make_simple_sampler(voltage)); @@ -268,7 +267,7 @@ int main(int argc, char** argv) { } // Write the samples to a json file. - if (root) { write_trace_json(voltage.at(0)); } + if (root) { write_trace_json(voltage); } auto profile = arb::profile::profiler_summary(); std::cout << profile << "\n"; @@ -284,27 +283,6 @@ int main(int argc, char** argv) { return 0; } -void write_trace_json(const arb::trace_data& trace) { - std::string path = "./voltages.json"; - - nlohmann::json json; - json["name"] = "ring demo"; - json["units"] = "mV"; - json["cell"] = "0.0"; - json["probe"] = "0"; - - auto& jt = json["data"]["time"]; - auto& jy = json["data"]["voltage"]; - - for (const auto& sample: trace) { - jt.push_back(sample.t); - jy.push_back(sample.v); - } - - std::ofstream file(path); - file << std::setw(1) << json << "\n"; -} - ring_params read_options(int argc, char** argv) { using sup::param_from_json; @@ -341,3 +319,21 @@ ring_params read_options(int argc, char** argv) { return params; } + +void write_trace_json(const sample_result& result) { + std::string path = "./voltages.json"; + + nlohmann::json json; + json["name"] = "network demo"; + json["units"] = "mV"; + json["cell"] = "0"; + json["probe"] = "Um"; + std::stringstream loc; + loc << result.metadata.at(0); + json["location"] = loc.str(); + json["data"]["time"] = result.time; + json["data"]["voltage"] = result.values.at(0); + + std::ofstream file(path); + file << std::setw(1) << json << "\n"; +} diff --git a/example/ornstein_uhlenbeck/ou.cpp b/example/ornstein_uhlenbeck/ou.cpp index 0496ad1396..7fbfe351a0 100644 --- a/example/ornstein_uhlenbeck/ou.cpp +++ b/example/ornstein_uhlenbeck/ou.cpp @@ -13,8 +13,7 @@ #include // a single-cell recipe with probes -class recipe: public arb::recipe { -public: +struct recipe: public arb::recipe { recipe(unsigned ncvs) { using namespace arb; using namespace arborio::literals; @@ -61,17 +60,12 @@ struct sampler { data_.resize(n_cvs*n_steps); } - void operator()(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { - const auto* m = arb::util::any_cast(pm.meta); - arb_assert(n_cvs_ == m->size()); - arb_assert(n_steps_ == n); - - for (std::size_t i=0; i(samples[i].data); - auto [lo, hi] = *data; - arb_assert(static_cast(hi-lo) == n_cvs_); - for (std::size_t j=0; j(pm.meta, samples); + auto width = reader.n_column(); + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + for (std::size_t iy = 0; iy < reader.n_column(); ++iy) { + data_[ix*width + iy] = reader.value(ix, iy); } } } diff --git a/example/plasticity/CMakeLists.txt b/example/plasticity/CMakeLists.txt index 904ea0943a..b64c52ce8a 100644 --- a/example/plasticity/CMakeLists.txt +++ b/example/plasticity/CMakeLists.txt @@ -1,4 +1,4 @@ add_executable(plasticity EXCLUDE_FROM_ALL plasticity.cpp) add_dependencies(examples plasticity) -target_link_libraries(plasticity PRIVATE arbor arborio arborenv) +target_link_libraries(plasticity PRIVATE arbor arborio arborenv fmt::fmt-header-only) diff --git a/example/plasticity/plasticity.cpp b/example/plasticity/plasticity.cpp index 214b640463..dd58bbd95e 100644 --- a/example/plasticity/plasticity.cpp +++ b/example/plasticity/plasticity.cpp @@ -1,6 +1,5 @@ #include #include -#include #include #include @@ -11,6 +10,8 @@ #include #include +#include + using namespace arborio::literals; namespace U = arb::units; @@ -77,14 +78,15 @@ struct recipe: public arb::recipe { // NEVER do this in HPC!!! std::mutex mtx; -void sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { - auto* loc = arb::util::any_cast(pm.meta); - - for (std::size_t i = 0; i lock{mtx}; - auto* value = arb::util::any_cast(samples[i].data); - std::cout << std::fixed << std::setprecision(4) - << "| " << samples[i].time << " | " << loc->pos << " | " << *value << " |\n"; +void sampler(arb::probe_metadata pm, const arb::sample_records& samples) { + using probe_t = arb::cable_probe_membrane_voltage; + auto reader = arb::sample_reader(pm.meta, samples); + std::lock_guard lock{mtx}; + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + auto time = reader.time(ix); + auto value = reader.value(ix); + arb::mlocation meta = reader.metadata(0); + std::cout << fmt::format("| {:7.4f} | {:11.1f} | {:9.2f} |\n", time, meta.pos, value); } } diff --git a/example/probe-demo/CMakeLists.txt b/example/probe-demo/CMakeLists.txt index e7cf7ad751..734b6218bd 100644 --- a/example/probe-demo/CMakeLists.txt +++ b/example/probe-demo/CMakeLists.txt @@ -1,3 +1,3 @@ add_executable(probe-demo EXCLUDE_FROM_ALL probe-demo.cpp) add_dependencies(examples probe-demo) -target_link_libraries(probe-demo PRIVATE arbor ext-tinyopt) +target_link_libraries(probe-demo PRIVATE arbor ext-tinyopt fmt::fmt-header-only) diff --git a/example/probe-demo/probe-demo.cpp b/example/probe-demo/probe-demo.cpp index 2a108ab28d..00d6efd2b0 100644 --- a/example/probe-demo/probe-demo.cpp +++ b/example/probe-demo/probe-demo.cpp @@ -1,12 +1,14 @@ #include #include -#include #include #include #include -#include #include +#include +#include + +#include #include #include #include @@ -19,15 +21,14 @@ // Simulate a cell modelled as a simple cable with HH dynamics, // emitting the results of a user specified probe over time. - -using std::any; -using arb::util::any_cast; +namespace U = arb::units; +using namespace arb::units::literals; const char* help_msg = "[OPTION]... PROBE\n" "\n" - " --dt=TIME set simulation dt to TIME [ms]\n" - " --until=TIME simulate until TIME [ms]\n" + " -d, --dt=TIME set simulation dt to TIME [ms]\n" + " -T, --until=TIME simulate until TIME [ms]\n" " -n, --n-cv=N discretize with N CVs\n" " -t, --sample=TIME take a sample every TIME [ms]\n" " -x, --at=X take sample at relative position X along cable or index of synapse\n" @@ -69,45 +70,73 @@ const char* help_msg = " all_hh_n HH state variable n in each CV\n" " all_expsyn_g expsyn state variable g for all synapses\n"; +enum probe_kind { invalid, state, point, cell }; + struct options { double sim_end = 100.0; // [ms] double sim_dt = 0.025; // [ms] double sample_dt = 1.0; // [ms] unsigned n_cv = 10; - - bool scalar_probe = true; - any probe_addr; + std::any probe_addr; std::string value_name; + probe_kind kind = probe_kind::invalid; }; bool parse_options(options&, int& argc, char** argv); -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*); + +template +std::string show_location(const M& where) { + if constexpr (std::is_same_v, arb::mcable>) { + return fmt::format("(cable {:1d} {:3.1f} {:3.1f})", where.branch, where.prox_pos, where.dist_pos); + } + else if constexpr (std::is_same_v, arb::mlocation>) { + return fmt::format("(location {:1d} {:3.1f})", where.branch, where.pos); + } + else if constexpr (std::is_same_v, arb::cable_probe_point_info>) { + return where.target; + } + else { + throw std::runtime_error{"Unexpected metadata type"}; + } +} + +// Do this once +static std::atomic printed_header = 0; + +template +void sampler(arb::probe_metadata pm, const arb::sample_records& samples) { + auto reader = arb::sample_reader(pm.meta, samples); + // Print CSV header for sample output + if (0 == printed_header.fetch_add(1)) { + std::cout << fmt::format("t", ""); + for (std::size_t iy = 0; iy < reader.n_column(); ++iy) std::cout << ", " << show_location(reader.metadata(iy)); + std::cout << '\n'; + } + + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + std::cout << fmt::format("{:7.3f}", reader.time(ix)); + for (unsigned iy = 0; iy < reader.n_column(); ++iy) { + std::cout << fmt::format(", {:-8.4f}", reader.value(ix, iy)); + } + std::cout << '\n'; + } +} + struct cable_recipe: public arb::recipe { arb::cable_cell_global_properties gprop; - any probe_addr; + std::any probe_addr; - explicit cable_recipe(any probe_addr, unsigned n_cv): - probe_addr(std::move(probe_addr)) - { + explicit cable_recipe(std::any probe_addr, unsigned n_cv): + probe_addr(std::move(probe_addr)) { gprop.default_parameters = arb::neuron_parameter_defaults; gprop.default_parameters.discretization = arb::cv_policy_fixed_per_branch(n_cv); } arb::cell_size_type num_cells() const override { return 1; } - - std::vector get_probes(arb::cell_gid_type) const override { - return {{probe_addr, "probe"}}; - } - - arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { - return arb::cell_kind::cable; - } - - any get_global_properties(arb::cell_kind) const override { - return gprop; - } + std::vector get_probes(arb::cell_gid_type) const override { return {{probe_addr, "probe"}}; } + arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { return arb::cell_kind::cable; } + std::any get_global_properties(arb::cell_kind) const override { return gprop; } arb::util::unique_any get_cell_description(arb::cell_gid_type) const override { const double length = 1000; // [µm] @@ -118,38 +147,51 @@ struct cable_recipe: public arb::recipe { auto decor = arb::decor{} .paint(arb::reg::all(), arb::density("hh")) // HH mechanism over whole cell. - .place(arb::mlocation{0, 0.}, arb::i_clamp{1.*arb::units::nA}, "iclamp") // Inject a 1 nA current indefinitely. - .place(arb::mlocation{0, 0.}, arb::synapse("expsyn"), "synapse1") // a synapse + .place(arb::mlocation{0, 0.0}, arb::i_clamp{1._nA}, "iclamp") // Inject a 1 nA current indefinitely. + .place(arb::mlocation{0, 0.0}, arb::synapse("expsyn"), "synapse1") // a synapse .place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "synapse2"); // another synapse return arb::cable_cell(tree, decor); } virtual std::vector event_generators(arb::cell_gid_type) const override { - return {arb::poisson_generator({"synapse1"}, .005, 0.*arb::units::ms, 0.1*arb::units::kHz), - arb::poisson_generator({"synapse2"}, .1, 0.*arb::units::ms, 0.1*arb::units::kHz)}; + return {arb::poisson_generator({"synapse1"}, .005, 0._ms, 0.1_kHz), + arb::poisson_generator({"synapse2"}, .1, 0._ms, 0.1_kHz)}; } - }; int main(int argc, char** argv) { try { options opt; - if (!parse_options(opt, argc, argv)) { - return 0; - } + if (!parse_options(opt, argc, argv)) return -1; cable_recipe R(opt.probe_addr, opt.n_cv); arb::simulation sim(R); - sim.add_sampler(arb::all_probes, - arb::regular_schedule(opt.sample_dt*arb::units::ms), - opt.scalar_probe? scalar_sampler: vector_sampler); - - // CSV header for sample output: - std::cout << "t, " << (opt.scalar_probe? "x, ": "x0, x1, ") << opt.value_name << '\n'; - - sim.run(opt.sim_end*arb::units::ms, opt.sim_dt*arb::units::ms); + switch (opt.kind) { + case probe_kind::cell: + sim.add_sampler(arb::all_probes, + arb::regular_schedule(opt.sample_dt*U::ms), + sampler); + break; + case probe_kind::state: + sim.add_sampler(arb::all_probes, + arb::regular_schedule(opt.sample_dt*U::ms), + sampler); + break; + case probe_kind::point: + sim.add_sampler(arb::all_probes, + arb::regular_schedule(opt.sample_dt*U::ms), + sampler); + break; + default: + std::cerr << "Invalid probe kind\n"; + return -1; + } + std::cout << "Samples for '" << opt.value_name << "'\n" + << std::string(opt.value_name.size() + 14, '=') + << "\n\n"; + sim.run(opt.sim_end*U::ms, opt.sim_dt*U::ms); } catch (to::option_error& e) { to::usage_error(argv[0], "[OPTIONS]... PROBE\nTry '--help' for more information.", e.what()); @@ -157,97 +199,68 @@ int main(int argc, char** argv) { } catch (std::exception& e) { std::cerr << "caught exception: " << e.what() << "\n"; - return 2; + return -2; } } -void scalar_sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { - auto* loc = any_cast(pm.meta); - auto* point_info = any_cast(pm.meta); - assert((loc != nullptr) || (point_info != nullptr)); - - loc = loc ? loc : &(point_info->loc); - - std::cout << std::fixed << std::setprecision(4); - for (std::size_t i = 0; i(samples[i].data); - assert(value); - - std::cout << samples[i].time << ", " << loc->pos << ", " << *value << '\n'; +static auto any2loc(std::any a) -> arb::mlocation { + auto pos = 0.5; + try { + pos = std::any_cast(a); + } + catch (const std::bad_any_cast& e) { + std::cerr << "Invalid position specification, using default\n"; } + return arb::mlocation { .branch=0, .pos=pos }; } -void vector_sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { - auto* cables_ptr = any_cast(pm.meta); - auto* point_info_ptr = any_cast*>(pm.meta); - - assert((cables_ptr != nullptr) || (point_info_ptr != nullptr)); - - unsigned n_entities = cables_ptr ? cables_ptr->size() : point_info_ptr->size(); - - std::cout << std::fixed << std::setprecision(4); - for (std::size_t i = 0; i(samples[i].data); - assert(value_range); - const auto& [lo, hi] = *value_range; - assert(n_entities==hi-lo); - - for (unsigned j = 0; j arb::cell_tag_type { + arb::cell_tag_type pos = "synapse1"; + try { + pos = std::any_cast(a); + } + catch (const std::bad_any_cast& e) { + std::cerr << "Invalid synapse specification, using default\n"; } + return pos; } + bool parse_options(options& opt, int& argc, char** argv) { using std::get; using namespace to; auto do_help = [&]() { usage(argv[0], 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. - - using probe_spec_t = std::tuple>; - + // Map probe argument to output variable name and a lambda that makes specific probe address from a location. + using probe_spec_t = std::tuple>; std::pair probe_tbl[] { // located probes - {"v", {"v", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_membrane_voltage{L{0, x}}; }}}, - {"i_axial", {"i_axial", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_axial_current{L{0, x}}; }}}, - {"j_ion", {"j_ion", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_total_ion_current_density{L{0, x}}; }}}, - {"j_na", {"j_na", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_ion_current_density{L{0, x}, "na"}; }}}, - {"j_k", {"j_k", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_ion_current_density{L{0, x}, "k"}; }}}, - {"c_na", {"c_na", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_ion_int_concentration{L{0, x}, "na"}; }}}, - {"c_k", {"c_k", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_ion_int_concentration{L{0, x}, "k"}; }}}, - {"hh_m", {"hh_m", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_density_state{L{0, x}, "hh", "m"}; }}}, - {"hh_h", {"hh_h", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_density_state{L{0, x}, "hh", "h"}; }}}, - {"hh_n", {"hh_n", true, [](std::any a) { auto x = std::any_cast(a); return arb::cable_probe_density_state{L{0, x}, "hh", "n"}; }}}, - {"expsyn_g", {"expsyn_ g", true, [](std::any a) { auto t = std::any_cast(a); return arb::cable_probe_point_state{t, "expsyn", "g"}; }}}, + {"v", {"v", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_membrane_voltage{any2loc(a)}; }}}, + {"i_axial", {"i_axial", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_axial_current{any2loc(a)}; }}}, + {"j_ion", {"j_ion", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_total_ion_current_density{any2loc(a)}; }}}, + {"j_na", {"j_na", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_ion_current_density{any2loc(a), "na"}; }}}, + {"j_k", {"j_k", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_ion_current_density{any2loc(a), "k"}; }}}, + {"c_na", {"c_na", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_ion_int_concentration{any2loc(a), "na"}; }}}, + {"c_k", {"c_k", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_ion_int_concentration{any2loc(a), "k"}; }}}, + {"hh_m", {"hh_m", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_density_state{any2loc(a), "hh", "m"}; }}}, + {"hh_h", {"hh_h", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_density_state{any2loc(a), "hh", "h"}; }}}, + {"hh_n", {"hh_n", probe_kind::state, [](std::any a) -> std::any { return arb::cable_probe_density_state{any2loc(a), "hh", "n"}; }}}, + {"expsyn_g", {"expsyn_g", probe_kind::point, [](std::any a) -> std::any { return arb::cable_probe_point_state{any2syn(a), "expsyn", "g"}; }}}, // all-of-cell probes - {"all_v", {"v", false, [](std::any) { return arb::cable_probe_membrane_voltage_cell{}; }}}, - {"all_i_ion", {"i_ion", false, [](std::any) { return arb::cable_probe_total_ion_current_cell{}; }}}, - {"all_i_na", {"i_na", false, [](std::any) { return arb::cable_probe_ion_current_cell{"na"}; }}}, - {"all_i_k", {"i_k", false, [](std::any) { return arb::cable_probe_ion_current_cell{"k"}; }}}, - {"all_i", {"i", false, [](std::any) { return arb::cable_probe_total_current_cell{}; }}}, - {"all_c_na", {"c_na", false, [](std::any) { return arb::cable_probe_ion_int_concentration_cell{"na"}; }}}, - {"all_c_k", {"c_k", false, [](std::any) { return arb::cable_probe_ion_int_concentration_cell{"k"}; }}}, - {"all_hh_m", {"hh_m", false, [](std::any) { return arb::cable_probe_density_state_cell{"hh", "m"}; }}}, - {"all_hh_h", {"hh_h", false, [](std::any) { return arb::cable_probe_density_state_cell{"hh", "h"}; }}}, - {"all_hh_n", {"hh_n", false, [](std::any) { return arb::cable_probe_density_state_cell{"hh", "n"}; }}}, - {"all_expsyn_g", {"expsyn_ g", false, [](std::any) { return arb::cable_probe_point_state_cell{"expsyn", "g"}; }}}, + {"all_v", {"v", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_membrane_voltage_cell{}; }}}, + {"all_i_ion", {"i_ion", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_total_ion_current_cell{}; }}}, + {"all_i_na", {"i_na", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_ion_current_cell{"na"}; }}}, + {"all_i_k", {"i_k", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_ion_current_cell{"k"}; }}}, + {"all_i", {"i", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_total_current_cell{}; }}}, + {"all_c_na", {"c_na", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_ion_int_concentration_cell{"na"}; }}}, + {"all_c_k", {"c_k", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_ion_int_concentration_cell{"k"}; }}}, + {"all_hh_m", {"hh_m", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_density_state_cell{"hh", "m"}; }}}, + {"all_hh_h", {"hh_h", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_density_state_cell{"hh", "h"}; }}}, + {"all_hh_n", {"hh_n", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_density_state_cell{"hh", "n"}; }}}, + {"all_expsyn_g", {"expsyn_g", probe_kind::cell, [](std::any) -> std::any { return arb::cable_probe_point_state_cell{"expsyn", "g"}; }}}, }; - probe_spec_t probe_spec; - std::any p_pos; - auto double_or_string = [](const char* arg) -> to::maybe { try { return {{std::stod(arg)}}; @@ -257,30 +270,26 @@ bool parse_options(options& opt, int& argc, char** argv) { } }; + probe_spec_t probe_spec; + std::any p_pos; to::option cli_opts[] = { { to::action(do_help), to::flag, to::exit, "-h", "--help" }, - { opt.sim_dt, "--dt" }, - { opt.sim_end, "--until" }, + { opt.sim_dt, "-d", "--dt" }, + { opt.sim_end, "-T", "--until" }, { opt.sample_dt, "-t", "--sample" }, { to::sink(p_pos, double_or_string), "-x", "--at" }, { opt.n_cv, "-n", "--n-cv" }, { {probe_spec, to::keywords(probe_tbl)}, to::single }, }; - const auto& [p_name, p_scalar, p_addr] = probe_spec; - if (!p_pos.has_value() && (p_name == "exp_syn_g")) { - p_pos = "synapse0"; - } - else { - p_pos = 0.5; - } + const auto& [p_name, p_kind, p_addr] = probe_spec; if (!to::run(cli_opts, argc, argv+1)) return false; if (!p_addr) throw to::user_option_error("missing PROBE"); if (argv[1]) throw to::user_option_error("unrecognized option"); - opt.value_name = p_name; - opt.scalar_probe = p_scalar; - opt.probe_addr = p_addr(p_pos); + opt.value_name = p_name; + opt.probe_addr = p_addr(p_pos); + opt.kind = p_kind; return true; } diff --git a/example/remote/remote.cpp b/example/remote/remote.cpp index a966d72a90..9b476c8327 100644 --- a/example/remote/remote.cpp +++ b/example/remote/remote.cpp @@ -62,7 +62,7 @@ struct mpi_handle { mpi_handle setup_mpi(); -void sampler(arb::probe_metadata, std::size_t, const arb::sample_record*); +void sampler(arb::probe_metadata, const arb::sample_records&); std::vector> trace; @@ -147,13 +147,13 @@ mpi_handle setup_mpi() { return result; } -void sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { - for (std::size_t ix = 0; ix < n; ++ix) { - if (pm.id.gid != 0) continue; - if (pm.id.tag != "Um") continue; - const auto& sample = samples[ix]; - auto value = *arb::util::any_cast(sample.data); - auto time = sample.time; +void sampler(arb::probe_metadata pm, const arb::sample_records& samples) { + if (pm.id.gid != 0) return; + if (pm.id.tag != "Um") return; + auto reader = arb::sample_reader(pm.meta, samples); + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + double time = reader.time(ix); + double value = reader.value(ix); trace.emplace_back(time, value); } } diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp index db488dfcd4..99b23a81b1 100644 --- a/example/ring/ring.cpp +++ b/example/ring/ring.cpp @@ -52,62 +52,46 @@ using arb::cell_size_type; using arb::cell_kind; using arb::time_type; +// result of simple sampler for probe type +using sample_result = arb::simple_sampler_result; + // Writes voltage trace as a json file. -void write_trace_json(const arb::trace_data& trace); +void write_trace_json(const sample_result&); // Generate a cell. arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params); -class ring_recipe: public arb::recipe { -public: +struct ring_recipe: public arb::recipe { ring_recipe(unsigned num_cells, cell_parameters params, unsigned min_delay): num_cells_(num_cells), cell_params_(params), - min_delay_(min_delay) - { + min_delay_(min_delay) { gprop_.default_parameters = arb::neuron_parameter_defaults; } - cell_size_type num_cells() const override { - return num_cells_; - } + cell_size_type num_cells() const override { return num_cells_; } - arb::util::unique_any get_cell_description(cell_gid_type gid) const override { - return branch_cell(gid, cell_params_); - } + arb::util::unique_any get_cell_description(cell_gid_type gid) const override { return branch_cell(gid, cell_params_); } - cell_kind get_cell_kind(cell_gid_type gid) const override { - return cell_kind::cable; - } + cell_kind get_cell_kind(cell_gid_type gid) const override { return cell_kind::cable; } // Each cell has one incoming connection, from cell with gid-1. std::vector connections_on(cell_gid_type gid) const override { - std::vector cons; - cell_gid_type src = gid? gid-1: num_cells_-1; - cons.push_back(arb::cell_connection({src, "detector"}, {"primary_syn"}, event_weight_, min_delay_*U::ms)); - return cons; + cell_gid_type src = gid ? gid - 1: num_cells_ - 1; + return { arb::cell_connection({src, "detector"}, {"primary_syn"}, event_weight_, min_delay_*U::ms) }; } - // Return one event generator on gid 0. This generates a single event that will - // kick start the spiking. + // Return one event generator on gid 0. This generates a single event that + // will kick start the spiking. std::vector event_generators(cell_gid_type gid) const override { - std::vector gens; - if (!gid) { - gens.push_back(arb::explicit_generator_from_milliseconds({"primary_syn"}, event_weight_, std::vector{1.0})); - } - return gens; - } - - std::vector get_probes(cell_gid_type gid) const override { - // Measure membrane voltage at end of soma. - arb::mlocation loc{0, 0.0}; - return {{arb::cable_probe_membrane_voltage{loc}, "Um"}}; + if (!gid) return { arb::explicit_generator_from_milliseconds({"primary_syn"}, event_weight_, std::vector{1.0}) }; + return {}; } - std::any get_global_properties(arb::cell_kind) const override { - return gprop_; - } + // Measure membrane voltage at end of soma. + std::vector get_probes(cell_gid_type gid) const override { return {{arb::cable_probe_membrane_voltage{arb::mlocation {0, 0.0}}, "Um"}}; } + std::any get_global_properties(arb::cell_kind) const override { return gprop_; } private: cell_size_type num_cells_; @@ -159,13 +143,12 @@ int main(int argc, char** argv) { arb::simulation sim(recipe, context, decomposition); // Set up the probe that will measure voltage in the cell. - // The id of the only probe on the cell: the cell_member type points to (cell 0, probe 0) auto probeset_id = arb::cell_address_type{0, "Um"}; // The schedule for sampling every 1 ms. auto sched = arb::regular_schedule(1*arb::units::ms); // This is where the voltage samples will be stored as (time, value) pairs - arb::trace_vector voltage; + sample_result voltage; // Now attach the sampler at probeset_id, with sampling schedule sched, writing to voltage sim.add_sampler(arb::one_probe(probeset_id), sched, arb::make_simple_sampler(voltage)); @@ -211,9 +194,7 @@ int main(int argc, char** argv) { } // Write the samples to a json file. - if (root) { - write_trace_json(voltage.at(0)); - } + if (root) write_trace_json(voltage); auto profile = arb::profile::profiler_summary(); std::cout << profile << "\n"; @@ -229,22 +210,19 @@ int main(int argc, char** argv) { return 0; } -void write_trace_json(const arb::trace_data& trace) { +void write_trace_json(const sample_result& result) { std::string path = "./voltages.json"; nlohmann::json json; - json["name"] = "ring demo"; + json["name"] = "ring_demo"; json["units"] = "mV"; - json["cell"] = "0.0"; - json["probe"] = "0"; - - auto& jt = json["data"]["time"]; - auto& jy = json["data"]["voltage"]; - - for (const auto& sample: trace) { - jt.push_back(sample.t); - jy.push_back(sample.v); - } + json["cell"] = "0"; + json["probe"] = "Um"; + std::stringstream loc; + loc << result.metadata.at(0); + json["location"] = loc.str(); + json["data"]["time"] = result.time; + json["data"]["voltage"] = result.values.at(0); std::ofstream file(path); file << std::setw(1) << json << "\n"; diff --git a/example/single/single.cpp b/example/single/single.cpp index 53b05855ef..a24733ffb8 100644 --- a/example/single/single.cpp +++ b/example/single/single.cpp @@ -1,11 +1,10 @@ #include -#include #include -#include #include #include #include +#include #include #include @@ -14,12 +13,14 @@ #include #include -#include - +#include #include using namespace arborio::literals; +// result of simple sampler for probe type +using sample_result = arb::simple_sampler_result; + struct options { std::string swc_file; double t_end = 20; @@ -95,20 +96,32 @@ int main(int argc, char** argv) { arb::simulation sim(R); // Attach a sampler to the probe described in the recipe, sampling every 0.1 ms. - arb::trace_vector traces; + sample_result traces; sim.add_sampler(arb::all_probes, arb::regular_schedule(0.1*arb::units::ms), arb::make_simple_sampler(traces)); sim.run(opt.t_end*arb::units::ms, opt.dt*arb::units::ms); - for (auto entry: traces.at(0)) { - std::cout << entry.t << ", " << entry.v << "\n"; + std::cout << "Membrane Potential\n" + "==================\n\n" + " t"; + for (std::size_t iy = 0; iy < traces.width; ++iy) { + arb::mlocation loc = traces.metadata.at(iy); + std::cout << fmt::format(", Um@(location {} {})\n", loc.branch, loc.pos); + } + + for (std::size_t ix = 0; ix < traces.n_sample; ++ix) {{} + std::cout << fmt::format("{:7.3f}", traces.time.at(ix)); + for (std::size_t iy = 0; iy < traces.width; ++iy) { + std::cout << fmt::format(", {:-8.4f}", traces.values[iy][ix]); + } + std::cout << '\n'; } } catch (std::exception& e) { std::cerr << "caught exception: " << e.what() << "\n"; - return 2; + return -2; } } diff --git a/example/v_clamp/v-clamp.cpp b/example/v_clamp/v-clamp.cpp index 86f4441c3f..e3147af275 100644 --- a/example/v_clamp/v-clamp.cpp +++ b/example/v_clamp/v-clamp.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include @@ -12,8 +13,7 @@ #include #include -#include - +#include #include using namespace arborio::literals; @@ -27,6 +27,9 @@ struct options { arb::cv_policy policy = arb::default_cv_policy(); }; +// result of simple sampler for probe type +using sample_result = arb::simple_sampler_result; + options parse_options(int argc, char** argv); arborio::loaded_morphology default_morphology(); @@ -96,14 +99,25 @@ int main(int argc, char** argv) { arb::simulation sim(R); // Attach a sampler to the probe described in the recipe, sampling every 0.1 ms. - - arb::trace_vector traces; + sample_result traces; sim.add_sampler(arb::all_probes, arb::regular_schedule(0.1*arb::units::ms), arb::make_simple_sampler(traces)); sim.run(opt.t_end*arb::units::ms, opt.dt*arb::units::ms); - for (auto entry: traces.at(0)) { - std::cout << entry.t << ", " << entry.v << "\n"; + std::cout << "Membrane Potential\n" + "==================\n\n" + " t"; + for (std::size_t iy = 0; iy < traces.width; ++iy) { + arb::mlocation loc = traces.metadata.at(iy); + std::cout << fmt::format(", Um@(location {} {})\n", loc.branch, loc.pos); + } + + for (std::size_t ix = 0; ix < traces.n_sample; ++ix) {{} + std::cout << fmt::format("{:7.3f}", traces.time.at(ix)); + for (std::size_t iy = 0; iy < traces.width; ++iy) { + std::cout << fmt::format(", {:-8.4f}", traces.values[iy][ix]); + } + std::cout << '\n'; } } catch (std::exception& e) { diff --git a/pyproject.toml b/pyproject.toml index 2551d531d1..c492a25fe9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,6 +57,7 @@ exclude = [ "build", "dist", "ext", + "install-*", "doc/scripts/inputs.py", "doc/scripts/make_images.py", ".*", diff --git a/python/example/connectivity/brunel.py b/python/example/connectivity/brunel.py index 1426880994..784a61598f 100644 --- a/python/example/connectivity/brunel.py +++ b/python/example/connectivity/brunel.py @@ -36,7 +36,7 @@ def network_description(self): (random {self.seed} {self.p}))""" inh = f"(gid-range {self.n_exc} {self.N})" weight = f"""(if-else (source-cell {inh}) - (scalar {self.g*self.weight}) + (scalar {self.g * self.weight}) (scalar {self.weight}))""" delay = "(scalar 0.5)" return A.network_description(rand, weight, delay, {}) diff --git a/python/example/diffusion.py b/python/example/diffusion.py index b4281815ae..8c3f61eeb9 100644 --- a/python/example/diffusion.py +++ b/python/example/diffusion.py @@ -59,14 +59,15 @@ def event_generators(self, _): cel = A.cable_cell(tree, dec, discretization=A.cv_policy("(max-extent 5)")) rec = recipe(cel, prb) sim = A.simulation(rec) -hdl = (sim.sample((0, "nad"), A.regular_schedule(0.1 * U.ms)),) +hdl = (sim.sample((0, "nad"), A.regular_schedule(0.01 * U.ms)),) -sim.run(tfinal=0.5 * U.ms) +sim.run(tfinal=0.5 * U.ms, dt=0.01 * U.ms) sns.set_theme() fg, ax = plt.subplots() for h in hdl: for d, m in sim.samples(h): + print(d) # Plot for lbl, ix in zip(m, range(1, d.shape[1])): ax.plot(d[:, 0], d[:, ix], label=lbl) diff --git a/python/probes.cpp b/python/probes.cpp index bb264688d4..424a7fddf4 100644 --- a/python/probes.cpp +++ b/python/probes.cpp @@ -41,57 +41,32 @@ struct recorder_base: sample_recorder { return py::cast(meta_); } - void reset() override { - sample_raw_.clear(); - } + void reset() override { sample_raw_.clear(); } protected: - Meta meta_; + Meta* meta_; std::vector sample_raw_; std::ptrdiff_t stride_; - recorder_base(const Meta* meta_ptr, std::ptrdiff_t width): - meta_(*meta_ptr), stride_(1+width) + recorder_base(const Meta* meta_ptr, std::size_t width): + meta_(const_cast(meta_ptr)), stride_(1 + width) {} }; -template -struct recorder_cable_scalar: recorder_base { - using recorder_base::sample_raw_; - - void record(any_ptr, std::size_t n_sample, const arb::sample_record* records) override { - for (std::size_t i = 0; i(records[i].data)) { - sample_raw_.push_back(records[i].time); - sample_raw_.push_back(*v_ptr); - } - else { - throw arb::arbor_internal_error("unexpected sample type"); - } - } - } - -protected: - recorder_cable_scalar(const Meta* meta_ptr): recorder_base(meta_ptr, 1) {} -}; - struct recorder_lif: recorder_base { using recorder_base::sample_raw_; - void record(any_ptr, std::size_t n_sample, const arb::sample_record* records) override { - for (std::size_t i = 0; i(records[i].data)) { - sample_raw_.push_back(records[i].time); - sample_raw_.push_back(*v_ptr); - } - else { - std::string ty = records[i].data.type().name(); - throw arb::arbor_internal_error("LIF recorder: unexpected sample type " + ty); - } + void record(any_ptr pm, const arb::sample_records& records) override { + auto reader = arb::sample_reader(pm, records); + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + auto t = reader.time(ix); + sample_raw_.push_back(t); + auto v = reader.value(ix); + sample_raw_.push_back(v); } } - recorder_lif(const arb::lif_probe_metadata* meta_ptr): recorder_base(meta_ptr, 1) {} + recorder_lif(const arb::lif_probe_metadata* meta_ptr, std::size_t): recorder_base(meta_ptr, 1) {} }; @@ -99,14 +74,14 @@ template struct recorder_cable_vector: recorder_base { using recorder_base::sample_raw_; - void record(any_ptr, std::size_t n_sample, const arb::sample_record* records) override { - for (std::size_t i = 0; i(records[i].data)) { - sample_raw_.push_back(records[i].time); - sample_raw_.insert(sample_raw_.end(), v_ptr->first, v_ptr->second); - } - else { - throw arb::arbor_internal_error("unexpected sample type"); + void record(any_ptr pm, const arb::sample_records& records) override { + auto reader = arb::sample_reader(pm, records); + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + auto t = reader.time(ix); + sample_raw_.push_back(t); + for (std::size_t iy = 0; iy < reader.n_column(); ++iy) { + auto v = reader.value(ix, iy); + sample_raw_.push_back(v); } } } @@ -117,39 +92,35 @@ struct recorder_cable_vector: recorder_base { }; // Specific recorder classes: - -struct recorder_cable_scalar_mlocation: recorder_cable_scalar { - explicit recorder_cable_scalar_mlocation(const arb::mlocation* meta_ptr): - recorder_cable_scalar(meta_ptr) {} +struct recorder_cable_vector_mcable: recorder_cable_vector { + explicit recorder_cable_vector_mcable(const arb::mcable* ptr, std::size_t n): + recorder_cable_vector(ptr, n) {} }; -struct recorder_cable_scalar_point_info: recorder_cable_scalar { - explicit recorder_cable_scalar_point_info(const arb::cable_probe_point_info* meta_ptr): - recorder_cable_scalar(meta_ptr) {} +struct recorder_cable_vector_mlocation: recorder_cable_vector { + explicit recorder_cable_vector_mlocation(const arb::mlocation* ptr, std::size_t n): + recorder_cable_vector(ptr, n) {} }; -struct recorder_cable_vector_mcable: recorder_cable_vector { - explicit recorder_cable_vector_mcable(const arb::mcable_list* meta_ptr): - recorder_cable_vector(meta_ptr, std::ptrdiff_t(meta_ptr->size())) {} -}; -struct recorder_cable_vector_point_info: recorder_cable_vector> { - explicit recorder_cable_vector_point_info(const std::vector* meta_ptr): - recorder_cable_vector(meta_ptr, std::ptrdiff_t(meta_ptr->size())) {} +struct recorder_cable_vector_point_info: recorder_cable_vector { + explicit recorder_cable_vector_point_info(const arb::cable_probe_point_info* ptr, std::size_t n): + recorder_cable_vector(ptr, n) {} }; // Helper for registering sample recorder factories and (trivial) metadata conversions. - template void register_probe_meta_maps(pyarb_global_ptr g) { g->recorder_factories.assign( - [](any_ptr meta_ptr) -> std::unique_ptr { - return std::unique_ptr(new Recorder(any_cast(meta_ptr))); + [](any_ptr meta_ptr, std::size_t n) -> std::unique_ptr { + return std::make_unique(any_cast(meta_ptr), n); }); g->probe_meta_converters.assign( - [](any_ptr meta_ptr) -> py::object { - return py::cast(*any_cast(meta_ptr)); + [](any_ptr ptr, std::size_t n) -> py::object { + auto raw = any_cast(ptr); + std::vector result(raw, raw + n); + return py::cast(result); }); } @@ -270,7 +241,6 @@ void register_cable_probes(pybind11::module& m, pyarb_global_ptr global_ptr) { return pprintf("", m.target, m.lid, m.multiplicity, m.loc);}); // Probe address constructors: - m.def("lif_probe_voltage", &lif_probe_voltage, "Probe specification for LIF cell membrane voltage.", "tag"_a); @@ -348,11 +318,9 @@ void register_cable_probes(pybind11::module& m, pyarb_global_ptr global_ptr) { "ion"_a, "tag"_a); // Add probe metadata to maps for converters and recorders. - - register_probe_meta_maps(global_ptr); - register_probe_meta_maps(global_ptr); - register_probe_meta_maps(global_ptr); - register_probe_meta_maps, recorder_cable_vector_point_info>(global_ptr); + register_probe_meta_maps(global_ptr); + register_probe_meta_maps(global_ptr); + register_probe_meta_maps(global_ptr); register_probe_meta_maps(global_ptr); } diff --git a/python/pyarb.hpp b/python/pyarb.hpp index 1b4241c113..14907af73c 100644 --- a/python/pyarb.hpp +++ b/python/pyarb.hpp @@ -30,7 +30,7 @@ namespace pyarb { // Sample recorder object interface. struct sample_recorder { - virtual void record(arb::util::any_ptr meta, std::size_t n_sample, const arb::sample_record* records) = 0; + virtual void record(arb::util::any_ptr meta, const arb::sample_records& records) = 0; virtual pybind11::object samples() const = 0; virtual pybind11::object meta() const = 0; virtual void reset() = 0; @@ -40,7 +40,7 @@ struct sample_recorder { // Recorder 'factory' type: given an any_ptr to probe metadata of a specific subset of types, // return a corresponding sample_recorder instance. -using sample_recorder_factory = std::function (arb::util::any_ptr)>; +using sample_recorder_factory = std::function (arb::util::any_ptr, std::size_t)>; // Holds map: probe metadata pointer type → recorder object factory. @@ -52,9 +52,9 @@ struct recorder_factory_map { map_[typeid(const Meta*)] = std::move(rf); } - std::unique_ptr make_recorder(arb::util::any_ptr meta) const { + std::unique_ptr make_recorder(arb::util::any_ptr meta, std::size_t n) const { try { - return map_.at(meta.type())(meta); + return map_.at(meta.type())(meta, n); } catch (std::out_of_range&) { std::string ty = meta.type().name(); @@ -64,24 +64,19 @@ struct recorder_factory_map { }; // Probe metadata to Python object converter. - -using probe_meta_converter = std::function; +using probe_meta_converter = std::function; struct probe_meta_cvt_map { std::unordered_map map_; template - void assign(probe_meta_converter cvt) { - map_[typeid(const Meta*)] = std::move(cvt); - } + void assign(probe_meta_converter cvt) { map_[typeid(const Meta*)] = std::move(cvt); } - pybind11::object convert(arb::util::any_ptr meta) const { + pybind11::object convert(arb::util::any_ptr meta, std::size_t n) const { if (auto iter = map_.find(meta.type()); iter!=map_.end()) { - return iter->second(meta); - } - else { - return pybind11::none(); + return iter->second(meta, n); } + return pybind11::none(); } }; diff --git a/python/simulation.cpp b/python/simulation.cpp index d7ea8a9ad2..6a16b0ec42 100644 --- a/python/simulation.cpp +++ b/python/simulation.cpp @@ -42,9 +42,7 @@ class simulation_shim { struct sampler_callback { std::shared_ptr recorders; - void operator()(arb::probe_metadata pm, std::size_t n_record, const arb::sample_record* records) { - recorders->at(pm.index)->record(pm.meta, n_record, records); - } + void operator()(arb::probe_metadata pm, const arb::sample_records& records) { recorders->at(pm.index)->record(pm.meta, records); } py::list samples() const { std::size_t size = recorders->size(); @@ -61,8 +59,7 @@ class simulation_shim { public: simulation_shim(std::shared_ptr& rec, const context_shim& ctx, const arb::domain_decomposition& decomp, std::uint64_t seed, pyarb_global_ptr global_ptr): - global_ptr_(global_ptr) - { + global_ptr_(global_ptr) { try { sim_.reset(new arb::simulation(recipe_shim(rec), ctx.context, decomp, seed)); } @@ -129,9 +126,9 @@ class simulation_shim { spike_record_.insert(spike_record_.end(), spikes.begin(), spikes.end()); // Sort the newly appended spikes. std::sort(spike_record_.begin()+old_size, spike_record_.end(), - [](const auto& lhs, const auto& rhs) { - return std::tie(lhs.time, lhs.source.gid, lhs.source.index)get_probe_metadata(probeset_id)) { - result.append(global_ptr_->probe_meta_converters.convert(pm.meta)); + result.append(global_ptr_->probe_meta_converters.convert(pm.meta, pm.width)); } return result; } @@ -166,7 +163,7 @@ class simulation_shim { std::shared_ptr recorders{new sample_recorder_vec}; for (const arb::probe_metadata& pm: sim_->get_probe_metadata(probeset_id)) { - recorders->push_back(global_ptr_->recorder_factories.make_recorder(pm.meta)); + recorders->push_back(global_ptr_->recorder_factories.make_recorder(pm.meta, pm.width)); } // Constructed callbacks are passed to the underlying simulator object, _and_ a copy diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp index 6bd0f7b001..14faf7e708 100644 --- a/python/single_cell_model.cpp +++ b/python/single_cell_model.cpp @@ -13,15 +13,12 @@ #include #include #include -#include #include "event_generator.hpp" #include "error.hpp" #include "strprintf.hpp" #include "label_dict.hpp" -using arb::util::any_cast; - namespace pyarb { // @@ -55,18 +52,18 @@ struct trace_callback { trace_callback(std::vector& ts, const std::unordered_map& ls): traces_(ts), locmap_(ls) {} - void operator()(arb::probe_metadata md, std::size_t n, const arb::sample_record* recs) { - // Push each (time, value) pair from the last epoch into trace_. - auto* loc = any_cast(md.meta); - if (locmap_.count(*loc)) { - auto& trace = traces_[locmap_.at(*loc)]; - for (std::size_t i=0; i(recs[i].data)) { - trace.t.push_back(recs[i].time); - trace.v.push_back(*p); - } - else { - throw std::runtime_error("unexpected sample type"); + void operator()(arb::probe_metadata pm, const arb::sample_records& recs) { + auto reader = arb::sample_reader(pm.meta, recs); + + for (std::size_t j = 0; j < reader.n_column(); ++j) { + const auto& loc = reader.metadata(j); + if (locmap_.count(loc)) { + auto& trace = traces_[locmap_.at(loc)]; + for (std::size_t i = 0; i < reader.n_row(); ++i) { + auto time = reader.time(i); + auto value = reader.value(i, j); + trace.t.push_back(time); + trace.v.push_back(value); } } } diff --git a/python/test/unit/test_multiple_connections.py b/python/test/unit/test_multiple_connections.py index b15a5d9236..081a4b2d43 100644 --- a/python/test/unit/test_multiple_connections.py +++ b/python/test/unit/test_multiple_connections.py @@ -55,9 +55,10 @@ def evaluate_outcome(self, sim, handle_mem): # membrane potential should temporarily be above the spiking threshold at around 1.0 ms (only testing this if the current node keeps the data, cf. GitHub issue #1892) if len(sim.samples(handle_mem)) > 0: data_mem, _ = sim.samples(handle_mem)[0] - # print(data_mem[(data_mem[:, 0] >= 1.0), 1]) - self.assertGreater(data_mem[(np.round(data_mem[:, 0], 2) == 1.02), 1], -10) - self.assertLess(data_mem[(np.round(data_mem[:, 0], 2) == 1.05), 1], -10) + um_pre = data_mem[(np.round(data_mem[:, 0], 2) == 1.02), 1][0] + um_pst = data_mem[(np.round(data_mem[:, 0], 2) == 1.05), 1][0] + self.assertGreater(um_pre, -10) + self.assertLess(um_pst, -10) # neuron 3 should spike at around 1.0 ms, when the added input from all connections will cause threshold crossing spike_times = sim.spikes()["time"] @@ -143,7 +144,7 @@ def cell_description(self, gid): sim.record(A.spike_recording.all) # create schedule and handle to record the membrane potential of neuron 3 - reg_sched = A.regular_schedule(0 * U.ms, self.dt * U.ms, self.runtime * U.ms) + reg_sched = A.regular_schedule(self.dt * U.ms) handle_mem = sim.sample((3, "Um"), reg_sched) # run the simulation @@ -372,7 +373,7 @@ def cell_description(self, gid): sim.record(A.spike_recording.all) # create schedule and handle to record the membrane potential of neuron 3 - reg_sched = A.regular_schedule(0 * U.ms, self.dt * U.ms, self.runtime * U.ms) + reg_sched = A.regular_schedule(self.dt * U.ms) handle_mem = sim.sample((3, "Um"), reg_sched) # run the simulation diff --git a/python/test/unit/test_probes.py b/python/test/unit/test_probes.py index d38d4ce3dd..afa405e325 100644 --- a/python/test/unit/test_probes.py +++ b/python/test/unit/test_probes.py @@ -96,7 +96,7 @@ def test_probe_addr_metadata(self): m = sim.probe_metadata((0, "Um")) self.assertEqual(1, len(m)) - self.assertEqual(A.location(0, 0.0), m[0]) + self.assertEqual([A.location(0, 0.0)], m[0]) m = sim.probe_metadata((0, "Um-all")) self.assertEqual(1, len(m)) @@ -104,11 +104,11 @@ def test_probe_addr_metadata(self): m = sim.probe_metadata((0, "Iax")) self.assertEqual(1, len(m)) - self.assertEqual(A.location(0, 0.02), m[0]) + self.assertEqual([A.location(0, 0.02)], m[0]) m = sim.probe_metadata((0, "Iion")) self.assertEqual(1, len(m)) - self.assertEqual(A.location(0, 0.03), m[0]) + self.assertEqual([A.location(0, 0.03)], m[0]) m = sim.probe_metadata((0, "Iion-all")) self.assertEqual(1, len(m)) @@ -120,7 +120,7 @@ def test_probe_addr_metadata(self): m = sim.probe_metadata((0, "hh-m")) self.assertEqual(1, len(m)) - self.assertEqual(A.location(0, 0.06), m[0]) + self.assertEqual([A.location(0, 0.06)], m[0]) m = sim.probe_metadata((0, "hh-n-all")) self.assertEqual(1, len(m)) @@ -128,9 +128,9 @@ def test_probe_addr_metadata(self): m = sim.probe_metadata((0, "expsyn-g")) self.assertEqual(1, len(m)) - self.assertEqual(A.location(0, 0.08), m[0].location) - self.assertEqual(1, m[0].multiplicity) - self.assertEqual("syn0", m[0].target) + self.assertEqual(A.location(0, 0.08), m[0][0].location) + self.assertEqual(1, m[0][0].multiplicity) + self.assertEqual("syn0", m[0][0].target) m = sim.probe_metadata((0, "expsyn-B-all")) self.assertEqual(1, len(m)) @@ -141,7 +141,7 @@ def test_probe_addr_metadata(self): m = sim.probe_metadata((0, "ina")) self.assertEqual(1, len(m)) - self.assertEqual(A.location(0, 0.10), m[0]) + self.assertEqual([A.location(0, 0.10)], m[0]) m = sim.probe_metadata((0, "ina-all")) self.assertEqual(1, len(m)) @@ -149,7 +149,7 @@ def test_probe_addr_metadata(self): m = sim.probe_metadata((0, "nai")) self.assertEqual(1, len(m)) - self.assertEqual(A.location(0, 0.12), m[0]) + self.assertEqual([A.location(0, 0.12)], m[0]) m = sim.probe_metadata((0, "nai-all")) self.assertEqual(1, len(m)) @@ -157,7 +157,7 @@ def test_probe_addr_metadata(self): m = sim.probe_metadata(0, "nao") self.assertEqual(1, len(m)) - self.assertEqual(A.location(0, 0.14), m[0]) + self.assertEqual([A.location(0, 0.14)], m[0]) m = sim.probe_metadata((0, "nao-all")) self.assertEqual(1, len(m)) @@ -198,10 +198,9 @@ class TestLifProbes(unittest.TestCase): def test_probe_addr_metadata(self): rec = lif_recipe() sim = A.simulation(rec) - - m = sim.probe_metadata((0, "Um")) - self.assertEqual(1, len(m)) - self.assertTrue(all(isinstance(i, A.lif_probe_metadata) for i in m)) + ms = sim.probe_metadata((0, "Um")) + self.assertEqual(1, len(ms)) + self.assertTrue(all(isinstance(i, A.lif_probe_metadata) for m in ms for i in m)) def test_probe_result(self): rec = lif_recipe() diff --git a/test/simple_recipes.hpp b/test/simple_recipes.hpp index 52b6638f41..233a189993 100644 --- a/test/simple_recipes.hpp +++ b/test/simple_recipes.hpp @@ -21,8 +21,7 @@ namespace U = units; // Common functionality: maintain an unordered map of probe data // per gid, built with `add_probe()`. -class simple_recipe_base: public recipe { -public: +struct simple_recipe_base: public recipe { simple_recipe_base() { cell_gprop_.catalogue = global_default_catalogue(); cell_gprop_.default_parameters = neuron_parameter_defaults; @@ -94,8 +93,7 @@ class homogeneous_recipe: public simple_recipe_base { // // Cell descriptions passed to the constructor are cloned. -class cable1d_recipe: public simple_recipe_base { -public: +struct cable1d_recipe: public simple_recipe_base { template explicit cable1d_recipe(const Seq& cells, bool coalesce = true) { for (const auto& c: cells) { diff --git a/test/unit/test_diffusion.cpp b/test/unit/test_diffusion.cpp index 8f51e38faa..99e6b8c44a 100644 --- a/test/unit/test_diffusion.cpp +++ b/test/unit/test_diffusion.cpp @@ -104,16 +104,15 @@ testing::AssertionResult all_near(const result_t& a, const result_t& b, double e testing::AssertionResult run(const linear& rec, const result_t exp) { result_t sample_values; - auto sampler = [&sample_values](arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { + auto sampler = [&sample_values](arb::probe_metadata pm, const arb::sample_records& recs) { sample_values.clear(); - auto ptr = arb::util::any_cast(pm.meta); - ASSERT_NE(ptr, nullptr); - auto n_cable = ptr->size(); - for (std::size_t i = 0; i(samples[i].data); - for (unsigned j = 0; j(pm.meta, recs); + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + auto time = reader.time(ix); + for (std::size_t iy = 0; iy < reader.n_column(); ++iy) { + auto cable = reader.metadata(iy); + auto value = reader.value(ix, iy); + sample_values.push_back({time, cable.prox_pos, value}); } } }; diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp index 3772f223b0..3f616af246 100644 --- a/test/unit/test_fvm_lowered.cpp +++ b/test/unit/test_fvm_lowered.cpp @@ -464,9 +464,10 @@ TEST(fvm_lowered, derived_mechs) { std::vector samples[3]; sampler_function sampler = - [&](probe_metadata pm, std::size_t n, const sample_record* records) { - for (std::size_t i = 0; i(records[i].data); + [&](probe_metadata pm, const sample_records& records) { + auto reader = sample_reader(pm.meta, records); + for (std::size_t i = 0; i < reader.n_row(); ++i) { + double v = reader.value(i); samples[pm.id.gid].push_back(v); } }; diff --git a/test/unit/test_lif_cell_group.cpp b/test/unit/test_lif_cell_group.cpp index 25829ad95c..0d21558cae 100644 --- a/test/unit/test_lif_cell_group.cpp +++ b/test/unit/test_lif_cell_group.cpp @@ -246,13 +246,14 @@ struct Um_type { TEST(lif_cell_group, probe) { auto ums = std::unordered_map>{}; - auto fun = [&ums](probe_metadata pm, - std::size_t n, - const sample_record* samples) { - for (std::size_t ix = 0; ix < n; ++ix) { - const auto& [t, v] = samples[ix]; - double u = *util::any_cast(v); - ums[pm.id].push_back({t, u}); + auto fun = [&ums](const probe_metadata& pm, + const sample_records& samples) { + using probe_t = arb::lif_probe_voltage; + auto reader = arb::sample_reader(pm.meta, samples); + for (std::size_t ix = 0ul; ix < reader.n_row(); ++ix) { + auto t = reader.time(ix); + auto v = reader.value(ix); + ums[pm.id].push_back({t, v}); } }; auto rec = probe_recipe{}; @@ -683,12 +684,13 @@ TEST(lif_cell_group, probe) { TEST(lif_cell_group, probe_with_connections) { auto ums = std::unordered_map>{}; auto fun = [&ums](probe_metadata pm, - std::size_t n, - const sample_record* samples) { - for (std::size_t ix = 0; ix < n; ++ix) { - const auto& [t, v] = samples[ix]; - double u = *util::any_cast(v); - ums[pm.id].push_back({t, u}); + const sample_records& recs) { + using meta_t = lif_probe_voltage::meta_type; + auto reader = arb::sample_reader(pm.meta, recs); + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + auto time = reader.time(ix); + auto value = reader.value(ix); + ums[pm.id].push_back({time, value}); } }; auto rec = probe_recipe{5}; diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp index 489f9adaa5..ee9d59ab29 100644 --- a/test/unit/test_probe.cpp +++ b/test/unit/test_probe.cpp @@ -139,9 +139,9 @@ void run_v_i_probe_test(context ctx) { // implemented as an interpolated probe too, in order to account // for stimulus contributions. - ASSERT_TRUE(std::get_if(&probe_map.data_on({0, "Um-l0"}).front()->info)); - ASSERT_TRUE(std::get_if(&probe_map.data_on({0, "Um-l1"}).front()->info)); - ASSERT_TRUE(std::get_if(&probe_map.data_on({0, "Ii-l2"}).front()->info)); + ASSERT_TRUE(std::get_if(&probe_map.data_on({0, "Um-l0"}).front()->info)); + ASSERT_TRUE(std::get_if(&probe_map.data_on({0, "Um-l1"}).front()->info)); + ASSERT_TRUE(std::get_if(&probe_map.data_on({0, "Ii-l2"}).front()->info)); probe_handle p0a = get_probe_raw_handle({0, "Um-l0"}, 0); probe_handle p0b = get_probe_raw_handle({0, "Um-l0"}, 1); @@ -576,11 +576,10 @@ void run_ion_density_probe_test(context ctx) { const auto& probe_map = fvm_info.probe_map; - // Should be no sodium ion instantiated on CV 0, so probe (0, 6) should + // Should be no sodium ion instantiated on CV 0, so probe (0, nai-l0) should // have been silently discared. Similarly, write_ca2 is not instantiated on - // CV 0, and so probe (0, 8) should have been discarded. All other probes - // should be in the map. - + // CV 0, and so probe (0, state-s-l0) should have been discarded. All other + // probes should be in the map. EXPECT_EQ(13u, rec.get_probes(0).size()); EXPECT_EQ(11u, probe_map.size()); @@ -754,7 +753,6 @@ void run_partial_density_probe_test(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.get_probes(0).size()); EXPECT_EQ(10u, rec.get_probes(1).size()); EXPECT_EQ(10u, probe_map.size()); @@ -849,45 +847,32 @@ void run_axial_and_ion_current_sampled_probe_test(context ctx) { std::vector i_memb(n_cv), i_stim(n_cv); sim.add_sampler(all_probes, explicit_schedule(std::vector{20*tau}), - [&](probe_metadata pm, - std::size_t n_sample, - const sample_record* samples) { + [&](probe_metadata pm, const sample_records& recs) { // Expect exactly one sample. - ASSERT_EQ(1u, n_sample); + ASSERT_EQ(1u, recs.n_sample); if (pm.id.tag == "I-total") { - auto m = any_cast(pm.meta); - ASSERT_NE(nullptr, m); + using probe_t = cable_probe_ion_current_cell; + auto reader = arb::sample_reader(pm.meta, recs); // Metadata should comprise one cable per CV. - ASSERT_EQ(n_cv, m->size()); - const cable_sample_range* s = any_cast(samples->data); - ASSERT_NE(nullptr, s); - auto [s_beg, s_end] = *s; - ASSERT_EQ(s_end - s_beg, n_cv); - for (unsigned ix = 0; ix < n_cv; ++ix) i_memb[ix] = *s_beg++; + ASSERT_EQ(n_cv, reader.n_column()); + for (unsigned iy = 0; iy < n_cv; ++iy) i_memb[iy] = reader.value(0, iy); } else if (pm.id.tag == "I-stimulus") { - auto m = any_cast(pm.meta); - ASSERT_NE(nullptr, m); + using probe_t = cable_probe_stimulus_current_cell; + auto reader = arb::sample_reader(pm.meta, recs); // Metadata should comprise one cable per CV. - ASSERT_EQ(n_cv, m->size()); - const cable_sample_range* s = any_cast(samples->data); - ASSERT_NE(nullptr, s); - auto [s_beg, s_end] = *s; - ASSERT_EQ(s_end - s_beg, n_cv); - for (unsigned ix = 0; ix < n_cv; ++ix) i_stim[ix] = *s_beg++; + ASSERT_EQ(n_cv, reader.n_column()); + for (unsigned iy = 0; iy < n_cv; ++iy) i_stim[iy] = reader.value(0, iy); } else { + ASSERT_EQ(recs.width, 1u); // Probe id tells us which axial current this is. ASSERT_LT(indices.count(pm.id.tag), n_axial_probe); auto idx = indices.at(pm.id.tag); - auto m = any_cast(pm.meta); - ASSERT_NE(nullptr, m); - - auto s = any_cast(samples->data); - ASSERT_NE(nullptr, s); - - i_axial.at(idx) = *s; + using probe_t = cable_probe_axial_current; + auto reader = arb::sample_reader(pm.meta, recs); + i_axial.at(idx) = reader.value(0, 0); } }); @@ -916,41 +901,28 @@ void run_axial_and_ion_current_sampled_probe_test(context ctx) { // // Use the default mechanism catalogue augmented by unit test specific mechanisms. // (Timestep fixed at 0.025 ms). - -template -auto run_simple_samplers(const arb::context& ctx, +template +auto run_simple_sampler(const arb::context& ctx, U::quantity t_end, const std::vector& cells, - const cell_address_type& probe, - const std::vector& probe_addrs, + const cell_address_type& where, + const P& probe, const std::vector& when) { cable1d_recipe rec(cells, false); rec.catalogue() = make_unit_test_catalogue(global_default_catalogue()); - unsigned n_probe = probe_addrs.size(); - for (auto& addr: probe_addrs) rec.add_probe(probe.gid, probe.tag, addr); + + rec.add_probe(where.gid, where.tag, probe); partition_hint_map phints = {{cell_kind::cable, {partition_hint::max_size, partition_hint::max_size, true}}}; simulation sim(rec, ctx, partition_load_balance(rec, ctx, phints)); - std::vector> traces(n_probe); - for (unsigned i = 0; i trace; + sim.add_sampler(one_probe(where), + explicit_schedule(when), + make_simple_sampler(trace)); sim.run(t_end, 0.025*U::ms); - return traces; -} - -template -auto run_simple_sampler(const arb::context& ctx, - U::quantity t_end, - const std::vector& cells, - const cell_address_type& probe, - const std::any& probe_addr, - const std::vector& when) { - return run_simple_samplers(ctx, t_end, cells, probe, {probe_addr}, when).at(0); + return trace; } template @@ -960,27 +932,25 @@ void run_multi_probe_test(context ctx) { // m_mlt_b6 has terminal branches 1, 2, 4, and 5. auto m = common_morphology::m_mlt_b6; - decor d; - // Paint mechanism on branches 1, 2, and 5, omitting branch 4. - d.paint(reg::branch(1), density("param_as_state", {{"p", 10.}})); - d.paint(reg::branch(2), density("param_as_state", {{"p", 20.}})); - d.paint(reg::branch(5), density("param_as_state", {{"p", 50.}})); - - auto tracev = run_simple_sampler(ctx, 0.1*U::ms, - {cable_cell{m, d}}, - {0, "probe"}, - cable_probe_density_state{ls::terminal(), "param_as_state", "s"}, - {0.0*U::ms}); + auto d = decor{} + .paint(reg::branch(1), density("param_as_state", {{"p", 10.}})) + .paint(reg::branch(2), density("param_as_state", {{"p", 20.}})) + .paint(reg::branch(5), density("param_as_state", {{"p", 50.}})); + + auto trace = run_simple_sampler(ctx, 0.1*U::ms, {cable_cell{m, d}}, + {0, "probe"}, cable_probe_density_state{ ls::terminal(), "param_as_state", "s"}, + {0.0*U::ms}); + ASSERT_EQ(3u, trace.values.size()); + ASSERT_EQ(3u, trace.metadata.size()); + for (const auto& val: trace.values) ASSERT_EQ(1u, val.size()); // Expect to have received a sample on each of the terminals of branches 1, 2, and 5. - ASSERT_EQ(3u, tracev.size()); - - // Expect a single sample per terminal std::vector> vals; - for (auto& trace: tracev) { - // ASSERT_EQ(1u, trace.size()); - vals.push_back({trace.meta, trace[0].v}); + for (size_t ix = 0; ix < trace.n_sample; ++ix) { + for (size_t iy = 0; iy < trace.width; ++iy) { + vals.emplace_back(trace.metadata.at(iy), trace.values[iy][ix]); + } } util::sort(vals); @@ -1016,27 +986,23 @@ void run_v_sampled_probe_test(context ctx) { const auto t_end = 1.*U::ms; // [ms] std::vector when = {0.3*U::ms, 0.6*U::ms}; // Sample at 0.3 and 0.6 ms. - auto trace0 = run_simple_sampler(ctx, t_end, cells, {0, "Um-loc"}, - cable_probe_membrane_voltage{probe_loc}, - when).at(0); - ASSERT_TRUE(trace0); - - EXPECT_EQ(probe_loc, trace0.meta); - EXPECT_EQ(2u, trace0.size()); + auto probe = cable_probe_membrane_voltage{ls::location(probe_loc)}; - auto trace1 = run_simple_sampler(ctx, t_end, cells, {1, "Um-loc"}, - cable_probe_membrane_voltage{probe_loc}, - when).at(0); - ASSERT_TRUE(trace1); + auto trace0 = run_simple_sampler(ctx, t_end, cells, {0, "Um-loc"}, probe, when); + EXPECT_EQ(probe_loc, trace0.metadata.at(0)); + EXPECT_EQ(1u, trace0.values.size()); + EXPECT_EQ(2u, trace0.values[0].size()); - EXPECT_EQ(probe_loc, trace1.meta); - EXPECT_EQ(2u, trace1.size()); + auto trace1 = run_simple_sampler(ctx, t_end, cells, {1, "Um-loc"}, probe, when); + EXPECT_EQ(probe_loc, trace1.metadata.at(0)); + EXPECT_EQ(1u, trace1.values.size()); + EXPECT_EQ(2u, trace1.values[0].size()); - EXPECT_EQ(trace0[0].t, trace1[0].t); - EXPECT_EQ(trace0[0].v, trace1[0].v); + EXPECT_EQ(trace0.time[0], trace1.time[0]); + EXPECT_EQ(trace0.values[0][0], trace1.values[0][0]); - EXPECT_EQ(trace0[1].t, trace1[1].t); - EXPECT_NE(trace0[1].v, trace1[1].v); + EXPECT_EQ(trace0.time[1], trace1.time[1]); + EXPECT_NE(trace0.values[0][1], trace1.values[0][1]); } @@ -1074,41 +1040,42 @@ void run_total_current_probe_test(context ctx) { // at the fork points, and once with a non-trivial CV centred on the fork // point. - trace_data, mcable_list> traces[2]; - trace_data, mcable_list> ion_traces[2]; - trace_data, mcable_list> stim_traces[2]; + std::vector> traces(2); + std::vector> ion_traces(2); + std::vector> stim_traces(2); // Run the cells sampling at τ and 20τ for both total membrane // current and total membrane ionic current. - auto run_cells = [&](bool interior_forks) { auto flags = interior_forks? cv_policy_flag::interior_forks: cv_policy_flag::none; cv_policy policy = cv_policy_fixed_per_branch(n_cv_per_branch, flags); std::vector cells = {{m, d0, {}, policy}, {m, d1, {}, policy}}; - for (unsigned i = 0; i<2; ++i) { + for (unsigned i = 0; i < 2; ++i) { SCOPED_TRACE(i); + auto& trace = traces[i]; + auto& ion_trace = ion_traces[i]; + auto& stim_trace = stim_traces[i]; + auto t_end = 21*tau; // [ms] - traces[i] = run_simple_sampler, mcable_list>(ctx, t_end, cells, - {i, "Itotal"}, - cable_probe_total_current_cell{}, - {tau, 20*tau}).at(0); + trace = run_simple_sampler(ctx, t_end, cells, + {i, "Itotal"}, cable_probe_total_current_cell{}, + {tau, 20*tau}); - ion_traces[i] = run_simple_sampler, mcable_list>(ctx, t_end, cells, - {i, "Iion"}, - cable_probe_total_ion_current_cell{}, - {tau, 20*tau}).at(0); + ion_trace = run_simple_sampler(ctx, t_end, cells, + {i, "Iion"}, cable_probe_total_ion_current_cell{}, + {tau, 20*tau}); - stim_traces[i] = run_simple_sampler, mcable_list>(ctx, t_end, cells, - {i, "Istim"}, - cable_probe_stimulus_current_cell{}, - {tau, 20*tau}).at(0); + stim_trace = run_simple_sampler(ctx, t_end, cells, + {i, "Istim"}, cable_probe_stimulus_current_cell{}, + {tau, 20*tau}); - ASSERT_EQ(2u, traces[i].size()); - ASSERT_EQ(2u, ion_traces[i].size()); + ASSERT_EQ(2u, trace.n_sample); + ASSERT_EQ(2u, ion_trace.n_sample); + ASSERT_EQ(2u, stim_trace.n_sample); // Check metadata size: // * With trivial forks, should have n_cv_per_branch*n_branch cables; zero-length cables @@ -1118,23 +1085,24 @@ void run_total_current_probe_test(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].width); - ASSERT_EQ((n_cv_per_branch+(int)interior_forks)*n_branch, traces[i].meta.size()); - ASSERT_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()); + ASSERT_EQ(ion_trace.width, trace.width); + ASSERT_EQ(ion_trace.metadata, trace.metadata); + for (std::size_t ix = 0; ix < trace.width; ++ix) { + EXPECT_EQ(ion_trace.values[ix].size(), trace.values[ix].size()); + } // Check total membrane currents + stimulus currents are individually non-zero, but sum is, // both at t=τ (j=0) and t=20τ (j=1). - ASSERT_EQ(ion_traces[i].meta, stim_traces[i].meta); - - for (unsigned j: {0u, 1u}) { + ASSERT_EQ(ion_traces[i].metadata, stim_traces[i].metadata); + for (unsigned j: util::make_span(2)) { double max_abs_current = 0; double sum_current = 0; - for (auto k: util::count_along(traces[i].meta)) { - double current = traces[i][j].v[k] + stim_traces[i][j].v[k]; - + + for (auto k: util::make_span(trace.width)) { + double current = trace.values[k][j] + stim_trace.values[k][j]; EXPECT_NE(0.0, current); max_abs_current = std::max(max_abs_current, std::abs(current)); sum_current += current; @@ -1144,23 +1112,22 @@ void run_total_current_probe_test(context ctx) { } // Confirm that total (non-stim) and ion currents differ at τ but are close at 20τ. - - for (unsigned k = 0; k, mcable_list> traces[2]; + + using probe_t = cable_probe_stimulus_current_cell; + simple_sampler_result traces[2]; for (unsigned i: {0u, 1u}) { - traces[i] = run_simple_sampler, mcable_list>(ctx, 2.5*stim_until, cells, {i, "Istim"}, - cable_probe_stimulus_current_cell{}, - {0.5*stim_until, 2*stim_until}).at(0); + traces[i] = run_simple_sampler(ctx, 2.5*stim_until, cells, + {i, "Istim"}, probe_t{}, + {0.5*stim_until, 2*stim_until}); - ASSERT_EQ(3u, traces[i].meta.size()); + ASSERT_EQ(3u, traces[i].width); for ([[maybe_unused]] unsigned cv: {0u, 1u, 2u}) { - ASSERT_EQ(2u, traces[i].size()); + ASSERT_EQ(2u, traces[i].values[3].size()); } } // Every sample in each trace should be zero _except_ the first sample for cell 0, cv 1 // and the first sample for cell 1, cv 2. - - EXPECT_EQ((std::vector{0, expected_stim0, 0}), traces[0][0]); - EXPECT_EQ((std::vector{0, 0, expected_stim1}), traces[1][0]); - EXPECT_EQ((std::vector(3)), traces[0][1]); - EXPECT_EQ((std::vector(3)), traces[1][1]); + EXPECT_EQ((std::vector{0, expected_stim0, 0}), traces[0].values[0]); + EXPECT_EQ((std::vector{0, 0, expected_stim1}), traces[1].values[0]); + EXPECT_EQ((std::vector(3)), traces[0].values[1]); + EXPECT_EQ((std::vector(3)), traces[1].values[1]); } /* @@ -1384,20 +1352,18 @@ ARB_PP_FOREACH(RUN_GPU, PROBE_TESTS) // Test simulator `get_probe_metadata` interface. // (No need to run this on GPU back-end as well.) - TEST(probe, get_probe_metadata) { // Reuse multiprobe test set-up to confirm simulator::get_probe_metadata returns // correct vector of metadata. - auto m = common_morphology::m_mlt_b6; - decor d; - - // Paint mechanism on branches 1, 2, and 5, omitting branch 4. - d.paint(reg::branch(1), density("param_as_state", {{"p", 10.}})); - d.paint(reg::branch(2), density("param_as_state", {{"p", 20.}})); - d.paint(reg::branch(5), density("param_as_state", {{"p", 50.}})); + auto mrf = common_morphology::m_mlt_b6; + auto dec = decor{} + // Paint mechanism on branches 1, 2, and 5, omitting branch 4. + .paint(reg::branch(1), density("param_as_state", {{"p", 10.}})) + .paint(reg::branch(2), density("param_as_state", {{"p", 20.}})) + .paint(reg::branch(5), density("param_as_state", {{"p", 50.}})); - cable1d_recipe rec(cable_cell{m, d}, false); + cable1d_recipe rec(cable_cell{mrf, dec}, false); rec.catalogue() = make_unit_test_catalogue(global_default_catalogue()); rec.add_probe(0, "param-as-state", cable_probe_density_state{ls::terminal(), "param_as_state", "s"}); @@ -1408,25 +1374,15 @@ TEST(probe, get_probe_metadata) { simulation sim(rec, ctx, partition_load_balance(rec, ctx, phints)); std::vector mm = sim.get_probe_metadata({0, "param-as-state"}); - ASSERT_EQ(3u, mm.size()); + ASSERT_EQ(1u, mm.size()); EXPECT_EQ((cell_address_type{0, "param-as-state"}), mm[0].id); - EXPECT_EQ((cell_address_type{0, "param-as-state"}), mm[1].id); - EXPECT_EQ((cell_address_type{0, "param-as-state"}), mm[2].id); EXPECT_EQ(0u, mm[0].index); - EXPECT_EQ(1u, mm[1].index); - EXPECT_EQ(2u, mm[2].index); - - const mlocation* l0 = any_cast(mm[0].meta); - const mlocation* l1 = any_cast(mm[1].meta); - const mlocation* l2 = any_cast(mm[2].meta); - - ASSERT_TRUE(l0); - ASSERT_TRUE(l1); - ASSERT_TRUE(l2); - std::vector locs = {*l0, *l1, *l2}; + // TODO This isn't as nice as we'd like. + auto ptr = any_cast(mm[0].meta); + std::vector locs(ptr, ptr + 3); util::sort(locs); EXPECT_EQ((mlocation{1, 1.}), locs[0]); EXPECT_EQ((mlocation{2, 1.}), locs[1]); diff --git a/test/unit/test_sde.cpp b/test/unit/test_sde.cpp index e0f097db9f..f49992f7e0 100644 --- a/test/unit/test_sde.cpp +++ b/test/unit/test_sde.cpp @@ -426,7 +426,7 @@ TEST(sde, reproducibility) { // - 1 density processes with 1 random variable std::size_t n_rv_densities = ncells*ncvs; std::size_t n_rv_per_dt = n_rv_densities; - std::size_t n_rv = n_rv_per_dt*(nsteps); + std::size_t n_rv = n_rv_per_dt*nsteps; // setup storage std::vector data; @@ -520,7 +520,7 @@ TEST(sde, normality) { std::size_t n_rv_synapses = ncells*nsynapses*(1+1+2+2); std::size_t n_rv_densities = ncells*ncvs*(1+1+2+2); std::size_t n_rv_per_dt = n_rv_synapses + n_rv_densities; - std::size_t n_rv = n_rv_per_dt*(nsteps); + std::size_t n_rv = n_rv_per_dt*nsteps; // setup storage std::vector data; @@ -667,48 +667,45 @@ TEST(sde, solver) { dec.place(*labels.locset("locs"), synapse(m4), "m4"); // a basic sampler: stores result in a vector - auto sampler_ = [nsteps] (std::vector& results, unsigned count, - probe_metadata pm, std::size_t n, sample_record const * samples) { - - auto* point_info_ptr = arb::util::any_cast*>(pm.meta); - assert(point_info_ptr != nullptr); - - unsigned n_entities = point_info_ptr->size(); - assert(n_entities == count); - - unsigned offset = pm.id.gid*(nsteps)*n_entities; - unsigned stride = n_entities; - assert(n == nsteps); - for (std::size_t i = 0; i(samples[i].data); - assert(value_range); - const auto& [lo, hi] = *value_range; - assert(n_entities==hi-lo); - for (unsigned j = 0; j& results, const probe_metadata& pm, const sample_records& samples) { + std::size_t n_entities = samples.width; + std::size_t offset = pm.id.gid*nsteps*n_entities; + std::size_t stride = n_entities; + assert(samples.n_sample == nsteps); + + using probe_t = arb::cable_probe_point_state_cell; + auto reader = arb::sample_reader(pm.meta, samples); + for (std::size_t ix = 0; ix results_m1(ncells*(nsteps)*(nsynapses)); - auto sampler_m1 = [&] (probe_metadata pm, std::size_t n, sample_record const * samples) { - sampler_(results_m1, nsynapses, pm, n, samples); + std::vector results_m1(ncells*nsteps*nsynapses); + auto sampler_m1 = [&] (const probe_metadata& pm, const sample_records& samples) { + ASSERT_EQ(nsynapses, samples.width); + sampler_(results_m1, pm, samples); }; // concrete sampler for process m2 - std::vector results_m2(ncells*(nsteps)*(nsynapses)); - auto sampler_m2 = [&] (probe_metadata pm, std::size_t n, sample_record const * samples) { - sampler_(results_m2, nsynapses, pm, n, samples); + std::vector results_m2(ncells*nsteps*nsynapses); + auto sampler_m2 = [&] (const probe_metadata& pm, const sample_records& samples) { + ASSERT_EQ(nsynapses, samples.width); + sampler_(results_m2, pm, samples); }; // concrete sampler for process m3 - std::vector results_m3(ncells*(nsteps)*(nsynapses)); - auto sampler_m3 = [&] (probe_metadata pm, std::size_t n, sample_record const * samples) { - sampler_(results_m3, nsynapses, pm, n, samples); + std::vector results_m3(ncells*nsteps*nsynapses); + auto sampler_m3 = [&] (const probe_metadata& pm, const sample_records& samples) { + ASSERT_EQ(nsynapses, samples.width); + sampler_(results_m3, pm, samples); }; // concrete sampler for process m4 - std::vector results_m4(ncells*(nsteps)*(nsynapses)); - auto sampler_m4 = [&] (probe_metadata pm, std::size_t n, sample_record const * samples) { - sampler_(results_m4, nsynapses, pm, n, samples); + std::vector results_m4(ncells*nsteps*nsynapses); + auto sampler_m4 = [&] (const probe_metadata& pm, const sample_records& samples) { + ASSERT_EQ(nsynapses, samples.width); + sampler_(results_m4, pm, samples); }; // instantiate recipe @@ -747,11 +744,11 @@ TEST(sde, solver) { // accumulate statistics for sampled data for (unsigned int k=0; k& results, unsigned count, - probe_metadata pm, std::size_t n, sample_record const * samples) { - - auto* point_info_ptr = arb::util::any_cast*>(pm.meta); - assert(point_info_ptr != nullptr); - - unsigned n_entities = point_info_ptr->size(); - assert(n_entities == count); - - unsigned offset = pm.id.gid*(nsteps)*n_entities; - unsigned stride = n_entities; - assert(n == nsteps); - for (std::size_t i = 0; i(samples[i].data); - assert(value_range); - const auto& [lo, hi] = *value_range; - assert(n_entities==hi-lo); - for (unsigned j = 0; j& results, + const probe_metadata& pm, + const sample_records& samples) { + std::size_t n_entities = samples.width; + std::size_t offset = pm.id.gid*nsteps*n_entities; + std::size_t stride = n_entities; + assert(samples.n_sample == nsteps); + + using probe_t = arb::cable_probe_point_state_cell; + auto reader = arb::sample_reader(pm.meta, samples); + for (std::size_t ix = 0; ix results_P(ncells*(nsteps)*(nsynapses)); - auto sampler_P = [&] (probe_metadata pm, std::size_t n, sample_record const * samples) { - sampler_(results_P, nsynapses, pm, n, samples); + std::vector results_P(ncells*nsteps*nsynapses); + auto sampler_P = [&] (probe_metadata pm, const sample_records& samples) { + ASSERT_EQ(nsynapses, samples.width); + sampler_(results_P, pm, samples); }; // concrete sampler for sigma - std::vector results_sigma(ncells*(nsteps)*(nsynapses)); - auto sampler_sigma = [&] (probe_metadata pm, std::size_t n, sample_record const * samples) { - sampler_(results_sigma, nsynapses, pm, n, samples); + std::vector results_sigma(ncells*nsteps*nsynapses); + auto sampler_sigma = [&] (probe_metadata pm, const sample_records& samples) { + ASSERT_EQ(nsynapses, samples.width); + sampler_(results_sigma, pm, samples); }; // instantiate recipe @@ -881,9 +875,9 @@ TEST(sde, coupled) { // accumulate statistics for sampled data for (unsigned int k=0; k arch_cpu(n_rv); diff --git a/test/unit/test_serdes.cpp b/test/unit/test_serdes.cpp index 546362daac..25ca2f6e18 100644 --- a/test/unit/test_serdes.cpp +++ b/test/unit/test_serdes.cpp @@ -179,15 +179,13 @@ struct serdes_recipe: public arb::recipe { static std::vector* output; void sampler(arb::probe_metadata pm, - std::size_t n, - const arb::sample_record* samples) { - auto* loc = arb::util::any_cast(pm.meta); - auto* point_info = arb::util::any_cast(pm.meta); - loc = loc ? loc : &(point_info->loc); - - for (std::size_t i = 0; i(samples[i].data); - output->push_back(*value); + const arb::sample_records& samples) { + using probe_t = arb::cable_probe_membrane_voltage; + auto reader = arb::sample_reader(pm.meta, samples); + assert(reader.n_column() == 1); + for (std::size_t ix = 0; ix < reader.n_row(); ++ix) { + auto value = reader.value(ix); + output->push_back(value); } } diff --git a/test/unit/test_v_clamp.cpp b/test/unit/test_v_clamp.cpp index cb5dd6d460..9eb74b66fe 100644 --- a/test/unit/test_v_clamp.cpp +++ b/test/unit/test_v_clamp.cpp @@ -72,17 +72,18 @@ using um_s_type = std::vector; TEST(v_process, clamp) { auto u_soma = um_s_type{}; auto u_dend = um_s_type{}; - auto fun = [&u_soma, &u_dend](arb::probe_metadata pm, - std::size_t n, - const arb::sample_record* samples) { - for (std::size_t ix = 0ul; ix < n; ++ix) { - const auto& [t, v] = samples[ix]; - double u = *arb::util::any_cast(v); + auto fun = [&u_soma, &u_dend](const arb::probe_metadata& pm, + const arb::sample_records& samples) { + using probe_t = arb::cable_probe_membrane_voltage; + auto reader = arb::sample_reader(pm.meta, samples); + for (std::size_t ix = 0ul; ix < reader.n_row(); ++ix) { + auto t = reader.time(ix); + auto v = reader.value(ix); if (pm.id.tag == "Um-(0, 0.125)") { - u_soma.push_back({t, u}); + u_soma.push_back({t, v}); } else if (pm.id.tag == "Um-(0, 0.625)") { - u_dend.push_back({t, u}); + u_dend.push_back({t, v}); } else { throw std::runtime_error{"Unexpected probe id"}; @@ -140,17 +141,18 @@ TEST(v_process, clamp) { TEST(v_process, limit) { auto u_soma = um_s_type{}; auto u_dend = um_s_type{}; - auto fun = [&u_soma, &u_dend](arb::probe_metadata pm, - std::size_t n, - const arb::sample_record* samples) { - for (std::size_t ix = 0ul; ix < n; ++ix) { - const auto& [t, v] = samples[ix]; - double u = *arb::util::any_cast(v); + auto fun = [&u_soma, &u_dend](const arb::probe_metadata& pm, + const arb::sample_records& samples) { + using probe_t = arb::cable_probe_membrane_voltage; + auto reader = arb::sample_reader(pm.meta, samples); + for (std::size_t ix = 0ul; ix < reader.n_row(); ++ix) { + auto t = reader.time(ix); + auto v = reader.value(ix); if (pm.id.tag == "Um-(0, 0.125)") { - u_soma.push_back({t, u}); + u_soma.push_back({t, v}); } else if (pm.id.tag == "Um-(0, 0.625)") { - u_dend.push_back({t, u}); + u_dend.push_back({t, v}); } else { throw std::runtime_error{"Unexpected probe id"}; @@ -208,17 +210,18 @@ TEST(v_process, limit) { TEST(v_process, clamp_fine) { auto u_soma = um_s_type{}; auto u_dend = um_s_type{}; - auto fun = [&u_soma, &u_dend](arb::probe_metadata pm, - std::size_t n, - const arb::sample_record* samples) { - for (std::size_t ix = 0ul; ix < n; ++ix) { - const auto& [t, v] = samples[ix]; - double u = *arb::util::any_cast(v); + auto fun = [&u_soma, &u_dend](const arb::probe_metadata& pm, + const arb::sample_records& samples) { + using probe_t = arb::cable_probe_membrane_voltage; + auto reader = arb::sample_reader(pm.meta, samples); + for (std::size_t ix = 0ul; ix < reader.n_row(); ++ix) { + auto t = reader.time(ix); + auto v = reader.value(ix); if (pm.id.tag == "Um-(0, 0.125)") { - u_soma.push_back({t, u}); + u_soma.push_back({t, v}); } else if (pm.id.tag == "Um-(0, 0.625)") { - u_dend.push_back({t, u}); + u_dend.push_back({t, v}); } else { throw std::runtime_error{"Unexpected probe id"};