FEATURE: fixed dt (#2053)

Motivation
----------------

So far, multi-compartment cable cells use a variable time step: the
user-specified maximum time step $dt$ may be shortened based on incoming
events, such as spikes, generators, and (exact) sampling, so that the
event delivery is performed exactly at the requested times.

All events with time $t_s$ are now delivered at simulation
time $t_i = t_0 + i ~dt$ if
$t_s  \in \left[ t_i , t_i + dt\right).$

Removing exact delivery increases the speed of the simulation due to
elimination of small steps, makes the numerics independent of presence
of sampling, and also leads to a number of code simplifications. In
particular, we have no more need for

- integration domains,
- event binning,
- sampling policies.

Main Changes
---------------------

The new event delivery logic is mainly implemented in
`mc_cell_group::advance` and
`fvm_lowered_cell_impl<Backend>::integrate`.

Incidental Changes
----------------------------
- `arb_mechanism_ppack` now exposes a scalar time and scalar time step:
`t` and `dt`
- All cells advance with the same increments, i.e. no separate
integration domains
- Deliverable events now organized in vector of vector of events.
Outermost vector has one element per `mechanism_id`, inner vector has
one element per time step interval.
- This vector is then used to initialize a `event_stream` for each
`shared_state::mech_storage`, while earlier, we had only one
`multi_event_stream` per `shared_state`.
- Helps with simpler and faster event delivery, and exposes more
parallelism on the GPU
- GPU `event_stream` is now nearly identical to multi-core
implementation (marking etc. is done on CPU) except for optional stable
sorting according to `mechanism_index`, which is done in parallel using
a `thread_group`.
  - Less maintenance
- `shared_state` gains some responsibility to handle events:
`register_events`, `mark_events`, `deliver_events`
  - Better encapsulation

---------

Co-authored-by: bcumming <[email protected]>
Co-authored-by: Thorsten Hater <[email protected]>

arbor/CMakeLists.txt

@@ -17,7 +17,6 @@ set(arbor_sources
     domain_decomposition.cpp
     execution_context.cpp
     gpu_context.cpp
-    event_binner.cpp
     fvm_layout.cpp
     fvm_lowered_cell_impl.cpp
     hardware/memory.cpp
@@ -76,8 +75,6 @@ if(ARB_WITH_GPU)
         backends/gpu/diffusion.cu
         backends/gpu/fine.cu
         backends/gpu/matrix_solve.cu
-        backends/gpu/multi_event_stream.cpp
-        backends/gpu/multi_event_stream.cu
         backends/gpu/shared_state.cu
         backends/gpu/forest.cpp
         backends/gpu/stimulus.cu

arbor/backends/event.hpp

@@ -2,6 +2,8 @@
 
 #include <arbor/common_types.hpp>
 #include <arbor/fvm_types.hpp>
+#include <arbor/mechanism_abi.h>
+#include <arbor/generic_event.hpp>
 
 // Structures for the representation of event delivery targets and
 // staged events.
@@ -13,11 +15,10 @@ namespace arb {
 struct target_handle {
     cell_local_size_type mech_id;    // mechanism type identifier (per cell group).
     cell_local_size_type mech_index; // instance of the mechanism
-    cell_size_type intdom_index;     // which integration domain (acts as index into arrays)
 
     target_handle() = default;
-    target_handle(cell_local_size_type mech_id, cell_local_size_type mech_index, cell_size_type intdom_index):
-        mech_id(mech_id), mech_index(mech_index), intdom_index(intdom_index) {}
+    target_handle(cell_local_size_type mech_id, cell_local_size_type mech_index):
+        mech_id(mech_id), mech_index(mech_index) {}
 };
 
 struct deliverable_event {
@@ -30,23 +31,23 @@ struct deliverable_event {
         time(time), weight(weight), handle(handle) {}
 };
 
+template<>
+struct has_event_index<deliverable_event> : public std::true_type {};
+
 // Stream index accessor function for multi_event_stream:
-inline cell_size_type event_index(const deliverable_event& ev) {
-    return ev.handle.intdom_index;
+inline cell_local_size_type event_index(const arb_deliverable_event_data& ed) {
+    return ed.mech_index;
 }
 
-// Subset of event information required for mechanism delivery.
-struct deliverable_event_data {
-    cell_local_size_type mech_id;    // same as target_handle::mech_id
-    cell_local_size_type mech_index; // same as target_handle::mech_index
-    float weight;
-};
-
 // Delivery data accessor function for multi_event_stream:
-inline deliverable_event_data event_data(const deliverable_event& ev) {
-    return {ev.handle.mech_id, ev.handle.mech_index, ev.weight};
+inline arb_deliverable_event_data event_data(const deliverable_event& ev) {
+    return {ev.handle.mech_index, ev.weight};
 }
 
+inline arb_deliverable_event_stream make_event_stream_state(arb_deliverable_event_data* begin,
+                                                            arb_deliverable_event_data* end) {
+    return {begin, end};
+}
 
 // Sample events (raw values from back-end state).
 
@@ -59,17 +60,11 @@ struct raw_probe_info {
 
 struct sample_event {
     time_type time;
-    cell_size_type intdom_index;  // which integration domain probe is on
     raw_probe_info raw;           // event payload: what gets put where on sample
 };
 
 inline raw_probe_info event_data(const sample_event& ev) {
     return ev.raw;
 }
 
-inline cell_size_type event_index(const sample_event& ev) {
-    return ev.intdom_index;
-}
-
-
 } // namespace arb

arbor/backends/event_delivery.md

@@ -18,9 +18,7 @@ the lowered cell, and is passed by reference as a parameter to the mechanism
 
 Target handles are used by the lowered cell implementation to identify a particular mechanism
 instance that can receive events — ultimately via `net_receive` — and corresponding simulated
-cell. The cell information is given as an index into the cell group collection of cells,
-and is used to group events by integration domain (we have one domain per cell in each cell
-group).
+cell. The cell information is given as an index into the cell group collection of cells.
 
 Target handles are represented by the `target_handle` struct, opaque to `mc_cell_group`,
 but created in the `fvm_multicell` for each point mechanism (synapse) in the cell group.
@@ -33,34 +31,38 @@ a `target_handle` describing their destination, and a weight.
 
 ### Back-end event streams
 
-`backend::multi_event_stream` represents a set (one per cell/integration domain)
-of event streams. There is one `multi_event_stream` per lowered cell.
-
+`backend::multi_event_stream` represents a set of event streams. There is one `multi_event_stream`
+per mechanism storage `backend::shared_state::mech_storage`.
 From the perspective of the lowered cell, it must support the methods below.
-In the following, `time` is a `view` or `const_view` of an array with one
-element per stream.
 
 *  `void init(const std::vector<deliverable_event>& staged_events)`
 
-   Take a copy of the staged events (which must be ordered by increasing event time)
-   and initialize the streams by gathering the events by cell.
+   Take a copy of the staged events (which must be partitioned by the event index (a.k.a. the stream
+   index), and ordered by increasing event time within each partition) and initialize the streams.
 
 *  `bool empty() const`
 
    Return true if and only if there are no un-retired events left in any stream.
 
-*  `void clear()`
+*  `size_type n_streams() const`
 
-   Retire all events, leaving the `multi_event_stream` in an empty state.
+   Number of partitions/streams.
+
+*  `size_type n_remaining() const`
 
-*  `void mark_until_after(const_view time)`
+   Number of remaining un-retired events among all streams.
+
+*  `size_type n_marked() const`
+
+   Number of marked events among all streams.
+
+*  `void clear()`
 
-   For all streams, mark events for delivery in the _i_ th stream with event time ≤ _time[i]_.
+   Retire all events, leaving the `multi_event_stream` in an empty state.
 
-*  `void event_times_if_before(view time) const`
+*  `void mark_until_after(arb_value_type t_until)`
 
-   For each stream, set _time[i]_ to the time of the next event time in the _i_ th stream
-   if such an event exists and has time less than _time[i]_.
+   For all streams, mark events for delivery with event time ≤ `t_until`.
 
 *  `void drop_marked_events()`
 
@@ -72,52 +74,33 @@ element per stream.
 Event delivery is performed as part of the integration loop within the lowered
 cell implementation. The interface is provided by the `multi_event_stream`
 described above, together with the mechanism method that handles the delivery proper,
-`mechanism::deliver_events` and a `backend` static method that computes the
-integration step end time before considering any pending events.
+`mechanism::deliver_events` and `backend` methods
+`shared_state::register_events`, `shared_state::mark_events` and `shared_state::deliver_events`.
+Events are considered as pending when their event time is within one time step of the current time:
+`event_time > time - dt/2 && event_time <= time + dt/2`.
 
 For `fvm_multicell` one integration step comprises:
 
 1.  Events for each cell that are due at that cell's corresponding time are
-    gathered with `events_.mark_events(time_)` where `time_` is a
-    `const_view` of the cell times and `events_` is a reference to the
-    `backend::multi_event_stream` object.
+    gathered with `state_>mark_events(step_midpoint)` where `step_midpoint` is
+    the current time plus half a time step (upper bound for pending event times). This method, in
+    turn, calls the `multi_event_stream::mark_until_after(step_midpoint)`.
 
 2.  Each mechanism is requested to deliver to itself any marked events that
-    are associated with that mechanism, via the virtual
-    `mechanism::apply_events(backend::multi_event_stream&)` method.
+    are associated with that mechanism, via the
+    `shared_state::deliver_events(const mechanism&)` method. This method eventually calls the virtual
+    `mechanism::apply_events(backend::multi_event_stream&)` method and retires the deliverd events
+    using `multi_event_stream::drop_marked_events`.
 
     This action must precede the computation of mechanism current contributions
     with `mechanism::compute_currents()`.
 
-3.  Marked events are discarded with `events_.drop_marked_events()`.
-
-4.  The upper bound on the integration step stop time `time_to_` is
-    computed via `backend::update_time_to(time_to_, time_, dt_max_, tfinal_)`,
-    as the minimum of the per-cell time `time_` plus `dt_max_` and
-    the final integration stop time `tfinal_`.
-
-5.  The integration step stop time `time_to_` is reduced to match any
-    pending events on each cell with `events_.event_times_if_before(time_to)`.
-
-6.  The solver matrix is assembled and solved to compute the voltages, using the
+3.  The solver matrix is assembled and solved to compute the voltages, using the
     newly computed currents and integration step times.
 
-7.  The mechanism states are updated with `mechanism::advance_state()`.
+4.  The mechanism states are updated with `mechanism::advance_state()`.
 
-8.  The cell times `time_` are set to the integration step stop times `time_to_`.
+5.  The cell times `time` are set to the integration step stop times `time_to`.
 
-9.  Spike detection for the last integration step is performed via the
+6.  Spike detection for the last integration step is performed via the
     `threshold_watcher_` object.
-
-## Consequences for the integrator
-
-Towards the end of the integration period, an integration step may have a zero _dt_
-for one or more cells within the group, and this needs to be handled correctly:
-
-*   Generated mechanism `advance_state()` methods should be numerically correct with
-    zero _dt_; a possibility is to guard the integration step with a _dt_ check.
-
-*   Matrix assemble and solve must check for zero _dt_. In the FVM `multicore`
-    matrix implementation, zero _dt_ sets the diagonal element to zero and the
-    rhs to the voltage; the solve procedure then ignores cells with a zero
-    diagonal element.

arbor/backends/event_stream_base.hpp

@@ -0,0 +1,64 @@
+#pragma once
+
+#include <type_traits>
+#include <vector>
+
+#include <arbor/generic_event.hpp>
+
+#include "backends/event.hpp"
+#include "backends/event_stream_state.hpp"
+
+namespace arb {
+
+template <typename Event, typename Span>
+class event_stream_base {
+public: // member types
+    using size_type = std::size_t;
+    using event_type = Event;
+    using event_time_type = ::arb::event_time_type<Event>;
+    using event_data_type = ::arb::event_data_type<Event>;
+
+protected: // private member types
+    using span_type = Span;
+
+    static_assert(std::is_same<decltype(std::declval<span_type>().begin()), event_data_type*>::value);
+    static_assert(std::is_same<decltype(std::declval<span_type>().end()), event_data_type*>::value);
+
+protected: // members
+    std::vector<event_data_type> ev_data_;
+    std::vector<span_type> ev_spans_;
+    size_type index_ = 0;
+
+public:
+    event_stream_base() = default;
+
+    // returns true if the currently marked time step has no events
+    bool empty() const {
+        return ev_spans_.empty() || ev_data_.empty() || !index_ || index_ > ev_spans_.size() ||
+            !ev_spans_[index_-1].size();
+    }
+
+    void mark() {
+        index_ += (index_ <= ev_spans_.size() ? 1 : 0);
+    }
+
+    auto marked_events() {
+        using std::begin;
+        using std::end;
+        if (empty()) {
+            return make_event_stream_state((event_data_type*)nullptr, (event_data_type*)nullptr);
+        } else {
+            return make_event_stream_state(begin(ev_spans_[index_-1]), end(ev_spans_[index_-1]));
+        }
+    }
+
+    // clear all previous data
+    void clear() {
+        ev_data_.clear();
+        ev_spans_.clear();
+        index_ = 0;
+    }
+
+};
+
+} // namespace arb

arbor/backends/event_stream_state.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <arbor/fvm_types.hpp>
+
+// Pointer representation of event stream marked event state,
+// common across CPU and GPU backends.
+
+namespace arb {
+
+template <typename EvData>
+struct event_stream_state {
+    using value_type = EvData;
+
+    const value_type* begin_marked;   // offset to beginning of marked events
+    const value_type* end_marked;     // offset to one-past-end of marked events
+
+    std::size_t size() const noexcept {
+        return end_marked - begin_marked;
+    }
+
+    bool empty() const noexcept {
+        return (size() == 0u);
+    }
+};
+
+template <typename EvData>
+inline event_stream_state<EvData> make_event_stream_state(EvData* begin, EvData* end) {
+    return {begin, end};
+}
+
+} // namespace arb