Merge pull request #3305 from heshpdx/iaf_psc_model

Optimize computation in iaf_psc_alpha_ps model
nest · Sep 27, 2024 · f0da7ab · f0da7ab
2 parents e61f34f + e63b1e9
commit f0da7ab
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 8 deletions.
diff --git a/models/iaf_psc_alpha_ps.cpp b/models/iaf_psc_alpha_ps.cpp
@@ -268,8 +268,14 @@ 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_;
+  // 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_ = 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_ / P_.tau_m_ );
@@ -500,9 +506,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 * V_.inv_tau_m_ );
 
-    const double ps_P30 = -P_.tau_m_ / P_.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;
@@ -518,8 +524,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 * 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_;
@@ -578,9 +584,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 * V_.inv_tau_m_ );
 
-  const double ps_P30 = -P_.tau_m_ / P_.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
@@ -407,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 -------------------------------