Description

This model demonstrates how to tests the first and second order derivatives
of a given function by comparing the (exact/analytical) derivatives with
derivatives calculated by finite differences.

GAMS provides through the function suffix .grad/.hess the (exact) derivatives
and the finite difference derivatives through function suffix gradn/.hessn.
These facilities work for intrinisc and extrinsic functions.

Keywords: GAMS language features, accuracy test, analytical derivatives,
          finite difference

Small Model of Type : GAMS

Category : GAMS Model library

Main file : derivtst.gms

$title How to test Derivatives of Functions (DERIVTST,SEQ=406)

$onText
This model demonstrates how to tests the first and second order derivatives
of a given function by comparing the (exact/analytical) derivatives with
derivatives calculated by finite differences.

GAMS provides through the function suffix .grad/.hess the (exact) derivatives
and the finite difference derivatives through function suffix gradn/.hessn.
These facilities work for intrinisc and extrinsic functions.

Keywords: GAMS language features, accuracy test, analytical derivatives,
          finite difference
$offText

$if not set sampleSize $set sampleSize 100
$if not set tol        $set tol        1e-6
$if not set func       $set func       power

Set
   s  'sample'        / s1*s%sampleSize%  /
   ai 'argument info' / lo, up, int, endo /;

* Function data
$ifThenI %func%==pdfnormal
$  funcLibIn mylib stodclib
   function %func% / mylib.%func% /;
$  set args sp(s,'a1'),sp(s,'a2'),sp(s,'a3')
   Set a 'function arguments' / a1*a3 /;

   Table fi(a,ai) 'function argument information'
            lo  up  int  endo
      a1   -10  10          1
      a2    -5   5
      a3   0.5   1;
$elseIfI %func%==power
$  set args sp(s,'a1'),sp(s,'a2')
   Set a 'function arguments' / a1*a2 /;

   Table fi(a,ai) 'function argument information'
            lo   up  int  endo
      a1  -100  100          1
      a2    -5    5    1;
$else
$  abort 'no function data for %func%
$endIf

Parameter
   sp(s,a)      'sample points for argument'
   grad(s,a)    'exact gradient'
   gradn(s,a)   'numerical gradient'
   hess(s,a,a)  'exact gradient'
   hessn(s,a,a) 'numerical gradient';

Set
   ae(a)    'endogenos arguments'
   badsp(s) 'sample points that trigger an execerror';

Alias (ae,aep);

ae(a)   = fi(a,'endo');
sp(s,a) = uniformInt(fi(a,'lo'),fi(a,'up'))$(fi(a,'int'))
        + uniform(fi(a,'lo'),fi(a,'up'))$(not fi(a,'int'));

loop(s,
   grad (s,ae)     = %func%.grad (ae.pos:%args%);
   gradn(s,ae)     = %func%.gradn(ae.pos:%args%);
   hess (s,ae,aep) = %func%.hess (ae.pos:aep.pos:%args%);
   hessn(s,ae,aep) = %func%.hessn(ae.pos:aep.pos:%args%);
   if(execError,
      badsp(s)  = yes;
      execError = 0;
   );
);

* Compare gradient and hessian
Set csp(s) 'sample points to compare';

Parameter
   badargs(s,a)     'sample points that triggered an execution error'
   gradError(s,a)   'error in gradient calculation'
   HessError(s,a,a) 'error in hessian  calculation';

csp(s) = not badsp(s);
badargs(badsp,a) = sp(badsp,a);
$macro Error(a,b) abs(a-b)/abs(a+1)
gradError(csp,ae)     = Error(grad(csp,ae),gradn(csp,ae));
hessError(csp,ae,aep) = Error(hess(csp,ae,aep),hessn(csp,ae,aep));
gradError(csp,ae    )$(gradError(csp,ae)     < %tol%) = 0;
hessError(csp,ae,aep)$(hessError(csp,ae,aep) < %tol%) = 0;

abort$(card(badsp) + card(gradError) + card(hessError)) badargs, gradError, hessError;