Add a Transformer to support multiple backends with different preferred prefixes a.k.a Non-dimensionalisation-transformer

doc/running/index.rst

@@ -4,8 +4,18 @@ Running NESTML
 Running NESTML causes several processing steps to occur:
 
 1. The model is parsed from file and checked (syntax, consistent physical units, and so on).
-2. Code is generated from the model by one of the "code generators" selected when NESTML was invoked.
-3. If necessary, the code is compiled and built by the "builder" that belongs to the selected code generator.
+2. A series of transformation steps are optionally carried out on the model. Which set of transformers is used is automatically determined based on the target platform selected when NESTML was invoked. (See the section "Available transformers" below.)
+3. Code is generated from the model by one of the "code generators" selected when NESTML was invoked. If necessary, the code is compiled and built by the "builder" that belongs to the selected code generator. (See the section "Supported target platforms" below.)
+
+
+Available transformers
+----------------------
+
+- **Non-dimensionalisation transformer**
+
+  :class:`pynestml.transformers.non_dimensionalisation_transformer.NonDimensionalisationTransformer`
+
+  Eliminates physical units from models, converting numerical values to real numbers.
 
 
 Supported target platforms

models/neurons/aeif_cond_exp_neuron.nestml

@@ -1,140 +1,7 @@
-# aeif_cond_exp - Conductance based exponential integrate-and-fire neuron model
-# #############################################################################
-# 
-# Description
-# +++++++++++
-# 
-# aeif_cond_exp is the adaptive exponential integrate and fire neuron
-# according to Brette and Gerstner (2005), with post-synaptic
-# conductances in the form of truncated exponentials.
-# 
-# The membrane potential is given by the following differential equation:
-# 
-# .. math::
-# 
-#    C_m \frac{dV_m}{dt} =
-#    -g_L(V_m-E_L)+g_L\Delta_T\exp\left(\frac{V_m-V_{th}}{\Delta_T}\right) - g_e(t)(V_m-E_e) \\
-#    -g_i(t)(V_m-E_i)-w +I_e
-# 
-# and
-# 
-# .. math::
-# 
-#    \tau_w \frac{dw}{dt} = a(V_m-E_L) - w
-# 
-# Note that the membrane potential can diverge to positive infinity due to the exponential term. To avoid numerical instabilities, instead of :math:`V_m`, the value :math:`\min(V_m,V_{peak})` is used in the dynamical equations.
-# 
-# .. note::
-# 
-#    The default refractory period for ``aeif`` models is zero, consistent with the model definition in
-#    Brette & Gerstner [1]_.  Thus, an ``aeif`` neuron with default parameters can fire multiple spikes in a single
-#    time step, which can lead to exploding spike numbers and extreme slow-down of simulations.
-#    To avoid such unphysiological behavior, you should set a refractory time ``refr_t > 0``.
-# 
-# 
-# References
-# ++++++++++
-# 
-# .. [1] Brette R and Gerstner W (2005). Adaptive exponential
-#        integrate-and-fire model as an effective description of neuronal
-#        activity. Journal of Neurophysiology. 943637-3642
-#        DOI: https://doi.org/10.1152/jn.00686.2005
-# 
-# 
-# See also
-# ++++++++
-# 
-# iaf_cond_exp, aeif_cond_alpha
-#
-#
-# Copyright statement
-# +++++++++++++++++++
-#
-# 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/>.
-#
-#
 model aeif_cond_exp_neuron:
-
     state:
-        V_m mV = E_L   # Membrane potential
-        w pA = 0 pA    # Spike-adaptation current
-        refr_t ms = 0 ms    # Refractory period timer
-
-    equations:
-        inline V_bounded mV = min(V_m, V_peak)    # prevent exponential divergence
-        kernel g_inh = exp(-t / tau_syn_inh)
-        kernel g_exc = exp(-t / tau_syn_exc)
-
-        # Add inlines to simplify the equation definition of V_m
-        inline exp_arg real = (V_bounded - V_th) / Delta_T
-        inline I_spike pA = g_L * Delta_T * exp(exp_arg)
-        inline I_syn_exc pA = convolve(g_exc, exc_spikes) * nS * (V_bounded - E_exc)
-        inline I_syn_inh pA = convolve(g_inh, inh_spikes) * nS * (V_bounded - E_inh)
-
-        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)
+        V_m V = E_L   # Membrane potential
 
     parameters:
