Refactor FVM Lowered

arbor/fvm_layout.cpp

@@ -1503,6 +1503,8 @@ make_point_mechanism_config(const std::unordered_map<std::string, mlocation_map<
              });
 
         auto config = make_mechanism_config(info, arb_mechanism_kind_point);
+        config.has_post_event = info.post_events;
+
         // Do coalesce?
         if (!info.random_variables.size() && info.linear && data.coalesce) {
             for (auto& [k, _v]: parameters) {

arbor/fvm_layout.hpp

@@ -147,7 +147,7 @@ struct fvm_diffusion_info {
     fvm_diffusion_info(value_type d): default_value(d) {}
     fvm_diffusion_info(): fvm_diffusion_info{0.0} {}
 };
-    
+
 struct fvm_cv_discretization {
     using size_type = arb_size_type;
     using index_type = arb_index_type;
@@ -209,6 +209,8 @@ struct fvm_mechanism_config {
 
     arb_mechanism_kind kind;
 
+    bool has_post_event = false;
+
     // Ordered CV indices where mechanism is present; may contain
     // duplicates for point mechanisms.
     std::vector<index_type> cv;

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<mechanism_ptr>& mechanisms() { return mechanisms_; }
-
     ARB_SERDES_ENABLE(fvm_lowered_cell_impl<Backend>, 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<shared_state> state_; // Cell state shared across mechanisms.
 
-    std::vector<mechanism_ptr> mechanisms_; // excludes reversal potential calculators.
+    std::vector<mechanism_ptr> density_mechanisms_;
+    std::vector<mechanism_ptr> junction_mechanisms_;
+    std::vector<mechanism_ptr> point_mechanisms_;
     std::vector<mechanism_ptr> revpot_mechanisms_;
     std::vector<mechanism_ptr> voltage_mechanisms_;
 
+    // mechs with POST_EVENTS, contains subset(s) from all above.
+    // NOTE: these are non-owning!
+    std::vector<mechanism*> post_event_mechanisms_;
+
     // Handles for accessing event targets.
     std::vector<target_handle> 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 <typename Backend>
 void fvm_lowered_cell_impl<Backend>::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<Backend>::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<Backend>::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,20 +228,15 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::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);
         state_->test_thresholds();
         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<Backend>::integrate(const timestep_
 template <typename Backend>
 void fvm_lowered_cell_impl<Backend>::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 <typename Backend>
@@ -398,10 +407,9 @@ fvm_lowered_cell_impl<Backend>::initialize(const std::vector<cell_gid_type>& 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<arb_index_type> 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<Backend>::initialize(const std::vector<cell_gid_type>& 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<Backend>::initialize(const std::vector<cell_gid_type>& 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<Backend>::initialize(const std::vector<cell_gid_type>& 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<Backend>::initialize(const std::vector<cell_gid_type>& 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));

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

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<backend>;
 
-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 foo::*, secret_mptr, &foo::secret>;
+//
+//     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 <typename V, V& store, V value>
+    struct bind {
+        static struct binder {
+            binder() { store = value; }
+        } init;
+    };
+
+    template <typename V, V& store, V value>
+    typename bind<V, store, value>::binder bind<V, store, value>::init;
+} // namespace access
+
+#define ACCESS_BIND(type, global, value)\
+namespace { using global ## _type_ = type; global ## _type_ global; }\
+template struct arb_access::bind<type, global, value>;
+
+ACCESS_BIND(std::vector<arb::mechanism_ptr> fvm_cell::*, private_pp_mechanisms_ptr, &fvm_cell::point_mechanisms_)
+ACCESS_BIND(std::vector<arb::mechanism_ptr> 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);