Merge pull request #1006 from clinssen/inhomogeneous_poisson

Add inhomogeneous Poisson generator NESTML model

doc/nestml_language/nestml_language_concepts.rst

@@ -517,6 +517,9 @@ The following functions are predefined in NESTML and can be used out of the box.
    * - ``random_normal``
      - mean, std
      - Returns a sample from a normal (Gaussian) distribution with parameters "mean" and "standard deviation"
+   * - ``random_poisson``
+     - rate
+     - Returns a sample from a Poissonian distribution with rate parameter (expected value) "rate".
    * - ``random_uniform``
      - offset, scale
      - Returns a sample from a uniform distribution in the interval [offset, offset + scale)

doc/tutorials/inhomogeneous_poisson/inhomogeneous_poisson.nestml

@@ -0,0 +1,51 @@
+"""
+inhomogeneous_poisson - Inhomogeneous Poisson process model
+###########################################################
+
+Description
++++++++++++
+
+Inhomogeneous Poisson process model.
+
+The rate of the model is piecewise constant and is defined by an array containing desired rates (in units of 1/s) and an array of equal length containing the corresponding times (in units of ms). Please see the documentation for the NEST built-in inhomogeneous_poisson_generator for more details [2].
+
+
+See also
+++++++++
+
+See the inhomogeneous Poisson generator NESTML tutorial for a usage example.
+
+
+References
+++++++++++
+
+[1] Wikipedia contributors. 'Poisson Point Process.' Wikipedia, The Free Encyclopedia. Accessed on February 23, 2024. https://en.wikipedia.org/wiki/Poisson_point_process.
+
+[2] https://nest-simulator.readthedocs.io/en/stable/models/inhomogeneous_poisson_generator.html
+"""
+model inhomogeneous_poisson_neuron:
+
+    state:
+        idx integer = 0
+
+    parameters:
+        N integer = 1
+        rate_times[N] ms = 0 ms
+        rate_values[N] ms**-1 = 0 s**-1
+
+    output:
+        spike
+
+    update:
+        # update the rate
+        if N > 0:
+            while idx < N - 1 and t >= rate_times[idx + 1]:
+                idx += 1
+
+        rate ms**-1 = rate_values[idx]
+
+        n_spikes_in_this_timestep integer = random_poisson(rate * resolution() * 1E-3)
+        if n_spikes_in_this_timestep > 0:
+            emit_spike()
+            # XXX: TODO: support multiplicity in spikes (see https://github.com/nest/nestml/issues/946)
+

doc/tutorials/tutorials_list.rst

@@ -17,6 +17,10 @@ Creating neuron models
 
   Implement the Ornstein-Uhlenbeck process in NESTML and use it to inject a noise current into a neuron.
 
+* :doc:`Inhomogeneous Poisson generator </tutorials/inhomogeneous_poisson/inhomogeneous_poisson>`
+
+  Create a model that emits spikes according to an inhomogeneous Poisson distribution.
+
 
 Creating synapse models
 -----------------------

pynestml/cocos/co_co_vector_declaration_right_size.py

@@ -53,45 +53,42 @@ def visit_variable(self, node: ASTVariable):
             if isinstance(self._neuron.get_parent(node), ASTDeclaration):
                 # node is being declared: size should be >= 1
                 min_index = 1
-            else:
-                # node is being indexed; index should be >= 0
-                min_index = 0
 
