/*
 *  test_spike_poisson_ps_base2_f.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/>.
 *
 */

%  THIS FILE IS A DESIGN STUDY IT WILL PROBABLY NEVER RUN

/*BeginDocumentation
Name: testsuite::test_spike_poisson_ps_base2_f.sli

Synopsis: (test_spike_poisson_ps_base2_f) run

Description:
In NEST spike times are represented by an integer part in units of h
and an offset in units of milliseconds. This helps to maintain a
uniform absolute accuracy of spike times even for long simulation
times.

The user has access to the two components of the spike if the spike 
detector is set to
   /time_in_steps true  /precise_times true .
In thie case the spike detector returns events with the properties
/times and /offsets, where /times are the integer parts s in units of h
and offsets are the fractional parts o in milliseconds. According to
NEST's definition of a grid-constrained spike as a spike occuring 
somewhere in (t-h,t], the precise spike time is
   t = s*h - o

Access to spike times with a uniform absolute accuracy is benefitial
when testing and comparing integrators for neuron models, see appendix
A.2 of [1] for details.

This script tests whether the accuracy of spike times is maintained
independent of the choice of computation step size h.

This assumes that also the poisson generator is capable of maintaining
the accuracy independent of computation step size.
If this test fails go back to
 test_spike_poisson_ps.sli
to check whether poisson_generator_ps can emit spike times at double
precision or whether spike times are limited to the precision of a
tic.

References:
[1] 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

Author: May 2010, adapted to NEST2, Diesmann
SeeAlso: testsuite::test_spike_poisson_ps
*/


%`  % <---- functional expression for immediate compilation and execution
   %       is enclosed by this symbol or similar
   % maybe ` is the better character.

      % --------- the prime indicates a literal ('x translated to /x)
      % |
      % v

(
require('unittest, "(6688)")
using('unittest)

setverbosity(M_ERROR)

min_exponent            = -20;
spike_absolute_accuracy = inv(1e15);
T                       = 4.0;  
h_min                   = 2^min_exponent;


Transmission = Function(
 {
  ResetKernel()

  SetStatus(0, << 'tics_per_ms : 2^-min_exponent, 'resolution :  h >>)

  sp = Create('spike_detector,  
   << 'withgid : false, 'to_memory : true, 'withtime : true, 'time_in_steps : true, 'precise_times : true >> 
  );

  pn = Create('poisson_generator_ps, << 'rate : 16384.0 >>);

  Connect(pn,sp)

  Simulate(T)
                          
  sp // {get(#,['events, ['times, 'offsets]])} // cva

 },

 'h
);



Map( Range([0, min_exponent, -2]),
 Function( 
  {
   Transmission(2^e) // 
     { [ 2^(e-min_exponent)*First(#), Last(#)] }    
   },
  'e 
 )
)


// 
 Function(
  {                                  
   l=Last(r);                        
   Most(r) // {Map(#, {join(#,l)})}
  },
  'r
 )


//
{
 Map(#,
  {
   Apply(#, Function(            
    {(s-sr)*h_min - o + or},     
    's,'o,'sr,'or                
   )) 
  }
 )
}

// Flatten
// {Select(#,{ abs(#) > spike_absolute_accuracy})}
// { # == []}
// assert_or_die

)

CompileMath exec

% all time steps in units of smallest h

% combine results with the one at smallest h as reference
% to 4-tuples [s o sr or]


% s o sr or, we compute (s*hmin-o)-(sr*hmin-or)
% trying to keep the highest possible accuracy: 
% (s-sr)*hmin - o + or

% comment these two last lines to inspect the 
% differences causing the failure