Add CoCo to ensure defining expressions only refer from parameters ->internals -> state blocks

.gitignore

@@ -1,3 +1,4 @@
+*.orig
 .settings/
 target/
 .idea/

models/neurons/hh_cond_exp_destexhe_neuron.nestml

@@ -54,9 +54,9 @@ model hh_cond_exp_destexhe_neuron:
         g_noise_exc uS = g_noise_exc0
         g_noise_inh uS = g_noise_inh0
 
-        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
+        V_m mV = V_m_init        # Membrane potential
+        V_m_old mV = V_m_init    # Membrane potential at the previous timestep
+        refr_t ms = 0 ms         # Refractory period timer
 
         Act_m real = alpha_m_init / (alpha_m_init + beta_m_init)
         Act_h real = alpha_h_init / (alpha_h_init + beta_h_init)
@@ -109,6 +109,7 @@ model hh_cond_exp_destexhe_neuron:
         E_Na mV = 60 mV            # Reversal potentials
         E_K mV = -90 mV            # Potassium reversal potential
         E_L mV = -80 mV            # Leak reversal potential (a.k.a. resting potential)
+        V_m_init mV = E_L          # Initial membrane potential
         V_T mV = -58 mV            # Voltage offset that controls dynamics. For default
                                    # parameters, V_T = -63mV results in a threshold around -50mV.
         tau_syn_exc ms = 2.7 ms    # Synaptic time constant for excitatory synapse

models/neurons/traub_cond_multisyn_neuron.nestml

@@ -44,7 +44,7 @@
 #
 model traub_cond_multisyn_neuron:
     state:
-        V_m mV = -70. mV    # Membrane potential
+        V_m mV = V_m_init   # Membrane potential
         V_m_old mV = E_L    # Membrane potential at previous timestep for threshold check
         refr_t ms = 0 ms    # Refractory period timer
 
@@ -107,6 +107,7 @@ model traub_cond_multisyn_neuron:
                g_GABAB$' = -g_GABAB$ / tau_GABAB_1
 
     parameters:
+        V_m_init mV = -70 mV    # Initial membrane potential
         C_m pF = 100.0 pF       # Membrane capacitance
         g_Na nS = 10000.0 nS    # Sodium peak conductance
         g_K nS = 8000.0 nS      # Potassium peak conductance
@@ -149,12 +150,12 @@ model traub_cond_multisyn_neuron:
         GABA_AInitialValue real = compute_synapse_constant(tau_GABAA_1, tau_GABAA_2, GABA_A_g_peak)
         GABA_BInitialValue real = compute_synapse_constant(tau_GABAB_1, tau_GABAB_2, GABA_B_g_peak)
 
-        alpha_n_init real = 0.032 * (V_m / mV + 52.) / (1. - exp(-(V_m / mV + 52.) / 5.))
-        beta_n_init real = 0.5 * exp(-(V_m / mV + 57.) / 40.)
-        alpha_m_init real = 0.32 * (V_m / mV + 54.) / (1.0 - exp(-(V_m / mV + 54.) / 4.))
-        beta_m_init real = 0.28 * (V_m / mV + 27.) / (exp((V_m / mV + 27.) / 5.) - 1.)
-        alpha_h_init real = 0.128 * exp(-(V_m / mV + 50.0) / 18.0)
-        beta_h_init real = 4.0 / (1.0 + exp(-(V_m / mV + 27.) / 5.))
+        alpha_n_init real = 0.032 * (V_m_init / mV + 52.) / (1. - exp(-(V_m_init / mV + 52.) / 5.))
+        beta_n_init real = 0.5 * exp(-(V_m_init / mV + 57.) / 40.)
+        alpha_m_init real = 0.32 * (V_m_init / mV + 54.) / (1.0 - exp(-(V_m_init / mV + 54.) / 4.))
+        beta_m_init real = 0.28 * (V_m_init / mV + 27.) / (exp((V_m_init / mV + 27.) / 5.) - 1.)
+        alpha_h_init real = 0.128 * exp(-(V_m_init / mV + 50.0) / 18.0)
+        beta_h_init real = 4.0 / (1.0 + exp(-(V_m_init / mV + 27.) / 5.))
 
     input:
         AMPA <- spike

models/neurons/traub_psc_alpha_neuron.nestml

