Merge pull request #2986 from janskaar/sli2py_iaf_psc_exp_multisynapse

bug fix iaf_psc_exp_multisynapse / sli2py test_iaf_psc_exp_multisynapse

models/iaf_psc_alpha_multisynapse.cpp

@@ -110,7 +110,6 @@ iaf_psc_alpha_multisynapse::Parameters_::Parameters_()
 iaf_psc_alpha_multisynapse::State_::State_()
   : I_const_( 0.0 )
   , V_m_( 0.0 )
-  , current_( 0.0 )
   , refractory_steps_( 0 )
 {
   y1_syn_.clear();
@@ -338,11 +337,9 @@ iaf_psc_alpha_multisynapse::update( Time const& origin, const long from, const l
       // neuron not refractory
       S_.V_m_ = V_.P30_ * ( S_.I_const_ + P_.I_e_ ) + V_.P33_ * S_.V_m_;
 
-      S_.current_ = 0.0;
       for ( size_t i = 0; i < P_.n_receptors_(); i++ )
       {
         S_.V_m_ += V_.P31_syn_[ i ] * S_.y1_syn_[ i ] + V_.P32_syn_[ i ] * S_.y2_syn_[ i ];
-        S_.current_ += S_.y2_syn_[ i ];
       }
 
       // lower bound of membrane potential

models/iaf_psc_alpha_multisynapse.h

@@ -214,8 +214,6 @@ class iaf_psc_alpha_multisynapse : public ArchivingNode
     std::vector< double > y2_syn_;
     //! This is the membrane potential RELATIVE TO RESTING POTENTIAL.
     double V_m_;
-    double current_; //! This is the current in a time step. This is only here
-                     //! to allow logging
 
     int refractory_steps_; //!< Number of refractory steps remaining
 
@@ -300,7 +298,7 @@ class iaf_psc_alpha_multisynapse : public ArchivingNode
     }
     else if ( elem == State_::I )
     {
-      return S_.current_;
+      return std::accumulate( S_.y2_syn_.begin(), S_.y2_syn_.end(), 0.0 );
     }
     else
     {

models/iaf_psc_exp_multisynapse.cpp

@@ -22,7 +22,6 @@
 
 #include "iaf_psc_exp_multisynapse.h"
 
-
 // Includes from libnestutil:
 #include "dict_util.h"
 #include "exceptions.h"
@@ -107,7 +106,6 @@ iaf_psc_exp_multisynapse::Parameters_::Parameters_()
 iaf_psc_exp_multisynapse::State_::State_()
   : I_const_( 0.0 )
   , V_m_( 0.0 )
-  , current_( 0.0 )
   , refractory_steps_( 0 )
 {
   i_syn_.clear();
@@ -313,17 +311,16 @@ iaf_psc_exp_multisynapse::update( const Time& origin, const long from, const lon
     {
       S_.V_m_ = S_.V_m_ * V_.P22_ + ( P_.I_e_ + S_.I_const_ ) * V_.P20_; // not sure about this
 
-      S_.current_ = 0.0;
       for ( size_t i = 0; i < P_.n_receptors_(); i++ )
       {
         S_.V_m_ += V_.P21_syn_[ i ] * S_.i_syn_[ i ];
-        S_.current_ += S_.i_syn_[ i ]; // not sure about this
       }
     }
     else
     {
       --S_.refractory_steps_; // neuron is absolute refractory
     }
+
     for ( size_t i = 0; i < P_.n_receptors_(); i++ )
     {
       // exponential decaying PSCs

models/iaf_psc_exp_multisynapse.h

@@ -211,9 +211,7 @@ class iaf_psc_exp_multisynapse : public ArchivingNode
 
     double I_const_; //!< synaptic dc input current, variable 0
     std::vector< double > i_syn_;
-    double V_m_;     //!< membrane potential, variable 2
-    double current_; //!< This is the current in a time step. This is only
-                     //!< here to allow logging
+    double V_m_; //!< membrane potential, variable 2
 
     //! absolute refractory counter (no membrane potential propagation)
     int refractory_steps_;
@@ -300,7 +298,7 @@ class iaf_psc_exp_multisynapse : public ArchivingNode
     }
     else if ( elem == State_::I )
     {
-      return S_.current_;
+      return std::accumulate( S_.i_syn_.begin(), S_.i_syn_.end(), 0.0 );
     }
     else
     {

testsuite/pytests/test_iaf_psc_alpha_multisynapse.py → testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_multisynapse.py

@@ -28,6 +28,7 @@
 import numpy as np
 import numpy.testing as nptest
 import pytest
+from scipy.linalg import expm
 
 
 @pytest.fixture(autouse=True)
@@ -37,8 +38,7 @@ def reset():
 
 def alpha_fn(t, tau_syn):
     vals = np.zeros_like(t)
-    zero_inds = t <= 0.0
-    nonzero_inds = ~zero_inds
+    nonzero_inds = t > 0.0
     vals[nonzero_inds] = np.e / tau_syn * t[nonzero_inds] * np.exp(-t[nonzero_inds] / tau_syn)
     return vals
 
@@ -50,6 +50,16 @@ def test_I_syn_1_in_recordables():
     assert "I_syn_1" in nrn.get("recordables")
 
 
+def alpha_psc_voltage_response(t, tau_syn, tau_m, C_m, w):
+    vals = np.zeros_like(t)
+    nonzero_inds = t > 0.0
+    A = np.array([[-1.0 / tau_syn, 0.0, 0.0], [1.0, -1.0 / tau_syn, 0.0], [0.0, 1.0 / C_m, -1.0 / tau_m]])
+
+    expAt = expm(A[None, ...] * t[nonzero_inds, None, None])  # shape (t, 3, 3)
+    vals[nonzero_inds] = expAt[:, 2, 0] * w * np.e / tau_syn  # first two state variables are 0
+    return vals
+
+
 def test_resize_recordables():
     """
     Test resizing of recordables.
@@ -80,40 +90,57 @@ def test_simulation_against_analytical_soln():
     from multiple different synaptic ports are the same as the analytical solution.
     """
 
-    tau_syn = [2.0, 20.0, 60.0, 100.0]
-    delays = [100.0, 200.0, 500.0, 1200.0]
-    weight = 1.0
-    spike_time = 10.0
-    simtime = 2500.0
+    tau_syns = [2.0, 20.0, 60.0, 100.0]
+    delays = [7.0, 5.0, 2.0, 1.0]
+    weights = [30.0, 50.0, 20.0, 10.0]
+    C_m = 250.0
+    tau_m = 15.0
+    spike_time = 1.0
+    simtime = 100.0
+    dt = 1.0
+
+    nest.set(resolution=dt)
 
     nrn = nest.Create(
         "iaf_psc_alpha_multisynapse",
         params={
-            "C_m": 250.0,
+            "C_m": C_m,
             "E_L": 0.0,
             "V_m": 0.0,
             "V_th": 1500.0,
             "I_e": 0.0,
-            "tau_m": 15.0,
-            "tau_syn": tau_syn,
+            "tau_m": tau_m,
+            "tau_syn": tau_syns,
         },
     )
-    sg = nest.Create("spike_generator", params={"spike_times": [spike_time]})
 
-    for i, syn_id in enumerate(range(1, 5)):
-        syn_spec = {"synapse_model": "static_synapse", "delay": delays[i], "weight": weight, "receptor_type": syn_id}
+    sg = nest.Create("spike_generator", params={"spike_times": [spike_time]})
 
+    for syn_idx, (delay, weight) in enumerate(zip(delays, weights)):
+        syn_spec = {
+            "synapse_model": "static_synapse",
+            "delay": delay,
+            "weight": weight,
+            "receptor_type": syn_idx + 1,
+        }
         nest.Connect(sg, nrn, conn_spec="one_to_one", syn_spec=syn_spec)
 
-    mm = nest.Create("multimeter", params={"record_from": ["I_syn_1", "I_syn_2", "I_syn_3", "I_syn_4"]})
+    mm = nest.Create(
+        "multimeter",
+        params={"record_from": ["I_syn_1", "I_syn_2", "I_syn_3", "I_syn_4", "V_m", "I_syn"], "interval": dt},
+    )
 
     nest.Connect(mm, nrn)
     nest.Simulate(simtime)
     times = mm.get("events", "times")
-    I_syn = np.sum([mm.get("events", f"I_syn_{i}") for i in range(1, 5)], axis=0)
 
-    I_syn_analytical = np.zeros_like(times, dtype=np.float64)
-    for i in range(4):
-        I_syn_analytical += alpha_fn(times - delays[i] - spike_time, tau_syn[i])
+    I_syns_analytical = []
+    V_m_analytical = np.zeros_like(times)
+    for weight, delay, tau_s in zip(weights, delays, tau_syns):
+        I_syns_analytical.append(alpha_fn(times - delay - spike_time, tau_s) * weight)
+        V_m_analytical += alpha_psc_voltage_response(times - delay - spike_time, tau_s, tau_m, C_m, weight)
+
+    for idx, I_syn_analytical in enumerate(I_syns_analytical):
+        nptest.assert_array_almost_equal(mm.get("events", f"I_syn_{idx+1}"), I_syn_analytical)
 
-    nptest.assert_array_almost_equal(I_syn, I_syn_analytical)
+    nptest.assert_array_almost_equal(mm.get("events", "V_m"), V_m_analytical)

testsuite/pytests/sli2py_neurons/test_iaf_psc_exp_multisynapse.py

@@ -0,0 +1,157 @@
+# -*- coding: utf-8 -*-
+#
+# test_iaf_psc_exp_multisynapse.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/>.
+
+"""
+Test ``iaf_psc_exp_multisynapse`` recordables and simulated PSCs against expectation.
+"""
+
+
+import nest
+import numpy as np
+import numpy.testing as nptest
+import pytest
+
+
+@pytest.fixture(autouse=True)
+def reset():
+    nest.ResetKernel()
+
+
+def exp_psc_fn(t, tau_syn):
+    vals = np.zeros_like(t)
+    nonzero_inds = t > 0.0
+    vals[nonzero_inds] = np.exp(-t[nonzero_inds] / tau_syn)
+    return vals
+
+
+def exp_psc_voltage_response(t, tau_syn, tau_m, C_m, w):
+    vals = np.zeros_like(t)
+    nonzero_inds = t > 0.0
+    delta_e = np.exp(-t[nonzero_inds] / tau_m) - np.exp(-t[nonzero_inds] / tau_syn)
+    vals[nonzero_inds] = w / (C_m * (1.0 / tau_syn - 1.0 / tau_m)) * delta_e
+    return vals
+
+
+def test_set_synaptic_time_constants():
+    """Tests that synaptic time constants can be set correctly"""
+    taus = [2.0, 20.0, 60.0, 100.0]
+    nrn = nest.Create("iaf_psc_exp_multisynapse")
+    nrn.set(tau_syn=taus)
+    nptest.assert_array_almost_equal(nrn.get("tau_syn"), taus)
+
+
+def test_simulation_against_analytical_solution():
+    """
+    Test simulated PSCs against analytical expectation.
+
+    This test checks that the integration of the exponential currents of inputs
+    from multiple different synaptic ports are the same as the analytical solution.
+    """
+
+    tau_syns = [2.0, 20.0, 60.0, 100.0]
+    delays = [7.0, 5.0, 2.0, 1.0]
+    weights = [30.0, 50.0, 20.0, 10.0]
+    C_m = 250.0
+    tau_m = 15.0
+    spike_time = 0.1
+    simtime = 8.0
+    dt = 0.1
+
+    nest.set(resolution=dt)
+
+    nrn = nest.Create(
+        "iaf_psc_exp_multisynapse",
+        params={
+            "C_m": C_m,
+            "E_L": 0.0,
+            "V_m": 0.0,
+            "V_th": 1500.0,
+            "I_e": 0.0,
+            "tau_m": tau_m,
+            "tau_syn": tau_syns,
+        },
+    )
+
+    sg = nest.Create("spike_generator", params={"spike_times": [spike_time]})
+
+    for syn_idx, (delay, weight) in enumerate(zip(delays, weights)):
+        syn_spec = {
+            "synapse_model": "static_synapse",
+            "delay": delay,
+            "weight": weight,
+            "receptor_type": syn_idx + 1,
+        }
+        nest.Connect(sg, nrn, conn_spec="one_to_one", syn_spec=syn_spec)
+
+    mm = nest.Create(
+        "multimeter",
+        params={"record_from": ["I_syn_1", "I_syn_2", "I_syn_3", "I_syn_4", "V_m", "I_syn"], "interval": dt},
+    )
+
+    nest.Connect(mm, nrn, syn_spec={"delay": 0.1})
+    nest.Simulate(simtime)
+    times = mm.get("events", "times")
+
+    I_syns_analytical = []
+    V_m_analytical = np.zeros_like(times)
+    for weight, delay, tau_s in zip(weights, delays, tau_syns):
+        I_syns_analytical.append(exp_psc_fn(times - delay - spike_time, tau_s) * weight)
+        V_m_analytical += exp_psc_voltage_response(times - delay - spike_time, tau_s, tau_m, C_m, weight)
+
+    for idx, I_syn_analytical in enumerate(I_syns_analytical):
+        nptest.assert_array_almost_equal(mm.get("events", f"I_syn_{idx+1}"), I_syn_analytical)
+    nptest.assert_array_almost_equal(mm.get("events", "V_m"), V_m_analytical)
+
+
+# The following tests address #800
+# - Test that the default recordables are V_m, w and I_syn_1
+# - Test that the recordable I_syn's change when changing the number of receptor ports
+
+
+def test_default_recordables():
+    nrn = nest.Create("iaf_psc_exp_multisynapse")
+    recordables = nrn.get("recordables")
+    assert len(recordables) == 3
+    assert "I_syn" in recordables
+    assert "I_syn_1" in recordables
+    assert "V_m" in recordables
+
+
+def test_resize_recordables():
+    """
+    Test resizing of recordables.
+
+    This test ensures that recordables are updated correctly when the number
+    of synaptic ports are changed.
+    """
+
+    tau_syn1 = [5.0, 1.0, 25.0]
+    tau_syn2 = [5.0, 1.0]
+    tau_syn3 = [5.0, 1.0, 25.0, 50.0]
+
+    nrn = nest.Create("iaf_psc_alpha_multisynapse", params={"tau_syn": tau_syn1})
+    assert len(nrn.recordables) == 5
+
+    nrn.set(tau_syn=tau_syn2)
+    assert len(nrn.recordables) == 4
+
+    nrn.set(tau_syn=tau_syn3)
+    assert len(nrn.recordables) == 6