add timestep() predefined function akin to resolution()

doc/nestml_language/nestml_language_concepts.rst

@@ -559,6 +559,9 @@ The following functions are predefined in NESTML and can be used out of the box.
    * - ``steps``
      - t
      - Convert a time into a number of simulation steps. See the section :ref:`Handling of time` for more information.
+   * - ``resolution``
+     -
+     - In the ``update`` block, or in initialising expressions, returns the current timestep taken in milliseconds. See the section :ref:`Handling of time` for more information.
    * - ``timestep``
      -
      - In the ``update`` block, returns the current timestep taken in milliseconds. See the section :ref:`Handling of time` for more information.
@@ -1104,6 +1107,7 @@ Retrieving simulation timing parameters
 
 To retrieve timing parameters from the simulator kernel, two special functions are built into NESTML:
 
+- ``resolution`` returns the current timestep taken. Can be used only inside the ``update`` block and in intialising expressions. The use of this function assumes that the simulator uses fixed resolution steps, therefore it is recommended to use ``timestep()`` instead in order to make the models more generic.
 - ``timestep`` returns the current timestep taken. Can be used only inside the ``update`` block.
 - ``steps`` takes one parameter of type ``ms`` and returns the number of simulation steps in the current simulation resolution. This only makes sense in case of a fixed simulation resolution (such as in NEST); hence, use of this function is not recommended, because it precludes the models from being compatible with other simulation platforms where a non-constant simulation timestep is used.
 

models/neurons/aeif_cond_alpha_neuron.nestml

@@ -45,7 +45,6 @@ model aeif_cond_alpha_neuron:
         V_m mV = E_L       # Membrane potential
         w pA = 0 pA        # Spike-adaptation current
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
     equations:
         inline V_bounded mV = min(V_m, V_peak) # prevent exponential divergence
@@ -61,6 +60,8 @@ model aeif_cond_alpha_neuron:
         V_m' = (-g_L * (V_bounded - E_L) + I_spike - I_syn_exc - I_syn_inh - w + I_e + I_stim) / C_m
         w' = (a * (V_bounded - E_L) - w) / tau_w
 
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
+
     parameters:
         # membrane parameters
 
@@ -105,23 +106,16 @@ model aeif_cond_alpha_neuron:
         spike
 
     update:
-        if is_refractory:
+        if refr_t > 0 ms:
             # neuron is absolute refractory, do not evolve V_m