@@ -94,12 +94,12 @@ model traub_psc_alpha_neuron:
         I_e pA = 0 pA
 
     internals:
-        alpha_n_init real = 0.032 * (V_m / mV + 52.) / (1. - exp(-(V_m / mV + 52.) / 5.))
-        beta_n_init real = 0.5 * exp(-(V_m / mV + 57.) / 40.)
-        alpha_m_init real = 0.32 * (V_m / mV + 54.) / (1.0 - exp(-(V_m / mV + 54.) / 4.))
-        beta_m_init real = 0.28 * (V_m / mV + 27.) / (exp((V_m / mV + 27.) / 5.) - 1.)
-        alpha_h_init real = 0.128 * exp(-(V_m / mV + 50.0) / 18.0)
-        beta_h_init real = 4.0 / (1.0 + exp(-(V_m / mV + 27.) / 5.))
+        alpha_n_init real = 0.032 * (V_m_init / mV + 52.) / (1. - exp(-(V_m_init / mV + 52.) / 5.))
+        beta_n_init real = 0.5 * exp(-(V_m_init / mV + 57.) / 40.)
+        alpha_m_init real = 0.32 * (V_m_init / mV + 54.) / (1.0 - exp(-(V_m_init / mV + 54.) / 4.))
+        beta_m_init real = 0.28 * (V_m_init / mV + 27.) / (exp((V_m_init / mV + 27.) / 5.) - 1.)
+        alpha_h_init real = 0.128 * exp(-(V_m_init / mV + 50.0) / 18.0)
+        beta_h_init real = 4.0 / (1.0 + exp(-(V_m_init / mV + 27.) / 5.))
 
     input:
         exc_spikes <- excitatory spike

models/neurons/wb_cond_multisyn_neuron.nestml

@@ -47,7 +47,7 @@
 #
 model wb_cond_multisyn_neuron:
     state:
-        V_m mV = -65. mV    # Membrane potential
+        V_m mV = V_m_init    # Membrane potential
         V_m_old mV = E_L    # Membrane potential at previous timestep for threshold check
         refr_t ms = 0 ms    # Refractory period timer
 
@@ -104,6 +104,7 @@ model wb_cond_multisyn_neuron:
         E_Na mV = 55.0 mV      # Sodium reversal potential
         E_K mV = -90.0 mV      # Potassium reversal potential
         E_L mV = -65.0 mV      # Leak reversal potential (a.k.a. resting potential)
+        V_m_init mV = E_L      # Initial membrane potential
         V_Tr mV = -55.0 mV     # Spike threshold
         refr_T ms = 2 ms       # Duration of refractory period
 
@@ -140,12 +141,12 @@ model wb_cond_multisyn_neuron:
         GABA_AInitialValue real = compute_synapse_constant(GABA_A_Tau_1, GABA_A_Tau_2, GABA_A_g_peak)
         GABA_BInitialValue real = compute_synapse_constant(GABA_B_Tau_1, GABA_B_Tau_2, GABA_B_g_peak)
 
-        alpha_n_init real = -0.05 * (V_m / mV + 34.0) / (exp(-0.1 * (V_m / mV + 34.0)) - 1.0)
-        beta_n_init real = 0.625 * exp(-(V_m / mV + 44.0) / 80.0)
-        alpha_m_init real = 0.1 * (V_m / mV + 35.0) / (1.0 - exp(-0.1 * (V_m / mV + 35.)))
-        beta_m_init  real = 4.0 * exp(-(V_m / mV + 60.0) / 18.0)
-        alpha_h_init real = 0.35 * exp(-(V_m / mV + 58.0) / 20.0)
-        beta_h_init  real = 5.0 / (exp(-0.1 * (V_m / mV + 28.0)) + 1.0)
+        alpha_n_init real = -0.05 * (V_m_init / mV + 34.0) / (exp(-0.1 * (V_m_init / mV + 34.0)) - 1.0)
+        beta_n_init real = 0.625 * exp(-(V_m_init / mV + 44.0) / 80.0)
+        alpha_m_init real = 0.1 * (V_m_init / mV + 35.0) / (1.0 - exp(-0.1 * (V_m_init / mV + 35.)))
+        beta_m_init  real = 4.0 * exp(-(V_m_init / mV + 60.0) / 18.0)
+        alpha_h_init real = 0.35 * exp(-(V_m_init / mV + 58.0) / 20.0)
+        beta_h_init  real = 5.0 / (exp(-0.1 * (V_m_init / mV + 28.0)) + 1.0)
 
     input:
         AMPA <- spike