-            # account for negative numbers as an index, which is an expression (e.g. "-42")
-            if isinstance(vector_parameter, ASTExpression) and vector_parameter.is_unary_operator() and vector_parameter.get_unary_operator().is_unary_minus and vector_parameter.get_expression().is_numeric_literal():
-                code, message = Messages.get_vector_parameter_wrong_size(node.get_complete_name() + "[" + str(vector_parameter) + "]", str(vector_parameter))
-                Logger.log_message(error_position=node.get_source_position(), log_level=LoggingLevel.ERROR,
-                                   code=code, message=message)
-                return
-
-            # otherwise, ``vector_parameter`` is a SimpleExpression
-            assert vector_parameter.is_variable() or vector_parameter.is_numeric_literal()
-
-            vector_parameter_val = None
-
-            if vector_parameter.is_numeric_literal():
-                if not isinstance(vector_parameter.type, IntegerTypeSymbol):
-                    code, message = Messages.get_vector_parameter_wrong_type(node.get_complete_name() + "[" + str(vector_parameter) + "]")
+                # account for negative numbers as an index, which is an expression (e.g. "-42")
+                if isinstance(vector_parameter, ASTExpression) and vector_parameter.is_unary_operator() and vector_parameter.get_unary_operator().is_unary_minus and vector_parameter.get_expression().is_numeric_literal():
+                    code, message = Messages.get_vector_parameter_wrong_size(node.get_complete_name() + "[" + str(vector_parameter) + "]", str(vector_parameter))
                     Logger.log_message(error_position=node.get_source_position(), log_level=LoggingLevel.ERROR,
                                        code=code, message=message)
                     return
 
-                vector_parameter_val = int(vector_parameter.get_numeric_literal())
+                # otherwise, ``vector_parameter`` is a SimpleExpression
+                assert vector_parameter.is_variable() or vector_parameter.is_numeric_literal()
 
-            if vector_parameter.is_variable():
-                symbol = vector_parameter.get_scope().resolve_to_symbol(vector_parameter.get_variable().get_complete_name(),
-                                                                        SymbolKind.VARIABLE)
+                vector_parameter_val = None
 
-                if symbol is None or not isinstance(symbol.get_type_symbol(), IntegerTypeSymbol):
-                    code, message = Messages.get_vector_parameter_wrong_type(node.get_complete_name() + "[" + str(vector_parameter) + "]")
-                    Logger.log_message(error_position=node.get_source_position(), log_level=LoggingLevel.ERROR,
-                                       code=code, message=message)
-                    return
+                if vector_parameter.is_numeric_literal():
+                    if not isinstance(vector_parameter.type, IntegerTypeSymbol):
+                        code, message = Messages.get_vector_parameter_wrong_type(node.get_complete_name() + "[" + str(vector_parameter) + "]")
+                        Logger.log_message(error_position=node.get_source_position(), log_level=LoggingLevel.ERROR,
+                                           code=code, message=message)
+                        return
 
-                vector_parameter_val = int(str(symbol.get_declaring_expression()))
+                    vector_parameter_val = int(vector_parameter.get_numeric_literal())
 
-            if vector_parameter_val is not None and vector_parameter_val < min_index:
-                code, message = Messages.get_vector_parameter_wrong_size(node.get_complete_name(),
-                                                                         str(vector_parameter_val))
-                Logger.log_message(error_position=node.get_source_position(), log_level=LoggingLevel.ERROR,
-                                   code=code, message=message)
+                if vector_parameter.is_variable():
+                    symbol = vector_parameter.get_scope().resolve_to_symbol(vector_parameter.get_variable().get_complete_name(),
+                                                                            SymbolKind.VARIABLE)
+
+                    if symbol is None or not isinstance(symbol.get_type_symbol(), IntegerTypeSymbol):
+                        code, message = Messages.get_vector_parameter_wrong_type(node.get_complete_name() + "[" + str(vector_parameter) + "]")
+                        Logger.log_message(error_position=node.get_source_position(), log_level=LoggingLevel.ERROR,
+                                           code=code, message=message)
+                        return
+
+                    vector_parameter_val = int(str(symbol.get_declaring_expression()))
+
+                if vector_parameter_val is not None and vector_parameter_val < min_index:
+                    code, message = Messages.get_vector_parameter_wrong_size(node.get_complete_name(),
+                                                                             str(vector_parameter_val))
+                    Logger.log_message(error_position=node.get_source_position(), log_level=LoggingLevel.ERROR,
+                                       code=code, message=message)

pynestml/cocos/co_co_vector_variable_in_non_vector_declaration.py

