Merge pull request #3304 from heshpdx/master

Optimize computation in GrowthCurveGaussian
nest · Sep 5, 2024 · a343246 · a343246
2 parents ef0074a + 2a8c7f8
commit a343246
Show file tree

Hide file tree

Showing 4 changed files with 21 additions and 6 deletions.
diff --git a/libnestutil/numerics.cpp b/libnestutil/numerics.cpp
@@ -93,6 +93,8 @@ const double numerics::nan = NAN;
 const double numerics::nan = 0.0 / 0.0;
 #endif
 
+const double numerics::sqrt_log_two = std::sqrt( std::log( 2.0 ) );
+
 // later also in namespace
 long
 ld_round( double x )

diff --git a/libnestutil/numerics.h b/libnestutil/numerics.h
@@ -47,6 +47,7 @@ namespace numerics
 extern const double e;
 extern const double pi;
 extern const double nan;
+extern const double sqrt_log_two;
 
 /** Supply expm1() function independent of system.
  *  @note Implemented inline for efficiency.

diff --git a/nestkernel/growth_curve.cpp b/nestkernel/growth_curve.cpp
@@ -78,6 +78,7 @@ nest::GrowthCurveGaussian::GrowthCurveGaussian()
   , eta_( 0.1 )
   , eps_( 0.7 )
 {
+  compute_local_();
 }
 
 void
@@ -93,6 +94,7 @@ nest::GrowthCurveGaussian::set( const DictionaryDatum& d )
 {
   updateValue< double >( d, names::eps, eps_ );
   updateValue< double >( d, names::eta, eta_ );
+  compute_local_();
 }
 
 double
@@ -106,22 +108,28 @@ nest::GrowthCurveGaussian::update( double t,
   // Numerical integration from t_minus to t
   // use standard forward Euler numerics
   const double h = Time::get_resolution().get_ms();
-  const double zeta = ( eta_ - eps_ ) / ( 2.0 * sqrt( log( 2.0 ) ) );
-  const double xi = ( eta_ + eps_ ) / 2.0;
+  const double inv_tau_Ca = 1.0 / tau_Ca;
 
   double z_value = z_minus;
   double Ca = Ca_minus;
 
-  for ( double lag = t_minus; lag < ( t - h / 2.0 ); lag += h )
+  for ( double lag = t_minus; lag < ( t - h * 0.5 ); lag += h )
   {
-    Ca = Ca - ( ( Ca / tau_Ca ) * h );
-    const double dz = h * growth_rate * ( 2.0 * exp( -pow( ( Ca - xi ) / zeta, 2 ) ) - 1.0 );
-    z_value = z_value + dz;
+    Ca = Ca - ( ( Ca * inv_tau_Ca ) * h );
+    const double dz = h * growth_rate * ( 2.0 * std::exp( -std::pow( ( Ca - xi_ ) * inv_zeta_, 2 ) ) - 1.0 );
+    z_value += dz;
   }
 
   return std::max( z_value, 0.0 );
 }
 
+void
+nest::GrowthCurveGaussian::compute_local_()
+{
+  inv_zeta_ = 2.0 * numerics::sqrt_log_two / ( eta_ - eps_ );
+  xi_ = ( eta_ + eps_ ) * 0.5;
+}
+
 /* ----------------------------------------------------------------
  * GrowthCurveSigmoid
  * ---------------------------------------------------------------- */

diff --git a/nestkernel/growth_curve.h b/nestkernel/growth_curve.h
@@ -227,8 +227,12 @@ class GrowthCurveGaussian : public GrowthCurve
   update( double t, double t_minus, double Ca_minus, double z, double tau_Ca, double growth_rate ) const override;
 
 private:
+  void compute_local_();
+
   double eta_;
   double eps_;
+  double inv_zeta_;
+  double xi_;
 };
 
 /** @BeginDocumentation