fix convolutions integration and enhance integrate_odes() to integrat…

…e specific ODEs

pynestml/codegeneration/nest_code_generator.py

@@ -302,12 +302,12 @@ def analyse_neuron(self, neuron: ASTModel) -> Tuple[Dict[str, ASTAssignment], Di
 
         equations_block = neuron.get_equations_blocks()[0]
 
-        delta_factors = ASTUtils.get_delta_factors_(neuron, equations_block)
         kernel_buffers = ASTUtils.generate_kernel_buffers(neuron, equations_block)
-        ASTUtils.replace_convolve_calls_with_buffers_(neuron, equations_block)
         ASTUtils.make_inline_expressions_self_contained(equations_block.get_inline_expressions())
         ASTUtils.replace_inline_expressions_through_defining_expressions(
             equations_block.get_ode_equations(), equations_block.get_inline_expressions())
+        delta_factors = ASTUtils.get_delta_factors_(neuron, equations_block)
+        ASTUtils.replace_convolve_calls_with_buffers_(neuron, equations_block)
 
         # Collect all equations with delay variables and replace ASTFunctionCall to ASTVariable wherever necessary
         equations_with_delay_vars_visitor = ASTEquationsWithDelayVarsVisitor()
@@ -382,12 +382,12 @@ def analyse_synapse(self, synapse: ASTModel) -> Dict[str, ASTAssignment]:
 
             equations_block = synapse.get_equations_blocks()[0]
 
-            delta_factors = ASTUtils.get_delta_factors_(synapse, equations_block)
             kernel_buffers = ASTUtils.generate_kernel_buffers(synapse, equations_block)
-            ASTUtils.replace_convolve_calls_with_buffers_(synapse, equations_block)
             ASTUtils.make_inline_expressions_self_contained(equations_block.get_inline_expressions())
             ASTUtils.replace_inline_expressions_through_defining_expressions(
                 equations_block.get_ode_equations(), equations_block.get_inline_expressions())
+            delta_factors = ASTUtils.get_delta_factors_(synapse, equations_block)
+            ASTUtils.replace_convolve_calls_with_buffers_(synapse, equations_block)
 
             analytic_solver, numeric_solver = self.ode_toolbox_analysis(synapse, kernel_buffers)
             self.analytic_solver[synapse.get_name()] = analytic_solver

pynestml/utils/ast_utils.py

@@ -1928,6 +1928,7 @@ def get_delta_factors_(cls, neuron: ASTModel, equations_block: ASTEquationsBlock
         For every occurrence of a convolution of the form `x^(n) = a * convolve(kernel, inport) + ...` where `kernel` is a delta function, add the element `(x^(n), inport) --> a` to the set.
         """
         delta_factors = {}
+
         for ode_eq in equations_block.get_ode_equations():
             var = ode_eq.get_lhs()
             expr = ode_eq.get_rhs()

tests/nest_tests/resources/iaf_psc_exp_multisynapse.nestml

@@ -16,23 +16,19 @@ model iaf_psc_exp_multisynapse_neuron:
         V_m mV = E_L  # membrane potential
         refr_t ms = 0 ms    # Refractory period timer
         is_refractory boolean = false
-        I_syn1 pA = 0 pA
-        I_syn2 pA = 0 pA
-        I_syn3 pA = 0 pA
 
     equations:
-        kernel I_kernel1 = exp(-t / tau_syn1)
+        kernel I_kernel1 = delta(t)
         kernel I_kernel2 = (e / tau_syn2) * t * exp(-t / tau_syn2)
         kernel I_kernel3 = -exp(-t / tau_syn3)
 
-        inline I_syn pA = (convolve(I_kernel1, spikes1) - convolve(I_kernel2, spikes2) + convolve(I_kernel3, spikes3)) * pA
+        recordable inline I_syn pA = (convolve(I_kernel1, spikes1) - convolve(I_kernel2, spikes2) + convolve(I_kernel3, spikes3)) * pA
 
         V_m' = -(V_m - E_L) / tau_m + (I_syn + I_e + I_stim) / C_m
 
     parameters:
         C_m pF = 250 pF           # Capacitance of the membrane
         tau_m ms = 10 ms          # Membrane time constant
-        tau_syn1 ms = .2ms        # Time constant of synaptic current.
         tau_syn2 ms = 2ms         # Time constant of synaptic current.
         tau_syn3 ms = 20ms        # Time constant of synaptic current.
         refr_T ms = 2 ms          # Duration of refractory period

tests/nest_tests/test_multisynapse.py

@@ -43,7 +43,7 @@ def test_multisynapse(self):
         input_path = os.path.join(os.path.realpath(os.path.join(
             os.path.dirname(__file__), "resources", "iaf_psc_exp_multisynapse.nestml")))
         target_path = "target"
-        logging_level = "INFO"
+        logging_level = "DEBUG"
         module_name = "nestmlmodule"
         suffix = "_nestml"
 
@@ -52,13 +52,13 @@ def test_multisynapse(self):
                              logging_level=logging_level,
                              module_name=module_name,
                              suffix=suffix)
-        nest.set_verbosity("M_ALL")
+        nest.Install(module_name)
 
         nest.ResetKernel()
-        nest.Install(module_name)
+        nest.set_verbosity("M_ALL")
+        nest.resolution = 0.1
 
         # network construction
-
         neuron = nest.Create("iaf_psc_exp_multisynapse_neuron_nestml")
 
         receptor_types = nest.GetStatus(neuron, "receptor_types")[0]
@@ -73,10 +73,10 @@ def test_multisynapse(self):
         nest.Connect(sg3, neuron, syn_spec={"receptor_type": receptor_types["SPIKES3"], "weight": 500., "delay": 0.1})
 
         mm = nest.Create("multimeter", params={"record_from": [
-                         "I_kernel1__X__spikes1", "I_kernel2__X__spikes2", "I_kernel3__X__spikes3"], "interval": 0.1})
+                         "I_syn", "I_kernel2__X__spikes2", "I_kernel3__X__spikes3"], "interval": nest.resolution})
         nest.Connect(mm, neuron)
 
-        vm_1 = nest.Create("voltmeter")
+        vm_1 = nest.Create("voltmeter", params={"interval": nest.resolution})
         nest.Connect(vm_1, neuron)
 
         # simulate
@@ -93,7 +93,7 @@ def test_multisynapse(self):
             ax[0].plot(V_m_timevec, V_m, label="V_m")
             ax[0].set_ylabel("voltage")
 
-            ax[1].plot(mm["times"], mm["I_kernel1__X__spikes1"], label="I_kernel1")
+            ax[1].plot(mm["times"], mm["I_syn"], label="I_syn")
             ax[1].set_ylabel("current")
 
             ax[2].plot(mm["times"], mm["I_kernel2__X__spikes2"], label="I_kernel2")
@@ -115,13 +115,13 @@ def test_multisynapse(self):
             fig.savefig("/tmp/test_multisynapse.png")
 
         # testing
-        np.testing.assert_almost_equal(V_m[-1], -72.93610279332957)
+        np.testing.assert_almost_equal(V_m[-1], -72.77625579314515)
 
     def test_multisynapse_with_vector_input_ports(self):
         input_path = os.path.join(os.path.realpath(os.path.join(
             os.path.dirname(__file__), "resources", "iaf_psc_exp_multisynapse_vectors.nestml")))
         target_path = "target"
-        logging_level = "INFO"
+        logging_level = "DEBUG"
         module_name = "nestml_module"
         suffix = "_nestml"
 
@@ -134,6 +134,7 @@ def test_multisynapse_with_vector_input_ports(self):
 
         nest.ResetKernel()
         nest.Install(module_name)
+        nest.resolution = 0.1
 
         # network construction
         neuron = nest.Create("iaf_psc_exp_multisynapse_vectors_neuron_nestml")
@@ -151,10 +152,10 @@ def test_multisynapse_with_vector_input_ports(self):
         nest.Connect(sg3, neuron, syn_spec={"receptor_type": receptor_types["SPIKES_3"], "weight": 500., "delay": 0.1})
 
         mm = nest.Create("multimeter", params={"record_from": [
-            "I_kernel1__X__spikes_1", "I_kernel2__X__spikes_2", "I_kernel3__X__spikes_3"], "interval": 0.1})
+            "I_kernel1__X__spikes_1", "I_kernel2__X__spikes_2", "I_kernel3__X__spikes_3"], "interval": nest.resolution})
         nest.Connect(mm, neuron)
 
-        vm_1 = nest.Create("voltmeter")
+        vm_1 = nest.Create("voltmeter", params={"interval": nest.resolution})
         nest.Connect(vm_1, neuron)
 
         # simulate
@@ -192,4 +193,4 @@ def test_multisynapse_with_vector_input_ports(self):
             fig.savefig("/tmp/test_multisynapse_vector.png")
 
         # testing
-        np.testing.assert_almost_equal(V_m[-1], -72.89041451202327)
+        np.testing.assert_almost_equal(V_m[-1], -72.77067117288824)