eprop_iaf_adapt_bsshslm_2020.h

Go to the documentation of this file.

/*
 *  eprop_iaf_adapt_bsshslm_2020.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef EPROP_IAF_ADAPT_BSSHSLM_2020_H
#define EPROP_IAF_ADAPT_BSSHSLM_2020_H
 
// nestkernel
#include "connection.h"
#include "eprop_archiving_node_impl.h"
#include "eprop_archiving_node_recurrent.h"
#include "event.h"
#include "nest_types.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"
 
namespace nest
{
 
/* BeginUserDocs: neuron, e-prop plasticity, current-based, integrate-and-fire, adaptation, hard threshold, Bellec
 
Short description
+++++++++++++++++
 
Current-based leaky integrate-and-fire neuron model with delta-shaped or exponentially filtered
postsynaptic currents and threshold adaptation for e-prop plasticity
 
Description
+++++++++++
 
``eprop_iaf_adapt_bsshslm_2020`` is an implementation of a leaky integrate-and-fire
neuron model with delta-shaped postsynaptic currents and threshold adaptation
used for eligibility propagation (e-prop) plasticity.
 
E-prop plasticity was originally introduced and implemented in TensorFlow in :footcite:p:`Bellec2020`.
 
The suffix ``_bsshslm_2020`` follows the NEST convention to indicate in the
model name the paper that introduced it by the first letter of the authors' last
names and the publication year.
 
 .. note::
   The neuron dynamics of the ``eprop_iaf_adapt_bsshslm_2020`` model (excluding
   e-prop plasticity and the threshold adaptation) are similar to the neuron
   dynamics of the ``iaf_psc_delta`` model, with minor differences, such as the
   propagator of the post-synaptic current and the voltage reset upon a spike.
 
The membrane voltage time course :math:`v_j^t` of the neuron :math:`j` is given by:
 
.. math::
  v_j^t &= \alpha v_j^{t-1} + \zeta \sum_{i \neq j} W_{ji}^\text{rec} z_i^{t-1}
    + \zeta \sum_i W_{ji}^\text{in} x_i^t - z_j^{t-1} v_\text{th} \,, \\
  \alpha &= e^{ -\frac{ \Delta t }{ \tau_\text{m} } } \,, \\
  \zeta &=
    \begin{cases}
      1 \\
      1 - \alpha
    \end{cases} \,, \\
 
where :math:`W_{ji}^\text{rec}` and :math:`W_{ji}^\text{in}` are the recurrent and
input synaptic weight matrices, and :math:`z_i^{t-1}` is the recurrent presynaptic
state variable, while :math:`x_i^t` represents the input at time :math:`t`.
 
Descriptions of further parameters and variables can be found in the table below.
 
The threshold adaptation is given by:
 
.. math::
  A_j^t &= v_\text{th} + \beta a_j^t \,, \\
  a_j^t &= \rho a_j^{t-1} + z_j^{t-1} \,, \\
  \rho &= e^{-\frac{ \Delta t }{ \tau_\text{a} }} \,. \\
 
The spike state variable is expressed by a Heaviside function:
 
.. math::
  z_j^t = H \left( v_j^t - A_j^t \right) \,. \\
 
If the membrane voltage crosses the adaptive threshold voltage :math:`A_j^t`, a spike is
emitted and the membrane voltage is reduced by :math:`v_\text{th}` in the next
time step. After the time step of the spike emission, the neuron is not
able to spike for an absolute refractory period :math:`t_\text{ref}`.
 
An additional state variable and the corresponding differential equation
represents a piecewise constant external current.
 
See the documentation on the :doc:`iaf_psc_delta<../models/iaf_psc_delta/>` neuron model
for more information on the integration of the subthreshold dynamics.
 
The change of the synaptic weight is calculated from the gradient :math:`g` of
the loss :math:`E` with respect to the synaptic weight :math:`W_{ji}`:
:math:`\frac{ \text{d}E }{ \text{d} W_{ij} }`
which depends on the presynaptic
spikes :math:`z_i^{t-1}`, the surrogate gradient or pseudo-derivative
of the spike state variable with respect to the postsynaptic membrane
voltage :math:`\psi_j^t` (the product of which forms the eligibility
trace :math:`e_{ji}^t`), and the learning signal :math:`L_j^t` emitted
by the readout neurons.
 
.. math::
  \frac{ \text{d} E }{ \text{d} W_{ji} } &= \sum_t L_j^t \bar{e}_{ji}^t \,, \\
  e_{ji}^t &= \psi_j^t \left( \bar{z}_i^{t-1} - \beta \epsilon_{ji,a}^{t-1} \right) \,, \\
  \epsilon^{t-1}_{ji,\text{a}} &= \psi_j^{t-1} \bar{z}_i^{t-2} + \left( \rho - \psi_j^{t-1} \beta \right)
    \epsilon^{t-2}_{ji,a} \,. \\
 
.. include:: ../models/eprop_iaf.rst
   :start-after: .. start_surrogate-gradient-functions
   :end-before: .. end_surrogate-gradient-functions
 
The eligibility trace and the presynaptic spike trains are low-pass filtered
with the following exponential kernels:
 
.. math::
  \bar{e}_{ji}^t &= \mathcal{F}_\kappa \left( e_{ji}^t \right) \,, \\
  \kappa &= e^{ -\frac{\Delta t }{ \tau_\text{m,out} }} \,, \\
  \bar{z}_i^t &= \mathcal{F}_\alpha(z_i^t) \,, \\
  \mathcal{F}_\alpha \left( z_i^t \right) &= \alpha \mathcal{F}_\alpha \left( z_i^{t-1} \right) + z_i^t \,, \\
  \mathcal{F}_\alpha \left( z_i^0 \right) &= z_i^0 \,, \\
 
where :math:`\tau_\text{m,out}` is the membrane time constant of the readout neuron.
 
Furthermore, a firing rate regularization mechanism keeps the average firing
rate :math:`f^\text{av}_j` of the postsynaptic neuron close to a target firing rate
:math:`f^\text{target}`. The gradient :math:`g_\text{reg}` of the regularization loss :math:`E_\text{reg}`
with respect to the synaptic weight :math:`W_{ji}` is given by:
 
.. math::
  \frac{ \text{d} E_\text{reg} }{ \text{d} W_{ji} }
    = c_\text{reg} \sum_t \frac{ 1 }{ T n_\text{trial} }
    \left( f^\text{target} - f^\text{av}_j \right) e_{ji}^t \,, \\
 
where :math:`c_\text{reg}` is a constant scaling factor and the average
is taken over the time that passed since the previous update, that is, the number of
trials :math:`n_\text{trial}` times the duration of an update interval :math:`T`.
 
The overall gradient is given by the addition of the two gradients.
 
For more information on e-prop plasticity, see the documentation on the other e-prop models:
 
 * :doc:`eprop_iaf_bsshslm_2020<../models/eprop_iaf_bsshslm_2020/>`
 * :doc:`eprop_readout_bsshslm_2020<../models/eprop_readout_bsshslm_2020/>`
 * :doc:`eprop_synapse_bsshslm_2020<../models/eprop_synapse_bsshslm_2020/>`
 * :doc:`eprop_learning_signal_connection_bsshslm_2020<../models/eprop_learning_signal_connection_bsshslm_2020/>`
 
Details on the event-based NEST implementation of e-prop can be found in :footcite:p:`KorcsakGorzo2025`.
 
Parameters
++++++++++
 
The following parameters can be set in the status dictionary.
 
=========================== ======= ======================= ================ ===================================
**Neuron parameters**
----------------------------------------------------------------------------------------------------------------
Parameter                   Unit    Math equivalent         Default          Description
=========================== ======= ======================= ================ ===================================
``adapt_beta``                      :math:`\beta`                        1.0 Prefactor of the threshold
                                                                             adaptation
``adapt_tau``               ms      :math:`\tau_\text{a}`               10.0 Time constant of the threshold
                                                                             adaptation
``C_m``                     pF      :math:`C_\text{m}`                 250.0 Capacitance of the membrane
``E_L``                     mV      :math:`E_\text{L}`                 -70.0 Leak / resting membrane potential
``I_e``                     pA      :math:`I_\text{e}`                   0.0 Constant external input current
``regular_spike_arrival``   Boolean                                 ``True`` If ``True``, the input spikes
                                                                             arrive at the end of the time step,
                                                                             if ``False`` at the beginning
                                                                             (determines PSC scale)
``t_ref``                   ms      :math:`t_\text{ref}`                 2.0 Duration of the refractory period
``tau_m``                   ms      :math:`\tau_\text{m}`               10.0 Time constant of the membrane
``V_min``                   mV      :math:`v_\text{min}`    negative maximum Absolute lower bound of the
                                                            value            membrane voltage
                                                            representable by
                                                            ``double``
                                                            type in C++
``V_th``                    mV      :math:`v_\text{th}`                -55.0 Spike threshold voltage
=========================== ======= ======================= ================ ===================================
 
=============================== ======= ======================= ================== =============================
**E-prop parameters**
----------------------------------------------------------------------------------------------------------------
Parameter                       Unit    Math equivalent         Default            Description
=============================== ======= ======================= ================== =============================
``flush_event_send_interval``   ms                              maximum value      Interval since previous event
                                                                representable by   after which a flush event is
                                                                ``double`` type in sent
                                                                C++
``c_reg``                               :math:`c_\text{reg}`                 0.0   Coefficient of firing rate
                                                                                   regularization
``f_target``                    Hz      :math:`f^\text{target}`             10.0   Target firing rate of rate
                                                                                   regularization
``surrogate_gradient_function``         :math:`\psi`            "piecewise_linear" Surrogate gradient /
                                                                                   pseudo-derivative function
                                                                                   ["piecewise_linear",
                                                                                   "exponential",
                                                                                   "fast_sigmoid_derivative",
                                                                                   "arctan_derivative"]
``surrogate_gradient_height``           :math:`\gamma`                        0.3  Height scaling of surrogate
                                                                                   gradient / pseudo-derivative
                                                                                   of membrane voltage
``surrogate_gradient_width``            :math:`1/\beta`                       1.0  Width scaling of surrogate
                                                                                   gradient / pseudo-derivative
                                                                                   of membrane voltage
=============================== ======= ======================= ================== =============================
 
Recordables
+++++++++++
 
The following state variables evolve during simulation and can be recorded.
 
================== ==== =============== ============= ========================
**Neuron state variables and recordables**
------------------------------------------------------------------------------
State variable     Unit Math equivalent Initial value Description
================== ==== =============== ============= ========================
``adaptation``          :math:`a_j`               0.0 Adaptation variable
``V_m``              mV :math:`v_j`             -70.0 Membrane voltage
``V_th_adapt``       mV :math:`A_j`             -55.0 Adapting spike threshold
================== ==== =============== ============= ========================
 
====================== ==== =============== ============= =========================================
**E-prop state variables and recordables**
---------------------------------------------------------------------------------------------------
State variable         Unit Math equivalent Initial value Description
====================== ==== =============== ============= =========================================
``learning_signal``      pA :math:`L_j`               0.0 Learning signal
``surrogate_gradient``      :math:`\psi_j`            0.0 Surrogate gradient / pseudo-derivative of
                                                          membrane voltage
====================== ==== =============== ============= =========================================
 
Usage
+++++
 
This model can only be used in combination with the other e-prop models
and the network architecture requires specific wiring, input, and output.
The usage is demonstrated in several
:doc:`supervised regression and classification tasks <../auto_examples/eprop_plasticity/index>`
reproducing among others the original proof-of-concept tasks in :footcite:p:`Bellec2020`.
 
References
++++++++++
 
.. footbibliography::
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, LearningSignalConnectionEvent, DataLoggingRequest
 
See also
++++++++
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: eprop_iaf_adapt_bsshslm_2020
 
EndUserDocs */
 
