{$L ./gcvspl.o}
{$H+}          // ansi strings
{$C+}          // assertions on

unit gcv09;

interface

uses mtx09;

var // kludge to simplify prototype see gcvspl_
    // used in avatar 3 and 4 ...
  xweights: double;
  mode: longint;
  valr: double;
  verbose: boolean = false; // if verbose dump c matrix in gcv3.txt

// (1) and (2) use static c mtx internally whereas (3) returns this
// c matrix for further use in splder() to interp. later

// cross val. unweighted data single data set mode 2 not sure to need this ? may be
// take mode=2 if ese=0 in *2*
// in this mode I think that derivatives are useless
function gcvspl(x,y,u:mtx; var s0, s1:mtx):boolean; // *1* m5

// idem with ese of yi assuming same for all mode 3
// this seems to be my prefered mode using ese=25 the actual ese !
// to get a proper smoothing i.e. a decent first derivative s1
function gcvspl(x,y,u:mtx; ese:double; var s0, s1:mtx):boolean;  // *2* mmmdmm todo with s2 mmmdmmm*

// smooth x y result in c to use with splder...
function gcvspl(x,y:mtx; ese:double; var c:mtx):boolean;  // *3* mmdm

// Ignores xweights but uses mode and valr
// caller should take care that w makes sense with mode and valr...
// intended use: mode=3 valr=1.0 w_j = fact/ese^2 of j^eme data
// NOTE mode=1 valr=0 should do no smoothing same as fmm spline.f
function gcvspl(x,y,w:mtx; var cc:mtx):boolean;  // *4* mmmM

// after calling gcvspl() *3* one can calculate interpolated values
// v for the vector u along with the firts and second derivative.
function splder(x: mtx; cc:mtx; u:mtx; var v, d1, d2:mtx):boolean;
// function splder(x: mtx; c:mtx; u:double; var v, d1, d2:double):boolean;
// TODO scalar version

implementation

{$I gcvspl.ppi} // Pascal prototypes of gcvspl fortran function

// ===========================================================================

var
  // NOTE mtx should be static so that e field is init with NIL
  wx: mtx; // TODO get rid of these global needed by interim avatars
  wy: mtx; // final should be *4*
  c: mtx;
  wk: mtx;           // worspace of gcvspl_
  q: mtx;            // workspace of splder_
  ider: longint = 0; // static var for splder retains previous index

// ===========================================================================

function gcvspl(x,y,u:mtx; var s0, s1:mtx):boolean; // *1* mmmMM
var
  n: longint = 0;
  nc: longint = 0;
  ny: longint = 0;
  m: longint = 2;
  k: longint = 1;
  md: longint = 2;
  val: double = 0.0;
  ier: longint = -999;  // should return 0 if all ok or >0 if err. in fortran
  ider0: longint = 0;   // splder() calc. function
  ider1: longint = 1;   // splder() calc. first derivative
  l: longint;
  i: longint;
  xi, yi0, yi1: double;
begin
  zer(s0,lgn(x),col(x));
  zer(s1,lgn(x),col(x));
  one(wy,1,1);
  n := lgn(x); // TODO lgn(x) or col(x)
  ny := n;
  nc := n; // I use 3 var with same value because fortran uses var args
  // TODO sanity check of args dims
  one(c,n,1);
  zer(wk,n+6*(n*m+1),1);
  one(wx,n,1);
  gcvspl_(ptr(x),ptr(y),ny,ptr(wx),ptr(wy),m,n,k,md,val,ptr(c),nc,ptr(wk),ier);
  //if verbose then mwr(c,'gcv1.txt');
  zer(q,2*m,1); // wrk space for splder_()
  for i := 1 to n do begin
    xi := elm(x,i,1);
    yi0  := splder_(ider0,m,n,xi,ptr(x),ptr(c),l,ptr(q));
    elm(s0,i,1,yi0);
    yi1  := splder_(ider1,m,n,xi,ptr(x),ptr(c),l,ptr(q));
    elm(s1,i,1,yi1);
  end;
  gcvspl := (ier=0)
end;

// ===========================================================================

// only 3 lines changed from previous: function: md=3, val is ese
// and init wx = 1/ese^2 the last seems redundant with val ?
function gcvspl(x,y,u:mtx; ese:double; var s0, s1:mtx):boolean; // *2* mmmdMM
var
  n: longint = 0;
  nc: longint = 0;
  ny: longint = 0;
  m: longint = 2;
  k: longint = 1;
  md: longint = 3; // HERE yi ese provided -> build wx = 1/ese^2
  //val: double = 0.0;
  ier: longint = -999;  // should return 0 if all ok or >0 if err. in fortran
  ider0: longint = 0;   // splder() calc. function
  ider1: longint = 1;   // splder() calc. first derivative
  l: longint;
  i: longint;
  xi, yi0, yi1: double;
