diff --git a/models/iaf_psc_alpha_hom.cpp b/models/iaf_psc_alpha_hom.cpp index 68135383d7..f9e5b03a01 100644 --- a/models/iaf_psc_alpha_hom.cpp +++ b/models/iaf_psc_alpha_hom.cpp @@ -1,24 +1,24 @@ /* -* iaf_psc_alpha_hom.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 . -* -*/ + * iaf_psc_alpha_hom.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 "iaf_psc_alpha_hom.h" @@ -45,364 +45,364 @@ namespace nest void register_iaf_psc_alpha_hom( const std::string& name ) { - register_node_model< iaf_psc_alpha_hom >( name ); + register_node_model< iaf_psc_alpha_hom >( name ); } /* -* Override the create() method with one call to RecordablesMap::insert_() -* for each quantity to be recorded. -*/ + * Override the create() method with one call to RecordablesMap::insert_() + * for each quantity to be recorded. + */ template <> void RecordablesMap< iaf_psc_alpha_hom >::create() { - // use standard names wherever you can for consistency! - insert_( names::V_m, &iaf_psc_alpha_hom::get_V_m_ ); - insert_( names::I_syn_ex, &iaf_psc_alpha_hom::get_I_syn_ex_ ); - insert_( names::I_syn_in, &iaf_psc_alpha_hom::get_I_syn_in_ ); + // use standard names wherever you can for consistency! + insert_( names::V_m, &iaf_psc_alpha_hom::get_V_m_ ); + insert_( names::I_syn_ex, &iaf_psc_alpha_hom::get_I_syn_ex_ ); + insert_( names::I_syn_in, &iaf_psc_alpha_hom::get_I_syn_in_ ); } /* ---------------------------------------------------------------- -* Default constructors defining default parameters and state -* ---------------------------------------------------------------- */ + * Default constructors defining default parameters and state + * ---------------------------------------------------------------- */ iaf_psc_alpha_hom::Parameters_::Parameters_() - : Tau_( 10.0 ) // ms - , C_( 250.0 ) // pF - , TauR_( 2.0 ) // ms - , E_L_( -70.0 ) // mV - , I_e_( 0.0 ) // pA - , V_reset_( -70.0 - E_L_ ) // mV, rel to E_L_ - , Theta_( -55.0 - E_L_ ) // mV, rel to E_L_ - , LowerBound_( -std::numeric_limits< double >::infinity() ) - , tau_ex_( 2.0 ) // ms - , tau_in_( 2.0 ) // ms + : Tau_( 10.0 ) // ms + , C_( 250.0 ) // pF + , TauR_( 2.0 ) // ms + , E_L_( -70.0 ) // mV + , I_e_( 0.0 ) // pA + , V_reset_( -70.0 - E_L_ ) // mV, rel to E_L_ + , Theta_( -55.0 - E_L_ ) // mV, rel to E_L_ + , LowerBound_( -std::numeric_limits< double >::infinity() ) + , tau_ex_( 2.0 ) // ms + , tau_in_( 2.0 ) // ms { } iaf_psc_alpha_hom::State_::State_() - : y0_( 0.0 ) - , dI_ex_( 0.0 ) - , I_ex_( 0.0 ) - , dI_in_( 0.0 ) - , I_in_( 0.0 ) - , y3_( 0.0 ) - , r_( 0 ) + : y0_( 0.0 ) + , dI_ex_( 0.0 ) + , I_ex_( 0.0 ) + , dI_in_( 0.0 ) + , I_in_( 0.0 ) + , y3_( 0.0 ) + , r_( 0 ) { } /* ---------------------------------------------------------------- -* Parameter and state extractions and manipulation functions -* ---------------------------------------------------------------- */ + * Parameter and state extractions and manipulation functions + * ---------------------------------------------------------------- */ void iaf_psc_alpha_hom::Parameters_::get( DictionaryDatum& d ) const { - def< double >( d, names::E_L, E_L_ ); // Resting potential - def< double >( d, names::I_e, I_e_ ); - def< double >( d, names::V_th, Theta_ + E_L_ ); // threshold value - def< double >( d, names::V_reset, V_reset_ + E_L_ ); - def< double >( d, names::V_min, LowerBound_ + E_L_ ); - def< double >( d, names::C_m, C_ ); - def< double >( d, names::tau_m, Tau_ ); - def< double >( d, names::t_ref, TauR_ ); - def< double >( d, names::tau_syn_ex, tau_ex_ ); - def< double >( d, names::tau_syn_in, tau_in_ ); + def< double >( d, names::E_L, E_L_ ); // Resting potential + def< double >( d, names::I_e, I_e_ ); + def< double >( d, names::V_th, Theta_ + E_L_ ); // threshold value + def< double >( d, names::V_reset, V_reset_ + E_L_ ); + def< double >( d, names::V_min, LowerBound_ + E_L_ ); + def< double >( d, names::C_m, C_ ); + def< double >( d, names::tau_m, Tau_ ); + def< double >( d, names::t_ref, TauR_ ); + def< double >( d, names::tau_syn_ex, tau_ex_ ); + def< double >( d, names::tau_syn_in, tau_in_ ); } double iaf_psc_alpha_hom::Parameters_::set( const DictionaryDatum& d, Node* node ) { - // if E_L_ is changed, we need to adjust all variables defined relative to - // E_L_ - const double ELold = E_L_; - updateValueParam< double >( d, names::E_L, E_L_, node ); - const double delta_EL = E_L_ - ELold; - - if ( updateValueParam< double >( d, names::V_reset, V_reset_, node ) ) - { - V_reset_ -= E_L_; - } - else - { - V_reset_ -= delta_EL; - } - - if ( updateValueParam< double >( d, names::V_th, Theta_, node ) ) - { - Theta_ -= E_L_; - } - else - { - Theta_ -= delta_EL; - } - - if ( updateValueParam< double >( d, names::V_min, LowerBound_, node ) ) - { - LowerBound_ -= E_L_; - } - else - { - LowerBound_ -= delta_EL; - } - - updateValueParam< double >( d, names::I_e, I_e_, node ); - updateValueParam< double >( d, names::C_m, C_, node ); - updateValueParam< double >( d, names::tau_m, Tau_, node ); - updateValueParam< double >( d, names::tau_syn_ex, tau_ex_, node ); - updateValueParam< double >( d, names::tau_syn_in, tau_in_, node ); - updateValueParam< double >( d, names::t_ref, TauR_, node ); - - if ( C_ <= 0.0 ) - { - throw BadProperty( "Capacitance must be > 0." ); - } - - if ( Tau_ <= 0.0 ) - { - throw BadProperty( "Membrane time constant must be > 0." ); - } - - if ( tau_ex_ <= 0.0 or tau_in_ <= 0.0 ) - { - throw BadProperty( "All synaptic time constants must be > 0." ); - } - - if ( TauR_ < 0.0 ) - { - throw BadProperty( "The refractory time t_ref can't be negative." ); - } - if ( V_reset_ >= Theta_ ) - { - throw BadProperty( "Reset potential must be smaller than threshold." ); - } - - return delta_EL; + // if E_L_ is changed, we need to adjust all variables defined relative to + // E_L_ + const double ELold = E_L_; + updateValueParam< double >( d, names::E_L, E_L_, node ); + const double delta_EL = E_L_ - ELold; + + if ( updateValueParam< double >( d, names::V_reset, V_reset_, node ) ) + { + V_reset_ -= E_L_; + } + else + { + V_reset_ -= delta_EL; + } + + if ( updateValueParam< double >( d, names::V_th, Theta_, node ) ) + { + Theta_ -= E_L_; + } + else + { + Theta_ -= delta_EL; + } + + if ( updateValueParam< double >( d, names::V_min, LowerBound_, node ) ) + { + LowerBound_ -= E_L_; + } + else + { + LowerBound_ -= delta_EL; + } + + updateValueParam< double >( d, names::I_e, I_e_, node ); + updateValueParam< double >( d, names::C_m, C_, node ); + updateValueParam< double >( d, names::tau_m, Tau_, node ); + updateValueParam< double >( d, names::tau_syn_ex, tau_ex_, node ); + updateValueParam< double >( d, names::tau_syn_in, tau_in_, node ); + updateValueParam< double >( d, names::t_ref, TauR_, node ); + + if ( C_ <= 0.0 ) + { + throw BadProperty( "Capacitance must be > 0." ); + } + + if ( Tau_ <= 0.0 ) + { + throw BadProperty( "Membrane time constant must be > 0." ); + } + + if ( tau_ex_ <= 0.0 or tau_in_ <= 0.0 ) + { + throw BadProperty( "All synaptic time constants must be > 0." ); + } + + if ( TauR_ < 0.0 ) + { + throw BadProperty( "The refractory time t_ref can't be negative." ); + } + if ( V_reset_ >= Theta_ ) + { + throw BadProperty( "Reset potential must be smaller than threshold." ); + } + + return delta_EL; } void iaf_psc_alpha_hom::State_::get( DictionaryDatum& d, const Parameters_& p ) const { - def< double >( d, names::V_m, y3_ + p.E_L_ ); // Membrane potential + def< double >( d, names::V_m, y3_ + p.E_L_ ); // Membrane potential } void iaf_psc_alpha_hom::State_::set( const DictionaryDatum& d, const Parameters_& p, double delta_EL, Node* node ) { - if ( updateValueParam< double >( d, names::V_m, y3_, node ) ) - { - y3_ -= p.E_L_; - } - else - { - y3_ -= delta_EL; - } + if ( updateValueParam< double >( d, names::V_m, y3_, node ) ) + { + y3_ -= p.E_L_; + } + else + { + y3_ -= delta_EL; + } } iaf_psc_alpha_hom::Buffers_::Buffers_( iaf_psc_alpha_hom& n ) - : logger_( n ) + : logger_( n ) { } iaf_psc_alpha_hom::Buffers_::Buffers_( const Buffers_&, iaf_psc_alpha_hom& n ) - : logger_( n ) + : logger_( n ) { } /* ---------------------------------------------------------------- -* Default and copy constructor for node -* ---------------------------------------------------------------- */ + * Default and copy constructor for node + * ---------------------------------------------------------------- */ iaf_psc_alpha_hom::iaf_psc_alpha_hom() - : ArchivingNodeHom() - , P_() - , S_() - , B_( *this ) + : ArchivingNodeHom() + , P_() + , S_() + , B_( *this ) { - recordablesMap_.create(); + recordablesMap_.create(); } iaf_psc_alpha_hom::iaf_psc_alpha_hom( const iaf_psc_alpha_hom& n ) - : ArchivingNodeHom( n ) - , P_( n.P_ ) - , S_( n.S_ ) - , B_( n.B_, *this ) + : ArchivingNodeHom( n ) + , P_( n.P_ ) + , S_( n.S_ ) + , B_( n.B_, *this ) { } /* ---------------------------------------------------------------- -* Node initialization functions -* ---------------------------------------------------------------- */ + * Node initialization functions + * ---------------------------------------------------------------- */ void iaf_psc_alpha_hom::init_buffers_() { - B_.input_buffer_.clear(); // includes resize + B_.input_buffer_.clear(); // includes resize - B_.logger_.reset(); + B_.logger_.reset(); - ArchivingNodeHom::clear_history(); + ArchivingNodeHom::clear_history(); } void iaf_psc_alpha_hom::pre_run_hook() { - // ensures initialization in case mm connected after Simulate - B_.logger_.init(); - - const double h = Time::get_resolution().get_ms(); - - // these P are independent - V_.P11_ex_ = V_.P22_ex_ = std::exp( -h / P_.tau_ex_ ); - V_.P11_in_ = V_.P22_in_ = std::exp( -h / P_.tau_in_ ); - - V_.P33_ = std::exp( -h / P_.Tau_ ); - - V_.expm1_tau_m_ = numerics::expm1( -h / P_.Tau_ ); - - // these depend on the above. Please do not change the order. - V_.P30_ = -P_.Tau_ / P_.C_ * numerics::expm1( -h / P_.Tau_ ); - V_.P21_ex_ = h * V_.P11_ex_; - V_.P21_in_ = h * V_.P11_in_; - - // these are determined according to a numeric stability criterion - std::tie( V_.P31_ex_, V_.P32_ex_ ) = IAFPropagatorAlpha( P_.tau_ex_, P_.Tau_, P_.C_ ).evaluate( h ); - std::tie( V_.P31_in_, V_.P32_in_ ) = IAFPropagatorAlpha( P_.tau_in_, P_.Tau_, P_.C_ ).evaluate( h ); - - V_.EPSCInitialValue_ = 1.0 * numerics::e / P_.tau_ex_; - V_.IPSCInitialValue_ = 1.0 * numerics::e / P_.tau_in_; - - // TauR specifies the length of the absolute refractory period as - // a double in ms. The grid based iaf_psc_alpha can only handle refractory - // periods that are integer multiples of the computation step size (h). - // To ensure consistency with the overall simulation scheme such conversion - // should be carried out via objects of class nest::Time. The conversion - // requires 2 steps: - // 1. A time object is constructed defining representation of - // TauR in tics. This representation is then converted to computation - // time steps again by a strategy defined by class nest::Time. - // 2. The refractory time in units of steps is read out get_steps(), a - // member function of class nest::Time. - // - // The definition of the refractory period of the iaf_psc_alpha is consistent - // the one of iaf_psc_alpha_ps. - // - // Choosing a TauR that is not an integer multiple of the computation time - // step h will lead to accurate (up to the resolution h) and self-consistent - // results. However, a neuron model capable of operating with real valued - // spike time may exhibit a different effective refractory time. - - V_.RefractoryCounts_ = Time( Time::ms( P_.TauR_ ) ).get_steps(); - // since t_ref_ >= 0, this can only fail in error - assert( V_.RefractoryCounts_ >= 0 ); + // ensures initialization in case mm connected after Simulate + B_.logger_.init(); + + const double h = Time::get_resolution().get_ms(); + + // these P are independent + V_.P11_ex_ = V_.P22_ex_ = std::exp( -h / P_.tau_ex_ ); + V_.P11_in_ = V_.P22_in_ = std::exp( -h / P_.tau_in_ ); + + V_.P33_ = std::exp( -h / P_.Tau_ ); + + V_.expm1_tau_m_ = numerics::expm1( -h / P_.Tau_ ); + + // these depend on the above. Please do not change the order. + V_.P30_ = -P_.Tau_ / P_.C_ * numerics::expm1( -h / P_.Tau_ ); + V_.P21_ex_ = h * V_.P11_ex_; + V_.P21_in_ = h * V_.P11_in_; + + // these are determined according to a numeric stability criterion + std::tie( V_.P31_ex_, V_.P32_ex_ ) = IAFPropagatorAlpha( P_.tau_ex_, P_.Tau_, P_.C_ ).evaluate( h ); + std::tie( V_.P31_in_, V_.P32_in_ ) = IAFPropagatorAlpha( P_.tau_in_, P_.Tau_, P_.C_ ).evaluate( h ); + + V_.EPSCInitialValue_ = 1.0 * numerics::e / P_.tau_ex_; + V_.IPSCInitialValue_ = 1.0 * numerics::e / P_.tau_in_; + + // TauR specifies the length of the absolute refractory period as + // a double in ms. The grid based iaf_psc_alpha can only handle refractory + // periods that are integer multiples of the computation step size (h). + // To ensure consistency with the overall simulation scheme such conversion + // should be carried out via objects of class nest::Time. The conversion + // requires 2 steps: + // 1. A time object is constructed defining representation of + // TauR in tics. This representation is then converted to computation + // time steps again by a strategy defined by class nest::Time. + // 2. The refractory time in units of steps is read out get_steps(), a + // member function of class nest::Time. + // + // The definition of the refractory period of the iaf_psc_alpha is consistent + // the one of iaf_psc_alpha_ps. + // + // Choosing a TauR that is not an integer multiple of the computation time + // step h will lead to accurate (up to the resolution h) and self-consistent + // results. However, a neuron model capable of operating with real valued + // spike time may exhibit a different effective refractory time. + + V_.RefractoryCounts_ = Time( Time::ms( P_.TauR_ ) ).get_steps(); + // since t_ref_ >= 0, this can only fail in error + assert( V_.RefractoryCounts_ >= 0 ); } /* ---------------------------------------------------------------- -* Update and spike handling functions -*/ + * Update and spike handling functions + */ void iaf_psc_alpha_hom::update( Time const& origin, const long from, const long to ) { - for ( long lag = from; lag < to; ++lag ) - { - if ( S_.r_ == 0 ) - { - // neuron not refractory - S_.y3_ = V_.P30_ * ( S_.y0_ + P_.I_e_ ) + V_.P31_ex_ * S_.dI_ex_ + V_.P32_ex_ * S_.I_ex_ + V_.P31_in_ * S_.dI_in_ - + V_.P32_in_ * S_.I_in_ + V_.expm1_tau_m_ * S_.y3_ + S_.y3_; - - // lower bound of membrane potential - S_.y3_ = ( S_.y3_ < P_.LowerBound_ ? P_.LowerBound_ : S_.y3_ ); - } - else - { - // neuron is absolute refractory - --S_.r_; - } - - // alpha shape EPSCs - S_.I_ex_ = V_.P21_ex_ * S_.dI_ex_ + V_.P22_ex_ * S_.I_ex_; - S_.dI_ex_ *= V_.P11_ex_; - - // get read access to the correct input-buffer slot - const size_t input_buffer_slot = kernel().event_delivery_manager.get_modulo( lag ); - auto& input = B_.input_buffer_.get_values_all_channels( input_buffer_slot ); - - // Apply spikes delivered in this step; spikes arriving at T+1 have - // an immediate effect on the state of the neuron - V_.weighted_spikes_ex_ = input[ Buffers_::SYN_EX ]; - S_.dI_ex_ += V_.EPSCInitialValue_ * V_.weighted_spikes_ex_; - - // alpha shape EPSCs - S_.I_in_ = V_.P21_in_ * S_.dI_in_ + V_.P22_in_ * S_.I_in_; - S_.dI_in_ *= V_.P11_in_; - - // Apply spikes delivered in this step; spikes arriving at T+1 have - // an immediate effect on the state of the neuron - V_.weighted_spikes_in_ = input[ Buffers_::SYN_IN ]; - S_.dI_in_ += V_.IPSCInitialValue_ * V_.weighted_spikes_in_; - - // threshold crossing - if ( S_.y3_ >= P_.Theta_ ) - { - S_.r_ = V_.RefractoryCounts_; - S_.y3_ = P_.V_reset_; - // A supra-threshold membrane potential should never be observable. - // The reset at the time of threshold crossing enables accurate - // integration independent of the computation step size, see [2,3] for - // details. - - set_spiketime( Time::step( origin.get_steps() + lag + 1 ) ); - SpikeEvent se; - kernel().event_delivery_manager.send( *this, se, lag ); - } - - // set new input current - S_.y0_ = input[ Buffers_::I0 ]; - - // reset all values in the currently processed input-buffer slot - B_.input_buffer_.reset_values_all_channels( input_buffer_slot ); - - // log state data - B_.logger_.record_data( origin.get_steps() + lag ); - } + for ( long lag = from; lag < to; ++lag ) + { + if ( S_.r_ == 0 ) + { + // neuron not refractory + S_.y3_ = V_.P30_ * ( S_.y0_ + P_.I_e_ ) + V_.P31_ex_ * S_.dI_ex_ + V_.P32_ex_ * S_.I_ex_ + V_.P31_in_ * S_.dI_in_ + + V_.P32_in_ * S_.I_in_ + V_.expm1_tau_m_ * S_.y3_ + S_.y3_; + + // lower bound of membrane potential + S_.y3_ = ( S_.y3_ < P_.LowerBound_ ? P_.LowerBound_ : S_.y3_ ); + } + else + { + // neuron is absolute refractory + --S_.r_; + } + + // alpha shape EPSCs + S_.I_ex_ = V_.P21_ex_ * S_.dI_ex_ + V_.P22_ex_ * S_.I_ex_; + S_.dI_ex_ *= V_.P11_ex_; + + // get read access to the correct input-buffer slot + const size_t input_buffer_slot = kernel().event_delivery_manager.get_modulo( lag ); + auto& input = B_.input_buffer_.get_values_all_channels( input_buffer_slot ); + + // Apply spikes delivered in this step; spikes arriving at T+1 have + // an immediate effect on the state of the neuron + V_.weighted_spikes_ex_ = input[ Buffers_::SYN_EX ]; + S_.dI_ex_ += V_.EPSCInitialValue_ * V_.weighted_spikes_ex_; + + // alpha shape EPSCs + S_.I_in_ = V_.P21_in_ * S_.dI_in_ + V_.P22_in_ * S_.I_in_; + S_.dI_in_ *= V_.P11_in_; + + // Apply spikes delivered in this step; spikes arriving at T+1 have + // an immediate effect on the state of the neuron + V_.weighted_spikes_in_ = input[ Buffers_::SYN_IN ]; + S_.dI_in_ += V_.IPSCInitialValue_ * V_.weighted_spikes_in_; + + // threshold crossing + if ( S_.y3_ >= P_.Theta_ ) + { + S_.r_ = V_.RefractoryCounts_; + S_.y3_ = P_.V_reset_; + // A supra-threshold membrane potential should never be observable. + // The reset at the time of threshold crossing enables accurate + // integration independent of the computation step size, see [2,3] for + // details. + + set_spiketime( Time::step( origin.get_steps() + lag + 1 ) ); + SpikeEvent se; + kernel().event_delivery_manager.send( *this, se, lag ); + } + + // set new input current + S_.y0_ = input[ Buffers_::I0 ]; + + // reset all values in the currently processed input-buffer slot + B_.input_buffer_.reset_values_all_channels( input_buffer_slot ); + + // log state data + B_.logger_.record_data( origin.get_steps() + lag ); + } } void iaf_psc_alpha_hom::handle( SpikeEvent& e ) { - assert( e.get_delay_steps() > 0 ); + assert( e.get_delay_steps() > 0 ); - const size_t input_buffer_slot = kernel().event_delivery_manager.get_modulo( - e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ) ); + const size_t input_buffer_slot = kernel().event_delivery_manager.get_modulo( + e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ) ); - const double s = e.get_weight() * e.get_multiplicity(); + const double s = e.get_weight() * e.get_multiplicity(); - // separate buffer channels for excitatory and inhibitory inputs - B_.input_buffer_.add_value( input_buffer_slot, s > 0 ? Buffers_::SYN_EX : Buffers_::SYN_IN, s ); + // separate buffer channels for excitatory and inhibitory inputs + B_.input_buffer_.add_value( input_buffer_slot, s > 0 ? Buffers_::SYN_EX : Buffers_::SYN_IN, s ); } void iaf_psc_alpha_hom::handle( CurrentEvent& e ) { - assert( e.get_delay_steps() > 0 ); + assert( e.get_delay_steps() > 0 ); - const size_t input_buffer_slot = kernel().event_delivery_manager.get_modulo( - e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ) ); + const size_t input_buffer_slot = kernel().event_delivery_manager.get_modulo( + e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ) ); - const double I = e.get_current(); - const double w = e.get_weight(); + const double I = e.get_current(); + const double w = e.get_weight(); - B_.input_buffer_.add_value( input_buffer_slot, Buffers_::I0, w * I ); + B_.input_buffer_.add_value( input_buffer_slot, Buffers_::I0, w * I ); } void iaf_psc_alpha_hom::handle( DataLoggingRequest& e ) { - B_.logger_.handle( e ); + B_.logger_.handle( e ); } } // namespace diff --git a/models/iaf_psc_alpha_hom.h b/models/iaf_psc_alpha_hom.h index 04ff3e27d0..b2d387ee91 100644 --- a/models/iaf_psc_alpha_hom.h +++ b/models/iaf_psc_alpha_hom.h @@ -1,24 +1,24 @@ /* -* iaf_psc_alpha_hom.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 . -* -*/ + * iaf_psc_alpha_hom.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 IAF_PSC_ALPHA_HOM_H #define IAF_PSC_ALPHA_HOM_H @@ -215,271 +215,271 @@ class iaf_psc_alpha_hom : public ArchivingNodeHom { public: - iaf_psc_alpha_hom(); - iaf_psc_alpha_hom( const iaf_psc_alpha_hom& ); + iaf_psc_alpha_hom(); + iaf_psc_alpha_hom( const iaf_psc_alpha_hom& ); - /** - * Import sets of overloaded virtual functions. - * @see Technical Issues / Virtual Functions: Overriding, Overloading, and - * Hiding - */ - using Node::handle; - using Node::handles_test_event; + /** + * Import sets of overloaded virtual functions. + * @see Technical Issues / Virtual Functions: Overriding, Overloading, and + * Hiding + */ + using Node::handle; + using Node::handles_test_event; - size_t send_test_event( Node&, size_t, synindex, bool ) override; + size_t send_test_event( Node&, size_t, synindex, bool ) override; - void handle( SpikeEvent& ) override; - void handle( CurrentEvent& ) override; - void handle( DataLoggingRequest& ) override; + void handle( SpikeEvent& ) override; + void handle( CurrentEvent& ) override; + void handle( DataLoggingRequest& ) override; - size_t handles_test_event( SpikeEvent&, size_t ) override; - size_t handles_test_event( CurrentEvent&, size_t ) override; - size_t handles_test_event( DataLoggingRequest&, size_t ) override; + size_t handles_test_event( SpikeEvent&, size_t ) override; + size_t handles_test_event( CurrentEvent&, size_t ) override; + size_t handles_test_event( DataLoggingRequest&, size_t ) override; - void get_status( DictionaryDatum& ) const override; - void set_status( const DictionaryDatum& ) override; + void get_status( DictionaryDatum& ) const override; + void set_status( const DictionaryDatum& ) override; private: - void init_buffers_() override; - void pre_run_hook() override; + void init_buffers_() override; + void pre_run_hook() override; - void update( Time const&, const long, const long ) override; + void update( Time const&, const long, const long ) override; - // The next two classes need to be friends to access the State_ class/member - friend class RecordablesMap< iaf_psc_alpha_hom >; - friend class UniversalDataLogger< iaf_psc_alpha_hom >; + // The next two classes need to be friends to access the State_ class/member + friend class RecordablesMap< iaf_psc_alpha_hom >; + friend class UniversalDataLogger< iaf_psc_alpha_hom >; - // ---------------------------------------------------------------- + // ---------------------------------------------------------------- - struct Parameters_ - { - /** Membrane time constant in ms. */ - double Tau_; + struct Parameters_ + { + /** Membrane time constant in ms. */ + double Tau_; - /** Membrane capacitance in pF. */ - double C_; + /** Membrane capacitance in pF. */ + double C_; - /** Refractory period in ms. */ - double TauR_; + /** Refractory period in ms. */ + double TauR_; - /** Resting potential in mV. */ - double E_L_; + /** Resting potential in mV. */ + double E_L_; - /** External current in pA */ - double I_e_; + /** External current in pA */ + double I_e_; - /** Reset value of the membrane potential */ - double V_reset_; + /** Reset value of the membrane potential */ + double V_reset_; - /** Threshold, RELATIVE TO RESTING POTENTIAL(!). - I.e. the real threshold is (E_L_+Theta_). */ - double Theta_; + /** Threshold, RELATIVE TO RESTING POTENTIAL(!). + I.e. the real threshold is (E_L_+Theta_). */ + double Theta_; - /** Lower bound, RELATIVE TO RESTING POTENTIAL(!). - I.e. the real lower bound is (LowerBound_+E_L_). */ - double LowerBound_; + /** Lower bound, RELATIVE TO RESTING POTENTIAL(!). + I.e. the real lower bound is (LowerBound_+E_L_). */ + double LowerBound_; - /** Time constant of excitatory synaptic current in ms. */ - double tau_ex_; + /** Time constant of excitatory synaptic current in ms. */ + double tau_ex_; - /** Time constant of inhibitory synaptic current in ms. */ - double tau_in_; + /** Time constant of inhibitory synaptic current in ms. */ + double tau_in_; - Parameters_(); //!< Sets default parameter values + Parameters_(); //!< Sets default parameter values - void get( DictionaryDatum& ) const; //!< Store current values in dictionary - - /** Set values from dictionary. - * @returns Change in reversal potential E_L, to be passed to State_::set() - */ - double set( const DictionaryDatum&, Node* node ); - }; - - // ---------------------------------------------------------------- - - struct State_ - { - double y0_; //!< Constant current - double dI_ex_; - double I_ex_; - double dI_in_; - double I_in_; - //! This is the membrane potential RELATIVE TO RESTING POTENTIAL. - double y3_; - - int r_; //!< Number of refractory steps remaining - - State_(); //!< Default initialization - - void get( DictionaryDatum&, const Parameters_& ) const; - - /** Set values from dictionary. - * @param dictionary to take data from - * @param current parameters - * @param Change in reversal potential E_L specified by this dict - */ - void set( const DictionaryDatum&, const Parameters_&, double, Node* node ); - }; - - // ---------------------------------------------------------------- - - struct Buffers_ - { - - Buffers_( iaf_psc_alpha_hom& ); - Buffers_( const Buffers_&, iaf_psc_alpha_hom& ); - - //! Indices for access to different channels of input_buffer_ - enum - { - SYN_IN = 0, - SYN_EX, - I0, - NUM_INPUT_CHANNELS - }; - - /** buffers and sums up incoming spikes/currents */ - MultiChannelInputBuffer< NUM_INPUT_CHANNELS > input_buffer_; - - //! Logger for all analog data - UniversalDataLogger< iaf_psc_alpha_hom > logger_; - }; - - // ---------------------------------------------------------------- - - struct Variables_ - { - - /** Amplitude of the synaptic current. - This value is chosen such that a postsynaptic potential with - weight one has an amplitude of 1 mV. - */ - double EPSCInitialValue_; - double IPSCInitialValue_; - int RefractoryCounts_; - - double P11_ex_; - double P21_ex_; - double P22_ex_; - double P31_ex_; - double P32_ex_; - double P11_in_; - double P21_in_; - double P22_in_; - double P31_in_; - double P32_in_; - double P30_; - double P33_; - double expm1_tau_m_; - - double weighted_spikes_ex_; - double weighted_spikes_in_; - }; - - // Access functions for UniversalDataLogger ------------------------------- - - //! Read out the real membrane potential - inline double - get_V_m_() const - { - return S_.y3_ + P_.E_L_; - } - - inline double - get_I_syn_ex_() const - { - return S_.I_ex_; - } - - inline double - get_I_syn_in_() const - { - return S_.I_in_; - } - - // Data members ----------------------------------------------------------- - - /** - * Instances of private data structures for the different types - * of data pertaining to the model. - * @note The order of definitions is important for speed. - * @{ - */ - Parameters_ P_; - State_ S_; - Variables_ V_; - Buffers_ B_; - /** @} */ - - //! Mapping of recordables names to access functions - static RecordablesMap< iaf_psc_alpha_hom > recordablesMap_; + void get( DictionaryDatum& ) const; //!< Store current values in dictionary + + /** Set values from dictionary. + * @returns Change in reversal potential E_L, to be passed to State_::set() + */ + double set( const DictionaryDatum&, Node* node ); + }; + + // ---------------------------------------------------------------- + + struct State_ + { + double y0_; //!< Constant current + double dI_ex_; + double I_ex_; + double dI_in_; + double I_in_; + //! This is the membrane potential RELATIVE TO RESTING POTENTIAL. + double y3_; + + int r_; //!< Number of refractory steps remaining + + State_(); //!< Default initialization + + void get( DictionaryDatum&, const Parameters_& ) const; + + /** Set values from dictionary. + * @param dictionary to take data from + * @param current parameters + * @param Change in reversal potential E_L specified by this dict + */ + void set( const DictionaryDatum&, const Parameters_&, double, Node* node ); + }; + + // ---------------------------------------------------------------- + + struct Buffers_ + { + + Buffers_( iaf_psc_alpha_hom& ); + Buffers_( const Buffers_&, iaf_psc_alpha_hom& ); + + //! Indices for access to different channels of input_buffer_ + enum + { + SYN_IN = 0, + SYN_EX, + I0, + NUM_INPUT_CHANNELS + }; + + /** buffers and sums up incoming spikes/currents */ + MultiChannelInputBuffer< NUM_INPUT_CHANNELS > input_buffer_; + + //! Logger for all analog data + UniversalDataLogger< iaf_psc_alpha_hom > logger_; + }; + + // ---------------------------------------------------------------- + + struct Variables_ + { + + /** Amplitude of the synaptic current. + This value is chosen such that a postsynaptic potential with + weight one has an amplitude of 1 mV. + */ + double EPSCInitialValue_; + double IPSCInitialValue_; + int RefractoryCounts_; + + double P11_ex_; + double P21_ex_; + double P22_ex_; + double P31_ex_; + double P32_ex_; + double P11_in_; + double P21_in_; + double P22_in_; + double P31_in_; + double P32_in_; + double P30_; + double P33_; + double expm1_tau_m_; + + double weighted_spikes_ex_; + double weighted_spikes_in_; + }; + + // Access functions for UniversalDataLogger ------------------------------- + + //! Read out the real membrane potential + inline double + get_V_m_() const + { + return S_.y3_ + P_.E_L_; + } + + inline double + get_I_syn_ex_() const + { + return S_.I_ex_; + } + + inline double + get_I_syn_in_() const + { + return S_.I_in_; + } + + // Data members ----------------------------------------------------------- + + /** + * Instances of private data structures for the different types + * of data pertaining to the model. + * @note The order of definitions is important for speed. + * @{ + */ + Parameters_ P_; + State_ S_; + Variables_ V_; + Buffers_ B_; + /** @} */ + + //! Mapping of recordables names to access functions + static RecordablesMap< iaf_psc_alpha_hom > recordablesMap_; }; inline size_t nest::iaf_psc_alpha_hom::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 ); + SpikeEvent e; + e.set_sender( *this ); + return target.handles_test_event( e, receptor_type ); } inline size_t iaf_psc_alpha_hom::handles_test_event( SpikeEvent&, size_t receptor_type ) { - if ( receptor_type != 0 ) - { - throw UnknownReceptorType( receptor_type, get_name() ); - } - return 0; + if ( receptor_type != 0 ) + { + throw UnknownReceptorType( receptor_type, get_name() ); + } + return 0; } inline size_t iaf_psc_alpha_hom::handles_test_event( CurrentEvent&, size_t receptor_type ) { - if ( receptor_type != 0 ) - { - throw UnknownReceptorType( receptor_type, get_name() ); - } - return 0; + if ( receptor_type != 0 ) + { + throw UnknownReceptorType( receptor_type, get_name() ); + } + return 0; } inline size_t iaf_psc_alpha_hom::handles_test_event( DataLoggingRequest& dlr, size_t receptor_type ) { - if ( receptor_type != 0 ) - { - throw UnknownReceptorType( receptor_type, get_name() ); - } - return B_.logger_.connect_logging_device( dlr, recordablesMap_ ); + if ( receptor_type != 0 ) + { + throw UnknownReceptorType( receptor_type, get_name() ); + } + return B_.logger_.connect_logging_device( dlr, recordablesMap_ ); } inline void iaf_psc_alpha_hom::get_status( DictionaryDatum& d ) const { - P_.get( d ); - S_.get( d, P_ ); - ArchivingNodeHom::get_status( d ); + P_.get( d ); + S_.get( d, P_ ); + ArchivingNodeHom::get_status( d ); - ( *d )[ names::recordables ] = recordablesMap_.get_list(); + ( *d )[ names::recordables ] = recordablesMap_.get_list(); } inline void iaf_psc_alpha_hom::set_status( const DictionaryDatum& d ) { - Parameters_ ptmp = P_; // temporary copy in case of errors - const double delta_EL = ptmp.set( d, this ); // throws if BadProperty - State_ stmp = S_; // temporary copy in case of errors - stmp.set( d, ptmp, delta_EL, this ); // throws if BadProperty - - // We now know that (ptmp, stmp) are consistent. We do not - // write them back to (P_, S_) before we are also sure that - // the properties to be set in the parent class are internally - // consistent. - ArchivingNodeHom::set_status( d ); - - // if we get here, temporaries contain consistent set of properties - P_ = ptmp; - S_ = stmp; + Parameters_ ptmp = P_; // temporary copy in case of errors + const double delta_EL = ptmp.set( d, this ); // throws if BadProperty + State_ stmp = S_; // temporary copy in case of errors + stmp.set( d, ptmp, delta_EL, this ); // throws if BadProperty + + // We now know that (ptmp, stmp) are consistent. We do not + // write them back to (P_, S_) before we are also sure that + // the properties to be set in the parent class are internally + // consistent. + ArchivingNodeHom::set_status( d ); + + // if we get here, temporaries contain consistent set of properties + P_ = ptmp; + S_ = stmp; } } // namespace diff --git a/nestkernel/archiving_node.cpp b/nestkernel/archiving_node.cpp index 03e7105fd3..7ad511eb3c 100644 --- a/nestkernel/archiving_node.cpp +++ b/nestkernel/archiving_node.cpp @@ -1,24 +1,24 @@ /* -* archiving_node.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 . -* -*/ + * archiving_node.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 "archiving_node.h" @@ -34,243 +34,243 @@ namespace nest // member functions for ArchivingNode nest::ArchivingNode::ArchivingNode() - : n_incoming_( 0 ) - , Kminus_( 0.0 ) - , Kminus_triplet_( 0.0 ) - , tau_minus_( 20.0 ) - , tau_minus_inv_( 1. / tau_minus_ ) - , tau_minus_triplet_( 110.0 ) - , tau_minus_triplet_inv_( 1. / tau_minus_triplet_ ) - , max_delay_( 0 ) - , trace_( 0.0 ) - , last_spike_( -1.0 ) + : n_incoming_( 0 ) + , Kminus_( 0.0 ) + , Kminus_triplet_( 0.0 ) + , tau_minus_( 20.0 ) + , tau_minus_inv_( 1. / tau_minus_ ) + , tau_minus_triplet_( 110.0 ) + , tau_minus_triplet_inv_( 1. / tau_minus_triplet_ ) + , max_delay_( 0 ) + , trace_( 0.0 ) + , last_spike_( -1.0 ) { } nest::ArchivingNode::ArchivingNode( const ArchivingNode& n ) - : StructuralPlasticityNode( n ) - , n_incoming_( n.n_incoming_ ) - , Kminus_( n.Kminus_ ) - , Kminus_triplet_( n.Kminus_triplet_ ) - , tau_minus_( n.tau_minus_ ) - , tau_minus_inv_( n.tau_minus_inv_ ) - , tau_minus_triplet_( n.tau_minus_triplet_ ) - , tau_minus_triplet_inv_( n.tau_minus_triplet_inv_ ) - , max_delay_( n.max_delay_ ) - , trace_( n.trace_ ) - , last_spike_( n.last_spike_ ) + : StructuralPlasticityNode( n ) + , n_incoming_( n.n_incoming_ ) + , Kminus_( n.Kminus_ ) + , Kminus_triplet_( n.Kminus_triplet_ ) + , tau_minus_( n.tau_minus_ ) + , tau_minus_inv_( n.tau_minus_inv_ ) + , tau_minus_triplet_( n.tau_minus_triplet_ ) + , tau_minus_triplet_inv_( n.tau_minus_triplet_inv_ ) + , max_delay_( n.max_delay_ ) + , trace_( n.trace_ ) + , last_spike_( n.last_spike_ ) { } void ArchivingNode::register_stdp_connection( double t_first_read, double delay ) { - // Mark all entries in the deque, which we will not read in future as read by - // this input, so that we safely increment the incoming number of - // connections afterwards without leaving spikes in the history. - // For details see bug #218. MH 08-04-22 + // Mark all entries in the deque, which we will not read in future as read by + // this input, so that we safely increment the incoming number of + // connections afterwards without leaving spikes in the history. + // For details see bug #218. MH 08-04-22 - for ( std::deque< histentry >::iterator runner = history_.begin(); - runner != history_.end() and ( t_first_read - runner->t_ > -1.0 * kernel().connection_manager.get_stdp_eps() ); - ++runner ) - { - ( runner->access_counter_ )++; - } + for ( std::deque< histentry >::iterator runner = history_.begin(); + runner != history_.end() and ( t_first_read - runner->t_ > -1.0 * kernel().connection_manager.get_stdp_eps() ); + ++runner ) + { + ( runner->access_counter_ )++; + } - n_incoming_++; + n_incoming_++; - max_delay_ = std::max( delay, max_delay_ ); + max_delay_ = std::max( delay, max_delay_ ); } double nest::ArchivingNode::get_K_value( double t ) { - // case when the neuron has not yet spiked - if ( history_.empty() ) - { - trace_ = 0.; - return trace_; - } + // case when the neuron has not yet spiked + if ( history_.empty() ) + { + trace_ = 0.; + return trace_; + } - // search for the latest post spike in the history buffer that came strictly - // before `t` - int i = history_.size() - 1; - while ( i >= 0 ) - { - if ( t - history_[ i ].t_ > kernel().connection_manager.get_stdp_eps() ) - { - trace_ = ( history_[ i ].Kminus_ * std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ) ); - return trace_; - } - --i; - } + // search for the latest post spike in the history buffer that came strictly + // before `t` + int i = history_.size() - 1; + while ( i >= 0 ) + { + if ( t - history_[ i ].t_ > kernel().connection_manager.get_stdp_eps() ) + { + trace_ = ( history_[ i ].Kminus_ * std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ) ); + return trace_; + } + --i; + } - // this case occurs when the trace was requested at a time precisely at or - // before the first spike in the history - trace_ = 0.; - return trace_; + // this case occurs when the trace was requested at a time precisely at or + // before the first spike in the history + trace_ = 0.; + return trace_; } void nest::ArchivingNode::get_K_values( double t, - double& K_value, - double& nearest_neighbor_K_value, - double& K_triplet_value ) + double& K_value, + double& nearest_neighbor_K_value, + double& K_triplet_value ) { - // case when the neuron has not yet spiked - if ( history_.empty() ) - { - K_triplet_value = Kminus_triplet_; - nearest_neighbor_K_value = Kminus_; - K_value = Kminus_; - return; - } + // case when the neuron has not yet spiked + if ( history_.empty() ) + { + K_triplet_value = Kminus_triplet_; + nearest_neighbor_K_value = Kminus_; + K_value = Kminus_; + return; + } - // search for the latest post spike in the history buffer that came strictly - // before `t` - int i = history_.size() - 1; - while ( i >= 0 ) - { - if ( t - history_[ i ].t_ > kernel().connection_manager.get_stdp_eps() ) - { - K_triplet_value = - ( history_[ i ].Kminus_triplet_ * std::exp( ( history_[ i ].t_ - t ) * tau_minus_triplet_inv_ ) ); - K_value = ( history_[ i ].Kminus_ * std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ) ); - nearest_neighbor_K_value = std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ); - return; - } - --i; - } + // search for the latest post spike in the history buffer that came strictly + // before `t` + int i = history_.size() - 1; + while ( i >= 0 ) + { + if ( t - history_[ i ].t_ > kernel().connection_manager.get_stdp_eps() ) + { + K_triplet_value = + ( history_[ i ].Kminus_triplet_ * std::exp( ( history_[ i ].t_ - t ) * tau_minus_triplet_inv_ ) ); + K_value = ( history_[ i ].Kminus_ * std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ) ); + nearest_neighbor_K_value = std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ); + return; + } + --i; + } - // this case occurs when the trace was requested at a time precisely at or - // before the first spike in the history - K_triplet_value = 0.0; - nearest_neighbor_K_value = 0.0; - K_value = 0.0; + // this case occurs when the trace was requested at a time precisely at or + // before the first spike in the history + K_triplet_value = 0.0; + nearest_neighbor_K_value = 0.0; + K_value = 0.0; } void nest::ArchivingNode::get_history( double t1, - double t2, - std::deque< histentry >::iterator* start, - std::deque< histentry >::iterator* finish ) + double t2, + std::deque< histentry >::iterator* start, + std::deque< histentry >::iterator* finish ) { - *finish = history_.end(); - if ( history_.empty() ) - { - *start = *finish; - return; - } - std::deque< histentry >::reverse_iterator runner = history_.rbegin(); - const double t2_lim = t2 + kernel().connection_manager.get_stdp_eps(); - const double t1_lim = t1 + kernel().connection_manager.get_stdp_eps(); - while ( runner != history_.rend() and runner->t_ >= t2_lim ) - { - ++runner; - } - *finish = runner.base(); - while ( runner != history_.rend() and runner->t_ >= t1_lim ) - { - runner->access_counter_++; - ++runner; - } - *start = runner.base(); + *finish = history_.end(); + if ( history_.empty() ) + { + *start = *finish; + return; + } + std::deque< histentry >::reverse_iterator runner = history_.rbegin(); + const double t2_lim = t2 + kernel().connection_manager.get_stdp_eps(); + const double t1_lim = t1 + kernel().connection_manager.get_stdp_eps(); + while ( runner != history_.rend() and runner->t_ >= t2_lim ) + { + ++runner; + } + *finish = runner.base(); + while ( runner != history_.rend() and runner->t_ >= t1_lim ) + { + runner->access_counter_++; + ++runner; + } + *start = runner.base(); } void nest::ArchivingNode::set_spiketime( Time const& t_sp, double offset ) { - StructuralPlasticityNode::set_spiketime( t_sp, offset ); + StructuralPlasticityNode::set_spiketime( t_sp, offset ); - const double t_sp_ms = t_sp.get_ms() - offset; + const double t_sp_ms = t_sp.get_ms() - offset; - if ( n_incoming_ ) - { - // prune all spikes from history which are no longer needed - // only remove a spike if: - // - its access counter indicates it has been read out by all connected - // STDP synapses, and - // - there is another, later spike, that is strictly more than - // (min_global_delay + max_local_delay + eps) away from the new spike (at t_sp_ms) - while ( history_.size() > 1 ) - { - const double next_t_sp = history_[ 1 ].t_; - if ( history_.front().access_counter_ >= n_incoming_ - and t_sp_ms - next_t_sp > max_delay_ + Time::delay_steps_to_ms( kernel().connection_manager.get_min_delay() ) - + kernel().connection_manager.get_stdp_eps() ) - { - history_.pop_front(); - } - else - { - break; - } - } - // update spiking history - Kminus_ = Kminus_ * std::exp( ( last_spike_ - t_sp_ms ) * tau_minus_inv_ ) + 1.0; - Kminus_triplet_ = Kminus_triplet_ * std::exp( ( last_spike_ - t_sp_ms ) * tau_minus_triplet_inv_ ) + 1.0; - last_spike_ = t_sp_ms; - history_.push_back( histentry( last_spike_, Kminus_, Kminus_triplet_, 0 ) ); - } - else - { - last_spike_ = t_sp_ms; - } + if ( n_incoming_ ) + { + // prune all spikes from history which are no longer needed + // only remove a spike if: + // - its access counter indicates it has been read out by all connected + // STDP synapses, and + // - there is another, later spike, that is strictly more than + // (min_global_delay + max_local_delay + eps) away from the new spike (at t_sp_ms) + while ( history_.size() > 1 ) + { + const double next_t_sp = history_[ 1 ].t_; + if ( history_.front().access_counter_ >= n_incoming_ + and t_sp_ms - next_t_sp > max_delay_ + Time::delay_steps_to_ms( kernel().connection_manager.get_min_delay() ) + + kernel().connection_manager.get_stdp_eps() ) + { + history_.pop_front(); + } + else + { + break; + } + } + // update spiking history + Kminus_ = Kminus_ * std::exp( ( last_spike_ - t_sp_ms ) * tau_minus_inv_ ) + 1.0; + Kminus_triplet_ = Kminus_triplet_ * std::exp( ( last_spike_ - t_sp_ms ) * tau_minus_triplet_inv_ ) + 1.0; + last_spike_ = t_sp_ms; + history_.push_back( histentry( last_spike_, Kminus_, Kminus_triplet_, 0 ) ); + } + else + { + last_spike_ = t_sp_ms; + } } void nest::ArchivingNode::get_status( DictionaryDatum& d ) const { - def< double >( d, names::t_spike, get_spiketime_ms() ); - def< double >( d, names::tau_minus, tau_minus_ ); - def< double >( d, names::tau_minus_triplet, tau_minus_triplet_ ); - def< double >( d, names::post_trace, trace_ ); + def< double >( d, names::t_spike, get_spiketime_ms() ); + def< double >( d, names::tau_minus, tau_minus_ ); + def< double >( d, names::tau_minus_triplet, tau_minus_triplet_ ); + def< double >( d, names::post_trace, trace_ ); #ifdef DEBUG_ARCHIVER - def< int >( d, names::archiver_length, history_.size() ); + def< int >( d, names::archiver_length, history_.size() ); #endif - // add status dict items from the parent class - StructuralPlasticityNode::get_status( d ); + // add status dict items from the parent class + StructuralPlasticityNode::get_status( d ); } void nest::ArchivingNode::set_status( const DictionaryDatum& d ) { - // We need to preserve values in case invalid values are set - double new_tau_minus = tau_minus_; - double new_tau_minus_triplet = tau_minus_triplet_; - updateValue< double >( d, names::tau_minus, new_tau_minus ); - updateValue< double >( d, names::tau_minus_triplet, new_tau_minus_triplet ); + // We need to preserve values in case invalid values are set + double new_tau_minus = tau_minus_; + double new_tau_minus_triplet = tau_minus_triplet_; + updateValue< double >( d, names::tau_minus, new_tau_minus ); + updateValue< double >( d, names::tau_minus_triplet, new_tau_minus_triplet ); - if ( new_tau_minus <= 0.0 or new_tau_minus_triplet <= 0.0 ) - { - throw BadProperty( "All time constants must be strictly positive." ); - } + if ( new_tau_minus <= 0.0 or new_tau_minus_triplet <= 0.0 ) + { + throw BadProperty( "All time constants must be strictly positive." ); + } - StructuralPlasticityNode::set_status( d ); + StructuralPlasticityNode::set_status( d ); - // do the actual update - tau_minus_ = new_tau_minus; - tau_minus_triplet_ = new_tau_minus_triplet; - tau_minus_inv_ = 1. / tau_minus_; - tau_minus_triplet_inv_ = 1. / tau_minus_triplet_; + // do the actual update + tau_minus_ = new_tau_minus; + tau_minus_triplet_ = new_tau_minus_triplet; + tau_minus_inv_ = 1. / tau_minus_; + tau_minus_triplet_inv_ = 1. / tau_minus_triplet_; - // check, if to clear spike history and K_minus - bool clear = false; - updateValue< bool >( d, names::clear, clear ); - if ( clear ) - { - clear_history(); - } + // check, if to clear spike history and K_minus + bool clear = false; + updateValue< bool >( d, names::clear, clear ); + if ( clear ) + { + clear_history(); + } } void nest::ArchivingNode::clear_history() { - last_spike_ = -1.0; - Kminus_ = 0.0; - Kminus_triplet_ = 0.0; - history_.clear(); + last_spike_ = -1.0; + Kminus_ = 0.0; + Kminus_triplet_ = 0.0; + history_.clear(); } diff --git a/nestkernel/archiving_node.h b/nestkernel/archiving_node.h index 825e77553a..86b402ed1a 100644 --- a/nestkernel/archiving_node.h +++ b/nestkernel/archiving_node.h @@ -1,24 +1,24 @@ /* -* archiving_node.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 . -* -*/ + * archiving_node.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 ARCHIVING_NODE_H #define ARCHIVING_NODE_H @@ -43,111 +43,111 @@ namespace nest { /** -* A node which archives spike history for the purposes of spike-timing -* dependent plasticity (STDP) -*/ + * A node which archives spike history for the purposes of spike-timing + * dependent plasticity (STDP) + */ class ArchivingNode : public StructuralPlasticityNode { public: - ArchivingNode(); - - ArchivingNode( const ArchivingNode& ); - - /** - * Return the Kminus (synaptic trace) value at time t (given in ms). - * - * When the trace is requested at the exact same time that the neuron emits a spike, - * the trace value as it was just before the spike is returned. - */ - double get_K_value( double t ) override; - - /** - * Write the different STDP K values at time t (in ms) to the provided locations. - * - * @param Kminus the eligibility trace for STDP - * @param nearest_neighbour_Kminus eligibility trace for nearest-neighbour STDP, like Kminus, - but increased to 1, rather than by 1, when a spike occurs - * @param Kminus_triplet eligibility trace for triplet STDP - * - * @throws UnexpectedEvent - */ - void get_K_values( double t, double& Kminus, double& nearest_neighbor_Kminus, double& Kminus_triplet ) override; - - /** - * Return the triplet Kminus value for the associated iterator. - */ - double get_K_triplet_value( const std::deque< histentry >::iterator& iter ); - - /** - * Return the spike times (in steps) of spikes which occurred in the range [t1,t2]. - */ - void get_history( double t1, - double t2, - std::deque< histentry >::iterator* start, - std::deque< histentry >::iterator* finish ) override; - - /** - * Register a new incoming STDP connection. - * - * t_first_read: The newly registered synapse will read the history entries - * with t > t_first_read. - */ - void register_stdp_connection( double t_first_read, double delay ) override; - - void get_status( DictionaryDatum& d ) const override; - void set_status( const DictionaryDatum& d ) override; + ArchivingNode(); + + ArchivingNode( const ArchivingNode& ); + + /** + * Return the Kminus (synaptic trace) value at time t (given in ms). + * + * When the trace is requested at the exact same time that the neuron emits a spike, + * the trace value as it was just before the spike is returned. + */ + double get_K_value( double t ) override; + + /** + * Write the different STDP K values at time t (in ms) to the provided locations. + * + * @param Kminus the eligibility trace for STDP + * @param nearest_neighbour_Kminus eligibility trace for nearest-neighbour STDP, like Kminus, + but increased to 1, rather than by 1, when a spike occurs + * @param Kminus_triplet eligibility trace for triplet STDP + * + * @throws UnexpectedEvent + */ + void get_K_values( double t, double& Kminus, double& nearest_neighbor_Kminus, double& Kminus_triplet ) override; + + /** + * Return the triplet Kminus value for the associated iterator. + */ + double get_K_triplet_value( const std::deque< histentry >::iterator& iter ); + + /** + * Return the spike times (in steps) of spikes which occurred in the range [t1,t2]. + */ + void get_history( double t1, + double t2, + std::deque< histentry >::iterator* start, + std::deque< histentry >::iterator* finish ) override; + + /** + * Register a new incoming STDP connection. + * + * t_first_read: The newly registered synapse will read the history entries + * with t > t_first_read. + */ + void register_stdp_connection( double t_first_read, double delay ) override; + + void get_status( DictionaryDatum& d ) const override; + void set_status( const DictionaryDatum& d ) override; protected: - /** - * Record spike history - */ - void set_spiketime( Time const& t_sp, double offset = 0.0 ); - - /** - * Return most recent spike time in ms - */ - inline double get_spiketime_ms() const; - - /** - * Clear spike history - */ - void clear_history(); - - /** - * Number of incoming connections from STDP connectors. - * - * This variable is needed to determine if every incoming connection has - * read the spikehistory for a given point in time - */ - size_t n_incoming_; + /** + * Record spike history + */ + void set_spiketime( Time const& t_sp, double offset = 0.0 ); + + /** + * Return most recent spike time in ms + */ + inline double get_spiketime_ms() const; + + /** + * Clear spike history + */ + void clear_history(); + + /** + * Number of incoming connections from STDP connectors. + * + * This variable is needed to determine if every incoming connection has + * read the spikehistory for a given point in time + */ + size_t n_incoming_; private: - // sum exp(-(t-ti)/tau_minus) - double Kminus_; + // sum exp(-(t-ti)/tau_minus) + double Kminus_; - // sum exp(-(t-ti)/tau_minus_triplet) - double Kminus_triplet_; + // sum exp(-(t-ti)/tau_minus_triplet) + double Kminus_triplet_; - double tau_minus_; - double tau_minus_inv_; + double tau_minus_; + double tau_minus_inv_; - // time constant for triplet low pass filtering of "post" spike train - double tau_minus_triplet_; - double tau_minus_triplet_inv_; + // time constant for triplet low pass filtering of "post" spike train + double tau_minus_triplet_; + double tau_minus_triplet_inv_; - double max_delay_; - double trace_; + double max_delay_; + double trace_; - double last_spike_; + double last_spike_; - // spiking history needed by stdp synapses - std::deque< histentry > history_; + // spiking history needed by stdp synapses + std::deque< histentry > history_; }; inline double ArchivingNode::get_spiketime_ms() const { - return last_spike_; + return last_spike_; } } // of namespace diff --git a/nestkernel/archiving_node_hom.cpp b/nestkernel/archiving_node_hom.cpp index e73a117c7f..cc653d4ec3 100644 --- a/nestkernel/archiving_node_hom.cpp +++ b/nestkernel/archiving_node_hom.cpp @@ -219,7 +219,8 @@ nest::ArchivingNodeHom::set_spiketime( Time const& t_sp, double offset ) // update spiking history Kminus_ = Kminus_ * std::exp( Time( Time::step( last_spike_ - t_sp_steps ) ).get_ms() * tau_minus_inv_ ) + 1.0; Kminus_triplet_ = - Kminus_triplet_ * std::exp( Time( Time::step( last_spike_ - t_sp_steps ) ).get_ms() * tau_minus_triplet_inv_ ) + 1.0; + Kminus_triplet_ * std::exp( Time( Time::step( last_spike_ - t_sp_steps ) ).get_ms() * tau_minus_triplet_inv_ ) + + 1.0; last_spike_ = t_sp_steps; history_.push_back( histentry_step( last_spike_, Kminus_, Kminus_triplet_, 0 ) ); } diff --git a/nestkernel/node.h b/nestkernel/node.h index 157b9cfca9..535aa825d0 100644 --- a/nestkernel/node.h +++ b/nestkernel/node.h @@ -765,7 +765,7 @@ class Node */ virtual double get_K_value( long t, size_t& dt_steps ); - virtual double get_K_value (double ); + virtual double get_K_value( double ); virtual double get_LTD_value( double t ); diff --git a/testsuite/pytests/test_stdp_pl_synapse_hom.py b/testsuite/pytests/test_stdp_pl_synapse_hom.py index fba655e902..56f6ac1203 100644 --- a/testsuite/pytests/test_stdp_pl_synapse_hom.py +++ b/testsuite/pytests/test_stdp_pl_synapse_hom.py @@ -59,10 +59,13 @@ def init_params(self): "alpha": 1.0, "mu": 0.4, "tau_plus": self.tau_pre, - "tau_minus": self.tau_post + "tau_minus": self.tau_post, + } + self.synapse_parameters = { + "synapse_model": self.synapse_model, + "receptor_type": 0, + "weight": self.init_weight, } - self.synapse_parameters = {"synapse_model": self.synapse_model, "receptor_type": 0, - "weight": self.init_weight, } self.neuron_parameters = {} self.nest_neuron_model = "iaf_psc_alpha_hom" @@ -224,7 +227,7 @@ def facilitate(w, Kplus): def depress(w, Kminus): new_weight = ( - w - self.synapse_common_properties["alpha"] * self.synapse_common_properties["lambda"] * w * Kminus + w - self.synapse_common_properties["alpha"] * self.synapse_common_properties["lambda"] * w * Kminus ) return new_weight if new_weight > 0.0 else 0.0 @@ -329,16 +332,16 @@ def depress(w, Kminus): return np.array(t_log), np.array(w_log), Kplus_log, Kminus_log def plot_weight_evolution( - self, - pre_spikes, - post_spikes, - t_log, - w_log, - Kpre_log=None, - Kpost_log=None, - pre_indices=slice(-1), - fname_snip="", - title_snip="", + self, + pre_spikes, + post_spikes, + t_log, + w_log, + Kpre_log=None, + Kpost_log=None, + pre_indices=slice(-1), + fname_snip="", + title_snip="", ): if not DEBUG_PLOTS: # make pylint happy if no matplotlib return