fnlse.gms : Rough correctness test for LSE max/min intrinsics

Description

```LSE max is defined as:    LSEMax(x) := log(exp(x1)+exp(x2)+...)
LSE min is defined as:    LSEMin(x) := -LSEMax(-x) = -log(exp(-x1)+exp(-x2)+...)

scaled LSE max is defined as:    LSEMaxSc(t,x) := LSEMax(tx) / t
scaled LSE min is defined as:    LSEMinSc(t,x) := LSEMin(tx) / t

This is not a very precise test: we simply verify that the direct/naive
computation of the function is pretty close to what the intrinsics produce,
and that the numerical gradients and Hessians are pretty close to
what the intrinsics produce.
```

Small Model of Type : GAMS

Category : GAMS Test library

Main file : fnlse.gms

``````\$title 'Rough correctness test for LSE max/min intrinsics' (FNLSE,SEQ=912)

\$ontext
LSE max is defined as:    LSEMax(x) := log(exp(x1)+exp(x2)+...)
LSE min is defined as:    LSEMin(x) := -LSEMax(-x) = -log(exp(-x1)+exp(-x2)+...)

scaled LSE max is defined as:    LSEMaxSc(t,x) := LSEMax(tx) / t
scaled LSE min is defined as:    LSEMinSc(t,x) := LSEMin(tx) / t

This is not a very precise test: we simply verify that the direct/naive
computation of the function is pretty close to what the intrinsics produce,
and that the numerical gradients and Hessians are pretty close to
what the intrinsics produce.
\$offtext

Sets
arg / n00*n19 /
P   / p1*p1000 /;
alias (arg,argp);
singleton set PP(P);
scalars
t
aErr0
aTol0  / 3e-14 /
aTol1  / 4e-7 /
aTol2  / 1e-4 /
rTol2 'numerical Hessians are not very precise' / .9 /
;
Parameters
data(arg)
gotFunc
wantFunc
aErr1(arg)
gotHess(arg,arg)
wantHess(arg,arg)
aErr2(arg,arg)
mx2(arg,argp)
rErr2(arg,argp)
;

