/*
* test_iaf_ps_dc_accuracy.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 <http://www.gnu.org/licenses/>.
*
*/
/* BeginDocumentation
Name: testsuite::test_iaf_ps_dc_accuracy - test of accuracy of neuron models subject to DC current
Synopsis: (test_iaf_ps_dc_accuracy) run -> comparison with analytical solution
Description:
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.
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.
Remarks:
The script checks whether the kernel can be appropriately configured.
The script can be used to check whether the accuracy has survived optimization
by the C++ compiler.
In case the accuracy is higher than specified by IEEE arithmetics this might
be due to the use of processor registers. The compiler option
-ffloat-store of the gcc compiler ensures that doubles are not stored in
registers.
FirstVersion: May 2005
Author: March 2009, Diesmann
References:
[1] Rotter S & Diesmann M (1999) Exact simulation of time-invariant linear
systems with applications to neuronal modeling. Biologial Cybernetics
81:381-402.
[2] Morrison A, Straube S, Plesser H E, & Diesmann M (2007) Exact Subthreshold
Integration with Continuous Spike Times in Discrete Time Neural Network
Simulations. Neural Computation 19:47--79
SeeAlso: iaf_psc_alpha_canon, iaf_psc_alpha_presc, iaf_psc_delta_canon, testsuite::test_iaf_ps_dc_t_accuracy
*/
/unittest (8093) require
/unittest using
M_ERROR setverbosity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Parameters of simulation schedule.
%
-14 /min_exponent Set
[0 min_exponent -1] Range /hlist Set
[ % time [ms] tolerated error [mV]
[ 5.0 1e-13 ]
[ 500.0 1e-9 ] % error larger because of accumulation
] /Tlist Set % at very small computation step sizes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Check if kernel accepts high precision
%
0 <<
/tics_per_ms min_exponent neg dexp
/resolution 0 dexp % 1 ms default
>> SetStatus
0 /ms_per_tic get frexp
exch
{0.5 eq} assert_or_die % base 2 tic size?
{1 sub min_exponent leq} assert_or_die % sufficient resolution?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Parameters of neuron model.
%
<<
/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
>> /params Set
params begin userdict begin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Simulation at given resolution returns membrane potential
% at end of simulation time and difference to exact value
% for both implementations.
%
/SimAtResolution
{
dup /i Set
dexp /h Set
ResetKernel % resets tic base and computation time step
0 << /tics_per_ms min_exponent neg dexp /resolution h >> SetStatus
[
/iaf_psc_alpha_canon
/iaf_psc_alpha_presc
/iaf_psc_delta_canon % this list can be extended
]
{Create dup params SetStatus} Map /neurons Set
T Simulate
neurons
{[exch /V_m get dup V sub abs]} Map Flatten
i prepend
} def
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% do all simulation times
%
{
Tlist
{
[/T /tolerance] Set
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Reference value
%
% V is the exact value of the membrane potential at the end
% of the simulation.
%
(I_e * tau_m/C_m * (1. - exp(-T/tau_m)) ) ExecMath /V Set
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Perform simulations at all resolutions and collect results
%
hlist {SimAtResolution} Map
%dup print_details % uncomment for debugging
{Rest 2 Partition [/All 2] Part} Map % select columns containing
% the voltage errors
Flatten {tolerance lt} Map
} Map
Flatten true exch {and} Fold
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The output below was generated on a Pentium 4M with GSL-1.5, gcc-4.0.0 -O2.
%
%
% Theory for the fixed point of the iteration (case T=500ms)
% ----------------------------------------------------------
%
% The difference between the theoretical value reached by the membrane potential
% at large t and the fixed point of the iteration is neither determined by
% the accuracy exponentials can be computed nor by the accumulation of round-off
% errors. The only relevant term is the smallest increment to the membrane potential
% that can be represented in the vicinity the fixed point. The increment in
% the iteration is the product of the distance from the limiting value and a term
% decreasing with decreasing h. Thus, for smaller h a larger error is
% to be tolerated.
%
%
% The iteration for the membrane potential is of the form
%
% y_{i+1} = p * (y_i - y^*) + y_i
%
% where y^* is the theoretical fixed point
%
% y^* = Tau/C * I0
%
% and p element in (-1,0] an h-dependent coefficient (i.e. p=expm1(-h/Tau) ).
% The iteration
%
% for(int j=0; j<(500/ldexp(1,i)); j++)
% {
% y = p * (y - ystar) + y;
% }
%
% exactly reproduces the values of the fourth column in the second table below.
%
% The fixed point y of the iteration is reached when p*(y-y^*) is so small
% that compared to y in double representation the sum of the two is
% indistinguishable from y.
% Let us assume that y is close to the theoretical value. The smallest
% representable difference is then given by
%
% s = 2^floor( log2(y^*) ) * eps
%
% where eps is the machine epsilon. In C++ this reads
%
% x=std::frexp(ystar,&p); s=std::ldexp(1,p -1)*std::numeric_limits<double>::epsilon();
%
% The difference that can just not be represented is s/2.
% The condition for the fixed point of the iteration therefore is
%
% s/2 = p * (y - y^*)
%
% solving for the difference dy between the fixed point and its
% theoretical value we obtain
%
% dy = s/2/p
% or
% y = y^* + s/2/p
%
% The fixed points predicted by these considerations are:
%
% h in ms y^* + s/2/p
% -------------------------------------------
% 0 40.0000000000000
% -1 39.9999999999999
% -2 39.9999999999999
% -3 39.9999999999997
% -4 39.9999999999994
% -5 39.9999999999989
% -6 39.9999999999977
% -7 39.9999999999955
% -8 39.9999999999909
% -9 39.9999999999818
% -10 39.9999999999636
% -11 39.9999999999272
% -12 39.9999999998545
% -13 39.9999999997090
% -14 39.9999999994179
%
%
%
% Exact value of membrane potential after 5 ms is 15.7387736114947 mV.
%
% h in ms exp*y [mV] error [mV] expm1*y + y [mV] error [mV]
% ---------------------------------------------------------------------------------------------------------------
% 0 15.7387736114947 7.105427357601e-15 15.7387736114947 0
% -1 15.7387736114947 0 15.7387736114947 1.77635683940025e-15
% -2 15.7387736114947 2.8421709430404e-14 15.7387736114947 1.77635683940025e-15
% -3 15.7387736114947 1.24344978758018e-14 15.7387736114947 0
% -4 15.7387736114947 5.50670620214078e-14 15.7387736114947 5.32907051820075e-15
% -5 15.7387736114947 7.105427357601e-15 15.7387736114947 1.77635683940025e-15
% -6 15.7387736114948 1.35003119794419e-13 15.7387736114947 3.5527136788005e-15
% -7 15.7387736114945 1.15463194561016e-13 15.7387736114946 1.4210854715202e-14
% -8 15.7387736114961 1.4228618283596e-12 15.7387736114947 3.5527136788005e-15
% -9 15.7387736114945 1.49213974509621e-13 15.7387736114947 0
% -10 15.738773611496 1.30384592011978e-12 15.7387736114946 1.4210854715202e-14
% -11 15.7387736114969 2.22577511976851e-12 15.7387736114947 1.77635683940025e-15
% -12 15.7387736114679 2.67803557107982e-11 15.7387736114947 5.32907051820075e-15
% -13 15.7387736115287 3.40296679723906e-11 15.7387736114947 7.46069872548105e-14
% -14 15.7387736114987 4.08384437378118e-12 15.7387736114946 3.01980662698043e-14
%
%
% Exact value of membrane potential after 500 ms is 40 mV.
%
% h in ms exp*y [mV] error [mV] expm1*y + y [mV] error [mV]
% ---------------------------------------------------------------------------------------------------------------
% 0 40 3.5527136788005e-14 40 3.5527136788005e-14
% -1 39.9999999999999 7.105427357601e-14 39.9999999999999 7.105427357601e-14
% -2 39.9999999999999 1.4210854715202e-13 39.9999999999999 1.35003119794419e-13
% -3 39.9999999999997 2.8421709430404e-13 39.9999999999997 2.8421709430404e-13
% -4 39.9999999999994 5.61328761250479e-13 39.9999999999994 5.61328761250479e-13
% -5 39.9999999999989 1.13686837721616e-12 39.9999999999989 1.13686837721616e-12
% -6 39.9999999999977 2.26663132707472e-12 39.9999999999977 2.27373675443232e-12
% -7 39.9999999999955 4.54036808150704e-12 39.9999999999955 4.54747350886464e-12
% -8 39.9999999999909 9.09494701772928e-12 39.9999999999909 9.09494701772928e-12
% -9 39.9999999999818 1.81898940354586e-11 39.9999999999818 1.81898940354586e-11
% -10 39.9999999999636 3.63939989256323e-11 39.9999999999636 3.63939989256323e-11
% -11 39.9999999999272 7.27879978512647e-11 39.9999999999272 7.27879978512647e-11
% -12 39.9999999998544 1.45590206557245e-10 39.9999999998544 1.45590206557245e-10
% -13 39.9999999997088 2.91180413114489e-10 39.9999999997088 2.91180413114489e-10
% -14 39.9999999994176 5.82360826228978e-10 39.9999999994176 5.82360826228978e-10
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% helper function for debugging,
% prints detailed table of results
%
/print_details
{
cout default 15 setprecision
endl
endl
(Exact value of membrane potential after ) <-
T <- ( ms is ) <-
V <- ( mV.) <- endl
endl
( h in ms ) <-
( alpha_canon [mV]) <-
( error [mV]) <-
( alpha_presc [mV]) <-
( error [mV]) <-
( delta_canon [mV]) <-
( error [mV]) <-
endl
(--------------------------) <-
(--------------------------) <-
(--------------------------) <-
(--------------------------) <-
(--------------------------) <-
(--------------------------) <-
(------------------------) <-
endl
exch
{
{
exch 24 setw exch <- ( ) <-
}
forall
endl
}
forall
;
}
def
% executes the overall test
assert_or_die