From 146164efad0335776b117ee58aac1e97dab8bc45 Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Sat, 3 Dec 2016 14:58:48 +0100 Subject: [PATCH 01/11] corrected refractory implementation --- models/hh_cond_exp_traub.cpp | 28 ++++++++++++++++++++-------- models/hh_cond_exp_traub.h | 4 +++- 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/models/hh_cond_exp_traub.cpp b/models/hh_cond_exp_traub.cpp index 6c56f83066..c88e16af72 100644 --- a/models/hh_cond_exp_traub.cpp +++ b/models/hh_cond_exp_traub.cpp @@ -136,18 +136,19 @@ hh_cond_exp_traub_dynamics( double, const double y[], double f[], void* pnode ) * ---------------------------------------------------------------- */ nest::hh_cond_exp_traub::Parameters_::Parameters_() - : g_Na( 20000.0 ) // Sodium Conductance (nS) - , g_K( 6000.0 ) // K Conductance (nS) - , g_L( 10.0 ) // Leak Conductance (nS) - , C_m( 200.0 ) // Membrane Capacitance (pF) - , E_Na( 50.0 ) // Reversal potentials (mV) + : g_Na( 20000.0 ) // Sodium Conductance (nS) + , g_K( 6000.0 ) // K Conductance (nS) + , g_L( 10.0 ) // Leak Conductance (nS) + , C_m( 200.0 ) // Membrane Capacitance (pF) + , E_Na( 50.0 ) // Reversal potentials (mV) , E_K( -90.0 ) , E_L( -60.0 ) - , V_T( -63.0 ) // adjusts threshold to around -50 mV + , V_T( -63.0 ) // adjusts threshold to around -50 mV , E_ex( 0.0 ) , E_in( -80.0 ) , tau_synE( 5.0 ) // Synaptic Time Constant Excitatory Synapse (ms) , tau_synI( 10.0 ) // Synaptic Time Constant Excitatory Synapse (ms) + , t_ref_( 0.0 ) // Refractory time in ms , I_e( 0.0 ) // Stimulus Current (pA) { } @@ -213,6 +214,7 @@ nest::hh_cond_exp_traub::Parameters_::get( DictionaryDatum& d ) const def< double >( d, names::E_in, E_in ); def< double >( d, names::tau_syn_ex, tau_synE ); def< double >( d, names::tau_syn_in, tau_synI ); + def< double >( d, names::t_ref, t_ref_ ); def< double >( d, names::I_e, I_e ); } @@ -231,13 +233,23 @@ nest::hh_cond_exp_traub::Parameters_::set( const DictionaryDatum& d ) updateValue< double >( d, names::E_in, E_in ); updateValue< double >( d, names::tau_syn_ex, tau_synE ); updateValue< double >( d, names::tau_syn_in, tau_synI ); + updateValue< double >( d, names::t_ref, t_ref_ ); updateValue< double >( d, names::I_e, I_e ); if ( C_m <= 0 ) + { throw BadProperty( "Capacitance must be strictly positive." ); + } if ( tau_synE <= 0 || tau_synI <= 0 ) + { throw BadProperty( "All time constants must be strictly positive." ); + } + + if ( t_ref_ < 0 ) + { + throw BadProperty( "Refractory time cannot be negative." ); + } } void @@ -368,7 +380,7 @@ nest::hh_cond_exp_traub::calibrate() { // ensures initialization in case mm connected after Simulate B_.logger_.init(); - V_.RefractoryCounts_ = 20; + V_.refractory_counts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); V_.U_old_ = S_.y_[ State_::V_M ]; } @@ -422,7 +434,7 @@ nest::hh_cond_exp_traub::update( Time const& origin, if ( S_.y_[ State_::V_M ] >= P_.V_T + 30. && V_.U_old_ > S_.y_[ State_::V_M ] ) { - S_.r_ = V_.RefractoryCounts_; + S_.r_ = V_.refractory_counts_; set_spiketime( Time::step( origin.get_steps() + lag + 1 ) ); diff --git a/models/hh_cond_exp_traub.h b/models/hh_cond_exp_traub.h index d06663a23f..e520ecfce3 100644 --- a/models/hh_cond_exp_traub.h +++ b/models/hh_cond_exp_traub.h @@ -92,6 +92,7 @@ contains variables called y,s,r,q and \chi. function in ms. tau_syn_in double - Time constant of the inhibitory synaptic exponential function in ms. + t_ref double - Duration of refractory period in ms. E_ex double - Excitatory synaptic reversal potential in mV. E_in double - Inhibitory synaptic reversal potential in mV. E_Na double - Sodium reversal potential in mV. @@ -186,6 +187,7 @@ class hh_cond_exp_traub : public Archiving_Node double E_in; //!< Inhibitory reversal Potential in mV double tau_synE; //!< Synaptic Time Constant Excitatory Synapse in ms double tau_synI; //!< Synaptic Time Constant Inhibitory Synapse in ms + double t_ref_; //!< Refractory time in ms double I_e; //!< External Current in pA Parameters_(); @@ -235,7 +237,7 @@ class hh_cond_exp_traub : public Archiving_Node */ struct Variables_ { - int RefractoryCounts_; + int refractory_counts_; double U_old_; // for spike-detection }; From 1ec48815e23576ddbb40871934560f756f504b40 Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Mon, 12 Dec 2016 23:40:10 +0100 Subject: [PATCH 02/11] added refractory test --- pynest/nest/tests/test_refractory.py | 200 +++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) create mode 100644 pynest/nest/tests/test_refractory.py diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py new file mode 100644 index 0000000000..54f97a1c32 --- /dev/null +++ b/pynest/nest/tests/test_refractory.py @@ -0,0 +1,200 @@ +# -*- coding: utf-8 -*- +# +# test_refractory.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 unittest + +import numpy as np + +import nest + +""" +Assert that all neuronal models implement refractory period correctly. + +Details +------- +Submit the neuron to a constant excitatory current so that it spikes in the +[0, 50] ms. +A ``spike_detector`` is used to detect the time at which the neuron spikes and +a ``voltmeter`` is then used to make sure the voltage is clamped to ``V_reset`` +during exactly ``t_ref``. +""" + + +# --------------------------------------------------------------------------- # +# Models, specific parameters +# ------------------------- +# + +ignore_model = [ + "amat2_psc_exp", + "ginzburg_neuron", + "hh_cond_exp_traub", + "hh_psc_alpha", + "hh_psc_alpha_gap", + "ht_neuron", + "iaf_chs_2007", + "iaf_chxk_2008", + "iaf_cond_alpha_mc", + "iaf_psc_alpha_canon", + "iaf_psc_alpha_presc", + "iaf_psc_delta_canon", + "iaf_psc_exp_ps", + "iaf_tum_2000", + "izhikevich", + "mat2_psc_exp", + "mcculloch_pitts_neuron", + "parrot_neuron", + "parrot_neuron_ps", + "pp_pop_psc_delta", + "pp_psc_delta", + "sli_neuron", +] + +nodes = nest.Models("nodes") +# list of all neuronal models to be tested +neurons = [m for m in nodes if nest.GetDefaults(m, "element_type") == "neuron" + and m not in ignore_model] + +# additional parameters for the recorder +#~ add_vm_param = { + #~ "iaf_cond_alpha_mc": {"record_from": "V_m.s"}, +#~ } + + +# --------------------------------------------------------------------------- # +# Simulation and refractory time +# ------------------------- +# + +simtime = 100 +resolution = 0.1 +min_steps = 1 +max_steps = 200 + + +# --------------------------------------------------------------------------- # +# Test class +# ------------------------- +# + +def foreach_neuron(func): + ''' + Decorator that automatically does the test for all neurons. + ''' + def wrapper(*args, **kwargs): + self = args[0] + msd = 123456 + N_vp = nest.GetKernelStatus(['total_num_virtual_procs'])[0] + pyrngs = [np.random.RandomState(s) for s in range(msd, msd + N_vp)] + for name in neurons: + nest.ResetKernel() + nest.SetKernelStatus({ + 'resolution': resolution, 'grng_seed': msd + N_vp, + 'rng_seeds': range(msd + N_vp + 1, msd + 2 * N_vp + 1)}) + func(self, name, **kwargs) + return wrapper + + +class RefractoryTestCase(unittest.TestCase): + """ + Check the correct implementation of refractory time in all neuronal models. + """ + + def compute_reftime(self, model, sd, vm, Vr): + ''' + Compute the refractory time of the neuron. + + Parameters + ---------- + sd : tuple + GID of the spike detector. + vm : tuple + GID of the voltmeter. + Vr : double + Value of the reset potential. + + Returns + ------- + t_ref_sim : double + Value of the simulated refractory period. + ''' + spike_times = nest.GetStatus(sd, "events")[0]["times"] + times = nest.GetStatus(vm, "events")[0]["times"] + idx_max = (np.argwhere(np.isclose(times, spike_times[1]))[0] + if len(spike_times) > 1 else -1) + Vs = nest.GetStatus(vm, "events")[0]["V_m"] + # get the index at which the spike occured + idx_spike = np.argwhere(times == spike_times[0])[0] + idx_end = np.where(np.isclose(Vs[idx_spike:idx_max], Vr, 1e-6))[0][-1] + t_ref_sim = 0. + if "gif" in model or "iaf" in model: + t_ref_sim = idx_end * resolution + else: + t_ref_sim = (idx_end) * resolution + return t_ref_sim + + + @foreach_neuron + def test_refractory_time(self, model): + ''' + Check that refractory time implementation is correct. + ''' + # randomly set a refractory period + t_ref = resolution * np.random.randint(min_steps, max_steps) + # create the neuron and devices + nparams = {"t_ref": t_ref} + neuron = nest.Create(model, params=nparams) + vm_params = {"interval": resolution} + #~ vm_params.update(add_vm_param.get(model, {})) + vm = nest.Create("voltmeter", params=vm_params) + sd = nest.Create("spike_detector") + cg = nest.Create("dc_generator", params={"amplitude": 600.}) + # connect them and simulate + nest.Connect(vm, neuron) + nest.Connect(cg, neuron) + nest.Connect(neuron, sd) + + nest.Simulate(simtime) + + # get and compare t_ref + Vr = nest.GetStatus(neuron, "V_reset")[0] + t_ref_sim = self.compute_reftime(model, sd, vm, Vr) + + self.assertEqual(t_ref, t_ref_sim, '''Error in model {}: + {} != {}'''.format(model, t_ref, t_ref_sim)) + + +# --------------------------------------------------------------------------- # +# Run the comparisons +# ------------------------ +# + +def suite(): + return unittest.makeSuite(RefractoryTestCase, "test") + + +def run(): + runner = unittest.TextTestRunner(verbosity=2) + runner.run(suite()) + + +if __name__ == '__main__': + run() From c2e0238237d290c00decc7bd29d53e607b90d01a Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Tue, 13 Dec 2016 10:27:45 +0100 Subject: [PATCH 03/11] support for multicompartment iaf --- pynest/nest/tests/test_refractory.py | 29 ++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index 54f97a1c32..c9416f0c02 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -44,6 +44,9 @@ # ignore_model = [ + "aeif_cond_alpha_RK5", + "aeif_cond_alpha_multisynapse", + "aeif_cond_beta_multisynapse", "amat2_psc_exp", "ginzburg_neuron", "hh_cond_exp_traub", @@ -52,7 +55,6 @@ "ht_neuron", "iaf_chs_2007", "iaf_chxk_2008", - "iaf_cond_alpha_mc", "iaf_psc_alpha_canon", "iaf_psc_alpha_presc", "iaf_psc_delta_canon", @@ -73,10 +75,10 @@ neurons = [m for m in nodes if nest.GetDefaults(m, "element_type") == "neuron" and m not in ignore_model] -# additional parameters for the recorder -#~ add_vm_param = { - #~ "iaf_cond_alpha_mc": {"record_from": "V_m.s"}, -#~ } +# additional parameters for the connector +add_connect_param = { + "iaf_cond_alpha_mc": {"receptor_type": 7}, +} # --------------------------------------------------------------------------- # @@ -124,6 +126,8 @@ def compute_reftime(self, model, sd, vm, Vr): Parameters ---------- + model : str + Name of the neuronal model. sd : tuple GID of the spike detector. vm : tuple @@ -140,15 +144,12 @@ def compute_reftime(self, model, sd, vm, Vr): times = nest.GetStatus(vm, "events")[0]["times"] idx_max = (np.argwhere(np.isclose(times, spike_times[1]))[0] if len(spike_times) > 1 else -1) - Vs = nest.GetStatus(vm, "events")[0]["V_m"] + name_Vm = "V_m.s" if model == "iaf_cond_alpha_mc" else "V_m" + Vs = nest.GetStatus(vm, "events")[0][name_Vm] # get the index at which the spike occured idx_spike = np.argwhere(times == spike_times[0])[0] idx_end = np.where(np.isclose(Vs[idx_spike:idx_max], Vr, 1e-6))[0][-1] - t_ref_sim = 0. - if "gif" in model or "iaf" in model: - t_ref_sim = idx_end * resolution - else: - t_ref_sim = (idx_end) * resolution + t_ref_sim = idx_end * resolution return t_ref_sim @@ -162,14 +163,14 @@ def test_refractory_time(self, model): # create the neuron and devices nparams = {"t_ref": t_ref} neuron = nest.Create(model, params=nparams) - vm_params = {"interval": resolution} - #~ vm_params.update(add_vm_param.get(model, {})) + name_Vm = "V_m.s" if model == "iaf_cond_alpha_mc" else "V_m" + vm_params = {"interval": resolution, "record_from": [name_Vm]} vm = nest.Create("voltmeter", params=vm_params) sd = nest.Create("spike_detector") cg = nest.Create("dc_generator", params={"amplitude": 600.}) # connect them and simulate nest.Connect(vm, neuron) - nest.Connect(cg, neuron) + nest.Connect(cg, neuron, syn_spec=add_connect_param.get(model, {})) nest.Connect(neuron, sd) nest.Simulate(simtime) From aa068bc5104344936d0ea19e78069b788d8fb60f Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Tue, 13 Dec 2016 10:52:50 +0100 Subject: [PATCH 04/11] fixed formatting --- models/hh_cond_exp_traub.cpp | 12 ++++++------ pynest/nest/tests/test_refractory.py | 11 +++++------ 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/models/hh_cond_exp_traub.cpp b/models/hh_cond_exp_traub.cpp index c88e16af72..f607a091be 100644 --- a/models/hh_cond_exp_traub.cpp +++ b/models/hh_cond_exp_traub.cpp @@ -136,14 +136,14 @@ hh_cond_exp_traub_dynamics( double, const double y[], double f[], void* pnode ) * ---------------------------------------------------------------- */ nest::hh_cond_exp_traub::Parameters_::Parameters_() - : g_Na( 20000.0 ) // Sodium Conductance (nS) - , g_K( 6000.0 ) // K Conductance (nS) - , g_L( 10.0 ) // Leak Conductance (nS) - , C_m( 200.0 ) // Membrane Capacitance (pF) - , E_Na( 50.0 ) // Reversal potentials (mV) + : g_Na( 20000.0 ) // Sodium Conductance (nS) + , g_K( 6000.0 ) // K Conductance (nS) + , g_L( 10.0 ) // Leak Conductance (nS) + , C_m( 200.0 ) // Membrane Capacitance (pF) + , E_Na( 50.0 ) // Reversal potentials (mV) , E_K( -90.0 ) , E_L( -60.0 ) - , V_T( -63.0 ) // adjusts threshold to around -50 mV + , V_T( -63.0 ) // adjusts threshold to around -50 mV , E_ex( 0.0 ) , E_in( -80.0 ) , tau_synE( 5.0 ) // Synaptic Time Constant Excitatory Synapse (ms) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index c9416f0c02..186a559969 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -72,8 +72,8 @@ nodes = nest.Models("nodes") # list of all neuronal models to be tested -neurons = [m for m in nodes if nest.GetDefaults(m, "element_type") == "neuron" - and m not in ignore_model] +neurons = [m for m in nodes if (nest.GetDefaults( + m, "element_type") == "neuron" and m not in ignore_model)] # additional parameters for the connector add_connect_param = { @@ -100,7 +100,7 @@ def foreach_neuron(func): ''' Decorator that automatically does the test for all neurons. - ''' + ''' def wrapper(*args, **kwargs): self = args[0] msd = 123456 @@ -152,7 +152,6 @@ def compute_reftime(self, model, sd, vm, Vr): t_ref_sim = idx_end * resolution return t_ref_sim - @foreach_neuron def test_refractory_time(self, model): ''' @@ -172,13 +171,13 @@ def test_refractory_time(self, model): nest.Connect(vm, neuron) nest.Connect(cg, neuron, syn_spec=add_connect_param.get(model, {})) nest.Connect(neuron, sd) - + nest.Simulate(simtime) # get and compare t_ref Vr = nest.GetStatus(neuron, "V_reset")[0] t_ref_sim = self.compute_reftime(model, sd, vm, Vr) - + self.assertEqual(t_ref, t_ref_sim, '''Error in model {}: {} != {}'''.format(model, t_ref, t_ref_sim)) From 3d5a56371fb19194e30915e2d55dd16efa2d785f Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Wed, 14 Dec 2016 15:32:27 +0100 Subject: [PATCH 05/11] changed test to check more models --- pynest/nest/tests/test_refractory.py | 138 ++++++++++++++++++--------- 1 file changed, 93 insertions(+), 45 deletions(-) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index 186a559969..1694b1b885 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -26,7 +26,8 @@ import nest """ -Assert that all neuronal models implement refractory period correctly. +Assert that all neuronal models that have a refractory period implement it +correctly (except for Hodgkin-Huxley models which cannot be tested). Details ------- @@ -35,6 +36,28 @@ A ``spike_detector`` is used to detect the time at which the neuron spikes and a ``voltmeter`` is then used to make sure the voltage is clamped to ``V_reset`` during exactly ``t_ref``. + +For neurons that do not clamp the potential, use a very large current to +trigger immediate spiking + +Untested models +--------------- +* ``aeif_cond_alpha_RK5`` +* ``ginzburg_neuron`` +* ``hh_cond_exp_traub`` +* ``hh_psc_alpha`` +* ``hh_psc_alpha_gap`` +* ``ht_neuron`` +* ``iaf_chs_2007`` +* ``iaf_chxk_2008`` +* ``iaf_tum_2000`` +* ``izhikevich`` +* ``mcculloch_pitts_neuron`` +* ``parrot_neuron`` +* ``parrot_neuron_ps`` +* ``pp_pop_psc_delta`` +* ``pp_psc_delta`` +* ``sli_neuron`` """ @@ -43,37 +66,45 @@ # ------------------------- # -ignore_model = [ - "aeif_cond_alpha_RK5", - "aeif_cond_alpha_multisynapse", - "aeif_cond_beta_multisynapse", +# list of all neuronal models that can be tested by looking at clamped V +neurons_V_clamped = [ + 'aeif_cond_alpha', + 'aeif_cond_alpha_multisynapse', + 'aeif_cond_beta_multisynapse', + 'aeif_cond_exp', + 'aeif_psc_alpha', + 'aeif_psc_exp', + 'gif_cond_exp', + 'gif_cond_exp_multisynapse', + 'gif_psc_exp', + 'gif_psc_exp_multisynapse', + 'iaf_cond_alpha', + 'iaf_cond_alpha_mc', + 'iaf_cond_exp', + 'iaf_cond_exp_sfa_rr', + 'iaf_neuron', + 'iaf_psc_alpha', + 'iaf_psc_alpha_multisynapse', + 'iaf_psc_delta', + 'iaf_psc_exp', + 'iaf_psc_exp_multisynapse' +] + +# neurons that must be tested through a high current to spike immediately +# (t_ref = interspike) +neurons_interspike = [ "amat2_psc_exp", - "ginzburg_neuron", - "hh_cond_exp_traub", - "hh_psc_alpha", - "hh_psc_alpha_gap", - "ht_neuron", - "iaf_chs_2007", - "iaf_chxk_2008", + "mat2_psc_exp" +] + +neurons_interspike_ps = [ "iaf_psc_alpha_canon", "iaf_psc_alpha_presc", "iaf_psc_delta_canon", "iaf_psc_exp_ps", - "iaf_tum_2000", - "izhikevich", - "mat2_psc_exp", - "mcculloch_pitts_neuron", - "parrot_neuron", - "parrot_neuron_ps", - "pp_pop_psc_delta", - "pp_psc_delta", - "sli_neuron", ] -nodes = nest.Models("nodes") -# list of all neuronal models to be tested -neurons = [m for m in nodes if (nest.GetDefaults( - m, "element_type") == "neuron" and m not in ignore_model)] +tested_models = neurons_V_clamped + neurons_interspike + neurons_interspike_ps # additional parameters for the connector add_connect_param = { @@ -106,7 +137,7 @@ def wrapper(*args, **kwargs): msd = 123456 N_vp = nest.GetKernelStatus(['total_num_virtual_procs'])[0] pyrngs = [np.random.RandomState(s) for s in range(msd, msd + N_vp)] - for name in neurons: + for name in tested_models: nest.ResetKernel() nest.SetKernelStatus({ 'resolution': resolution, 'grng_seed': msd + N_vp, @@ -120,7 +151,7 @@ class RefractoryTestCase(unittest.TestCase): Check the correct implementation of refractory time in all neuronal models. """ - def compute_reftime(self, model, sd, vm, Vr): + def compute_reftime(self, model, sd, vm, neuron): ''' Compute the refractory time of the neuron. @@ -132,8 +163,8 @@ def compute_reftime(self, model, sd, vm, Vr): GID of the spike detector. vm : tuple GID of the voltmeter. - Vr : double - Value of the reset potential. + neuron : tuple + GID of the recorded neuron. Returns ------- @@ -141,16 +172,24 @@ def compute_reftime(self, model, sd, vm, Vr): Value of the simulated refractory period. ''' spike_times = nest.GetStatus(sd, "events")[0]["times"] - times = nest.GetStatus(vm, "events")[0]["times"] - idx_max = (np.argwhere(np.isclose(times, spike_times[1]))[0] - if len(spike_times) > 1 else -1) - name_Vm = "V_m.s" if model == "iaf_cond_alpha_mc" else "V_m" - Vs = nest.GetStatus(vm, "events")[0][name_Vm] - # get the index at which the spike occured - idx_spike = np.argwhere(times == spike_times[0])[0] - idx_end = np.where(np.isclose(Vs[idx_spike:idx_max], Vr, 1e-6))[0][-1] - t_ref_sim = idx_end * resolution - return t_ref_sim + if model in neurons_interspike: + # spike emitted at next timestep so substract resolution + return spike_times[1]-spike_times[0]-resolution + elif model in neurons_interspike_ps: + return spike_times[1]-spike_times[0] + else: + Vr = nest.GetStatus(neuron, "V_reset")[0] + times = nest.GetStatus(vm, "events")[0]["times"] + idx_max = (np.argwhere(np.isclose(times, spike_times[1]))[0][0] + if len(spike_times) > 1 else -1) + name_Vm = "V_m.s" if model == "iaf_cond_alpha_mc" else "V_m" + Vs = nest.GetStatus(vm, "events")[0][name_Vm] + # get the index at which the spike occured + idx_spike = np.argwhere(times == spike_times[0])[0][0] + idx_end = np.where( + np.isclose(Vs[idx_spike:idx_max], Vr, 1e-6))[0][-1] + t_ref_sim = idx_end * resolution + return t_ref_sim @foreach_neuron def test_refractory_time(self, model): @@ -165,8 +204,12 @@ def test_refractory_time(self, model): name_Vm = "V_m.s" if model == "iaf_cond_alpha_mc" else "V_m" vm_params = {"interval": resolution, "record_from": [name_Vm]} vm = nest.Create("voltmeter", params=vm_params) - sd = nest.Create("spike_detector") + sd = nest.Create("spike_detector", params={'precise_times': True}) cg = nest.Create("dc_generator", params={"amplitude": 600.}) + # for models that do not clamp V_m, use very large current to trigger + # almost immediate spiking => t_ref almost equals interspike + if model not in neurons_V_clamped: + nest.SetStatus(cg, "amplitude", 10000000.) # connect them and simulate nest.Connect(vm, neuron) nest.Connect(cg, neuron, syn_spec=add_connect_param.get(model, {})) @@ -175,11 +218,16 @@ def test_refractory_time(self, model): nest.Simulate(simtime) # get and compare t_ref - Vr = nest.GetStatus(neuron, "V_reset")[0] - t_ref_sim = self.compute_reftime(model, sd, vm, Vr) - - self.assertEqual(t_ref, t_ref_sim, '''Error in model {}: - {} != {}'''.format(model, t_ref, t_ref_sim)) + t_ref_sim = self.compute_reftime(model, sd, vm, neuron) + + # approximate result for precise spikes (interpolation error) + if model in neurons_interspike_ps: + self.assertAlmostEqual(t_ref, t_ref_sim, places=3, + msg='''Error in model {}: + {} != {}'''.format(model, t_ref, t_ref_sim)) + else: + self.assertAlmostEqual(t_ref, t_ref_sim, msg='''Error in model {}: + {} != {}'''.format(model, t_ref, t_ref_sim)) # --------------------------------------------------------------------------- # From b76cb4a09b16245ab69d12411e41c2dff37b668e Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Wed, 14 Dec 2016 15:39:09 +0100 Subject: [PATCH 06/11] corrected tested_neurons list --- pynest/nest/tests/test_refractory.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index 1694b1b885..2c9fda7a4f 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -104,7 +104,32 @@ "iaf_psc_exp_ps", ] -tested_models = neurons_V_clamped + neurons_interspike + neurons_interspike_ps +# models that cannot be tested +ignore_model = [ + "aeif_cond_alpha_RK5", # this one is known to be faulty and will be removed + "ginzburg_neuron", + "hh_cond_exp_traub", + "hh_psc_alpha", + "hh_psc_alpha_gap", + "ht_neuron", + "iaf_chs_2007", + "iaf_chxk_2008", + "iaf_psc_alpha_canon", + "iaf_psc_alpha_presc", + "iaf_psc_delta_canon", + "iaf_psc_exp_ps", + "iaf_tum_2000", + "izhikevich", + "mcculloch_pitts_neuron", + "parrot_neuron", + "parrot_neuron_ps", + "pp_pop_psc_delta", + "pp_psc_delta", + "sli_neuron", +] + +tested_models = [m for m in nest.Models("nodes") if (nest.GetDefaults( + m, "element_type") == "neuron" and m not in ignore_model)] # additional parameters for the connector add_connect_param = { From f3b4b79b55de36b67170250f585078eaec152882 Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Wed, 14 Dec 2016 15:42:07 +0100 Subject: [PATCH 07/11] corrected tested_neurons list 2 --- pynest/nest/tests/test_refractory.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index 2c9fda7a4f..d888fe6cbc 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -106,7 +106,7 @@ # models that cannot be tested ignore_model = [ - "aeif_cond_alpha_RK5", # this one is known to be faulty and will be removed + "aeif_cond_alpha_RK5", # this one is faulty and will be removed "ginzburg_neuron", "hh_cond_exp_traub", "hh_psc_alpha", @@ -114,10 +114,6 @@ "ht_neuron", "iaf_chs_2007", "iaf_chxk_2008", - "iaf_psc_alpha_canon", - "iaf_psc_alpha_presc", - "iaf_psc_delta_canon", - "iaf_psc_exp_ps", "iaf_tum_2000", "izhikevich", "mcculloch_pitts_neuron", From bcdced1c5b454c2c2f0885c478eb30f6cb5ff663 Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Fri, 16 Dec 2016 11:01:02 +0100 Subject: [PATCH 08/11] added ht_neuron --- pynest/nest/tests/test_refractory.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index d888fe6cbc..948bb41d90 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -87,14 +87,15 @@ 'iaf_psc_alpha_multisynapse', 'iaf_psc_delta', 'iaf_psc_exp', - 'iaf_psc_exp_multisynapse' + 'iaf_psc_exp_multisynapse', ] # neurons that must be tested through a high current to spike immediately # (t_ref = interspike) neurons_interspike = [ "amat2_psc_exp", - "mat2_psc_exp" + "mat2_psc_exp", + "ht_neuron", ] neurons_interspike_ps = [ @@ -111,7 +112,6 @@ "hh_cond_exp_traub", "hh_psc_alpha", "hh_psc_alpha_gap", - "ht_neuron", "iaf_chs_2007", "iaf_chxk_2008", "iaf_tum_2000", @@ -229,8 +229,10 @@ def test_refractory_time(self, model): cg = nest.Create("dc_generator", params={"amplitude": 600.}) # for models that do not clamp V_m, use very large current to trigger # almost immediate spiking => t_ref almost equals interspike - if model not in neurons_V_clamped: + if model in neurons_interspike_ps: nest.SetStatus(cg, "amplitude", 10000000.) + elif model in neurons_interspike: + nest.SetStatus(cg, "amplitude", 2000.) # connect them and simulate nest.Connect(vm, neuron) nest.Connect(cg, neuron, syn_spec=add_connect_param.get(model, {})) From 4a84dc180c39e923b91b14cbafbc641b9ed56dbb Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Wed, 11 Jan 2017 09:45:35 +0100 Subject: [PATCH 09/11] corrections in response to @jschuecker's review --- pynest/nest/tests/test_refractory.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index 948bb41d90..12fbabdbdb 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -134,14 +134,14 @@ # --------------------------------------------------------------------------- # -# Simulation and refractory time +# Simulation and refractory time limits # ------------------------- # simtime = 100 resolution = 0.1 -min_steps = 1 -max_steps = 200 +min_steps = 1 # minimal number of refractory steps (t_ref = resolution in ms) +max_steps = 200 # maximal number of steps (t_ref = 200 * resolution in ms) # --------------------------------------------------------------------------- # @@ -201,12 +201,13 @@ def compute_reftime(self, model, sd, vm, neuron): else: Vr = nest.GetStatus(neuron, "V_reset")[0] times = nest.GetStatus(vm, "events")[0]["times"] - idx_max = (np.argwhere(np.isclose(times, spike_times[1]))[0][0] - if len(spike_times) > 1 else -1) + # index of the 2nd spike + idx_max = np.argwhere(times == spike_times[1])[0][0] name_Vm = "V_m.s" if model == "iaf_cond_alpha_mc" else "V_m" Vs = nest.GetStatus(vm, "events")[0][name_Vm] # get the index at which the spike occured idx_spike = np.argwhere(times == spike_times[0])[0][0] + # find end of refractory period between 1st and 2nd spike idx_end = np.where( np.isclose(Vs[idx_spike:idx_max], Vr, 1e-6))[0][-1] t_ref_sim = idx_end * resolution @@ -226,7 +227,7 @@ def test_refractory_time(self, model): vm_params = {"interval": resolution, "record_from": [name_Vm]} vm = nest.Create("voltmeter", params=vm_params) sd = nest.Create("spike_detector", params={'precise_times': True}) - cg = nest.Create("dc_generator", params={"amplitude": 600.}) + cg = nest.Create("dc_generator", params={"amplitude": 900.}) # for models that do not clamp V_m, use very large current to trigger # almost immediate spiking => t_ref almost equals interspike if model in neurons_interspike_ps: From 0b583243c4effbc7c324ae7d28bf483bab170448 Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Wed, 11 Jan 2017 09:52:15 +0100 Subject: [PATCH 10/11] correction 2 --- pynest/nest/tests/test_refractory.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index 12fbabdbdb..662e3f7e07 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -134,7 +134,7 @@ # --------------------------------------------------------------------------- # -# Simulation and refractory time limits +# Simulation time and refractory time limits # ------------------------- # From 964b6141f18845caf5903de65b6374c5ccd257db Mon Sep 17 00:00:00 2001 From: Silmathoron Date: Wed, 11 Jan 2017 10:17:10 +0100 Subject: [PATCH 11/11] corrected pep8 --- pynest/nest/tests/test_refractory.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index 662e3f7e07..24016a05fc 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -140,8 +140,8 @@ simtime = 100 resolution = 0.1 -min_steps = 1 # minimal number of refractory steps (t_ref = resolution in ms) -max_steps = 200 # maximal number of steps (t_ref = 200 * resolution in ms) +min_steps = 1 # minimal number of refractory steps (t_ref = resolution) +max_steps = 200 # maximal number of steps (t_ref = 200 * resolution) # --------------------------------------------------------------------------- #