execseed = 20021215;
option fdDelta = 1.0e-5;
* ------------------------ LSEMax -------------------------
loop {P,
PP(P) = yes;
data(arg) = uniform(-10,10);
gotFunc = LSEMax.value(data('n00'),data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
wantFunc = log(sum(arg, exp(data(arg))));

loop{arg,
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
};
loop{(arg,argp),
gotHess(arg,argp) = LSEMax.hess(ord(arg):ord(argp): data('n00'),data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
wantHess(arg,argp) = LSEMax.hessn(ord(arg):ord(argp): data('n00'),data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
};

aErr0 = abs(gotFunc-wantFunc);
abort\$[aErr0 > aTol0] 'wrong function value in LSEMax', PP, data, gotFunc, wantFunc, aErr0;

aErr1(arg)\$[aErr1(arg) <= aTol1] = 0;

aErr2(arg,argp) = abs(gotHess(arg,argp)-wantHess(arg,argp));
aErr2(arg,argp)\$[aErr2(arg,argp) <= aTol2] = 0;
repeat {
break\$[card(aErr2) = 0];
mx2(arg,argp) = max(abs(gotHess(arg,argp)),abs(wantHess(arg,argp)));
rErr2(arg,argp) = [aErr2(arg,argp) / mx2(arg,argp)]\$aErr2(arg,argp);
rErr2(arg,argp)\$[rErr2(arg,argp) <= rTol2] = 0;
break\$[card(rErr2) = 0];
abort  'wrong Hessian value in LSEMax', PP, data, gothess, wantHess, aErr2, rErr2;
until yes };
};

gotHess(arg,argp) = 0;
wantHess(arg,argp) = 0;
* ------------------------ LSEMaxSc -------------------------
loop {P,
PP(P) = yes;
t = uniform(-1,2);  t = 10**t;
data(arg) = uniform(-15,7);  data('n00') = 0;
gotFunc = LSEMaxSc.value(        t,data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
wantFunc = log(sum(arg\$[ord(arg) > 1], exp(t*data(arg)))) / t;

loop{arg\$[ord(arg) > 1],
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
};
loop{(arg,argp)\$[(ord(arg) > 1) and (ord(argp) > 1)],
gotHess(arg,argp) = LSEMaxSc.hess(ord(arg):ord(argp):         t,data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
wantHess(arg,argp) = LSEMaxSc.hessn(ord(arg):ord(argp):         t,data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
};

aErr0 = abs(gotFunc-wantFunc);
abort\$[aErr0 > aTol0] 'wrong function value in LSEMaxSc', PP, data, gotFunc, wantFunc, aErr0;

aErr1(arg)\$[aErr1(arg) <= aTol1] = 0;

aErr2(arg,argp) = abs(gotHess(arg,argp)-wantHess(arg,argp));
aErr2(arg,argp)\$[aErr2(arg,argp) <= aTol2] = 0;
repeat {
break\$[card(aErr2) = 0];
mx2(arg,argp) = max(abs(gotHess(arg,argp)),abs(wantHess(arg,argp)));
rErr2(arg,argp) = [aErr2(arg,argp) / mx2(arg,argp)]\$aErr2(arg,argp);
rErr2(arg,argp)\$[rErr2(arg,argp) <= rTol2] = 0;
break\$[card(rErr2) = 0];
abort  'wrong Hessian value in LSEMaxSc', PP, data, gothess, wantHess, aErr2, rErr2;
until yes };
};

* ------------------------ LSEMin -------------------------
loop {P,
PP(P) = yes;
data(arg) = uniform(-10,10);
gotFunc = LSEMin.value(data('n00'),data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
wantFunc = -log(sum(arg, exp(-data(arg))));

loop{arg,
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
};
loop{(arg,argp),
gotHess(arg,argp) = LSEMin.hess(ord(arg):ord(argp): data('n00'),data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
wantHess(arg,argp) = LSEMin.hessn(ord(arg):ord(argp): data('n00'),data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
};

aErr0 = abs(gotFunc-wantFunc);
abort\$[aErr0 > aTol0] 'wrong function value in LSEMin', PP, data, gotFunc, wantFunc, aErr0;

aErr1(arg)\$[aErr1(arg) <= aTol1] = 0;

aErr2(arg,argp) = abs(gotHess(arg,argp)-wantHess(arg,argp));
aErr2(arg,argp)\$[aErr2(arg,argp) <= aTol2] = 0;
repeat {
break\$[card(aErr2) = 0];
mx2(arg,argp) = max(abs(gotHess(arg,argp)),abs(wantHess(arg,argp)));
rErr2(arg,argp) = [aErr2(arg,argp) / mx2(arg,argp)]\$aErr2(arg,argp);
rErr2(arg,argp)\$[rErr2(arg,argp) <= rTol2] = 0;
break\$[card(rErr2) = 0];
abort  'wrong Hessian value in LSEMin', PP, data, gothess, wantHess, aErr2, rErr2;
until yes };
};

gotHess(arg,argp) = 0;
wantHess(arg,argp) = 0;
* ------------------------ LSEMinSc -------------------------
loop {P,
PP(P) = yes;
t = uniform(-1,2);  t = 10**t;
data(arg) = uniform(-7,15);  data('n00') = 0;
gotFunc = LSEMinSc.value(        t,data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
wantFunc = -log(sum(arg\$[ord(arg) > 1], exp(-t*data(arg)))) / t;

loop{arg\$[ord(arg) > 1],
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
};
loop{(arg,argp)\$[(ord(arg) > 1) and (ord(argp) > 1)],
gotHess(arg,argp) = LSEMinSc.hess(ord(arg):ord(argp):         t,data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
wantHess(arg,argp) = LSEMinSc.hessn(ord(arg):ord(argp):         t,data('n01'),data('n02'),data('n03'),data('n04'),
data('n05'),data('n06'),data('n07'),data('n08'),data('n09'),
data('n10'),data('n11'),data('n12'),data('n13'),data('n14'),
data('n15'),data('n16'),data('n17'),data('n18'),data('n19') );
};

aErr0 = abs(gotFunc-wantFunc);
abort\$[aErr0 > aTol0] 'wrong function value in LSEMinSc', PP, data, gotFunc, wantFunc, aErr0;

aErr1(arg)\$[aErr1(arg) <= aTol1] = 0;

aErr2(arg,argp) = abs(gotHess(arg,argp)-wantHess(arg,argp));
aErr2(arg,argp)\$[aErr2(arg,argp) <= aTol2] = 0;
repeat {
break\$[card(aErr2) = 0];
mx2(arg,argp) = max(abs(gotHess(arg,argp)),abs(wantHess(arg,argp)));
rErr2(arg,argp) = [aErr2(arg,argp) / mx2(arg,argp)]\$aErr2(arg,argp);
rErr2(arg,argp)\$[rErr2(arg,argp) <= rTol2] = 0;
break\$[card(rErr2) = 0];
abort  'wrong Hessian value in LSEMinSc', PP, data, gothess, wantHess, aErr2, rErr2;
until yes };
};
``````
GAMS Development Corp.
GAMS Software GmbH

General Information and Sales
U.S. (+1) 202 342-0180
Europe: (+49) 221 949-9170