@@ -63,7 +63,7 @@ def visit_declaration(self, node):
             for variable in variables:
                 if variable is not None:
                     symbol = node.get_scope().resolve_to_symbol(variable.get_complete_name(), SymbolKind.VARIABLE)
-                    if symbol is not None and symbol.has_vector_parameter() and not node.has_size_parameter():
+                    if symbol is not None and symbol.has_vector_parameter() and not variable.has_vector_parameter():
                         code, message = Messages.get_vector_in_non_vector(vector=symbol.get_symbol_name(),
                                                                           non_vector=list(var.get_complete_name() for
                                                                                           var in

pynestml/codegeneration/printers/latex_function_call_printer.py

@@ -62,10 +62,13 @@ def _print_function_call(self, node: ASTFunctionCall) -> str:
             return r'\text{steps}'
 
         if function_name == PredefinedFunctions.RANDOM_NORMAL:
-            return r'\mathcal{N}({!s}, {!s}'
+            return r'\mathcal{N}({!s}, {!s})'
+
+        if function_name == PredefinedFunctions.RANDOM_POISSON:
+            return r'\mathrm{Poisson}({!s})'
 
         if function_name == PredefinedFunctions.RANDOM_UNIFORM:
-            return r'\mathcal{U}({!s}, {!s}'
+            return r'\mathcal{U}({!s}, {!s})'
 
         if function_name == PredefinedFunctions.EMIT_SPIKE:
             return r'\text{spike}'

pynestml/codegeneration/printers/nest2_cpp_function_call_printer.py

@@ -47,6 +47,9 @@ def _print_function_call_format_string(self, function_call: ASTFunctionCall) ->
         if function_name == PredefinedFunctions.RANDOM_NORMAL:
             return '(({!s}) + ({!s}) * ' + 'normal_dev_( nest::kernel().rng_manager.get_rng( ' + 'get_thread() ) ))'
 
+        if function_name == PredefinedFunctions.RANDOM_POISSON:
+            return '(poisson_dev_.set_lambda({!s}), poisson_dev_.ldev( nest::get_vp_specific_rng( ' + 'get_thread() ) ))'
+
         if function_name == PredefinedFunctions.RANDOM_UNIFORM:
             return '(({!s}) + ({!s}) * nest::kernel().rng_manager.get_rng( ' + 'get_thread() )->drand())'
 

pynestml/codegeneration/printers/nest_cpp_function_call_printer.py

@@ -55,6 +55,9 @@ def _print_function_call_format_string(self, function_call: ASTFunctionCall) ->
         if function_name == PredefinedFunctions.RANDOM_NORMAL:
             return '(({!s}) + ({!s}) * ' + 'normal_dev_( nest::get_vp_specific_rng( ' + 'get_thread() ) ))'
 
+        if function_name == PredefinedFunctions.RANDOM_POISSON:
+            return '([&]() -> int {{ nest::poisson_distribution::param_type poisson_params({!s}); int sample = poisson_dev_( nest::get_vp_specific_rng( get_thread() ), poisson_params); return sample; }})()'   # double curly braces {{ }} due to passing through str.format() later
+
         if function_name == PredefinedFunctions.RANDOM_UNIFORM:
             return '(({!s}) + ({!s}) * nest::get_vp_specific_rng( ' + 'get_thread() )->drand())'
 

pynestml/codegeneration/printers/nest_variable_printer.py

@@ -85,7 +85,7 @@ def print_variable(self, variable: ASTVariable) -> str:
 
         vector_param = ""
         if self.with_vector_parameter and symbol.has_vector_parameter():
-            vector_param = "[" + self._print_vector_parameter_name_reference(variable) + "]"
+            vector_param = "[" + self._expression_printer.print(variable.get_vector_parameter()) + "]"
 
         if symbol.is_buffer():
             if isinstance(symbol.get_type_symbol(), UnitTypeSymbol):
@@ -127,24 +127,6 @@ def _print_delay_variable(self, variable: ASTVariable) -> str:
 
         return ""
 
-    def _print_vector_parameter_name_reference(self, variable: ASTVariable) -> str:
-        """
-        Converts the vector parameter into NEST processable format
-        :param variable:
-        :return:
-        """
-        vector_parameter = variable.get_vector_parameter()
-        vector_parameter_var = vector_parameter.get_variable()
-        if vector_parameter_var:
-            vector_parameter_var.scope = variable.get_scope()
-
-            symbol = vector_parameter_var.get_scope().resolve_to_symbol(vector_parameter_var.get_complete_name(),
-                                                                        SymbolKind.VARIABLE)
-            if symbol and symbol.block_type == BlockType.LOCAL:
-                return symbol.get_symbol_name()
-
-        return self._expression_printer.print(vector_parameter)
-
     def _print_buffer_value(self, variable: ASTVariable) -> str:
         """
         Converts for a handed over symbol the corresponding name of the buffer to a nest processable format.

pynestml/codegeneration/printers/python_variable_printer.py

@@ -87,7 +87,7 @@ def print_variable(self, variable: ASTVariable) -> str:
 
         vector_param = ""
         if self.with_vector_parameter and symbol.has_vector_parameter():
-            vector_param = "[" + self._print_vector_parameter_name_reference(variable) + "]"
+            vector_param = "[" + self._expression_printer.print(variable.get_vector_parameter()) + "]"
 
         if symbol.is_buffer():
             if isinstance(symbol.get_type_symbol(), UnitTypeSymbol):
@@ -126,24 +126,6 @@ def _print_delay_variable(self, variable: ASTVariable):
                 return "get_delayed_" + variable.get_name() + "()"
         return ""
 
-    def _print_vector_parameter_name_reference(self, variable: ASTVariable) -> str:
-        """
-        Converts the vector parameter into Python format
-        :param variable:
-        :return:
-        """
-        vector_parameter = variable.get_vector_parameter()
-        vector_parameter_var = vector_parameter.get_variable()
-        if vector_parameter_var:
-            vector_parameter_var.scope = variable.get_scope()
-
-            symbol = vector_parameter_var.get_scope().resolve_to_symbol(vector_parameter_var.get_complete_name(),
-                                                                        SymbolKind.VARIABLE)
-            if symbol and symbol.block_type == BlockType.LOCAL:
-                return symbol.get_symbol_name()
-
-        return self._expression_printer.print(vector_parameter)
-
     def _print(self, variable, symbol, with_origin: bool = True) -> str:
         variable_name = PythonVariablePrinter._print_python_name(variable.get_complete_name())
 

pynestml/codegeneration/printers/spinnaker_c_variable_printer.py

@@ -76,7 +76,7 @@ def print_variable(self, variable: ASTVariable) -> str:
 
         vector_param = ""
         if self.with_vector_parameter and symbol.has_vector_parameter():
-            vector_param = "[" + self._print_vector_parameter_name_reference(variable) + "]"
+            vector_param = "[" + self._expression_printer.print(variable.get_vector_parameter()) + "]"
 
         if symbol.is_buffer():
             if isinstance(symbol.get_type_symbol(), UnitTypeSymbol):
@@ -114,24 +114,6 @@ def _print_delay_variable(self, variable: ASTVariable) -> str:
 
         return ""
 
-    def _print_vector_parameter_name_reference(self, variable: ASTVariable) -> str:
-        """
-        Converts the vector parameter into SPINNAKER processable format
-        :param variable:
-        :return:
-        """
-        vector_parameter = variable.get_vector_parameter()
-        vector_parameter_var = vector_parameter.get_variable()
-        if vector_parameter_var:
-            vector_parameter_var.scope = variable.get_scope()
-
-            symbol = vector_parameter_var.get_scope().resolve_to_symbol(vector_parameter_var.get_complete_name(),
-                                                                        SymbolKind.VARIABLE)
-            if symbol and symbol.block_type == BlockType.LOCAL:
-                return symbol.get_symbol_name()
-
-        return self._expression_printer.print(vector_parameter)
-
     def _print_buffer_value(self, variable: ASTVariable) -> str:
         """
         Converts for a handed over symbol the corresponding name of the buffer to a SPINNAKER processable format.

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

@@ -63,6 +63,7 @@ along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 // Includes for random number generator
 {%- if nest_version.startswith("v2") %}
 #include "normal_randomdev.h"
+#include "poisson_randomdev.h"
 #include "uniform_randomdev.h"
 {%- else %}
 #include <random>
@@ -873,8 +874,10 @@ private:
 
 {%- if nest_version.startswith("v2") %}
   librandom::NormalRandomDev normal_dev_; //!< random deviate generator
+  librandom::PoissonRandomDev poisson_dev_; //!< random deviate generator
 {%- else %}
   nest::normal_distribution normal_dev_; //!< random deviate generator
+  nest::poisson_distribution poisson_dev_; //!< random deviate generator
 {%- endif %}
 {%- endif %}
 {%- if has_delay_variables %}

pynestml/symbols/predefined_functions.py

@@ -45,6 +45,7 @@ class PredefinedFunctions:
         LOGGER_INFO           The callee name of the logger-info function.
         LOGGER_WARNING        The callee name of the logger-warning function.
         RANDOM_NORMAL         The callee name of the function used to generate a random normal (Gaussian) distributed variable with parameters `mean` and `var` (variance).
+        RANDOM_POISSON        The callee name of the function used to generate a random Poissonian distributed variable with the rate parameter `lmbda` (expected value).
         RANDOM_UNIFORM        The callee name of the function used to generate a random sample from a uniform distribution in the interval `[offset, offset + scale)`.
         EXPM1                 The callee name of the exponent (alternative) function.
         DELTA                 The callee name of the delta function.
@@ -76,6 +77,7 @@ class PredefinedFunctions:
     LOGGER_INFO = "info"
     LOGGER_WARNING = "warning"
     RANDOM_NORMAL = "random_normal"
+    RANDOM_POISSON = "random_poisson"
     RANDOM_UNIFORM = "random_uniform"
     EXPM1 = "expm1"
     CLIP = "clip"
@@ -113,6 +115,7 @@ def register_functions(cls):
         cls.__register_logger_info_function()
         cls.__register_logger_warning_function()
         cls.__register_random_normal_function()
+        cls.__register_random_poisson_function()
         cls.__register_random_uniform_function()
         cls.__register_exp1_function()
         cls.__register_delta_function()
@@ -310,6 +313,16 @@ def __register_random_normal_function(cls):
                                 element_reference=None, is_predefined=True)
         cls.name2function[cls.RANDOM_NORMAL] = symbol
 
+    @classmethod
+    def __register_random_poisson_function(cls):
+        """
+        Registers the random method as used to generate a random Poissonian distributed variable with the rate parameter `lmbda` (expected value).
+        """
+        symbol = FunctionSymbol(name=cls.RANDOM_POISSON, param_types=[PredefinedTypes.get_real_type()],
+                                return_type=PredefinedTypes.get_integer_type(),
+                                element_reference=None, is_predefined=True)
+        cls.name2function[cls.RANDOM_POISSON] = symbol
+
     @classmethod
     def __register_random_uniform_function(cls):
         """

pynestml/visitors/ast_random_number_generator_visitor.py

@@ -59,6 +59,6 @@ def visit_simple_expression(self, node):
             node.type = ErrorTypeSymbol()
             return
 
-        if function_name == PredefinedFunctions.RANDOM_NORMAL:
+        if function_name in [PredefinedFunctions.RANDOM_NORMAL, PredefinedFunctions.RANDOM_POISSON, PredefinedFunctions.RANDOM_UNIFORM]:
             self._norm_rng_is_used = True
             return

pynestml/visitors/ast_set_vector_parameter_in_update_expressions.py

@@ -38,3 +38,5 @@ def visit_variable(self, node: ASTVariable):
         """
         if node.get_name() == self.var.get_name():
             node.set_vector_parameter(self.var.get_vector_parameter())
+            if self.var.get_vector_parameter():
+                self.var.get_vector_parameter().update_scope(node.get_scope())