Add absolute and relative GSL solver parameters (#1086)

Co-authored-by: C.A.P. Linssen <charl@turingbirds.com>

pynestml/codegeneration/nest_code_generator.py

@@ -110,6 +110,7 @@ class NESTCodeGenerator(CodeGenerator):
         - **module_templates**: A list of the jinja templates or a relative path to a directory containing the templates related to generating the NEST module.
     - **nest_version**: A string identifying the version of NEST Simulator to generate code for. The string corresponds to the NEST Simulator git repository tag or git branch name, for instance, ``"v2.20.2"`` or ``"master"``. The default is the empty string, which causes the NEST version to be automatically identified from the ``nest`` Python module.
     - **solver**: A string identifying the preferred ODE solver. ``"analytic"`` for propagator solver preferred; fallback to numeric solver in case ODEs are not analytically solvable. Use ``"numeric"`` to disable analytic solver.
+    - **gsl_adaptive_step_size_controller**: For the numeric (GSL) solver: how to interpret the absolute and relative tolerance values. Can be changed to trade off integration accuracy with numerical stability. The default value is ``"with_respect_to_solution"``. Can also be set to ``"with_respect_to_derivative"``. (Tolerance values can be specified at runtime as parameters of the model instance.) For further details, see https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control.
     - **numeric_solver**: A string identifying the preferred numeric ODE solver. Supported are ``"rk45"`` and ``"forward-Euler"``.
     - **delay_variable**: A mapping identifying, for each synapse (the name of which is given as a key), the variable or parameter in the model that corresponds with the NEST ``Connection`` class delay property.
     - **weight_variable**: Like ``delay_variable``, but for synaptic weight.
@@ -141,6 +142,7 @@ class NESTCodeGenerator(CodeGenerator):
         },
         "nest_version": "",
         "solver": "analytic",
+        "gsl_adaptive_step_size_controller": "with_respect_to_solution",
         "numeric_solver": "rk45",
         "delay_variable": {},
         "weight_variable": {}
@@ -489,6 +491,7 @@ def _get_model_namespace(self, astnode: ASTModel) -> Dict:
 
         if namespace["uses_numeric_solver"]:
             namespace["numeric_solver"] = self.get_option("numeric_solver")
+            namespace["gsl_adaptive_step_size_controller"] = self.get_option("gsl_adaptive_step_size_controller")
 
         return namespace
 

pynestml/codegeneration/resources_nest/point_neuron/common/NeuronClass.jinja2

@@ -353,9 +353,9 @@ void {{ neuronName }}::init_state_internal_()
   const double __resolution = nest::Time::get_resolution().get_ms();  // do not remove, this is necessary for the resolution() function
 
 {%- if numeric_solver == "rk45" %}
-  // by default, integrate all variables
-  // use a default "good enough" value for the absolute error. It can be adjusted via `node.set()`
-  P_.__gsl_error_tol = 1e-3;
+  // by default, integrate all variables with a conservative tolerance, in the sense that we err on the side of integrating very precisely at the expense of extra computation
+  P_.__gsl_abs_error_tol = 1e-6;
+  P_.__gsl_rel_error_tol = 1e-6;
 {%- endif %}
 
 {%- if parameter_vars_with_iv|length > 0 %}
@@ -463,11 +463,24 @@ void {{ neuronName }}::init_buffers_()
 
   if ( not B_.__c )
   {
-    B_.__c = gsl_odeiv_control_y_new( P_.__gsl_error_tol, 0.0 );
+{%- if gsl_adaptive_step_size_controller == "with_respect_to_solution" %}
+    B_.__c = gsl_odeiv_control_y_new( P_.__gsl_abs_error_tol, P_.__gsl_rel_error_tol );
+{%- elif gsl_adaptive_step_size_controller == "with_respect_to_derivative" %}
+    B_.__c = gsl_odeiv_control_yp_new( P_.__gsl_abs_error_tol, P_.__gsl_rel_error_tol );
+{%- else %}
+{{ raise('Unknown GSL adaptive step size controller value') }}
+{%- endif %}
   }
   else
   {
-    gsl_odeiv_control_init( B_.__c, P_.__gsl_error_tol, 0.0, 1.0, 0.0 );
+{%- if gsl_adaptive_step_size_controller == "with_respect_to_solution" %}
+    gsl_odeiv_control_init( B_.__c, P_.__gsl_abs_error_tol, P_.__gsl_rel_error_tol, 1.0, 0.0 );
+{%- elif gsl_adaptive_step_size_controller == "with_respect_to_derivative" %}
+    gsl_odeiv_control_init( B_.__c, P_.__gsl_abs_error_tol, P_.__gsl_rel_error_tol, 0.0, 1.0 );
+{%- else %}
+{{ raise('Unknown GSL adaptive step size controller value') }}
+{%- endif %}
+
   }
 
   if ( not B_.__e )

pynestml/codegeneration/resources_nest/point_neuron/common/NeuronHeader.jinja2

@@ -119,6 +119,9 @@ namespace {{names_namespace}}
     const Name _{{sym.get_symbol_name()}}( "{{sym.get_symbol_name()}}" );
 {%- endfor %}
 {%- endif %}