pynestml/cocos/co_co_parameters_internals_state_declaration_order.py

@@ -0,0 +1,153 @@
+# -*- coding: utf-8 -*-
+#
+# co_co_parameters_internals_state_declaration_order.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 typing import Optional
+
+from pynestml.cocos.co_co import CoCo
+from pynestml.meta_model.ast_block_with_variables import ASTBlockWithVariables
+from pynestml.meta_model.ast_declaration import ASTDeclaration
+from pynestml.meta_model.ast_model import ASTModel
+from pynestml.meta_model.ast_variable import ASTVariable
+from pynestml.symbols.predefined_variables import PredefinedVariables
+from pynestml.symbols.symbol import SymbolKind
+from pynestml.symbols.variable_symbol import BlockType, VariableSymbol
+from pynestml.utils.logger import Logger, LoggingLevel
+from pynestml.utils.messages import Messages
+from pynestml.visitors.ast_higher_order_visitor import ASTHigherOrderVisitor
+from pynestml.visitors.ast_visitor import ASTVisitor
+
+
+def get_variable_block_type(var: ASTVariable) -> Optional[BlockType]:
+    r"""
+    Determines in which block (state, parameters, or internals) a variable was declared.
+
+    :param var: an ASTVariable node
+    :return: the BlockType in which the variable was declared, or None if not found
+    """
+    # Get the symbol from the variable
+    if var.get_scope() is None:
+        return None
+
+    symbol = var.get_scope().resolve_to_symbol(var.get_name(), SymbolKind.VARIABLE)
+
+    if symbol is None or not isinstance(symbol, VariableSymbol):
+        return None
+
+    # Return the block_type from the symbol
+    return symbol.block_type
+
+
+class CoCoParametersInternalsStateDeclarationOrder:
+    r"""
+    Variables in the state block are only allowed to refer back to internals and parameters, internals are only allowed to refer back to parameters, and parameters can refer back to neither but need to be defined in a self-contained manner.
+
+    Allowed:
+
+    .. ::
+
+          parameters:
+              bar real = 123
+
+          state:
+              foo real = bar
+
+    Not allowed:
+
+    .. ::
+
+          parameters:
+              bar real = foo
+
+          state:
+              foo real = 123
+
+    Not allowed:
+
+    .. ::
+
+          state:
+              foo real = 123
+
+          parameters:
+              bar real = foo
+
+    """
+
+    @classmethod
+    def check_co_co(cls, node: ASTModel):
+        """
+        Ensures the coco for the handed over neuron.
+        :param node: a single neuron instance.
+        """
+        node.accept(CoCoParametersDefinedBeforeInternalsDefinedBeforeStateVisitor())
+
+
+class CoCoParametersDefinedBeforeInternalsDefinedBeforeStateVisitor(ASTVisitor):
+
+    def visit_declaration(self, node: ASTDeclaration):
+        """
+        Ensures the coco.
+        :param node: a single equation object.
+        """
+        if not isinstance(node.get_parent(), ASTBlockWithVariables):
+            # this is a local declaration -- skip
+            return
+
+        if node.get_parent().is_parameters:
+            block_containing_declaration = BlockType.PARAMETERS
+        elif node.get_parent().is_internals:
+            block_containing_declaration = BlockType.INTERNALS
+        elif node.get_parent().is_state:
+            block_containing_declaration = BlockType.STATE
+        else:
+            # This should not happen, as declarations should only be allowed in parameters, internals, or state blocks
+            raise Exception(f"Declaration is not contained in parameters, internals, or state block.")
+
+        # get defining expression
+        expr = node.get_expression()
+        assert expr is not None, "Could not find defining expression for declaration"
+
+        # grab variables in ``expr``
+        variables = []
+
+        def collect_vars(node):
+            if isinstance(node, ASTVariable):
+                variables.append(node)
+
+        expr.accept(ASTHigherOrderVisitor(visit_funcs=collect_vars))
+
+        # for each variable, check where it is declared
+        for var in variables:
+            if var.name in PredefinedVariables.get_variables():
+                # skip predefined variables
+                continue
+
+            block_type = get_variable_block_type(var)
+
+            if block_containing_declaration == BlockType.PARAMETERS and block_type == BlockType.STATE:
+                code, message = Messages.get_variable_definition_in_parameters_block_not_allowed_to_refer_to_state_block(str(var))
+                Logger.log_message(code=code, message=message, error_position=node.get_source_position(), log_level=LoggingLevel.ERROR)
+            elif block_containing_declaration == BlockType.PARAMETERS and block_type == BlockType.INTERNALS:
+                code, message = Messages.get_variable_definition_in_parameters_block_not_allowed_to_refer_to_internals_block(str(var))
+                Logger.log_message(code=code, message=message, error_position=node.get_source_position(), log_level=LoggingLevel.ERROR)
+            elif block_containing_declaration == BlockType.INTERNALS and block_type == BlockType.STATE:
+                code, message = Messages.get_variable_definition_in_internals_block_not_allowed_to_refer_to_state_block(str(var))
+                Logger.log_message(code=code, message=message, error_position=node.get_source_position(), log_level=LoggingLevel.ERROR)

