$title 'Test correctness of LSEMax intrinsic' (FNLSEMAX,SEQ=913)

$onText
  LSE max LSEMax is defined as

  f := log(exp(x1)+exp(x2)+...)
$offText

Sets
   P
   V
   arg 'expected content of set V' / v00*v19 /
   tst(*)
   ;
alias (V,Vp);
alias (u,*);
singleton set PP(P);
scalars
   aErr0
   mx0
   rErr0
   ;
parameters
   xIn(P,V)
   x(V)

   wantFunc(P)
   wntFunc
   gotFunc

   wantGrad(P,V)
   wntGrad(V)
   gotGrad(V)

   wantHess(P,V,V)
   wntHess(V,V)
   gotHess(V,V)

   aErr1(V)
   aErr1Show(V)
   mx1(V)
   rErr1(V)
   rErr1Show(V)

   aErr2(V,V)
   aErr2Show(V,V)
   mx2(V,V)
   rErr2(V,V)
   rErr2Show(V,V)
   ;

$gdxIn fnlsemax.gdx
$load P V xIn wantFunc wantGrad wantHess
$gdxIn

abort$[card(arg) <> card(V)] 'sets arg and V are not equal';
tst(arg) = yes;
tst(V) = no;
abort$[card(tst)] 'sets arg and V are not equal', tst;

* ------------------------ LSEMax -------------------------
loop {P,
  PP(P) = yes;
  x(V) = xIn(P,V);
  wntFunc = wantFunc(P);
  gotFunc = LSEMax.value(x('v00'),x('v01'),x('v02'),x('v03'),x('v04'),
                         x('v05'),x('v06'),x('v07'),x('v08'),x('v09'),
                         x('v10'),x('v11'),x('v12'),x('v13'),x('v14'),
                         x('v15'),x('v16'),x('v17'),x('v18'),x('v19') );

  wntGrad(V) = wantGrad(P,V);
  loop{V,
    gotGrad(V) = LSEMax.grad(ord(V): x('v00'),x('v01'),x('v02'),x('v03'),x('v04'),
                                     x('v05'),x('v06'),x('v07'),x('v08'),x('v09'),
                                     x('v10'),x('v11'),x('v12'),x('v13'),x('v14'),
                                     x('v15'),x('v16'),x('v17'),x('v18'),x('v19') );
  };

  wntHess(V,Vp) = wantHess(P,V,Vp);
  loop{(V,Vp),
    gotHess(V,Vp) = LSEMax.hess(ord(V):ord(Vp): x('v00'),x('v01'),x('v02'),x('v03'),x('v04'),
                                                x('v05'),x('v06'),x('v07'),x('v08'),x('v09'),
                                                x('v10'),x('v11'),x('v12'),x('v13'),x('v14'),
                                                x('v15'),x('v16'),x('v17'),x('v18'),x('v19') );
  };

  aErr0 = abs(gotFunc-wntFunc);
  repeat {
    break$[aErr0 <= 1e-200];
    mx0 = max(abs(gotFunc),abs(wntFunc));
    rErr0 = aErr0 / mx0;
    break$[rErr0 <= 5e-16];
    abort  'wrong function value in LSEMax', PP, x, gotFunc, wntFunc, aErr0, rErr0;
  until yes };

  aErr1(V) = abs(gotGrad(V)-wntGrad(V));
  aErr1Show(V) = aErr1(V)$[aErr1(V) > 1e-200];
  repeat {
    break$[card(aErr1Show) = 0];
    mx1(V) = max(abs(gotGrad(V)),abs(wntGrad(V)));
    rErr1(V) = [aErr1(V) / mx1(V)]$aErr1Show(V);
    rErr1Show(V) = rErr1(V)$[rErr1(V) > 5e-16];
    break$[card(rErr1Show) = 0];
    abort  'wrong gradient value in LSEMax', PP, x, gotGrad, wntGrad, aErr1Show, rErr1Show;
  until yes };

  aErr2(V,Vp) = abs(gotHess(V,Vp)-wntHess(V,Vp));
  aErr2Show(V,Vp) = aErr2(V,Vp)$[aErr2(V,Vp) > 1e-200];
  repeat {
    break$[card(aErr2Show) = 0];
    mx2(V,Vp) = max(abs(gotHess(V,Vp)),abs(wntHess(V,Vp)));
    rErr2(V,Vp) = [aErr2(V,Vp) / mx2(V,Vp)]$aErr2Show(V,Vp);
    rErr2Show(V,Vp) = rErr2(V,Vp)$[rErr2(V,Vp) > 5e-16];
    break$[card(rErr2Show) = 0];
    abort  'wrong Hessian value in LSEMax', PP, x, gotHess, wntHess, aErr2Show, rErr2Show;
  until yes };
};