$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 gotGrad(arg) wantGrad(arg) 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, gotGrad(arg) = LSEMax.grad(ord(arg): 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') ); wantGrad(arg) = LSEMax.gradn(ord(arg): 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') ); }; 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) = abs(gotGrad(arg)-wantGrad(arg)); aErr1(arg)$[aErr1(arg) <= aTol1] = 0; abort$[card(aErr1) > 0] 'wrong gradient value in LSEMax', PP, data, gotGrad, wantGrad, aErr1; 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 }; }; gotGrad(arg) = 0; wantGrad(arg) = 0; 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], gotGrad(arg) = LSEMaxSc.grad(ord(arg): 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') ); wantGrad(arg) = LSEMaxSc.gradn(ord(arg): 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') ); }; 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) = abs(gotGrad(arg)-wantGrad(arg)); aErr1(arg)$[aErr1(arg) <= aTol1] = 0; abort$[card(aErr1) > 0] 'wrong gradient value in LSEMaxSc', PP, data, gotGrad, wantGrad, aErr1; 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, gotGrad(arg) = LSEMin.grad(ord(arg): 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') ); wantGrad(arg) = LSEMin.gradn(ord(arg): 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') ); }; 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) = abs(gotGrad(arg)-wantGrad(arg)); aErr1(arg)$[aErr1(arg) <= aTol1] = 0; abort$[card(aErr1) > 0] 'wrong gradient value in LSEMin', PP, data, gotGrad, wantGrad, aErr1; 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 }; }; gotGrad(arg) = 0; wantGrad(arg) = 0; 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], gotGrad(arg) = LSEMinSc.grad(ord(arg): 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') ); wantGrad(arg) = LSEMinSc.gradn(ord(arg): 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') ); }; 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) = abs(gotGrad(arg)-wantGrad(arg)); aErr1(arg)$[aErr1(arg) <= aTol1] = 0; abort$[card(aErr1) > 0] 'wrong gradient value in LSEMinSc', PP, data, gotGrad, wantGrad, aErr1; 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 }; };