/* splinef.i
 * $Id: splinef.i,v 1.1 2007/03/04 05:13:15 dhmunro Exp $
 * piecewise cubic interpolation functions
 */
/* Copyright (c) 2007, The Regents of the University of California.
 * All rights reserved.
 * This file is part of yorick (http://yorick.sourceforge.net).
 * Read the accompanying LICENSE file for details.
 */

/* This package extends the original spline.i functions, concentrating
 * on how arbitrary piecewise cubic functions may be used as a natural
 * extension of the piecewise linear interp() built-in function.
 * Given a sequence of abcissa values x(i), the corresponding function
 * values y(i) uniquely determine a piecewise linear function.  Knowing
 * both the function and its first derivatives, y(i) and dydx(i),
 * uniquely determines a piecewise cubic function.
 *
 * Unlike a true "spline", a general piecewise cubic curve has a
 * discontinuous second derivative.  True spline curves depend only on
 * the function values y(i) at x(i); the derivatives dydx(i) are
 * computed to force continuity of the second derivative at x(i).
 * The less continuous general piecewise cubic permits you to make a
 * more robust fit, which is still smoother than piecewise linear, at
 * the cost of having to provide dydx(i) values by some other means.
 * The splinelsq function can select y(i) and dydx(i) to have least
 * square residuals of the piecewise cubic from a larger set of data,
 * like the fitlsq function in fitlsq.i does for piecewise linear fits.
 *
 * Note that a pair of piecewise cubic curves is a Bezier curve; that is,
 * if x(t) and y(t) are piecewise cubic, then (x,y) is a (cubic) Bezier.
 * These are the curves generated by the Postscript arcto operator.
 * If (x0,y0) is an initial point, (x1,y1) and (x2,y2) are the Bezier
 * control points, and (x3,y3) is the final point, then the spline
 * with 0<t<1 has dxdt0=3*(x1-x0), dydt0=3*(y1-y0), dxdt3=3*(x3-x2),
 * and dydt3=3*(y3-y2).
 */

func splinef(dydx, y, x, xp)
/* DOCUMENT yp = splinef(dydx, y, x, xp)
 *       or yp = splinef(x_y_dydx, xp)
 *   returns piecewise cubic function specified by DYDX, Y, X at
 *   the points XP.  Extrapolation beyond the extreme endpoints of X
 *   is linear, with slope equal to the final value of DYDX.  The
 *   return value dimensions are the same as the dimensions of XP.
 *
 *   In the second form, X_Y_DYDX is a 3-by-nknots array of [x,y,dydx]
 *   values.  The values of X in either case must either increase or
 *   decrease monotonically.
 *
 * SEE ALSO: interp, splined, splinei, splinelsq
 */
{
  local c0, c1, c2, c3, cm1;
  spline_coef;
  return poly(x, c0, c1, c2, c3);
}

func splined(dydx, y, x, xp)
/* DOCUMENT yp = splined(dydx, y, x, xp)
 *       or yp = splined(x_y_dydx, xp)
 *   returns the derivative of the piecewise cubic function specified
 *   by DYDX, Y, X at the points XP.  Extrapolation beyond the extreme
 *   endpoints of X is linear, so splined gives the final value of DYDX.
 *   The return value dimensions are the same as the dimensions of XP.
 *
 *   In the second form, X_Y_DYDX is a 3-by-nknots array of [x,y,dydx]
 *   values.  The values of X in either case must either increase or
 *   decrease monotonically.
 *
 * SEE ALSO: splined, splinei
 */
{
  local c0, c1, c2, c3, cm1;
  spline_coef;
  return poly(x, c1, c2+c2, 3.*c3);
}

func splinei(dydx, y, x, xp)
/* DOCUMENT yp = splinei(dydx, y, x, xp)
 *       or yp = splinei(x_y_dydx, xp)
 *   returns the integral of the piecewise cubic function specified
 *   by DYDX, Y, X at the points XP.  The integral is quadratic beyond
 *   the extreme endpoints of X, and zero at X(1).  The dimensions of
 *   the return value are the same as the dimensions of XP.
 *   This is the cubic analog of the integ function.
 *
 *   In the second form, X_Y_DYDX is a 3-by-nknots array of [x,y,dydx]
 *   values.  The values of X in either case must either increase or
 *   decrease monotonically.
 *
 * SEE ALSO: integ, splinef, splined
 */
{
  local c0, c1, c2, c3, cm1;
  cm1 = 1;
  spline_coef;
  return poly(x, cm1, c0, 0.5*c1, (1./3.)*c2, 0.25*c3);
}

func spline_coef
/* DOCUMENT spline_coef
 *   is the worker for the splinef, splined, and splinei functions.
 *   If you need to compute both function and derivative or integral,
 *   you will improve performance using spline_coef.  See the source
 *   code for those functions for usage.
 *
 * SEE ALSO: splinef, splined, splinei
 */
{
  extern c0, c1, c2, c3, cm1, dydx, y, x;
  if (is_void(xp)) {
    xp = y;
    x = dydx(1,);
    y = dydx(2,);
    dydx = dydx(3,);
  }
  u = digitize(xp, x);
  l = max(u - 1, 1);
  u = min(u, numberof(x));

  xl = x(l);
  c0 = y(l);
  c1 = dydx(l);

  if (cm1) {
    c2 = x(dif);
    cm1 = ((y(zcen) - (1./12.)*dydx(dif)*c2)*c2)(cum)(l);
  }

  x = x(u) - xl;     /* 0 outside endpoints of x */
  c2 = double(!x);
  x = (1. - c2)/(x + c2);
  c2 = (y(u) - c0) * x;
  c3 = dydx(u) + c1 - c2 - c2;
  c2 = (c2 - c1 - c3) * x;
  c3 *= x * x;

  x = xp - xl;
}

