Add support for axonal delays

doc/htmldoc/developer_space/axonal_delays.rst

@@ -0,0 +1,86 @@
+.. _axonal_delays_dev:
+
+Axonal Delays
+=============
+
+Adding axonal delays to NEST is non-trivial when it comes to their interaction with spike-timing dependent plasticity (STDP).
+Axonal delays lower than their dendritic counterpart are non-problematic, however larger axonal delays cause causality
+issues due to the way how and when pre- and post-synaptic spikes are processed in NEST by synapses implementing STDP weight dynamics.
+
+If a pre-synaptic spike is processed at a synapse, it will also process all post-synaptic spikes that reached the synapse
+between the last and current pre-synaptic spike. Weight changes due facilitation (post-synaptic spike following a
+pre-synaptic one) or depression (pre-synaptic spike following a post-synaptic one) are only relevant at the time when pre-synaptic spikes
+reach the synapse, as this is the only point in time when the exact weight is of importance. Post-synaptic spikes can
+therefore be archived in the post-synaptic neuron until the next pre-synaptic spike is processed by the synapse.
+As all pre-synaptic spikes are delivered to their target synapse and neuron right after they have been communicated,
+they might be processed before they would actually reach the synapse when taking axonal delays into account.
+If the axonal delay is now larger than the dendritic delay, post-synaptic
+spikes occurring at time `t` will reach the synapse before pre-synaptic spikes occurring before `t`,
+but might not be taken into account by the pre-synaptic spike, if it was already communicated,
+and thus delivered, before `t`. Each pre-synaptic spike sent over a connection
+with a predominant axonal delay must therefore also process post-synaptic spikes which have not yet occurred,
+but could be emitted in the future. Multiple implementations were implemented and
+benchmarked before coming to the conclusion that the implementation at hand should be used inside NEST.
+
+The main idea of the axonal delays implementation in NEST is based on the fact, that neurons only emit few spikes per second.
+It should thus be rare that a post-synaptic spike occurs right after a pre-synaptic one in the critical region before
+the pre-synaptic spike reaches the synapse, but has already been processed. In typical networks, there will most likely
+only be few occurrences where causality becomes an issue. In order to still guarantee correct synaptic weights,
+incorrect STDP weight changes are rolled back, re-calculated, and the weight of pre-synaptic spike, which already reached
+the target neuron's ring buffer, is corrected. Undoing the STDP weight changes and re-calculating them obviously comes
+with a cost, however as only few such occurrences are to be expected, this solution is more efficient than restructuring
+the kernel to make sure axonal delays are always handled correctly (see Alternative implementations).
+
+Changes to the kernel and neuron models
+---------------------------------------
+
+Introducing axonal delays changes the way the min- and max-delays must be calculated, as they are now a combination of
+dendritic and axonal delays. The default value for the delay which is now referring to the dendritic delay remains 1,
+while the default value for axonal_delay is set to 0. In the default case, purely dendritic delay is assumed.
+
+The ``ArchivingNode`` was made axonal-delay-aware. Each pre-synaptic spike after which a correction could potentially follow,
+will be archived in the post-synaptic neuron in a dynamic ring-buffer-like structure. Post-synaptic spikes will then
+trigger a correction for all relevant pre-synaptic spikes in this buffer. The way spikes are received at a neuron is
+model-dependent, as the implementation of spike accumulation and buffering until being processed might vary between
+neuron models. Neurons models will therefore also have to handle correction of previously handled spikes differently.
+In the simplest case, all incoming spikes to a neuron are simply accumulated in a single scalar value per time slot.
+A correction of a previously handled spike would therefore just subtract the previous, wrong weight and add the new,
+corrected weight. Therefore, simply sending another spike with the difference of the old and new weight would be
+sufficient in this case. However, some neurons might have different buffers for spikes being sent over inhibitory and
+excitatory connections, which could be distinguished by the sign of the weight. If a new spike is now sent to correct
+an old one, the sign might be negative even though both the old and new weight were originally positive, the new weight
+is just smaller. In such a case, the spike would be accumulated in the wrong buffer.
+
+Instead of sending a regular ``SpikeEvent`` to signal a correction, a ``CorrectionSpikeEvent`` is sent. Overloading the handle
+function now allows handling the correction in the correct way, depending on the model implementation.
+Furthermore, neuron models must now call ``ArchivingNode::pre_run_hook_()`` in their derived pre_run_hook implementation
+and call ``reset_correction_entries_stdp_ax_delay_()`` at the end of their update implementation.
+Currently, only the ``iaf_psc_alpha`` neuron model supports STDP with axonal delays.
+All other neurons will act as if the delay of incoming connections was purely dendritic.
+
+Synapse models only support dendritic delay by default. If axonal delays are required, the synapse model must be derived
+from ``AxonalDelayConnection`` instead of ``Connection``. The ``AxonalDelayConnection`` is derived from ``Connection`` and adds a single
+double-precision member for the axonal delay. The main differences compared to synapses with purely dendritic delays are
+different handling of delays inside the send function and the addition of the ``correct_synapse_stdp_ax_delay`` which is
+called by the ``ConnectionManager`` when a synapse needs to re-calculate its weight given a new post-synaptic spike and a previous pre-synaptic one.
+Currently, only the ``stdp_pl_synapse_hom_ax_delay`` synapse model supports axonal delays.
+
+Changes to the python interface
+-------------------------------
+
+In general, the kernel was made axonal-delay-aware and this is reflected in the user interface, as it is now possible
+to set the ``names::dendritic_delay`` and ``names::axonal_delay`` for each synapse (given that the synapse model is
+derived from ``AxonalDelayConnection``).
+
+Remaining work
+---------------
+
+
+Currently, only one neuron and synapse model are supporting axonal delays. All neuron models that support STDP could
+also support axonal delays, without sacrificing performance, changing their behavior, or requiring more memory, but need
+to be adapted slightly (i.e., implement handle for ``CorrectionSpikeEvent``, call ``ArchivingNode::pre_run_hook_`` and call
+``reset_correction_entries_stdp_ax_delay_``).
+
+Existing STDP synapse models need one version with and one without axonal delays. Alternatively, synapse models could
+be templatized to either use only dendritic or dendritic and axonal delays. However, this branching should be resolved
+at compile time to not negatively impact performance.

