fnsqlog102.gms : Test correctness of sqlog10 intrinsic

Description

sqlog10(x,S) = log_10(x)                         if x >= S
               (1/log 10) (ln(S) + r/S - 1/S^2)  otherwise, r = x - S
for S >= 1e-150, default S = 1e-150.

Contributor: Steve, Mar 2017


Small Model of Type : GAMS


Category : GAMS Test library


Main file : fnsqlog102.gms   includes :   fnset_xy.inc [htmlfntest_xy.inc [html]

$title 'Test correctness of sqlog10 intrinsic' (FNSQLOG102,SEQ=724)

$onText
sqlog10(x,S) = log_10(x)                         if x >= S
               (1/log 10) (ln(S) + r/S - 1/S^2)  otherwise, r = x - S
for S >= 1e-150, default S = 1e-150.

Contributor: Steve, Mar 2017
$offText

$include fnset_xy.inc

reps = 8e-16;
sets
  T /
    t1  * t6,
    t11 * t16
    t21 * t26
  /
  torig(t) 'unsmoothed inputs'
  ;

* doing anything with y = 1e-150 is pretty useless: the 2nd deriv is
* so large
table data(T,V)
        x          y       f_      fx_     fxx_    rc_     ec_
t1      1       1e-140
t2   1e-136     1e-140
t3   1e-139     1e-140
t4   1e-160     1e-140
t5     EPS      1e-140
t6   -1e1       1e-140

t11     1       1e-20
t12  1e-19      1e-20
t13  1e-20      1e-20
t14  1e-120     1e-20
t15    EPS      1e-20
t16  -1e40      1e-20

t21   4000       2000
t22   2001       2000
t23   2000       2000
t24  1999.999    2000
t25    EPS       2000
t26  -1e100      2000
;

scalars
  r
  L10
  lnS
  ;
* problematic to test exactly at the breakpoint: should 2nd deriv be 0
* or taken from the smooth part?
loop{T$[data(T,'x') = data(T,'y')],
  data(T,'x') = (1+4e-16) * data(T,'x');
};
torig(t) = [data(T,'x') > data(T,'y')];

L10 = log(10);
loop{torig(T),
  data(T,  'f_') = log(data(T,'x')) / L10;
  data(T, 'fx_') = 1 / (data(T,'x') * L10);
  data(T,'fxx_') = (-1 / data(T,'x')) / (data(T,'x')*L10);
};
loop{t$[not torig(t)],
  lnS = log(data(t,'y'));
  r = data(t,'x') - data(t,'y');
  data(t,  'f_') = (lnS + r/data(t,'y') - 0.5*sqr(r/data(t,'y'))) / L10;
  data(t, 'fx_') = (1 - r/data(t,'y')) / (data(t,'y')*L10);
  data(T,'fxx_') = (-1  / data(t,'y')) / (data(t,'y')*L10);
};

loop {T,
  data(T,'fyx_') = data(T,'fxy_');

  data(T,  'f')  = sqlog10.value(   data(T,'x'),data(T,'y'));
  data(T, 'fx')  = sqlog10.grad(1:  data(T,'x'),data(T,'y'));
  data(T, 'fy')  = sqlog10.grad(2:  data(T,'x'),data(T,'y'));
  data(T,'fxx')  = sqlog10.hess(1:1:data(T,'x'),data(T,'y'));
  data(T,'fxy')  = sqlog10.hess(1:2:data(T,'x'),data(T,'y'));
  data(T,'fyx')  = sqlog10.hess(2:1:data(T,'x'),data(T,'y'));
  data(T,'fyy')  = sqlog10.hess(2:2:data(T,'x'),data(T,'y'));
  data(T, 'rc')  = mathlastrc;
  data(T, 'ec')  = mathlastec;
};

display data;

$include fntest_xy.inc