void register_eprop_iaf_adapt_bsshslm_2020( const std::string& name );
 
class eprop_iaf_adapt_bsshslm_2020 : public EpropArchivingNodeRecurrent< true >
{
 
public:
  eprop_iaf_adapt_bsshslm_2020();
 
  eprop_iaf_adapt_bsshslm_2020( const eprop_iaf_adapt_bsshslm_2020& );
 
  using Node::handle;
  using Node::handles_test_event;
 
  size_t send_test_event( Node&, size_t, synindex, bool ) override;
 
  void handle( SpikeEvent& ) override;
  void handle( CurrentEvent& ) override;
  void handle( LearningSignalConnectionEvent& ) override;
  void handle( DataLoggingRequest& ) override;
 
  size_t handles_test_event( SpikeEvent&, size_t ) override;
  size_t handles_test_event( CurrentEvent&, size_t ) override;
  size_t handles_test_event( LearningSignalConnectionEvent&, size_t ) override;
  size_t handles_test_event( DataLoggingRequest&, size_t ) override;
 
  void get_status( Dictionary& ) const override;
  void set_status( const Dictionary& ) override;
 
private:
  void init_buffers_() override;
  void pre_run_hook() override;
 
  void update( Time const&, const long, const long ) override;
 
  double compute_gradient( std::vector< long >&, const long, const long, const double, const bool ) override;
 