+
+    const Name gsl_abs_error_tol("gsl_abs_error_tol");
+    const Name gsl_rel_error_tol("gsl_rel_error_tol");
 }
 }
 
@@ -582,7 +585,8 @@ private:
 {%- if uses_numeric_solver %}
 {%-     if numeric_solver == "rk45" %}
 
-    double __gsl_error_tol;
+    double __gsl_abs_error_tol;
+    double __gsl_rel_error_tol;
 {%-     endif %}
 {%- endif %}
 
@@ -1003,9 +1007,13 @@ inline void {{neuronName}}::get_status(DictionaryDatum &__d) const
   (*__d)[nest::names::recordables] = recordablesMap_.get_list();
 {%- if uses_numeric_solver %}
 {%-     if numeric_solver == "rk45" %}
-  def< double >(__d, nest::names::gsl_error_tol, P_.__gsl_error_tol);
-  if ( P_.__gsl_error_tol <= 0. ){
-    throw nest::BadProperty( "The gsl_error_tol must be strictly positive." );
+  def< double >(__d, nest::{{ neuronName }}_names::gsl_abs_error_tol, P_.__gsl_abs_error_tol);
+  if ( P_.__gsl_abs_error_tol <= 0. ){
+    throw nest::BadProperty( "The gsl_abs_error_tol must be strictly positive." );
+  }
+  def< double >(__d, nest::{{ neuronName }}_names::gsl_rel_error_tol, P_.__gsl_rel_error_tol);
+  if ( P_.__gsl_rel_error_tol < 0. ){
+    throw nest::BadProperty( "The gsl_rel_error_tol must be zero or positive." );
   }
 {%-     endif %}
 {%- endif %}
@@ -1063,10 +1071,15 @@ inline void {{neuronName}}::set_status(const DictionaryDatum &__d)
 
 {% if uses_numeric_solver %}
 {%-     if numeric_solver == "rk45" %}
-  updateValue< double >(__d, nest::names::gsl_error_tol, P_.__gsl_error_tol);
-  if ( P_.__gsl_error_tol <= 0. )
+  updateValue< double >(__d, nest::{{ neuronName }}_names::gsl_abs_error_tol, P_.__gsl_abs_error_tol);
+  if ( P_.__gsl_abs_error_tol <= 0. )
+  {
+    throw nest::BadProperty( "The gsl_abs_error_tol must be strictly positive." );
+  }
+  updateValue< double >(__d, nest::{{ neuronName }}_names::gsl_rel_error_tol, P_.__gsl_rel_error_tol);
+  if ( P_.__gsl_rel_error_tol < 0. )
   {
-    throw nest::BadProperty( "The gsl_error_tol must be strictly positive." );
+    throw nest::BadProperty( "The gsl_rel_error_tol must be zero or positive." );
   }
 {%-     endif %}
 {%- endif %}

pynestml/codegeneration/resources_nest/point_neuron/common/SynapseHeader.h.jinja2

@@ -339,12 +339,6 @@ private:
 {%- endfor %}
 {%- endfilter %}
 
-{% if uses_numeric_solver %}
-{%-     if numeric_solver == "rk45" %}
-    double __gsl_error_tol;
-{%     endif %}
-{% endif %}
-
     /** Initialize parameters to their default values. */
     Parameters_() {};
   };
@@ -1082,17 +1076,6 @@ set_delay(tmp_{{ nest_codegen_opt_delay_variable }});
 {%- endfor %}
 {%- endif %}
 
-{%  if uses_numeric_solver %}
-{%-     if numeric_solver == "rk45" %}
-
-  updateValue< double >(__d, nest::names::gsl_error_tol, P_.__gsl_error_tol);
-  if ( P_.__gsl_error_tol <= 0. )
-  {
-    throw nest::BadProperty( "The gsl_error_tol must be strictly positive." );
-  }
-{%-    endif %}
-{%  endif %}
-
   // recompute internal variables in case they are dependent on parameters or state that might have been updated in this call to set_status()
   V_.__h = nest::Time::get_resolution().get_ms();
   recompute_internal_variables();

pynestml/codegeneration/resources_spinnaker/@NEURON_NAME@_impl.h.jinja2