/* question: any simple way to get continuous 2nd derivative fit? */
func splinelsq(y, x, xfit, weight=, y0=, dydx0=, y1=, dydx1=, constrain=)
/* DOCUMENT x_y_dydx = splinelsq(y, x, xfit)
            ...
            yp = splinef(x_y_dydx, xp)
     performs a least squares fit to the data points (X, Y).  The input
     XFIT are the abcissas of the piecewise cubic function with knot
     points XFIT which is the least squares best fit to the data (X,Y).
     The XFIT must be strictly increase or decrease.

     Any points in XFIT with no data points in the intervals on
     either side will be removed.

     A weight= keyword of the same length as X and Y may be supplied in
     order to weight the various data points differently; a typical
     WEIGHT function is 1/sigma^2 where sigma are the standard deviations
     associated with the Y values.

     You can specify y0=, dydx0=, y1=, and dydx1= keywords to fix the
     value of the function or its derivative at the first (0) or last (1)
     endpoint.  Be sure there is at least one point in the final
     interval so that the XFIT at the endpoint is not removed.

     More generally, you can specify a constrain= keyword.  The value
     of constrain is a hook function which will be called just before
     the matrix solve.  Your constrain subroutine will be passed no
     arguments, but it can access and modify the mat and rhs variables.

   SEE ALSO: splinef
 */
{
  xfit = double(xfit);
  np1 = numberof(xfit)+1;

  /* bin the input data into the xfit */
  l = digitize(x, xfit);
  /* remove any points between adjacent empty bins */
  list = histogram(l, top=np1);
  mask = list(2:0) + list(1:-1);
  list = where(mask);
  nfit = numberof(xfit);
  if (numberof(list) < nfit) {
    if (numberof(list) < 2) error, "all data points outside xfit";
    xfit = xfit(list);
    np1 = numberof(xfit)+1;
    l = digitize(x, xfit);
  }
  /* create extended xfit for which l and l+1 may be used directly */
  dx = xfit(0) - xfit(1);
  xl = min(min(x), min(xfit)) - abs(dx);
  xu = max(max(x), max(xfit)) + abs(dx);
  xx = (dx>0.)? grow(xl, xfit, xu) : grow(xu, xfit, xl);

  xl = xx(l);
  xu = xx(l+1);
  dx = xu - xl;
  rdx = 1./dx;

  ld = x-xl;
  lx = ld * rdx;
  ud = lx*lx;
  lx = 1.-lx;
  ld *= lx*lx;
  uf = (1.+lx+lx)*ud;
  ud *= -lx*dx;
  lf = 1.-uf;

  if (is_void(weight)) weight = 1.0;
  lfy = histogram(l, lf*y*weight, top=np1);
  ldy = histogram(l, ld*y*weight, top=np1);
  ufy = histogram(l, uf*y*weight, top=np1);
  udy = histogram(l, ud*y*weight, top=np1);

  lflf = histogram(l, lf*lf*weight, top=np1);
  lfld = histogram(l, lf*ld*weight, top=np1);
  lfuf = histogram(l, lf*uf*weight, top=np1);
  lfud = histogram(l, lf*ud*weight, top=np1);
  ldld = histogram(l, ld*ld*weight, top=np1);
  lduf = histogram(l, ld*uf*weight, top=np1);
  ldud = histogram(l, ld*ud*weight, top=np1);
  ufuf = histogram(l, uf*uf*weight, top=np1);
  ufud = histogram(l, uf*ud*weight, top=np1);
  udud = histogram(l, ud*ud*weight, top=np1);

  n = 2*numberof(xfit);
  rhs = array(0., n);
  mat = array(rhs, n);
  n2 = n*n;
  dn = 2*(n+1);
  mat(1:n2:dn) = lflf(2:0) + ufuf(1:-1);
  mat(2:n2:dn) = mat(n+1:n2:dn) = lfld(2:0) + ufud(1:-1);
  mat(n+2:n2:dn) = ldld(2:0) + udud(1:-1);
  mat(3:n2-n:dn) = mat(dn-1:n2:dn) = lfuf(2:-1);
  mat(4:n2-n:dn) = mat(n+dn-1:n2:dn) = lfud(2:-1);
  mat(n+3:n2:dn) = mat(dn:n2:dn) = lduf(2:-1);
  mat(n+4:n2:dn) = mat(n+dn:n2:dn) = ldud(2:-1);
  rhs(1:n:2) = lfy(2:0) + ufy(1:-1);
  rhs(2:n:2) = ldy(2:0) + udy(1:-1);

  if (!is_void(y0)) {
    if (list(1)!=1) error, "y0= keyword with no data in first interval";
    rhs(1) = y0;
    mat(1,) = 0.0;
    mat(1,1) = 1.0;
  }
  if (!is_void(dydx0)) {
    if (list(1)!=1) error, "dydx0= keyword with no data in first interval";
    rhs(2) = dydx0;
    mat(2,) = 0.0;
    mat(2,2) = 1.0;
  }
  if (!is_void(y1)) {
    if (list(0)!=nfit) error, "y1= keyword with no data in last interval";
    rhs(-1) = y1;
    mat(-1,) = 0.0;
    mat(-1,-1) = 1.0;
  }
  if (!is_void(dydx1)) {
    if (list(0)!=nfit) error, "dydx1= keyword with no data in last interval";
    rhs(0) = dydx1;
    mat(0,) = 0.0;
    mat(0,0) = 1.0;
  }
  if (!is_void(constrain)) constrain;

  fd = LUsolve(mat, rhs);
  xfd = xfit(-:1:3,);
  xfd(2,) = fd(1:0:2);
  xfd(3,) = fd(2:0:2);
  return xfd;
}