diff --git a/testsuite/pytests/sli2py_neurons/test_iaf_psc_exp_ps.py b/testsuite/pytests/sli2py_neurons/test_iaf_psc_exp_ps.py new file mode 100644 index 0000000000..88d69ec325 --- /dev/null +++ b/testsuite/pytests/sli2py_neurons/test_iaf_psc_exp_ps.py @@ -0,0 +1,97 @@ +# -*- coding: utf-8 -*- +# +# test_iaf_psc_exp_ps.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 . + +""" +Name: testsuite::test_iaf_psc_exp_ps + +Synopsis: (test_iaf_psc_exp_ps) run -> compares response to current step with analytical solution + +Description: +A DC current is injected into the neuron using a current generator +device. The membrane potential is recorder, and compared to the expected +analytical solution, with a spike happening off-grid. +""" + +import nest +import numpy as np +import numpy.testing as nptest +import pytest + + +def test_iaf_psc_exp_ps_dc_input(): + dt = 0.1 + dc_amp = 1010.0 + + nest.ResetKernel() + nest.set(resolution=dt, local_num_threads=1) + + dc_gen = nest.Create("dc_generator", {"amplitude": dc_amp}) + nrn = nest.Create("iaf_psc_exp_ps", 1) + vm = nest.Create("voltmeter", {"interval": 0.1}) + + syn_spec = {"synapse_model": "static_synapse", "weight": 1.0, "delay": dt} + nest.Connect(dc_gen, nrn, syn_spec=syn_spec) + nest.Connect(vm, nrn, syn_spec=syn_spec) + + nest.Simulate(8.0) + + times = vm.get("events", "times") + times -= dt # account for delay to multimeter + + tau_m = nrn.get("tau_m") + R = tau_m / nrn.get("C_m") + theta = nrn.get("V_th") + E_L = nrn.get("E_L") + + # array for analytical solution + V_m_analytical = np.empty_like(times) + V_m_analytical[:] = nrn.get("E_L") + + # first index for which the DC current is received by neuron. + # DC current will be integrated from this time step + start_index = 1 + + # analytical solution without delay and threshold on a grid + vm_soln = E_L + (1 - np.exp(-times / tau_m)) * R * dc_amp + + # exact time at which the neuron will spike + exact_spiketime = -tau_m * np.log(1 - (theta - E_L) / (R * dc_amp)) + + # offset from grid point + time_offset = exact_spiketime - (exact_spiketime // dt) * dt + + # solution calculated on the grid, with t0 being the exact spike time + vm_soln_offset = E_L + (1 - np.exp(-(times - time_offset + dt) / tau_m)) * R * dc_amp + + # rise until threshold + V_m_analytical[start_index:] = vm_soln[:-start_index] + + # set refractory potential + # first index after spike, offset by time at which DC current arrives + crossing_ind = int(exact_spiketime // dt + 1) + start_index + num_ref = int(nrn.get("t_ref") / dt) + V_m_analytical[crossing_ind : crossing_ind + num_ref] = nrn.get("V_reset") + + # rise after refractory period + num_inds = len(times) - crossing_ind - num_ref + V_m_analytical[crossing_ind + num_ref :] = vm_soln_offset[:num_inds] + + nptest.assert_array_almost_equal(V_m_analytical, vm.get("events", "V_m")) diff --git a/testsuite/unittests/test_iaf_psc_exp_ps.sli b/testsuite/unittests/test_iaf_psc_exp_ps.sli deleted file mode 100644 index 14479a98de..0000000000 --- a/testsuite/unittests/test_iaf_psc_exp_ps.sli +++ /dev/null @@ -1,172 +0,0 @@ -/* - * test_iaf_psc_exp_ps.sli - * - * 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 . - * - */ - - /** @BeginDocumentation - Name: testsuite::test_iaf_psc_exp_ps - sli script for overall test of iaf_psc_exp_ps model - - Synopsis: (test_iaf_psc_exp_ps) run -> compares response to current step with reference data - - Description: - test_iaf_psc_exp_ps.sli is an overall test of the iaf_psc_exp_ps model connected - to some useful devices. - - A DC current is injected into the neuron using a current generator - device. The membrane potential as well as the spiking activity are - recorded by corresponding devices. - - It can be observed how the current charges the membrane, a spike - is emitted, the neuron becomes absolute refractory, and finally - starts to recover. - - The timing of the various events on the simulation grid is of - particular interest and crucial for the consistency of the - simulation scheme. - - Although 0.1 cannot be represented in the IEEE double data type, it - is safe to simulate with a resolution (computation step size) of 0.1 - ms because by default nest is built with a timebase enabling exact - representation of 0.1 ms. - - The expected output is documented and briefly commented at the end of - the script. - - Other test programs discuss the various aspects of this script in detail, - see the SeeAlso key below. - - Author: Jeyashree Krishnan, 2017 - SeeAlso: iaf_psc_exp, testsuite::test_iaf_i0, testsuite::test_iaf_i0_refractory, testsuite::test_iaf_dc, testsuite::test_iaf_psc_exp_ps_lossless -*/ - - -(unittest) run -/unittest using - -0.1 /h Set - -ResetKernel - -<< - /local_num_threads 1 - /resolution h ->> SetKernelStatus - -/iaf_psc_exp_ps Create /neuron Set - -/dc_generator Create /dc_gen Set -dc_gen << /amplitude 1000. >> SetStatus - -/voltmeter << /time_in_steps true /interval h >> Create /vm Set - -/spike_recorder Create /sp_det Set - - -dc_gen neuron 1.0 h Connect -vm neuron 1.0 h Connect -neuron sp_det 1.0 h Connect - -8 Simulate - -sp_det /events get /times get First stack % prints spike time - -{ % reference data - dup Transpose First /test_times Set % times of reference - - vm [/events [/times /V_m]] get cva % array of recorded data - 6 ToUnitTestPrecision % to precision of reference - Transpose dup == % all recorded tuples - {First test_times exch MemberQ } Select % those with reference - eq % compare -} - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% Expected output of this program: -% -% The output send to std::cout is a superposition of the output of -% the voltmeter and the spike recorder. Both, voltmeter and spike -% recorder are connected to the same neuron. -% -% time (in steps) voltage (in mV) -[ -[ 1 -70 ] %<----- The earliest time dc_gen can be switched on. -[ 2 -70 ] %<----- The DC current arrives at the neuron, it is -[ 3 -69.602 ] %<- reflected in the neuron's state variable y0, -[ 4 -69.2079 ] % | (initial condition) but has not yet affected -[ 5 -68.8178 ] % | the membrane potential. -[ 6 -68.4316 ] % | -[ 7 -68.0492 ] % --- the effect of the DC current is visible in the -[ 8 -67.6706 ] % membrane potential -[ 9 -67.2958 ] % -[ 10 -66.9247 ] % -% -% ... -% -[ 45 -56.0204 ] % -[ 46 -55.7615 ] % -[ 47 -55.5051 ] % -[ 48 -55.2513 ] % -[ 49 -55.0001 ] % -[ 50 -70 ] % <---- The membrane potential crossed threshold in the -[ 51 -70 ] % step 4.9 ms -> 5.0 ms. The membrane potential is -[ 52 -70 ] % reset (no super-threshold values can be observed). -[ 53 -70 ] % The precise spike time is reported at 4.90004 ms. -[ 54 -70 ] % -[ 55 -70 ] % -[ 56 -70 ] % -[ 57 -70 ] % -[ 58 -70 ] % -[ 59 -70 ] % -[ 60 -70 ] % -[ 61 -70 ] % -[ 62 -70 ] % -[ 63 -70 ] % -[ 64 -70 ] % -[ 65 -70 ] % -[ 66 -70 ] % -[ 67 -70 ] % -[ 68 -70 ] % -[ 69 -70 ] % -[ 70 -69.6021 ] % <---- Since the neuron uses precise spike times that are not -[ 71 -69.2081 ] % locked to the grid, the refractory period ended after -[ 72 -68.818 ] % 2.0 ms during the timestep 6.9 ms -> 7.0 ms. The membrane -[ 73 -68.4317 ] % potential has already started to increase again. -[ 74 -68.0493 ] % -[ 75 -67.6707 ] % -[ 76 -67.2959 ] % -[ 77 -66.9248 ] % -[ 78 -66.5574 ] % -[ 79 -66.1936 ] % <-- -] % | -% | -% - The simulation was run for 8.0 ms. However, in the step -% 7.9 ms -> 8.0 ms the voltmeter necessarily receives the -% voltages that occurred at time 7.9 ms (delay h). This -% results in different end times of the recorded voltage -% traces at different resolutions. In the current -% simulation kernel there is no general cure for this -% problem. One workaround is to end the simulation script -% with "h Simulate", thereby making the script resolution -% dependent. -% - -exch assert_or_die -