pynestml/cocos/co_cos_manager.py

@@ -56,6 +56,7 @@
 from pynestml.cocos.co_co_ode_functions_have_consistent_units import CoCoOdeFunctionsHaveConsistentUnits
 from pynestml.cocos.co_co_output_port_defined_if_emit_call import CoCoOutputPortDefinedIfEmitCall
 from pynestml.cocos.co_co_parameters_assigned_only_in_parameter_block import CoCoParametersAssignedOnlyInParameterBlock
+from pynestml.cocos.co_co_parameters_internals_state_declaration_order import CoCoParametersInternalsStateDeclarationOrder
 from pynestml.cocos.co_co_priorities_correctly_specified import CoCoPrioritiesCorrectlySpecified
 from pynestml.cocos.co_co_resolution_func_legally_used import CoCoResolutionFuncLegallyUsed
 from pynestml.cocos.co_co_resolution_func_used import CoCoResolutionOrStepsFuncUsed
@@ -109,6 +110,14 @@ def check_function_declared_and_correctly_typed(cls, model: ASTModel):
         """
         CoCoFunctionCallsConsistent.check_co_co(model)
 
+    @classmethod
+    def check_parameters_internals_state_declaration_order(cls, model: ASTModel):
+        """
+        Checks that all variables have been declared at most once per scope.
+        :param model: a single model instance
+        """
+        CoCoParametersInternalsStateDeclarationOrder.check_co_co(model)
+
     @classmethod
     def check_variables_unique_in_scope(cls, model: ASTModel):
         """
@@ -436,6 +445,7 @@ def check_cocos(cls, model: ASTModel, after_ast_rewrite: bool = False):
         Logger.set_current_node(model)
 
         cls.check_each_block_defined_at_most_once(model)
+        cls.check_parameters_internals_state_declaration_order(model)
         cls.check_function_defined(model)
         cls.check_variables_unique_in_scope(model)
         cls.check_inline_expression_not_assigned_to(model)

pynestml/symbols/predefined_variables.py

@@ -120,6 +120,15 @@ def get_variable(cls, name):
 
     @classmethod
     def get_variables(cls):