begin
  gcvspl := false;
  zer(s0,lgn(x),col(x));
  zer(s1,lgn(x),col(x));
  one(wy,1,1);
  n := lgn(x); // TODO lgn(x) or col(x)
  ny := n;
  nc := n; // I use 3 var with same value because fortran uses var args
  // TODO sanity check of args dims
  one(c,n,1); // GLOBAL static mtx c
  zer(wk,n+6*(n*m+1),1);
  elm(wx,n,1,1/ese/ese); // was dim() in mtx3
  gcvspl_(ptr(x),ptr(y),ny,ptr(wx),ptr(wy),m,n,k,md,ese,ptr(c),nc,ptr(wk),ier); // HERE
  if ier>0 then begin
    writeln('ier=',ier);
    exit
  end;
  zer(q,2*m,1); // wrk space for splder_()
  for i := 1 to n do begin
    xi := elm(x,i,1);
    yi0  := splder_(ider0,m,n,xi,ptr(x),ptr(c),l,ptr(q));
    elm(s0,i,1,yi0);
    yi1  := splder_(ider1,m,n,xi,ptr(x),ptr(c),l,ptr(q));
    elm(s1,i,1,yi1);
  end;
  //if verbose then mwr(c,'gcv2.txt');
  gcvspl := (ier=0)
end;

// ===========================================================================

// exactly the same as the previous one but c is passed as var arg.
// and there are no evaluation in s0 and s1... result is c mtx
// so that one can call the eval function ( splder() ) using several
// different x and c mtx in the calling program
// TODO sanity check of args dims
function gcvspl(x,y:mtx; ese:double; var c:mtx):boolean;  // *3* mmdm
var
  n: longint = 0;
  nc: longint = 0;
  ny: longint = 0;
  m: longint = 2;
  k: longint = 1;
  md: longint = 3; // HERE yi ese provided -> build wx = 1/ese^2
  ier: longint = -999;  // should return 0 if all ok or >0 if err. in fortran
  evar: double; // ese^2
begin
  evar := ese*ese;
  one(wy,1,1);
  n := lgn(x); // TODO lgn(x) or col(x)
  ny := n;
  nc := n; // I use 3 var with same value because fortran uses var args
  one(c,n,1); // GLOBAL static mtx c
  zer(wk,n+6*(n*m+1),1);
  elm(wx,n,1,xweights); // was dim() in mtx3
  gcvspl_(ptr(x),ptr(y),ny,ptr(wx),ptr(wy),m,n,k,md,evar,ptr(c),nc,ptr(wk),ier);
  // if verbose then mwr(c,'gcv3.txt');
  gcvspl := (ier=0)
end;

// ===========================================================================

// Ignores xweights but uses mode and valr
// caller should take care that w makes sense with mode and valr...
// intended use: mode=3 valr=1.0 w_j = fact/ese^2 of j^eme data
function gcvspl(x,y,w:mtx; var cc:mtx):boolean;  // *4* mmmM
var
  // valr, mode: interface global var
  n: longint = 0;
  nc: longint = 0;
  ny: longint = 0;
  m: longint = 2;
  k: longint = 1;
  ier: longint = -999;  // should return 0 if all ok or >0 if err. in fortran
begin
  one(wy,1,1);
  n := lgn(x); // TODO lgn(x) or col(x)
  ny := n;
  nc := n; // I use 3 var with same value because fortran uses var args
  one(cc,n,1); // achtung global c used by other gcvspl avatars
  zer(wk,n+6*(n*m+1),1);
  gcvspl_(ptr(x),ptr(y),ny,ptr(w),ptr(wy),m,n,k,mode,valr,ptr(cc),nc,ptr(wk),ier);
  // if verbose then mwr(cc,'gcv4.txt');
  gcvspl := (ier=0)
end;

// ===========================================================================
// calc. first and second derivative...
// I may well use spl sev instead...
// which are in mtx09... and get integral as well...
//
function splder(x: mtx; cc:mtx; u:mtx; var v, d1, d2:mtx):boolean;
var
  i: longint;
  n, nu: longint;
  m: longint = 2; // MUST be same as in gcvspl()
  xi, yi0, yi1, yi2: double;
  ider: longint;
  l: longint = 0;
begin
  splder := false;
  // TODO sanity check on u limits between xmin..xmax
  nu := lgn(u);
  n := lgn(x);
  zer(q,2*m,1); // wrk space for splder_()
  zer(v,nu,1);
  zer(d1,nu,1);
  zer(d2,nu,1);
  for i := 1 to nu do begin
    ider := 0;
    xi := elm(u,i,1);
    yi0  := splder_(ider,m,n,xi,ptr(x),ptr(cc),l,ptr(q));
    // writeln('#splder, i, l, xi, yi0: ',i:5,l:5,xi:10:2,yi0:10:4);
    elm(v,i,1,yi0);
    ider := 1;
    yi1  := splder_(ider,m,n,xi,ptr(x),ptr(cc),l,ptr(q));
    //yi1 := 0;
    elm(d1,i,1,yi1);
    ider := 2;
    yi2  := splder_(ider,m,n,xi,ptr(x),ptr(cc),l,ptr(q));
    //yi2 := 0;
    elm(d2,i,1,yi2)
  end;
  splder := true
end;

end.
