diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 000000000..9e5281ed7 --- /dev/null +++ b/.editorconfig @@ -0,0 +1,14 @@ +root = true + +# Unix-style newlines with a newline ending every file +[*] +end_of_line = lf +insert_final_newline = true + +[*.py] +charset = utf-8 + +# 4 space indentation +[*{.py,.cpp,.cu,.hpp}] +indent_style = space +indent_size = 4 diff --git a/.github/workflows/test-matrix.yml b/.github/workflows/test-matrix.yml index 21b93c3c4..264e96b86 100644 --- a/.github/workflows/test-matrix.yml +++ b/.github/workflows/test-matrix.yml @@ -29,8 +29,8 @@ jobs: - { name: "Linux Min Clang", os: "ubuntu-22.04", - cc: "clang-13", - cxx: "clang++-13", + cc: "clang-14", + cxx: "clang++-14", py: "3.10", cmake: "4.0.x", mpi: "ON", diff --git a/arbor/adex_cell_group.cpp b/arbor/adex_cell_group.cpp index 7553eac76..55d07854c 100644 --- a/arbor/adex_cell_group.cpp +++ b/arbor/adex_cell_group.cpp @@ -3,6 +3,7 @@ #include #include "arbor/math.hpp" +#include "util/maputil.hpp" #include "util/rangeutil.hpp" #include "util/span.hpp" #include "label_resolution.hpp" @@ -99,6 +100,49 @@ void adex_cell_group::reset() { spikes_.clear(); } +void +adex_cell_group::edit_cell(cell_gid_type gid, std::any cell_edit) { + try { + auto adex_edit = std::any_cast(cell_edit); + auto lid = util::binary_search_index(gids_, gid); + if (!lid) throw arb::arbor_internal_error{"gid " + std::to_string(gid) + " erroneuosly dispatched to cell group."}; + auto& lowered = cells_[*lid]; + auto tmp = adex_cell { + .source = lowered.source, // Label of source + .target = lowered.target, // Label of target + .delta = lowered.delta * U::mV, + .V_th = lowered.V_th * U::mV, + .C_m = lowered.C_m * U::nF, + .E_L = lowered.E_L * U::mV, + .E_R = lowered.E_R * U::mV, + .V_m = lowered.V_m * U::mV, + .t_ref = lowered.t_ref * U::ms, + .g = lowered.g * U::uS, + .tau = lowered.tau * U::ms, + .w = lowered.w * U::pA, + .a = lowered.a * U::uS, + .b = lowered.b * U::nA, + }; + adex_edit(tmp); + // NOTE: we forbid writing to V_m? Reasons + // * the cell might be in the refractory period which causes semantic issues + // - return to normal or not? + // - what should probes return + // * V_m is the _initial state_ only + if (tmp.V_m.value_as(U::mV) != lowered.V_m) throw bad_cell_edit(gid, "Initial voltage is not editable."); + if (tmp.w.value_as(U::pA) != lowered.w) throw bad_cell_edit(gid, "Adaption parameter is not editable."); + if (tmp.source != lowered.source) throw bad_cell_edit(gid, "Source is not editable."); + if (tmp.target != lowered.target) throw bad_cell_edit(gid, "Target is not editable."); + // Write back + lowered = adex_lowered_cell{tmp}; + + } + catch (const std::bad_any_cast& ){ + throw bad_cell_edit(gid, "Not an AdEx editor (C++ type-id: '" + std::string{cell_edit.type().name()} + "')"); + } +} + + // integrate a single cell's state from current time `cur` tos final time `end`. // Extra parameters // * the cell cannot be updated until time `nxt`, which might be in the past or future. @@ -118,7 +162,7 @@ void integrate_until(adex_lowered_cell& cell, const time_type end, const time_ty auto delta = end - cur; // membrane potential deviation from resting value auto dE = cell.V_m - cell.E_L; - // leak current + // leak current auto il = cell.g*dE; // spike current auto is = cell.g*cell.delta*exp((cell.V_m - cell.V_th)/cell.delta); diff --git a/arbor/adex_cell_group.hpp b/arbor/adex_cell_group.hpp index c0b509acc..68df87352 100644 --- a/arbor/adex_cell_group.hpp +++ b/arbor/adex_cell_group.hpp @@ -89,6 +89,8 @@ struct ARB_ARBOR_API adex_cell_group: public cell_group { static bool backend_supported(backend_kind kind) { return kind == backend_kind::multicore; } + void edit_cell(cell_gid_type gid, std::any edit) override; + private: enum class adex_probe_kind { voltage, adaption }; diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp index a40108b53..4069a5991 100644 --- a/arbor/backends/gpu/shared_state.cpp +++ b/arbor/backends/gpu/shared_state.cpp @@ -381,6 +381,24 @@ std::pair shared_state::voltage_bounds() const { return minmax_value_impl(n_cv, voltage.data()); } +void shared_state::update_range_parameter(unsigned mid, + arb_index_type lid, + arb_mechanism_ppack& ppack, + cell_gid_type pid, + const std::vector& vals) { + auto off = 0; + auto& store = storage[mid]; + auto& param_d = store.parameters_[pid]; + auto param_h = memory::on_host(memory::device_view(param_d, ppack.width)); + for (auto idx = 0ul; idx < ppack.width; ++idx) { + if (lid != ppack.vec_ci[ppack.node_index[idx]]) continue; + param_h[idx] = vals[off]; + ++off; + } + // conservative copy of the _full_ array... + memory::copy(param_h.data(), param_d); +} + void shared_state::take_samples() { sample_events.mark(); if (!sample_events.empty()) { diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp index 56237e00e..998a621ae 100644 --- a/arbor/backends/gpu/shared_state.hpp +++ b/arbor/backends/gpu/shared_state.hpp @@ -233,6 +233,12 @@ struct ARB_ARBOR_API shared_state: shared_state_base cell_id + cell_gid_type pid, // index of parameter + const std::vector& val); // values we are writing }; // For debugging only diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp index a624cbd3b..e5e2bce4e 100644 --- a/arbor/backends/multicore/shared_state.cpp +++ b/arbor/backends/multicore/shared_state.cpp @@ -23,7 +23,6 @@ #include "multicore_common.hpp" #include "shared_state.hpp" -#include "fvm.hpp" namespace arb { namespace multicore { @@ -183,7 +182,6 @@ void istim_state::add_current(const arb_value_type time, array& current_density) } // shared_state methods: - shared_state::shared_state(task_system_handle, // ignored in mc backend arb_size_type n_cell, arb_size_type n_cv_, @@ -403,10 +401,8 @@ unsigned shared_state::instantiate(arb::mechanism& m, bool peer_indices = !pos_data.peer_cv.empty(); // store indices for random number generation - if (m.mech_.n_random_variables) { - store.gid_ = pos_data.gid; - store.idx_ = pos_data.idx; - } + store.gid_ = pos_data.gid; + if (m.mech_.n_random_variables) store.idx_ = pos_data.idx; // Allocate view pointers (except globals!) store.state_vars_.resize(m.mech_.n_state_vars); m.ppack_.state_vars = store.state_vars_.data(); @@ -537,5 +533,19 @@ unsigned shared_state::instantiate(arb::mechanism& m, return id; } +void shared_state::update_range_parameter(unsigned mid, + arb_index_type lid, + arb_mechanism_ppack& ppack, + cell_gid_type pid, + const std::vector& vals) { + auto off = 0; + auto param = storage[mid].parameters_[pid]; + for (auto idx = 0ul; idx < ppack.width; ++idx) { + if (lid != ppack.vec_ci[ppack.node_index[idx]]) continue; + param[idx] = vals[off]; + ++off; + } +} + } // namespace multicore } // namespace arb diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp index 9adf99937..da609e96f 100644 --- a/arbor/backends/multicore/shared_state.hpp +++ b/arbor/backends/multicore/shared_state.hpp @@ -238,6 +238,14 @@ struct ARB_ARBOR_API shared_state: sample_time_host = util::range_pointer_view(sample_time); sample_value_host = util::range_pointer_view(sample_value); } + + bool mechanism_matches(); + + void update_range_parameter(unsigned mid, // mechanism id + arb_index_type cell_id, // to filter CVs of the cell we are editing + arb_mechanism_ppack& ppack, // ppack, used to read from CV -> cell_id + cell_gid_type pid, // index of parameter + const std::vector& val); // values we are writing }; // For debugging only: diff --git a/arbor/backends/shared_state_base.hpp b/arbor/backends/shared_state_base.hpp index e40f628c5..fbae163ec 100644 --- a/arbor/backends/shared_state_base.hpp +++ b/arbor/backends/shared_state_base.hpp @@ -42,7 +42,7 @@ struct shared_state_base { // samples auto n_samples = util::sum_by(samples, [] (const auto& s) {return s.size();}); if (d->sample_time.size() < n_samples) { - d->sample_time = array(n_samples); + d->sample_time = array(n_samples); d->sample_value = array(n_samples); } initialize(samples, d->sample_events); @@ -80,6 +80,12 @@ struct shared_state_base { } } + // overwrite a RANGE-type parameter in the mechanism `ppack` is attached to. + void update_range_parameter(cell_gid_type lid, arb_mechanism_ppack& ppack, cell_gid_type pid, const std::vector& vals) { + auto d = static_cast(this); + d->update_range_parameter(lid, ppack, pid, vals); + } + arb_value_type* mechanism_state_data(const mechanism& m, const std::string& key) { auto d = static_cast(this); diff --git a/arbor/benchmark_cell_group.cpp b/arbor/benchmark_cell_group.cpp index b1976250e..4222dd67c 100644 --- a/arbor/benchmark_cell_group.cpp +++ b/arbor/benchmark_cell_group.cpp @@ -11,6 +11,7 @@ #include "profile/profiler_macro.hpp" #include "util/span.hpp" +#include "util/maputil.hpp" template void serialize(arb::serializer& s, const K& k, const arb::benchmark_cell_group&); @@ -23,8 +24,7 @@ benchmark_cell_group::benchmark_cell_group(const std::vector& gid const recipe& rec, cell_label_range& cg_sources, cell_label_range& cg_targets): - gids_(gids) -{ + gids_(gids) { for (auto gid: gids_) { if (!rec.get_probes(gid).empty()) { throw bad_cell_probe(cell_kind::benchmark, gid); @@ -47,28 +47,39 @@ benchmark_cell_group::benchmark_cell_group(const std::vector& gid } void benchmark_cell_group::reset() { - for (auto& c: cells_) { - c.time_sequence.reset(); - } - + for (auto& c: cells_) c.time_sequence.reset(); clear_spikes(); } -void benchmark_cell_group::t_serialize(serializer& ser, const std::string& k) const { - serialize(ser, k, *this); -} -void benchmark_cell_group::t_deserialize(serializer& ser, const std::string& k) { - deserialize(ser, k, *this); +void +benchmark_cell_group::edit_cell(cell_gid_type gid, std::any cell_edit) { + try { + auto bench_edit = std::any_cast(cell_edit); + auto lid = util::binary_search_index(gids_, gid); + if (!lid) throw arb::arbor_internal_error{"gid " + std::to_string(gid) + " erroneuosly dispatched to cell group."}; + benchmark_cell& lowered = cells_[*lid]; + auto tmp = benchmark_cell{.source=lowered.source, .target=lowered.target, .time_sequence=std::move(lowered.time_sequence), .realtime_ratio=lowered.realtime_ratio}; + bench_edit(tmp); + if (tmp.source != lowered.source) throw bad_cell_edit(gid, "Source is not editable."); + if (tmp.target != lowered.target) throw bad_cell_edit(gid, "Target is not editable."); + // Write back + lowered.time_sequence = std::move(tmp.time_sequence); + lowered.realtime_ratio = tmp.realtime_ratio; + } + catch (const std::bad_any_cast&) { + throw bad_cell_edit(gid, "Not a Benchmark editor (C++ type-id: '" + std::string{cell_edit.type().name()} + "')"); + } } -cell_kind benchmark_cell_group::get_cell_kind() const { - return cell_kind::benchmark; -} +void benchmark_cell_group::t_serialize(serializer& ser, const std::string& k) const { serialize(ser, k, *this); } + +void benchmark_cell_group::t_deserialize(serializer& ser, const std::string& k) { deserialize(ser, k, *this); } + +cell_kind benchmark_cell_group::get_cell_kind() const { return cell_kind::benchmark; } void benchmark_cell_group::advance(epoch ep, time_type dt, - const event_lane_subrange& event_lanes) -{ + const event_lane_subrange& event_lanes) { using std::chrono::high_resolution_clock; using duration_type = std::chrono::duration; @@ -97,13 +108,9 @@ void benchmark_cell_group::advance(epoch ep, PL(cell); }; -const std::vector& benchmark_cell_group::spikes() const { - return spikes_; -} +const std::vector& benchmark_cell_group::spikes() const { return spikes_; } -void benchmark_cell_group::clear_spikes() { - spikes_.clear(); -} +void benchmark_cell_group::clear_spikes() { spikes_.clear(); } void benchmark_cell_group::add_sampler(sampler_association_handle h, cell_member_predicate probeset_ids, diff --git a/arbor/benchmark_cell_group.hpp b/arbor/benchmark_cell_group.hpp index fa46a8a5c..ad949672f 100644 --- a/arbor/benchmark_cell_group.hpp +++ b/arbor/benchmark_cell_group.hpp @@ -32,6 +32,8 @@ class benchmark_cell_group: public cell_group { void remove_all_samplers() override {} + void edit_cell(cell_gid_type gid, std::any edit) override; + ARB_SERDES_ENABLE(benchmark_cell_group, cells_, spikes_, gids_); void t_serialize(serializer& ser, const std::string& k) const override; diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp index 77afa11de..274d5ec72 100644 --- a/arbor/cable_cell.cpp +++ b/arbor/cable_cell.cpp @@ -81,6 +81,8 @@ struct cable_cell_impl { using index_type = cable_cell::index_type; using size_type = cable_cell::size_type; + bool is_mut_ = false; + // The label dictionary. label_dict dictionary; @@ -117,7 +119,8 @@ struct cable_cell_impl { // Discretization std::optional discretization_; - cable_cell_impl(const arb::morphology& m, const label_dict& labels, const decor& decorations, const std::optional& cvp): + cable_cell_impl(const arb::morphology& m, const label_dict& labels, const decor& decorations, const std::optional& cvp, cable_cell_mutability is_mut=cable_cell_mutability::disabled): + is_mut_(is_mut == cable_cell_mutability::enabled), dictionary(labels), provider(m, dictionary), decorations(decorations), @@ -239,8 +242,12 @@ void cable_cell_impl::init() { } } -cable_cell::cable_cell(const arb::morphology& m, const decor& decorations, const label_dict& dictionary, const std::optional& cvp): - impl_(make_impl(new cable_cell_impl(m, dictionary, decorations, cvp))) +cable_cell::cable_cell(const arb::morphology& m, + const decor& decorations, + const label_dict& dictionary, + const std::optional& cvp, + cable_cell_mutability is_mut): + impl_(make_impl(new cable_cell_impl(m, dictionary, decorations, cvp, is_mut))) {} cable_cell::cable_cell(): impl_(make_impl(new cable_cell_impl())) {} @@ -249,6 +256,8 @@ cable_cell::cable_cell(const cable_cell& other): impl_(make_impl(new cable_cell_impl(*other.impl_))) {} +bool cable_cell::is_editable() const { return impl_->is_mut_; } + const label_dict& cable_cell::labels() const { return impl_->dictionary; } const concrete_embedding& cable_cell::embedding() const { return impl_->provider.embedding(); } const arb::morphology& cable_cell::morphology() const { return impl_->provider.morphology(); } diff --git a/arbor/cable_cell_group.cpp b/arbor/cable_cell_group.cpp index 3cd501894..7c5a0902e 100644 --- a/arbor/cable_cell_group.cpp +++ b/arbor/cable_cell_group.cpp @@ -18,6 +18,7 @@ #include "util/partition.hpp" #include "util/range.hpp" #include "util/span.hpp" +#include "util/maputil.hpp" namespace arb { @@ -491,4 +492,22 @@ std::vector cable_cell_group::get_probe_metadata(const cell_addr } return result; } + +void +cable_cell_group::edit_cell(cell_gid_type gid, std::any cell_edit) { + auto lid = util::binary_search_index(gids_, gid); + if (!lid) throw arb::arbor_internal_error{"gid " + std::to_string(gid) + " erroneuosly dispatched to cell group."}; + try { + auto cc_edit = std::any_cast(cell_edit); + lowered_->edit_cell(gid, *lid, cc_edit); + } + catch (std::bad_any_cast& ex) { + throw bad_cell_edit(gid, "Not a Cable Cell editor (C++ type-id: '" + + std::string{cell_edit.type().name()} + + " ./. " + + std::string{typeid(cable_cell_editor).name()} + + "')"); + } +} + } // namespace arb diff --git a/arbor/cable_cell_group.hpp b/arbor/cable_cell_group.hpp index 4c9f981f3..ba174811f 100644 --- a/arbor/cable_cell_group.hpp +++ b/arbor/cable_cell_group.hpp @@ -44,6 +44,8 @@ struct ARB_ARBOR_API cable_cell_group: public cell_group { void remove_all_samplers() override; + void edit_cell(cell_gid_type gid, std::any edit) override; + std::vector get_probe_metadata(const cell_address_type&) const override; ARB_SERDES_ENABLE(cable_cell_group, gids_, spikes_, lowered_); @@ -52,10 +54,11 @@ struct ARB_ARBOR_API cable_cell_group: public cell_group { void t_deserialize(serializer& ser, const std::string& k) override; static bool backend_supported(backend_kind kind) { return kind == backend_kind::multicore || kind == backend_kind::gpu; } -private: + private: // List of the gids of the cells in the group. std::vector gids_; + // The lowered cell state (e.g. FVM) of the cell. fvm_lowered_cell_ptr lowered_; diff --git a/arbor/cell_group.hpp b/arbor/cell_group.hpp index f7689dd97..6ac757a83 100644 --- a/arbor/cell_group.hpp +++ b/arbor/cell_group.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include "epoch.hpp" #include "event_lane.hpp" @@ -20,8 +21,7 @@ // ranges are needed to map (gid, label) pairs to their corresponding lid sets. namespace arb { -class cell_group { -public: +struct cell_group { virtual ~cell_group() = default; virtual cell_kind get_cell_kind() const = 0; @@ -32,20 +32,23 @@ class cell_group { virtual const std::vector& spikes() const = 0; virtual void clear_spikes() = 0; - // Sampler association methods below should be thread-safe, as they might be invoked - // from a sampler call back called from a different cell group running on a different thread. - + // Sampler association methods below should be thread-safe, as they might be + // invoked from a sampler call back called from a different cell group + // running on a different thread. virtual void add_sampler(sampler_association_handle, cell_member_predicate, schedule, sampler_function) = 0; virtual void remove_sampler(sampler_association_handle) = 0; virtual void remove_all_samplers() = 0; - // Probe metadata queries might also be called while a simulation is running, and so should - // also be thread-safe. + // allow editing of certain cell properties + virtual void edit_cell(cell_gid_type gid, std::any edit) = 0; + // Probe metadata queries might also be called while a simulation is + // running, and so should also be thread-safe. virtual std::vector get_probe_metadata(const cell_address_type&) const { return {}; } + // trampolines for serialization virtual void t_serialize(serializer& s, const std::string&) const = 0; - virtual void t_deserialize(serializer& s, const std::string&) = 0; + virtual void t_deserialize(serializer& s, const std::string&) = 0; }; using cell_group_ptr = std::unique_ptr; diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index 3d2f9813c..d38aa1811 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -455,12 +455,12 @@ fvm_cv_discretize(const cable_cell& cell, } ARB_ARBOR_API fvm_cv_discretization fvm_cv_discretize(const std::vector& cells, - const cable_cell_parameter_set& global_defaults, - const arb::execution_context& ctx) -{ + const cable_cell_parameter_set& global_defaults, + const arb::execution_context& ctx) { std::vector cell_disc(cells.size()); - threading::parallel_for::apply(0, cells.size(), ctx.thread_pool.get(), - [&] (int i) { cell_disc[i]=fvm_cv_discretize(cells[i], global_defaults);}); + threading::parallel_for::apply(0, cells.size(), + ctx.thread_pool.get(), + [&] (int i) { cell_disc[i] = fvm_cv_discretize(cells[i], global_defaults);}); fvm_cv_discretization combined; for (auto cell_idx: count_along(cells)) { @@ -927,7 +927,7 @@ make_revpot_mechanism_config(const std::unordered_map& gj_conns, diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp index 9b8ec3bd7..61ed211d5 100644 --- a/arbor/fvm_layout.hpp +++ b/arbor/fvm_layout.hpp @@ -320,6 +320,13 @@ struct fvm_mechanism_data { bool post_events = false; }; +ARB_ARBOR_API fvm_mechanism_data +fvm_build_mechanism_data(const cable_cell_global_properties& gprop, + const cable_cell& cell, + const std::vector& gj_conns, + const fvm_cv_discretization& D, + arb_size_type cell_idx); + ARB_ARBOR_API fvm_mechanism_data fvm_build_mechanism_data( const cable_cell_global_properties& gprop, const std::vector& cells, diff --git a/arbor/fvm_lowered_cell.hpp b/arbor/fvm_lowered_cell.hpp index 050713b11..8d0544c46 100644 --- a/arbor/fvm_lowered_cell.hpp +++ b/arbor/fvm_lowered_cell.hpp @@ -230,8 +230,22 @@ struct fvm_initialization_data { } }; -// Common base class for FVM implementation on host or gpu back-end. +struct fvm_mutable_data { + // densities :: (where, derived-mech-name) -> parameters + std::vector> density_overrides_; + // synapses :: (where, derived-mech-name, tag) -> parameters + std::vector> synapse_overrides_; + // morphology and CVP (needed for re-discretisation) + morphology morph; + std::optional cvp; + decor dec; + label_dict lbl; + // local cell index + arb_index_type lid; +}; + +// Common base class for FVM implementation on host or gpu back-end. struct fvm_lowered_cell { virtual void reset() = 0; @@ -242,18 +256,20 @@ struct fvm_lowered_cell { const event_lane_subrange& event_lanes, const std::vector>& staged_samples) = 0; + virtual void edit_cell(cell_gid_type gid, cell_gid_type lid, const cable_cell_editor& edit) = 0; + virtual arb_value_type time() const = 0; virtual ~fvm_lowered_cell() {} virtual void t_serialize(serializer& ser, const std::string& k) const = 0; virtual void t_deserialize(serializer& ser, const std::string& k) = 0; + }; using fvm_lowered_cell_ptr = std::unique_ptr; -ARB_ARBOR_API fvm_lowered_cell_ptr make_fvm_lowered_cell(backend_kind p, const execution_context& ctx, - std::uint64_t seed = 0); +ARB_ARBOR_API fvm_lowered_cell_ptr make_fvm_lowered_cell(backend_kind p, const execution_context& ctx, std::uint64_t seed = 0); inline void serialize(serializer& s, const std::string& k, const fvm_lowered_cell& v) { v.t_serialize(s, k); } diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp index e7fa2ccb3..a5be20aac 100644 --- a/arbor/fvm_lowered_cell_impl.hpp +++ b/arbor/fvm_lowered_cell_impl.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -32,6 +33,7 @@ #include "util/transform.hpp" namespace arb { + template struct fvm_lowered_cell_impl: public fvm_lowered_cell { using backend = Backend; @@ -55,14 +57,13 @@ struct fvm_lowered_cell_impl: public fvm_lowered_cell { value_type time() const override { return state_->time; } - //Exposed for testing purposes - std::vector& mechanisms() { return mechanisms_; } - ARB_SERDES_ENABLE(fvm_lowered_cell_impl, seed_, state_); void t_serialize(serializer& ser, const std::string& k) const override { serialize(ser, k, *this); } void t_deserialize(serializer& ser, const std::string& k) override { deserialize(ser, k, *this); } + void edit_cell(cell_gid_type gid, cell_gid_type lid, const cable_cell_editor& edit) override; + // Host or GPU-side back-end dependent storage. using array = typename backend::array; using shared_state = typename backend::shared_state; @@ -71,10 +72,17 @@ struct fvm_lowered_cell_impl: public fvm_lowered_cell { std::unique_ptr state_; // Cell state shared across mechanisms. + std::unordered_map mechptr_by_name_; + std::vector mechanisms_; // excludes reversal potential calculators. std::vector revpot_mechanisms_; std::vector voltage_mechanisms_; + // mutable cells. + std::unordered_map mutable_cell_data_; + // _needed_ for mutable cells + cable_cell_global_properties gprop; + // Handles for accessing event targets. std::vector target_handles_; // Lookup table for target ids -> local target handle indices. @@ -116,7 +124,6 @@ struct fvm_lowered_cell_impl: public fvm_lowered_cell { 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); @@ -297,7 +304,6 @@ fvm_lowered_cell_impl::add_probes(const std::vector& gid 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) { @@ -308,7 +314,7 @@ fvm_lowered_cell_impl::add_probes(const std::vector& gid cell_gid_type gid = gids[cell_idx]; const auto& rec_probes = rec.get_probes(gid); for (const auto& pi: rec_probes) { - resolve_probe_address(probe_data, cells, cell_idx, pi.address, D, mech_data, target_handles, mechptr_by_name); + resolve_probe_address(probe_data, cells, cell_idx, pi.address, D, mech_data, target_handles, mechptr_by_name_); 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); @@ -318,6 +324,72 @@ fvm_lowered_cell_impl::add_probes(const std::vector& gid } } + +template void +fvm_lowered_cell_impl::edit_cell(cell_gid_type gid, + cell_gid_type lid, + const cable_cell_editor& edit) { + if (!mutable_cell_data_.contains(gid)) throw bad_cell_edit{gid, "Cable cell is not editable. Did you enable editing at creation?"}; + auto& mut_data = mutable_cell_data_[gid]; + auto new_dec = decor{}; + // Handle density changes + if (auto fun = edit.on_density) { + for (auto& [reg, mech, params]: mut_data.density_overrides_) { + auto new_params = fun(reg, mech, params); + for (const auto& [nk, nv]: new_params) { + bool found = false; + for (auto& [ok, ov]: params) { + if (ok == nk) { + ov = nv; + found = true; + break; + } + } + if (!found) throw bad_cell_edit{gid, "Unknown parameter '" + std::string{nk} + "' for mechanism '" + mech + "'." }; + } + // found all new values. add mech to decor + new_dec.paint(reg, density(mech.c_str(), params)); + } + } + if (auto fun = edit.on_synapse) { + for (auto& [loc, mech, tag, params]: mut_data.synapse_overrides_) { + auto new_params = fun(loc, mech, params); + for (const auto& [nk, nv]: new_params) { + // TODO: Using STL patterns here trips up older compilers + // due to capturing from a structured binding. + bool found = false; + for (auto& [ok, ov]: params) { + if (ok == nk) { + ov = nv; + found = true; + break; + } + } + if (!found) throw bad_cell_edit{gid, "Unknown parameter '" + std::string{nk} + "' for mechanism '" + mech + "'." }; + } + // found all new values. add mech to decor + new_dec.place(loc, synapse(mech.c_str(), params), tag); + } + } + + auto new_cell = cable_cell{mut_data.morph, new_dec, mut_data.lbl, mut_data.cvp}; + + auto D = fvm_cv_discretize(new_cell, gprop.default_parameters); + // TODO do we need the (correct) gap junctions here? Only if we do GJs? + auto M = fvm_build_mechanism_data(gprop, new_cell, {}, D, 0); + for (const auto& [mech, data]: M.mechanisms) { + auto& ptr = mechptr_by_name_.at(mech); + auto& type = ptr->mech_; + for (auto off: util::make_span(type.n_parameters)) { + auto pname = type.parameters[off].name; + auto pset = std::find_if(data.param_values.begin(), data.param_values.end(), + [&pname] (const auto& it) { return it.first == pname; }); + if (pset == data.param_values.end()) arbor_internal_error{"Cannot find parameter by name? Expected " + std::string{pname}}; + state_->update_range_parameter(ptr->mechanism_id(), mut_data.lid, ptr->ppack_, off, pset->second); + } + } +} + inline auto get_cable_cell_global_properties(const recipe& rec) { try { @@ -378,6 +450,9 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid } fvm_info.shrink_to_fit(); + // TODO think about (how) making this optional. + gprop = std::any_cast(rec.get_global_properties(cell_kind::cable)); + // extract and verify global settings auto global_props = get_cable_cell_global_properties(rec); // fetch backend specific mechanism data @@ -400,7 +475,7 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid // Fill src_to_spike and cv_to_cell vectors only if mechanisms with post_events implemented are present. post_events_ = mech_data.post_events; auto max_detector = 0; - std::vector src_to_spike, cv_to_cell; + std::vector src_to_spike; if (post_events_) { max_detector = util::max_element_by(fvm_info.num_sources, [](const auto& elem) { return util::second(elem); })->second; for (auto cell_idx: util::make_span(ncell)) { @@ -409,9 +484,10 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid } } src_to_spike.shrink_to_fit(); - cv_to_cell = D.geometry.cv_to_cell; } + auto cv_to_cell = D.geometry.cv_to_cell; + // map control volume ids to global cell ids std::vector cv_to_gid(D.geometry.cv_to_cell.size()); std::transform(D.geometry.cv_to_cell.begin(), D.geometry.cv_to_cell.end(), @@ -437,8 +513,8 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid data_alignment, seed_); + target_handles_.resize(mech_data.n_target); // Keep track of mechanisms by name for probe lookup. - std::unordered_map mechptr_by_name; target_handles_.resize(mech_data.n_target); for (const auto& [name, config]: mech_data.mechanisms) { auto n_cv = config.cv.size(); @@ -486,7 +562,7 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid } } } - mechptr_by_name[name] = mech.get(); + mechptr_by_name_[name] = mech.get(); mechanisms_.emplace_back(mech.release()); break; } @@ -499,7 +575,7 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid } auto [mech, over] = mech_instance(name); state_->instantiate(*mech, over, layout, config.param_values); - mechptr_by_name[name] = mech.get(); + mechptr_by_name_[name] = mech.get(); mechanisms_.emplace_back(mech.release()); break; } @@ -515,7 +591,7 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid } auto [mech, over] = mech_instance(name); state_->instantiate(*mech, over, layout, config.param_values); - mechptr_by_name[name] = mech.get(); + mechptr_by_name_[name] = mech.get(); voltage_mechanisms_.emplace_back(mech.release()); break; } @@ -531,7 +607,7 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid } auto [mech, over] = mech_instance(name); state_->instantiate(*mech, over, layout, config.param_values); - mechptr_by_name[name] = mech.get(); + mechptr_by_name_[name] = mech.get(); mechanisms_.emplace_back(mech.release()); break; } @@ -540,7 +616,7 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid // to any currents, so leave weights as zero. auto [mech, over] = mech_instance(name); state_->instantiate(*mech, over, layout, config.param_values); - mechptr_by_name[name] = mech.get(); + mechptr_by_name_[name] = mech.get(); revpot_mechanisms_.emplace_back(mech.release()); break; } @@ -548,14 +624,73 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid } } - add_probes(gids, cells, rec, D, mechptr_by_name, mech_data, target_handles_, fvm_info.probe_map); + // prepare data for editable cells + for (auto ix: util::make_span(ncell)) { + const auto& cell = cells[ix]; + const auto gid = gids[ix]; + if (cell.is_editable()) { + mutable_cell_data_[gid] = fvm_mutable_data { + .morph = cell.morphology(), + .cvp = cell.discretization(), + .lbl = cell.labels(), + .lid = static_cast(ix), + }; + + const auto& dec = cell.decorations(); + + for (const auto& [reg, it]: dec.paintings()) { + if (std::holds_alternative(it)) { + const auto& mech = std::get(it); + const auto& name = mech.mech.name(); + const auto& ptr = mechptr_by_name_.at(name); + const auto& type = ptr->mech_; + parameter_map params; + for (arb_size_type ix = 0; ix < type.n_parameters; ++ix) { + const auto& param = type.parameters[ix]; + const auto& key = param.name; + arb_value_type val = param.default_value; + for (const auto& [k, v]: mech.mech.values()) { + if (k == param.name) val = v; + } + params.emplace_back(key, val); + } + mutable_cell_data_[gid].density_overrides_.emplace_back(reg, name, params); + } + } + + for (const auto& [loc, it, hash]: dec.placements()) { + if (std::holds_alternative(it)) { + const auto& mech = std::get(it); + const auto& name = mech.mech.name(); + const auto& ptr = mechptr_by_name_.at(name); + const auto& type = ptr->mech_; + const auto& tag = dec.tag_of(hash); + parameter_map params; + for (arb_size_type ix = 0; ix < type.n_parameters; ++ix) { + const auto& param = type.parameters[ix]; + const auto& key = param.name; + arb_value_type val = param.default_value; + for (const auto& [k, v]: mech.mech.values()) { + if (k == param.name) val = v; + } + params.emplace_back(key, val); + } + mutable_cell_data_[gid].synapse_overrides_.emplace_back(loc, name, tag, params); + } + } + } + } + + + + add_probes(gids, cells, rec, D, mech_data, target_handles_, fvm_info.probe_map); // Create lookup structure for target ids. util::make_partition(target_handle_divisions_, util::transform_view(gids, [&](cell_gid_type i) { return fvm_info.num_targets[i]; })); - + reset(); return fvm_info; } @@ -950,7 +1085,7 @@ void resolve_probe(const cable_probe_point_state_cell& p, probe_resolution_data< ++lid; } - + result.metadata = std::move(metadata); result.shrink_to_fit(); R.result.push_back(std::move(result)); diff --git a/arbor/include/arbor/adex_cell.hpp b/arbor/include/arbor/adex_cell.hpp index 48c9a7602..4893fd2d2 100644 --- a/arbor/include/arbor/adex_cell.hpp +++ b/arbor/include/arbor/adex_cell.hpp @@ -49,4 +49,6 @@ struct ARB_SYMBOL_VISIBLE adex_probe_voltage {}; // Sample value type: `double` struct ARB_SYMBOL_VISIBLE adex_probe_adaption {}; +using adex_cell_editor = std::function; + } // namespace arb diff --git a/arbor/include/arbor/arbexcept.hpp b/arbor/include/arbor/arbexcept.hpp index 22ced920b..d03ff85e1 100644 --- a/arbor/include/arbor/arbexcept.hpp +++ b/arbor/include/arbor/arbexcept.hpp @@ -47,6 +47,12 @@ struct ARB_SYMBOL_VISIBLE resolution_disabled: arbor_exception { {} }; +struct ARB_SYMBOL_VISIBLE bad_cell_edit: arbor_exception { + bad_cell_edit(cell_gid_type gid, std::string why): + arbor_exception("Cannot edit cell gid=" + std::to_string(gid) + ": " + why) + {} + std::string why; +}; struct ARB_SYMBOL_VISIBLE dup_cell_probe: arbor_exception { dup_cell_probe(cell_kind kind, cell_gid_type gid, cell_tag_type tag); diff --git a/arbor/include/arbor/benchmark_cell.hpp b/arbor/include/arbor/benchmark_cell.hpp index f8569fe31..c15ac25f1 100644 --- a/arbor/include/arbor/benchmark_cell.hpp +++ b/arbor/include/arbor/benchmark_cell.hpp @@ -22,15 +22,11 @@ struct ARB_SYMBOL_VISIBLE benchmark_cell { // Time taken in ms to advance the cell one ms of simulation time. // If equal to 1, then a single cell can be advanced in realtime - double realtime_ratio; - - benchmark_cell() = default; - benchmark_cell(cell_tag_type source, cell_tag_type target, schedule seq, double ratio): - source(source), target(target), time_sequence(seq), realtime_ratio(ratio) {}; + double realtime_ratio = 1.0; ARB_SERDES_ENABLE(benchmark_cell, source, target, time_sequence, realtime_ratio); }; -} // namespace arb - +using benchmark_cell_editor = std::function; +} // namespace arb diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp index 4c15c2ee9..a91ecc741 100644 --- a/arbor/include/arbor/cable_cell.hpp +++ b/arbor/include/arbor/cable_cell.hpp @@ -24,6 +24,12 @@ namespace arb { +enum class cable_cell_mutability { + disabled, + enabled, +}; + + // `cable_sample_range` is simply a pair of `const double*` pointers describing the sequence // of double values associated with the cell-wide sample. @@ -259,6 +265,18 @@ using location_assignment = std::conditional_t || std std::unordered_map>, mlocation_map>; +// Allowed edits on cable cells +using parameter_map = std::vector>; + +using density_editor = std::function; +using synapse_editor = std::function; + +// Overwrite a list of named parameters on a given mechanism +struct ARB_SYMBOL_VISIBLE cable_cell_editor { + density_editor on_density; + synapse_editor on_synapse; +}; + // High-level abstract representation of a cell. struct ARB_SYMBOL_VISIBLE cable_cell { using lid_range_map = std::unordered_multimap; @@ -281,7 +299,8 @@ struct ARB_SYMBOL_VISIBLE cable_cell { cable_cell(const class morphology& m, const decor& d, const label_dict& l={}, - const std::optional& = {}); + const std::optional& = {}, + const cable_cell_mutability is_mut=cable_cell_mutability::disabled); /// Access to labels const label_dict& labels() const; @@ -292,10 +311,10 @@ struct ARB_SYMBOL_VISIBLE cable_cell { const mprovider& provider() const; // Convenience access to placed items. - const std::unordered_map>& synapses() const; - const std::unordered_map>& junctions() const; - const mlocation_map& detectors() const; - const mlocation_map& stimuli() const; + const location_assignment& synapses() const; + const location_assignment& junctions() const; + const location_assignment& detectors() const; + const location_assignment& stimuli() const; // Convenience access to painted items. const region_assignment densities() const; @@ -329,6 +348,8 @@ struct ARB_SYMBOL_VISIBLE cable_cell { const lid_range_map& synapse_ranges() const; const lid_range_map& junction_ranges() const; + bool is_editable() const; + private: std::unique_ptr impl_; }; diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp index 44a496e8f..74f4c2df9 100644 --- a/arbor/include/arbor/cable_cell_param.hpp +++ b/arbor/include/arbor/cable_cell_param.hpp @@ -427,9 +427,9 @@ class ARB_ARBOR_API decor { std::unordered_map hashes_; public: - const auto& paintings() const {return paintings_; } - const auto& placements() const {return placements_; } - const auto& defaults() const {return defaults_; } + const auto& paintings() const { return paintings_; } + const auto& placements() const { return placements_; } + const auto& defaults() const { return defaults_; } decor& paint(region, paintable); decor& place(locset, placeable, cell_tag_type); diff --git a/arbor/include/arbor/lif_cell.hpp b/arbor/include/arbor/lif_cell.hpp index 2af3c2bf7..bc29fd945 100644 --- a/arbor/include/arbor/lif_cell.hpp +++ b/arbor/include/arbor/lif_cell.hpp @@ -1,5 +1,7 @@ #pragma once +#include + #include #include #include @@ -32,4 +34,6 @@ struct ARB_SYMBOL_VISIBLE lif_probe_metadata {}; // Sample value type: `double` struct ARB_SYMBOL_VISIBLE lif_probe_voltage {}; +using lif_cell_editor = std::function; + } // namespace arb diff --git a/arbor/include/arbor/mechanism.hpp b/arbor/include/arbor/mechanism.hpp index 691c5a856..2b13b59fc 100644 --- a/arbor/include/arbor/mechanism.hpp +++ b/arbor/include/arbor/mechanism.hpp @@ -95,7 +95,7 @@ class mechanism { // Per-cell group identifier for an instantiated mechanism. unsigned mechanism_id() const { return ppack_.mechanism_id; } - arb_mechanism_type mech_; + arb_mechanism_type mech_; arb_mechanism_interface iface_; arb_mechanism_ppack ppack_; diff --git a/arbor/include/arbor/simulation.hpp b/arbor/include/arbor/simulation.hpp index 824a2e731..3044b13f9 100644 --- a/arbor/include/arbor/simulation.hpp +++ b/arbor/include/arbor/simulation.hpp @@ -3,7 +3,6 @@ #include #include #include -#include #include #include @@ -28,7 +27,7 @@ using balancer_function = std::function seqs; + std::vector schedules; spike_source_cell() = delete; template - spike_source_cell(cell_tag_type source, Seqs&&... seqs): source(std::move(source)), seqs{std::forward(seqs)...} {} - spike_source_cell(cell_tag_type source, std::vector seqs): source(std::move(source)), seqs(std::move(seqs)) {} + spike_source_cell(cell_tag_type source, Seqs&&... seqs): source(std::move(source)), schedules{std::forward(seqs)...} {} + spike_source_cell(cell_tag_type source, std::vector seqs): source(std::move(source)), schedules(std::move(seqs)) {} }; +using spike_source_cell_editor = std::function; + } // namespace arb diff --git a/arbor/lif_cell_group.cpp b/arbor/lif_cell_group.cpp index 4b9dd05b5..113c9ecd3 100644 --- a/arbor/lif_cell_group.cpp +++ b/arbor/lif_cell_group.cpp @@ -4,6 +4,7 @@ #include "lif_cell_group.hpp" #include "profile/profiler_macro.hpp" #include "util/rangeutil.hpp" +#include "util/maputil.hpp" #include "util/span.hpp" using namespace arb; @@ -39,12 +40,47 @@ lif_cell_group::lif_cell_group(const std::vector& gids, } } } + + if (!util::is_sorted(gids_)) throw arb::arbor_internal_error{"gids must be sorted?!"}; } -cell_kind lif_cell_group::get_cell_kind() const { - return cell_kind::lif; +void +lif_cell_group::edit_cell(cell_gid_type gid, std::any cell_edit) { + try { + auto lif_edit = std::any_cast(cell_edit); + auto lid = util::binary_search_index(gids_, gid); + if (!lid) throw arb::arbor_internal_error{"gid " + std::to_string(gid) + " erroneuosly dispatched to cell group."}; + auto& lowered = cells_[*lid]; + auto tmp = lif_cell{ + .source = lowered.source, + .target = lowered.target, + .tau_m = lowered.tau_m * U::ms, + .V_th = lowered.V_th * U::mV, + .C_m = lowered.C_m * U::pF, + .E_L = lowered.E_L * U::mV, + .E_R = lowered.E_R * U::mV, + .V_m = lowered.V_m * U::mV, + .t_ref = lowered.t_ref * U::ms + }; + lif_edit(tmp); + // NOTE: we forbid writing to V_m? Reasons + // * the cell might be in the refractory period which causes semantic issues + // - return to normal or not? + // - what should probes return + // * V_m is the _initial state_ only + if (tmp.V_m.value_as(U::mV) != lowered.V_m) throw bad_cell_edit(gid, "Initial voltage is not editable."); + if (tmp.source != lowered.source) throw bad_cell_edit(gid, "Source is not editable."); + if (tmp.target != lowered.target) throw bad_cell_edit(gid, "Target is not editable."); + // Write back + lowered = lif_lowered_cell{tmp}; + } + catch (const std::bad_any_cast&) { + throw bad_cell_edit(gid, "Not a LIF editor (C++ type-id: '" + std::string{cell_edit.type().name()} + "')"); + } } +cell_kind lif_cell_group::get_cell_kind() const { return cell_kind::lif; } + void lif_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& event_lanes) { PE(lif); for (auto lid: util::make_span(gids_.size())) { diff --git a/arbor/lif_cell_group.hpp b/arbor/lif_cell_group.hpp index de8342755..3448d0db7 100644 --- a/arbor/lif_cell_group.hpp +++ b/arbor/lif_cell_group.hpp @@ -77,6 +77,8 @@ struct ARB_ARBOR_API lif_cell_group: public cell_group { virtual void t_serialize(serializer& ser, const std::string& k) const override; virtual void t_deserialize(serializer& ser, const std::string& k) override; + void edit_cell(cell_gid_type gid, std::any edit) override; + static bool backend_supported(backend_kind kind) { return kind == backend_kind::multicore; } private: diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index 229786052..ef4fa723f 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include #include #include @@ -74,6 +76,8 @@ struct simulation_state { void update(const recipe& rec); + void edit_cell(cell_gid_type gid, std::any edit); + void reset(); time_type run(time_type tfinal, time_type dt); @@ -565,9 +569,7 @@ std::vector simulation_state::get_probe_metadata(const cell_addr if (auto lidx = util::value_by_key(gid_to_cell_index_, probeset_id.gid)) { return cell_groups_.at(*lidx)->get_probe_metadata(probeset_id); } - else { - return {}; - } + return {}; } // Simulation class implementations forward to implementation class. @@ -581,12 +583,19 @@ simulation::simulation(const recipe& rec, impl_.reset(new simulation_state(rec, decomp, ctx, seed)); } -void simulation::reset() { - impl_->reset(); -} +void simulation::reset() { impl_->reset(); } void simulation::update(const recipe& rec) { impl_->update(rec); } +// facilitate cell editig +void simulation::edit_cell(cell_gid_type gid, std::any edit) { impl_->edit_cell(gid, edit); } +void simulation_state::edit_cell(cell_gid_type gid, std::any edit) { + if (gid >= ddc_->num_global_cells()) throw std::range_error{"Not a valid gid: " + std::to_string(gid)}; + if (auto gidx = util::value_by_key(gid_to_cell_index_, gid)) { + cell_groups_[*gidx]->edit_cell(gid, edit); + } +} + time_type simulation::run(const units::quantity& tfinal, const units::quantity& dt) { auto dt_ms = dt.value_as(units::ms); if (dt_ms <= 0.0 || !std::isfinite(dt_ms)) throw domain_error("Finite time-step must be supplied."); diff --git a/arbor/spike_source_cell_group.cpp b/arbor/spike_source_cell_group.cpp index 8c8dc4f97..49dbd4891 100644 --- a/arbor/spike_source_cell_group.cpp +++ b/arbor/spike_source_cell_group.cpp @@ -8,16 +8,15 @@ #include "profile/profiler_macro.hpp" #include "spike_source_cell_group.hpp" #include "util/span.hpp" +#include "util/maputil.hpp" namespace arb { -spike_source_cell_group::spike_source_cell_group( - const std::vector& gids, - const recipe& rec, - cell_label_range& cg_sources, - cell_label_range& cg_targets): - gids_(gids) -{ +spike_source_cell_group::spike_source_cell_group(const std::vector& gids, + const recipe& rec, + cell_label_range& cg_sources, + cell_label_range& cg_targets): + gids_(gids) { for (auto gid: gids_) { if (!rec.get_probes(gid).empty()) { throw bad_cell_probe(cell_kind::spike_source, gid); @@ -30,8 +29,9 @@ spike_source_cell_group::spike_source_cell_group( cg_targets.add_cell(); try { auto cell = util::any_cast(rec.get_cell_description(gid)); - time_sequences_.emplace_back(cell.seqs); + time_sequences_.emplace_back(cell.schedules); cg_sources.add_label(hash_value(cell.source), {0, 1}); + sources_.push_back(cell.source); } catch (std::bad_any_cast& e) { throw bad_cell_description(cell_kind::spike_source, gid); @@ -39,9 +39,7 @@ spike_source_cell_group::spike_source_cell_group( } } -cell_kind spike_source_cell_group::get_cell_kind() const { - return cell_kind::spike_source; -} +cell_kind spike_source_cell_group::get_cell_kind() const { return cell_kind::spike_source; } void spike_source_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& event_lanes) { PE(sscell); @@ -68,22 +66,36 @@ void spike_source_cell_group::reset() { clear_spikes(); } -const std::vector& spike_source_cell_group::spikes() const { - return spikes_; +void +spike_source_cell_group::edit_cell(cell_gid_type gid, std::any cell_edit) { + try { + auto source_edit = std::any_cast(cell_edit); + auto lid = util::binary_search_index(gids_, gid); + if (!lid) throw arb::arbor_internal_error{"gid " + std::to_string(gid) + " erroneuosly dispatched to cell group."}; + auto idx = *lid; + auto tmp = spike_source_cell{sources_[idx], std::move(time_sequences_[idx])}; + source_edit(tmp); + // NOTE: we forbid writing to V_m? Reasons + // * the cell might be in the refractory period which causes semantic issues + // - return to normal or not? + // - what should probes return + // * V_m is the _initial state_ only + if (tmp.source != sources_[idx]) throw bad_cell_edit(gid, "Source is not editable."); + // Write back + time_sequences_[idx] = std::move(tmp.schedules); + } + catch (const std::bad_any_cast&) { + throw bad_cell_edit(gid, "Not a source cell editor (C++ typid: '" + std::string{cell_edit.type().name()} + "')"); + } } -void spike_source_cell_group::clear_spikes() { - spikes_.clear(); -} +const std::vector& spike_source_cell_group::spikes() const { return spikes_; } -void spike_source_cell_group::t_serialize(serializer& ser, const std::string& k) const { - serialize(ser, k, *this); -} +void spike_source_cell_group::clear_spikes() { spikes_.clear(); } -void spike_source_cell_group::t_deserialize(serializer& ser, const std::string& k) { - deserialize(ser, k, *this); -} +void spike_source_cell_group::t_serialize(serializer& ser, const std::string& k) const { serialize(ser, k, *this); } +void spike_source_cell_group::t_deserialize(serializer& ser, const std::string& k) { deserialize(ser, k, *this); } void spike_source_cell_group::add_sampler(sampler_association_handle, cell_member_predicate, schedule, sampler_function) {} diff --git a/arbor/spike_source_cell_group.hpp b/arbor/spike_source_cell_group.hpp index 148fd2238..6446f70e3 100644 --- a/arbor/spike_source_cell_group.hpp +++ b/arbor/spike_source_cell_group.hpp @@ -35,14 +35,17 @@ class ARB_ARBOR_API spike_source_cell_group: public cell_group { void remove_all_samplers() override {} - ARB_SERDES_ENABLE(spike_source_cell_group, spikes_, gids_, time_sequences_); + ARB_SERDES_ENABLE(spike_source_cell_group, sources_, spikes_, gids_, time_sequences_); - virtual void t_serialize(serializer& ser, const std::string& k) const override; - virtual void t_deserialize(serializer& ser, const std::string& k) override; + void t_serialize(serializer& ser, const std::string& k) const override; + void t_deserialize(serializer& ser, const std::string& k) override; + + void edit_cell(cell_gid_type gid, std::any edit) override; static bool backend_supported(backend_kind kind) { return kind == backend_kind::multicore; } private: + std::vector sources_; std::vector spikes_; std::vector gids_; std::vector> time_sequences_; diff --git a/doc/concepts/benchmark_cell.rst b/doc/concepts/benchmark_cell.rst index d443e180e..1b84170ea 100644 --- a/doc/concepts/benchmark_cell.rst +++ b/doc/concepts/benchmark_cell.rst @@ -3,10 +3,12 @@ Benchmark cells =============== -The description of a benchmark cell is used to determine the spiking schedule of the cell and manipulate its -performance efficiency. This cell is mainly used by developers. +The description of a benchmark cell is used to determine the spiking schedule of +the cell and manipulate its performance efficiency. This cell is mainly used by +developers. API --- +* :ref:`C++ ` * :ref:`Python ` diff --git a/doc/concepts/cable_cell.rst b/doc/concepts/cable_cell.rst index 9a0b6e225..01b33b85b 100644 --- a/doc/concepts/cable_cell.rst +++ b/doc/concepts/cable_cell.rst @@ -49,6 +49,25 @@ Once constructed, the cable cell can be queried for specific information about t * 3 decors (1 for each of purkinje, granule and pyramidal). * 1 label dictionary that defines the region types. +Certain parameters on cable cells can be updated ('edited') in between calls to +``simulation::run``. Currently, this covers + +* Density mechanism paramters (*not* ``iexpr``) +* Point mechanism paramters (*not* ``iexpr``) + +and needs to be enabled per cell. + +.. Note:: + + 1. We do offer this functionality, however, it comes at comes at a cost. + Largely, this concerns runtime performance and memory usage as storing and + updating the information incurs non-trivial bookkeeping cost. If you want + to adapt ion channels over time, consider adding a ``time`` parameter to + the appropriate ``.mod`` files and dispatch behaviour based on that. + 2. This is an experimental feature. More types of parameters may be covered + in the future, please feel encouraged to make requests accordingly to your + needs. + .. toctree:: :maxdepth: 2 diff --git a/doc/concepts/simulation.rst b/doc/concepts/simulation.rst index dcfaeab9a..dbf8e711e 100644 --- a/doc/concepts/simulation.rst +++ b/doc/concepts/simulation.rst @@ -33,11 +33,18 @@ Simulation execution and interaction Simulations provide an interface for executing and interacting with the model: -* The simulation is executed/*run* by advancing the model state from the current simulation time to another - with maximum time step size. -* The model state can be *reset* to its initial state before the simulation was started. -* *Sampling* of the simulation state can be performed during execution with samplers and probes - and spike output with the total number of spikes generated since either construction or reset. +* The simulation is executed/*run* by advancing the model state from the current + simulation time to another with maximum time step size. +* The model state can be *reset* to its initial state before the simulation was + started. +* *Sampling* of the simulation state can be performed during execution with + samplers and probes and spike output with the total number of spikes generated + since either construction or reset. +* *Updating the (biological) network*, whenever control is passed to the user, + i.e. when the current execution step is finished, the connections between cells + can be altered. +* *Tweaking cell settings*, similarly, at synchronisation points, some cell + parameters can be changed. API --- diff --git a/doc/cpp/adex_cell.rst b/doc/cpp/adex_cell.rst index e43681ff7..51e53ab82 100644 --- a/doc/cpp/adex_cell.rst +++ b/doc/cpp/adex_cell.rst @@ -70,3 +70,10 @@ AdEx cells .. cpp:member:: double b Increase in :math:`w` after emitted spike :math:`b = 0.08\,nA` [nA] + +.. cpp:type:: adex_cell_editor = std::function + + Callback function to update setting of AdEx cells in place via the + ``simulation::edit_cell`` interace. All values are changeable, except: + ``source``, ``target``, ``w``, and ``V_m``. + diff --git a/doc/cpp/benchmark_cell.rst b/doc/cpp/benchmark_cell.rst new file mode 100644 index 000000000..c843022c2 --- /dev/null +++ b/doc/cpp/benchmark_cell.rst @@ -0,0 +1,51 @@ +.. _cppbenchcell: + +Benchmark cells +=============== + +.. cpp:namespace:: arb + +.. cpp:class:: lif_cell + + A benchmarking cell, used by Arbor developers to test communication performance. + + .. cpp:function:: benchmark_cell(const cell_tag_type& source, const cell_tag_type& target, schedule, double realtime_ratio) + + Construct a benchmark cell with a single built-in source with label + ``source``; and a single built-in target with label ``target``. The + labels can be used for forming connections from/to the cell in the + :cpp:class:`arb::recipe` by creating a :cpp:class:`arb::connection`. + + A benchmark cell generates spikes at a user-defined sequence of time points: + + - at regular intervals (using an :cpp:class:`arb::regular_schedule`) + - at a sequence of user-defined times (using an :cpp:class:`arb::explicit_schedule`) + - at times defined by a Poisson sequence (using an :cpp:class:`arb::poisson_schedule`) + + and the time taken to integrate a cell can be tuned by setting the parameter ``realtime_ratio``. + + .. cpp:member:: cell_tag_type source + + Label of the source on the cell. + + .. cpp:member:: cell_tag_type target + + Label of the target on the cell. + + .. cpp:member:: schedule time_sequence + + User-defined sequence of time points, e.g. + :cpp:class:`arb::regular_schedule`, + :cpp:class:`arb::explicit_schedule`, or + :cpp:class:`arb::poisson_schedule`. + + .. cpp:member:: double ratio + + Time taken to integrate a cell; for example, if ``realtime_ratio`` = + 2, a cell will take 2 seconds of CPU time to simulate 1 second. + +.. cpp:type:: benchmakr_cell_editor = std::function + + Callback function to update setting of LIF cells in place via the + ``simulation::edit_cell`` interace. All values are changeable, except: + ``source`` and ``target``. diff --git a/doc/cpp/cable_cell.rst b/doc/cpp/cable_cell.rst index 1b564c337..d86357959 100644 --- a/doc/cpp/cable_cell.rst +++ b/doc/cpp/cable_cell.rst @@ -151,13 +151,13 @@ Electrical properties and ion values ------------------------------------- On each cell segment, electrical and ion properties can be specified by the -:cpp:expr:`parameters` field, of type :cpp:type:`cable_cell_local_parameter_set`. +:cpp:expr:`parameters` field, of type :cpp:type:`cable_cell_parameter_set`. -The :cpp:type:`cable_cell_local_parameter_set` has the following members, +The :cpp:type:`cable_cell_parameter_set` has the following members, where an empty optional value or missing map key indicates that the corresponding value should be taken from the cell or global parameter set. -.. cpp:class:: cable_cell_local_parameter_set +.. cpp:class:: cable_cell_parameter_set .. cpp:member:: std::unordered_map ion_data @@ -192,10 +192,10 @@ value should be taken from the cell or global parameter set. Default parameters for a cell are given by the :cpp:expr:`default_parameters` field in the :cpp:type:`cable_cell` object. This is a value of type :cpp:type:`cable_cell_parameter_set`, -which extends :cpp:type:`cable_cell_local_parameter_set` by adding an additional +which extends :cpp:type:`cable_cell_parameter_set` by adding an additional field describing reversal potential computation: -.. cpp:class:: cable_cell_parameter_set: public cable_cell_local_parameter_set +.. cpp:class:: cable_cell_parameter_set: public cable_cell_parameter_set .. cpp:member:: std::unordered_map reversal_potential_method @@ -254,10 +254,11 @@ Global properties for reversal potential calculation. -For convenience, :cpp:expr:`neuron_parameter_defaults` is a predefined :cpp:type:`cable_cell_local_parameter_set` -value that holds values that correspond to NEURON defaults. To use these values, -assign them to the :cpp:expr:`default_parameters` field of the global properties -object returned in the recipe. +For convenience, :cpp:expr:`neuron_parameter_defaults` is a predefined +:cpp:type:`cable_cell_parameter_set` value that holds values that correspond to +NEURON defaults. To use these values, assign them to the +:cpp:expr:`default_parameters` field of the global properties object returned in +the recipe. .. _cppcablecell-revpot: @@ -316,3 +317,36 @@ Overriding properties locally TODO: using ``paint`` to specify electrical properties on subsections of the morphology. + +Editing Cell Parameters +----------------------- + +Whenever control is given to the user, i.e. between calls to +``simulation::run``, certain parameters oon the cell model may be altered. +Thisis done using ``simulation::edit_cell``, passing the cell's ``gid`` and a +callback object. For cable cells, this is a structure + +.. cpp:class:: cable_cell_editor + + .. cpp:type:: map = std::vector> + + .. cpp:member:: std::function + + Callback function to update setting of LIF cells in place via the + ``simulation::edit_cell`` interace. All values are changeable, except: + ``source``, ``target``, and ``V_m``. diff --git a/doc/cpp/simulation.rst b/doc/cpp/simulation.rst index 78fe0ada4..8807b181a 100644 --- a/doc/cpp/simulation.rst +++ b/doc/cpp/simulation.rst @@ -130,6 +130,13 @@ Class documentation :ref:`sampling_api` documentation. + .. cpp:function:: void edit_cell(cell_gid_type gid, std::any editor) + + Execute an editing operation ``editor`` against the cell with the given + ``gid``. The operation is a function operating on the interface type of + the cell in most cases, see the individual cell types' documentation for + the options and manner of changes. + .. cpp:function:: void remove_sampler(sampler_association_handle) Remove a sampler. diff --git a/doc/cpp/spike_source_cell.rst b/doc/cpp/spike_source_cell.rst index b083bc03c..90724e096 100644 --- a/doc/cpp/spike_source_cell.rst +++ b/doc/cpp/spike_source_cell.rst @@ -18,3 +18,9 @@ Spike source cells :param source: label of the source on the cell. :param schedule: User-defined sequence of time points + +.. cpp:type:: spike_source_cell_editor = std::function + + Callback function to update setting of LIF cells in place via the + ``simulation::edit_cell`` interace. All values are changeable, except: + ``source``. diff --git a/doc/python/simulation.rst b/doc/python/simulation.rst index 09b497509..b2cdc5e46 100644 --- a/doc/python/simulation.rst +++ b/doc/python/simulation.rst @@ -78,10 +78,54 @@ over the local and distributed hardware resources (see :ref:`pydomdec`). Then, t .. function:: update_connections(recipe) - Rebuild the connection table as described by - :py:class:`arbor.recipe::connections_on` The recipe must differ **only** - in the return value of its :py:func:`connections_on` when compared to - the original recipe used to construct the simulation object. + Rebuild the connection table as described by :py:class:`arbor.recipe` + and :py:func:`connections_on` The recipe must differ **only** in the + return value of its :py:func:`connections_on` when compared to the + original recipe used to construct the simulation object. + + .. function:: edit_cell(gid, cell_kind, editor) + + Execute an editing operation ``editor`` against the cell with the given + ``gid`` which must be of type ``cell_kind``. The operation is a function + operating on the interface type of the cell in most cases, see the + individual cell types' documentation for the options and manner of + changes. + + Cell editing is done through *functions* on the original cell mode + updating it in-place. Example + + .. code-block:: pythonic + + # assume sim is constructed s.t. gid=0 is an adex cell + sim = A.simulation(...) + + # a lambda doesn't work, we cannot assign in it + def adex_update(cell): + cell.C_m += 2 + + # pass in the kind here, since we cannot infer it from + # the Python type + sim.edit_cell(0, A.cell_kind.adex, adex_update) + + Exception here is the ``cable`` kind, which requires a ``dict`` with + specific keys + + .. code-block:: pythonic + + # assume sim is constructed s.t. gid=0 is a cable cell + @ containing the ``HH`` mechanism + sim = A.simulation(...) + + # tweak the leak conductance only on HH + def hh_update(region, mech, params): + if mech != 'hh': + return [] + return [('gl', 42)] + + # pass in the kind here, since we cannot infer it from + # the Python type + sim.edit_cell(0, A.cell_kind.cable, + {'on_density': hh_update}) .. function:: reset() @@ -160,10 +204,10 @@ over the local and distributed hardware resources (see :ref:`pydomdec`). Then, t .. function:: samples(handle) Retrieve a list of sample data associated with the given :term:`handle`. - There will be one entry in the list per probe associated with the :term:`probeset id` used when the sampling was set up. + There will be one entry in the list per probe associated with the :term:`probeset id` used when the sampling was set up. Each entry is a pair ``(samples, meta)`` where ``meta`` is the probe metadata as would be returned by ``probe_metadata(probeset_id)``, and ``samples`` contains the recorded values. - + For example, if a probe was placed on a locset describing three positions, ``samples(handle)`` will return a list of length `3`. In each element, each corresponding to a location in the locset, you'll find a tuple containing ``metadata`` and ``data``, where ``metadata`` will be a string describing the location, @@ -210,15 +254,16 @@ cells, but by default, spikes are not recorded. Recording is enabled with the the simulation object to record no spikes, all locally generated spikes, or all spikes generated by any MPI rank. -Spikes recorded during a simulation are returned as a NumPy structured datatype with two fields, -``source`` and ``time``. The ``source`` field itself is a structured datatype with two fields, -``gid`` and ``index``, identifying the threshold detector that generated the spike. +Spikes recorded during a simulation are returned as a NumPy structured datatype +with two fields, ``source`` and ``time``. The ``source`` field itself is a +structured datatype with two fields, ``gid`` and ``index``, identifying the +threshold detector that generated the spike. .. Note:: - The spikes returned by :py:func:`simulation.record` are sorted in ascending order of spike time. - Spikes that have the same spike time are sorted in ascending order of gid and local index of the - spike source. + The spikes returned by :py:func:`simulation.record` are sorted in ascending + order of spike time. Spikes that have the same spike time are sorted in + ascending order of gid and local index of the spike source. .. container:: example-code @@ -239,7 +284,7 @@ Spikes recorded during a simulation are returned as a NumPy structured datatype dt = 0.025 * U.ms sim.run(tSim, dt) - # Print the spikes: time and source + # Print the spikes: time and source for s in sim.spikes(): print(s) @@ -347,4 +392,3 @@ Example >>> [ 2.7 -62.53349949] >>> [ 2.8 -69.22068995] >>> [ 2.9 -73.41691825]] - diff --git a/example/bench/bench.cpp b/example/bench/bench.cpp index 78b685069..25c701453 100644 --- a/example/bench/bench.cpp +++ b/example/bench/bench.cpp @@ -95,7 +95,7 @@ class bench_recipe: public arb::recipe { // different MPI ranks and threads. auto sched = arb::poisson_schedule(params_.cell.spike_freq_hz*arb::units::Hz, gid); - return arb::benchmark_cell("src", "tgt", sched, params_.cell.realtime_ratio); + return arb::benchmark_cell{.source="src", .target="tgt", .time_sequence=sched, .realtime_ratio=params_.cell.realtime_ratio}; } arb::cell_kind get_cell_kind(arb::cell_gid_type gid) const override { diff --git a/example/busyring/ring.cpp b/example/busyring/ring.cpp index 1718dc5eb..3cb445b22 100644 --- a/example/busyring/ring.cpp +++ b/example/busyring/ring.cpp @@ -541,4 +541,4 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param // Make a CV between every sample in the sample tree. return {arb::morphology(tree), decor, {}, arb::cv_policy_every_segment()}; -} \ No newline at end of file +} diff --git a/example/drybench/drybench.cpp b/example/drybench/drybench.cpp index 554f36518..31da30b82 100644 --- a/example/drybench/drybench.cpp +++ b/example/drybench/drybench.cpp @@ -97,7 +97,7 @@ class tile_desc: public arb::tile { arb::util::unique_any get_cell_description(cell_gid_type gid) const override { auto gen = arb::poisson_schedule(params_.cell.spike_freq_hz*arb::units::Hz, gid); - return arb::benchmark_cell("src", "tgt", std::move(gen), params_.cell.realtime_ratio); + return arb::benchmark_cell{.source="src", .target="tgt", .time_sequence=std::move(gen), .realtime_ratio=params_.cell.realtime_ratio}; } cell_kind get_cell_kind(cell_gid_type gid) const override { diff --git a/python/cells.cpp b/python/cells.cpp index b023b47f7..d08cae54a 100644 --- a/python/cells.cpp +++ b/python/cells.cpp @@ -28,7 +28,6 @@ #include #include -#include "arbor/recipe.hpp" #include "conversion.hpp" #include "error.hpp" #include "label_dict.hpp" @@ -39,7 +38,6 @@ namespace pyarb { namespace U = arb::units; - namespace py = pybind11; template @@ -179,8 +177,7 @@ std::tuple value_and_scale(const paintable_arg& arg) { void register_cells(py::module& m) { using namespace py::literals; - using std::optional; - + using std::optional; py::class_ spike_source_cell(m, "spike_source_cell", "A spike source cell, that generates a user-defined sequence of spikes that act as inputs for other cells in the network."); py::class_ cell_cv_data(m, "cell_cv_data", @@ -987,22 +984,28 @@ void register_cells(py::module& m) { "The group of spike detectors has the label 'label', used for forming connections between cells."); cable_cell .def(py::init( - [](const arb::morphology& m, const arb::decor& d, const std::optional<::pyarb::label_dict>& l, const std::optional& p) { - if (l) return arb::cable_cell(m, d, l->dict, p); - return arb::cable_cell(m, d, {}, p); - }), - "morphology"_a, "decor"_a, "labels"_a=py::none(), "discretization"_a=py::none(), - "Construct with a morphology, decor, label dictionary, and cv policy.") + [](const arb::morphology& m, const arb::decor& d, const std::optional<::pyarb::label_dict>& l, const std::optional& p, bool editable) { + auto is_mut = editable ? arb::cable_cell_mutability::enabled : arb::cable_cell_mutability::disabled; + if (l) return arb::cable_cell(m, d, l->dict, p, is_mut); + return arb::cable_cell(m, d, {}, p, is_mut); + }), + "morphology"_a, "decor"_a, "labels"_a=py::none(), "discretization"_a=py::none(), "editable"_a=false, + "Construct with a morphology, decor, label dictionary, and cv policy.") .def(py::init( - [](const arb::segment_tree& t, const arb::decor& d, const std::optional<::pyarb::label_dict>& l, const std::optional& p) { - if (l) return arb::cable_cell({t}, d, l->dict, p); - return arb::cable_cell({t}, d, {}, p); - }), - "segment_tree"_a, "decor"_a, "labels"_a=py::none(), "discretization"_a=py::none(), - "Construct with a morphology derived from a segment tree, decor, label dictionary, and cv policy.") + [](const arb::segment_tree& t, const arb::decor& d, const std::optional<::pyarb::label_dict>& l, const std::optional& p, bool editable) { + auto is_mut = editable ? arb::cable_cell_mutability::enabled : arb::cable_cell_mutability::disabled; + if (l) return arb::cable_cell({t}, d, l->dict, p, is_mut); + return arb::cable_cell({t}, d, {}, p, is_mut); + }), + "segment_tree"_a, "decor"_a, "labels"_a=py::none(), "discretization"_a=py::none(), "editable"_a=false, + "Construct with a morphology derived from a segment tree, decor, label dictionary, and cv policy.") + .def_property_readonly("is_editable", + &arb::cable_cell::is_editable, + "Editing allowed?.") .def_property_readonly("num_branches", [](const arb::cable_cell& c) {return c.morphology().num_branches();}, "The number of unbranched cable sections in the morphology.") + // Get locations associated with a locset label. .def("locations", [](arb::cable_cell& c, const char* label) {return c.concrete_locset(arborio::parse_locset_expression(label).unwrap());}, diff --git a/python/simulation.cpp b/python/simulation.cpp index ea23429fa..739369f9d 100644 --- a/python/simulation.cpp +++ b/python/simulation.cpp @@ -1,4 +1,6 @@ #include +#include + #include #include #include @@ -7,6 +9,12 @@ #include #include #include +#include +#include +#include +#include +#include +#include #include "context.hpp" #include "error.hpp" @@ -22,14 +30,10 @@ namespace py = pybind11; namespace pyarb { // Argument type for simulation_shim::record() (see below). - -enum class spike_recording { - off, local, all -}; +enum class spike_recording { off, local, all }; // Wraps an arb::simulation object and in addition manages a set of // sampler callbacks for retrieving probe data. - class simulation_shim { std::unique_ptr sim_; std::vector spike_record_; @@ -122,6 +126,70 @@ class simulation_shim { return sim_->run(tfinal, dt); } + void edit_cell(arb::cell_gid_type gid, arb::cell_kind kind, py::object fun) { + switch (kind) { + case arb::cell_kind::adex: { + sim_->edit_cell(gid, + arb::adex_cell_editor( + [fun](arb::adex_cell& cell) -> void { + py::gil_scoped_acquire gil; + fun(py::cast(&cell, py::return_value_policy::reference)); + })); + break; + } + case arb::cell_kind::cable: { + if (!py::isinstance(fun)) throw std::runtime_error{"Wrong argument type: need dict, got " + std::string(py::str(py::type::of(fun)))}; + // this is now ok + auto dict = fun.cast(); + auto edit = arb::cable_cell_editor {}; + if (dict.contains("on_density")) { + edit.on_density = arb::density_editor( + [dict] (const arb::region&, const std::string&, const arb::parameter_map&) { + return arb::parameter_map{}; + }); + } + sim_->edit_cell(gid, edit); + break; + } + case arb::cell_kind::lif: { + sim_->edit_cell(gid, + arb::lif_cell_editor( + [fun](arb::lif_cell& x) { + py::gil_scoped_acquire gil; + fun(py::cast(&x, py::return_value_policy::reference)); + })); + + break; + } + case arb::cell_kind::spike_source: { + sim_->edit_cell(gid, + arb::spike_source_cell_editor( + [fun](arb::spike_source_cell& x) { + py::gil_scoped_acquire gil; + fun(py::cast(&x, py::return_value_policy::reference)); + })); + + break; + } + + case arb::cell_kind::benchmark: { + sim_->edit_cell(gid, + arb::benchmark_cell_editor( + [fun](arb::benchmark_cell& x) { + py::gil_scoped_acquire gil; + fun(py::cast(&x, py::return_value_policy::reference)); + })); + break; + } + default: + ARB_UNREACHABLE; + } + } + void edit_adex_cell(arb::cell_gid_type gid, std::function fun) { sim_->edit_cell(gid, fun); } + void edit_lif_cell(arb::cell_gid_type gid, arb::lif_cell_editor fun) { sim_->edit_cell(gid, fun); } + void edit_source_cell(arb::cell_gid_type gid, arb::spike_source_cell_editor fun) { sim_->edit_cell(gid, fun); } + void edit_benchmark_cell(arb::cell_gid_type gid, arb::benchmark_cell_editor fun) { sim_->edit_cell(gid, fun); } + void record(spike_recording policy) { auto spike_recorder = [this](const std::vector& spikes) { auto old_size = spike_record_.size(); @@ -282,7 +350,15 @@ void register_simulation(py::module& m, pyarb_global_ptr global_ptr) { return sim.get_probe_metadata({std::get<0>(addr), std::get<1>(addr)}); }, "Retrieve metadata associated with given probe id.", - "addr"_a) + "addr"_a) + .def("edit_cell", + &simulation_shim::edit_cell, + "Edit cell with gid and kind, passing a function (lif, adex, ...) that mutates its argument.\n" + "Cable cells expect a dictionary with function values:\n" + " {'on_density': f, 'on_synapse': g}\n" + "where each function `f`, `g` receives a location expression (region|locset), a mechanism name,\n" + "and the list of parameters [(key, value)] to edit. It is expected to return the _updated_ parameters.", + "gid"_a, "kind"_a, "lambda"_a) .def("probe_metadata", [](const simulation_shim& sim, arb::cell_gid_type gid, const arb::cell_tag_type& tag) { @@ -321,4 +397,4 @@ void register_simulation(py::module& m, pyarb_global_ptr global_ptr) { } -} // namespace pyarb +} // namespace pyar diff --git a/test/ubench/CMakeLists.txt b/test/ubench/CMakeLists.txt index 679f737e4..fd10c1c60 100644 --- a/test/ubench/CMakeLists.txt +++ b/test/ubench/CMakeLists.txt @@ -1,7 +1,6 @@ include(ExternalProject) # List of micro benchmarks to build. - set(bench_sources accumulate_functor_values.cpp default_construct.cpp @@ -10,18 +9,15 @@ set(bench_sources fvm_discretize.cpp mech_vec.cpp task_system.cpp - merge.cpp -) + merge.cpp) if(ARB_WITH_GPU) -list(APPEND bench_sources - cuda_compare_and_reduce.cu - cuda_reduce_by_key.cu -) + list(APPEND bench_sources + cuda_compare_and_reduce.cu + cuda_reduce_by_key.cu) endif() # Build benches. - foreach(bench_src ${bench_sources}) string(REGEX REPLACE "\\.[^.]*$" "" bench_exe ${bench_src}) add_executable(${bench_exe} EXCLUDE_FROM_ALL "${bench_src}") diff --git a/test/ubench/mech_vec.cpp b/test/ubench/mech_vec.cpp index b61512076..c1d954364 100644 --- a/test/ubench/mech_vec.cpp +++ b/test/ubench/mech_vec.cpp @@ -2,9 +2,6 @@ // // Start with pas (passive dendrite) mechanism -// NOTE: This targets an earlier version of the Arbor API and -// will need to be reworked in order to compile. - #include #include @@ -21,8 +18,28 @@ using namespace arb; using backend = arb::multicore::backend; using fvm_cell = arb::fvm_lowered_cell_impl; +namespace arb_bench_access { + template + struct bind { + static struct binder { + binder() { store = value; } + } init; + }; + + template + typename bind::binder bind::init; +} // namespace access + +#define ARB_BENCH_ACCESS_BIND(type, global, value)\ +namespace { using global ## _type_ = type; global ## _type_ global; }\ +template struct arb_bench_access::bind; + +ARB_BENCH_ACCESS_BIND(std::vector fvm_cell::*, private_mechanisms_ptr, &fvm_cell::mechanisms_) + +#undef ARB_BENCH_ACCESS_BIND + mechanism_ptr& find_mechanism(const std::string& name, fvm_cell& cell) { - auto &mechs = cell.mechanisms(); + auto &mechs = cell.*private_mechanisms_ptr; auto it = std::find_if(mechs.begin(), mechs.end(), [&](mechanism_ptr& m){return m->internal_name()==name;}); diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index f72dd530a..ffc034bd9 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -64,6 +64,7 @@ set(unit_sources test_any_visitor.cpp test_backend.cpp test_cable_cell.cpp + test_cell_editor.cpp test_counter.cpp test_cv_geom.cpp test_cv_layout.cpp diff --git a/test/unit/test_cell_editor.cpp b/test/unit/test_cell_editor.cpp new file mode 100644 index 000000000..66e8d0258 --- /dev/null +++ b/test/unit/test_cell_editor.cpp @@ -0,0 +1,979 @@ +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include "util/span.hpp" + +#ifdef ARB_GPU_ENABLED +constexpr int with_gpu = 0; +#else +constexpr int with_gpu = -1; +#endif + + +using namespace arb::units::literals; + +struct lif_recipe: arb::recipe { + + struct param_t { + double weight = 0; + double cm_pF = 0; + size_t n_100 = 0; + size_t n_200 = 0; + }; + + lif_recipe(double w, double cm_pf): weight(w), C_m(cm_pf*arb::units::pF) {} + + arb::cell_size_type num_cells() const override { return N; } + arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { return arb::cell_kind::lif; } + arb::util::unique_any get_cell_description(arb::cell_gid_type gid) const override { + auto cell = arb::lif_cell{"src", "tgt"}; + cell.C_m = C_m; + return cell; + } + + std::vector event_generators(arb::cell_gid_type gid) const override { + return {arb::regular_generator({"tgt"}, weight, 0_ms, 0.5_ms)}; + } + + arb::cell_size_type N = 10; + + double weight = 100; + arb::units::quantity C_m = 20_pF; +}; + +TEST(edit_lif, no_edit) { + using param_t = lif_recipe::param_t; + // check base case at 20pF + // weight c_m t=100 200 + for (const auto& param: {param_t{ 10.0, 20.0, 20, 50}, + param_t{ 100.0, 20.0, 330, 670}, + param_t{1000.0, 20.0, 500, 1000}}) { + auto rec = lif_recipe{param.weight, param.cm_pF}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + // check base case at 40pF + // weight c_m t=100 200 + for (const auto& param: {param_t{ 10.0, 40.0, 00, 0}, + param_t{ 100.0, 40.0, 250, 500}, + param_t{1000.0, 40.0, 500, 1000}}) { + auto rec = lif_recipe{param.weight, param.cm_pF}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } +} + +TEST(edit_lif, edit) { + using param_t = lif_recipe::param_t; + + arb::lif_cell_editor edit = [](arb::lif_cell& cell) { cell.C_m = 40_pF; }; + + auto ctx = arb::make_context(); + // scan group sizes + for (auto group: arb::util::make_span(1, 10)) { + auto phm = arb::partition_hint_map{ + {arb::cell_kind::lif, + arb::partition_hint{ + .cpu_group_size=std::size_t(group), + .gpu_group_size=std::size_t(group), + .prefer_gpu=true, + } + } + }; + // check transition from 20pF -> 40pF for cell gid=0 + // weight c_m t=100 200 + for (const auto& param: {param_t{ 10.0, 20.0, 20, 47}, + param_t{ 100.0, 20.0, 330, 661}, + param_t{1000.0, 20.0, 500, 1000}}) { + auto rec = lif_recipe{param.weight, param.cm_pF}; + auto ddc = arb::partition_load_balance(rec, ctx, phm); + auto sim = arb::simulation{rec, ctx, ddc}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.edit_cell(0, edit); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + // check transition from 20pF -> 40pF for half of cells + // weight c_m t=100 200 + for (const auto& param: {param_t{ 10.0, 20.0, 20, 35}, + param_t{ 100.0, 20.0, 330, 625}, + param_t{1000.0, 20.0, 500, 1000}}) { + auto rec = lif_recipe{param.weight, param.cm_pF}; + auto ddc = arb::partition_load_balance(rec, ctx, phm); + auto sim = arb::simulation{rec, ctx, ddc}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + for (arb::cell_gid_type gid = 0; gid < rec.num_cells(); gid += 2) sim.edit_cell(gid, edit); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + // check transition from 20pF -> 40pF for all cells + // weight c_m t=100 200 + for (const auto& param: {param_t{ 10.0, 20.0, 20, 20}, + param_t{ 100.0, 20.0, 330, 580}, + param_t{1000.0, 20.0, 500, 1000}}) { + auto rec = lif_recipe{param.weight, param.cm_pF}; + auto ddc = arb::partition_load_balance(rec, ctx, phm); + auto sim = arb::simulation{rec, ctx, ddc}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + for (arb::cell_gid_type gid = 0; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + // edits are idempotent + // check transition from 20pF -> 40pF for all cells + // weight c_m t=100 200 + for (const auto& param: {param_t{ 10.0, 20.0, 20, 20}, + param_t{ 100.0, 20.0, 330, 580}, + param_t{1000.0, 20.0, 500, 1000}}) { + auto rec = lif_recipe{param.weight, param.cm_pF}; + auto ddc = arb::partition_load_balance(rec, ctx, phm); + auto sim = arb::simulation{rec, ctx, ddc}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + for (arb::cell_gid_type gid = 0; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + for (arb::cell_gid_type gid = 0; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + + } +} + +TEST(edit_lif, errors) { + auto rec = lif_recipe{0, 0}; + auto sim = arb::simulation{rec}; + // Check that errors are actually thrown. + EXPECT_THROW(sim.edit_cell( 0, arb::lif_cell_editor([](auto& cell) { cell.V_m = 42_mV; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, arb::lif_cell_editor([](auto& cell) { cell.source = "foo"; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, arb::lif_cell_editor([](auto& cell) { cell.target = "foo"; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, 42), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell(42, arb::lif_cell_editor([](arb::lif_cell& cell) { cell.C_m = 40_pF; })), std::range_error); +} + +struct adex_recipe: arb::recipe { + + struct param_t { + double weight = 0; + double cm_pF = 0; + size_t n_100 = 0; + size_t n_200 = 0; + }; + + adex_recipe(double w, double cm_pf): weight(w), C_m(cm_pf*arb::units::pF) {} + + arb::cell_size_type num_cells() const override { return N; } + arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { return arb::cell_kind::adex; } + arb::util::unique_any get_cell_description(arb::cell_gid_type gid) const override { + auto cell = arb::adex_cell{"src", "tgt"}; + cell.C_m = C_m; + return cell; + } + + std::vector event_generators(arb::cell_gid_type gid) const override { + return {arb::regular_generator({"tgt"}, weight, 0_ms, 0.5_ms)}; + } + + arb::cell_size_type N = 10; + + double weight = 100; + arb::units::quantity C_m = 20_pF; +}; + +TEST(edit_adex, no_edit) { + using param_t = adex_recipe::param_t; + // check base case at 20pF + // weight c_m t=100 200 + for (const auto& param: {param_t{ 1.5, 20.0, 210, 300 }, + param_t{ 2.5, 20.0, 360, 630 }, + param_t{ 3.0, 20.0, 370, 710 }}) { + auto rec = adex_recipe{param.weight, param.cm_pF}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + // check base case at 40pF + // weight c_m t=100 200 + for (const auto& param: {param_t{ 1.5, 40.0, 130, 170 }, + param_t{ 2.5, 40.0, 310, 520 }, + param_t{ 3.0, 40.0, 340, 620 }}) { + auto rec = adex_recipe{param.weight, param.cm_pF}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } +} + +TEST(edit_adex, edit) { + using param_t = adex_recipe::param_t; + + arb::adex_cell_editor edit = [](arb::adex_cell& cell) { cell.C_m = 40_pF; }; + + auto ctx = arb::make_context(); + // scan group sizes + for (auto group: arb::util::make_span(1, 10)) { + auto phm = arb::partition_hint_map{ + {arb::cell_kind::lif, + arb::partition_hint{ + .cpu_group_size=std::size_t(group), + .gpu_group_size=std::size_t(group), + .prefer_gpu=true, + } + } + }; + // check transition from 20pF -> 40pF for cell gid=0 + // weight c_m t=100 200 + for (const auto& param: {param_t{ 1.5, 20.0, 210, 304}, + param_t{ 2.5, 20.0, 360, 634}, + param_t{ 3.0, 20.0, 370, 709}}) { + auto rec = adex_recipe{param.weight, param.cm_pF}; + auto ddc = arb::partition_load_balance(rec, ctx, phm); + auto sim = arb::simulation{rec, ctx, ddc}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.edit_cell(0, edit); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + // check transition from 20pF -> 40pF for half of cells + // weight c_m t=100 200 + for (const auto& param: {param_t{ 1.5, 20.0, 210, 320}, + param_t{ 2.5, 20.0, 360, 650}, + param_t{ 3.0, 20.0, 370, 705}}) { + auto rec = adex_recipe{param.weight, param.cm_pF}; + auto ddc = arb::partition_load_balance(rec, ctx, phm); + auto sim = arb::simulation{rec, ctx, ddc}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + for (arb::cell_gid_type gid = 0; gid < rec.num_cells(); gid += 2) sim.edit_cell(gid, edit); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + // check transition from 20pF -> 40pF for all cells + // weight c_m t=100 200 + for (const auto& param: {param_t{ 1.5, 20.0, 210, 340}, + param_t{ 2.5, 20.0, 360, 670}, + param_t{ 3.0, 20.0, 370, 700}}) { + auto rec = adex_recipe{param.weight, param.cm_pF}; + auto ddc = arb::partition_load_balance(rec, ctx, phm); + auto sim = arb::simulation{rec, ctx, ddc}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + for (arb::cell_gid_type gid = 0; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + // edits are idempotent + // check transition from 20pF -> 40pF for all cells + // weight c_m t=100 200 + for (const auto& param: {param_t{ 1.5, 20.0, 210, 340}, + param_t{ 2.5, 20.0, 360, 670}, + param_t{ 3.0, 20.0, 370, 700}}) { + auto rec = adex_recipe{param.weight, param.cm_pF}; + auto ddc = arb::partition_load_balance(rec, ctx, phm); + auto sim = arb::simulation{rec, ctx, ddc}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + for (arb::cell_gid_type gid = 0; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + for (arb::cell_gid_type gid = 0; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + break; + } +} + +TEST(edit_adex, errors) { + auto rec = adex_recipe{0, 0}; + auto sim = arb::simulation{rec}; + // Check that errors are actually thrown. + EXPECT_THROW(sim.edit_cell( 0, arb::adex_cell_editor([](auto& cell) { cell.V_m = 42_mV; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, arb::adex_cell_editor([](auto& cell) { cell.source = "foo"; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, arb::adex_cell_editor([](auto& cell) { cell.target = "foo"; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, 42), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell(42, arb::adex_cell_editor([](auto&& cell) { cell.C_m = 40_pF; })), std::range_error); +} + +struct bench_recipe: arb::recipe { + + struct param_t { + double rtr = 1.0; // real time + double nu_kHz = 1.0; // freq + size_t n_100 = 0; // spikes after N ms + size_t n_200 = 0; + }; + + bench_recipe(double r, double f): ratio(r), freq(f*arb::units::kHz) {} + + arb::cell_size_type num_cells() const override { return N; } + arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { return arb::cell_kind::benchmark; } + arb::util::unique_any get_cell_description(arb::cell_gid_type gid) const override { + return arb::benchmark_cell{.source="src", .target="tgt", .time_sequence=arb::regular_schedule(1/freq), .realtime_ratio=ratio}; + } + + arb::cell_size_type N = 10; + + double ratio = 1.0; + arb::units::quantity freq = 1_kHz; +}; + +TEST(edit_bench, no_edit) { + using param_t = bench_recipe::param_t; + // rtr nu 100 200 + for (const auto& param: {param_t{1e-4, 1.0, 1000, 2000}, // 10 cells x 100ms x 1kHz = 1000 spikes + param_t{1e-4, 2.0, 2000, 4000}, + param_t{1e-4, 4.0, 4000, 8000}}) { + auto rec = bench_recipe{param.rtr, param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } +} + +TEST(edit_bench, edit_rate) { + using param_t = bench_recipe::param_t; + // rtr nu 100 200 + for (const auto& param: {param_t{1e-4, 1.0, 1000, 2100}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{1e-4, 2.0, 2000, 4200}, + param_t{1e-4, 4.0, 4000, 8400}}) { + arb::benchmark_cell_editor edit = [&](auto& cell) { + cell.time_sequence = arb::regular_schedule(1.0/(2.0_kHz * param.nu_kHz)); + }; + auto rec = bench_recipe{param.rtr, param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + sim.edit_cell(5, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + + // rtr nu 100 200 + for (const auto& param: {param_t{1e-4, 1.0, 1000, 2500}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{1e-4, 2.0, 2000, 5000}, + param_t{1e-4, 4.0, 4000, 10000}}) { + arb::benchmark_cell_editor edit = [&](auto& cell) { cell.time_sequence = arb::regular_schedule(1.0/(2.0_kHz * param.nu_kHz)); }; + auto rec = bench_recipe{param.rtr, param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); gid += 2) sim.edit_cell(gid, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + + // rtr nu 100 200 + for (const auto& param: {param_t{1e-4, 1.0, 1000, 3000}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{1e-4, 2.0, 2000, 6000}, + param_t{1e-4, 4.0, 4000, 12000}}) { + arb::benchmark_cell_editor edit = [&](auto& cell) { cell.time_sequence = arb::regular_schedule(1.0/(2.0_kHz * param.nu_kHz)); }; + auto rec = bench_recipe{param.rtr, param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + +} + +TEST(edit_bench, edit_schedule) { + // allow -- due to the stochastic nature -- up to + auto eps = 0.15; + + using param_t = bench_recipe::param_t; + // rtr nu 100 200 + for (const auto& param: {param_t{1e-4, 1.0, 1000, 2000}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{1e-4, 2.0, 2000, 4000}, + param_t{1e-4, 4.0, 4000, 8000}}) { + arb::benchmark_cell_editor edit = [&](auto& cell) { cell.time_sequence = arb::poisson_schedule(1.0_ms/param.nu_kHz); }; + auto rec = bench_recipe{param.rtr, param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + sim.edit_cell(5, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_GE(sim.num_spikes(), param.n_200*(1.0 - eps)); + EXPECT_LE(sim.num_spikes(), param.n_200*(1.0 + eps)); + + } + + // rtr nu 100 200 + for (const auto& param: {param_t{1e-4, 1.0, 1000, 2000}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{1e-4, 2.0, 2000, 4000}, + param_t{1e-4, 4.0, 4000, 8000}}) { + arb::benchmark_cell_editor edit = [&](auto& cell) { cell.time_sequence = arb::poisson_schedule(1.0_ms/param.nu_kHz); }; + auto rec = bench_recipe{param.rtr, param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); gid += 2) sim.edit_cell(gid, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_GE(sim.num_spikes(), param.n_200*(1.0 - eps)); + EXPECT_LE(sim.num_spikes(), param.n_200*(1.0 + eps)); + } + + // rtr nu 100 200 + for (const auto& param: {param_t{1e-4, 1.0, 1000, 2000}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{1e-4, 2.0, 2000, 4000}, + param_t{1e-4, 4.0, 4000, 8000}}) { + arb::benchmark_cell_editor edit = [&](auto& cell) { cell.time_sequence = arb::poisson_schedule(1.0_ms/param.nu_kHz); }; + auto rec = bench_recipe{param.rtr, param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_GE(sim.num_spikes(), param.n_200*(1.0 - eps)); + EXPECT_LE(sim.num_spikes(), param.n_200*(1.0 + eps)); + } +} + +TEST(edit_benchmark, errors) { + auto rec = bench_recipe{1, 1}; + auto sim = arb::simulation{rec}; + // Check that errors are actually thrown. + EXPECT_THROW(sim.edit_cell( 0, arb::benchmark_cell_editor([](auto& cell) { cell.source = "foo"; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, arb::benchmark_cell_editor([](auto& cell) { cell.target = "foo"; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, 42), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell(42, arb::benchmark_cell_editor([](auto& cell) { cell.realtime_ratio = 42; })), std::range_error); +} + +TEST(edit_benchmark, do_nothing_does_nothing) { + arb::benchmark_cell_editor edit = [](auto& cell) { cell.time_sequence = arb::poisson_schedule(1.0_kHz, 42);}; + arb::benchmark_cell_editor noop = [](auto& cell) {}; + + size_t n_noop = 0; + { + auto rec = bench_recipe{1e-4, 1.0}; + auto sim = arb::simulation{rec}; + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, noop); + sim.run(200_ms, 0.1_ms); + n_noop = sim.num_spikes(); + } + size_t n_expt = 0; + { + auto rec = bench_recipe{1e-4, 1.0}; + auto sim = arb::simulation{rec}; + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + sim.run(200_ms, 0.1_ms); + n_expt = sim.num_spikes(); + } + EXPECT_EQ(n_expt, n_noop); + EXPECT_GE(n_noop, 2000); + EXPECT_LE(n_noop, 2100); +} + + +struct source_recipe: arb::recipe { + + struct param_t { + double nu_kHz = 1.0; // freq + size_t n_100 = 0; // spikes after N ms + size_t n_200 = 0; + }; + + source_recipe(double f): freq(f*arb::units::kHz) {} + + arb::cell_size_type num_cells() const override { return N; } + arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { return arb::cell_kind::spike_source; } + arb::util::unique_any get_cell_description(arb::cell_gid_type gid) const override { + return arb::spike_source_cell{"tgt", arb::regular_schedule(1/freq)}; + } + + arb::cell_size_type N = 10; + arb::units::quantity freq = 1_kHz; +}; + +TEST(edit_source, no_edit) { + using param_t = source_recipe::param_t; + // rtr nu 100 200 + for (const auto& param: {param_t{1.0, 1000, 2000}, // 10 cells x 100ms x 1kHz = 1000 spikes + param_t{2.0, 2000, 4000}, + param_t{4.0, 4000, 8000}}) { + auto rec = source_recipe{param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } +} + +TEST(edit_source, edit_rate) { + using param_t = source_recipe::param_t; + // rtr nu 100 200 + for (const auto& param: {param_t{1.0, 1000, 2100}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{2.0, 2000, 4200}, + param_t{4.0, 4000, 8400}}) { + arb::spike_source_cell_editor edit = [&](auto& cell) { cell.schedules = {arb::regular_schedule(1.0/(2.0_kHz * param.nu_kHz))}; }; + auto rec = source_recipe{param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + sim.edit_cell(5, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + + // rtr nu 100 200 + for (const auto& param: {param_t{1.0, 1000, 2500}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{2.0, 2000, 5000}, + param_t{4.0, 4000, 10000}}) { + arb::spike_source_cell_editor edit = [&](auto& cell) { cell.schedules = {arb::regular_schedule(1.0/(2.0_kHz * param.nu_kHz))}; }; + auto rec = source_recipe{param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); gid += 2) sim.edit_cell(gid, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + + // rtr nu 100 200 + for (const auto& param: {param_t{1.0, 1000, 3000}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{2.0, 2000, 6000}, + param_t{4.0, 4000, 12000}}) { + arb::spike_source_cell_editor edit = [&](auto& cell) { cell.schedules = {arb::regular_schedule(1.0/(2.0_kHz * param.nu_kHz))}; }; + auto rec = source_recipe{param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_EQ(sim.num_spikes(), param.n_200); + } + +} + +TEST(edit_source, do_nothing_does_nothing) { + arb::spike_source_cell_editor edit = [](auto& cell) { cell.schedules = {arb::poisson_schedule(1.0_kHz, 42) };}; + arb::spike_source_cell_editor noop = [](arb::spike_source_cell& cell) {}; + + size_t n_noop = 0; + { + auto rec = source_recipe{1.0}; + auto sim = arb::simulation{rec}; + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, noop); + sim.run(200_ms, 0.1_ms); + n_noop = sim.num_spikes(); + } + size_t n_expt = 0; + { + auto rec = source_recipe{1.0}; + auto sim = arb::simulation{rec}; + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + sim.run(200_ms, 0.1_ms); + n_expt = sim.num_spikes(); + } + EXPECT_EQ(n_expt, n_noop); + EXPECT_GE(n_noop, 2000); + EXPECT_LE(n_noop, 2100); +} + +TEST(edit_source, edit_schedule) { + auto eps = 0.15; + + using param_t = source_recipe::param_t; + // rtr nu 100 200 + for (const auto& param: {param_t{1.0, 1000, 2100}, // one cell adds 100ms x 2kHz, the others stay at 1Khz => 100 spikes extra + param_t{2.0, 2000, 4100}, + param_t{4.0, 4000, 8100}}) { + arb::spike_source_cell_editor edit = [&](auto& cell) { cell.schedules.push_back(arb::poisson_schedule(1.0_ms/param.nu_kHz)); }; + auto rec = source_recipe{param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + sim.edit_cell(5, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_GE(sim.num_spikes(), param.n_200*(1.0 - eps)); + EXPECT_LE(sim.num_spikes(), param.n_200*(1.0 + eps)); + } + + // rtr nu 100 200 + for (const auto& param: {param_t{1.0, 1000, 2500}, + param_t{2.0, 2000, 5000}, + param_t{4.0, 4000, 10000}}) { + arb::spike_source_cell_editor edit = [&](auto& cell) { cell.schedules.push_back(arb::poisson_schedule(1.0_ms/param.nu_kHz)); }; + auto rec = source_recipe{param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); gid += 2) sim.edit_cell(gid, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_GE(sim.num_spikes(), param.n_200*(1.0 - eps)); + EXPECT_LE(sim.num_spikes(), param.n_200*(1.0 + eps)); + } + + // rtr nu 100 200 + for (const auto& param: {param_t{1.0, 1000, 3100}, // 0-100: 1000 spikes 100-200: 1000 + ~1000 + param_t{2.0, 2000, 6100}, + param_t{4.0, 4000, 12100}}) { + arb::spike_source_cell_editor edit = [&](auto& cell) { cell.schedules.push_back(arb::poisson_schedule(1.0_ms/param.nu_kHz)); }; + auto rec = source_recipe{param.nu_kHz}; + auto sim = arb::simulation{rec}; + sim.run(100_ms, 0.1_ms); + for (auto gid = 0u; gid < rec.num_cells(); ++gid) sim.edit_cell(gid, edit); + EXPECT_EQ(sim.num_spikes(), param.n_100); + sim.run(200_ms, 0.1_ms); + EXPECT_GE(sim.num_spikes(), param.n_200*(1.0 - eps)); + EXPECT_LE(sim.num_spikes(), param.n_200*(1.0 + eps)); + } +} + +TEST(edit_spike_source, errors) { + auto rec = source_recipe{1}; + auto sim = arb::simulation{rec}; + // Check that errors are actually thrown. + EXPECT_THROW(sim.edit_cell( 0, arb::spike_source_cell_editor([](auto& cell) { cell.source = "foo"; })), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell( 0, 42), arb::bad_cell_edit); + EXPECT_THROW(sim.edit_cell(42, arb::spike_source_cell_editor([](auto& cell) {})), std::range_error); +} + +constexpr size_t N = 4; +constexpr double eps = 1e-6; +constexpr double T = 40; +constexpr double dt = 1; +constexpr size_t n_step = T/dt; +using result_t = std::vector>; + +struct cable_recipe: arb::recipe { + cable_recipe(arb::decor dec_): dec{std::move(dec_)} { + props.default_parameters = arb::neuron_parameter_defaults; + } + + arb::cell_size_type num_cells() const override { return N; } + arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { return arb::cell_kind::cable; } + arb::util::unique_any get_cell_description(arb::cell_gid_type) const override { + // Create a cable cell + // + // +------+ + // | hh |=== pas === + // +------+ + // + auto par = arb::mnpos; + auto seg = arb::segment_tree{}; + par = seg.append(par, { 0, 0, 0, 42}, {10, 0, 0, 42}, 1); // soma + par = seg.append(par, {10, 0, 0, 23}, {20, 0, 0, 23}, 2); // dendrite + auto mrf = arb::morphology{seg}; + auto lbl = arb::label_dict{}; + return arb::cable_cell{mrf, dec, lbl, cvp, arb::cable_cell_mutability::enabled}; + } + + virtual std::vector get_probes(arb::cell_gid_type) const override { return {{arb::cable_probe_membrane_voltage{arb::ls::location(0, 0.5)}, "Um"}}; } + std::any get_global_properties(arb::cell_kind) const override { return props; } + + arb::cable_cell_global_properties props; + arb::decor dec; + arb::cv_policy cvp = arb::cv_policy_max_extent(1.0_um); +}; + +testing::AssertionResult all_near(const std::vector& a, const result_t& b, int iy, double eps) { + if (a.size() != b.size()) return testing::AssertionFailure() << "sequences differ in length" + << " #expected=" << b.size() + << " #received=" << a.size(); + std::stringstream res; + res << std::setprecision(9); + for (size_t ix = 0; ix < a.size(); ++ix) { + // printf("%9.6f, ", b[ix][iy]); + auto ax = a[ix]; + auto bx = b[ix][iy]; + if (fabs(ax - bx) > eps) { + res << " elements " << ax << " and " << bx << " differ at index " << ix << ", " << iy; + break; + } + } + std::string str = res.str(); + std::cerr << res.str(); + if (str.empty()) return testing::AssertionSuccess(); + else return testing::AssertionFailure() << str; +} + +TEST(edit_cable, errors) { + auto dec = arb::decor{} + .paint(arb::reg::tagged(1), arb::density("hh", {{"gkbar", 0.036}})) + .paint(arb::reg::tagged(2), arb::density("pas")) + .place(arb::ls::location(0, 0.5), arb::i_clamp::box(10_ms, 20_ms, 100_pA), "ic1") + ; + + auto rec = cable_recipe{dec}; + auto sim = arb::simulation{rec}; + // wrong editor + EXPECT_THROW(sim.edit_cell( 0, arb::spike_source_cell_editor([](auto& cell) { cell.source = "foo"; })), arb::bad_cell_edit); + // non-existant gid + EXPECT_THROW(sim.edit_cell(42, + arb::cable_cell_editor { + .on_density = [](const arb::region& where, const std::string& what, const arb::parameter_map& old) { + return old; + }, + }), + std::range_error); +} + +TEST(edit_cable, hh) { + result_t sample_values; + sample_values.resize(n_step); + auto sampler = [&sample_values](arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { + auto gid = pm.id.gid; + for (std::size_t ix = 0; ix < n; ++ix) { + sample_values[ix][gid] = *arb::util::any_cast(samples[ix].data); + } + }; + + std::vector unedited = {-65.000000, -65.976650, -66.650927, -67.003375, -67.167843, -67.211650, -67.190473, -67.136324, -67.069902, -67.002777, -66.941313, -65.466084, -64.416894, -63.769651, -63.388147, -63.232318, -63.250644, -63.392280, -63.602123, -63.830035, -64.038928, -64.208130, -64.331221, -64.411239, -64.455972, -64.474584, -64.475653, -64.466276, -64.451815, -64.436002, -64.421192, -65.752674, -66.673119, -67.139232, -67.345711, -67.391215, -67.353235, -67.274783, -67.182795, -67.091783 }; + std::vector edited = {-65.000000, -64.712831, -64.4678048, -64.284311, -64.1428594, -64.0383144, -63.9637155, -63.9138841, -63.8840799, -63.8701359, -63.8683377, -62.3633429, -60.8617346, -59.3328255, -57.403936, -54.2563308, -46.4880246, -13.1916514, 36.3749165, -2.63933723, -32.8848184, -51.0183638, -65.6607921, -69.4867711, -69.6836194, -68.922367, -67.7889208, -66.6558301, -65.5674181, -64.5892848, -63.701057, -64.349117, -64.8180212, -65.033465, -65.1135912, -65.1018904, -65.0391091, -64.9484026, -64.8455287, -64.740033 }; + + auto dec_gkbar_0036 = arb::decor{} + .paint(arb::reg::tagged(1), arb::density("hh", {{"gkbar", 0.036}})) + .paint(arb::reg::tagged(2), arb::density("pas")) + .place(arb::ls::location(0, 0.5), arb::i_clamp::box(10_ms, 20_ms, 100_pA), "ic1") + ; + auto dec_gkbar_0008 = arb::decor{} + .paint(arb::reg::tagged(1), arb::density("hh", {{"gkbar", 0.008}})) + .paint(arb::reg::tagged(2), arb::density("pas")) + .place(arb::ls::location(0, 0.5), arb::i_clamp::box(10_ms, 20_ms, 100_pA), "ic1") + ; + + auto ctx = arb::make_context({arbenv::default_concurrency(), with_gpu}); + + // results must be invariant under the group size, even if it doesn't divide into N + for (size_t g_size = N; g_size <= N; ++g_size) { + // ... and the gid targeted + for (size_t gid = 1; gid < N; ++gid) { + auto rec_0008 = cable_recipe{dec_gkbar_0008}; + auto sim_0008 = arb::simulation{rec_0008, + ctx, + partition_load_balance(rec_0008, + ctx, + {{arb::cell_kind::cable, arb::partition_hint{.cpu_group_size=g_size}}})}; + sim_0008.add_sampler(arb::all_probes, arb::regular_schedule(dt*arb::units::ms), sampler); + // now re-write the gkbar parameter + sim_0008.edit_cell(gid, + arb::cable_cell_editor { + .on_density = arb::density_editor([] (const auto&, const auto& what, const auto&) -> arb::parameter_map { + if (what == "hh") return {{"gkbar", 0.036}}; + return {}; + })}); + sim_0008.run(T*arb::units::ms, dt*arb::units::ms); + auto samples_0008 = sample_values; + for (auto& row: sample_values) std::fill(row.begin(), row.end(), 0.0); + + + auto rec_0036 = cable_recipe{dec_gkbar_0036}; + auto sim_0036 = arb::simulation{rec_0036, + ctx, + partition_load_balance(rec_0036, + ctx, + {{arb::cell_kind::cable, arb::partition_hint{.cpu_group_size=g_size}}})}; + sim_0036.add_sampler(arb::all_probes, arb::regular_schedule(dt*arb::units::ms), sampler); + // now re-write the gkbar parameter + sim_0036.edit_cell(gid, + arb::cable_cell_editor { + .on_density = arb::density_editor([] (const auto&, const auto& what, const auto&) -> arb::parameter_map { + if (what == "hh") return {{"gkbar", 0.008}}; + return {}; + })}); + sim_0036.run(T*arb::units::ms, dt*arb::units::ms); + auto samples_0036 = sample_values; + for (auto& row: sample_values) std::fill(row.begin(), row.end(), 0.0); + + // check each unaltered cell against each other + for (unsigned col = 0; col < N; ++col) { + if (col == gid) continue; + EXPECT_TRUE(all_near(unedited, samples_0036, col, eps)); + EXPECT_TRUE(all_near( edited, samples_0008, col, eps)); + } + // now check the edited ones + EXPECT_TRUE(all_near(unedited, samples_0008, gid, eps)); + EXPECT_TRUE(all_near( edited, samples_0036, gid, eps)); + } + } +} + +TEST(edit_cable, pas) { + result_t sample_values; + sample_values.resize(n_step); + auto sampler = [&sample_values](arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { + auto gid = pm.id.gid; + for (std::size_t ix = 0; ix < n; ++ix) { + sample_values[ix][gid] = *arb::util::any_cast(samples[ix].data); + } + }; + + std::vector unedited = { -65.000000, -65.976650, -66.650927, -67.003375, -67.167843, -67.211650, -67.190473, -67.136324, -67.069902, -67.002777, -66.941313, -65.466084, -64.416894, -63.769651, -63.388147, -63.232318, -63.250644, -63.392280, -63.602123, -63.830035, -64.038928, -64.208130, -64.331221, -64.411239, -64.455972, -64.474584, -64.475653, -64.466276, -64.451815, -64.436002, -64.421192, -65.752674, -66.673119, -67.139232, -67.345711, -67.391215, -67.353235, -67.274783, -67.182795, -67.091783}; + std::vector edited = { -65.000000, -68.3113158, -69.2069298, -69.3547088, -69.3686045, -69.338811, -69.3119782, -69.2868814, -69.2672765, -69.2510752, -69.2381611, -68.6364915, -68.4719094, -68.4352472, -68.4243892, -68.4229709, -68.4226988, -68.4231265, -68.4234685, -68.423814, -68.4240821, -68.4243069, -68.4244849, -68.4246288, -68.4247431, -68.4248342, -68.4249062, -68.4249629, -68.4250073, -68.4250417, -68.4250682, -69.0158919, -69.172337, -69.2018512, -69.2068637, -69.2035358, -69.2000019, -69.1964224, -69.1935063, -69.1910237 }; + + auto dec = arb::decor{} + .paint(arb::reg::tagged(1), arb::density("hh", {{"gkbar", 0.036}})) + .paint(arb::reg::tagged(2), arb::density("pas")) + .place(arb::ls::location(0, 0.5), arb::i_clamp::box(10_ms, 20_ms, 100_pA), "ic1") + ; + auto ctx = arb::make_context({arbenv::default_concurrency(), with_gpu}); + auto rec = cable_recipe{dec}; + + // results must be invariant under the group size, even if it doesn't divide into N + for (size_t g_size = 1; g_size <= N; ++g_size) { + // ... and the gid targeted + for (size_t gid = 0; gid < N; ++gid) { + auto sim = arb::simulation{rec, + ctx, + partition_load_balance(rec, + ctx, + {{arb::cell_kind::cable, arb::partition_hint{.cpu_group_size=g_size}}})}; + sim.add_sampler(arb::all_probes, arb::regular_schedule(dt*arb::units::ms), sampler); + sim.edit_cell(gid, + arb::cable_cell_editor { + .on_density = [] (const auto&, const auto& what, const auto&) -> arb::parameter_map { + if (what == "pas") return {{"g", 0.008}}; + return {}; + }}); + + sim.run(T*arb::units::ms, dt*arb::units::ms); + // all gids present 'unedited' traces, except the one gid we targeted + for (size_t col = 0; col < N; ++col) { + if (col == gid) { + EXPECT_TRUE(all_near(edited, sample_values, col, eps)); + } + else { + EXPECT_TRUE(all_near(unedited, sample_values, col, eps)); + } + } + } + } +} + +struct cable_gen_recipe: cable_recipe { + cable_gen_recipe(const arb::decor& dec): cable_recipe(dec) {} + + std::vector event_generators(arb::cell_gid_type) const override { + return {arb::poisson_generator({"syn"}, 0.005, 0_ms, 1.0_kHz, 42)}; + } +}; + +TEST(edit_cable, expsyn) { + result_t sample_values; + sample_values.resize(n_step); + auto sampler = [&sample_values](arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { + auto gid = pm.id.gid; + for (std::size_t ix = 0; ix < n; ++ix) { + sample_values[ix][gid] = *arb::util::any_cast(samples[ix].data); + } + }; + + std::vector exp_20 = { -65.0, -54.7514995, -42.2686825, 0.137012342, 21.8149167, -45.364232, -67.9129598, -72.1109645, -68.779846, -64.7865025, -64.9626734, -61.1617474, -59.1857019, -60.6333545, -57.8448538, -55.8318071, -56.8708576, -55.3365449, -49.8457899, -41.2636889, -20.9274785, -2.98642657, -39.2739903, -64.2217407, -71.2949142, -71.7633758, -69.3955239, -66.2353185, -60.7635546, -55.6521287, -50.6940555, -41.8170467, -21.6652507, -1.49099787, -38.0130863, -63.7076861, -70.3662875, -70.8823108, -71.0628905, -65.465745, }; + std::vector exp_40 = { -65.0, -54.7514995, -41.1813592, 3.78404154, 20.0808805, -46.28961, -67.2793905, -69.8274219, -65.9098412, -61.2754815, -60.0784287, -56.4210018, -53.7311027, -53.206293, -51.0969989, -48.7490513, -47.3695081, -45.168154, -41.1002217, -36.1433987, -33.5465413, -39.3917815, -52.7034955, -64.0943558, -64.1922785, -63.7424232, -61.3787013, -58.6864932, -54.1679162, -48.8244327, -40.8628796, -24.7089839, -11.9162012, -37.5491484, -59.5190358, -68.2175105, -65.258128, -64.1614843, -63.660979, -58.9115015, }; + + auto dec_tau_20 = arb::decor{} + .paint(arb::reg::tagged(1), arb::density("hh", {{"gkbar", 0.036}})) + .paint(arb::reg::tagged(2), arb::density("pas")) + .place(arb::ls::location(0, 0.5), arb::synapse("expsyn", {{"tau", 2.0}}), "syn") + ; + auto dec_tau_40 = arb::decor{} + .paint(arb::reg::tagged(1), arb::density("hh", {{"gkbar", 0.036}})) + .paint(arb::reg::tagged(2), arb::density("pas")) + .place(arb::ls::location(0, 0.5), arb::synapse("expsyn", {{"tau", 4.0}}), "syn") + ; + + auto ctx = arb::make_context({arbenv::default_concurrency(), with_gpu}); + + // results must be invariant under the group size, even if it doesn't divide into N + for (size_t g_size = N; g_size > 1; --g_size) { + // ... and the gid targeted + for (size_t gid = 0; gid < N; ++gid) { + auto rec_20 = cable_gen_recipe{dec_tau_20}; + auto sim_20 = arb::simulation{rec_20, + ctx, + partition_load_balance(rec_20, + ctx, + {{arb::cell_kind::cable, arb::partition_hint{.cpu_group_size=g_size}}})}; + sim_20.edit_cell(gid, + arb::cable_cell_editor { + .on_synapse = [] (const auto&, const auto& what, const auto&) -> arb::parameter_map { + if (what == "expsyn") return {{"tau", 4.0}}; + return {}; + }}); + sim_20.add_sampler(arb::all_probes, arb::regular_schedule(dt*arb::units::ms), sampler); + sim_20.run(T*arb::units::ms, dt*arb::units::ms); + auto samples_20 = sample_values; + for (auto& row: sample_values) std::fill(row.begin(), row.end(), 0.0); + + auto rec_40 = cable_gen_recipe{dec_tau_40}; + auto sim_40 = arb::simulation{rec_40, + ctx, + partition_load_balance(rec_40, + ctx, + {{arb::cell_kind::cable, arb::partition_hint{.cpu_group_size=g_size}}})}; + sim_40.edit_cell(gid, + arb::cable_cell_editor { + .on_synapse = [] (const auto&, const auto& what, const auto&) -> arb::parameter_map { + if (what == "expsyn") return {{"tau", 2.0}}; + return {}; + }}); + sim_40.add_sampler(arb::all_probes, arb::regular_schedule(dt*arb::units::ms), sampler); + sim_40.run(T*arb::units::ms, dt*arb::units::ms); + auto samples_40 = sample_values; + for (auto& row: sample_values) std::fill(row.begin(), row.end(), 0.0); + + // check each unaltered cell against each other + for (unsigned col = 0; col < N; ++col) { + if (col == gid) continue; + EXPECT_TRUE(all_near(exp_20, samples_20, col, eps)); + EXPECT_TRUE(all_near(exp_40, samples_40, col, eps)); + } + // now check the edited ones + EXPECT_TRUE(all_near(exp_20, samples_40, gid, eps)); + EXPECT_TRUE(all_near(exp_40, samples_20, gid, eps)); + } + } +} + +// NOTE just tests _that_ we can, not correctness. Should probably change that. +TEST(edit_cable, can_edit_derived) { + auto dec = arb::decor{} + .paint(arb::reg::tagged(1), arb::density("hh", {{"gkbar", 0.036}})) + .paint(arb::reg::tagged(2), arb::density("pas/e=-80", {{"g", 0.1}})) + .paint(arb::reg::tagged(2), arb::density("pas/e=-70", {{"g", 0.2}})) + .paint(arb::reg::tagged(2), arb::density("pas/e=-60", {{"g", 0.3}})) + .place(arb::ls::location(0, 0.5), arb::i_clamp::box(10_ms, 20_ms, 100_pA), "ic1") + ; + auto rec = cable_recipe{dec}; + auto sim = arb::simulation{rec}; + sim.edit_cell(1, + arb::cable_cell_editor { + .on_density = [] (const auto&, const auto& what, const auto&) -> arb::parameter_map { + if (what == "pas/e=-70") return {{"g", 0.008}}; + return {}; + }}); + sim.run(T*arb::units::ms, dt*arb::units::ms); +}