-            refr_t -= resolution()
+            integrate_odes(w, refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs (including V_m)
-            integrate_odes()
+            # neuron not refractory
+            integrate_odes(V_m, w)
 
-    onCondition(V_m >= V_peak):
+    onCondition(V_m >= V_peak and refr_t <= 0 ms):
         # threshold crossing
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         V_m = V_reset
         w += b
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false
-

models/neurons/aeif_cond_exp_neuron.nestml

@@ -46,7 +46,6 @@ model aeif_cond_exp_neuron:
         V_m mV = E_L   # Membrane potential
         w pA = 0 pA    # Spike-adaptation current
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
     equations:
         inline V_bounded mV = min(V_m, V_peak)    # prevent exponential divergence
@@ -62,6 +61,8 @@ model aeif_cond_exp_neuron:
         V_m' = (-g_L * (V_bounded - E_L) + I_spike - I_syn_exc - I_syn_inh - w + I_e + I_stim) / C_m
         w' = (a * (V_bounded - E_L) - w) / tau_w
 
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
+
     parameters:
         # membrane parameters
         C_m pF = 281.0 pF       # Membrane Capacitance
@@ -96,23 +97,16 @@ model aeif_cond_exp_neuron:
         spike
 
     update:
-        if is_refractory:
+        if refr_t > 0 ms:
             # neuron is absolute refractory, do not evolve V_m
-            refr_t -= resolution()
-            integrate_odes(w)
+            integrate_odes(w, refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs
+            # neuron not refractory
             integrate_odes(w, V_m)
 
-    onCondition(V_m >= V_peak):
+    onCondition(V_m >= V_th and refr_t <= 0 ms):
         # threshold crossing
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         V_m = V_reset # clamp potential
         w += b
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false

models/neurons/hh_cond_exp_destexhe_neuron.nestml

@@ -39,7 +39,6 @@ model hh_cond_exp_destexhe_neuron:
         V_m mV = E_L #  Membrane potential
         V_m_old mV = E_L #  Membrane potential at the previous timestep
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
         Act_m real =  alpha_m_init / ( alpha_m_init + beta_m_init )
         Act_h real = alpha_h_init / ( alpha_h_init + beta_h_init )
@@ -64,6 +63,7 @@ model hh_cond_exp_destexhe_neuron:
         inline I_syn_inh pA = convolve(g_inh, inh_spikes) * nS * ( V_m - E_inh )
 
         V_m' =( -I_Na - I_K - I_M - I_L - I_syn_exc - I_syn_inh + I_e + I_stim - I_noise) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
         # channel dynamics
         inline V_rel mV = V_m - V_T
@@ -124,8 +124,6 @@ model hh_cond_exp_destexhe_neuron:
     internals:
         D_exc uS**2/ms = 2 * sigma_noise_exc**2 / tau_syn_exc
         D_inh uS**2/ms = 2 * sigma_noise_inh**2 / tau_syn_inh
-        A_exc uS = ((D_exc * tau_syn_exc / 2) * (1 - exp(-2 * resolution() / tau_syn_exc )))**.5
-        A_inh uS = ((D_inh * tau_syn_inh / 2) * (1 - exp(-2 * resolution() / tau_syn_inh )))**.5
 
     input:
         inh_spikes <- inhibitory spike
@@ -137,22 +135,20 @@ model hh_cond_exp_destexhe_neuron:
 
     update:
         # Hodgkin-Huxley type model: ODEs are always integrated, regardless of refractory state
-        g_noise_exc = g_noise_exc0 + (g_noise_exc - g_noise_exc0) * exp(-resolution() / tau_syn_exc) + A_exc * random_normal(0, 1)
-        g_noise_inh = g_noise_inh0 + (g_noise_inh - g_noise_inh0) * exp(-resolution() / tau_syn_inh) + A_inh * random_normal(0, 1)
+        A_exc uS = ((D_exc * tau_syn_exc / 2) * (1 - exp(-2 * timestep() / tau_syn_exc )))**.5
+        A_inh uS = ((D_inh * tau_syn_inh / 2) * (1 - exp(-2 * timestep() / tau_syn_inh )))**.5
+        g_noise_exc = g_noise_exc0 + (g_noise_exc - g_noise_exc0) * exp(-timestep() / tau_syn_exc) + A_exc * random_normal(0, 1)
+        g_noise_inh = g_noise_inh0 + (g_noise_inh - g_noise_inh0) * exp(-timestep() / tau_syn_inh) + A_inh * random_normal(0, 1)
 
         V_m_old mV = V_m
-        integrate_odes()
 
-        if is_refractory:
-            # neuron is absolute refractory, decrease refractoriness timer
-            refr_t -= resolution()
+        if refr_t > 0 ms:
+            # neuron is absolute refractory
+            integrate_odes(Act_m, Act_h, Inact_n, Noninact_p, V_m, refr_t)
+        else:
+            # neuron not refractory
+            integrate_odes(Act_m, Act_h, Inact_n, Noninact_p, V_m)
 
-    onCondition(not is_refractory and V_m > V_T + 30mV and V_m_old > V_m):
+    onCondition(refr_t <= 0 ms and V_m > V_T + 30mV and V_m_old > V_m):
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false

models/neurons/iaf_psc_delta_neuron.nestml

@@ -44,11 +44,11 @@ model iaf_psc_delta_neuron:
     state:
         V_m mV = E_L     # Membrane potential
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
     equations:
         kernel K_delta = delta(t)
         V_m' = -(V_m - E_L) / tau_m + convolve(K_delta, spikes) * (mV / ms) + (I_e + I_stim) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
     parameters:
         tau_m   ms = 10 ms              # Membrane time constant
@@ -71,25 +71,15 @@ model iaf_psc_delta_neuron:
         spike
 
     update:
-        if is_refractory:
+        if refr_t > 0 ms:
             # neuron is absolute refractory, do not evolve V_m
-            refr_t -= resolution()
+            integrate_odes(refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs (including V_m)
-            integrate_odes()
+            # neuron not refractory
+            integrate_odes(V_m)
 
-    onCondition(not is_refractory and V_m >= V_th):
+    onCondition(V_m >= V_th and refr_t <= 0 ms):
         # threshold crossing
         V_m = V_reset
-
-        if refr_T > 0 ms:
-            refr_t = refr_T    # start of the refractory period
-            is_refractory = true
-
+        refr_t = refr_T    # start of the refractory period
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false
-

models/neurons/iaf_psc_exp_neuron.nestml

@@ -57,14 +57,14 @@ model iaf_psc_exp_neuron:
     state:
         V_m mV = E_L    # Membrane potential
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
         I_syn_exc pA = 0 pA
         I_syn_inh pA = 0 pA
 
     equations:
         I_syn_exc' = -I_syn_exc / tau_syn_exc
         I_syn_inh' = -I_syn_inh / tau_syn_inh
         V_m' = -(V_m - E_L) / tau_m + (I_syn_exc - I_syn_inh + I_e + I_stim) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
     parameters:
         C_m pF = 250 pF           # Capacitance of the membrane
@@ -88,28 +88,21 @@ model iaf_psc_exp_neuron:
         spike
 
     update:
-        if is_refractory:
+        if refr_t > 0 ms:
             # neuron is absolute refractory, do not evolve V_m
-            refr_t -= resolution()
-            integrate_odes(I_syn_exc, I_syn_inh)
+            integrate_odes(I_syn_exc, I_syn_inh, refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs (including V_m)
-            integrate_odes()
+            # neuron not refractory
+            integrate_odes(I_syn_exc, I_syn_inh, V_m)
 
     onReceive(exc_spikes):
         I_syn_exc += exc_spikes * pA * s
 
     onReceive(inh_spikes):
         I_syn_inh += inh_spikes * pA * s
 
-    onCondition(not is_refractory and V_m >= V_th):
+    onCondition(V_m >= V_th and refr_t <= 0 ms):
         # threshold crossing
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         V_m = V_reset
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false

models/neurons/wb_cond_multisyn_neuron.nestml

@@ -31,7 +31,6 @@ model wb_cond_multisyn_neuron:
         V_m mV = -65. mV # Membrane potential
         V_m_old mV = E_L    # Membrane potential at previous timestep for threshold check
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
         Inact_h real = alpha_h_init / ( alpha_h_init + beta_h_init )    # Inactivation variable h for Na
         Act_n real = alpha_n_init / (alpha_n_init + beta_n_init)  # Activation variable n for K
@@ -138,23 +137,13 @@ model wb_cond_multisyn_neuron:
         spike
 
     update:
-        if is_refractory:
-            refr_t -= resolution()
-
         V_m_old = V_m
-        integrate_odes()
+        integrate_odes()    # in this model, V_m is always integrated, regardless of refractory state
 
-    onCondition(not is_refractory and V_m > V_Tr and V_m_old > V_m):
-        # threshold crossing and maximum
+    onCondition(refr_t <= 0 ms and V_m > V_Tr and V_m_old > V_m):        # if neuron is not refractory, threshold crossing and maximum
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         emit_spike()
 
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false
-
     function compute_synapse_constant(Tau_1 ms, Tau_2 ms, g_peak nS) real:
         # Factor used to account for the missing 1/((1/Tau_2)-(1/Tau_1)) term
         # in the ht_neuron_dynamics integration of the synapse terms.

pynestml/cocos/co_co_resolution_func_used.py

@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+#
+# co_co_resolution_func_used.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/>.
+
+from pynestml.cocos.co_co import CoCo
+from pynestml.meta_model.ast_model import ASTModel
+from pynestml.meta_model.ast_simple_expression import ASTSimpleExpression
+from pynestml.symbols.predefined_functions import PredefinedFunctions
+from pynestml.utils.logger import LoggingLevel, Logger
+from pynestml.utils.messages import Messages
+from pynestml.visitors.ast_visitor import ASTVisitor
+
+
+class CoCoResolutionFuncUsed(CoCo):
+    r"""
+    This Coco emits a warning in case the ``resolution()`` predefined function is used.
+    """
+
+    @classmethod
+    def check_co_co(cls, model: ASTModel):
+        """
+        Checks the coco.
+        :param model: a single neuron
+        """
+        visitor = CoCoResolutionFuncUsedVisitor()
+        visitor.neuron = model
+        model.accept(visitor)
+
+
+class CoCoResolutionFuncUsedVisitor(ASTVisitor):
+    def visit_simple_expression(self, node):
+        if node.get_function_call() is None:
+            return
+
+        function_name = node.get_function_call().get_name()
+        if function_name == PredefinedFunctions.TIME_RESOLUTION:
+            code, message = Messages.get_resolution_func_used()
+            Logger.log_message(code=code, message=message, error_position=node.get_source_position(), log_level=LoggingLevel.WARNING)

pynestml/cocos/co_cos_manager.py

@@ -52,6 +52,7 @@
 from pynestml.cocos.co_co_internals_assigned_only_in_internals_block import CoCoInternalsAssignedOnlyInInternalsBlock
 from pynestml.cocos.co_co_parameters_assigned_only_in_parameter_block import CoCoParametersAssignedOnlyInParameterBlock
 from pynestml.cocos.co_co_resolution_func_legally_used import CoCoResolutionFuncLegallyUsed
+from pynestml.cocos.co_co_resolution_func_used import CoCoResolutionFuncUsed
 from pynestml.cocos.co_co_simple_delta_function import CoCoSimpleDeltaFunction
 from pynestml.cocos.co_co_state_variables_initialized import CoCoStateVariablesInitialized
 from pynestml.cocos.co_co_convolve_has_correct_parameter import CoCoConvolveHasCorrectParameter
@@ -328,6 +329,14 @@ def check_invariant_type_correct(cls, model: ASTModel):
         """
         CoCoInvariantIsBoolean.check_co_co(model)
 
+    @classmethod
+    def check_resolution_func_used(cls, model: ASTModel):
+        """
+        Checks if all invariants are of type boolean.
+        :param model: a single model object.
+        """
+        CoCoResolutionFuncUsed.check_co_co(model)
+
     @classmethod
     def check_vector_in_non_vector_declaration_detected(cls, model: ASTModel):
         """
@@ -451,4 +460,5 @@ def post_symbol_table_builder_checks(cls, model: ASTModel, after_ast_rewrite: bo
         cls.check_vector_declaration_size(model)
         cls.check_co_co_priorities_correctly_specified(model)
         cls.check_resolution_func_legally_used(model)
+        cls.check_resolution_func_used(model)
         cls.check_input_port_size_type(model)

pynestml/utils/messages.py

@@ -130,6 +130,7 @@ class MessageCode(Enum):
     SYNS_BAD_BUFFER_COUNT = 107
     CM_NO_V_COMP = 108
     MECHS_DICTIONARY_INFO = 109
+    RESOLUTION_FUNC_USED = 110
 
 
 class Messages:
@@ -1312,3 +1313,8 @@ def get_mechs_dictionary_info(cls, chan_info, syns_info, conc_info, con_in_info)
         message += "con_in_info:\n" + con_in_info + "\n"
 
         return MessageCode.MECHS_DICTIONARY_INFO, message
+
+    @classmethod
+    def get_resolution_func_used(cls):
+        message = "Model contains a call to the ``resolution()`` function. This restricts the model to being compatible only with fixed-timestep simulators. Consider using ``timestep()`` instead."
+        return MessageCode.RESOLUTION_FUNC_USED, message