From b6e73ce102d0d82e24be9f8fe2f50d3ba0da9dd0 Mon Sep 17 00:00:00 2001 From: "C.A.P. Linssen" Date: Fri, 20 Sep 2024 20:21:03 +0200 Subject: [PATCH] add timestep() predefined function akin to resolution() --- .../nestml_language_concepts.rst | 4 ++ models/neurons/aeif_cond_alpha_neuron.nestml | 20 +++---- models/neurons/aeif_cond_exp_neuron.nestml | 18 ++---- .../hh_cond_exp_destexhe_neuron.nestml | 28 ++++------ models/neurons/iaf_psc_delta_neuron.nestml | 24 +++----- models/neurons/iaf_psc_exp_neuron.nestml | 19 ++----- models/neurons/wb_cond_multisyn_neuron.nestml | 15 +---- pynestml/cocos/co_co_resolution_func_used.py | 55 +++++++++++++++++++ pynestml/cocos/co_cos_manager.py | 10 ++++ pynestml/utils/messages.py | 6 ++ 10 files changed, 115 insertions(+), 84 deletions(-) create mode 100644 pynestml/cocos/co_co_resolution_func_used.py diff --git a/doc/nestml_language/nestml_language_concepts.rst b/doc/nestml_language/nestml_language_concepts.rst index 226eb8133..e25824dff 100644 --- a/doc/nestml_language/nestml_language_concepts.rst +++ b/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. diff --git a/models/neurons/aeif_cond_alpha_neuron.nestml b/models/neurons/aeif_cond_alpha_neuron.nestml index 4b35dbd6b..6640a38f4 100644 --- a/models/neurons/aeif_cond_alpha_neuron.nestml +++ b/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 - diff --git a/models/neurons/aeif_cond_exp_neuron.nestml b/models/neurons/aeif_cond_exp_neuron.nestml index 0d2a6e6b3..b3dd70734 100644 --- a/models/neurons/aeif_cond_exp_neuron.nestml +++ b/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 diff --git a/models/neurons/hh_cond_exp_destexhe_neuron.nestml b/models/neurons/hh_cond_exp_destexhe_neuron.nestml index e00e2c401..c2bfced5a 100644 --- a/models/neurons/hh_cond_exp_destexhe_neuron.nestml +++ b/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 diff --git a/models/neurons/iaf_psc_delta_neuron.nestml b/models/neurons/iaf_psc_delta_neuron.nestml index f3d4005fd..ca0897d0d 100644 --- a/models/neurons/iaf_psc_delta_neuron.nestml +++ b/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 - diff --git a/models/neurons/iaf_psc_exp_neuron.nestml b/models/neurons/iaf_psc_exp_neuron.nestml index 882cbf14f..f20ed3d77 100644 --- a/models/neurons/iaf_psc_exp_neuron.nestml +++ b/models/neurons/iaf_psc_exp_neuron.nestml @@ -57,7 +57,6 @@ 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 @@ -65,6 +64,7 @@ model iaf_psc_exp_neuron: 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,13 +88,12 @@ 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 @@ -102,14 +101,8 @@ model iaf_psc_exp_neuron: 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 diff --git a/models/neurons/wb_cond_multisyn_neuron.nestml b/models/neurons/wb_cond_multisyn_neuron.nestml index 6f6886557..de8de886b 100644 --- a/models/neurons/wb_cond_multisyn_neuron.nestml +++ b/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. diff --git a/pynestml/cocos/co_co_resolution_func_used.py b/pynestml/cocos/co_co_resolution_func_used.py new file mode 100644 index 000000000..9127cc96c --- /dev/null +++ b/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 . + +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) diff --git a/pynestml/cocos/co_cos_manager.py b/pynestml/cocos/co_cos_manager.py index 01d008890..8cfcae7c6 100644 --- a/pynestml/cocos/co_cos_manager.py +++ b/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) diff --git a/pynestml/utils/messages.py b/pynestml/utils/messages.py index 3d5673081..dec92bea1 100644 --- a/pynestml/utils/messages.py +++ b/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