From 5023910d36531be3914415f7680625d7cd052d22 Mon Sep 17 00:00:00 2001 From: Mahesh Madhav Date: Tue, 3 Sep 2024 21:24:35 +0000 Subject: [PATCH 1/7] Optimize computation in iaf_psc_alpha_ps model Replace FP divides with FP multiplies and pre-compute loop invariants. Similar to https://github.com/nest/nest-simulator/pull/3304 --- models/iaf_psc_alpha_ps.cpp | 37 +++++++++++++++++++++++++------------ models/iaf_psc_alpha_ps.h | 7 +++++++ 2 files changed, 32 insertions(+), 12 deletions(-) diff --git a/models/iaf_psc_alpha_ps.cpp b/models/iaf_psc_alpha_ps.cpp index acc3e06d5a..24b544b8c0 100644 --- a/models/iaf_psc_alpha_ps.cpp +++ b/models/iaf_psc_alpha_ps.cpp @@ -85,6 +85,7 @@ nest::iaf_psc_alpha_ps::Parameters_::Parameters_() , U_min_( -std::numeric_limits< double >::infinity() ) // mV , U_reset_( -70.0 - E_L_ ) // mV, rel to E_L_ { + compute_local_(); } nest::iaf_psc_alpha_ps::State_::State_() @@ -187,9 +188,21 @@ nest::iaf_psc_alpha_ps::Parameters_::set( const DictionaryDatum& d, Node* node ) throw BadProperty( "All time constants must be strictly positive." ); } + compute_local_(); + return delta_EL; } +void +nest::iaf_psc_alpha_ps::Parameters_::compute_local_() +{ + // pre-compute inverse member variables for speed + inv_c_m_ = 1.0 / c_m_; + inv_tau_m_ = 1.0 / tau_m_; + inv_tau_syn_ex_ = 1.0 / tau_syn_ex_; + inv_tau_syn_in_ = 1.0 / tau_syn_in_; +} + void nest::iaf_psc_alpha_ps::State_::get( DictionaryDatum& d, const Parameters_& p ) const { @@ -268,15 +281,15 @@ nest::iaf_psc_alpha_ps::pre_run_hook() V_.h_ms_ = Time::get_resolution().get_ms(); - V_.psc_norm_ex_ = 1.0 * numerics::e / P_.tau_syn_ex_; - V_.psc_norm_in_ = 1.0 * numerics::e / P_.tau_syn_in_; + V_.psc_norm_ex_ = 1.0 * numerics::e * P_.inv_tau_syn_ex_; + V_.psc_norm_in_ = 1.0 * numerics::e * P_.inv_tau_syn_in_; // pre-compute matrix for full time step - V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ / P_.tau_m_ ); - V_.exp_tau_syn_ex_ = std::exp( -V_.h_ms_ / P_.tau_syn_ex_ ); - V_.exp_tau_syn_in_ = std::exp( -V_.h_ms_ / P_.tau_syn_in_ ); + V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ * P_.inv_tau_m_ ); + V_.exp_tau_syn_ex_ = std::exp( -V_.h_ms_ * P_.inv_tau_syn_ex_ ); + V_.exp_tau_syn_in_ = std::exp( -V_.h_ms_ * P_.inv_tau_syn_in_ ); - V_.P30_ = -P_.tau_m_ / P_.c_m_ * V_.expm1_tau_m_; + V_.P30_ = -P_.tau_m_ * P_.inv_c_m_ * V_.expm1_tau_m_; // these are determined according to a numeric stability criterion propagator_ex_ = IAFPropagatorAlpha( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_ ); std::tie( V_.P31_ex_, V_.P32_ex_ ) = propagator_ex_.evaluate( V_.h_ms_ ); @@ -500,9 +513,9 @@ nest::iaf_psc_alpha_ps::propagate_( const double dt ) // V_m_ remains unchanged at 0.0 while neuron is refractory if ( not S_.is_refractory_ ) { - const double expm1_tau_m = numerics::expm1( -dt / P_.tau_m_ ); + const double expm1_tau_m = numerics::expm1( -dt * P_.inv_tau_m_ ); - const double ps_P30 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; + const double ps_P30 = -P_.tau_m_ * P_.inv_c_m_ * expm1_tau_m; double ps_P31_ex; double ps_P32_ex; @@ -518,8 +531,8 @@ nest::iaf_psc_alpha_ps::propagate_( const double dt ) S_.V_m_ = ( S_.V_m_ < P_.U_min_ ? P_.U_min_ : S_.V_m_ ); } - const double ps_e_TauSyn_ex = std::exp( -dt / P_.tau_syn_ex_ ); - const double ps_e_TauSyn_in = std::exp( -dt / P_.tau_syn_in_ ); + const double ps_e_TauSyn_ex = std::exp( -dt * P_.inv_tau_syn_ex_ ); + const double ps_e_TauSyn_in = std::exp( -dt * P_.inv_tau_syn_in_ ); // now the synaptic components S_.I_ex_ = ps_e_TauSyn_ex * dt * S_.dI_ex_ + ps_e_TauSyn_ex * S_.I_ex_; @@ -578,9 +591,9 @@ nest::iaf_psc_alpha_ps::emit_instant_spike_( Time const& origin, const long lag, double nest::iaf_psc_alpha_ps::threshold_distance( double t_step ) const { - const double expm1_tau_m = numerics::expm1( -t_step / P_.tau_m_ ); + const double expm1_tau_m = numerics::expm1( -t_step * P_.inv_tau_m_ ); - const double ps_P30 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; + const double ps_P30 = -P_.tau_m_ * P_.inv_c_m_ * expm1_tau_m; double ps_P31_ex; double ps_P32_ex; diff --git a/models/iaf_psc_alpha_ps.h b/models/iaf_psc_alpha_ps.h index 9c358ebeaa..4ec9ecba17 100644 --- a/models/iaf_psc_alpha_ps.h +++ b/models/iaf_psc_alpha_ps.h @@ -291,13 +291,17 @@ class iaf_psc_alpha_ps : public ArchivingNode /** Membrane time constant in ms. */ double tau_m_; + double inv_tau_m_; /* inverse of tau_m_ */ /** Time constant of synaptic current in ms. */ double tau_syn_ex_; double tau_syn_in_; + double inv_tau_syn_ex_; /* inverse of tau_syn_ex_ */ + double inv_tau_syn_in_; /* inverse of tau_syn_in_ */ /** Membrane capacitance in pF. */ double c_m_; + double inv_c_m_; /* inverse of c_m_ */ /** Refractory period in ms. */ double t_ref_; @@ -330,6 +334,9 @@ class iaf_psc_alpha_ps : public ArchivingNode * @returns Change in reversal potential E_L, to be passed to State_::set() */ double set( const DictionaryDatum&, Node* ); + + /** Used for pre-computing member variables for performance */ + void compute_local_(); }; // ---------------------------------------------------------------- From de67e0abbf2d609ef1b75642cdf53dbfe316ae62 Mon Sep 17 00:00:00 2001 From: Mahesh Madhav Date: Tue, 3 Sep 2024 21:36:42 +0000 Subject: [PATCH 2/7] Fix formatting --- models/iaf_psc_alpha_ps.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/iaf_psc_alpha_ps.cpp b/models/iaf_psc_alpha_ps.cpp index 24b544b8c0..bc875cce88 100644 --- a/models/iaf_psc_alpha_ps.cpp +++ b/models/iaf_psc_alpha_ps.cpp @@ -197,8 +197,8 @@ void nest::iaf_psc_alpha_ps::Parameters_::compute_local_() { // pre-compute inverse member variables for speed - inv_c_m_ = 1.0 / c_m_; - inv_tau_m_ = 1.0 / tau_m_; + inv_c_m_ = 1.0 / c_m_; + inv_tau_m_ = 1.0 / tau_m_; inv_tau_syn_ex_ = 1.0 / tau_syn_ex_; inv_tau_syn_in_ = 1.0 / tau_syn_in_; } From b7c5e5b0f13b1cc478658ac22664d5d8e91a8c8e Mon Sep 17 00:00:00 2001 From: Mahesh Madhav Date: Fri, 6 Sep 2024 05:22:27 +0000 Subject: [PATCH 3/7] Transfer code from P_ (parameters) to V_ (variables) --- models/iaf_psc_alpha_ps.cpp | 44 ++++++++++++++++--------------------- models/iaf_psc_alpha_ps.h | 11 ++++------ 2 files changed, 23 insertions(+), 32 deletions(-) diff --git a/models/iaf_psc_alpha_ps.cpp b/models/iaf_psc_alpha_ps.cpp index bc875cce88..760383eafe 100644 --- a/models/iaf_psc_alpha_ps.cpp +++ b/models/iaf_psc_alpha_ps.cpp @@ -85,7 +85,6 @@ nest::iaf_psc_alpha_ps::Parameters_::Parameters_() , U_min_( -std::numeric_limits< double >::infinity() ) // mV , U_reset_( -70.0 - E_L_ ) // mV, rel to E_L_ { - compute_local_(); } nest::iaf_psc_alpha_ps::State_::State_() @@ -188,21 +187,9 @@ nest::iaf_psc_alpha_ps::Parameters_::set( const DictionaryDatum& d, Node* node ) throw BadProperty( "All time constants must be strictly positive." ); } - compute_local_(); - return delta_EL; } -void -nest::iaf_psc_alpha_ps::Parameters_::compute_local_() -{ - // pre-compute inverse member variables for speed - inv_c_m_ = 1.0 / c_m_; - inv_tau_m_ = 1.0 / tau_m_; - inv_tau_syn_ex_ = 1.0 / tau_syn_ex_; - inv_tau_syn_in_ = 1.0 / tau_syn_in_; -} - void nest::iaf_psc_alpha_ps::State_::get( DictionaryDatum& d, const Parameters_& p ) const { @@ -281,15 +268,21 @@ nest::iaf_psc_alpha_ps::pre_run_hook() V_.h_ms_ = Time::get_resolution().get_ms(); - V_.psc_norm_ex_ = 1.0 * numerics::e * P_.inv_tau_syn_ex_; - V_.psc_norm_in_ = 1.0 * numerics::e * P_.inv_tau_syn_in_; + // pre-compute inverse member variables + V_.inv_c_m_ = 1.0 / P_.c_m_; + V_.inv_tau_m_ = 1.0 / P_.tau_m_; + V_.inv_tau_syn_ex_ = 1.0 / P_.tau_syn_ex_; + V_.inv_tau_syn_in_ = 1.0 / P_.tau_syn_in_; + + V_.psc_norm_ex_ = 1.0 * numerics::e * V_.inv_tau_syn_ex_; + V_.psc_norm_in_ = 1.0 * numerics::e * V_.inv_tau_syn_in_; // pre-compute matrix for full time step - V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ * P_.inv_tau_m_ ); - V_.exp_tau_syn_ex_ = std::exp( -V_.h_ms_ * P_.inv_tau_syn_ex_ ); - V_.exp_tau_syn_in_ = std::exp( -V_.h_ms_ * P_.inv_tau_syn_in_ ); + V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ * V_.inv_tau_m_ ); + V_.exp_tau_syn_ex_ = std::exp( -V_.h_ms_ * V_.inv_tau_syn_ex_ ); + V_.exp_tau_syn_in_ = std::exp( -V_.h_ms_ * V_.inv_tau_syn_in_ ); - V_.P30_ = -P_.tau_m_ * P_.inv_c_m_ * V_.expm1_tau_m_; + V_.P30_ = -P_.tau_m_ * V_.inv_c_m_ * V_.expm1_tau_m_; // these are determined according to a numeric stability criterion propagator_ex_ = IAFPropagatorAlpha( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_ ); std::tie( V_.P31_ex_, V_.P32_ex_ ) = propagator_ex_.evaluate( V_.h_ms_ ); @@ -303,6 +296,7 @@ nest::iaf_psc_alpha_ps::pre_run_hook() V_.refractory_steps_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); // since t_ref_ >= sim step size, this can only fail in error assert( V_.refractory_steps_ >= 1 ); + } /* ---------------------------------------------------------------- @@ -513,9 +507,9 @@ nest::iaf_psc_alpha_ps::propagate_( const double dt ) // V_m_ remains unchanged at 0.0 while neuron is refractory if ( not S_.is_refractory_ ) { - const double expm1_tau_m = numerics::expm1( -dt * P_.inv_tau_m_ ); + const double expm1_tau_m = numerics::expm1( -dt * V_.inv_tau_m_ ); - const double ps_P30 = -P_.tau_m_ * P_.inv_c_m_ * expm1_tau_m; + const double ps_P30 = -P_.tau_m_ * V_.inv_c_m_ * expm1_tau_m; double ps_P31_ex; double ps_P32_ex; @@ -531,8 +525,8 @@ nest::iaf_psc_alpha_ps::propagate_( const double dt ) S_.V_m_ = ( S_.V_m_ < P_.U_min_ ? P_.U_min_ : S_.V_m_ ); } - const double ps_e_TauSyn_ex = std::exp( -dt * P_.inv_tau_syn_ex_ ); - const double ps_e_TauSyn_in = std::exp( -dt * P_.inv_tau_syn_in_ ); + const double ps_e_TauSyn_ex = std::exp( -dt * V_.inv_tau_syn_ex_ ); + const double ps_e_TauSyn_in = std::exp( -dt * V_.inv_tau_syn_in_ ); // now the synaptic components S_.I_ex_ = ps_e_TauSyn_ex * dt * S_.dI_ex_ + ps_e_TauSyn_ex * S_.I_ex_; @@ -591,9 +585,9 @@ nest::iaf_psc_alpha_ps::emit_instant_spike_( Time const& origin, const long lag, double nest::iaf_psc_alpha_ps::threshold_distance( double t_step ) const { - const double expm1_tau_m = numerics::expm1( -t_step * P_.inv_tau_m_ ); + const double expm1_tau_m = numerics::expm1( -t_step * V_.inv_tau_m_ ); - const double ps_P30 = -P_.tau_m_ * P_.inv_c_m_ * expm1_tau_m; + const double ps_P30 = -P_.tau_m_ * V_.inv_c_m_ * expm1_tau_m; double ps_P31_ex; double ps_P32_ex; diff --git a/models/iaf_psc_alpha_ps.h b/models/iaf_psc_alpha_ps.h index 4ec9ecba17..0f813d7c32 100644 --- a/models/iaf_psc_alpha_ps.h +++ b/models/iaf_psc_alpha_ps.h @@ -291,17 +291,13 @@ class iaf_psc_alpha_ps : public ArchivingNode /** Membrane time constant in ms. */ double tau_m_; - double inv_tau_m_; /* inverse of tau_m_ */ /** Time constant of synaptic current in ms. */ double tau_syn_ex_; double tau_syn_in_; - double inv_tau_syn_ex_; /* inverse of tau_syn_ex_ */ - double inv_tau_syn_in_; /* inverse of tau_syn_in_ */ /** Membrane capacitance in pF. */ double c_m_; - double inv_c_m_; /* inverse of c_m_ */ /** Refractory period in ms. */ double t_ref_; @@ -334,9 +330,6 @@ class iaf_psc_alpha_ps : public ArchivingNode * @returns Change in reversal potential E_L, to be passed to State_::set() */ double set( const DictionaryDatum&, Node* ); - - /** Used for pre-computing member variables for performance */ - void compute_local_(); }; // ---------------------------------------------------------------- @@ -414,6 +407,10 @@ class iaf_psc_alpha_ps : public ArchivingNode double dI_ex_before_; //!< at beginning of mini-step double dI_in_before_; //!< at beginning of mini-step double V_m_before_; //!< at beginning of mini-step + double inv_tau_m_; //!< 1 / tau_m + double inv_tau_syn_ex_; //!< 1 / tau_syn_ex + double inv_tau_syn_in_; //!< 1 / tau_syn_in + double inv_c_m_; //!< 1 / c_m }; // Access functions for UniversalDataLogger ------------------------------- From 186e72ce52099d48af83805a2bd107063b51dd6b Mon Sep 17 00:00:00 2001 From: Mahesh Madhav <67384846+heshpdx@users.noreply.github.com> Date: Thu, 5 Sep 2024 22:40:36 -0700 Subject: [PATCH 4/7] Fix formatting --- models/iaf_psc_alpha_ps.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/models/iaf_psc_alpha_ps.cpp b/models/iaf_psc_alpha_ps.cpp index 760383eafe..467b25bcd9 100644 --- a/models/iaf_psc_alpha_ps.cpp +++ b/models/iaf_psc_alpha_ps.cpp @@ -296,7 +296,6 @@ nest::iaf_psc_alpha_ps::pre_run_hook() V_.refractory_steps_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); // since t_ref_ >= sim step size, this can only fail in error assert( V_.refractory_steps_ >= 1 ); - } /* ---------------------------------------------------------------- From 9701ed29dbdc34105400963591ab7c0e7bdf2716 Mon Sep 17 00:00:00 2001 From: Mahesh Madhav <67384846+heshpdx@users.noreply.github.com> Date: Tue, 10 Sep 2024 19:37:38 -0700 Subject: [PATCH 5/7] Update models/iaf_psc_alpha_ps.cpp --- models/iaf_psc_alpha_ps.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/iaf_psc_alpha_ps.cpp b/models/iaf_psc_alpha_ps.cpp index 467b25bcd9..6768e5d8f3 100644 --- a/models/iaf_psc_alpha_ps.cpp +++ b/models/iaf_psc_alpha_ps.cpp @@ -274,8 +274,8 @@ nest::iaf_psc_alpha_ps::pre_run_hook() V_.inv_tau_syn_ex_ = 1.0 / P_.tau_syn_ex_; V_.inv_tau_syn_in_ = 1.0 / P_.tau_syn_in_; - V_.psc_norm_ex_ = 1.0 * numerics::e * V_.inv_tau_syn_ex_; - V_.psc_norm_in_ = 1.0 * numerics::e * V_.inv_tau_syn_in_; + V_.psc_norm_ex_ = numerics::e * V_.inv_tau_syn_ex_; + V_.psc_norm_in_ = numerics::e * V_.inv_tau_syn_in_; // pre-compute matrix for full time step V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ * V_.inv_tau_m_ ); From c00ebaccbf2d7f1410885f60258beab84756b36e Mon Sep 17 00:00:00 2001 From: Mahesh Madhav Date: Thu, 12 Sep 2024 17:17:08 +0000 Subject: [PATCH 6/7] Revert pre-compute matrix to original computation --- models/iaf_psc_alpha_ps.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/models/iaf_psc_alpha_ps.cpp b/models/iaf_psc_alpha_ps.cpp index 6768e5d8f3..34f57bd4c4 100644 --- a/models/iaf_psc_alpha_ps.cpp +++ b/models/iaf_psc_alpha_ps.cpp @@ -278,9 +278,9 @@ nest::iaf_psc_alpha_ps::pre_run_hook() V_.psc_norm_in_ = numerics::e * V_.inv_tau_syn_in_; // pre-compute matrix for full time step - V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ * V_.inv_tau_m_ ); - V_.exp_tau_syn_ex_ = std::exp( -V_.h_ms_ * V_.inv_tau_syn_ex_ ); - V_.exp_tau_syn_in_ = std::exp( -V_.h_ms_ * V_.inv_tau_syn_in_ ); + V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ / P_.tau_m_ ); + V_.exp_tau_syn_ex_ = std::exp( -V_.h_ms_ / P_.tau_syn_ex_ ); + V_.exp_tau_syn_in_ = std::exp( -V_.h_ms_ / P_.tau_syn_in_ ); V_.P30_ = -P_.tau_m_ * V_.inv_c_m_ * V_.expm1_tau_m_; // these are determined according to a numeric stability criterion From e63b1e9733ececfa125745ccd3bef4a113f960fa Mon Sep 17 00:00:00 2001 From: Mahesh Madhav Date: Thu, 12 Sep 2024 17:18:58 +0000 Subject: [PATCH 7/7] One more rollback --- models/iaf_psc_alpha_ps.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/iaf_psc_alpha_ps.cpp b/models/iaf_psc_alpha_ps.cpp index 34f57bd4c4..885d1f2051 100644 --- a/models/iaf_psc_alpha_ps.cpp +++ b/models/iaf_psc_alpha_ps.cpp @@ -282,7 +282,7 @@ nest::iaf_psc_alpha_ps::pre_run_hook() V_.exp_tau_syn_ex_ = std::exp( -V_.h_ms_ / P_.tau_syn_ex_ ); V_.exp_tau_syn_in_ = std::exp( -V_.h_ms_ / P_.tau_syn_in_ ); - V_.P30_ = -P_.tau_m_ * V_.inv_c_m_ * V_.expm1_tau_m_; + V_.P30_ = -P_.tau_m_ / P_.c_m_ * V_.expm1_tau_m_; // these are determined according to a numeric stability criterion propagator_ex_ = IAFPropagatorAlpha( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_ ); std::tie( V_.P31_ex_, V_.P32_ex_ ) = propagator_ex_.evaluate( V_.h_ms_ );