/*
 *  test_multimeter_accu.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_multimeter_accu - test support for multimeter accumulation mode

Synopsis: (test_multimeter_accu.sli) run -> dies if assertion fails

Description:
1. Test that all recording flags incompatible with /to_accumulator are set to false
2. Simulate iaf_cond_alpha net with poisson input and record detailed and accumulated
   potentials, compare.
3. Run network for different thread numbers to check accumulation across threads.

Author: Plesser
FirstVersion: 2011-02-11
*/

% The following test needs the model iaf_cond_alpha,
% which is only compiled if we have the GSL
statusdict/have_gsl :: {

/unittest (8831) require
/unittest using

/arr (9161) require

M_ERROR setverbosity


/clear_error
{ 
  counttomark npop % pop all but mark
  errordict begin /newerror false def end
} def

% integer_data --- return array with only unique entries in ascending order
%                  data must be intergers
%           
/unique_ascending
{
  /d << >> def
  { cvs cvlit d exch 0 put } forall
  d keys { cvs cvi } Map Sort
}
def

% times data --- times and data have corresponding entries; entries in 
%                data with equalt times are summed.
/sum_by_times
{
  /data Set
  dup { cvs cvlit } Map /times Set
  unique_ascending { cvs cvlit } Map /utimes Set
  /res << >> def
  utimes { res exch 0 put } forall
  [ data times ] { /key Set res key get add res exch key exch put  } ScanThread
  [ ] utimes { res exch get append } forall
}
def


% first test:
% make sure that to_accumulator forces all conflicting options to false
{
  << >> begin
    /mm /multimeter << /to_accumulator true >> Create def
    /mms mm GetStatus def
    true [ /to_file /to_screen /to_memory /withgid /withweight ]
  { mms exch get not and } Fold    
  end
} assert_or_die

% second test:
% connect normal and accumulating multimeter, simulate, compare
{
  ResetKernel
  << >> begin
    /recvars [ /g_in ] def
    /N 100 def

        /nnet /subnet Create def
        nnet ChangeSubnet
	/iaf_cond_alpha N Create ;
        0 ChangeSubnet
        /nrns nnet GetGlobalLeaves def

        /multimeter << /record_from recvars /withgid false /withtime true 
                       /interval 0.7 /start 4.8 /stop 233.1 /time_in_steps true >> SetDefaults
        /mm /multimeter Create def	
        /ac /multimeter << /to_memory false /to_accumulator true >> Create def	
        mm nrns DivergentConnect
        ac nrns DivergentConnect

        /pg /poisson_generator << /rate 1000. >> Create def
        /static_synapse /exc << /weight  1.0 >> CopyModel 
        /static_synapse /inh << /weight -1.0 >> CopyModel  
        pg nrns /exc DivergentConnect
        pg nrns /inh DivergentConnect

        250. Simulate

        /detevs mm /events get def
        /dettimes detevs /times get cva def          
        /accevs ac /events get def	
        /acctimes accevs /times get cva def          
        
        recvars 
        { 
          /var Set
          % get data from individual recordings
	  detevs var get cva
          % compute mean for error estimate
          dup Mean /varmean Set
          % sum results for equal times
          dettimes exch sum_by_times

          % data from accumulated recording
          accevs var get cva

          % test difference relative to mean	  
          sub true exch { varmean div abs 1e-15 lt and } Fold 
	} Map 

        true exch { and } Fold 

        % check that times are equal  
        dettimes unique_ascending acctimes eq and  

  end
} assert_or_die


% third test:
% as second test, but with three threads (NB: thread number should NOT divide number of neurons)
% since we know from the previous test that accumulation works on a single thread, we can
% test against single-thread accumulated results. This saves us a lot of work in matching
% gids to threads.
/run_mma
{
  << >> begin
    /n_threads Set

    /N 100 def  % should not be divisible by thread number
  
    ResetKernel
    0 /num_processes get 1 eq assert % distributed setting not covered
    0 << /local_num_threads n_threads >> SetStatus

    % actual neurons, placed in subnet
    /nnet /subnet Create def
    nnet ChangeSubnet
    /iaf_cond_alpha N << /I_e 40. >> Create ;
    0 ChangeSubnet
    /nrns nnet GetGlobalLeaves def
    nrns { dup /global_id get 1 sub 30. N div mul -90. add 
               /foo << >> def foo exch /V_m exch put foo SetStatus } forall

    % multimeter for accumulated recording
    /ac /multimeter << /record_from [ /V_m ] /withgid false 
                       /to_memory false /to_accumulator true >> Create def	
    ac nrns DivergentConnect

    250. Simulate

    % obtain data 
    ac /events get /V_m get cva
  end
}
def

{
  << >> begin
    1 run_mma dup  % run with single thread
    Mean /mn Set   % store mean
    3 run_mma      % run with three threads

    sub mn div true exch { abs 1e-14 lt and } Fold
  end
} assert_or_die

endusing



} if % have_gsl