/*
 *  test_lambertw.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_lambertw - check if lambertw function works

Synopsis: (test_lambertw) run -> some known relations

Description:
The script tests whether the Lambert W-function [1] provided by the 
GNU Scientific Library (GSL) [2] fulfills some known relations [3].
In the absence of the GSL NEST uses a simple iterative scheme to
compute the Lambert-W function. In this case we apply less strict
criteria for the required accuracy.
The relationships tested are:
 (1) the principal branch crosses the origin
 (2) at -1/e both branches meet athe value W=-1
 (3) the principal branch fulfills the "golden ratio" of exponentials
 (4) the non-principal branch yields the same result as we find by
     bisectioning for the problem given as an example in the documentation
     of LambertWm1 with realistic parameters.

References:
  [1] Corless, R. M., Gonnet, G. H., Hare, D. E. G., Jeffrey, D. J., & Knuth, D. E.
      (1996). On the lambert w function. Advances in Computational Mathematics 
      5, 329{359.
  [2] Galassi, M., Davies, J., Theiler, J., Gough, B., Jungman, G., Booth, M.,
      & Rossi, F. (2006). GNU Scientific Library Reference Manual (2nd Ed.).
      Network Theory Limited.
  [3] Weisstein, E. W. (1999). CRC Concise Encyclopedia of Mathematics. 
      "Lambert W-Function" Boca Raton, London, New York, Washington D.C.: CRC Press.

Author:  July 2009, Diesmann
SeeAlso: LambertW, LambertW0, LambertWm1, testsuite::test_iaf_psp_peak, testsuite::test_iaf_psp_normalized
*/

/unittest (6335) require
/unittest using

M_ERROR setverbosity




% (1)
%  LambertW0 crosses the origin

0. LambertW0 0. eq                 assert_or_die


% (2)
% at -1/e the values of both branches coincide with W=-1

statusdict/have_gsl :: 
{
 -1. exp neg LambertW0  -1. eq      assert_or_die
 -1. exp neg LambertWm1 -1. eq      assert_or_die
}
{
 -1. exp neg LambertW0  -1. sub abs 1e-8 lt  assert_or_die
 -1. exp neg LambertWm1 -1. sub abs 1e-8 lt  assert_or_die
}
ifelse


% (3)
%  W(1)=0.567143.. is called the omega constant 
%  or the  "golden ratio" of exponentials  because of the relation 
%     exp -W(1) == W(1), 
%  see [2]

statusdict/have_gsl :: 
{
 1. LambertW0   1. LambertW0 neg exp eq                 assert_or_die
}
{
 1. LambertW0   1. LambertW0 neg exp sub abs 1e-12 lt   assert_or_die
}
ifelse 



% (4)
% In the example of finding the peak time of the postsynaptic potential (PSP)
% generated by an alpha-function shaped current (PSC) the Lambert-W function
% is evaluated for the argument -exp(-1/a)/a, where a is the ratio of the 
% membrane time constant tau_m and the time constant of the alpha-function
% tau_syn. The physically meaning full solution leading to a positive peak time
% is located on the non-principal branch of the Lambert-W function LambertWm1.
% Here, we evaluate LambertWm1 for the argument discussed above with realistic 
% choices of the parameters. The result is compared with the value obtained by
% a simple bisectioning algorithm.


20.  /tau_m  Set
0.5  /tau_syn Set
tau_m tau_syn div /a Set


(LambertWm1(-exp(-1./a)/a)) ExecMath 

[/w] (w*exp(w)+exp(-1./a)/a) Function  -10. -1.  1e-11 FindRoot

sub abs 1e-9 lt assert_or_die