eprop_readout_bsshslm_2020.h

Go to the documentation of this file.

/*
 *  eprop_readout_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_READOUT_BSSHSLM_2020_H
#define EPROP_READOUT_BSSHSLM_2020_H
 
// nestkernel
#include "connection.h"
#include "eprop_archiving_node_impl.h"
#include "eprop_archiving_node_readout.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, Bellec
 
Short description
+++++++++++++++++
 
Current-based leaky integrate readout neuron model with delta-shaped or exponentially filtered
postsynaptic currents for e-prop plasticity
 
Description
+++++++++++
 
``eprop_readout_bsshslm_2020`` is an implementation of an integrate-and-fire neuron model
with delta-shaped postsynaptic currents used as readout neuron 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.
 
The membrane voltage time course :math:`v_j^t` of the neuron :math:`j` is given by:
 
.. math::
  v_j^t &= \kappa v_j^{t-1} + \zeta \sum_{i \neq j} W_{ji}^\text{out} z_i^{t-1} \,, \\
  \kappa &= e^{ -\frac{ \Delta t }{ \tau_\text{m} } } \,, \\
  \zeta &=
    \begin{cases}
      1 \\
      1 - \kappa
    \end{cases} \,, \\
 
where :math:`W_{ji}^\text{out}` is the output synaptic weight matrix and
:math:`z_i^{t-1}` is the recurrent presynaptic spike state variable.
 
Descriptions of further parameters and variables can be found in the table below.
 
The spike state variable of a presynaptic neuron is expressed by a Heaviside function:
 
.. math::
  z_i^t = H \left( v_i^t - v_\text{th} \right) \,. \\
 
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}` 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{z}_i^{t-1} \,. \\
 
The presynaptic spike trains are low-pass filtered with the following exponential kernel:
 
.. math::
  \bar{z}_i^t &=\mathcal{F}_\kappa(z_i^t) \,, \\
  \mathcal{F}_\kappa(z_i^t) &= \kappa \mathcal{F}_\kappa \left( z_i^{t-1} \right) + z_i^t \,, \\
  \mathcal{F}_\kappa(z_i^0) &= z_i^0 \,. \\
 
Since readout neurons are leaky integrators without a spiking mechanism, the
formula for computing the gradient lacks the surrogate gradient /
pseudo-derivative and a firing regularization term.
 
The learning signal :math:`L_j^t` is given by the non-plastic feedback weight
matrix :math:`B_{jk}` and the continuous error signal :math:`e_k^t` emitted by
readout neuron :math:`k`:
 
.. math::
  L_j^t = B_{jk} e_k^t \,. \\
 
The error signal depends on the selected loss function.
If a mean squared error loss is selected, then:
 
.. math::
  e_k^t = y_k^t - y_k^{*,t} \,, \\
 
where the readout signal :math:`y_k^t` corresponds to the membrane voltage of
readout neuron :math:`k` and :math:`y_k^{*,t}` is the real-valued target signal.
 
If a cross-entropy loss is selected, then:
 
.. math::
  e^k_t &= \pi_k^t - \pi_k^{*,t} \,, \\
  \pi_k^t &= \text{softmax}_k \left( y_1^t, ..., y_K^t \right) =
    \frac{ \exp \left( y_k^t\right) }{ \sum_{k'} \exp \left( y_{k'}^t \right) } \,, \\
 
where the readout signal :math:`\pi_k^t` corresponds to the softmax of the
membrane voltage of readout neuron :math:`k` and :math:`\pi_k^{*,t}` is the
one-hot encoded target signal.
 
Furthermore, the readout and target signal are zero before the onset of the
learning window in each update interval.
 
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_iaf_adapt_bsshslm_2020<../models/eprop_iaf_adapt_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
========================= ======= ===================== ================== =====================================
``C_m``                   pF      :math:`C_\text{m}`                 250.0 Capacitance of the membrane
``E_L``                   mV      :math:`E_\text{L}`                   0.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)
``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 membrane
                                                        value              voltage
                                                        representable by
                                                        ``double`` type in
                                                        C++
========================= ======= ===================== ================== =====================================
 
========== ======= ===================== ==================== =========================================
**E-prop parameters**
-------------------------------------------------------------------------------------------------------
Parameter  Unit    Math equivalent       Default              Description
========== ======= ===================== ==================== =========================================
``loss``           :math:`E`             "mean_squared_error" Loss function
                                                              ["mean_squared_error", "cross_entropy"]
========== ======= ===================== ==================== =========================================
 
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
=============== ==== =============== ============= ================
``V_m``         mV   :math:`v_j`               0.0 Membrane voltage
=============== ==== =============== ============= ================
 
========================= ==== =============== ============= ===============================
**E-prop state variables and recordables**
--------------------------------------------------------------------------------------------
State variable            Unit Math equivalent Initial value Description
========================= ==== =============== ============= ===============================
``error_signal``          mV   :math:`L_j`               0.0 Error signal
``readout_signal``        mV   :math:`y_j`               0.0 Readout signal
``readout_signal_unnorm`` mV                             0.0 Unnormalized readout signal
``target_signal``         mV   :math:`y^*_j`             0.0 Target signal
========================= ==== =============== ============= ===============================
 
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
+++++
 
LearningSignalConnectionEvent, DelayedRateConnectionEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DelayedRateConnectionEvent, DataLoggingRequest
 
See also
++++++++
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: eprop_readout_bsshslm_2020
 
EndUserDocs */
 
