{$H+}  // ansi strings
{$C+}  // assertions on
program smo;

uses
  getopts, strutils, mtx09, gcv09;

var
  fi, fo: text;       // input output and parameters ( mbg ese )
  fnam: string = '';  //
  mnam: string = '*'; // input data file name and matrix name
  n: longint;         // number of segments of interp. table
  wfact: double;      // 0 -> unit weights > 0 uses col. 3 as Sy w = wfact/(Sy)
  trace: boolean = false;   // option -t
  bkgsub: boolean = false;   // option -b
  nodeFile: string;   // option -s
  ri, rx, r0: mtx;    // roi include/exclude and null matrices
  md, ud: mtx;        // md def.loop to build ud array
  deb, pas: double;   // see option -d x grid from deb by pas
  ul: mtx;            // list of valid values in ul array

procedure usage;
const
  help={$i hlp.pas}
begin
  writeln;
  writeln(help);
  halt
end;

procedure options;
var
  nw: integer;
  c: char;
  opti: string = '';     // input file
  opto: string = '';     // output file
  optn: string = '0';    // number of points in smoothing vector
  optm: string = '3';    // mode global var see gcvspl.f source OR 3
  optv: string = '1.0';  // valr global var ( val arg of gcvspl_ )
  optf: string = '1.0';  // multiply ese vector ( increase smoothing )
                         // if f=0 then no smoothing same as spline.f
  optx: string = '';     // exclude region from smoothing
  optr: string = '';     // def. empty roi means use all data !
  opts: string = 'smo.xy'; // spline coeff. output default output file
  optd: string = '';     // interp. array: beg step
  optl: string = '';     // interp. list
  optstr: string = 'hi:o:n:m:v:f:r:x:bts:d:l:';
  code: integer;
begin
  c := getopt(optstr);
  while c<>EndOfOptions do begin
    case c of
      'h': usage;
      'i': opti := optarg;
      'o': opto := optarg;
      'n': optn := optarg;
      'm': optm := optarg;
      'v': optv := optarg;
      'f': optf := optarg;
      'r': optr := optarg;
      'x': optx := optarg;
      'b': bkgsub := true;
      't': trace := true;
      's': opts := optarg;
      'd': optd := optarg;
      'l': optl := optarg;
    end;
    c := getopt(optstr);
  end;
// split opti into the file name and the matrix name with default
  nw := WordCount(opti,WordDelim);
  if nw>1 then begin
    mnam := ExtractWord(2,opti,WordDelim);
  end else mnam := '*';
  if nw>0 then begin
    fnam := ExtractWord(1,opti,WordDelim);
  end else fnam := '';
  assign(fi,fnam); reset(fi);
  assign(fo,opto); rewrite(fo);
  val(optn,n,code); assert( code=0 );
  val(optm,gcv09.mode,code); assert( code=0 );
  val(optv,gcv09.valr,code); assert( code=0 );
  val(optf,wfact,code); assert( code=0 );
  if wfact<=0 then begin
    // force simple interpolation with mode=1 and valr=0
    // ignoring any previous values of -m and -v options
    gcv09.mode := 1;
    gcv09.valr := 0.0;
  end;
  equ(ri,optr); // region of interest (smooth and output)
  equ(rx,optx); // region to exclude from smoothing ...
  equ(r0,'');   // null matrix
  equ(ul,optl); // interp. each value in ml mtx
  equ(md,optd); // u = md(1) to md(3) by md(2) ; def. md(3) is x(n)
  nodeFile := opts; // defaults smo.xy
end;

var
  d, dx: mtx; // d: original data dx: sub set after removing -x ranges
  x, y, s, w, c0: mtx;
  u, v, d1, d2, o: mtx;
  s5: mtx; // cubic spline coeff. s5(5,m)  x a b c d returned by spl()
  t5: mtx; // transposed s5 by cpy (plain mtx)
  i0, n0: longint;
  x1, xf: double;

begin
  options;
  mrd(d,fi,mnam);          // TODO use mnam set by opt? ... shift in option ?
  equ(dx,d);               // duplicate matrix d
  rix(d,ri,r0);            // note here rx is empty (optx not used yet)
  rix(dx,ri,rx);           // simple to understand... d all usefull data
  assert( col(d)>=2 );     // opt. 3d col. is e.s.e.
  equ(x,col(dx,1));        // cpy x1 and y1 from d: x1, y1 must have contiguous elem.
  equ(y,col(dx,2));        // whereas col(d,1) elem. are not contiguous
  one(w,x);                // init weight vector : w = ones(x)
  if col(d)>2 then begin
    equ(s,col(dx,3));      // std. dev. use below to calc relative inverse variance
  end else begin
    one(s,x)
  end;
  if (wfact>0) then begin  // use unit weights even if 3 cols. -f 0
    mul(s,wfact);          // if uniform ese "-f ese" and 2 col. data file
                           // should be the same -f ese -v 1 ? and -v ese^2 -f 1.0 ?
  end else one(s,x);
  dvd(w,w,s);              // w = w ./ (s*wfact)
  mul(w,w,w);              // w_j = (1/wfact/s_j)^2
  assert( gcvspl(x,y,w,c0) ); // avatar *4*
  // now do the interpolation as required by -n -d -l options...
  if not mty(md) then begin
    assert((lgn(md)>=1) and (col(md)>=2) ); // only md(1,1) and (1,2) used
    deb := elm(md,1,1);
    pas := elm(md,1,2);
    x1 := elm(x,1,1); // first
    xf := elm(x,lgn(x),1); // last
    assert( pas > 0 );
    while deb < x1 do deb := deb + pas;
    n0 := trunc((xf-deb)/pas);
    zer(ud,n0+1,1);             // TODO ud should it be row or col
    for i0 := 1 to n0+1 do begin  // vector or don't care TODO
      elm(ud,i0,1,deb+pred(i0)*pas);
    end
  end;
  stk(u,ud); // interpolation values calculated from options
  stk(u,ul); // -l and -d
  if (n=0) and mty(u) then begin
    equ(u,col(d,1)); // def. smooth eval on nodes
  end;
  if (n>0) then begin // -l and -d ignored if -n option is given
    // eval on a regularly spaced grid n segments i.e. n+1 values.
    // TODO there is a lsp avatar which may be nicer to use here
    lsp(u,elm(d,1,1),elm(d,lgn(d),1),n); // interpol. values linspace()
  end;
// TODO if we use spline/seval do we need splder ? in this case ?
// NOTE second derivative is useless... unless smo factor is high...
  splder(x,c0,u,v,d1,d2);  // v_j : smoothed interp on u_j values d1 and d2 1st & 2d deriv.
// calculate using fmm's spline.f the 3d degree spline coeff. table: x a b c d
// on the smoothed vector c0
  spl(s5,x,c0);      // use fmm spline to retrieve a b c d coeff. for smoothed nodes
  equ(t5,trp(s5));   //
  mwr(t5,nodeFile);      // spline coefficients output
//
  zer(o,lgn(u),5);   // build a table for output col 5 for data - bkg
  elm(col(o,1),u);   // fill table with calc. vectors
  elm(col(o,2),v);
  elm(col(o,3),d1);
  elm(col(o,4),d2);
// interim do substraction here if bkgsub ( see option -b )
// TODO assert option -x ??? to skip peaks area
// TODO exclude option -d -l or -n
  if bkgsub then begin
    elm(col(o,5),col(d,2));
    sub(col(o,5),col(o,2))
  end;
  mwr(o,fo);         // output
  flush(fo)
end.