-        # membrane parameters
-        C_m pF = 281.0 pF       # Membrane Capacitance
-        refr_T ms = 2 ms        # Duration of refractory period
-        V_reset mV = -60.0 mV   # Reset Potential
-        g_L nS = 30.0 nS        # Leak Conductance
-        E_L mV = -70.6 mV       # Leak reversal Potential (aka resting potential)
-
-        # spike adaptation parameters
-        a nS = 4 nS             # Subthreshold adaptation
-        b pA = 80.5 pA          # Spike-triggered adaptation
-        Delta_T mV = 2.0 mV     # Slope factor
-        tau_w ms = 144.0 ms     # Adaptation time constant
-        V_th mV = -50.4 mV      # Spike initiation threshold
-        V_peak mV = 0 mV        # Spike detection threshold
-
-        # synaptic parameters
-        E_exc mV = 0 mV            # Excitatory reversal Potential
-        tau_syn_exc ms = 0.2 ms    # Synaptic Time Constant Excitatory Synapse
-        E_inh mV = -85.0 mV        # Inhibitory reversal Potential
-        tau_syn_inh ms = 2.0 ms    # Synaptic Time Constant for Inhibitory Synapse
-
-        # constant external input current
-        I_e pA = 0 pA
-
-    input:
-        inh_spikes <- inhibitory spike
-        exc_spikes <- excitatory spike
-        I_stim pA <- continuous
-
-    output:
-        spike
-
-    update:
-        if refr_t > 0 ms:
-            # neuron is absolute refractory, do not evolve V_m
-            integrate_odes(w, refr_t)
-        else:
-            # neuron not refractory
-            integrate_odes(w, V_m)
-
-    onCondition(refr_t <= 0 ms and V_m >= V_peak):
-        # threshold crossing
-        refr_t = refr_T    # start of the refractory period
-        V_m = V_reset # clamp potential
-        w += b
-        emit_spike()
+        E_L V = -70.6E3 uV       # Leak reversal Potential (aka resting potential)
+        foo V = 1000 * V_m       # Leak reversal Potential (aka resting potential)

pynestml/cocos/co_co_function_calls_consistent.py

@@ -28,7 +28,6 @@
 from pynestml.utils.ast_utils import ASTUtils
 from pynestml.utils.logger import Logger, LoggingLevel
 from pynestml.utils.messages import Messages
-from pynestml.utils.type_caster import TypeCaster
 from pynestml.visitors.ast_visitor import ASTVisitor
 
 
@@ -109,5 +108,5 @@ def visit_function_call(self, node):
                 # variadic type symbol accepts anything
                 return
 
-            if not actual_type.equals(expected_type) and not isinstance(expected_type, TemplateTypeSymbol):
-                TypeCaster.try_to_recover_or_error(expected_type, actual_type, actual_arg)
+            # if not actual_type.equals(expected_type) and not isinstance(expected_type, TemplateTypeSymbol):
+            #     TypeCaster.try_to_recover_or_error(expected_type, actual_type, actual_arg)

pynestml/cocos/co_co_illegal_expression.py

@@ -28,7 +28,6 @@
 from pynestml.utils.logger import LoggingLevel, Logger
 from pynestml.utils.logging_helper import LoggingHelper
 from pynestml.utils.messages import Messages
-from pynestml.utils.type_caster import TypeCaster
 from pynestml.visitors.ast_visitor import ASTVisitor
 
 
@@ -68,9 +67,8 @@ def visit_declaration(self, node):
             if isinstance(rhs_type, ErrorTypeSymbol):
                 LoggingHelper.drop_missing_type_error(node)
                 return