doc/htmldoc/developer_space/index.rst

@@ -6,9 +6,9 @@ Developer space
 Here is all documentation pertaining to the development of NEST.
 It is documentation for anyone needing to touch the code or documentation.
 
-.. grid:: 3 
+.. grid:: 3
 
-  .. grid-item-card:: 
+  .. grid-item-card::
        :link-type: ref
        :link: dev_install
        :class-card: sd-bg-success sd-text-white
@@ -112,7 +112,7 @@ Developer guides
 
 .. toctree::
    :maxdepth: 1
-   :hidden: 
+   :hidden:
    :glob:
 
    workflows/*
@@ -122,4 +122,5 @@ Developer guides
    guidelines/styleguide/vim_support_sli
    templates/*
    sli_docs/index
+   axonal_delays
    cppcomments

doc/htmldoc/examples/index.rst

@@ -219,6 +219,10 @@ PyNEST examples
            * :doc:`/auto_examples/eprop_plasticity/eprop_supervised_regression_sine-waves`
            * :doc:`/auto_examples/eprop_plasticity/eprop_supervised_classification_neuromorphic_mnist`
 
+    .. grid-item-card:: Axonal delays
+           :img-top: ../static/img/nest_logo-faded.png
+
+           * :doc:`/auto_examples/axonal_delays`
 
 .. grid:: 1 1 2 3
 
@@ -375,3 +379,4 @@ PyNEST examples
    ../auto_examples/pong/run_simulations
    ../auto_examples/pong/pong
    ../auto_examples/pong/generate_gif
+   ../auto_examples/axonal_delays

doc/htmldoc/get-started_index.rst

@@ -159,6 +159,7 @@ More topics
 
        * :ref:`sim_gap_junctions`
        * :ref:`weight_normalization`
+       * :ref:`delays`
 
 
 

doc/htmldoc/synapses/delays.rst

@@ -0,0 +1,73 @@
+.. _delays:
+
+Delays
+======
+
+
+In NEST, transmission delays are specified with the ``delay`` parameter.
+Delays are considered fully dendritic by all built-in models and therefore, the ``delay`` parameter is still used by
+most models.
+
+Since NEST 3.9, it is also possible to specify both explicit, heterogeneous axonal and dendritic delays for models
+supporting this feature. This is useful for STDP models and other models relying on the timing of spike arrival at the
+synapse.
+
+Currently, only ``stdp_pl_synapse_hom_ax_delay`` supports explicitly specifying axonal and dendritic delays with the
+``axonal_delay`` and ``dendtritic_delay`` parameters. For STDP with predominant axonal delays, neuron models must be
+adjusted to correctly handle these delays. At this point, only ``iaf_psc_alpha`` supports STDP with predominant axonal
+delays.
+
+When using ``stdp_pl_synapse_hom_ax_delay``:
+
+- The parameter ``delay`` is no longer valid. This is to prevent ambiguity between the two types of delays.
+- The parameter names ``dendritic_delay`` and ``axonal_delay`` have to be used to specify delay.
+- If these parameters are not explicitly provided, then the default values are used:
+
+  ``dendritic_delay: 1.0``  and ``axonal_delay: 0.0``.
+- If only axonal delay is provided and no dendritic delay, the dendritic delay is assumed to be 0 and vice-versa.
+
+
+Use of ``axonal_delay`` and ``dendritic_delay`` is the same as ``delay``:
+
+
+**Using syn_spec**
+
+.. code-block:: python
+
+   nest.Create("iaf_psc_alpha")
+   nest.Connect(neuron, neuron, syn_spec={"synapse_model": "stdp_pl_synapse_hom", "delay": 1.0})
+
+.. code-block:: python
+
+   nest.Create("iaf_psc_alpha")
+   nest.Connect(neuron, neuron, syn_spec=
+               {"synapse_model": "stdp_pl_synapse_hom_ax_delay", "axonal_delay": 1.0, "dendritic_delay": 1.0})
+
+**Using SetStatus**
+
+.. code-block:: python
+
+   conn = nest.Connect(neuron, neuron, syn_spec={"synapse_model": "stdp_pl_synapse_hom"})
+   nest.SetStatus(conn, {"delay": 1.0})
+
+.. code-block:: python
+
+      conn = nest.Connect(neuron, neuron, syn_spec={"synapse_model": "stdp_pl_synapse_hom_ax_delay"})
+      nest.SetStatus(conn, {"axonal_delay": 1.0, "dendritic_delay": 1.0})
+
+**Using SetDefaults**
+
+.. code-block:: python
+
+   nest.SetDefaults("stdp_pl_synapse_hom", {"delay": 1.0})
+
+.. code-block:: python
+
+   nest.SetDefaults("stdp_pl_synapse_hom_ax_delay", {"axonal_delay": 1.0, "dendritic_delay": 1.0})
+
+
+.. seealso::
+
+   :doc:`Example using axonal delays </auto_examples/axonal_delays>`
+
+   For details on further developments see :ref:`axonal_delays_dev`.

libnestutil/nest_types.h

@@ -91,9 +91,15 @@ constexpr uint8_t NUM_BITS_LCID = 27U;
 constexpr uint8_t NUM_BITS_PROCESSED_FLAG = 1U;
 constexpr uint8_t NUM_BITS_MARKER_SPIKE_DATA = 2U;
 constexpr uint8_t NUM_BITS_LAG = 14U;
-constexpr uint8_t NUM_BITS_DELAY = 21U;
 constexpr uint8_t NUM_BITS_NODE_ID = 62U;
 
+// These types are used in delay_types.h and denote the space available for the dendritic and axonal portions of the
+// total transmission delay. The delay is only split into two parts for selected synapse types.
+// Given that axonal delays can be much larger than dendritic/backpropagation delays, they require more bits.
+constexpr uint8_t NUM_BITS_DENDRITIC_DELAY = 14U;
+constexpr uint8_t NUM_BITS_AXONAL_DELAY = sizeof( unsigned int ) * 8 - NUM_BITS_DENDRITIC_DELAY;
+
+
 // Maximally allowed values for bitfields
 
 constexpr uint64_t MAX_LCID = generate_max_value( NUM_BITS_LCID );

models/ac_generator.cpp

@@ -56,30 +56,29 @@ RecordablesMap< ac_generator >::create()
 {
   insert_( Name( names::I ), &ac_generator::get_I_ );
 }
-}
 
 /* ----------------------------------------------------------------
  * Default constructors defining default parameters and state
  * ---------------------------------------------------------------- */
 