void register_eprop_readout_bsshslm_2020( const std::string& name );
 
class eprop_readout_bsshslm_2020 : public EpropArchivingNodeReadout< true >
{
 
public:
  eprop_readout_bsshslm_2020();
 
  eprop_readout_bsshslm_2020( const eprop_readout_bsshslm_2020& );
 
  using Node::handle;
  using Node::handles_test_event;
 
  using Node::sends_secondary_event;
 
  void
  sends_secondary_event( LearningSignalConnectionEvent& ) override
  {
  }
 
  void
  sends_secondary_event( DelayedRateConnectionEvent& ) override
  {
  }
 
  void handle( SpikeEvent& ) override;
  void handle( CurrentEvent& ) override;
  void handle( DelayedRateConnectionEvent& ) 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( DelayedRateConnectionEvent&, 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 >& presyn_isis,
    const long t_previous_update,
    const long t_previous_trigger_spike,
    const double kappa,
    const bool average_gradient ) override;
 
  long get_shift() const override;
  bool is_eprop_recurrent_node() const override;
 
  void compute_error_signal_mean_squared_error();
 
  void compute_error_signal_cross_entropy();
 
  void ( eprop_readout_bsshslm_2020::*compute_error_signal )();
 
  friend class RecordablesMap< eprop_readout_bsshslm_2020 >;
 
  friend class UniversalDataLogger< eprop_readout_bsshslm_2020 >;
 
  struct Parameters_
  {
    double C_m_;
 
    double E_L_;
 
    double I_e_;
 
    std::string loss_;
 
    bool regular_spike_arrival_;
 
    double tau_m_;
 
    double V_min_;
 
    Parameters_();
 
    void get( Dictionary& ) const;
 
    double set( const Dictionary&, Node* );
  };
 
  struct State_
  {
    double error_signal_;
 
    double readout_signal_;
 
    double readout_signal_unnorm_;
 
    double target_signal_;
 
    double i_in_;
 
    double v_m_;
 
    double z_in_;
 
    State_();
 
    void get( Dictionary&, const Parameters_& ) const;
 
    void set( const Dictionary&, const Parameters_&, double, Node* );
  };
 
  struct Buffers_
  {
    Buffers_( eprop_readout_bsshslm_2020& );
 
    Buffers_( const Buffers_&, eprop_readout_bsshslm_2020& );
 
    double normalization_rate_;
 
    RingBuffer spikes_;
 
    RingBuffer currents_;
 
    UniversalDataLogger< eprop_readout_bsshslm_2020 > logger_;
  };
 
  struct Variables_
  {
    double P_v_m_;
 
    double P_z_in_;
 
    double P_i_in_;
 
    bool signal_to_other_readouts_;
  };
 
  static const size_t MIN_RATE_RECEPTOR = 1;
 
  enum RateSynapseTypes
  {
    READOUT_SIG = MIN_RATE_RECEPTOR,
    TARGET_SIG,
    SUP_RATE_RECEPTOR
  };
 
  double
  get_v_m_() const
  {
    return S_.v_m_ + P_.E_L_;
  }
 
  double
  get_readout_signal_() const
  {
    return S_.readout_signal_;
  }
 
  double
  get_readout_signal_unnorm_() const
  {
    return S_.readout_signal_unnorm_;
  }
 
  double
  get_target_signal_() const
  {
    return S_.target_signal_;
  }
 
  double
  get_error_signal_() const
  {
    return S_.error_signal_;
  }
 
  // the order in which the structure instances are defined is important for speed
 
  Parameters_ P_;
 
  State_ S_;
 
  Variables_ V_;
 
  Buffers_ B_;
 
  static RecordablesMap< eprop_readout_bsshslm_2020 > recordablesMap_;
};
 
inline long
eprop_readout_bsshslm_2020::get_shift() const
{
  return offset_gen_ + delay_in_rec_ + delay_rec_out_;
}
 
inline bool
eprop_readout_bsshslm_2020::is_eprop_recurrent_node() const
{
  return false;
}
 
inline size_t
eprop_readout_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_readout_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_readout_bsshslm_2020::handles_test_event( DelayedRateConnectionEvent& e, size_t receptor_type )
{
  size_t step_rate_model_id = kernel().model_manager.get_node_model_id( "step_rate_generator" );
  size_t model_id = e.get_sender().get_model_id();
 
  if ( step_rate_model_id == model_id and receptor_type != TARGET_SIG )
  {
    throw IllegalConnection(
      "eprop_readout_bsshslm_2020 neurons expect a connection with a step_rate_generator node through receptor_type "
      "2." );
  }
 
  if ( receptor_type < MIN_RATE_RECEPTOR or receptor_type >= SUP_RATE_RECEPTOR )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
 
  return receptor_type;
}
 
inline size_t
eprop_readout_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_readout_bsshslm_2020::get_status( Dictionary& d ) const
{
  EpropArchivingNodeReadout::get_status( d );
  P_.get( d );
  S_.get( d, P_ );
  d[ names::recordables ] = recordablesMap_.get_list();
 
  Dictionary receptor_dict_;
  receptor_dict_[ names::readout_signal ] = static_cast< long >( READOUT_SIG );
  receptor_dict_[ names::target_signal ] = static_cast< long >( TARGET_SIG );
 
  d[ names::receptor_types ] = receptor_dict_;
}
 
inline void
eprop_readout_bsshslm_2020::set_status( const Dictionary& d )
{
  EpropArchivingNodeReadout::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_READOUT_BSSHSLM_2020_H