@@ -133,8 +133,6 @@ typedef struct neuron_impl_t {
 // TODO: implement gsl solver support and find out where to put gsl
 /***
 {%- if uses_numeric_solver %}
-    REAL __gsl_error_tol;
-
 
     // -----------------------------------------------------------------------
     //   GSL ODE solver data structures

tests/nest_tests/nest_integration_test.py

@@ -129,12 +129,14 @@ def test_nest_integration(self):
         self._test_model_equivalence_spiking("izhikevich", "izhikevich_neuron_alt_int_nestml")
         self._test_model_equivalence_fI_curve("izhikevich", "izhikevich_neuron_alt_int_nestml")
 
-        self._test_model_equivalence_subthreshold("hh_psc_alpha", "hh_psc_alpha_neuron_nestml")
-        self._test_model_equivalence_spiking("hh_psc_alpha", "hh_psc_alpha_neuron_nestml", tolerance=1E-5)
-        self._test_model_equivalence_fI_curve("hh_psc_alpha", "hh_psc_alpha_neuron_nestml")
-
-        self._test_model_equivalence_subthreshold("hh_cond_exp_traub", "hh_cond_exp_traub_neuron_nestml")
-        self._test_model_equivalence_fI_curve("hh_cond_exp_traub", "hh_cond_exp_traub_neuron_nestml")
+        nestml_hh_psc_alpha_model_parameters = {"gsl_abs_error_tol": 1E-3, "gsl_rel_error_tol": 0.}  # matching the defaults in NEST
+        self._test_model_equivalence_subthreshold("hh_psc_alpha", "hh_psc_alpha_neuron_nestml", nestml_model_parameters=nestml_hh_psc_alpha_model_parameters)
+        self._test_model_equivalence_spiking("hh_psc_alpha", "hh_psc_alpha_neuron_nestml", tolerance=1E-5, nestml_model_parameters=nestml_hh_psc_alpha_model_parameters)
+        self._test_model_equivalence_fI_curve("hh_psc_alpha", "hh_psc_alpha_neuron_nestml", nestml_model_parameters=nestml_hh_psc_alpha_model_parameters)
+
+        nestml_hh_cond_exp_traub_model_parameters = {"gsl_abs_error_tol": 1E-3, "gsl_rel_error_tol": 0.}  # matching the defaults in NEST
+        self._test_model_equivalence_subthreshold("hh_cond_exp_traub", "hh_cond_exp_traub_neuron_nestml", nestml_model_parameters=nestml_hh_cond_exp_traub_model_parameters)
+        self._test_model_equivalence_fI_curve("hh_cond_exp_traub", "hh_cond_exp_traub_neuron_nestml", nestml_model_parameters=nestml_hh_cond_exp_traub_model_parameters)
 
         self._test_model_equivalence_subthreshold("aeif_cond_exp", "aeif_cond_exp_neuron_alt_nestml", kernel_opts={"resolution": .01})    # needs resolution 0.01 because the NEST model overrides this internally. Subthreshold only because threshold detection is inside the while...gsl_odeiv_evolve_apply() loop in NEST but outside the loop (strictly after gsl_odeiv_evolve_apply()) in NESTML, causing spike times to differ slightly
         self._test_model_equivalence_fI_curve("aeif_cond_exp", "aeif_cond_exp_neuron_alt_nestml")

tests/nest_tests/test_aeif_integrator.py

@@ -0,0 +1,84 @@
+# -*- coding: utf-8 -*-
+#
+# test_aeif_integrator.py
+#
+# 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 <http://www.gnu.org/licenses/>.
+
+import numpy as np
+import os
+import pytest
+
+import nest
+from nest.lib.hl_api_exceptions import NESTErrors
+
+from pynestml.frontend.pynestml_frontend import generate_nest_target
+from pynestml.codegeneration.nest_tools import NESTTools
+
+
+@pytest.mark.skipif(NESTTools.detect_nest_version().startswith("v2"),
+                    reason="This test does not support NEST 2")
+class TestAEIFIntegrator_with_respect_to_solution:
+    r"""Check that the integrator default tolerances are set such that the adaptive exponential models do not cause trouble for the numerical integrator"""
+
+    @pytest.mark.parametrize("gsl_adaptive_step_size_controller", ["with_respect_to_solution", "with_respect_to_derivative"])
+    def test_aeif_integrator(self, gsl_adaptive_step_size_controller: str):
+        # generate the code
+
+        input_path = [os.path.realpath(os.path.join(os.path.dirname(__file__), os.path.join(os.pardir, os.pardir, "models", "neurons", "aeif_cond_exp_neuron.nestml")))]
+        logging_level = "DEBUG"
+        suffix = "_nestml"
+
+        nest.set_verbosity("M_ALL")
+
+        generate_nest_target(input_path,
+                             logging_level=logging_level,
+                             suffix=suffix,
+                             codegen_opts={"gsl_adaptive_step_size_controller": gsl_adaptive_step_size_controller})
+
+        # run the simulation
+
+        nest.ResetKernel()
+        nest.resolution = 0.01
+        nest.Install("nestmlmodule")
+
+        params = {
+            'E_L': -80,
+            'V_reset': -60.0,
+            'refr_T': 2.0,
+            'g_L': 4,
+            'C_m': 32,
+            'E_exc': 0.0,
+            'E_inh': -80.0,
+            'tau_syn_exc': 1.5,
+            'tau_syn_inh': 4.2,
+            'a': 0.8,
+            'b': 80,
+            'Delta_T': 0.8,
+            'tau_w': 144.0,
+            'V_th': -57.0
+        }
+
+        neuron = nest.Create("aeif_cond_exp_neuron_nestml")
+        neuron.set(**params)
+        neuron.set({'V_m': -54.})
+
+        # the test succeeds if the integrator terminates
+
+        nest.Simulate(100.)
+
+        assert not np.isnan(neuron.V_m)