fnedist.gms : Test correctness of edist intrinsic

Description

Norm 2 function

f := sqrt(sqr(x1)+sqr(x2)+...)

g1.. = x1/f  x2/f  ...

hij = -xi*xj/f^3     for i<>j
      (f^2-xi^2)/f^3 for i=j


Small Model of Type : GAMS


Category : GAMS Test library


Main file : fnedist.gms

$title 'Test correctness of edist intrinsic' (FNEDIST,SEQ=197)
$ontext
  Norm 2 function

  f := sqrt(sqr(x1)+sqr(x2)+...)

  g1.. = x1/f  x2/f  ...

  hij = -xi*xj/f^3     for i<>j
        (f^2-xi^2)/f^3 for i=j
$offtext

Sets 
   arg / n0*n19 /
   T   / p1*p1000 /;
alias (arg,argp);
Parameters
   data(arg)
   testi(*,*)
   testm(*,*), a, b;

loop {T,
  data(arg) = uniform(-100,100);
  testi(''  ,''  )   = edist.value(data('n0') ,data('n1') ,data('n2') ,data('n3') ,data('n4') ,
                                   data('n5') ,data('n6') ,data('n7') ,data('n8') ,data('n9') ,
								   data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
                                   data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
  loop(arg, 
    testi(arg,''  )  = edist.grad (ord(arg): data('n0') ,data('n1') ,data('n2') ,data('n3') ,data('n4') ,
                                             data('n5') ,data('n6') ,data('n7') ,data('n8') ,data('n9') ,
											 data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
											 data('n15'),data('n16'),data('n17'),data('n18'),data('n19') ));
  loop((arg,argp), 
    testi(arg,argp)  = edist.hess (ord(arg):ord(argp): data('n0') ,data('n1') ,data('n2') ,data('n3') ,data('n4') ,
                                                       data('n5') ,data('n6') ,data('n7') ,data('n8') ,data('n9') ,
													   data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
													   data('n15'),data('n16'),data('n17'),data('n18'),data('n19') ));

                                                       
  testm('','') = sqrt(sum(arg, sqr(data(arg))));
  a = 1/testm('','');
  b = a*a*a;
  testm(arg,'')   = data(arg)*a;
  testm(arg,argp) = -data(arg)*data(argp)*b;
* Since we already calculated hess(xi,xi) with the formula for i<>j we can use: hess(xi,xi) = 1/f + hess(xi,xi)
  testm(arg,arg)  = a + testm(arg,arg);

  abort$(abs(                testi( '',  '')-testm( '',  '')) >1e-6) 'wrong function value', data, testi, testm;
  abort$(abs(sum(arg,        testi(arg,  '')-testm(arg,  '')))>1e-6) 'wrong gradient value', data, testi, testm;
  abort$(abs(sum((arg,argp), testi(arg,argp)-testm(arg,argp)))>1e-6) 'wrong hessian  value', data, testi, testm;
};