-nest::ac_generator::Parameters_::Parameters_()
+ac_generator::Parameters_::Parameters_()
   : amp_( 0.0 )     // pA
   , offset_( 0.0 )  // pA
   , freq_( 0.0 )    // Hz
   , phi_deg_( 0.0 ) // degree
 {
 }
 
-nest::ac_generator::Parameters_::Parameters_( const Parameters_& p )
+ac_generator::Parameters_::Parameters_( const Parameters_& p )
   : amp_( p.amp_ )
   , offset_( p.offset_ )
   , freq_( p.freq_ )
   , phi_deg_( p.phi_deg_ )
 {
 }
 
-nest::ac_generator::Parameters_&
-nest::ac_generator::Parameters_::operator=( const Parameters_& p )
+ac_generator::Parameters_&
+ac_generator::Parameters_::operator=( const Parameters_& p )
 {
   if ( this == &p )
   {
@@ -94,19 +93,19 @@ nest::ac_generator::Parameters_::operator=( const Parameters_& p )
   return *this;
 }
 
-nest::ac_generator::State_::State_()
+ac_generator::State_::State_()
   : y_0_( 0.0 )
   , y_1_( 0.0 ) // pA
   , I_( 0.0 )   // pA
 {
 }
 
-nest::ac_generator::Buffers_::Buffers_( ac_generator& n )
+ac_generator::Buffers_::Buffers_( ac_generator& n )
   : logger_( n )
 {
 }
 
-nest::ac_generator::Buffers_::Buffers_( const Buffers_&, ac_generator& n )
+ac_generator::Buffers_::Buffers_( const Buffers_&, ac_generator& n )
   : logger_( n )
 {
 }
@@ -116,7 +115,7 @@ nest::ac_generator::Buffers_::Buffers_( const Buffers_&, ac_generator& n )
  * ---------------------------------------------------------------- */
 
 void
-nest::ac_generator::Parameters_::get( DictionaryDatum& d ) const
+ac_generator::Parameters_::get( DictionaryDatum& d ) const
 {
   ( *d )[ names::amplitude ] = amp_;
   ( *d )[ names::offset ] = offset_;
@@ -125,14 +124,14 @@ nest::ac_generator::Parameters_::get( DictionaryDatum& d ) const
 }
 
 void
-nest::ac_generator::State_::get( DictionaryDatum& d ) const
+ac_generator::State_::get( DictionaryDatum& d ) const
 {
   ( *d )[ names::y_0 ] = y_0_;
   ( *d )[ names::y_1 ] = y_1_;
 }
 
 void
-nest::ac_generator::Parameters_::set( const DictionaryDatum& d, Node* node )
+ac_generator::Parameters_::set( const DictionaryDatum& d, Node* node )
 {
   updateValueParam< double >( d, names::amplitude, amp_, node );
   updateValueParam< double >( d, names::offset, offset_, node );
@@ -145,7 +144,7 @@ nest::ac_generator::Parameters_::set( const DictionaryDatum& d, Node* node )
  * Default and copy constructor for node
  * ---------------------------------------------------------------- */
 
-nest::ac_generator::ac_generator()
+ac_generator::ac_generator()
   : StimulationDevice()
   , P_()
   , S_()
@@ -154,7 +153,7 @@ nest::ac_generator::ac_generator()
   recordablesMap_.create();
 }
 
-nest::ac_generator::ac_generator( const ac_generator& n )
+ac_generator::ac_generator( const ac_generator& n )
   : StimulationDevice( n )
   , P_( n.P_ )
   , S_( n.S_ )
@@ -168,20 +167,20 @@ nest::ac_generator::ac_generator( const ac_generator& n )
  * ---------------------------------------------------------------- */
 
 void
-nest::ac_generator::init_state_()
+ac_generator::init_state_()
 {
   StimulationDevice::init_state();
 }
 
 void
-nest::ac_generator::init_buffers_()
+ac_generator::init_buffers_()
 {
   StimulationDevice::init_buffers();
   B_.logger_.reset();
 }
 
 void
-nest::ac_generator::pre_run_hook()
+ac_generator::pre_run_hook()
 {
   B_.logger_.init();
 
@@ -206,7 +205,7 @@ nest::ac_generator::pre_run_hook()
 }
 
 void
-nest::ac_generator::update( Time const& origin, const long from, const long to )
+ac_generator::update( Time const& origin, const long from, const long to )
 {
   long start = origin.get_steps();
 
@@ -231,7 +230,7 @@ nest::ac_generator::update( Time const& origin, const long from, const long to )
 }
 
 void
-nest::ac_generator::handle( DataLoggingRequest& e )
+ac_generator::handle( DataLoggingRequest& e )
 {
   B_.logger_.handle( e );
 }
@@ -241,7 +240,7 @@ nest::ac_generator::handle( DataLoggingRequest& e )
  * ---------------------------------------------------------------- */
 
 void
-nest::ac_generator::set_data_from_stimulation_backend( std::vector< double >& input_param )
+ac_generator::set_data_from_stimulation_backend( std::vector< double >& input_param )
 {
   Parameters_ ptmp = P_; // temporary copy in case of errors
 
@@ -264,3 +263,5 @@ nest::ac_generator::set_data_from_stimulation_backend( std::vector< double >& in
   // if we get here, temporary contains consistent set of properties
   P_ = ptmp;
 }
+
+} // namespace nest

models/aeif_cond_alpha.cpp

@@ -426,6 +426,8 @@ nest::aeif_cond_alpha::init_buffers_()
 void
 nest::aeif_cond_alpha::pre_run_hook()
 {
+  ArchivingNode::pre_run_hook_();
+
   // ensures initialization in case mm connected after Simulate
   B_.logger_.init();