-            if self.__types_do_not_match(lhs_type, rhs_type):
-                TypeCaster.try_to_recover_or_error(lhs_type, rhs_type, node.get_expression())
-        return
+            # if self.__types_do_not_match(lhs_type, rhs_type):
+            #     TypeCaster.try_to_recover_or_error(lhs_type, rhs_type, node.get_expression())
 
     def visit_inline_expression(self, node):
         """
@@ -82,8 +80,8 @@ def visit_inline_expression(self, node):
         if isinstance(rhs_type, ErrorTypeSymbol):
             LoggingHelper.drop_missing_type_error(node)
             return
-        if self.__types_do_not_match(lhs_type, rhs_type):
-            TypeCaster.try_to_recover_or_error(lhs_type, rhs_type, node.get_expression())
+        # if self.__types_do_not_match(lhs_type, rhs_type):
+        #     TypeCaster.try_to_recover_or_error(lhs_type, rhs_type, node.get_expression())
 
     def visit_assignment(self, node):
         """
@@ -120,23 +118,23 @@ def handle_compound_assignment(self, node):
         lhs_type_symbol = lhs_variable_symbol.get_type_symbol()
 
         if node.is_compound_product:
-            if self.__types_do_not_match(lhs_type_symbol, lhs_type_symbol * rhs_type_symbol):
-                TypeCaster.try_to_recover_or_error(lhs_type_symbol, lhs_type_symbol * rhs_type_symbol,
-                                                   node.get_expression())
-                return
+            # if self.__types_do_not_match(lhs_type_symbol, lhs_type_symbol * rhs_type_symbol):
+            #     TypeCaster.try_to_recover_or_error(lhs_type_symbol, lhs_type_symbol * rhs_type_symbol,
+            #                                        node.get_expression())
+            #     return
             return
 
         if node.is_compound_quotient:
-            if self.__types_do_not_match(lhs_type_symbol, lhs_type_symbol / rhs_type_symbol):
-                TypeCaster.try_to_recover_or_error(lhs_type_symbol, lhs_type_symbol / rhs_type_symbol,
-                                                   node.get_expression())
-                return
+            # if self.__types_do_not_match(lhs_type_symbol, lhs_type_symbol / rhs_type_symbol):
+            #     TypeCaster.try_to_recover_or_error(lhs_type_symbol, lhs_type_symbol / rhs_type_symbol,
+            #                                        node.get_expression())
+            #     return
             return
 
         assert node.is_compound_sum or node.is_compound_minus
-        if self.__types_do_not_match(lhs_type_symbol, rhs_type_symbol):
-            TypeCaster.try_to_recover_or_error(lhs_type_symbol, rhs_type_symbol,
-                                               node.get_expression())
+        # if self.__types_do_not_match(lhs_type_symbol, rhs_type_symbol):
+        #     TypeCaster.try_to_recover_or_error(lhs_type_symbol, rhs_type_symbol,
+        #                                        node.get_expression())
 
     @staticmethod
     def __types_do_not_match(lhs_type_symbol, rhs_type_symbol):
@@ -154,11 +152,10 @@ def handle_simple_assignment(self, node):
             LoggingHelper.drop_missing_type_error(node)
             return
 
