diff --git a/libnestutil/nest_types.h b/libnestutil/nest_types.h index ad55ddfb9b..24dca07381 100644 --- a/libnestutil/nest_types.h +++ b/libnestutil/nest_types.h @@ -90,6 +90,7 @@ constexpr uint8_t NUM_BITS_SYN_ID = 6U; 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_MARKER_ACTIVATION = 1U; constexpr uint8_t NUM_BITS_LAG = 14U; constexpr uint8_t NUM_BITS_DELAY = 21U; constexpr uint8_t NUM_BITS_NODE_ID = 61U; diff --git a/models/eprop_iaf.cpp b/models/eprop_iaf.cpp index fb336a40f4..4a2d2ac8d1 100644 --- a/models/eprop_iaf.cpp +++ b/models/eprop_iaf.cpp @@ -85,6 +85,7 @@ eprop_iaf::Parameters_::Parameters_() , kappa_( 0.97 ) , kappa_reg_( 0.97 ) , eprop_isi_trace_cutoff_( 1000.0 ) + , activation_interval_( 3000.0 ) { } @@ -131,6 +132,7 @@ eprop_iaf::Parameters_::get( DictionaryDatum& d ) const def< double >( d, names::kappa, kappa_ ); def< double >( d, names::kappa_reg, kappa_reg_ ); def< double >( d, names::eprop_isi_trace_cutoff, eprop_isi_trace_cutoff_ ); + def< double >( d, names::activation_interval, activation_interval_ ); } double @@ -168,6 +170,7 @@ eprop_iaf::Parameters_::set( const DictionaryDatum& d, Node* node ) updateValueParam< double >( d, names::kappa, kappa_, node ); updateValueParam< double >( d, names::kappa_reg, kappa_reg_, node ); updateValueParam< double >( d, names::eprop_isi_trace_cutoff, eprop_isi_trace_cutoff_, node ); + updateValueParam< double >( d, names::activation_interval, activation_interval_, node ); if ( C_m_ <= 0 ) { @@ -209,11 +212,16 @@ eprop_iaf::Parameters_::set( const DictionaryDatum& d, Node* node ) throw BadProperty( "Firing rate low-pass filter for regularization kappa_reg from range [0, 1] required." ); } - if ( eprop_isi_trace_cutoff_ < 0.0 ) + if ( activation_interval_ < 0 ) { - throw BadProperty( "Cutoff of integration of eprop trace between spikes eprop_isi_trace_cutoff ≥ 0 required." ); + throw BadProperty( "Interval between activations activation_interval ≥ 0 required." ); } + if ( eprop_isi_trace_cutoff_ < 0.0 or eprop_isi_trace_cutoff_ > activation_interval_ ) + { + throw BadProperty( + "Computation cutoff of eprop trace 0 ≤ eprop trace eprop_isi_trace_cutoff ≤ activation_interval required." ); + } return delta_EL; } @@ -271,6 +279,7 @@ eprop_iaf::pre_run_hook() V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); V_.eprop_isi_trace_cutoff_steps_ = Time( Time::ms( P_.eprop_isi_trace_cutoff_ ) ).get_steps(); + V_.activation_interval_steps_ = Time( Time::ms( P_.activation_interval_ ) ).get_steps(); // calculate the entries of the propagator matrix for the evolution of the state vector @@ -314,6 +323,14 @@ eprop_iaf::update( Time const& origin, const long from, const long to ) S_.z_ = 1.0; S_.v_m_ -= P_.V_th_ * S_.z_; S_.r_ = V_.RefractoryCounts_; + set_last_event_time( t ); + } + else if ( get_last_event_time() > 0 and t - get_last_event_time() >= V_.activation_interval_steps_ ) + { + SpikeEvent se; + se.set_activation(); + kernel().event_delivery_manager.send( *this, se, lag ); + set_last_event_time( t ); } append_new_eprop_history_entry( t ); @@ -380,61 +397,70 @@ eprop_iaf::compute_gradient( const long t_spike, double& epsilon, double& weight, const CommonSynapseProperties& cp, - WeightOptimizer* optimizer ) + WeightOptimizer* optimizer, + const bool activation, + const bool previous_event_was_activation, + double& sum_grad ) { - double e = 0.0; // eligibility trace - double z = 0.0; // spiking variable - double z_current_buffer = 1.0; // buffer containing the spike that triggered the current integration - double psi = 0.0; // surrogate gradient - double L = 0.0; // learning signal - double firing_rate_reg = 0.0; // firing rate regularization - double grad = 0.0; // gradient + const auto& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); + const auto& opt_cp = *ecp.optimizer_cp_; + const bool optimize_each_step = opt_cp.optimize_each_step_; - const EpropSynapseCommonProperties& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); - const auto optimize_each_step = ( *ecp.optimizer_cp_ ).optimize_each_step_; + if ( not previous_event_was_activation ) + { + sum_grad = 0.0; // sum of gradients + } auto eprop_hist_it = get_eprop_history( t_spike_previous - 1 ); - const long t_compute_until = std::min( t_spike_previous + V_.eprop_isi_trace_cutoff_steps_, t_spike ); + const long cutoff_end = t_spike_previous + V_.eprop_isi_trace_cutoff_steps_; + const long t_compute_until = std::min( cutoff_end, t_spike ); - for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) + if ( not previous_event_was_activation ) { - z = z_previous_buffer; - z_previous_buffer = z_current_buffer; - z_current_buffer = 0.0; - - psi = eprop_hist_it->surrogate_gradient_; - L = eprop_hist_it->learning_signal_; - firing_rate_reg = eprop_hist_it->firing_rate_reg_; + double z_current_buffer = 1.0; // spike that triggered current computation - z_bar = V_.P_v_m_ * z_bar + z; - e = psi * z_bar; - e_bar = P_.kappa_ * e_bar + e; - e_bar_reg = P_.kappa_reg_ * e_bar_reg + ( 1.0 - P_.kappa_reg_ ) * e; - - if ( optimize_each_step ) + for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) { - grad = L * e_bar + firing_rate_reg * e_bar_reg; - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t, grad, weight ); - } - else - { - grad += L * e_bar + firing_rate_reg * e_bar_reg; + const double z = z_previous_buffer; // spiking variable + z_previous_buffer = z_current_buffer; + z_current_buffer = 0.0; + + const double psi = eprop_hist_it->surrogate_gradient_; // surrogate gradient + const double L = eprop_hist_it->learning_signal_; // learning signal + const double firing_rate_reg = eprop_hist_it->firing_rate_reg_; // firing rate regularization + + z_bar = V_.P_v_m_ * z_bar + z; + const double e = psi * z_bar; // eligibility trace + e_bar = P_.kappa_ * e_bar + e; + e_bar_reg = P_.kappa_reg_ * e_bar_reg + ( 1.0 - P_.kappa_reg_ ) * e; + + const double grad = L * e_bar + firing_rate_reg * e_bar_reg; + + if ( optimize_each_step ) + { + sum_grad = grad; + weight = optimizer->optimized_weight( opt_cp, t, sum_grad, weight ); + } + else + { + sum_grad += grad; + } } } - if ( not optimize_each_step ) + const long trace_decay_interval = t_spike - ( previous_event_was_activation ? t_spike_previous : t_compute_until ); + + if ( trace_decay_interval > 0 ) { - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t_compute_until, grad, weight ); + z_bar *= std::pow( V_.P_v_m_, trace_decay_interval ); + e_bar *= std::pow( P_.kappa_, trace_decay_interval ); + e_bar_reg *= std::pow( P_.kappa_reg_, trace_decay_interval ); } - const long cutoff_to_spike_interval = t_spike - t_compute_until; - - if ( cutoff_to_spike_interval > 0 ) + if ( not( activation or optimize_each_step ) ) { - z_bar *= std::pow( V_.P_v_m_, cutoff_to_spike_interval ); - e_bar *= std::pow( P_.kappa_, cutoff_to_spike_interval ); - e_bar_reg *= std::pow( P_.kappa_reg_, cutoff_to_spike_interval ); + weight = optimizer->optimized_weight( opt_cp, t_compute_until, sum_grad, weight ); } } diff --git a/models/eprop_iaf.h b/models/eprop_iaf.h index 225de65333..6a7aa0e149 100644 --- a/models/eprop_iaf.h +++ b/models/eprop_iaf.h @@ -230,6 +230,8 @@ Parameter Unit Math equivalent Default Des ---------------------------------------------------------------------------------------------------------------- Parameter Unit Math equivalent Default Description =============================== ======= =========================== ================== ========================= +``activation_interval`` ms 3000.0 Interval between two + activations ``c_reg`` :math:`c_\text{reg}` 0.0 Coefficient of firing rate regularization ``eprop_isi_trace_cutoff`` ms :math:`{\Delta t}_\text{c}` maximum value Cutoff for integration of @@ -397,7 +399,10 @@ class eprop_iaf : public EpropArchivingNodeRecurrent< false > double&, double&, const CommonSynapseProperties&, - WeightOptimizer* ) override; + WeightOptimizer*, + const bool, + const bool, + double& ) override; long get_shift() const override; bool is_eprop_recurrent_node() const override; @@ -458,6 +463,9 @@ class eprop_iaf : public EpropArchivingNodeRecurrent< false > //! Time interval from the previous spike until the cutoff of e-prop update integration between two spikes (ms). double eprop_isi_trace_cutoff_; + //! Interval between two activations. + long activation_interval_; + //! Default constructor. Parameters_(); @@ -535,6 +543,9 @@ class eprop_iaf : public EpropArchivingNodeRecurrent< false > //! Time steps from the previous spike until the cutoff of e-prop update integration between two spikes. long eprop_isi_trace_cutoff_steps_; + + //! Time steps of activation interval. + long activation_interval_steps_; }; //! Get the current value of the membrane voltage. diff --git a/models/eprop_iaf_adapt.cpp b/models/eprop_iaf_adapt.cpp index 10cb9ff224..7a3a617645 100644 --- a/models/eprop_iaf_adapt.cpp +++ b/models/eprop_iaf_adapt.cpp @@ -89,6 +89,7 @@ eprop_iaf_adapt::Parameters_::Parameters_() , kappa_( 0.97 ) , kappa_reg_( 0.97 ) , eprop_isi_trace_cutoff_( 1000.0 ) + , activation_interval_( 3000.0 ) { } @@ -139,6 +140,7 @@ eprop_iaf_adapt::Parameters_::get( DictionaryDatum& d ) const def< double >( d, names::kappa, kappa_ ); def< double >( d, names::kappa_reg, kappa_reg_ ); def< double >( d, names::eprop_isi_trace_cutoff, eprop_isi_trace_cutoff_ ); + def< double >( d, names::activation_interval, activation_interval_ ); } double @@ -178,6 +180,7 @@ eprop_iaf_adapt::Parameters_::set( const DictionaryDatum& d, Node* node ) updateValueParam< double >( d, names::kappa, kappa_, node ); updateValueParam< double >( d, names::kappa_reg, kappa_reg_, node ); updateValueParam< double >( d, names::eprop_isi_trace_cutoff, eprop_isi_trace_cutoff_, node ); + updateValueParam< double >( d, names::activation_interval, activation_interval_, node ); if ( adapt_beta_ < 0 ) { @@ -229,11 +232,16 @@ eprop_iaf_adapt::Parameters_::set( const DictionaryDatum& d, Node* node ) throw BadProperty( "Firing rate low-pass filter for regularization kappa_reg from range [0, 1] required." ); } - if ( eprop_isi_trace_cutoff_ < 0.0 ) + if ( activation_interval_ < 0 ) { - throw BadProperty( "Cutoff of integration of eprop trace between spikes eprop_isi_trace_cutoff ≥ 0 required." ); + throw BadProperty( "Interval between activations activation_interval ≥ 0 required." ); } + if ( eprop_isi_trace_cutoff_ < 0.0 or eprop_isi_trace_cutoff_ > activation_interval_ ) + { + throw BadProperty( + "Computation cutoff of eprop trace 0 ≤ eprop trace eprop_isi_trace_cutoff ≤ activation_interval required." ); + } return delta_EL; } @@ -305,6 +313,7 @@ eprop_iaf_adapt::pre_run_hook() V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); V_.eprop_isi_trace_cutoff_steps_ = Time( Time::ms( P_.eprop_isi_trace_cutoff_ ) ).get_steps(); + V_.activation_interval_steps_ = Time( Time::ms( P_.activation_interval_ ) ).get_steps(); // calculate the entries of the propagator matrix for the evolution of the state vector @@ -353,6 +362,14 @@ eprop_iaf_adapt::update( Time const& origin, const long from, const long to ) S_.z_ = 1.0; S_.v_m_ -= P_.V_th_ * S_.z_; S_.r_ = V_.RefractoryCounts_; + set_last_event_time( t ); + } + else if ( get_last_event_time() > 0 and t - get_last_event_time() >= V_.activation_interval_steps_ ) + { + SpikeEvent se; + se.set_activation(); + kernel().event_delivery_manager.send( *this, se, lag ); + set_last_event_time( t ); } append_new_eprop_history_entry( t ); @@ -419,63 +436,72 @@ eprop_iaf_adapt::compute_gradient( const long t_spike, double& epsilon, double& weight, const CommonSynapseProperties& cp, - WeightOptimizer* optimizer ) + WeightOptimizer* optimizer, + const bool activation, + const bool previous_event_was_activation, + double& sum_grad ) { - double e = 0.0; // eligibility trace - double z = 0.0; // spiking variable - double z_current_buffer = 1.0; // buffer containing the spike that triggered the current integration - double psi = 0.0; // surrogate gradient - double L = 0.0; // learning signal - double firing_rate_reg = 0.0; // firing rate regularization - double grad = 0.0; // gradient + const auto& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); + const auto& opt_cp = *ecp.optimizer_cp_; + const bool optimize_each_step = opt_cp.optimize_each_step_; - const EpropSynapseCommonProperties& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); - const auto optimize_each_step = ( *ecp.optimizer_cp_ ).optimize_each_step_; + if ( not previous_event_was_activation ) + { + sum_grad = 0.0; // sum of gradients + } auto eprop_hist_it = get_eprop_history( t_spike_previous - 1 ); - const long t_compute_until = std::min( t_spike_previous + V_.eprop_isi_trace_cutoff_steps_, t_spike ); + const long cutoff_end = t_spike_previous + V_.eprop_isi_trace_cutoff_steps_; + const long t_compute_until = std::min( cutoff_end, t_spike ); - for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) + if ( not previous_event_was_activation ) { - z = z_previous_buffer; - z_previous_buffer = z_current_buffer; - z_current_buffer = 0.0; - - psi = eprop_hist_it->surrogate_gradient_; - L = eprop_hist_it->learning_signal_; - firing_rate_reg = eprop_hist_it->firing_rate_reg_; + double z_current_buffer = 1.0; // spike that triggered current computation - z_bar = V_.P_v_m_ * z_bar + z; - e = psi * ( z_bar - P_.adapt_beta_ * epsilon ); - epsilon = V_.P_adapt_ * epsilon + e; - e_bar = P_.kappa_ * e_bar + e; - e_bar_reg = P_.kappa_reg_ * e_bar_reg + ( 1.0 - P_.kappa_reg_ ) * e; - - if ( optimize_each_step ) + for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) { - grad = L * e_bar + firing_rate_reg * e_bar_reg; - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t, grad, weight ); - } - else - { - grad += L * e_bar + firing_rate_reg * e_bar_reg; + const double z = z_previous_buffer; // spiking variable + z_previous_buffer = z_current_buffer; + z_current_buffer = 0.0; + + const double psi = eprop_hist_it->surrogate_gradient_; // surrogate gradient + const double L = eprop_hist_it->learning_signal_; // learning signal + const double firing_rate_reg = eprop_hist_it->firing_rate_reg_; // firing rate regularization + + z_bar = V_.P_v_m_ * z_bar + z; + const double e = psi * ( z_bar - P_.adapt_beta_ * epsilon ); // eligibility trace + epsilon = V_.P_adapt_ * epsilon + e; + e_bar = P_.kappa_ * e_bar + e; + e_bar_reg = P_.kappa_reg_ * e_bar_reg + ( 1.0 - P_.kappa_reg_ ) * e; + + const double grad = L * e_bar + firing_rate_reg * e_bar_reg; + + if ( optimize_each_step ) + { + sum_grad = grad; + weight = optimizer->optimized_weight( opt_cp, t, sum_grad, weight ); + } + else + { + sum_grad += grad; + } } } - if ( not optimize_each_step ) + const long trace_decay_interval = t_spike - ( previous_event_was_activation ? t_spike_previous : t_compute_until ); + + if ( trace_decay_interval > 0 ) { - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t_compute_until, grad, weight ); + z_bar *= std::pow( V_.P_v_m_, trace_decay_interval ); + e_bar *= std::pow( P_.kappa_, trace_decay_interval ); + e_bar_reg *= std::pow( P_.kappa_reg_, trace_decay_interval ); + epsilon *= std::pow( V_.P_adapt_, trace_decay_interval ); } - const long cutoff_to_spike_interval = t_spike - t_compute_until; - - if ( cutoff_to_spike_interval > 0 ) + if ( not( activation or optimize_each_step ) ) { - z_bar *= std::pow( V_.P_v_m_, cutoff_to_spike_interval ); - e_bar *= std::pow( P_.kappa_, cutoff_to_spike_interval ); - e_bar_reg *= std::pow( P_.kappa_reg_, cutoff_to_spike_interval ); - epsilon *= std::pow( V_.P_adapt_, cutoff_to_spike_interval ); + weight = optimizer->optimized_weight( opt_cp, t_compute_until, sum_grad, weight ); } } diff --git a/models/eprop_iaf_adapt.h b/models/eprop_iaf_adapt.h index d404dcbd69..3799cbf3bc 100644 --- a/models/eprop_iaf_adapt.h +++ b/models/eprop_iaf_adapt.h @@ -213,6 +213,8 @@ Parameter Unit Math equivalent Default Des ---------------------------------------------------------------------------------------------------------------- Parameter Unit Math equivalent Default Description =============================== ======= =========================== ================== ========================= +``activation_interval`` ms 3000.0 Interval between two + activations ``c_reg`` :math:`c_\text{reg}` 0.0 Coefficient of firing rate regularization ``eprop_isi_trace_cutoff`` ms :math:`{\Delta t}_\text{c}` maximum value Cutoff for integration of @@ -365,7 +367,10 @@ class eprop_iaf_adapt : public EpropArchivingNodeRecurrent< false > double&, double&, const CommonSynapseProperties&, - WeightOptimizer* ) override; + WeightOptimizer*, + const bool, + const bool, + double& ) override; long get_shift() const override; bool is_eprop_recurrent_node() const override; @@ -432,6 +437,9 @@ class eprop_iaf_adapt : public EpropArchivingNodeRecurrent< false > //! Time interval from the previous spike until the cutoff of e-prop update integration between two spikes (ms). double eprop_isi_trace_cutoff_; + //! Interval between two activations. + long activation_interval_; + //! Default constructor. Parameters_(); @@ -518,6 +526,9 @@ class eprop_iaf_adapt : public EpropArchivingNodeRecurrent< false > //! Time steps from the previous spike until the cutoff of e-prop update integration between two spikes. long eprop_isi_trace_cutoff_steps_; + + //! Time steps of activation interval. + long activation_interval_steps_; }; //! Get the current value of the membrane voltage. diff --git a/models/eprop_iaf_adapt_bsshslm_2020.cpp b/models/eprop_iaf_adapt_bsshslm_2020.cpp index f37f2203a0..95d68c06e5 100644 --- a/models/eprop_iaf_adapt_bsshslm_2020.cpp +++ b/models/eprop_iaf_adapt_bsshslm_2020.cpp @@ -87,6 +87,7 @@ eprop_iaf_adapt_bsshslm_2020::Parameters_::Parameters_() , tau_m_( 10.0 ) , V_min_( -std::numeric_limits< double >::max() ) , V_th_( -55.0 - E_L_ ) + , activation_interval_( 3 * kernel().simulation_manager.get_eprop_update_interval().get_ms() ) { } @@ -135,6 +136,7 @@ eprop_iaf_adapt_bsshslm_2020::Parameters_::get( DictionaryDatum& d ) const def< double >( d, names::tau_m, tau_m_ ); def< double >( d, names::V_min, V_min_ + E_L_ ); def< double >( d, names::V_th, V_th_ + E_L_ ); + def< double >( d, names::activation_interval, activation_interval_ ); } double @@ -162,6 +164,7 @@ eprop_iaf_adapt_bsshslm_2020::Parameters_::set( const DictionaryDatum& d, Node* updateValueParam< double >( d, names::gamma, gamma_, node ); updateValueParam< double >( d, names::I_e, I_e_, node ); updateValueParam< bool >( d, names::regular_spike_arrival, regular_spike_arrival_, node ); + updateValueParam< double >( d, names::activation_interval, activation_interval_, node ); if ( updateValueParam< std::string >( d, names::surrogate_gradient_function, surrogate_gradient_function_, node ) ) { @@ -213,6 +216,19 @@ eprop_iaf_adapt_bsshslm_2020::Parameters_::set( const DictionaryDatum& d, Node* throw BadProperty( "Spike threshold voltage V_th ≥ minimal voltage V_min required." ); } + if ( activation_interval_ < 0 ) + { + throw BadProperty( "Interval between activations activation_interval ≥ 0 required." ); + } + + if ( Time( Time::ms( activation_interval_ ) ).get_steps() + % kernel().simulation_manager.get_eprop_update_interval().get_steps() + != 0 ) + { + throw BadProperty( + "Interval between activations activation_interval required to be an integer multiple of the e-prop update " + "interval." ); + } return delta_EL; } @@ -283,6 +299,7 @@ eprop_iaf_adapt_bsshslm_2020::pre_run_hook() B_.logger_.init(); // ensures initialization in case multimeter connected after Simulate V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); + V_.activation_interval_steps_ = Time( Time::ms( P_.activation_interval_ ) ).get_steps(); // calculate the entries of the propagator matrix for the evolution of the state vector @@ -353,6 +370,14 @@ eprop_iaf_adapt_bsshslm_2020::update( Time const& origin, const long from, const S_.z_ = 1.0; S_.r_ = V_.RefractoryCounts_; + set_last_event_time( t ); + } + else if ( get_last_event_time() > 0 and t - get_last_event_time() >= V_.activation_interval_steps_ ) + { + SpikeEvent se; + se.set_activation(); + kernel().event_delivery_manager.send( *this, se, lag ); + set_last_event_time( t ); } append_new_eprop_history_entry( t ); @@ -423,43 +448,39 @@ eprop_iaf_adapt_bsshslm_2020::compute_gradient( std::vector< long >& presyn_isis { auto eprop_hist_it = get_eprop_history( t_previous_trigger_spike ); - double e = 0.0; // eligibility trace double e_bar = 0.0; // low-pass filtered eligibility trace double epsilon = 0.0; // adaptive component of eligibility vector double grad = 0.0; // gradient value to be calculated - double L = 0.0; // learning signal - double psi = 0.0; // surrogate gradient double sum_e = 0.0; // sum of eligibility traces - double z = 0.0; // spiking variable double z_bar = 0.0; // low-pass filtered spiking variable - for ( long presyn_isi : presyn_isis ) + for ( const long presyn_isi : presyn_isis ) { - z = 1.0; // set spiking variable to 1 for each incoming spike + double z = 1.0; // set spiking variable to 1 for each incoming spike - for ( long t = 0; t < presyn_isi; ++t ) + for ( long t = 0; t < presyn_isi; ++t, ++eprop_hist_it ) { assert( eprop_hist_it != eprop_history_.end() ); - psi = eprop_hist_it->surrogate_gradient_; - L = eprop_hist_it->learning_signal_; + const double psi = eprop_hist_it->surrogate_gradient_; // surrogate gradient + const double L = eprop_hist_it->learning_signal_; // learning signal z_bar = V_.P_v_m_ * z_bar + V_.P_z_in_ * z; - e = psi * ( z_bar - P_.adapt_beta_ * epsilon ); + const double e = psi * ( z_bar - P_.adapt_beta_ * epsilon ); // eligibility trace epsilon = V_.P_adapt_ * epsilon + e; e_bar = kappa * e_bar + ( 1.0 - kappa ) * e; + grad += L * e_bar; sum_e += e; - z = 0.0; // set spiking variable to 0 between spikes - ++eprop_hist_it; + z = 0.0; // set spiking variable to 0 between spikes } } presyn_isis.clear(); const long update_interval = kernel().simulation_manager.get_eprop_update_interval().get_steps(); const long learning_window = kernel().simulation_manager.get_eprop_learning_window().get_steps(); - const auto firing_rate_reg = get_firing_rate_reg_history( t_previous_update + get_shift() + update_interval ); + const double firing_rate_reg = get_firing_rate_reg_history( t_previous_update + get_shift() + update_interval ); grad += firing_rate_reg * sum_e; diff --git a/models/eprop_iaf_adapt_bsshslm_2020.h b/models/eprop_iaf_adapt_bsshslm_2020.h index f8d8009ba2..3426bf082c 100644 --- a/models/eprop_iaf_adapt_bsshslm_2020.h +++ b/models/eprop_iaf_adapt_bsshslm_2020.h @@ -196,6 +196,8 @@ Parameter Unit Math equivalent Default Des ---------------------------------------------------------------------------------------------------------------- Parameter Unit Math equivalent Default Description =============================== ======= ======================= ================== ============================= +``activation_interval`` ms 3000.0 Interval between two + activations ``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 @@ -386,6 +388,9 @@ class eprop_iaf_adapt_bsshslm_2020 : public EpropArchivingNodeRecurrent< true > //! Spike threshold voltage relative to the leak membrane potential (mV). double V_th_; + //! Interval between two activations. + long activation_interval_; + //! Default constructor. Parameters_(); @@ -473,6 +478,9 @@ class eprop_iaf_adapt_bsshslm_2020 : public EpropArchivingNodeRecurrent< true > //! Total refractory steps. int RefractoryCounts_; + + //! Time steps of activation interval. + long activation_interval_steps_; }; //! Get the current value of the membrane voltage. diff --git a/models/eprop_iaf_bsshslm_2020.cpp b/models/eprop_iaf_bsshslm_2020.cpp index a4c52557ff..5959194f54 100644 --- a/models/eprop_iaf_bsshslm_2020.cpp +++ b/models/eprop_iaf_bsshslm_2020.cpp @@ -83,6 +83,7 @@ eprop_iaf_bsshslm_2020::Parameters_::Parameters_() , tau_m_( 10.0 ) , V_min_( -std::numeric_limits< double >::max() ) , V_th_( -55.0 - E_L_ ) + , activation_interval_( 3 * kernel().simulation_manager.get_eprop_update_interval().get_ms() ) { } @@ -127,6 +128,7 @@ eprop_iaf_bsshslm_2020::Parameters_::get( DictionaryDatum& d ) const def< double >( d, names::tau_m, tau_m_ ); def< double >( d, names::V_min, V_min_ + E_L_ ); def< double >( d, names::V_th, V_th_ + E_L_ ); + def< double >( d, names::activation_interval, activation_interval_ ); } double @@ -152,6 +154,7 @@ eprop_iaf_bsshslm_2020::Parameters_::set( const DictionaryDatum& d, Node* node ) updateValueParam< double >( d, names::gamma, gamma_, node ); updateValueParam< double >( d, names::I_e, I_e_, node ); updateValueParam< bool >( d, names::regular_spike_arrival, regular_spike_arrival_, node ); + updateValueParam< double >( d, names::activation_interval, activation_interval_, node ); if ( updateValueParam< std::string >( d, names::surrogate_gradient_function, surrogate_gradient_function_, node ) ) { @@ -193,6 +196,19 @@ eprop_iaf_bsshslm_2020::Parameters_::set( const DictionaryDatum& d, Node* node ) throw BadProperty( "Spike threshold voltage V_th ≥ minimal voltage V_min required." ); } + if ( activation_interval_ < 0 ) + { + throw BadProperty( "Interval between activations activation_interval ≥ 0 required." ); + } + + if ( Time( Time::ms( activation_interval_ ) ).get_steps() + % kernel().simulation_manager.get_eprop_update_interval().get_steps() + != 0 ) + { + throw BadProperty( + "Interval between activations activation_interval required to be an integer multiple of the e-prop update " + "interval." ); + } return delta_EL; } @@ -249,6 +265,7 @@ eprop_iaf_bsshslm_2020::pre_run_hook() B_.logger_.init(); // ensures initialization in case multimeter connected after Simulate V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); + V_.activation_interval_steps_ = Time( Time::ms( P_.activation_interval_ ) ).get_steps(); // calculate the entries of the propagator matrix for the evolution of the state vector @@ -313,6 +330,14 @@ eprop_iaf_bsshslm_2020::update( Time const& origin, const long from, const long S_.z_ = 1.0; S_.r_ = V_.RefractoryCounts_; + set_last_event_time( t ); + } + else if ( get_last_event_time() > 0 and t - get_last_event_time() >= V_.activation_interval_steps_ ) + { + SpikeEvent se; + se.set_activation(); + kernel().event_delivery_manager.send( *this, se, lag ); + set_last_event_time( t ); } append_new_eprop_history_entry( t ); @@ -383,41 +408,37 @@ eprop_iaf_bsshslm_2020::compute_gradient( std::vector< long >& presyn_isis, { auto eprop_hist_it = get_eprop_history( t_previous_trigger_spike ); - double e = 0.0; // eligibility trace double e_bar = 0.0; // low-pass filtered eligibility trace double grad = 0.0; // gradient value to be calculated - double L = 0.0; // learning signal - double psi = 0.0; // surrogate gradient double sum_e = 0.0; // sum of eligibility traces - double z = 0.0; // spiking variable double z_bar = 0.0; // low-pass filtered spiking variable - for ( long presyn_isi : presyn_isis ) + for ( const long presyn_isi : presyn_isis ) { - z = 1.0; // set spiking variable to 1 for each incoming spike + double z = 1.0; // set spiking variable to 1 for each incoming spike - for ( long t = 0; t < presyn_isi; ++t ) + for ( long t = 0; t < presyn_isi; ++t, ++eprop_hist_it ) { assert( eprop_hist_it != eprop_history_.end() ); - psi = eprop_hist_it->surrogate_gradient_; - L = eprop_hist_it->learning_signal_; + const double psi = eprop_hist_it->surrogate_gradient_; // surrogate gradient + const double L = eprop_hist_it->learning_signal_; // learning signal z_bar = V_.P_v_m_ * z_bar + V_.P_z_in_ * z; - e = psi * z_bar; + const double e = psi * z_bar; // eligibility trace e_bar = kappa * e_bar + ( 1.0 - kappa ) * e; + grad += L * e_bar; sum_e += e; - z = 0.0; // set spiking variable to 0 between spikes - ++eprop_hist_it; + z = 0.0; // set spiking variable to 0 between spikes } } presyn_isis.clear(); const long update_interval = kernel().simulation_manager.get_eprop_update_interval().get_steps(); const long learning_window = kernel().simulation_manager.get_eprop_learning_window().get_steps(); - const auto firing_rate_reg = get_firing_rate_reg_history( t_previous_update + get_shift() + update_interval ); + const double firing_rate_reg = get_firing_rate_reg_history( t_previous_update + get_shift() + update_interval ); grad += firing_rate_reg * sum_e; diff --git a/models/eprop_iaf_bsshslm_2020.h b/models/eprop_iaf_bsshslm_2020.h index 34cc9a575d..2abad4464e 100644 --- a/models/eprop_iaf_bsshslm_2020.h +++ b/models/eprop_iaf_bsshslm_2020.h @@ -183,6 +183,7 @@ Parameter Unit Math equivalent Default Des ---------------------------------------------------------------------------------------------------------------- Parameter Unit Math equivalent Default Description =============================== ==== ======================= ================== ================================ +``activation_interval`` ms 3000.0 Interval between two activations ``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 @@ -365,6 +366,9 @@ class eprop_iaf_bsshslm_2020 : public EpropArchivingNodeRecurrent< true > //! Spike threshold voltage relative to the leak membrane potential (mV). double V_th_; + //! Interval between two activations. + long activation_interval_; + //! Default constructor. Parameters_(); @@ -443,6 +447,9 @@ class eprop_iaf_bsshslm_2020 : public EpropArchivingNodeRecurrent< true > //! Total refractory steps. int RefractoryCounts_; + + //! Time steps of activation interval. + long activation_interval_steps_; }; //! Get the current value of the membrane voltage. diff --git a/models/eprop_iaf_psc_delta.cpp b/models/eprop_iaf_psc_delta.cpp index a8b13ac4e4..e17555fd05 100644 --- a/models/eprop_iaf_psc_delta.cpp +++ b/models/eprop_iaf_psc_delta.cpp @@ -87,6 +87,7 @@ eprop_iaf_psc_delta::Parameters_::Parameters_() , kappa_( 0.97 ) , kappa_reg_( 0.97 ) , eprop_isi_trace_cutoff_( 1000.0 ) + , activation_interval_( 3000.0 ) { } @@ -134,6 +135,7 @@ eprop_iaf_psc_delta::Parameters_::get( DictionaryDatum& d ) const def< double >( d, names::kappa, kappa_ ); def< double >( d, names::kappa_reg, kappa_reg_ ); def< double >( d, names::eprop_isi_trace_cutoff, eprop_isi_trace_cutoff_ ); + def< double >( d, names::activation_interval, activation_interval_ ); } double @@ -174,6 +176,7 @@ eprop_iaf_psc_delta::Parameters_::set( const DictionaryDatum& d, Node* node ) updateValueParam< double >( d, names::kappa, kappa_, node ); updateValueParam< double >( d, names::kappa_reg, kappa_reg_, node ); updateValueParam< double >( d, names::eprop_isi_trace_cutoff, eprop_isi_trace_cutoff_, node ); + updateValueParam< double >( d, names::activation_interval, activation_interval_, node ); if ( V_th_ < V_min_ ) { @@ -225,11 +228,16 @@ eprop_iaf_psc_delta::Parameters_::set( const DictionaryDatum& d, Node* node ) throw BadProperty( "Firing rate low-pass filter for regularization kappa_reg from range [0, 1] required." ); } - if ( eprop_isi_trace_cutoff_ < 0.0 ) + if ( activation_interval_ < 0 ) { - throw BadProperty( "Cutoff of integration of eprop trace between spikes eprop_isi_trace_cutoff ≥ 0 required." ); + throw BadProperty( "Interval between activations activation_interval ≥ 0 required." ); } + if ( eprop_isi_trace_cutoff_ < 0.0 or eprop_isi_trace_cutoff_ > activation_interval_ ) + { + throw BadProperty( + "Computation cutoff of eprop trace 0 ≤ eprop trace eprop_isi_trace_cutoff ≤ activation_interval required." ); + } return delta_EL; } @@ -287,6 +295,7 @@ eprop_iaf_psc_delta::pre_run_hook() V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); V_.eprop_isi_trace_cutoff_steps_ = Time( Time::ms( P_.eprop_isi_trace_cutoff_ ) ).get_steps(); + V_.activation_interval_steps_ = Time( Time::ms( P_.activation_interval_ ) ).get_steps(); // calculate the entries of the propagator matrix for the evolution of the state vector @@ -347,6 +356,14 @@ eprop_iaf_psc_delta::update( Time const& origin, const long from, const long to kernel().event_delivery_manager.send( *this, se, lag ); z = 1.0; + set_last_event_time( t ); + } + else if ( get_last_event_time() > 0 and t - get_last_event_time() >= V_.activation_interval_steps_ ) + { + SpikeEvent se; + se.set_activation(); + kernel().event_delivery_manager.send( *this, se, lag ); + set_last_event_time( t ); } append_new_eprop_history_entry( t ); @@ -413,61 +430,70 @@ eprop_iaf_psc_delta::compute_gradient( const long t_spike, double& epsilon, double& weight, const CommonSynapseProperties& cp, - WeightOptimizer* optimizer ) + WeightOptimizer* optimizer, + const bool activation, + const bool previous_event_was_activation, + double& sum_grad ) { - double e = 0.0; // eligibility trace - double z = 0.0; // spiking variable - double z_current_buffer = 1.0; // buffer containing the spike that triggered the current integration - double psi = 0.0; // surrogate gradient - double L = 0.0; // learning signal - double firing_rate_reg = 0.0; // firing rate regularization - double grad = 0.0; // gradient + const auto& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); + const auto& opt_cp = *ecp.optimizer_cp_; + const bool optimize_each_step = opt_cp.optimize_each_step_; - const EpropSynapseCommonProperties& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); - const auto optimize_each_step = ( *ecp.optimizer_cp_ ).optimize_each_step_; + if ( not previous_event_was_activation ) + { + sum_grad = 0.0; // sum of gradients + } auto eprop_hist_it = get_eprop_history( t_spike_previous - 1 ); - const long t_compute_until = std::min( t_spike_previous + V_.eprop_isi_trace_cutoff_steps_, t_spike ); + const long cutoff_end = t_spike_previous + V_.eprop_isi_trace_cutoff_steps_; + const long t_compute_until = std::min( cutoff_end, t_spike ); - for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) + if ( not previous_event_was_activation ) { - z = z_previous_buffer; - z_previous_buffer = z_current_buffer; - z_current_buffer = 0.0; + double z_current_buffer = 1.0; // spike that triggered current computation - psi = eprop_hist_it->surrogate_gradient_; - L = eprop_hist_it->learning_signal_; - firing_rate_reg = eprop_hist_it->firing_rate_reg_; + for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) + { + const double z = z_previous_buffer; // spiking variable + z_previous_buffer = z_current_buffer; + z_current_buffer = 0.0; - z_bar = V_.P_v_m_ * z_bar + z; - e = psi * z_bar; - e_bar = P_.kappa_ * e_bar + e; - e_bar_reg = P_.kappa_reg_ * e_bar_reg + ( 1.0 - P_.kappa_reg_ ) * e; + const double psi = eprop_hist_it->surrogate_gradient_; // surrogate gradient + const double L = eprop_hist_it->learning_signal_; // learning signal + const double firing_rate_reg = eprop_hist_it->firing_rate_reg_; // firing rate regularization - if ( optimize_each_step ) - { - grad = L * e_bar + firing_rate_reg * e_bar_reg; - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t, grad, weight ); - } - else - { - grad += L * e_bar + firing_rate_reg * e_bar_reg; + z_bar = V_.P_v_m_ * z_bar + z; + const double e = psi * z_bar; // eligibility trace + e_bar = P_.kappa_ * e_bar + e; + e_bar_reg = P_.kappa_reg_ * e_bar_reg + ( 1.0 - P_.kappa_reg_ ) * e; + + const double grad = L * e_bar + firing_rate_reg * e_bar_reg; + + if ( optimize_each_step ) + { + sum_grad = grad; + weight = optimizer->optimized_weight( opt_cp, t, sum_grad, weight ); + } + else + { + sum_grad += grad; + } } } - if ( not optimize_each_step ) + const long trace_decay_interval = t_spike - ( previous_event_was_activation ? t_spike_previous : t_compute_until ); + + if ( trace_decay_interval > 0 ) { - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t_compute_until, grad, weight ); + z_bar *= std::pow( V_.P_v_m_, trace_decay_interval ); + e_bar *= std::pow( P_.kappa_, trace_decay_interval ); + e_bar_reg *= std::pow( P_.kappa_reg_, trace_decay_interval ); } - const long cutoff_to_spike_interval = t_spike - t_compute_until; - - if ( cutoff_to_spike_interval > 0 ) + if ( not( activation or optimize_each_step ) ) { - z_bar *= std::pow( V_.P_v_m_, cutoff_to_spike_interval ); - e_bar *= std::pow( P_.kappa_, cutoff_to_spike_interval ); - e_bar_reg *= std::pow( P_.kappa_reg_, cutoff_to_spike_interval ); + weight = optimizer->optimized_weight( opt_cp, t_compute_until, sum_grad, weight ); } } diff --git a/models/eprop_iaf_psc_delta.h b/models/eprop_iaf_psc_delta.h index c36066a984..04738a7841 100644 --- a/models/eprop_iaf_psc_delta.h +++ b/models/eprop_iaf_psc_delta.h @@ -236,6 +236,8 @@ Parameter Unit Math equivalent Default Des ---------------------------------------------------------------------------------------------------------------- Parameter Unit Math equivalent Default Description =============================== ======= =========================== ================== ========================= +``activation_interval`` ms 3000.0 Interval between two + activations ``c_reg`` :math:`c_\text{reg}` 0.0 Coefficient of firing rate regularization ``eprop_isi_trace_cutoff`` ms :math:`{\Delta t}_\text{c}` maximum value Cutoff for integration of @@ -409,7 +411,10 @@ class eprop_iaf_psc_delta : public EpropArchivingNodeRecurrent< false > double&, double&, const CommonSynapseProperties&, - WeightOptimizer* ) override; + WeightOptimizer*, + const bool, + const bool, + double& ) override; long get_shift() const override; bool is_eprop_recurrent_node() const override; @@ -476,6 +481,9 @@ class eprop_iaf_psc_delta : public EpropArchivingNodeRecurrent< false > //! Time interval from the previous spike until the cutoff of e-prop update integration between two spikes (ms). double eprop_isi_trace_cutoff_; + //! Interval between two activations. + long activation_interval_; + //! Default constructor. Parameters_(); @@ -550,6 +558,9 @@ class eprop_iaf_psc_delta : public EpropArchivingNodeRecurrent< false > //! Time steps from the previous spike until the cutoff of e-prop update integration between two spikes. long eprop_isi_trace_cutoff_steps_; + + //! Time steps of activation interval. + long activation_interval_steps_; }; //! Get the current value of the membrane voltage. diff --git a/models/eprop_iaf_psc_delta_adapt.cpp b/models/eprop_iaf_psc_delta_adapt.cpp index 9ffe2a8b69..9ebf077b17 100644 --- a/models/eprop_iaf_psc_delta_adapt.cpp +++ b/models/eprop_iaf_psc_delta_adapt.cpp @@ -91,6 +91,7 @@ eprop_iaf_psc_delta_adapt::Parameters_::Parameters_() , kappa_( 0.97 ) , kappa_reg_( 0.97 ) , eprop_isi_trace_cutoff_( 1000.0 ) + , activation_interval_( 3000.0 ) { } @@ -143,6 +144,7 @@ eprop_iaf_psc_delta_adapt::Parameters_::get( DictionaryDatum& d ) const def< double >( d, names::kappa, kappa_ ); def< double >( d, names::kappa_reg, kappa_reg_ ); def< double >( d, names::eprop_isi_trace_cutoff, eprop_isi_trace_cutoff_ ); + def< double >( d, names::activation_interval, activation_interval_ ); } double @@ -184,6 +186,7 @@ eprop_iaf_psc_delta_adapt::Parameters_::set( const DictionaryDatum& d, Node* nod updateValueParam< double >( d, names::kappa, kappa_, node ); updateValueParam< double >( d, names::kappa_reg, kappa_reg_, node ); updateValueParam< double >( d, names::eprop_isi_trace_cutoff, eprop_isi_trace_cutoff_, node ); + updateValueParam< double >( d, names::activation_interval, activation_interval_, node ); if ( V_th_ < V_min_ ) { @@ -245,11 +248,16 @@ eprop_iaf_psc_delta_adapt::Parameters_::set( const DictionaryDatum& d, Node* nod throw BadProperty( "Firing rate low-pass filter for regularization kappa_reg from range [0, 1] required." ); } - if ( eprop_isi_trace_cutoff_ < 0.0 ) + if ( activation_interval_ < 0 ) { - throw BadProperty( "Cutoff of integration of eprop trace between spikes eprop_isi_trace_cutoff ≥ 0 required." ); + throw BadProperty( "Interval between activations activation_interval ≥ 0 required." ); } + if ( eprop_isi_trace_cutoff_ < 0.0 or eprop_isi_trace_cutoff_ > activation_interval_ ) + { + throw BadProperty( + "Computation cutoff of eprop trace 0 ≤ eprop trace eprop_isi_trace_cutoff ≤ activation_interval required." ); + } return delta_EL; } @@ -321,6 +329,7 @@ eprop_iaf_psc_delta_adapt::pre_run_hook() V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); V_.eprop_isi_trace_cutoff_steps_ = Time( Time::ms( P_.eprop_isi_trace_cutoff_ ) ).get_steps(); + V_.activation_interval_steps_ = Time( Time::ms( P_.activation_interval_ ) ).get_steps(); // calculate the entries of the propagator matrix for the evolution of the state vector @@ -386,6 +395,14 @@ eprop_iaf_psc_delta_adapt::update( Time const& origin, const long from, const lo kernel().event_delivery_manager.send( *this, se, lag ); S_.z_ = 1.0; + set_last_event_time( t ); + } + else if ( get_last_event_time() > 0 and t - get_last_event_time() >= V_.activation_interval_steps_ ) + { + SpikeEvent se; + se.set_activation(); + kernel().event_delivery_manager.send( *this, se, lag ); + set_last_event_time( t ); } append_new_eprop_history_entry( t ); @@ -452,63 +469,72 @@ eprop_iaf_psc_delta_adapt::compute_gradient( const long t_spike, double& epsilon, double& weight, const CommonSynapseProperties& cp, - WeightOptimizer* optimizer ) + WeightOptimizer* optimizer, + const bool activation, + const bool previous_event_was_activation, + double& sum_grad ) { - double e = 0.0; // eligibility trace - double z = 0.0; // spiking variable - double z_current_buffer = 1.0; // buffer containing the spike that triggered the current integration - double psi = 0.0; // surrogate gradient - double L = 0.0; // learning signal - double firing_rate_reg = 0.0; // firing rate regularization - double grad = 0.0; // gradient + const auto& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); + const auto& opt_cp = *ecp.optimizer_cp_; + const bool optimize_each_step = opt_cp.optimize_each_step_; - const EpropSynapseCommonProperties& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); - const auto optimize_each_step = ( *ecp.optimizer_cp_ ).optimize_each_step_; + if ( not previous_event_was_activation ) + { + sum_grad = 0.0; // sum of gradients + } auto eprop_hist_it = get_eprop_history( t_spike_previous - 1 ); - const long t_compute_until = std::min( t_spike_previous + V_.eprop_isi_trace_cutoff_steps_, t_spike ); + const long cutoff_end = t_spike_previous + V_.eprop_isi_trace_cutoff_steps_; + const long t_compute_until = std::min( cutoff_end, t_spike ); - for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) + if ( not previous_event_was_activation ) { - z = z_previous_buffer; - z_previous_buffer = z_current_buffer; - z_current_buffer = 0.0; + double z_current_buffer = 1.0; // spike that triggered current computation - psi = eprop_hist_it->surrogate_gradient_; - L = eprop_hist_it->learning_signal_; - firing_rate_reg = eprop_hist_it->firing_rate_reg_; + for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) + { + const double z = z_previous_buffer; // spiking variable + z_previous_buffer = z_current_buffer; + z_current_buffer = 0.0; - z_bar = V_.P_v_m_ * z_bar + z; - e = psi * ( z_bar - P_.adapt_beta_ * epsilon ); - epsilon = V_.P_adapt_ * epsilon + e; - e_bar = P_.kappa_ * e_bar + e; - e_bar_reg = P_.kappa_reg_ * e_bar_reg + ( 1.0 - P_.kappa_reg_ ) * e; + const double psi = eprop_hist_it->surrogate_gradient_; // surrogate gradient + const double L = eprop_hist_it->learning_signal_; // learning signal + const double firing_rate_reg = eprop_hist_it->firing_rate_reg_; // firing rate regularization - if ( optimize_each_step ) - { - grad = L * e_bar + firing_rate_reg * e_bar_reg; - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t, grad, weight ); - } - else - { - grad += L * e_bar + firing_rate_reg * e_bar_reg; + z_bar = V_.P_v_m_ * z_bar + z; + const double e = psi * ( z_bar - P_.adapt_beta_ * epsilon ); // eligibility trace + epsilon = V_.P_adapt_ * epsilon + e; + e_bar = P_.kappa_ * e_bar + e; + e_bar_reg = P_.kappa_reg_ * e_bar_reg + ( 1.0 - P_.kappa_reg_ ) * e; + + const double grad = L * e_bar + firing_rate_reg * e_bar_reg; + + if ( optimize_each_step ) + { + sum_grad = grad; + weight = optimizer->optimized_weight( opt_cp, t, sum_grad, weight ); + } + else + { + sum_grad += grad; + } } } - if ( not optimize_each_step ) + const long trace_decay_interval = t_spike - ( previous_event_was_activation ? t_spike_previous : t_compute_until ); + + if ( trace_decay_interval > 0 ) { - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t_compute_until, grad, weight ); + z_bar *= std::pow( V_.P_v_m_, trace_decay_interval ); + e_bar *= std::pow( P_.kappa_, trace_decay_interval ); + e_bar_reg *= std::pow( P_.kappa_reg_, trace_decay_interval ); + epsilon *= std::pow( V_.P_adapt_, trace_decay_interval ); } - const long cutoff_to_spike_interval = t_spike - t_compute_until; - - if ( cutoff_to_spike_interval > 0 ) + if ( not( activation or optimize_each_step ) ) { - z_bar *= std::pow( V_.P_v_m_, cutoff_to_spike_interval ); - e_bar *= std::pow( P_.kappa_, cutoff_to_spike_interval ); - e_bar_reg *= std::pow( P_.kappa_reg_, cutoff_to_spike_interval ); - epsilon *= std::pow( V_.P_adapt_, cutoff_to_spike_interval ); + weight = optimizer->optimized_weight( opt_cp, t_compute_until, sum_grad, weight ); } } diff --git a/models/eprop_iaf_psc_delta_adapt.h b/models/eprop_iaf_psc_delta_adapt.h index 3c0949de3a..77110ddb48 100644 --- a/models/eprop_iaf_psc_delta_adapt.h +++ b/models/eprop_iaf_psc_delta_adapt.h @@ -249,6 +249,8 @@ Parameter Unit Math equivalent Default Des ---------------------------------------------------------------------------------------------------------------- Parameter Unit Math equivalent Default Description =============================== ======= =========================== ================== ========================= +``activation_interval`` ms 3000.0 Interval between two + activations ``c_reg`` :math:`c_\text{reg}` 0.0 Coefficient of firing rate regularization ``eprop_isi_trace_cutoff`` ms :math:`{\Delta t}_\text{c}` maximum value Cutoff for integration of @@ -425,7 +427,10 @@ class eprop_iaf_psc_delta_adapt : public EpropArchivingNodeRecurrent< false > double&, double&, const CommonSynapseProperties&, - WeightOptimizer* ) override; + WeightOptimizer*, + const bool, + const bool, + double& ) override; long get_shift() const override; bool is_eprop_recurrent_node() const override; @@ -498,6 +503,9 @@ class eprop_iaf_psc_delta_adapt : public EpropArchivingNodeRecurrent< false > //! Time interval from the previous spike until the cutoff of e-prop update integration between two spikes (ms). double eprop_isi_trace_cutoff_; + //! Interval between two activations. + long activation_interval_; + //! Default constructor. Parameters_(); @@ -584,6 +592,9 @@ class eprop_iaf_psc_delta_adapt : public EpropArchivingNodeRecurrent< false > //! Time steps from the previous spike until the cutoff of e-prop update integration between two spikes. long eprop_isi_trace_cutoff_steps_; + + //! Time steps of activation interval. + long activation_interval_steps_; }; //! Get the current value of the membrane voltage. diff --git a/models/eprop_input_neuron.cpp b/models/eprop_input_neuron.cpp new file mode 100644 index 0000000000..b7b8aa5d2c --- /dev/null +++ b/models/eprop_input_neuron.cpp @@ -0,0 +1,113 @@ +/* + * eprop_input_neuron.cpp + * + * 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 . + * + */ + + +#include "eprop_input_neuron.h" + +// Includes from libnestutil: +#include "numerics.h" + +// Includes from nestkernel: +#include "event_delivery_manager_impl.h" +#include "exceptions.h" +#include "kernel_manager.h" +#include "nest_impl.h" + +// Includes from sli: +#include "dictutils.h" + +namespace nest +{ +void +register_eprop_input_neuron( const std::string& name ) +{ + register_node_model< eprop_input_neuron >( name ); +} + + +eprop_input_neuron::eprop_input_neuron() + : ArchivingNode() +{ +} + +void +eprop_input_neuron::init_buffers_() +{ + B_.n_spikes_.clear(); // includes resize + ArchivingNode::clear_history(); +} + +void +eprop_input_neuron::update( Time const& origin, const long from, const long to ) +{ + for ( long lag = from; lag < to; ++lag ) + { + const long t = origin.get_steps() + lag; + const unsigned long current_spikes_n = static_cast< unsigned long >( B_.n_spikes_.get_value( lag ) ); + if ( current_spikes_n > 0 ) + { + // create a new SpikeEvent, set its multiplicity and send it + SpikeEvent se; + se.set_multiplicity( current_spikes_n ); + kernel().event_delivery_manager.send( *this, se, lag ); + + // set the spike times, respecting the multiplicity + for ( unsigned long i = 0; i < current_spikes_n; i++ ) + { + set_spiketime( Time::step( t + 1 ) ); + } + last_event_time_ = t; + } + else if ( last_event_time_ > 0 and t - last_event_time_ >= activation_interval_ ) + { + SpikeEvent se; + se.set_activation(); + kernel().event_delivery_manager.send( *this, se, lag ); + last_event_time_ = t; + } + } +} + +void +eprop_input_neuron::get_status( DictionaryDatum& d ) const +{ + ArchivingNode::get_status( d ); +} + +void +eprop_input_neuron::set_status( const DictionaryDatum& d ) +{ + ArchivingNode::set_status( d ); +} + +void +eprop_input_neuron::handle( SpikeEvent& e ) +{ + // Repeat only spikes incoming on port 0, port 1 will be ignored + if ( 0 == e.get_rport() ) + { + B_.n_spikes_.add_value( e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ), + static_cast< double >( e.get_multiplicity() ) ); + } +} + +} // namespace diff --git a/models/eprop_input_neuron.h b/models/eprop_input_neuron.h new file mode 100644 index 0000000000..58971aed11 --- /dev/null +++ b/models/eprop_input_neuron.h @@ -0,0 +1,170 @@ +/* + * eprop_input_neuron.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 . + * + */ + +#ifndef EPROP_INPUT_NEURON_H +#define EPROP_INPUT_NEURON_H + +// Includes from nestkernel: +#include "archiving_node.h" +#include "connection.h" +#include "event.h" +#include "nest_types.h" +#include "ring_buffer.h" + +namespace nest +{ + +/* BeginUserDocs: neuron, parrot + +Short description ++++++++++++++++++ + +Neuron that repeats incoming spikes + +Description ++++++++++++ + +The parrot neuron simply emits one spike for every incoming spike. +An important application is to provide identical poisson spike +trains to a group of neurons. The ``poisson_generator`` sends a different +spike train to each of its target neurons. By connecting one +``poisson_generator`` to a ``eprop_input_neuron`` and then that ``eprop_input_neuron`` to +a group of neurons, all target neurons will receive the same poisson +spike train. + +Please note that weights of connections *to* the ``eprop_input_neuron`` +are ignored, while weights on connections *from* the ``eprop_input_neuron`` +to the target are handled as usual. Delays are honored on both +incoming and outgoing connections. + +Only spikes arriving on connections to port 0 will be repeated. +Connections onto port 1 will be accepted, but spikes incoming +through port 1 will be ignored. This allows setting exact pre- +and postsynaptic spike times for STDP protocols by connecting +two parrot neurons spiking at desired times by, for example, a +``stdp_synapse`` onto port 1 on the postsynaptic parrot neuron. + +Receives +++++++++ + +SpikeEvent + +Sends ++++++ + +SpikeEvent + +Examples using this model ++++++++++++++++++++++++++ + +.. listexamples:: eprop_input_neuron + +EndUserDocs */ + +void register_eprop_input_neuron( const std::string& name ); + +class eprop_input_neuron : public ArchivingNode +{ + +public: + eprop_input_neuron(); + + /** + * Import sets of overloaded virtual functions. + * @see Technical Issues / Virtual Functions: Overriding, + * Overloading, and Hiding + */ + using Node::handle; + using Node::handles_test_event; + using Node::receives_signal; + using Node::sends_signal; + + size_t send_test_event( Node&, size_t, synindex, bool ) override; + SignalType sends_signal() const override; + SignalType receives_signal() const override; + + void handle( SpikeEvent& ) override; + size_t handles_test_event( SpikeEvent&, size_t ) override; + + void get_status( DictionaryDatum& ) const override; + void set_status( const DictionaryDatum& ) override; + +private: + void init_buffers_() override; + void + pre_run_hook() override + { + } // no variables + + void update( Time const&, const long, const long ) override; + + /** + Buffers and accumulates the number of incoming spikes per time step; + RingBuffer stores doubles; for now the numbers are casted. + */ + struct Buffers_ + { + RingBuffer n_spikes_; + }; + + Buffers_ B_; +}; + +inline size_t +eprop_input_neuron::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_input_neuron::handles_test_event( SpikeEvent&, size_t receptor_type ) +{ + // Allow connections to port 0 (spikes to be repeated) + // and port 1 (spikes to be ignored). + if ( receptor_type == 0 or receptor_type == 1 ) + { + return receptor_type; + } + else + { + throw UnknownReceptorType( receptor_type, get_name() ); + } +} + +inline SignalType +eprop_input_neuron::sends_signal() const +{ + return ALL; +} + +inline SignalType +eprop_input_neuron::receives_signal() const +{ + return ALL; +} + +} // namespace + +#endif // EPROP_INPUT_NEURON_H diff --git a/models/eprop_readout.cpp b/models/eprop_readout.cpp index b2f740d70e..4a9507bd1c 100644 --- a/models/eprop_readout.cpp +++ b/models/eprop_readout.cpp @@ -314,51 +314,63 @@ eprop_readout::compute_gradient( const long t_spike, double& epsilon, double& weight, const CommonSynapseProperties& cp, - WeightOptimizer* optimizer ) + WeightOptimizer* optimizer, + const bool activation, + const bool previous_event_was_activation, + double& sum_grad ) { - double z = 0.0; // spiking variable - double z_current_buffer = 1.0; // buffer containing the spike that triggered the current integration - double L = 0.0; // error signal - double grad = 0.0; // gradient + const auto& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); + const auto& opt_cp = *ecp.optimizer_cp_; + const bool optimize_each_step = opt_cp.optimize_each_step_; - const EpropSynapseCommonProperties& ecp = static_cast< const EpropSynapseCommonProperties& >( cp ); - const auto optimize_each_step = ( *ecp.optimizer_cp_ ).optimize_each_step_; + if ( not previous_event_was_activation ) + { + sum_grad = 0.0; // sum of gradients + } auto eprop_hist_it = get_eprop_history( t_spike_previous - 1 ); - const long t_compute_until = std::min( t_spike_previous + V_.eprop_isi_trace_cutoff_steps_, t_spike ); + const long cutoff_end = t_spike_previous + V_.eprop_isi_trace_cutoff_steps_; + const long t_compute_until = std::min( cutoff_end, t_spike ); - for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) + if ( not previous_event_was_activation ) { - z = z_previous_buffer; - z_previous_buffer = z_current_buffer; - z_current_buffer = 0.0; + double z_current_buffer = 1.0; // spike that triggered current computation - L = eprop_hist_it->error_signal_; + for ( long t = t_spike_previous; t < t_compute_until; ++t, ++eprop_hist_it ) + { + const double z = z_previous_buffer; // spiking variable + z_previous_buffer = z_current_buffer; + z_current_buffer = 0.0; - z_bar = V_.P_v_m_ * z_bar + z; + const double E = eprop_hist_it->error_signal_; // error signal - if ( optimize_each_step ) - { - grad = L * z_bar; - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t, grad, weight ); - } - else - { - grad += L * z_bar; + z_bar = V_.P_v_m_ * z_bar + z; + + const double grad = E * z_bar; + + if ( optimize_each_step ) + { + sum_grad = grad; + weight = optimizer->optimized_weight( opt_cp, t, sum_grad, weight ); + } + else + { + sum_grad += grad; + } } } - if ( not optimize_each_step ) + const long trace_decay_interval = t_spike - ( previous_event_was_activation ? t_spike_previous : t_compute_until ); + + if ( trace_decay_interval > 0 ) { - weight = optimizer->optimized_weight( *ecp.optimizer_cp_, t_compute_until, grad, weight ); + z_bar *= std::pow( V_.P_v_m_, trace_decay_interval ); } - const long cutoff_to_spike_interval = t_spike - t_compute_until; - - if ( cutoff_to_spike_interval > 0 ) + if ( not( activation or optimize_each_step ) ) { - z_bar *= std::pow( V_.P_v_m_, cutoff_to_spike_interval ); + weight = optimizer->optimized_weight( opt_cp, t_compute_until, sum_grad, weight ); } } diff --git a/models/eprop_readout.h b/models/eprop_readout.h index b96c007542..e3e0492156 100644 --- a/models/eprop_readout.h +++ b/models/eprop_readout.h @@ -308,7 +308,10 @@ class eprop_readout : public EpropArchivingNodeReadout< false > double&, double&, const CommonSynapseProperties&, - WeightOptimizer* ) override; + WeightOptimizer*, + const bool, + const bool, + double& ) override; long get_shift() const override; bool is_eprop_recurrent_node() const override; diff --git a/models/eprop_readout_bsshslm_2020.cpp b/models/eprop_readout_bsshslm_2020.cpp index 79f059fb59..10d70f2d59 100644 --- a/models/eprop_readout_bsshslm_2020.cpp +++ b/models/eprop_readout_bsshslm_2020.cpp @@ -387,25 +387,22 @@ eprop_readout_bsshslm_2020::compute_gradient( std::vector< long >& presyn_isis, auto eprop_hist_it = get_eprop_history( t_previous_trigger_spike ); double grad = 0.0; // gradient value to be calculated - double L = 0.0; // error signal - double z = 0.0; // spiking variable double z_bar = 0.0; // low-pass filtered spiking variable - for ( long presyn_isi : presyn_isis ) + for ( const long presyn_isi : presyn_isis ) { - z = 1.0; // set spiking variable to 1 for each incoming spike + double z = 1.0; // set spiking variable to 1 for each incoming spike - for ( long t = 0; t < presyn_isi; ++t ) + for ( long t = 0; t < presyn_isi; ++t, ++eprop_hist_it ) { assert( eprop_hist_it != eprop_history_.end() ); - L = eprop_hist_it->error_signal_; + const double E = eprop_hist_it->error_signal_; // error signal z_bar = V_.P_v_m_ * z_bar + V_.P_z_in_ * z; - grad += L * z_bar; - z = 0.0; // set spiking variable to 0 between spikes + grad += E * z_bar; - ++eprop_hist_it; + z = 0.0; // set spiking variable to 0 between spikes } } presyn_isis.clear(); diff --git a/models/eprop_synapse.h b/models/eprop_synapse.h index 5509019296..b109ecc02b 100644 --- a/models/eprop_synapse.h +++ b/models/eprop_synapse.h @@ -317,6 +317,9 @@ class eprop_synapse : public Connection< targetidentifierT > //! The time step when the previous spike arrived. long t_spike_previous_ = 0; + //! Last event was an activation. + bool previous_event_was_activation_ = false; + //! The time step when the spike arrived that triggered the previous e-prop update. long t_previous_trigger_spike_ = 0; @@ -335,6 +338,9 @@ class eprop_synapse : public Connection< targetidentifierT > //! Value of spiking variable one time step before t_previous_spike_. double z_previous_buffer_ = 0.0; + //! Sum of gradients. + double sum_grad_ = 0.0; + /** * Optimizer * @@ -365,6 +371,7 @@ eprop_synapse< targetidentifierT >::eprop_synapse() : ConnectionBase() , weight_( 1.0 ) , t_spike_previous_( 0 ) + , previous_event_was_activation_( false ) , t_previous_trigger_spike_( 0 ) , optimizer_( nullptr ) { @@ -399,12 +406,14 @@ eprop_synapse< targetidentifierT >::operator=( const eprop_synapse& es ) weight_ = es.weight_; t_spike_previous_ = es.t_spike_previous_; + previous_event_was_activation_ = es.previous_event_was_activation_; t_previous_trigger_spike_ = es.t_previous_trigger_spike_; z_bar_ = es.z_bar_; e_bar_ = es.e_bar_; e_bar_reg_ = es.e_bar_reg_; epsilon_ = es.epsilon_; z_previous_buffer_ = es.z_previous_buffer_; + sum_grad_ = es.sum_grad_; optimizer_ = es.optimizer_; return *this; @@ -415,11 +424,14 @@ eprop_synapse< targetidentifierT >::eprop_synapse( eprop_synapse&& es ) : ConnectionBase( es ) , weight_( es.weight_ ) , t_spike_previous_( es.t_spike_previous_ ) + , previous_event_was_activation_( es.previous_event_was_activation_ ) , t_previous_trigger_spike_( es.t_previous_trigger_spike_ ) , z_bar_( es.z_bar_ ) , e_bar_( es.e_bar_ ) , e_bar_reg_( es.e_bar_reg_ ) , epsilon_( es.epsilon_ ) + , z_previous_buffer_( es.z_previous_buffer_ ) + , sum_grad_( es.sum_grad_ ) , optimizer_( es.optimizer_ ) { // Move operator, therefore we must null the optimizer pointer in the source of the move. @@ -440,13 +452,14 @@ eprop_synapse< targetidentifierT >::operator=( eprop_synapse&& es ) weight_ = es.weight_; t_spike_previous_ = es.t_spike_previous_; + previous_event_was_activation_ = es.previous_event_was_activation_; t_previous_trigger_spike_ = es.t_previous_trigger_spike_; z_bar_ = es.z_bar_; e_bar_ = es.e_bar_; e_bar_reg_ = es.e_bar_reg_; epsilon_ = es.epsilon_; z_previous_buffer_ = es.z_previous_buffer_; - + sum_grad_ = es.sum_grad_; optimizer_ = es.optimizer_; // Move assignment, therefore we must null the optimizer pointer in the source of the move. @@ -471,7 +484,7 @@ eprop_synapse< targetidentifierT >::check_connection( Node& s, ConnTestDummyNode dummy_target; ConnectionBase::check_connection_( dummy_target, s, t, receptor_type ); - t.register_eprop_connection(); + t.register_synapse(); optimizer_ = cp.optimizer_cp_->get_optimizer(); } @@ -492,24 +505,43 @@ eprop_synapse< targetidentifierT >::send( Event& e, size_t thread, const EpropSy assert( target ); const long t_spike = e.get_stamp().get_steps(); + const bool activation = e.get_activation(); + + if ( previous_event_was_activation_ and activation ) + { + return false; + } if ( t_spike_previous_ != 0 ) { - target->compute_gradient( - t_spike, t_spike_previous_, z_previous_buffer_, z_bar_, e_bar_, e_bar_reg_, epsilon_, weight_, cp, optimizer_ ); + target->compute_gradient( t_spike, + t_spike_previous_, + z_previous_buffer_, + z_bar_, + e_bar_, + e_bar_reg_, + epsilon_, + weight_, + cp, + optimizer_, + activation, + previous_event_was_activation_, + sum_grad_ ); } - const long eprop_isi_trace_cutoff = target->get_eprop_isi_trace_cutoff(); - target->write_update_to_history( t_spike_previous_, t_spike, eprop_isi_trace_cutoff ); + target->erase_used_eprop_history( t_spike, t_spike_previous_ ); t_spike_previous_ = t_spike; + previous_event_was_activation_ = activation; - e.set_receiver( *target ); - e.set_weight( weight_ ); - e.set_delay_steps( get_delay_steps() ); - e.set_rport( get_rport() ); - e(); - + if ( not activation ) + { + e.set_receiver( *target ); + e.set_weight( weight_ ); + e.set_delay_steps( get_delay_steps() ); + e.set_rport( get_rport() ); + e(); + } return true; } diff --git a/models/eprop_synapse_bsshslm_2020.h b/models/eprop_synapse_bsshslm_2020.h index 3ac9144827..e02fc8055f 100644 --- a/models/eprop_synapse_bsshslm_2020.h +++ b/models/eprop_synapse_bsshslm_2020.h @@ -331,9 +331,14 @@ class eprop_synapse_bsshslm_2020 : public Connection< targetidentifierT > //! Synaptic weight. double weight_; + double gradient_sum_; + //! The time step when the previous spike arrived. long t_spike_previous_; + //! Last event was an activation. + bool previous_event_was_activation_; + //! The time step when the previous e-prop update was. long t_previous_update_; @@ -384,7 +389,9 @@ template < typename targetidentifierT > eprop_synapse_bsshslm_2020< targetidentifierT >::eprop_synapse_bsshslm_2020() : ConnectionBase() , weight_( 1.0 ) + , gradient_sum_( 0.0 ) , t_spike_previous_( 0 ) + , previous_event_was_activation_( false ) , t_previous_update_( 0 ) , t_next_update_( 0 ) , t_previous_trigger_spike_( 0 ) @@ -406,7 +413,9 @@ template < typename targetidentifierT > eprop_synapse_bsshslm_2020< targetidentifierT >::eprop_synapse_bsshslm_2020( const eprop_synapse_bsshslm_2020& es ) : ConnectionBase( es ) , weight_( es.weight_ ) + , gradient_sum_( es.gradient_sum_ ) , t_spike_previous_( 0 ) + , previous_event_was_activation_( false ) , t_previous_update_( 0 ) , t_next_update_( kernel().simulation_manager.get_eprop_update_interval().get_steps() ) , t_previous_trigger_spike_( 0 ) @@ -430,7 +439,9 @@ eprop_synapse_bsshslm_2020< targetidentifierT >::operator=( const eprop_synapse_ ConnectionBase::operator=( es ); weight_ = es.weight_; + gradient_sum_ = es.gradient_sum_; t_spike_previous_ = es.t_spike_previous_; + previous_event_was_activation_ = es.previous_event_was_activation_; t_previous_update_ = es.t_previous_update_; t_next_update_ = es.t_next_update_; t_previous_trigger_spike_ = es.t_previous_trigger_spike_; @@ -446,7 +457,9 @@ template < typename targetidentifierT > eprop_synapse_bsshslm_2020< targetidentifierT >::eprop_synapse_bsshslm_2020( eprop_synapse_bsshslm_2020&& es ) : ConnectionBase( es ) , weight_( es.weight_ ) + , gradient_sum_( es.gradient_sum_ ) , t_spike_previous_( 0 ) + , previous_event_was_activation_( false ) , t_previous_update_( 0 ) , t_next_update_( es.t_next_update_ ) , t_previous_trigger_spike_( 0 ) @@ -471,7 +484,9 @@ eprop_synapse_bsshslm_2020< targetidentifierT >::operator=( eprop_synapse_bsshsl ConnectionBase::operator=( es ); weight_ = es.weight_; + gradient_sum_ = es.gradient_sum_; t_spike_previous_ = es.t_spike_previous_; + previous_event_was_activation_ = es.previous_event_was_activation_; t_previous_update_ = es.t_previous_update_; t_next_update_ = es.t_next_update_; t_previous_trigger_spike_ = es.t_previous_trigger_spike_; @@ -524,6 +539,19 @@ eprop_synapse_bsshslm_2020< targetidentifierT >::send( Event& e, assert( target ); const long t_spike = e.get_stamp().get_steps(); + const bool activation = e.get_activation(); + + if ( previous_event_was_activation_ ) + { + if ( activation ) + { + return false; + } + + t_spike_previous_ = t_spike; + t_previous_trigger_spike_ = t_spike; + } + const long update_interval = kernel().simulation_manager.get_eprop_update_interval().get_steps(); const long shift = target->get_shift(); @@ -541,8 +569,8 @@ eprop_synapse_bsshslm_2020< targetidentifierT >::send( Event& e, if ( t_spike_previous_ > 0 ) { - const long t = t_spike >= t_next_update_ + shift ? t_next_update_ + shift : t_spike; - presyn_isis_.push_back( t - t_spike_previous_ ); + presyn_isis_.push_back( + std::min( t_spike, t_next_update_ + shift ) - std::min( t_spike_previous_, t_next_update_ + shift ) ); } if ( t_spike > t_next_update_ + shift ) @@ -550,26 +578,50 @@ eprop_synapse_bsshslm_2020< targetidentifierT >::send( Event& e, const long idx_current_update = ( t_spike - shift ) / update_interval; const long t_current_update = idx_current_update * update_interval; - target->write_update_to_history( t_previous_update_, t_current_update ); + target->write_update_to_history( t_previous_update_, t_current_update, activation, previous_event_was_activation_ ); const double gradient = target->compute_gradient( presyn_isis_, t_previous_update_, t_previous_trigger_spike_, kappa_, cp.average_gradient_ ); - weight_ = optimizer_->optimized_weight( *cp.optimizer_cp_, idx_current_update, gradient, weight_ ); + gradient_sum_ += gradient; + if ( not activation ) + { + weight_ = optimizer_->optimized_weight( *cp.optimizer_cp_, idx_current_update, gradient_sum_, weight_ ); + gradient_sum_ = 0.0; + } t_previous_update_ = t_current_update; t_next_update_ = t_current_update + update_interval; t_previous_trigger_spike_ = t_spike; } + else + { + if ( not activation and previous_event_was_activation_ ) + { + const long idx_current_update = ( t_spike - shift ) / update_interval; + const long t_current_update = idx_current_update * update_interval; + + target->write_update_to_history( + t_previous_update_, t_current_update, activation, previous_event_was_activation_ ); + + weight_ = optimizer_->optimized_weight( *cp.optimizer_cp_, idx_current_update, gradient_sum_, weight_ ); + gradient_sum_ = 0.0; + } + } - t_spike_previous_ = t_spike; + if ( not activation ) + { + t_spike_previous_ = t_spike; + + e.set_receiver( *target ); + e.set_weight( weight_ ); + e.set_delay_steps( get_delay_steps() ); + e.set_rport( get_rport() ); + e(); + } - e.set_receiver( *target ); - e.set_weight( weight_ ); - e.set_delay_steps( get_delay_steps() ); - e.set_rport( get_rport() ); - e(); + previous_event_was_activation_ = activation; return true; } diff --git a/modelsets/eprop b/modelsets/eprop index 0a9a5ed696..6825e1d7d4 100644 --- a/modelsets/eprop +++ b/modelsets/eprop @@ -4,7 +4,6 @@ multimeter spike_recorder weight_recorder -parrot_neuron poisson_generator spike_generator step_rate_generator @@ -12,6 +11,7 @@ step_rate_generator rate_connection_delayed static_synapse +eprop_input_neuron eprop_iaf_bsshslm_2020 eprop_iaf_adapt_bsshslm_2020 eprop_readout_bsshslm_2020 diff --git a/modelsets/full b/modelsets/full index 46d6c6ee5a..67c021fb35 100644 --- a/modelsets/full +++ b/modelsets/full @@ -21,6 +21,7 @@ correlomatrix_detector correlospinmatrix_detector dc_generator diffusion_connection +eprop_input_neuron eprop_iaf_bsshslm_2020 eprop_iaf_adapt_bsshslm_2020 eprop_readout_bsshslm_2020 diff --git a/nestkernel/archiving_node.cpp b/nestkernel/archiving_node.cpp index 7ad511eb3c..f1d8e7c574 100644 --- a/nestkernel/archiving_node.cpp +++ b/nestkernel/archiving_node.cpp @@ -34,7 +34,9 @@ namespace nest // member functions for ArchivingNode nest::ArchivingNode::ArchivingNode() - : n_incoming_( 0 ) + : activation_interval_( 3 * kernel().simulation_manager.get_eprop_update_interval().get_steps() ) + , last_event_time_( 0 ) + , n_incoming_( 0 ) , Kminus_( 0.0 ) , Kminus_triplet_( 0.0 ) , tau_minus_( 20.0 ) @@ -49,6 +51,8 @@ nest::ArchivingNode::ArchivingNode() nest::ArchivingNode::ArchivingNode( const ArchivingNode& n ) : StructuralPlasticityNode( n ) + , activation_interval_( n.activation_interval_ ) + , last_event_time_( n.last_event_time_ ) , n_incoming_( n.n_incoming_ ) , Kminus_( n.Kminus_ ) , Kminus_triplet_( n.Kminus_triplet_ ) diff --git a/nestkernel/archiving_node.h b/nestkernel/archiving_node.h index 86b402ed1a..6eef397cf8 100644 --- a/nestkernel/archiving_node.h +++ b/nestkernel/archiving_node.h @@ -98,6 +98,12 @@ class ArchivingNode : public StructuralPlasticityNode void set_status( const DictionaryDatum& d ) override; protected: + //! Interval between two activations in steps. + long activation_interval_; + + //! Time of last spike or activation event in steps. + long last_event_time_; + /** * Record spike history */ diff --git a/nestkernel/eprop_archiving_node.h b/nestkernel/eprop_archiving_node.h index 04cfc2d3ba..0bbe4f4d19 100644 --- a/nestkernel/eprop_archiving_node.h +++ b/nestkernel/eprop_archiving_node.h @@ -66,10 +66,13 @@ class EpropArchivingNode : public Node */ EpropArchivingNode( const EpropArchivingNode& other ); + void register_synapse() override; void register_eprop_connection() override; void write_update_to_history( const long t_previous_update, const long t_current_update, + const bool activation, + const bool previous_event_was_activation, const long eprop_isi_trace_cutoff = 0 ) override; /** @@ -104,7 +107,7 @@ class EpropArchivingNode : public Node * * @param eprop_isi_trace_cutoff The cutoff value for the inter-spike integration of the eprop trace. */ - void erase_used_eprop_history( const long eprop_isi_trace_cutoff ); + void erase_used_eprop_history( const long t_spike, const long t_spike_previous ); /** * @brief Retrieves eprop history size. diff --git a/nestkernel/eprop_archiving_node_impl.h b/nestkernel/eprop_archiving_node_impl.h index 7c4c60a52e..f8de1345ab 100644 --- a/nestkernel/eprop_archiving_node_impl.h +++ b/nestkernel/eprop_archiving_node_impl.h @@ -48,6 +48,13 @@ EpropArchivingNode< HistEntryT >::EpropArchivingNode( const EpropArchivingNode& { } +template < typename HistEntryT > +void +EpropArchivingNode< HistEntryT >::register_synapse() +{ + ++eprop_indegree_; +} + template < typename HistEntryT > void EpropArchivingNode< HistEntryT >::register_eprop_connection() @@ -72,6 +79,8 @@ template < typename HistEntryT > void EpropArchivingNode< HistEntryT >::write_update_to_history( const long t_previous_update, const long t_current_update, + const bool activation, + const bool previous_event_was_activation, const long eprop_isi_trace_cutoff ) { if ( eprop_indegree_ == 0 ) @@ -80,32 +89,35 @@ EpropArchivingNode< HistEntryT >::write_update_to_history( const long t_previous } const long shift = model_dependent_history_shift_(); + const long t_curr_update_shifted = t_current_update + shift; + const long t_prev_update_shifted = t_previous_update + shift; - const auto it_hist_curr = get_update_history( t_current_update + shift ); - - if ( it_hist_curr != update_history_.end() and it_hist_curr->t_ == t_current_update + shift ) + if ( not activation ) { - ++it_hist_curr->access_counter_; - } - else - { - update_history_.insert( it_hist_curr, HistEntryEpropUpdate( t_current_update + shift, 1 ) ); + auto it_hist_curr = get_update_history( t_curr_update_shifted ); - if ( not history_shift_required_() ) + if ( it_hist_curr != update_history_.end() and it_hist_curr->t_ == t_curr_update_shifted ) { - erase_used_eprop_history( eprop_isi_trace_cutoff ); + ++it_hist_curr->access_counter_; + } + else + { + update_history_.insert( it_hist_curr, HistEntryEpropUpdate( t_curr_update_shifted, 1 ) ); } } - const auto it_hist_prev = get_update_history( t_previous_update + shift ); - - if ( it_hist_prev != update_history_.end() and it_hist_prev->t_ == t_previous_update + shift ) + if ( not previous_event_was_activation ) { - // If an entry exists for the previous update time, decrement its access counter - --it_hist_prev->access_counter_; - if ( it_hist_prev->access_counter_ == 0 ) + auto it_hist_prev = get_update_history( t_prev_update_shifted ); + + if ( it_hist_prev != update_history_.end() and it_hist_prev->t_ == t_prev_update_shifted ) { - update_history_.erase( it_hist_prev ); + // If an entry exists for the previous update time, decrement its access counter + --it_hist_prev->access_counter_; + if ( it_hist_prev->access_counter_ == 0 ) + { + update_history_.erase( it_hist_prev ); + } } } } @@ -159,27 +171,43 @@ EpropArchivingNode< HistEntryT >::erase_used_eprop_history() template < typename HistEntryT > void -EpropArchivingNode< HistEntryT >::erase_used_eprop_history( const long eprop_isi_trace_cutoff ) +EpropArchivingNode< HistEntryT >::erase_used_eprop_history( const long t_spike, const long t_spike_previous ) { - if ( eprop_history_.empty() // nothing to remove - or update_history_.size() < 2 // no time markers to check - ) + + if ( not update_history_.empty() and update_history_.back().t_ == t_spike ) { - return; + ++update_history_.back().access_counter_; + } + else + { + update_history_.emplace_back( t_spike, 1 ); } - const long t_prev = ( update_history_.end() - 2 )->t_; - const long t_curr = ( update_history_.end() - 1 )->t_; + if ( t_spike_previous != 0 ) + { + auto it_hist_prev = std::lower_bound( update_history_.begin(), update_history_.end(), t_spike_previous ); - if ( t_prev + eprop_isi_trace_cutoff < t_curr ) + if ( it_hist_prev != update_history_.end() and it_hist_prev->t_ == t_spike_previous ) + { + --it_hist_prev->access_counter_; + if ( it_hist_prev->access_counter_ == 0 ) + { + update_history_.erase( it_hist_prev ); + } + } + } + + if ( eprop_history_.empty() ) { - // erase no longer needed entries to be ignored by trace cutoff - eprop_history_.erase( get_eprop_history( t_prev + eprop_isi_trace_cutoff ), get_eprop_history( t_curr ) ); + return; } + const long time_end = update_history_.front().t_ - 1; + auto it_end = std::lower_bound( eprop_history_.begin(), eprop_history_.end(), time_end ); - // erase no longer needed entries before the earliest current update - eprop_history_.erase( - get_eprop_history( std::numeric_limits< long >::min() ), get_eprop_history( update_history_.begin()->t_ - 1 ) ); + if ( it_end != eprop_history_.end() and it_end->t_ == time_end ) + { + eprop_history_.erase( eprop_history_.begin(), it_end ); + } } template < typename HistEntryT > diff --git a/nestkernel/eprop_archiving_node_recurrent.h b/nestkernel/eprop_archiving_node_recurrent.h index a8872e3297..150c8c82d8 100644 --- a/nestkernel/eprop_archiving_node_recurrent.h +++ b/nestkernel/eprop_archiving_node_recurrent.h @@ -239,12 +239,25 @@ class EpropArchivingNodeRecurrent : public EpropArchivingNode< HistEntryEpropRec */ void reset_spike_count(); + /** + * Sets the time the neuron spiked. + */ + void set_last_event_time( const long last_event_time ); + + /** + * Gets the last time the neuron spiked. + */ + long get_last_event_time(); + //! Firing rate regularization. double firing_rate_reg_; //! Average firing rate. double f_av_; + //! Last time the neuron spiked. + long last_event_time_; + protected: long model_dependent_history_shift_() const override; bool history_shift_required_() const override; @@ -283,6 +296,20 @@ EpropArchivingNodeRecurrent< hist_shift_required >::reset_spike_count() n_spikes_ = 0; } +template < bool hist_shift_required > +inline void +EpropArchivingNodeRecurrent< hist_shift_required >::set_last_event_time( const long last_event_time ) +{ + last_event_time_ = last_event_time; +} + +template < bool hist_shift_required > +inline long +EpropArchivingNodeRecurrent< hist_shift_required >::get_last_event_time() +{ + return last_event_time_; +} + template < bool hist_shift_required > long EpropArchivingNodeRecurrent< hist_shift_required >::model_dependent_history_shift_() const diff --git a/nestkernel/eprop_archiving_node_recurrent_impl.h b/nestkernel/eprop_archiving_node_recurrent_impl.h index 58dde7a87e..be1df3e5d6 100644 --- a/nestkernel/eprop_archiving_node_recurrent_impl.h +++ b/nestkernel/eprop_archiving_node_recurrent_impl.h @@ -48,6 +48,7 @@ EpropArchivingNodeRecurrent< hist_shift_required >::EpropArchivingNodeRecurrent( : EpropArchivingNode() , firing_rate_reg_( 0.0 ) , f_av_( 0.0 ) + , last_event_time_( 0 ) , n_spikes_( 0 ) { } @@ -57,6 +58,7 @@ EpropArchivingNodeRecurrent< hist_shift_required >::EpropArchivingNodeRecurrent( : EpropArchivingNode( n ) , firing_rate_reg_( n.firing_rate_reg_ ) , f_av_( n.f_av_ ) + , last_event_time_( n.last_event_time_ ) , n_spikes_( n.n_spikes_ ) { } diff --git a/nestkernel/event.h b/nestkernel/event.h index 780234b4d5..d173b7c2fc 100644 --- a/nestkernel/event.h +++ b/nestkernel/event.h @@ -322,6 +322,21 @@ class Event */ void set_stamp( Time const& ); + /** + * Sets the activation flag. + */ + void set_activation(); + + /** + * Unsets the activation flag. + */ + void unset_activation(); + + /** + * Returns whether the activation flag is set. + */ + bool get_activation(); + protected: size_t sender_node_id_; //!< node ID of sender or 0 SpikeData sender_spike_data_; //!< spike data of sender node, in some cases required to retrieve node ID @@ -976,6 +991,24 @@ Event::set_stamp( Time const& s ) // get_rel_delivery_steps) } +inline void +Event::set_activation() +{ + sender_spike_data_.set_activation_marker(); +} + +inline void +Event::unset_activation() +{ + sender_spike_data_.unset_activation_marker(); +} + +inline bool +Event::get_activation() +{ + return sender_spike_data_.is_activation_marker(); +} + inline long Event::get_delay_steps() const { diff --git a/nestkernel/event_delivery_manager.cpp b/nestkernel/event_delivery_manager.cpp index 8ba790b740..ff656b3d03 100644 --- a/nestkernel/event_delivery_manager.cpp +++ b/nestkernel/event_delivery_manager.cpp @@ -664,6 +664,14 @@ EventDeliveryManager::deliver_events_( const size_t tid, const std::vector< Spik syn_id_batch[ j ] = spike_data.get_syn_id(); lcid_batch[ j ] = spike_data.get_lcid(); se_batch[ j ].set_sender_node_id_info( tid_batch[ j ], syn_id_batch[ j ], lcid_batch[ j ] ); + if ( spike_data.is_activation_marker() ) + { + se_batch[ j ].set_activation(); + } + else + { + se_batch[ j ].unset_activation(); + } } for ( size_t j = 0; j < SPIKES_PER_BATCH; ++j ) { @@ -685,6 +693,14 @@ EventDeliveryManager::deliver_events_( const size_t tid, const std::vector< Spik syn_id_batch[ j ] = spike_data.get_syn_id(); lcid_batch[ j ] = spike_data.get_lcid(); se_batch[ j ].set_sender_node_id_info( tid_batch[ j ], syn_id_batch[ j ], lcid_batch[ j ] ); + if ( spike_data.is_activation_marker() ) + { + se_batch[ j ].set_activation(); + } + else + { + se_batch[ j ].unset_activation(); + } } for ( size_t j = 0; j < num_remaining_entries; ++j ) { @@ -705,6 +721,14 @@ EventDeliveryManager::deliver_events_( const size_t tid, const std::vector< Spik se_batch[ j ].set_stamp( prepared_timestamps[ spike_data.get_lag() ] ); se_batch[ j ].set_offset( spike_data.get_offset() ); + if ( spike_data.is_activation_marker() ) + { + se_batch[ j ].set_activation(); + } + else + { + se_batch[ j ].unset_activation(); + } syn_id_batch[ j ] = spike_data.get_syn_id(); // for compressed spikes lcid holds the index in the // compressed_spike_data structure @@ -723,7 +747,16 @@ EventDeliveryManager::deliver_events_( const size_t tid, const std::vector< Spik { // non-local sender -> receiver retrieves ID of sender Node from SourceTable based on tid, syn_id, lcid // only if needed, as this is computationally costly + auto activation = se_batch[ j ].get_activation(); se_batch[ j ].set_sender_node_id_info( tid, syn_id_batch[ j ], lcid_batch[ j ] ); + if ( activation ) + { + se_batch[ j ].set_activation(); + } + else + { + se_batch[ j ].unset_activation(); + } } } for ( size_t j = 0; j < SPIKES_PER_BATCH; ++j ) @@ -742,6 +775,14 @@ EventDeliveryManager::deliver_events_( const size_t tid, const std::vector< Spik recv_buffer[ rank * spike_buffer_size_per_rank + num_batches * SPIKES_PER_BATCH + j ]; se_batch[ j ].set_stamp( prepared_timestamps[ spike_data.get_lag() ] ); se_batch[ j ].set_offset( spike_data.get_offset() ); + if ( spike_data.is_activation_marker() ) + { + se_batch[ j ].set_activation(); + } + else + { + se_batch[ j ].unset_activation(); + } syn_id_batch[ j ] = spike_data.get_syn_id(); // for compressed spikes lcid holds the index in the // compressed_spike_data structure @@ -760,7 +801,16 @@ EventDeliveryManager::deliver_events_( const size_t tid, const std::vector< Spik { // non-local sender -> receiver retrieves ID of sender Node from SourceTable based on tid, syn_id, lcid // only if needed, as this is computationally costly + auto activation = se_batch[ j ].get_activation(); se_batch[ j ].set_sender_node_id_info( tid, syn_id_batch[ j ], lcid_batch[ j ] ); + if ( activation ) + { + se_batch[ j ].set_activation(); + } + else + { + se_batch[ j ].unset_activation(); + } } } for ( size_t j = 0; j < num_remaining_entries; ++j ) diff --git a/nestkernel/event_delivery_manager_impl.h b/nestkernel/event_delivery_manager_impl.h index 55310a24bd..8d8b629443 100644 --- a/nestkernel/event_delivery_manager_impl.h +++ b/nestkernel/event_delivery_manager_impl.h @@ -112,7 +112,7 @@ EventDeliveryManager::send_remote( size_t tid, SpikeEvent& e, const long lag ) // Unroll spike multiplicity as plastic synapses only handle individual spikes. for ( size_t i = 0; i < e.get_multiplicity(); ++i ) { - ( *emitted_spikes_register_[ tid ] ).emplace_back( target, lag ); + ( *emitted_spikes_register_[ tid ] ).emplace_back( target, lag, e.get_activation() ); } } } diff --git a/nestkernel/histentry.h b/nestkernel/histentry.h index 7d8985f998..635072ca4c 100644 --- a/nestkernel/histentry.h +++ b/nestkernel/histentry.h @@ -82,13 +82,13 @@ class HistEntryEprop { } - friend bool operator<( const HistEntryEprop& he, long t ); + friend bool operator<( const HistEntryEprop& a, const HistEntryEprop& b ); }; inline bool -operator<( const HistEntryEprop& he, long t ) +operator<( const HistEntryEprop& a, const HistEntryEprop& b ) { - return he.t_ < t; + return a.t_ < b.t_; } /** diff --git a/nestkernel/nest_names.cpp b/nestkernel/nest_names.cpp index 2b1d98435e..22ede20823 100644 --- a/nestkernel/nest_names.cpp +++ b/nestkernel/nest_names.cpp @@ -102,6 +102,7 @@ const Name capacity( "capacity" ); const Name center( "center" ); const Name circular( "circular" ); const Name clear( "clear" ); +const Name activation_interval( "activation_interval" ); const Name comp_idx( "comp_idx" ); const Name comparator( "comparator" ); const Name compartments( "compartments" ); diff --git a/nestkernel/nest_names.h b/nestkernel/nest_names.h index 44a75f62d0..6b704e35fc 100644 --- a/nestkernel/nest_names.h +++ b/nestkernel/nest_names.h @@ -129,6 +129,7 @@ extern const Name capacity; extern const Name center; extern const Name circular; extern const Name clear; +extern const Name activation_interval; extern const Name comp_idx; extern const Name comparator; extern const Name compartments; diff --git a/nestkernel/node.cpp b/nestkernel/node.cpp index 6f54cc1075..bf74f8e0f6 100644 --- a/nestkernel/node.cpp +++ b/nestkernel/node.cpp @@ -218,6 +218,12 @@ Node::register_stdp_connection( double, double ) throw IllegalConnection( "The target node does not support STDP synapses." ); } +void +Node::register_synapse() +{ + throw IllegalConnection( "The target node does not support eprop synapses." ); +} + void Node::register_eprop_connection() { @@ -231,7 +237,13 @@ Node::get_shift() const } void -Node::write_update_to_history( const long, const long, const long ) +Node::write_update_to_history( const long, const long, const bool, const bool, const long ) +{ + throw IllegalConnection( "The target node is not an e-prop neuron." ); +} + +void +Node::erase_used_eprop_history( const long, const long ) { throw IllegalConnection( "The target node is not an e-prop neuron." ); } @@ -559,7 +571,10 @@ nest::Node::compute_gradient( const long, double&, double&, const CommonSynapseProperties&, - WeightOptimizer* ) + WeightOptimizer*, + bool, + bool, + double& ) { throw IllegalConnection( "The target node does not support compute_gradient()." ); } diff --git a/nestkernel/node.h b/nestkernel/node.h index 12742558cc..825869e3b9 100644 --- a/nestkernel/node.h +++ b/nestkernel/node.h @@ -487,6 +487,16 @@ class Node */ virtual void register_stdp_connection( double, double ); + /** + * @brief Registers an eprop synapse and initializes the update history. + * + * The time for the first entry of the update history is set to the neuron specific shift for `bsshslm_2020` + * models and to the negative transmission delay from the recurrent to the output layer otherwise. + * + * @throws IllegalConnection + */ + virtual void register_synapse(); + /** * @brief Registers an eprop synapse and initializes the update history. * @@ -525,8 +535,12 @@ class Node */ virtual void write_update_to_history( const long t_previous_update, const long t_current_update, + const bool activation, + const bool previous_event_was_activation, const long eprop_isi_trace_cutoff = 0 ); + virtual void erase_used_eprop_history( const long t_spike, const long t_spike_previous ); + /** * Retrieves the maximum number of time steps integrated between two consecutive spikes. * @@ -858,7 +872,10 @@ class Node double& epsilon, double& weight, const CommonSynapseProperties& cp, - WeightOptimizer* optimizer ); + WeightOptimizer* optimizer, + bool activation, + bool previous_event_was_activation, + double& sum_grad ); /** * Compute gradient change for eprop synapses. diff --git a/nestkernel/spike_data.h b/nestkernel/spike_data.h index 87b1205ba0..4a1f678acd 100644 --- a/nestkernel/spike_data.h +++ b/nestkernel/spike_data.h @@ -85,6 +85,12 @@ enum enum_status_spike_data_id SPIKE_DATA_ID_INVALID }; +enum enum_activation_id +{ + SPIKE_DATA_ID_ACTIVATION_DEFAULT, + SPIKE_DATA_ID_ACTIVATION +}; + /** * Used to communicate spikes. These are the elements of the MPI * buffers. @@ -96,16 +102,17 @@ class SpikeData protected: static constexpr int MAX_LAG = generate_max_value( NUM_BITS_LAG ); - size_t lcid_ : NUM_BITS_LCID; //!< local connection index - unsigned int marker_ : NUM_BITS_MARKER_SPIKE_DATA; //!< status flag - unsigned int lag_ : NUM_BITS_LAG; //!< lag in this min-delay interval - unsigned int tid_ : NUM_BITS_TID; //!< thread index - synindex syn_id_ : NUM_BITS_SYN_ID; //!< synapse-type index + size_t lcid_ : NUM_BITS_LCID; //!< local connection index + unsigned int marker_ : NUM_BITS_MARKER_SPIKE_DATA; //!< status flag + unsigned int activation_marker_ : NUM_BITS_MARKER_ACTIVATION; //!< activation flag + unsigned int lag_ : NUM_BITS_LAG; //!< lag in this min-delay interval + unsigned int tid_ : NUM_BITS_TID; //!< thread index + synindex syn_id_ : NUM_BITS_SYN_ID; //!< synapse-type index public: SpikeData(); SpikeData( const SpikeData& rhs ); - SpikeData( const Target& target, const size_t lag ); + SpikeData( const Target& target, const size_t lag, const bool activation ); SpikeData( const size_t tid, const synindex syn_id, const size_t lcid, const unsigned int lag ); SpikeData& operator=( const SpikeData& rhs ); @@ -183,10 +190,25 @@ class SpikeData */ bool is_invalid_marker() const; + /** + * Returns whether the marker is the default marker. + */ + bool is_activation_marker() const; + /** * Returns offset. */ double get_offset() const; + + /** + * Sets the activation marker. + */ + void set_activation_marker(); + + /** + * Unsets the activation marker. + */ + void unset_activation_marker(); }; //! check legal size @@ -195,6 +217,7 @@ using success_spike_data_size = StaticAssert< sizeof( SpikeData ) == 8 >::succes inline SpikeData::SpikeData() : lcid_( 0 ) , marker_( SPIKE_DATA_ID_DEFAULT ) + , activation_marker_( SPIKE_DATA_ID_ACTIVATION_DEFAULT ) , lag_( 0 ) , tid_( 0 ) , syn_id_( 0 ) @@ -204,24 +227,31 @@ inline SpikeData::SpikeData() inline SpikeData::SpikeData( const SpikeData& rhs ) : lcid_( rhs.lcid_ ) , marker_( rhs.marker_ ) + , activation_marker_( rhs.activation_marker_ ) , lag_( rhs.lag_ ) , tid_( rhs.tid_ ) , syn_id_( rhs.syn_id_ ) { } -inline SpikeData::SpikeData( const Target& target, const size_t lag ) +inline SpikeData::SpikeData( const Target& target, const size_t lag, const bool activation ) : lcid_( target.get_lcid() ) , marker_( SPIKE_DATA_ID_DEFAULT ) + , activation_marker_( SPIKE_DATA_ID_ACTIVATION_DEFAULT ) , lag_( lag ) , tid_( target.get_tid() ) , syn_id_( target.get_syn_id() ) { + if ( activation ) + { + activation_marker_ = SPIKE_DATA_ID_ACTIVATION; + } } inline SpikeData::SpikeData( const size_t tid, const synindex syn_id, const size_t lcid, const unsigned int lag ) : lcid_( lcid ) , marker_( SPIKE_DATA_ID_DEFAULT ) + , activation_marker_( SPIKE_DATA_ID_ACTIVATION_DEFAULT ) , lag_( lag ) , tid_( tid ) , syn_id_( syn_id ) @@ -233,6 +263,7 @@ SpikeData::operator=( const SpikeData& rhs ) { lcid_ = rhs.lcid_; marker_ = rhs.marker_; + activation_marker_ = rhs.activation_marker_; lag_ = rhs.lag_; tid_ = rhs.tid_; syn_id_ = rhs.syn_id_; @@ -249,6 +280,7 @@ SpikeData::set( const size_t tid, const synindex syn_id, const size_t lcid, cons lcid_ = lcid; marker_ = SPIKE_DATA_ID_DEFAULT; + activation_marker_ = SPIKE_DATA_ID_ACTIVATION_DEFAULT; lag_ = lag; tid_ = tid; syn_id_ = syn_id; @@ -263,6 +295,7 @@ SpikeData::set( const TargetT& target, const unsigned int lag ) assert( lag < MAX_LAG ); lcid_ = target.get_lcid(); marker_ = SPIKE_DATA_ID_DEFAULT; + activation_marker_ = SPIKE_DATA_ID_ACTIVATION_DEFAULT; lag_ = lag; tid_ = target.get_tid(); syn_id_ = target.get_syn_id(); @@ -347,12 +380,30 @@ SpikeData::is_invalid_marker() const return marker_ == SPIKE_DATA_ID_INVALID; } +inline bool +SpikeData::is_activation_marker() const +{ + return activation_marker_ == SPIKE_DATA_ID_ACTIVATION; +} + inline double SpikeData::get_offset() const { return 0; } +inline void +SpikeData::set_activation_marker() +{ + activation_marker_ = SPIKE_DATA_ID_ACTIVATION; +} + +inline void +SpikeData::unset_activation_marker() +{ + activation_marker_ = SPIKE_DATA_ID_ACTIVATION_DEFAULT; +} + class OffGridSpikeData : public SpikeData { private: @@ -386,7 +437,7 @@ inline OffGridSpikeData::OffGridSpikeData() } inline OffGridSpikeData::OffGridSpikeData( const Target& target, const size_t lag, const double offset ) - : SpikeData( target, lag ) + : SpikeData( target, lag, false ) , offset_( offset ) { } @@ -447,6 +498,7 @@ OffGridSpikeData::set( const size_t tid, lcid_ = lcid; marker_ = SPIKE_DATA_ID_DEFAULT; + activation_marker_ = SPIKE_DATA_ID_ACTIVATION_DEFAULT; lag_ = lag; tid_ = tid; syn_id_ = syn_id; @@ -476,15 +528,15 @@ OffGridSpikeData::get_offset() const */ struct SpikeDataWithRank { - SpikeDataWithRank( const Target& target, const size_t lag ); + SpikeDataWithRank( const Target& target, const size_t lag, const bool activation ); const size_t rank; //!< rank of target neuron const SpikeData spike_data; //! data on spike transmitted }; -inline SpikeDataWithRank::SpikeDataWithRank( const Target& target, const size_t lag ) +inline SpikeDataWithRank::SpikeDataWithRank( const Target& target, const size_t lag, const bool activation ) : rank( target.get_rank() ) - , spike_data( target, lag ) + , spike_data( target, lag, activation ) { } diff --git a/pynest/examples/eprop_plasticity/eprop_supervised_classification_evidence-accumulation_bsshslm_2020.py b/pynest/examples/eprop_plasticity/eprop_supervised_classification_evidence-accumulation_bsshslm_2020.py index fc734ed8d4..d0b9cc454c 100644 --- a/pynest/examples/eprop_plasticity/eprop_supervised_classification_evidence-accumulation_bsshslm_2020.py +++ b/pynest/examples/eprop_plasticity/eprop_supervised_classification_evidence-accumulation_bsshslm_2020.py @@ -263,7 +263,7 @@ # since devices cannot establish plastic synapses for technical reasons gen_spk_in = nest.Create("spike_generator", n_in) -nrns_in = nest.Create("parrot_neuron", n_in) +nrns_in = nest.Create("eprop_input_neuron", n_in) # 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. diff --git a/pynest/examples/eprop_plasticity/eprop_supervised_classification_neuromorphic_mnist.py b/pynest/examples/eprop_plasticity/eprop_supervised_classification_neuromorphic_mnist.py index 84a0b990bf..30a6c16c52 100644 --- a/pynest/examples/eprop_plasticity/eprop_supervised_classification_neuromorphic_mnist.py +++ b/pynest/examples/eprop_plasticity/eprop_supervised_classification_neuromorphic_mnist.py @@ -244,7 +244,7 @@ # since devices cannot establish plastic synapses for technical reasons gen_spk_in = nest.Create("spike_generator", n_in) -nrns_in = nest.Create("parrot_neuron", n_in) +nrns_in = nest.Create("eprop_input_neuron", n_in) nrns_rec = nest.Create(model_nrn_rec, n_rec, params_nrn_rec) nrns_out = nest.Create("eprop_readout", n_out, params_nrn_out) diff --git a/pynest/examples/eprop_plasticity/eprop_supervised_regression_handwriting_bsshslm_2020.py b/pynest/examples/eprop_plasticity/eprop_supervised_regression_handwriting_bsshslm_2020.py index ef2fa7a51b..4d854eb339 100644 --- a/pynest/examples/eprop_plasticity/eprop_supervised_regression_handwriting_bsshslm_2020.py +++ b/pynest/examples/eprop_plasticity/eprop_supervised_regression_handwriting_bsshslm_2020.py @@ -233,7 +233,7 @@ # since devices cannot establish plastic synapses for technical reasons gen_spk_in = nest.Create("spike_generator", n_in) -nrns_in = nest.Create("parrot_neuron", n_in) +nrns_in = nest.Create("eprop_input_neuron", n_in) # 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. diff --git a/pynest/examples/eprop_plasticity/eprop_supervised_regression_lemniscate_bsshslm_2020.py b/pynest/examples/eprop_plasticity/eprop_supervised_regression_lemniscate_bsshslm_2020.py index 3b42282ae3..53668561be 100644 --- a/pynest/examples/eprop_plasticity/eprop_supervised_regression_lemniscate_bsshslm_2020.py +++ b/pynest/examples/eprop_plasticity/eprop_supervised_regression_lemniscate_bsshslm_2020.py @@ -219,7 +219,7 @@ # since devices cannot establish plastic synapses for technical reasons gen_spk_in = nest.Create("spike_generator", n_in) -nrns_in = nest.Create("parrot_neuron", n_in) +nrns_in = nest.Create("eprop_input_neuron", n_in) # 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. diff --git a/pynest/examples/eprop_plasticity/eprop_supervised_regression_sine-waves.py b/pynest/examples/eprop_plasticity/eprop_supervised_regression_sine-waves.py index bf0a66e992..0e0235fdbf 100644 --- a/pynest/examples/eprop_plasticity/eprop_supervised_regression_sine-waves.py +++ b/pynest/examples/eprop_plasticity/eprop_supervised_regression_sine-waves.py @@ -217,7 +217,7 @@ # since devices cannot establish plastic synapses for technical reasons gen_spk_in = nest.Create("spike_generator", n_in) -nrns_in = nest.Create("parrot_neuron", n_in) +nrns_in = nest.Create("eprop_input_neuron", n_in) nrns_rec = nest.Create(model_nrn_rec, n_rec, params_nrn_rec) nrns_out = nest.Create("eprop_readout", n_out, params_nrn_out) diff --git a/pynest/examples/eprop_plasticity/eprop_supervised_regression_sine-waves_bsshslm_2020.py b/pynest/examples/eprop_plasticity/eprop_supervised_regression_sine-waves_bsshslm_2020.py index 46387e4df5..cf667efc93 100644 --- a/pynest/examples/eprop_plasticity/eprop_supervised_regression_sine-waves_bsshslm_2020.py +++ b/pynest/examples/eprop_plasticity/eprop_supervised_regression_sine-waves_bsshslm_2020.py @@ -208,7 +208,7 @@ # since devices cannot establish plastic synapses for technical reasons gen_spk_in = nest.Create("spike_generator", n_in) -nrns_in = nest.Create("parrot_neuron", n_in) +nrns_in = nest.Create("eprop_input_neuron", n_in) # 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. diff --git a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py index 0ed52cac1f..d5c92999bd 100644 --- a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py +++ b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py @@ -55,6 +55,8 @@ "rate_transformer_sigmoid_gg_1998", # rate transformer "parrot_neuron", "parrot_neuron_ps", + "eprop_input_neuron", + "eprop_input_neuron_ps", "spike_train_injector", # spike emitting neuron, does not support spike input "cm_default", # cannot readout V_m directly "iaf_cond_alpha_mc", # cannot readout V_m directly diff --git a/testsuite/pytests/test_eprop_bsshslm_2020_plasticity.py b/testsuite/pytests/test_eprop_bsshslm_2020_plasticity.py index 69ed7d8fff..f53d8a28c5 100644 --- a/testsuite/pytests/test_eprop_bsshslm_2020_plasticity.py +++ b/testsuite/pytests/test_eprop_bsshslm_2020_plasticity.py @@ -171,7 +171,7 @@ def test_eprop_regression(): params_nrn_rec["beta"] /= np.abs(params_nrn_rec["V_th"]) gen_spk_in = nest.Create("spike_generator", n_in) - nrns_in = nest.Create("parrot_neuron", n_in) + nrns_in = nest.Create("eprop_input_neuron", n_in) nrns_rec = nest.Create("eprop_iaf_bsshslm_2020", n_rec, params_nrn_rec) nrns_out = nest.Create("eprop_readout_bsshslm_2020", n_out, params_nrn_out) gen_rate_target = nest.Create("step_rate_generator", n_out) @@ -513,6 +513,7 @@ def test_eprop_classification(batch_size, loss_nest_reference): "tau_m": 20.0, "V_m": 0.0, "V_th": 0.6, + "activation_interval": 3 * duration["sequence"], } params_nrn_reg["gamma"] /= params_nrn_reg["V_th"] @@ -534,6 +535,7 @@ def test_eprop_classification(batch_size, loss_nest_reference): "tau_m": 20.0, "V_m": 0.0, "V_th": 0.6, + "activation_interval": 3 * duration["sequence"], } params_nrn_ad["gamma"] /= params_nrn_ad["V_th"] @@ -545,7 +547,7 @@ def test_eprop_classification(batch_size, loss_nest_reference): ) gen_spk_in = nest.Create("spike_generator", n_in) - nrns_in = nest.Create("parrot_neuron", n_in) + nrns_in = nest.Create("eprop_input_neuron", n_in) nrns_reg = nest.Create("eprop_iaf_bsshslm_2020", n_reg, params_nrn_reg) nrns_ad = nest.Create("eprop_iaf_adapt_bsshslm_2020", n_ad, params_nrn_ad) nrns_out = nest.Create("eprop_readout_bsshslm_2020", n_out, params_nrn_out) @@ -861,7 +863,7 @@ def test_eprop_history_cleaning(neuron_model, eprop_history_duration_reference): # Create neurons gen_spk_in = nest.Create("spike_generator", 3) - nrns_in = nest.Create("parrot_neuron", 3) + nrns_in = nest.Create("eprop_input_neuron", 3) nrns_rec = nest.Create(neuron_model, 1) # Create recorders diff --git a/testsuite/pytests/test_eprop_plasticity.py b/testsuite/pytests/test_eprop_plasticity.py index 6bb51436ef..98705bfbbf 100644 --- a/testsuite/pytests/test_eprop_plasticity.py +++ b/testsuite/pytests/test_eprop_plasticity.py @@ -214,7 +214,7 @@ def test_eprop_regression(neuron_model, optimizer, loss_nest_reference): params_nrn_rec["adaptation"] = 0.0 gen_spk_in = nest.Create("spike_generator", n_in) - nrns_in = nest.Create("parrot_neuron", n_in) + nrns_in = nest.Create("eprop_input_neuron", n_in) nrns_rec = nest.Create(neuron_model, n_rec, params_nrn_rec) nrns_out = nest.Create("eprop_readout", n_out, params_nrn_out) gen_rate_target = nest.Create("step_rate_generator", n_out) @@ -523,7 +523,7 @@ def test_eprop_surrogate_gradients(surrogate_gradient_type, surrogate_gradient_r } gen_spk_in = nest.Create("spike_generator", 1) - nrns_in = nest.Create("parrot_neuron", 1) + nrns_in = nest.Create("eprop_input_neuron", 1) nrns_rec = nest.Create("eprop_iaf", 1, params_nrn_rec) params_mm_rec = { @@ -558,23 +558,9 @@ def test_eprop_surrogate_gradients(surrogate_gradient_type, surrogate_gradient_r @pytest.mark.parametrize( - "neuron_model,eprop_isi_trace_cutoff,eprop_history_duration_reference", - [ - ( - "eprop_iaf", - 5.0, - np.hstack( - [ - np.arange(x, y) - for x, y in [[1, 12], [6, 16], [11, 21], [16, 26], [21, 31], [17, 27], [12, 42], [12, 30]] - ] - ), - ), - ("eprop_iaf", 100000.0, np.hstack([np.arange(x, y) for x, y in [[1, 52], [33, 43], [23, 53], [43, 61]]])), - ("eprop_readout", 100000.0, np.hstack([np.arange(x, y) for x, y in [[1, 52], [33, 43], [23, 53], [43, 61]]])), - ], + "neuron_model,eprop_isi_trace_cutoff", [("eprop_iaf", 5.0), ("eprop_iaf", 100000.0), ("eprop_readout", 100000.0)] ) -def test_eprop_history_cleaning(neuron_model, eprop_isi_trace_cutoff, eprop_history_duration_reference): +def test_eprop_history_cleaning(neuron_model, eprop_isi_trace_cutoff): """ Test the e-prop archiving mechanism's cleaning process by ensuring that the length of the `eprop_history` buffer matches the expected values based on a given input firing pattern. These reference length values @@ -599,13 +585,16 @@ def test_eprop_history_cleaning(neuron_model, eprop_isi_trace_cutoff, eprop_hist # Create neurons - params_nrn_rec = { + params_nrn = { "eprop_isi_trace_cutoff": eprop_isi_trace_cutoff, } + if neuron_model != "eprop_readout": + params_nrn["activation_interval"] = eprop_isi_trace_cutoff + gen_spk_in = nest.Create("spike_generator", 3) - nrns_in = nest.Create("parrot_neuron", 3) - nrns_rec = nest.Create(neuron_model, 1, params_nrn_rec) + nrns_in = nest.Create("eprop_input_neuron", 3) + nrns = nest.Create(neuron_model, 1, params_nrn) # Create recorders @@ -635,8 +624,8 @@ def test_eprop_history_cleaning(neuron_model, eprop_isi_trace_cutoff, eprop_hist params_syn_in = params_syn_base.copy() nest.Connect(gen_spk_in, nrns_in, params_conn_one_to_one, params_syn_static) - nest.Connect(nrns_in, nrns_rec, params_conn_all_to_all, params_syn_in) - nest.Connect(mm_rec, nrns_rec, params_conn_all_to_all, params_syn_static) + nest.Connect(nrns_in, nrns, params_conn_all_to_all, params_syn_in) + nest.Connect(mm_rec, nrns, params_conn_all_to_all, params_syn_static) # Create input @@ -659,7 +648,11 @@ def test_eprop_history_cleaning(neuron_model, eprop_isi_trace_cutoff, eprop_hist events_mm_rec = mm_rec.get("events") eprop_history_duration = events_mm_rec["eprop_history_duration"] - senders = events_mm_rec["senders"] - eprop_history_duration = np.array([eprop_history_duration[senders == i] for i in set(senders)])[0] + eprop_history_duration_reference = np.hstack( + [ + np.arange(x, y + 1) + for x, y in [[1.0, 11.0], [2.0, 21.0], [12.0, 31.0], [12.0, 21.0], [12.0, 41.0], [32.0, 49.0]] + ] + ) assert np.allclose(eprop_history_duration, eprop_history_duration_reference, rtol=1e-8)