diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index 3d2f9813c..2fd219990 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -1503,6 +1503,8 @@ make_point_mechanism_config(const std::unordered_map cv; diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp index e7fa2ccb3..a3b5f23d1 100644 --- a/arbor/fvm_lowered_cell_impl.hpp +++ b/arbor/fvm_lowered_cell_impl.hpp @@ -55,9 +55,6 @@ 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); } @@ -71,10 +68,16 @@ struct fvm_lowered_cell_impl: public fvm_lowered_cell { std::unique_ptr state_; // Cell state shared across mechanisms. - std::vector mechanisms_; // excludes reversal potential calculators. + std::vector density_mechanisms_; + std::vector junction_mechanisms_; + std::vector point_mechanisms_; std::vector revpot_mechanisms_; std::vector voltage_mechanisms_; + // mechs with POST_EVENTS, contains subset(s) from all above. + // NOTE: these are non-owning! + std::vector post_event_mechanisms_; + // Handles for accessing event targets. std::vector target_handles_; // Lookup table for target ids -> local target handle indices. @@ -86,9 +89,6 @@ struct fvm_lowered_cell_impl: public fvm_lowered_cell { // random number generator seed value arb_seed_type seed_; - // Flag indicating that at least one of the mechanisms implements the post_events procedure - bool post_events_ = false; - // Reset concentrations for timestep, then have all mechanisms run their update void update_ion_state(); @@ -126,18 +126,22 @@ template void fvm_lowered_cell_impl::reset() { state_->reset(); - for (auto& m: voltage_mechanisms_) m->initialize(); - for (auto& m: revpot_mechanisms_) m->initialize(); - for (auto& m: mechanisms_) m->initialize(); + for (auto& m: voltage_mechanisms_) m->initialize(); + for (auto& m: revpot_mechanisms_) m->initialize(); + for (auto& m: point_mechanisms_) m->initialize(); + for (auto& m: density_mechanisms_) m->initialize(); + for (auto& m: junction_mechanisms_) m->initialize(); update_ion_state(); state_->zero_currents(); // Note: mechanisms must be initialized again after the ion state is updated, // as mechanisms can read/write the ion_state within the initialize block - for (auto& m: revpot_mechanisms_) m->initialize(); - for (auto& m: mechanisms_) m->initialize(); - for (auto& m: voltage_mechanisms_) m->initialize(); + for (auto& m: revpot_mechanisms_) m->initialize(); + for (auto& m: point_mechanisms_) m->initialize(); + for (auto& m: voltage_mechanisms_) m->initialize(); + for (auto& m: density_mechanisms_) m->initialize(); + for (auto& m: junction_mechanisms_) m->initialize(); // NOTE: Threshold watcher reset must come after the voltage values are set, // as voltage is implicitly read by watcher to set initial state. @@ -164,28 +168,29 @@ fvm_integration_result fvm_lowered_cell_impl::integrate(const timestep_ arb_assert(state_->time == ts.t_begin()); // Update integration step time information visible to mechanisms. - for (auto& m: mechanisms_) m->set_dt(state_->dt); - for (auto& m: revpot_mechanisms_) m->set_dt(state_->dt); - for (auto& m: voltage_mechanisms_) m->set_dt(state_->dt); + const auto dt = state_->dt; + for (auto& m: density_mechanisms_) m->set_dt(dt); + for (auto& m: point_mechanisms_) m->set_dt(dt); + for (auto& m: revpot_mechanisms_) m->set_dt(dt); + for (auto& m: voltage_mechanisms_) m->set_dt(dt); + for (auto& m: junction_mechanisms_) m->set_dt(dt); // Update any required reversal potentials based on ionic concentrations - for (auto& m: revpot_mechanisms_) { - m->update_current(); - } + for (auto& m: revpot_mechanisms_) m->update_current(); PE(zero); state_->zero_currents(); PL(zero); // Deliver events and accumulate mechanism current contributions. - // Mark all events due before (but not including) the end of this time step (state_->time_to) for delivery state_->mark_events(); - for (auto& m: mechanisms_) { - // apply the events and drop them afterwards - state_->deliver_events(*m); - m->update_current(); - } + for (auto& m: point_mechanisms_) state_->deliver_events(*m); + + // Now we are fre to update currents. + for (auto& m: density_mechanisms_) m->update_current(); + for (auto& m: point_mechanisms_) m->update_current(); + for (auto& m: junction_mechanisms_) m->update_current(); // Add stimulus current contributions. // NOTE: performed after dt, time_to calculation, in case we want to @@ -204,11 +209,16 @@ fvm_integration_result fvm_lowered_cell_impl::integrate(const timestep_ state_->integrate_cable_state(); PL(cable); - // Integrate mechanism state for density - for (auto& m: mechanisms_) { - state_->update_prng_state(*m); - m->update_state(); - } + // advance PRNG states. + for (auto& m: density_mechanisms_) state_->update_prng_state(*m); + for (auto& m: point_mechanisms_) state_->update_prng_state(*m); + for (auto& m: voltage_mechanisms_) state_->update_prng_state(*m); + for (auto& m: junction_mechanisms_) state_->update_prng_state(*m); + + // Integrate mechanism state + for (auto& m: density_mechanisms_) m->update_state(); + for (auto& m: point_mechanisms_) m->update_state(); + for (auto& m: junction_mechanisms_) m->update_state(); // Update ion concentrations. PE(ionupdate); @@ -218,10 +228,7 @@ fvm_integration_result fvm_lowered_cell_impl::integrate(const timestep_ // voltage mechs run now; after the cable_solver, but before the // threshold test for (auto& m: voltage_mechanisms_) m->update_current(); - for (auto& m: voltage_mechanisms_) { - state_->update_prng_state(*m); - m->update_state(); - } + for (auto& m: voltage_mechanisms_) m->update_state(); // Update time and test for spike threshold crossings. PE(threshold); @@ -229,9 +236,7 @@ fvm_integration_result fvm_lowered_cell_impl::integrate(const timestep_ PL(threshold); PE(post); - if (post_events_) { - for (auto& m: mechanisms_) m->post_event(); - } + for (auto& m: post_event_mechanisms_) m->post_event(); PL(post); // Advance epoch @@ -251,7 +256,11 @@ fvm_integration_result fvm_lowered_cell_impl::integrate(const timestep_ template void fvm_lowered_cell_impl::update_ion_state() { state_->ions_init_concentration(); - for (auto& m: mechanisms_) m->update_ions(); + // NOTE neither voltage nor revpot mechanisms can alter ions + // TODO Can junction/point? + for (auto& m: point_mechanisms_) m->update_ions(); + for (auto& m: density_mechanisms_) m->update_ions(); + for (auto& m: junction_mechanisms_) m->update_ions(); } template @@ -398,10 +407,9 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid fvm_mechanism_data mech_data = fvm_build_mechanism_data(global_props, cells, gids, gj_conns, D, context_); // 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; - if (post_events_) { + if (mech_data.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)) { for (auto lid: util::make_span(fvm_info.num_sources[gids[cell_idx]])) { @@ -487,7 +495,8 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid } } mechptr_by_name[name] = mech.get(); - mechanisms_.emplace_back(mech.release()); + if (config.has_post_event) post_event_mechanisms_.push_back(mech.get()); + point_mechanisms_.emplace_back(mech.release()); break; } case arb_mechanism_kind_gap_junction: { @@ -500,7 +509,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(); - mechanisms_.emplace_back(mech.release()); + junction_mechanisms_.emplace_back(mech.release()); break; } case arb_mechanism_kind_voltage: { @@ -532,7 +541,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(); - mechanisms_.emplace_back(mech.release()); + density_mechanisms_.emplace_back(mech.release()); break; } case arb_mechanism_kind_reversal_potential: { @@ -555,7 +564,7 @@ fvm_lowered_cell_impl::initialize(const std::vector& gid util::transform_view(gids, [&](cell_gid_type i) { return fvm_info.num_targets[i]; })); - + reset(); return fvm_info; } @@ -950,7 +959,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/modcc/module.cpp b/modcc/module.cpp index 918099683..fec9b15c1 100644 --- a/modcc/module.cpp +++ b/modcc/module.cpp @@ -249,7 +249,6 @@ bool Module::semantic() { }; // ... except for write_ions(), which we construct by hand here. - expr_list_type ion_assignments; auto sym_to_id = [](Symbol* sym) -> expression_ptr { @@ -568,6 +567,9 @@ bool Module::semantic() { error(pprintf("NET_RECEIVE takes at most one argument (Arbor limitation!)"), net_rec_api.first->location()); } net_rec_api.first->body(net_rec_api.second->body()->clone()); + if (kind_ != moduleKind::point) { + error("NET_RECEIVE only allowed on POINT_PROCESS", net_rec_api.first->location()); + } if (net_rec_api.second) { for (auto &s: net_rec_api.second->body()->statements()) { if (s->is_assignment()) { @@ -600,6 +602,9 @@ bool Module::semantic() { if (post_events_) { auto post_events_api = make_empty_api_method("post_event_api", "post_event"); post_events_api.first->body(post_events_api.second->body()->clone()); + if (kind_ != moduleKind::point) { + error("POST_EVENT only allowed on POINT_PROCESS", post_events_api.first->location()); + } } // check voltage mechanisms before rev pot ... otherwise we are in trouble diff --git a/test/ubench/mech_vec.cpp b/test/ubench/mech_vec.cpp index b61512076..c047ed5ad 100644 --- a/test/ubench/mech_vec.cpp +++ b/test/ubench/mech_vec.cpp @@ -21,16 +21,59 @@ using namespace arb; using backend = arb::multicore::backend; using fvm_cell = arb::fvm_lowered_cell_impl; -mechanism_ptr& find_mechanism(const std::string& name, fvm_cell& cell) { - auto &mechs = cell.mechanisms(); - auto it = std::find_if(mechs.begin(), - mechs.end(), +// Subvert class access protections. Demo: +// +// class foo { +// int secret = 7; +// }; +// +// int foo::* secret_mptr; +// template class access::bind; +// +// int seven = foo{}.*secret_mptr; +// +// Or with shortcut define (places global in anonymous namespace): +// +// ACCESS_BIND(int foo::*, secret_mptr, &foo::secret) +// +// int seven = foo{}.*secret_mptr; +namespace arb_access { + template + struct bind { + static struct binder { + binder() { store = value; } + } init; + }; + + template + typename bind::binder bind::init; +} // namespace access + +#define ACCESS_BIND(type, global, value)\ +namespace { using global ## _type_ = type; global ## _type_ global; }\ +template struct arb_access::bind; + +ACCESS_BIND(std::vector fvm_cell::*, private_pp_mechanisms_ptr, &fvm_cell::point_mechanisms_) +ACCESS_BIND(std::vector fvm_cell::*, private_de_mechanisms_ptr, &fvm_cell::density_mechanisms_) + + +const mechanism_ptr& find_mechanism(const std::string& name, fvm_cell& cell) { + auto& pp = cell.*private_pp_mechanisms_ptr; + auto it = std::find_if(pp.begin(), + pp.end(), [&](mechanism_ptr& m){return m->internal_name()==name;}); - if (it==mechs.end()) { + auto& de = cell.*private_de_mechanisms_ptr; + if (it == pp.end()) { + it = std::find_if(de.begin(), + de.end(), + [&](mechanism_ptr& m){return m->internal_name()==name;}); + } + if (it==de.end()) { std::cerr << "couldn't find mechanism with name " << name << "\n"; exit(1); } return *it; + return *it; } class recipe_expsyn_1_branch: public recipe { @@ -103,7 +146,7 @@ class recipe_pas_1_branch: public recipe { // Add soma. auto s0 = tree.append(arb::mnpos, {0,0,-soma_radius,soma_radius}, {0,0,soma_radius,soma_radius}, 1); - // Add dendrite + // Add dendrite auto s1 = tree.append(s0, {0,0,soma_radius,dend_radius}, 3); tree.append(s1, {0,0,soma_radius+dend_length,dend_radius}, 3); diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp index ba34b1de3..dcb6dd000 100644 --- a/test/unit/test_fvm_lowered.cpp +++ b/test/unit/test_fvm_lowered.cpp @@ -47,22 +47,19 @@ ACCESS_BIND(std::unique_ptr fvm_cell::*, private_state_ptr, &fvm_c using matrix = arb::multicore::cable_solver; -ACCESS_BIND(std::vector fvm_cell::*, private_mechanisms_ptr, &fvm_cell::mechanisms_) +ACCESS_BIND(std::vector fvm_cell::*, private_pp_mechanisms_ptr, &fvm_cell::point_mechanisms_) +ACCESS_BIND(std::vector fvm_cell::*, private_de_mechanisms_ptr, &fvm_cell::density_mechanisms_) arb::mechanism* find_mechanism(fvm_cell& fvcell, const std::string& name) { - for (auto& mech: fvcell.*private_mechanisms_ptr) { - if (mech->internal_name()==name) { - return mech.get(); - } + for (auto& mech: fvcell.*private_de_mechanisms_ptr) { + if (mech->internal_name() == name) return mech.get(); + } + for (auto& mech: fvcell.*private_pp_mechanisms_ptr) { + if (mech->internal_name() == name) return mech.get(); } return nullptr; } -arb::mechanism* find_mechanism(fvm_cell& fvcell, int index) { - auto& mechs = fvcell.*private_mechanisms_ptr; - return index<(int)mechs.size()? mechs[index].get(): nullptr; -} - using namespace arb; class gap_recipe_0: public recipe { @@ -438,18 +435,15 @@ TEST(fvm_lowered, derived_mechs) { { // Test initialization and global parameter values. - fvm_cell fvcell(*context); fvcell.initialize({0, 1, 2}, rec); // Both mechanisms will have the same internal name, "test_kin1". - using fvec = std::vector; fvec tau_values; - for (auto& mech: fvcell.*private_mechanisms_ptr) { + for (auto& mech: fvcell.*private_de_mechanisms_ptr) { ASSERT_TRUE(mech); EXPECT_EQ("test_kin1"s, mech->internal_name()); - tau_values.push_back(mechanism_global(mech, "tau")); } util::sort(tau_values); @@ -550,7 +544,9 @@ TEST(fvm_lowered, read_valence) { fvm_cell fvcell(*context); fvcell.initialize({0}, rec); - auto cr_mech_ptr = find_mechanism(fvcell, 0); + // NOTE: the intternal name is not changed, so use the _base_ name here + auto cr_mech_ptr = find_mechanism(fvcell, "test_ca_read_valence"); + ASSERT_NE(cr_mech_ptr, nullptr); auto cr_record_z = mechanism_field(cr_mech_ptr, "record_z"); ASSERT_EQ(7.0, cr_record_z.at(0)); }