+        """
+        Returns the list of all defined variables.
+        :return: a list of variable symbols.
+        :rtype: list(VariableSymbol)
+        """
+        return cls.name2variable.keys()
+
+    @classmethod
+    def get_variables_to_symbols_dict(cls):
         """
         Returns the list of all defined variables.
         :return: a list of variable symbols.

pynestml/utils/ast_synapse_information_collector.py

@@ -148,7 +148,7 @@ def get_variable_names_of_synapse(self, synapse_inline: ASTInlineExpression, exc
 
     @classmethod
     def get_variable_names_to_ignore(cls):
-        return set(PredefinedVariables.get_variables().keys()).union({"v_comp"})
+        return set(PredefinedVariables.get_variables()).union({"v_comp"})
 
     def get_synapse_specific_internal_declarations(self, synapse_inline: ASTInlineExpression) -> defaultdict:
         synapse_variable_names = self.get_variable_names_of_synapse(

pynestml/utils/ast_utils.py

@@ -637,7 +637,7 @@ def replace_var(_expr=None):
             else:
                 return
 
-            if not var.get_name() in PredefinedVariables.get_variables().keys() \
+            if not var.get_name() in PredefinedVariables.get_variables() \
                and not var.get_name().endswith(suffix):
                 symbol = scope.resolve_to_symbol(var.get_name(), SymbolKind.VARIABLE)
                 if symbol:    # make sure it is not a unit (like "ms")
@@ -662,7 +662,7 @@ def replace_var(_expr=None):
             else:
                 return
 
-            if not var.get_name() in PredefinedVariables.get_variables().keys():
+            if not var.get_name() in PredefinedVariables.get_variables():
                 symbol = new_scope.resolve_to_symbol(var.get_name(), SymbolKind.VARIABLE)
                 if symbol:    # make sure it is not a unit (like "ms")
                     var.update_scope(new_scope)

pynestml/utils/messages.py

@@ -148,6 +148,9 @@ class MessageCode(Enum):
     WEIGHT_VARIABLE_NOT_SPECIFIED = 119
     DELAY_VARIABLE_NOT_FOUND = 120
     WEIGHT_VARIABLE_NOT_FOUND = 121
+    VARIABLE_DEFINITION_IN_PARAMETERS_BLOCK_NOT_ALLOWED_TO_REFER_TO_STATE_BLOCK = 122
+    VARIABLE_DEFINITION_IN_PARAMETERS_BLOCK_NOT_ALLOWED_TO_REFER_TO_INTERNALS_BLOCK = 123
+    VARIABLE_DEFINITION_IN_INTERNALS_BLOCK_NOT_ALLOWED_TO_REFER_TO_STATE_BLOCK = 124
 
 
 class Messages:
@@ -1181,3 +1184,21 @@ def get_weight_variable_not_found(cls, variable_name: str) -> Tuple[MessageCode,
         message = "Weight variable \"" + variable_name + "\" not found in synapse. Please see https://nestml.readthedocs.io/en/latest/running/running_nest.html#dendritic-delay-and-synaptic-weight"
 
         return MessageCode.WEIGHT_VARIABLE_NOT_FOUND, message
+
+    @classmethod
+    def get_variable_definition_in_parameters_block_not_allowed_to_refer_to_state_block(cls, variable_name: str) -> Tuple[MessageCode, str]:
+        message = "Definition of variable \"" + variable_name + "\" in ``parameters`` block is not allowed to refer to variables defined in the ``state`` block."
+
+        return MessageCode.VARIABLE_DEFINITION_IN_PARAMETERS_BLOCK_NOT_ALLOWED_TO_REFER_TO_STATE_BLOCK, message
+
+    @classmethod
+    def get_variable_definition_in_parameters_block_not_allowed_to_refer_to_internals_block(cls, variable_name: str) -> Tuple[MessageCode, str]:
+        message = "Definition of variable \"" + variable_name + "\" in ``parameters`` block is not allowed to refer to variables defined in the ``internals`` block."
+
+        return MessageCode.VARIABLE_DEFINITION_IN_PARAMETERS_BLOCK_NOT_ALLOWED_TO_REFER_TO_INTERNALS_BLOCK, message
+
+    @classmethod
+    def get_variable_definition_in_internals_block_not_allowed_to_refer_to_state_block(cls, variable_name: str) -> Tuple[MessageCode, str]:
+        message = "Definition of variable \"" + variable_name + "\" in ``internals`` block is not allowed to refer to variables defined in the ``state`` block."
+
+        return MessageCode.VARIABLE_DEFINITION_IN_INTERNALS_BLOCK_NOT_ALLOWED_TO_REFER_TO_STATE_BLOCK, message

pynestml/visitors/ast_symbol_table_visitor.py

@@ -63,7 +63,7 @@ def visit_model(self, node: ASTModel) -> None:
         node.update_scope(scope)
         node.get_body().update_scope(scope)
         # now first, we add all predefined elements to the scope
-        variables = PredefinedVariables.get_variables()
+        variables = PredefinedVariables.get_variables_to_symbols_dict()
         functions = PredefinedFunctions.get_function_symbols()
         types = PredefinedTypes.get_types()
         for symbol in variables.keys():