  long get_shift() const override;
  bool is_eprop_recurrent_node() const override;
 
  friend class RecordablesMap< eprop_iaf_adapt_bsshslm_2020 >;
 
  friend class UniversalDataLogger< eprop_iaf_adapt_bsshslm_2020 >;
 
  struct Parameters_
  {
    double adapt_beta_;
 
    double adapt_tau_;
 
    double C_m_;
 
    double c_reg_;
 
    double E_L_;
 
    double f_target_;
 
    double I_e_;
 
    bool regular_spike_arrival_;
 
    std::string surrogate_gradient_function_;
 
    double surrogate_gradient_height_;
 
    double surrogate_gradient_width_;
 
    double t_ref_;
 
    double tau_m_;
 
    double V_min_;
 
    double V_th_;
 
    Parameters_();
 
    void get( Dictionary& ) const;
 
    double set( const Dictionary&, Node* );
  };
 
  struct State_
  {
    double adapt_;
 
    double v_th_adapt_;
 
    double learning_signal_;
 
    int r_;
 
    double surrogate_gradient_;
 
    double i_in_;
 
    double v_m_;
 
    double z_;
 
    double z_in_;
 
    State_();
 
    void get( Dictionary&, const Parameters_& ) const;
 
    void set( const Dictionary&, const Parameters_&, double, Node* );
  };
 