-        if lhs_variable_symbol is not None and self.__types_do_not_match(lhs_variable_symbol.get_type_symbol(),
-                                                                         rhs_type_symbol):
-            TypeCaster.try_to_recover_or_error(lhs_variable_symbol.get_type_symbol(), rhs_type_symbol,
-                                               node.get_expression())
-        return
+        # if lhs_variable_symbol is not None and self.__types_do_not_match(lhs_variable_symbol.get_type_symbol(),
+        #                                                                  rhs_type_symbol):
+        #     TypeCaster.try_to_recover_or_error(lhs_variable_symbol.get_type_symbol(), rhs_type_symbol,
+        #                                        node.get_expression())
 
     def visit_if_clause(self, node):
         """

pynestml/cocos/co_co_output_port_defined_if_emit_call.py

@@ -104,6 +104,7 @@ def visit_function_call(self, node: ASTFunctionCall):
                     # types are not equal, but castable
                     code, message = Messages.get_implicit_cast_rhs_to_lhs(output_block_attr_type_sym.print_symbol(),
                                                                           emit_spike_arg_type_sym.print_symbol())
+                    # XXX: this should be removed
                     Logger.log_message(error_position=output_blocks[0].get_source_position(),
                                        code=code, message=message, log_level=LoggingLevel.WARNING)
                     continue

pynestml/cocos/co_co_user_defined_function_correctly_defined.py

@@ -28,7 +28,6 @@
 from pynestml.symbols.symbol import SymbolKind
 from pynestml.utils.logger import LoggingLevel, Logger
 from pynestml.utils.messages import Messages
-from pynestml.utils.type_caster import TypeCaster
 
 
 class CoCoUserDefinedFunctionCorrectlyDefined(CoCo):
@@ -128,9 +127,9 @@ def __check_return_recursively(cls, type_symbol=None, stmts=None, ret_defined=Fa
                         code, message = Messages.get_type_could_not_be_derived(cls.processed_function.get_name())
                         Logger.log_message(error_position=stmt.get_source_position(),
                                            code=code, message=message, log_level=LoggingLevel.ERROR)
-                    elif not type_of_return.equals(type_symbol):
-                        TypeCaster.try_to_recover_or_error(type_symbol, type_of_return,
-                                                           stmt.get_return_stmt().get_expression())
+                    # elif not type_of_return.equals(type_symbol):
+                    #     TypeCaster.try_to_recover_or_error(type_symbol, type_of_return,
+                    #                                        stmt.get_return_stmt().get_expression())
             elif isinstance(stmt, ASTCompoundStmt):
                 # otherwise it is a compound stmt, thus check recursively
                 if stmt.is_if_stmt():

pynestml/codegeneration/nest_code_generator.py

@@ -18,7 +18,6 @@
 #
 # 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 Any, Dict, List, Mapping, Optional, Sequence, Tuple
 
 import datetime
@@ -62,6 +61,7 @@
 from pynestml.symbols.unit_type_symbol import UnitTypeSymbol
 from pynestml.symbols.symbol import SymbolKind
 from pynestml.transformers.inline_expression_expansion_transformer import InlineExpressionExpansionTransformer
+from pynestml.transformers.non_dimensionalisation_transformer import NonDimensionalisationTransformer
 from pynestml.transformers.synapse_post_neuron_transformer import SynapsePostNeuronTransformer
 from pynestml.utils.ast_utils import ASTUtils
 from pynestml.utils.logger import Logger
@@ -413,6 +413,24 @@ def analyse_neuron(self, neuron: ASTModel) -> Tuple[Dict[str, ASTAssignment], Di
 
         spike_updates, post_spike_updates = self.get_spike_update_expressions(neuron, kernel_buffers, [analytic_solver, numeric_solver], delta_factors)
 
+        # update expressions are processed by the transformer
+        for var_name, spike_updates_for_var in spike_updates.items():
+            transformed_spike_update_exprs = []
+            for spike_update_expr in spike_updates_for_var:
+                transformed_spike_update_exprs.append(NonDimensionalisationTransformer({"quantity_to_preferred_prefix": {"electrical potential": "m",
+                                                         "electrical current": "p",
+                                                         "electrical capacitance": "p",
+                                                         "electrical resistance": "M",
+                                                         "electrical conductance": "n",
+                                                         "time": "m"}}).transform([spike_update_expr])[0])
+
+            spike_updates[var_name] = transformed_spike_update_exprs
+
+
+
+
+        # post_spike_updates = NonDimensionalisationTransformer().transform(post_spike_updates)
+
         return spike_updates, post_spike_updates, equations_with_delay_vars, equations_with_vector_vars
 
     def analyse_synapse(self, synapse: ASTModel) -> Dict[str, ASTAssignment]: