diff --git a/testsuite/pytests/conftest.py b/testsuite/pytests/conftest.py index 8e073adb9b..02440c1e21 100644 --- a/testsuite/pytests/conftest.py +++ b/testsuite/pytests/conftest.py @@ -31,7 +31,6 @@ def test_gsl(): pass """ -import dataclasses import os import pathlib import subprocess @@ -45,7 +44,6 @@ def test_gsl(): # Ignore it during test collection collect_ignore = ["utilities"] -import testsimulation # noqa import testutil # noqa @@ -179,24 +177,3 @@ def skipif_incompatible_mpi(request, subprocess_compatible_mpi): if not subprocess_compatible_mpi and request.node.get_closest_marker("skipif_incompatible_mpi"): pytest.skip("skipped because MPI is incompatible with subprocess") - - -@pytest.fixture(autouse=True) -def simulation_class(request): - return getattr(request, "param", testsimulation.Simulation) - - -@pytest.fixture -def simulation(request): - marker = request.node.get_closest_marker("simulation") - sim_cls = marker.args[0] if marker else testsimulation.Simulation - sim = sim_cls(*(request.getfixturevalue(field.name) for field in dataclasses.fields(sim_cls))) - nest.ResetKernel() - if getattr(sim, "set_resolution", True): - nest.resolution = sim.resolution - nest.local_num_threads = sim.local_num_threads - return sim - - -# Inject the root simulation fixtures into this module to be always available. -testutil.create_dataclass_fixtures(testsimulation.Simulation, __name__) diff --git a/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha.py b/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha.py deleted file mode 100644 index 8ef1815eb0..0000000000 --- a/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha.py +++ /dev/null @@ -1,431 +0,0 @@ -# -*- coding: utf-8 -*- -# -# test_iaf_psc_alpha.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 . - -import dataclasses -import math - -import nest -import numpy as np -import pytest -import testsimulation -import testutil -from scipy.special import lambertw - -# Notes: -# * copy docs -# * add docs & examples -# * add comments -# * restructure -# * test min delay stuff tests sim indep of mindelay, move out? -# * `test_iaf_ps_dc_accuracy` tests kernel precision, move out? - - -@dataclasses.dataclass -class IAFPSCAlphaSimulation(testsimulation.Simulation): - def setup(self): - self.neuron = nest.Create("iaf_psc_alpha") - vm = self.voltmeter = nest.Create("voltmeter") - vm.interval = self.resolution - sr = self.spike_recorder = nest.Create("spike_recorder") - nest.Connect(vm, self.neuron, syn_spec={"weight": 1.0, "delay": self.delay}) - nest.Connect(self.neuron, sr, syn_spec={"weight": 1.0, "delay": self.delay}) - - @property - def spikes(self): - return np.column_stack( - ( - self.spike_recorder.events["senders"], - self.spike_recorder.events["times"], - ) - ) - - -@dataclasses.dataclass -class MinDelaySimulation(IAFPSCAlphaSimulation): - amplitude: float = 1000.0 - min_delay: float = 0.0 - - def setup(self): - dc = self.dc_generator = nest.Create("dc_generator") - dc.amplitude = self.amplitude - - super().setup() - - nest.Connect( - dc, - self.neuron, - syn_spec={"weight": 1.0, "delay": self.delay}, - ) - - -@testutil.use_simulation(IAFPSCAlphaSimulation) -class TestIAFPSCAlpha: - def test_iaf_psc_alpha(self, simulation): - dc = simulation.dc_generator = nest.Create("dc_generator") - dc.amplitude = 1000 - - simulation.setup() - - nest.Connect( - dc, - simulation.neuron, - syn_spec={"weight": 1.0, "delay": simulation.resolution}, - ) - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expected_default) - assert actual == expected - - @pytest.mark.parametrize("duration", [20.0]) - def test_iaf_psc_alpha_fudge(self, simulation): - simulation.setup() - - tau_m = 20 - tau_syn = 0.5 - C_m = 250.0 - a = tau_m / tau_syn - b = 1.0 / tau_syn - 1.0 / tau_m - t_max = 1.0 / b * (-lambertw(-math.exp(-1.0 / a) / a, k=-1) - 1.0 / a).real - V_max = ( - math.exp(1) - / (tau_syn * C_m * b) - * ((math.exp(-t_max / tau_m) - math.exp(-t_max / tau_syn)) / b - t_max * math.exp(-t_max / tau_syn)) - ) - simulation.neuron.set(tau_m=tau_m, tau_syn_ex=tau_syn, tau_syn_in=tau_syn, C_m=C_m) - sg = nest.Create( - "spike_generator", - params={"precise_times": False, "spike_times": [simulation.resolution]}, - ) - nest.Connect( - sg, - simulation.neuron, - syn_spec={"weight": float(1.0 / V_max), "delay": simulation.resolution}, - ) - - results = simulation.simulate() - - actual_t_max = results[np.argmax(results[:, 1]), 0] - assert actual_t_max == pytest.approx(t_max + 0.2, abs=0.05) - - def test_iaf_psc_alpha_i0(self, simulation): - simulation.setup() - - simulation.neuron.I_e = 1000 - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expected_i0) - assert actual == expected - assert simulation.spikes == pytest.approx(expected_i0_t) - - @pytest.mark.parametrize("resolution", [0.1, 0.2, 0.5, 1.0]) - def test_iaf_psc_alpha_i0_refractory(self, simulation): - simulation.setup() - - simulation.neuron.I_e = 1450 - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expected_i0_refr) - assert actual == expected - - -@pytest.mark.parametrize("min_delay", [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1.0, 2.0]) -@pytest.mark.parametrize("delay, duration", [(2.0, 10.5)]) -@testutil.use_simulation(MinDelaySimulation) -class TestMinDelayUsingIAFPSCAlpha: - def test_iaf_psc_alpha_mindelay_create(self, simulation, min_delay): - simulation.setup() - - # Connect 2 throwaway neurons with `min_delay` to force `min_delay` - nest.Connect(*nest.Create("iaf_psc_alpha", 2), syn_spec={"delay": min_delay, "weight": 1.0}) - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expected_mindelay) - assert actual == expected - - def test_iaf_psc_alpha_mindelay_set(self, simulation, min_delay, delay): - nest.set(min_delay=min_delay, max_delay=delay) - nest.SetDefaults("static_synapse", {"delay": delay}) - - simulation.setup() - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expected_mindelay) - assert actual == expected - - def test_iaf_psc_alpha_mindelay_simblocks(self, simulation, min_delay, delay): - nest.set(min_delay=min_delay, max_delay=delay) - nest.SetDefaults("static_synapse", {"delay": delay}) - - simulation.setup() - - for _ in range(22): - nest.Simulate(0.5) - # duration=0 so that `simulation.simulate` is noop but - # still extracts results for us. - simulation.duration = 0 - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expected_mindelay) - assert actual == expected - - -def test_kernel_precision(): - nest.ResetKernel() - nest.set(tics_per_ms=2**14, resolution=2**0) - assert math.frexp(nest.ms_per_tic) == (0.5, -13) - - -@dataclasses.dataclass -class DCAccuracySimulation(testsimulation.Simulation): - # Don't autoset the resolution in the fixture, we do it in setup. - set_resolution = False - model: str = "iaf_psc_alpha" - params: dict = dataclasses.field(default_factory=dict) - - def setup(self): - nest.ResetKernel() - nest.set(tics_per_ms=2**14, resolution=self.resolution) - self.neuron = nest.Create(self.model, params=self.params) - - -@testutil.use_simulation(DCAccuracySimulation) -@pytest.mark.parametrize( - "model", - [ - "iaf_psc_alpha_ps", - "iaf_psc_delta_ps", - "iaf_psc_exp_ps", - "iaf_psc_exp_ps_lossless", - ], -) -@pytest.mark.parametrize("resolution", [2**i for i in range(0, -14, -1)]) -class TestIAFPSDCAccuracy: - @pytest.mark.parametrize( - "params", - [ - { - "E_L": 0.0, # resting potential in mV - "V_m": 0.0, # initial membrane potential in mV - "V_th": 2000.0, # spike threshold in mV - "I_e": 1000.0, # DC current in pA - "tau_m": 10.0, # membrane time constant in ms - "C_m": 250.0, # membrane capacity in pF - } - ], - ) - @pytest.mark.parametrize("duration, tolerance", [(5, 1e-13), (500.0, 1e-9)]) - def test_iaf_ps_dc_accuracy(self, simulation, duration, tolerance, params): - simulation.run() - # Analytical solution - V = params["I_e"] * params["tau_m"] / params["C_m"] * (1.0 - math.exp(-duration / params["tau_m"])) - # Check that membrane potential is within tolerance of analytical solution. - assert math.fabs(simulation.neuron.V_m - V) < tolerance - - @pytest.mark.parametrize( - "params", - [ - { - "E_L": 0.0, # resting potential in mV - "V_m": 0.0, # initial membrane potential in mV - "V_th": 15.0, # spike threshold in mV - "I_e": 1000.0, # DC current in pA - "tau_m": 10.0, # membrane time constant in ms - "C_m": 250.0, # membrane capacity in pF - } - ], - ) - @pytest.mark.parametrize("duration, tolerance", [(5, 1e-13)]) - def test_iaf_ps_dc_t_accuracy(self, simulation, params, tolerance): - simulation.run() - t = -params["tau_m"] * math.log(1.0 - (params["C_m"] * params["V_th"]) / (params["tau_m"] * params["I_e"])) - assert math.fabs(simulation.neuron.t_spike - t) < tolerance - - -expected_default = np.array( - [ - [0.1, -70], - [0.2, -70], - [0.3, -69.602], - [0.4, -69.2079], - [0.5, -68.8178], - [0.6, -68.4316], - [0.7, -68.0492], - [0.8, -67.6706], - [0.9, -67.2958], - [1.0, -66.9247], - [4.5, -56.0204], - [4.6, -55.7615], - [4.7, -55.5051], - [4.8, -55.2513], - [4.9, -55.0001], - [5.0, -70], - [5.1, -70], - [5.2, -70], - [5.3, -70], - [5.4, -70], - [5.5, -70], - [5.6, -70], - [5.7, -70], - [5.8, -70], - [5.9, -70], - [6.0, -70], - [6.1, -70], - [6.2, -70], - [6.3, -70], - [6.4, -70], - [6.5, -70], - [6.6, -70], - [6.7, -70], - [6.8, -70], - [6.9, -70], - [7.0, -70], - [7.1, -69.602], - [7.2, -69.2079], - [7.3, -68.8178], - [7.4, -68.4316], - [7.5, -68.0492], - [7.6, -67.6706], - [7.7, -67.2958], - [7.8, -66.9247], - [7.9, -66.5572], - ] -) - -expected_i0 = np.array( - [ - [0.1, -69.602], - [0.2, -69.2079], - [0.3, -68.8178], - [0.4, -68.4316], - [0.5, -68.0492], - [4.3, -56.0204], - [4.4, -55.7615], - [4.5, -55.5051], - [4.6, -55.2513], - [4.7, -55.0001], - [4.8, -70], - [4.9, -70], - [5.0, -70], - ] -) - -expected_i0_t = np.array( - [ - [1, 4.8], - ] -) - -expected_i0_refr = np.array( - [ - [0.1, -69.4229], - [0.2, -68.8515], - [0.3, -68.2858], - [0.4, -67.7258], - [0.5, -67.1713], - [0.6, -66.6223], - [0.7, -66.0788], - [0.8, -65.5407], - [0.9, -65.008], - [1.0, -64.4806], - [1.1, -63.9584], - [1.2, -63.4414], - [1.3, -62.9295], - [1.4, -62.4228], - [1.5, -61.9211], - [1.6, -61.4243], - [1.7, -60.9326], - [1.8, -60.4457], - [1.9, -59.9636], - [2.0, -59.4864], - [2.1, -59.0139], - [2.2, -58.5461], - [2.3, -58.0829], - [2.4, -57.6244], - [2.5, -57.1704], - [2.6, -56.721], - [2.7, -56.276], - [2.8, -55.8355], - [2.9, -55.3993], - [3.0, -70], - [3.1, -70], - [3.2, -70], - [3.3, -70], - [3.4, -70], - [3.5, -70], - [3.6, -70], - [3.7, -70], - [3.8, -70], - [3.9, -70], - [4.0, -70], - [4.1, -70], - [4.2, -70], - [4.3, -70], - [4.4, -70], - [4.5, -70], - [4.6, -70], - [4.7, -70], - [4.8, -70], - [4.9, -70], - [5.0, -70], - [ - 5.1, - -69.4229, - ], - [5.2, -68.8515], - [5.3, -68.2858], - [5.4, -67.7258], - [5.5, -67.1713], - [5.6, -66.6223], - [5.7, -66.0788], - [5.8, -65.5407], - [5.9, -65.008], - [6.0, -64.4806], - [6.1, -63.9584], - [6.2, -63.4414], - [6.3, -62.9295], - [6.4, -62.4228], - [6.5, -61.9211], - [6.6, -61.4243], - [6.7, -60.9326], - [6.8, -60.4457], - [6.9, -59.9636], - ] -) - -expected_mindelay = np.array( - [ - [1.000000e00, -7.000000e01], - [2.000000e00, -7.000000e01], - [3.000000e00, -6.655725e01], - [4.000000e00, -6.307837e01], - [5.000000e00, -5.993054e01], - [6.000000e00, -5.708227e01], - [7.000000e00, -7.000000e01], - [8.000000e00, -7.000000e01], - [9.000000e00, -6.960199e01], - [1.000000e01, -6.583337e01], - ] -) diff --git a/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha_1to2.py b/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha_1to2.py deleted file mode 100644 index ce7fcc9d8d..0000000000 --- a/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha_1to2.py +++ /dev/null @@ -1,155 +0,0 @@ -# -*- coding: utf-8 -*- -# -# test_iaf_psc_alpha_1to2.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 . - -import dataclasses - -import nest -import numpy as np -import pytest -import testsimulation -import testutil - - -@dataclasses.dataclass -class IAFPSCAlpha1to2Simulation(testsimulation.Simulation): - weight: float = 100.0 - delay: float = None - min_delay: float = None - - def __post_init__(self): - self.syn_spec = {"weight": self.weight} - if self.delay is not None: - self.syn_spec["delay"] = self.delay - - def setup(self): - n1, n2 = self.neurons = nest.Create("iaf_psc_alpha", 2) - n1.I_e = 1450.0 - vm = self.voltmeter = nest.Create("voltmeter") - vm.interval = self.resolution - vm_spec = {} - if self.delay is not None: - vm_spec["delay"] = self.delay - nest.Connect(vm, n2, syn_spec=vm_spec) - nest.Connect(n1, n2, syn_spec=self.syn_spec) - - -@pytest.mark.parametrize("resolution", [0.1, 0.2, 0.5, 1.0]) -@testutil.use_simulation(IAFPSCAlpha1to2Simulation) -class TestIAFPSCAlpha1to2WithMultiRes: - @pytest.mark.parametrize("delay", [1.0]) - def test_1to2(self, simulation): - simulation.setup() - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expect_default) - assert actual == expected - - def test_default_delay(self, simulation): - nest.SetDefaults("static_synapse", {"delay": 1.0}) - simulation.setup() - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expect_default) - assert actual == expected - - -@testutil.use_simulation(IAFPSCAlpha1to2Simulation) -@pytest.mark.parametrize("delay,resolution", [(2.0, 0.1)]) -@pytest.mark.parametrize("min_delay", [0.1, 0.5, 2.0]) -def test_mindelay_invariance(simulation): - assert simulation.min_delay <= simulation.delay - nest.set(min_delay=simulation.min_delay, max_delay=simulation.delay) - simulation.setup() - results = simulation.simulate() - actual, expected = testutil.get_comparable_timesamples(results, expect_inv) - assert actual == expected - - -expect_default = np.array( - [ - [2.5, -70], - [2.6, -70], - [2.7, -70], - [2.8, -70], - [2.9, -70], - [3.0, -70], - [3.1, -70], - [3.2, -70], - [3.3, -70], - [3.4, -70], - [3.5, -70], - [3.6, -70], - [3.7, -70], - [3.8, -70], - [3.9, -70], - [4.0, -70], - [4.1, -69.9974], - [4.2, -69.9899], - [4.3, -69.9781], - [4.4, -69.9624], - [4.5, -69.9434], - [4.6, -69.9213], - [4.7, -69.8967], - [4.8, -69.8699], - [4.9, -69.8411], - [5.0, -69.8108], - [5.1, -69.779], - [5.2, -69.7463], - [5.3, -69.7126], - [5.4, -69.6783], - [5.5, -69.6435], - [5.6, -69.6084], - [5.7, -69.5732], - ] -) - -expect_inv = np.array( - [ - [0.1, -70], - [0.2, -70], - [0.3, -70], - [0.4, -70], - [0.5, -70], - [2.8, -70], - [2.9, -70], - [3.0, -70], - [3.1, -70], - [3.2, -70], - [3.3, -70], - [3.4, -70], - [3.5, -70], - [4.8, -70], - [4.9, -70], - [5.0, -70], - [5.1, -69.9974], - [5.2, -69.9899], - [5.3, -69.9781], - [5.4, -69.9624], - [5.5, -69.9434], - [5.6, -69.9213], - [5.7, -69.8967], - [5.8, -69.8699], - [5.9, -69.8411], - [6.0, -69.8108], - ] -) diff --git a/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha_dc.py b/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha_dc.py deleted file mode 100644 index 666f339da9..0000000000 --- a/testsuite/pytests/sli2py_neurons/iaf_psc_alpha/test_iaf_psc_alpha_dc.py +++ /dev/null @@ -1,209 +0,0 @@ -# -*- coding: utf-8 -*- -# -# test_iaf_psc_alpha_dc.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 . - -import dataclasses - -import nest -import numpy as np -import pytest -import testsimulation -import testutil - - -@dataclasses.dataclass -class IAFPSCAlphaDCSimulation(testsimulation.Simulation): - amplitude: float = 1000.0 - origin: float = 0.0 - arrival: float = 3.0 - dc_delay: float = 1.0 - dc_visible: float = 3.0 - dc_duration: float = 2.0 - - def __post_init__(self): - self.dc_on = self.dc_visible - self.dc_delay - self.dc_off = self.dc_on + self.dc_duration - - def setup(self): - super().setup() - n1 = self.neuron = nest.Create("iaf_psc_alpha") - dc = self.dc_generator = nest.Create("dc_generator") - dc.amplitude = self.amplitude - vm = self.voltmeter = nest.Create("voltmeter") - vm.interval = self.resolution - nest.Connect(vm, n1) - - -@pytest.mark.parametrize("weight", [1.0]) -@testutil.use_simulation(IAFPSCAlphaDCSimulation) -class TestIAFPSCAlphaDC: - @pytest.mark.parametrize("delay", [0.1]) - def test_dc(self, simulation): - simulation.setup() - - dc_gen_spec = {"delay": simulation.delay, "weight": simulation.weight} - nest.Connect(simulation.dc_generator, simulation.neuron, syn_spec=dc_gen_spec) - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expect_default) - assert actual == expected - - @pytest.mark.parametrize("resolution,delay", [(0.1, 0.1), (0.2, 0.2), (0.5, 0.5), (1.0, 1.0)]) - def test_dc_aligned(self, simulation): - simulation.setup() - - simulation.dc_generator.set( - amplitude=simulation.amplitude, - origin=simulation.origin, - start=simulation.arrival - simulation.resolution, - ) - nest.Connect( - simulation.dc_generator, - simulation.neuron, - syn_spec={"delay": simulation.delay}, - ) - - results = simulation.simulate() - - actual, expected = testutil.get_comparable_timesamples(results, expect_aligned) - assert actual == expected - - @pytest.mark.parametrize("resolution,delay", [(0.1, 0.1), (0.2, 0.2), (0.5, 0.5), (1.0, 1.0)]) - def test_dc_aligned_auto(self, simulation): - simulation.setup() - - simulation.dc_generator.set( - amplitude=simulation.amplitude, - origin=simulation.origin, - start=simulation.dc_on, - ) - dc_gen_spec = {"delay": simulation.dc_delay, "weight": simulation.weight} - nest.Connect(simulation.dc_generator, simulation.neuron, syn_spec=dc_gen_spec) - - results = simulation.simulate() - actual, expected = testutil.get_comparable_timesamples(results, expect_aligned) - assert actual == expected - - @pytest.mark.parametrize("resolution,delay", [(0.1, 0.1), (0.2, 0.2), (0.5, 0.5), (1.0, 1.0)]) - @pytest.mark.parametrize("duration", [10.0]) - def test_dc_aligned_stop(self, simulation): - simulation.setup() - - simulation.dc_generator.set( - amplitude=simulation.amplitude, - origin=simulation.origin, - start=simulation.dc_on, - stop=simulation.dc_off, - ) - dc_gen_spec = {"delay": simulation.dc_delay, "weight": simulation.weight} - nest.Connect(simulation.dc_generator, simulation.neuron, syn_spec=dc_gen_spec) - - results = simulation.simulate() - actual, expected = testutil.get_comparable_timesamples(results, expect_stop) - assert actual == expected - - -expect_default = np.array( - [ - [0.1, -70], - [0.2, -70], - [0.3, -69.602], - [0.4, -69.2079], - [0.5, -68.8178], - [0.6, -68.4316], - [0.7, -68.0492], - [0.8, -67.6706], - [0.9, -67.2958], - [1.0, -66.9247], - [1.1, -66.5572], - [1.2, -66.1935], - [1.3, -65.8334], - [1.4, -65.4768], - [1.5, -65.1238], - [1.6, -64.7743], - ] -) - - -expect_aligned = np.array( - [ - [2.5, -70], - [2.6, -70], - [2.7, -70], - [2.8, -70], - [2.9, -70], - [3.0, -70], - [3.1, -69.602], - [3.2, -69.2079], - [3.3, -68.8178], - [3.4, -68.4316], - [3.5, -68.0492], - [3.6, -67.6706], - [3.7, -67.2958], - [3.8, -66.9247], - [3.9, -66.5572], - [4.0, -66.1935], - [4.1, -65.8334], - [4.2, -65.4768], - ] -) - - -expect_stop = np.array( - [ - [2.5, -70], - [2.6, -70], - [2.7, -70], - [2.8, -70], - [2.9, -70], - [3.0, -70], - [3.1, -69.602], - [3.2, -69.2079], - [3.3, -68.8178], - [3.4, -68.4316], - [3.5, -68.0492], - [3.6, -67.6706], - [3.7, -67.2958], - [3.8, -66.9247], - [3.9, -66.5572], - [4.0, -66.1935], - [4.1, -65.8334], - [4.2, -65.4768], - [4.3, -65.1238], - [4.4, -64.7743], - [4.5, -64.4283], - [4.6, -64.0858], - [4.7, -63.7466], - [4.8, -63.4108], - [4.9, -63.0784], - [5.0, -62.7492], - [5.1, -62.8214], - [5.2, -62.8928], - [5.3, -62.9635], - [5.4, -63.0335], - [5.5, -63.1029], - [5.6, -63.1715], - [5.7, -63.2394], - [5.8, -63.3067], - [5.9, -63.3733], - [6.0, -63.4392], - ] -) diff --git a/testsuite/pytests/sli2py_neurons/test_iaf_ps_dc_accuracy.py b/testsuite/pytests/sli2py_neurons/test_iaf_ps_dc_accuracy.py new file mode 100644 index 0000000000..904cf55dd8 --- /dev/null +++ b/testsuite/pytests/sli2py_neurons/test_iaf_ps_dc_accuracy.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +# +# test_iaf_ps_dc_accuracy.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 . + +import math + +import nest +import pytest + + +@pytest.mark.parametrize( + "model", + [ + "iaf_psc_alpha_ps", + "iaf_psc_delta_ps", + "iaf_psc_exp_ps", + "iaf_psc_exp_ps_lossless", + ], +) +@pytest.mark.parametrize("resolution", [2**i for i in range(0, -14, -1)]) +@pytest.mark.parametrize("duration, tolerance", [(5.0, 1e-13), (500.0, 1e-9)]) +def test_iaf_ps_dc_accuracy(model, resolution, duration, tolerance): + """ + A DC current is injected for a finite duration. The membrane potential at + the end of the simulated interval is compared to the theoretical value for + different computation step sizes. + + Computation step sizes are specified as base 2 values. + + Two different intervals are tested. At the end of the first interval the membrane + potential still steeply increases. At the end of the second, the membrane + potential has within double precision already reached the limit for large t. + + The high accuracy of the neuron models is achieved by the use of Exact Integration [1] + and an appropriate arrangement of the terms [2]. For small computation step sizes the + accuracy at large simulation time decreases because of the accumulation of errors. + + Reference output is documented at the end of the script. + + Individual simulation results can be inspected by uncommented the call + to function print_details. + """ + nest.ResetKernel() + nest.set(tics_per_ms=2**14, resolution=resolution) + + params = { + "E_L": 0.0, # resting potential in mV + "V_m": 0.0, # initial membrane potential in mV + "V_th": 2000.0, # spike threshold in mV + "I_e": 1000.0, # DC current in pA + "tau_m": 10.0, # membrane time constant in ms + "C_m": 250.0, # membrane capacity in pF + } + + neuron = nest.Create(model, params=params) + nest.Simulate(duration) + + expected_vm = params["I_e"] * params["tau_m"] / params["C_m"] * (1.0 - math.exp(-duration / params["tau_m"])) + + assert neuron.V_m == pytest.approx(expected_vm, abs=tolerance) + + +@pytest.mark.parametrize( + "model", + [ + "iaf_psc_alpha_ps", + "iaf_psc_delta_ps", + "iaf_psc_exp_ps", + "iaf_psc_exp_ps_lossless", + ], +) +@pytest.mark.parametrize("resolution", [2**i for i in range(0, -14, -1)]) +@pytest.mark.parametrize("duration, tolerance", [(5.0, 1e-13)]) +def test_iaf_ps_dc_t_accuracy(model, resolution, duration, tolerance): + """ + A DC current is injected for a finite duration. The time of the first + spike is compared to the theoretical value for different computation + step sizes. + + Computation step sizes are specified as base 2 values. + + The high accuracy of the neuron models is achieved by the use of + Exact Integration [1] and an appropriate arrangement of the terms + [2]. For small computation step sizes the accuracy at large + simulation time decreases because of the accumulation of errors. + + The expected output is documented at the end of the script. + Individual simulation results can be inspected by uncommented the + call to function print_details. + """ + nest.ResetKernel() + nest.set(tics_per_ms=2**14, resolution=resolution) + + params = { + "E_L": 0.0, # resting potential in mV + "V_m": 0.0, # initial membrane potential in mV + "V_th": 15.0, # spike threshold in mV + "I_e": 1000.0, # DC current in pA + "tau_m": 10.0, # membrane time constant in ms + "C_m": 250.0, # membrane capacity in pF + } + + neuron = nest.Create(model, params=params) + spike_recorder = nest.Create("spike_recorder") + nest.Connect(neuron, spike_recorder) + + nest.Simulate(duration) + + spike_times = spike_recorder.get("events", "times") + assert len(spike_times) == 1, "Neuron did not spike exactly once." + + t_spike = spike_times[0] + expected_t = -params["tau_m"] * math.log(1.0 - (params["C_m"] * params["V_th"]) / (params["tau_m"] * params["I_e"])) + + assert t_spike == pytest.approx(expected_t, abs=tolerance) diff --git a/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha.py b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha.py new file mode 100644 index 0000000000..0c75fff614 --- /dev/null +++ b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha.py @@ -0,0 +1,332 @@ +# -*- coding: utf-8 -*- +# +# test_iaf_psc_alpha.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 . + +import nest +import numpy as np +import pytest +import testutil + + +def test_iaf_psc_alpha_basic(): + """ + An overall test of the iaf_psc_alpha 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. + """ + nest.ResetKernel() + + # Simulation variables + duration = 8.0 + resolution = 0.1 + delay = 1.0 + + # Create neuron model + neuron = nest.Create("iaf_psc_alpha") + + # Create devices + voltmeter = nest.Create("voltmeter") + dc_generator = nest.Create("dc_generator") + + voltmeter.interval = resolution + dc_generator.amplitude = 1000.0 + + # Connect devices + nest.Connect(voltmeter, neuron, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(dc_generator, neuron, syn_spec={"weight": 1.0, "delay": resolution}) + + # Simulate + nest.Simulate(duration) + + # Collect spikes from spike_recorder + voltage = np.column_stack((voltmeter.events["times"], voltmeter.events["V_m"])) + + # Compare actual and expected spike times + actual, expected = testutil.get_comparable_timesamples( + resolution, + voltage, + np.array( + [ + [0.1, -70], + [0.2, -70], + [0.3, -69.602], + [0.4, -69.2079], + [0.5, -68.8178], + [0.6, -68.4316], + [0.7, -68.0492], + [0.8, -67.6706], + [0.9, -67.2958], + [1.0, -66.9247], + [4.5, -56.0204], + [4.6, -55.7615], + [4.7, -55.5051], + [4.8, -55.2513], + [4.9, -55.0001], + [5.0, -70], + [5.1, -70], + [5.2, -70], + [5.3, -70], + [5.4, -70], + [5.5, -70], + [5.6, -70], + [5.7, -70], + [5.8, -70], + [5.9, -70], + [6.0, -70], + [6.1, -70], + [6.2, -70], + [6.3, -70], + [6.4, -70], + [6.5, -70], + [6.6, -70], + [6.7, -70], + [6.8, -70], + [6.9, -70], + [7.0, -70], + [7.1, -69.602], + [7.2, -69.2079], + [7.3, -68.8178], + [7.4, -68.4316], + [7.5, -68.0492], + [7.6, -67.6706], + [7.7, -67.2958], + [7.8, -66.9247], + [7.9, -66.5572], + ] + ), + ) + + # Assert approximate equality, based on number of digits in reference data + np.testing.assert_allclose(actual, expected, rtol=1e-6) + + +def test_iaf_psc_alpha_i0(): + """ + Test of a specific feature of the iaf_psc_alpha + model. It is tested whether an internal DC current that is present + from the time of neuron initialization, correctly affects the membrane + potential. + + This is probably the simplest setup in which we can study how the + dynamics develops from an initial condition. + + When the DC current is supplied by a device external to the neuron + the situation is more complex because additional delays are introduced. + """ + nest.ResetKernel() + + # Simulation variables + resolution = 0.1 + delay = resolution + duration = 8.0 + + # Create neuron and devices + neuron = nest.Create("iaf_psc_alpha", params={"I_e": 1000.0}) + voltmeter = nest.Create("voltmeter") + voltmeter.interval = resolution + spike_recorder = nest.Create("spike_recorder") + + # Connect devices + nest.Connect(voltmeter, neuron, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(neuron, spike_recorder, syn_spec={"weight": 1.0, "delay": delay}) + + # Simulate + nest.Simulate(duration) + + # Get voltmeter output + results = np.column_stack((voltmeter.events["times"], voltmeter.events["V_m"])) + + # Get spike times + spikes = np.column_stack( + ( + spike_recorder.events["senders"], + spike_recorder.events["times"], + ) + ) + + # Compare voltmeter output to expected + actual, expected = testutil.get_comparable_timesamples( + resolution, + results, + np.array( + [ + [0.1, -69.602], + [0.2, -69.2079], + [0.3, -68.8178], + [0.4, -68.4316], + [0.5, -68.0492], + [4.3, -56.0204], + [4.4, -55.7615], + [4.5, -55.5051], + [4.6, -55.2513], + [4.7, -55.0001], + [4.8, -70], + [4.9, -70], + [5.0, -70], + ] + ), + ) + np.testing.assert_allclose(actual, expected, rtol=1e-6) + + # Compare spike times + assert spikes == pytest.approx( + np.array( + [ + [1, 4.8], + ] + ) + ) + + +@pytest.mark.parametrize("resolution", [0.1, 0.2, 0.5, 1.0]) +def test_iaf_psc_alpha_i0_refractory(resolution): + """ + Test a specific feature of the iaf_psc_alpha model. + + It is tested whether the voltage traces of simulations + carried out at different resolutions (computation step sizes) are well + aligned and identical when the neuron recovers from refractoriness. + + In grid based simulation a prerequisite is that the spike is reported at + a grid position shared by all the resolutions compared. + + Here, we compare resolutions 0.1, 0.2, 0.5, and 1.0 ms. Therefore, the + internal DC current is adjusted such (1450.0 pA) that the spike is + reported at time 3.0 ms, corresponding to computation step 30, 15, 6, + and 3, respectively. + + The results are consistent with those of iaf_psc_alpha_ps capable of + handling off-grid spike timing when the interpolation order is set to + 0. + """ + # Simulation variables + delay = resolution + duration = 8.0 + + nest.ResetKernel() + nest.resolution = resolution + + neuron = nest.Create("iaf_psc_alpha", params={"I_e": 1450.0}) + voltmeter = nest.Create("voltmeter", params={"interval": resolution}) + spike_recorder = nest.Create("spike_recorder") + + nest.Connect(voltmeter, neuron, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(neuron, spike_recorder, syn_spec={"weight": 1.0, "delay": delay}) + + nest.Simulate(duration) + + times = voltmeter.events["times"] + V_m = voltmeter.events["V_m"] + results = np.column_stack((times, V_m)) + + actual, expected = testutil.get_comparable_timesamples( + resolution, + results, + np.array( + [ + [0.1, -69.4229], + [0.2, -68.8515], + [0.3, -68.2858], + [0.4, -67.7258], + [0.5, -67.1713], + [0.6, -66.6223], + [0.7, -66.0788], + [0.8, -65.5407], + [0.9, -65.008], + [1.0, -64.4806], + [1.1, -63.9584], + [1.2, -63.4414], + [1.3, -62.9295], + [1.4, -62.4228], + [1.5, -61.9211], + [1.6, -61.4243], + [1.7, -60.9326], + [1.8, -60.4457], + [1.9, -59.9636], + [2.0, -59.4864], + [2.1, -59.0139], + [2.2, -58.5461], + [2.3, -58.0829], + [2.4, -57.6244], + [2.5, -57.1704], + [2.6, -56.721], + [2.7, -56.276], + [2.8, -55.8355], + [2.9, -55.3993], + [3.0, -70], + [3.1, -70], + [3.2, -70], + [3.3, -70], + [3.4, -70], + [3.5, -70], + [3.6, -70], + [3.7, -70], + [3.8, -70], + [3.9, -70], + [4.0, -70], + [4.1, -70], + [4.2, -70], + [4.3, -70], + [4.4, -70], + [4.5, -70], + [4.6, -70], + [4.7, -70], + [4.8, -70], + [4.9, -70], + [5.0, -70], + [5.1, -69.4229], + [5.2, -68.8515], + [5.3, -68.2858], + [5.4, -67.7258], + [5.5, -67.1713], + [5.6, -66.6223], + [5.7, -66.0788], + [5.8, -65.5407], + [5.9, -65.008], + [6.0, -64.4806], + [6.1, -63.9584], + [6.2, -63.4414], + [6.3, -62.9295], + [6.4, -62.4228], + [6.5, -61.9211], + [6.6, -61.4243], + [6.7, -60.9326], + [6.8, -60.4457], + [6.9, -59.9636], + ] + ), + ) + np.testing.assert_allclose(actual, expected, rtol=1e-6) diff --git a/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_1to2.py b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_1to2.py new file mode 100644 index 0000000000..43dcaecfce --- /dev/null +++ b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_1to2.py @@ -0,0 +1,183 @@ +# -*- coding: utf-8 -*- +# +# test_iaf_psc_alpha_1to2.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 . + + +import nest +import numpy as np +import pandas as pd +import pytest + + +def Vm_theory(t, nrn_params, syn_weight): + """ + Returns the value of the membrane potential at time t, assuming + alpha-shaped post-synaptic currents and an incoming spike at t=0. + """ + + assert (t >= 0).all() + + C_m = nrn_params["C_m"] + tau_m = nrn_params["tau_m"] + tau_syn = nrn_params["tau_syn_ex"] + E_L = nrn_params["E_L"] + + prefactor = syn_weight * np.exp(1) / (tau_syn * C_m) + term1 = (np.exp(-t / tau_m) - np.exp(-t / tau_syn)) / (1 / tau_syn - 1 / tau_m) ** 2 + term2 = t * np.exp(-t / tau_syn) / (1 / tau_syn - 1 / tau_m) + return E_L + prefactor * (term1 - term2) + + +@pytest.mark.parametrize( + "resolution, min_delay", [(0.1, 1.0), (0.2, 1.0), (0.5, 1.0), (1.0, 1.0), (0.1, 0.1), (0.1, 0.2), (0.1, 0.5)] +) +def test_1to2(resolution, min_delay): + """ + Checks the spike interaction of two iaf_psc_alpha model neurons. + + In order to obtain identical results for different computation step + sizes h, the SLI script needs to be independent of h. This is + achieved by specifying all time parameters in milliseconds (ms). In + particular the time of spike emission and the synaptic delay need to + be integer multiples of the computation step sizes to be + tested. test_iaf_dc_aligned_delay demonstrates the strategy for the + case of DC current input. + + A DC current in the pre-synaptic neuron is adjusted to cause a spike + at a grid position (t=3.0 ms) joined by all computation step sizes to + be tested. + + Note that in a neuron model where synaptic events are modeled by a + truncated exponential the effect of the incoming spike would be + visible at the time of impact (here, t=4.0 ms). This is because the + initial condition for the postsynaptic potential (PSP) has a + non-zero voltage component. For PSPs with finite rise time the + situation is different. In this case the voltage component of the + initial condition is zero (see documentation of + test_iaf_psp). Therefore, at the time of impact the PSP is only + visible in other components of the state vector. + + See end of file for commented table of expected output. + """ + + nest.ResetKernel() + nest.set(resolution=resolution, min_delay=min_delay, max_delay=1) + + # Parameters for network and some pre-calculated results + I_ext = 1450 + time_to_threshold = 3 # for iaf_psc_alpha with default params driven by I_ext + delay = 1 + weight = 100 + + # Create neurons + n1, n2 = nest.Create("iaf_psc_alpha", 2) + n1.I_e = I_ext + + # Spike recorder for neuron 1 + sr = nest.Create("spike_recorder", params={"time_in_steps": True}) + nest.Connect(n1, sr) + + # Voltmeter for n2 + vm = nest.Create("voltmeter", params={"interval": resolution, "time_in_steps": True}) + nest.Connect(vm, n2) + + nest.Connect(n1, n2, syn_spec={"weight": weight, "delay": delay}) + + # Expected spike times of first neuron for default parameters + t_expected = [time_to_threshold, 2 * time_to_threshold + n1.t_ref] + + # Run the simulation until time of second spike of n1, so it does not affect n2 + nest.Simulate(t_expected[1]) + + spikes = sr.events["times"] + voltages = pd.DataFrame.from_records(vm.events) + + assert list(spikes) == [round(t / resolution) for t in t_expected] + + # Effect of first spike on membrane potential of second. Begins with arrival of first spike of n1 at n2 + t_arrive = t_expected[0] + delay + t_arrive_steps = round(t_arrive / resolution) + + # Must have resting membrane potential up to spike arrival at n2 + v_pre_spike = voltages.loc[voltages.times <= t_arrive_steps].V_m + np.testing.assert_array_equal(v_pre_spike, n2.E_L) + + # Compare membrane potential from spike arrival onward to analytical solution + t_post_spike = voltages.loc[voltages.times >= t_arrive_steps].times * resolution - t_arrive + v_post_spike = voltages.loc[voltages.times >= t_arrive_steps].V_m + + v_expected = Vm_theory(t_post_spike, n2.get(), weight) + + np.testing.assert_allclose(v_post_spike, v_expected, rtol=1e-12) + + +# ------------------------------------------------------------------------------------------ +# +# Expected output of this simulation +# +# The output sent 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. +# +# +# h= (in ms) +# [ 0.1 0.2 0.5 1.0] +# +# time voltage +# [ +# ... +# [ 25 5 -70]% <-- Voltage trace of the postsynaptic neuron +# [ 26 13 -70]% (neuron2), at rest until a spike arrives. +# [ 27 -70] +# [ 28 14 -70] +# [ 29 -70] +# [ 30 15 6 3 -70] +# 1 30 % <-- The pre-synaptic neuron (neuron1) emits a +# [ 31 -70]% spike at t=3.0 ms. +# [ 32 16 -70] +# [ 33 -70] +# [ 34 17 -70] +# [ 35 7 -70]% <-- Synaptic delay of 1.0 ms. +# [ 36 18 -70] +# [ 37 -70] +# [ 38 19 -70] +# [ 39 -70] +# [ 40 20 8 4 -70]% <----------- Spike arrives at the postsynaptic neuron +# [ 41 -69.9974]% <- (neuron2) and changes the state vector of +# [ 42 21 -69.9899]% | the neuron, not visible in voltage because +# [ 43 -69.9781]% | voltage of PSP initial condition is 0. +# [ 44 22 -69.9624]% | +# [ 45 9 -69.9434]% - Arbitrarily close to the time of impact +# [ 46 23 -69.9213]% (t=4.0 ms) the effect of the spike (PSP) +# [ 47 -69.8967]% is visible in the voltage trace. +# [ 48 24 -69.8699] +# [ 49 -69.8411] +# [ 50 25 10 5 -69.8108] +# [ 51 -69.779 ] +# [ 52 26 -69.7463]% <--- The voltage trace is independent +# [ 53 -69.7126]% of the computation step size h. +# [ 54 27 -69.6783]% Larger step sizes only have fewer +# [ 55 11 -69.6435]% sample points. +# [ 56 28 -69.6084] +# [ 57 -69.5732] +# ... +# ] +# +# ------------------------------------------------------------------------------------------ diff --git a/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_fudge.py b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_fudge.py new file mode 100644 index 0000000000..e9db172a93 --- /dev/null +++ b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_fudge.py @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +# +# test_iaf_psc_alpha_fudge.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 . + +import math + +import nest +import numpy as np +import pytest +import testutil +from scipy.special import lambertw + + +def test_iaf_psc_alpha_fudge(): + """ + The peak time of the postsynaptic potential (PSP) is calculated using + the LambertW function. The theoretical peak voltage amplitude for a + postsynaptic current of amplitude 1pA is then used to adjust the + synaptic weight such that a PSP of amplitude 1mV is generated. The + success of this approach is verified by comparing the theoretical + value with the result of a simulation where a single spike is sent to + the neuron. + + The name of this test script has a historical explanation. Prior to + July 2009 the analytical expression for the peak time of the PSP was + not known to the NEST developers. Therefore the normalization factor + required to adjust the PSP amplitude was computed by root finding + outside of NEST. The factor was called "fudge" in examples and + application code. The root finding was not done in NEST because infix + mathematical notation only become available in SLI in January + 2009. The name "fudge" indicated that the origin of this value was not + obvious from the simulation scripts and usage was inherently dangerous + because a change of the time constants of the neuron model would + invalidate the value of "fudge". + """ + + # Simulation variables + resolution = 0.1 + delay = resolution + duration = 20.0 + + nest.resolution = resolution + + # Create neuron and devices + neuron = nest.Create("iaf_psc_alpha") + voltmeter = nest.Create("voltmeter", params={"interval": resolution}) + + # Connect voltmeter + nest.Connect(voltmeter, neuron) + + # Biophysical parameters + tau_m = 20.0 + tau_syn = 0.5 + C_m = 250.0 + + # Set neuron parameters + neuron.tau_m = tau_m + neuron.tau_syn_ex = tau_syn + neuron.tau_syn_in = tau_syn + neuron.C_m = C_m + + # Compute fudge factors + a = tau_m / tau_syn + b = 1.0 / tau_syn - 1.0 / tau_m + t_max = (1.0 / b) * (-lambertw(-math.exp(-1.0 / a) / a, k=-1) - 1.0 / a).real + + v_max = ( + math.exp(1) + / (tau_syn * C_m * b) + * ((math.exp(-t_max / tau_m) - math.exp(-t_max / tau_syn)) / b - t_max * math.exp(-t_max / tau_syn)) + ) + + # Create spike generator to fire once at resolution + spike_gen = nest.Create( + "spike_generator", + params={ + "precise_times": False, + "spike_times": [resolution], + }, + ) + + # Connect spike generator to neuron + nest.Connect(spike_gen, neuron, syn_spec={"weight": float(1.0 / v_max), "delay": delay}) + + # Simulate + nest.Simulate(duration) + + # Extract membrane potential trace + volt_data = voltmeter.events + times = volt_data["times"] + v_m = volt_data["V_m"] + results = np.column_stack((times, v_m)) + + # Find time of peak voltage and peak voltage and thus height of PSP + max_vm_ix = np.argmax(results[:, 1]) + actual_t_max = results[max_vm_ix, 0] + actual_vm_max = results[max_vm_ix, 1] + actual_psp_height = actual_vm_max - neuron.E_L + + expected_psp_height = 1 + expected_t_max = t_max + resolution + delay # spike is sent at t=resolution and arrives with delay + + assert actual_t_max == pytest.approx(expected_t_max, abs=resolution / 2) + assert actual_psp_height == pytest.approx(expected_psp_height, abs=1e-4) diff --git a/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_mindelay.py b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_mindelay.py new file mode 100644 index 0000000000..daecafd800 --- /dev/null +++ b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_mindelay.py @@ -0,0 +1,163 @@ +# -*- coding: utf-8 -*- +# +# test_iaf_psc_alpha_mindelay.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 . + + +import math + +import nest +import numpy as np +import pytest +import testutil + + +@pytest.mark.parametrize("min_delay", [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1.0, 2.0]) +def test_iaf_psc_alpha_mindelay_create(min_delay): + """ + Tests automatic adjustment of min_delay. + + The simulation is run with a range of different min_delays. All + should give identical results. This is achieved by sampling the + membrane potential at a fixed interval. + """ + nest.ResetKernel() + nest.resolution = 0.1 + + # Simulation variables + delay = 2.0 + duration = 10.5 + + neuron = nest.Create("iaf_psc_alpha") + voltmeter = nest.Create("voltmeter", params={"interval": 0.1}) + spike_recorder = nest.Create("spike_recorder") + dc_generator = nest.Create("dc_generator", params={"amplitude": 1000.0}) + + nest.Connect(voltmeter, neuron, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(neuron, spike_recorder, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(dc_generator, neuron, syn_spec={"weight": 1.0, "delay": delay}) + + # Force `min_delay` by connecting two throwaway neurons + throwaway = nest.Create("iaf_psc_alpha", 2) + nest.Connect(throwaway[0], throwaway[1], syn_spec={"weight": 1.0, "delay": min_delay}) + + nest.Simulate(duration) + + results = np.column_stack((voltmeter.events["times"], voltmeter.events["V_m"])) + + actual, expected = testutil.get_comparable_timesamples(nest.resolution, results, expected_mindelay) + np.testing.assert_allclose(actual, expected, rtol=1e-7) + + +@pytest.mark.parametrize("min_delay", [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1.0, 2.0]) +def test_iaf_psc_alpha_mindelay_set(min_delay): + """ + Tests explicitly setting min_delay. + + The simulation is run with a range of different min_delays. All + should give identical results. This is achieved by sampling the + membrane potential at a fixed interval. + """ + nest.ResetKernel() + nest.resolution = 0.1 + + # Simulation variables + delay = 2.0 + duration = 10.5 + + # Set up test min_delay conditions + nest.set(min_delay=min_delay, max_delay=delay) + nest.SetDefaults("static_synapse", {"delay": delay}) + + neuron = nest.Create("iaf_psc_alpha") + voltmeter = nest.Create("voltmeter", params={"interval": 0.1}) + spike_recorder = nest.Create("spike_recorder") + dc_generator = nest.Create("dc_generator", params={"amplitude": 1000.0}) + + nest.Connect(voltmeter, neuron, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(neuron, spike_recorder, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(dc_generator, neuron, syn_spec={"weight": 1.0, "delay": delay}) + + nest.Simulate(duration) + + results = np.column_stack((voltmeter.events["times"], voltmeter.events["V_m"])) + + actual, expected = testutil.get_comparable_timesamples(nest.resolution, results, expected_mindelay) + np.testing.assert_allclose(actual, expected, rtol=1e-7) + + +@pytest.mark.parametrize("min_delay", [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1.0, 2.0]) +def test_iaf_psc_alpha_mindelay_simblocks(min_delay): + """ + Tests explicitly setting min_delay across 21 simulation calls. + + The simulation is run with a range of different min_delays. All + should give identical results. This is achieved by sampling the + membrane potential at a fixed interval. + """ + nest.ResetKernel() + nest.resolution = 0.1 + + # Simulation variables + delay = 2.0 + + # Set up test min_delay conditions + nest.set(min_delay=min_delay, max_delay=delay) + nest.SetDefaults("static_synapse", {"delay": delay}) + + neuron = nest.Create("iaf_psc_alpha") + voltmeter = nest.Create("voltmeter", params={"interval": 0.1}) + spike_recorder = nest.Create("spike_recorder") + dc_generator = nest.Create("dc_generator", params={"amplitude": 1000.0}) + + nest.Connect(voltmeter, neuron, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(neuron, spike_recorder, syn_spec={"weight": 1.0, "delay": delay}) + nest.Connect(dc_generator, neuron, syn_spec={"weight": 1.0, "delay": delay}) + + for _ in range(21): + nest.Simulate(0.5) + + results = np.column_stack((voltmeter.events["times"], voltmeter.events["V_m"])) + + actual, expected = testutil.get_comparable_timesamples(nest.resolution, results, expected_mindelay) + + assert len(actual) == len(expected_mindelay) + np.testing.assert_allclose(actual, expected, rtol=1e-7) + + +def test_kernel_precision(): + nest.ResetKernel() + nest.set(tics_per_ms=2**14, resolution=2**0) + assert math.frexp(nest.ms_per_tic) == (0.5, -13) + + +expected_mindelay = np.array( + [ + [1.000000e00, -7.000000e01], + [2.000000e00, -7.000000e01], + [3.000000e00, -6.655725e01], + [4.000000e00, -6.307837e01], + [5.000000e00, -5.993054e01], + [6.000000e00, -5.708227e01], + [7.000000e00, -7.000000e01], + [8.000000e00, -7.000000e01], + [9.000000e00, -6.960199e01], + [1.000000e01, -6.583337e01], + ] +) diff --git a/testsuite/pytests/utilities/testsimulation.py b/testsuite/pytests/utilities/testsimulation.py deleted file mode 100644 index eb56136e7f..0000000000 --- a/testsuite/pytests/utilities/testsimulation.py +++ /dev/null @@ -1,51 +0,0 @@ -# -*- coding: utf-8 -*- -# -# testsimulation.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 . - -import dataclasses - -import nest -import numpy as np -import testutil - - -@dataclasses.dataclass -class Simulation: - local_num_threads: int = 1 - resolution: float = 0.1 - duration: float = 8.0 - weight: float = 100.0 - delay: float = 1.0 - - def setup(self): - pass - - def simulate(self): - nest.Simulate(self.duration) - if hasattr(self, "voltmeter"): - return np.column_stack((self.voltmeter.events["times"], self.voltmeter.events["V_m"])) - - def run(self): - self.setup() - return self.simulate() - - def __init_subclass__(cls, **kwargs): - super().__init_subclass__(**kwargs) - testutil.create_dataclass_fixtures(cls) diff --git a/testsuite/pytests/utilities/testutil.py b/testsuite/pytests/utilities/testutil.py index 5f8bed0cfd..63e44178e5 100644 --- a/testsuite/pytests/utilities/testutil.py +++ b/testsuite/pytests/utilities/testutil.py @@ -19,17 +19,13 @@ # You should have received a copy of the GNU General Public License # along with NEST. If not, see . -import dataclasses import sys -import numpy as np +import nest +import pandas as pd import pytest -def parameter_fixture(name, default_factory=lambda: None): - return pytest.fixture(autouse=True, name=name)(lambda request: getattr(request, "param", default_factory())) - - def dict_is_subset_of(small, big): """ Return true if dict `small` is subset of dict `big`. @@ -47,69 +43,26 @@ def dict_is_subset_of(small, big): return big == big | small -def isin_approx(A, B, tol=1e-06): - A = np.asarray(A) - B = np.asarray(B) - - Bs = np.sort(B) # skip if already sorted - idx = np.searchsorted(Bs, A) - - linvalid_mask = idx == len(B) - idx[linvalid_mask] = len(B) - 1 - lval = Bs[idx] - A - lval[linvalid_mask] *= -1 - - rinvalid_mask = idx == 0 - idx1 = idx - 1 - idx1[rinvalid_mask] = 0 - rval = A - Bs[idx1] - rval[rinvalid_mask] *= -1 - return np.minimum(lval, rval) <= tol - - -def get_comparable_timesamples(actual, expected): - simulated_points = isin_approx(actual[:, 0], expected[:, 0]) - expected_points = isin_approx(expected[:, 0], actual[:, 0]) - assert len(actual[simulated_points]) > 0, "The recorded data did not contain any relevant timesamples" - return actual[simulated_points], pytest.approx(expected[expected_points]) - - -def create_dataclass_fixtures(cls, module_name=None): - for field, type_ in getattr(cls, "__annotations__", {}).items(): - if isinstance(field, dataclasses.Field): - name = field.name - if field.default_factory is not dataclasses.MISSING: - default = field.default_factory - else: - - def default(d=field.default): - return d +def get_comparable_timesamples(resolution, actual, expected): + """ + Return result of inner join on time stamps for actual and expected given resolution. - else: - name = field - attr = getattr(cls, field) - # We may be receiving a mixture of literal default values, field defaults, - # and field default factories. - if isinstance(attr, dataclasses.Field): - if attr.default_factory is not dataclasses.MISSING: - default = attr.default_factory - else: + `actual` and `expected` must be arrays with columns representing times (in ms) and values. + Times will be converted to steps given the resolution and the data will then be inner- + joined on the steps, i.e., rows with equal steps will be extracted. - def default(d=attr.default): - return d + Returns two one-dimensional arrays containing the values at the joint points from + actual and expected, respectively. + """ - else: + tics = nest.tics_per_ms - def default(d=attr): - return d + actual = pd.DataFrame(actual, columns=["t", "val_a"]) + expected = pd.DataFrame(expected, columns=["t", "val_e"]) - setattr( - sys.modules[module_name or cls.__module__], - name, - parameter_fixture(name, default), - ) + actual["tics"] = (actual.t * tics).round().astype(int) + expected["tics"] = (expected.t * tics).round().astype(int) + common = pd.merge(actual, expected, how="inner", on="tics") -def use_simulation(cls): - # If `mark` receives one argument that is a class, it decorates that arg. - return pytest.mark.simulation(cls, "") + return common.val_a.values, common.val_e.values