  struct Buffers_
  {
    Buffers_( eprop_iaf_adapt_bsshslm_2020& );
 
    Buffers_( const Buffers_&, eprop_iaf_adapt_bsshslm_2020& );
 
    RingBuffer spikes_;
 
    RingBuffer currents_;
 
    UniversalDataLogger< eprop_iaf_adapt_bsshslm_2020 > logger_;
  };
 
  struct Variables_
  {
    double P_v_m_;
 
    double P_z_in_;
 
    double P_i_in_;
 
    double P_adapt_;
 
    int RefractoryCounts_;
  };
 
  double
  get_v_m_() const
  {
    return S_.v_m_ + P_.E_L_;
  }
 
  double
  get_surrogate_gradient_() const
  {
    return S_.surrogate_gradient_;
  }
 
  double
  get_learning_signal_() const
  {
    return S_.learning_signal_;
  }
 
  double
  get_v_th_adapt_() const
  {
    return S_.v_th_adapt_ + P_.E_L_;
  }
 
  double
  get_adaptation_() const
  {
    return S_.adapt_;
  }
 
  // the order in which the structure instances are defined is important for speed
 
  Parameters_ P_;
 
  State_ S_;
 
  Variables_ V_;
 
  Buffers_ B_;
 
  static RecordablesMap< eprop_iaf_adapt_bsshslm_2020 > recordablesMap_;
};
 
inline long
eprop_iaf_adapt_bsshslm_2020::get_shift() const
{
  return offset_gen_ + delay_in_rec_;
}
 
inline bool
eprop_iaf_adapt_bsshslm_2020::is_eprop_recurrent_node() const
{
  return true;
}
 
inline size_t
eprop_iaf_adapt_bsshslm_2020::send_test_event( Node& target, size_t receptor_type, synindex, bool )
{
  SpikeEvent e;
  e.set_sender( *this );
  return target.handles_test_event( e, receptor_type );
}
 
inline size_t
eprop_iaf_adapt_bsshslm_2020::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
 
  return 0;
}
 
inline size_t
eprop_iaf_adapt_bsshslm_2020::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
 
  return 0;
}
 
inline size_t
eprop_iaf_adapt_bsshslm_2020::handles_test_event( LearningSignalConnectionEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
 
  return 0;
}
 
inline size_t
eprop_iaf_adapt_bsshslm_2020::handles_test_event( DataLoggingRequest& dlr, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
 
  return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
}
 
inline void
eprop_iaf_adapt_bsshslm_2020::get_status( Dictionary& d ) const
{
  EpropArchivingNodeRecurrent::get_status( d );
  P_.get( d );
  S_.get( d, P_ );
  d[ names::recordables ] = recordablesMap_.get_list();
}
 
inline void
eprop_iaf_adapt_bsshslm_2020::set_status( const Dictionary& d )
{
  EpropArchivingNodeRecurrent::set_status( d );
  // temporary copies in case of errors
  Parameters_ ptmp = P_;
  State_ stmp = S_;
 
  // make sure that ptmp and stmp consistent - throw BadProperty if not
  const double delta_EL = ptmp.set( d, this );
  stmp.set( d, ptmp, delta_EL, this );
 
  P_ = ptmp;
  S_ = stmp;
}
 
}  // namespace nest
 
#endif  // EPROP_IAF_ADAPT_BSSHSLM_2020_H