$title 'Test correctness of sqexp intrinsic' (FNSQEXP2,SEQ=720)

$onText
sqexp(x,S) =   exp(x)                     if x <= S
               exp(S) * (1 + r + r^2/2)   otherwise, r = x - S
for S <= 150, default S = 150.

Contributor: Steve, Mar 2017
$offText

$include fnset_xy.inc

scalar xd;
reps = 8e-16;
sets
  T /
    t1  * t6,
    t11 * t16
    t21 * t26
  /
  tmap(t,t) /
    t3. (t4*t6)
    t13.(t14*t16)
    t23.(t24*t26)
  /
  tc(t)
  ;
alias (t,t2);
tc(t2) = sum{t, tmap(t,t2)};

table data(T,V)
        x       y       f_      fx_     fxx_    rc_     ec_
t1      1      150
t2     149     150
t3     150     150
t4   150.001   150
t5     151     150
t6    1e117    150

t11     1       20
t12    19       20
t13    20       20
t14  20.0001    20
t15    30       20
t16   1e145     20

t21   -50      -20
t22   -21      -20
t23   -20      -20
t24  -19.9999  -20
t25     0      -20
t26   1e148    -20
;

* loop{T$[data(T,'x') <= data(T,'y')],
loop{T$[not tc(t)],
  data(T,  'f_') = exp(data(T,'x'));
  data(T, 'fx_') = data(T,'f_');
  data(T,'fxx_') = data(T,'f_');
};
loop{tmap(t,t2),
  xd = data(t2,'x')-data(t,'x');
  data(t2,  'f_') = data(t,'f_') * (1 + xd + 0.5*sqr(xd));
  data(t2, 'fx_') = data(t,'f_') * (1 + xd);
  data(t2,'fxx_') = data(t,'f_');
};

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

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

display data;

$include fntest_xy.inc