#  File src/library/splines/R/splineClasses.R
#  Part of the R package, http://www.R-project.org
#  Copyright (C) 1998 Douglas M. Bates and William N. Venables.
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/

#### Classes and methods for determining and manipulating interpolation
#### splines.

### Major classes:
###   spline - a virtual class of representations of piecewise
###	       polynomial functions.  The join points of the polynomials
###	       are called "knots".  The order of the spline is the number
###	       of coefficients in the polynomial pieces.
###   bSpline - splines represented as linear combinations of B-splines
###   polySpline - splines represented as polynomials

### Minor classes:
###   nbSpline - "natural" bSplines.  That is, splines of order 4 with linear
###		 extrapolation beyond the limits of the knots.
###   npolySpline - polynomial representation of a natural spline
###   pbSpline - periodic bSplines
###   ppolySpline - periodic polynomial splines
###   backSpline - "splines" for inverse interpolation

splineDesign <-
    ## Creates the "design matrix" for a collection of B-splines.
    function(knots, x, ord = 4, derivs = integer(nx), outer.ok = FALSE,
             sparse = FALSE)
{
    knots <- sort(as.numeric(knots))
    if((nk <- length(knots)) <= 0) stop("must have at least 'ord' knots")
    x <- as.numeric(x)
    nx <- length(x)
    if(length(derivs) != nx)
	stop("length of 'derivs' must match length of 'x'")
    ord <- as.integer(ord)
    if(ord > nk || ord < 1)
	stop("'ord' must be positive integer, at most the number of knots")

    ## The x test w/ sorted knots assumes ord <= nk+1-ord, or nk >= 2*ord-1L:
    if(!outer.ok && nk < 2*ord-1)
        stop(gettextf("need at least %s (=%d) knots",
                      "2*ord -1", 2*ord -1),
             domain = NA)

    o1 <- ord - 1L
### FIXME: the 'outer.ok && need.outer' handling would more efficiently happen
###        in the underlying C code - with some programming effort though..
    if(need.outer <- any(out.x <- x < knots[ord] | knots[nk - o1] < x)) {
        if(outer.ok) { ## x[] is allowed to be 'anywhere'
            ## extend knots set "temporarily"
            in.x <- knots[1L] < x & x < knots[nk]
	    knots <- knots[c(rep.int(1L, o1), 1L:nk, rep.int(nk, o1))]
            if((x.out <- !all(in.x))) {
                x <- x[in.x]
                nnx <- length(x)
            }
	} else
	stop(gettextf("the 'x' data must be in the range %g to %g unless you set '%s'",
		      knots[ord],
                      knots[nk - o1],
                      "outer.ok = TRUE"),
             domain = NA)
    }
    temp <- .Call(C_spline_basis, knots, ord, x, derivs)
    ncoef <- nk - ord

    ii <- if(need.outer && x.out) { # only assign non-zero for x[]'s "inside" knots
        rep.int((1L:nx)[in.x], rep.int(ord, nnx))
    } else rep.int(1L:nx, rep.int(ord, nx))
    jj <- c(outer(1L:ord, attr(temp, "Offsets"), "+"))
    ## stopifnot(length(ii) == length(jj))

    if(sparse) {
	if(is.null(tryCatch(loadNamespace("Matrix"), error = function(e)NULL)))
	    stop(gettextf("%s needs package 'Matrix' correctly installed",
                          "splineDesign(*, sparse=TRUE)"),
                 domain = NA)

	if(need.outer) { ## shift column numbers and drop those "outside"
	    jj <- jj - o1 - 1L
	    ok <- 0 <= jj & jj < ncoef
	    as(new("dgTMatrix",
		   i = ii[ok] - 1L,
		   j = jj[ok],
		   x = as.double(temp[ok]),
		   Dim = c(nx, ncoef)), "CsparseMatrix")
	}
	else
	    as(new("dgTMatrix",
		   i = ii - 1L,
		   j = jj - 1L,
		   x = as.double(temp),
		   Dim = c(nx, ncoef)), "CsparseMatrix")

    } else { ## traditional (dense) matrix
	design <- matrix(double(nx * ncoef), nx, ncoef)
	if(need.outer) { ## shift column numbers and drop those "outside"
	    jj <- jj - o1
	    ok <- 1 <= jj & jj <= ncoef
	    design[cbind(ii, jj)[ok , , drop=FALSE]] <- temp[ok]
	}
	else
	    design[cbind(ii, jj)] <- temp
	design
    }
}

interpSpline <-
    ## Determine the natural interpolation spline.
    function(obj1, obj2, bSpline = FALSE, period = NULL, na.action = na.fail)
    UseMethod("interpSpline")

interpSpline.default <-
    function(obj1, obj2, bSpline = FALSE, period = NULL, na.action = na.fail)
{
    ord <- 4 # spline order -- future: want to use other 'order'/'degree'!
    deg <- ord - 1

    frm <- na.action(data.frame(x = as.numeric(obj1), y = as.numeric(obj2)))
    frm <- frm[order(frm$x), ]
    ndat <- nrow(frm)
    x <- frm$x
    if(anyDuplicated(x))
	stop("values of 'x' must be distinct")
    ## 'deg' extra knots (shifted) out on each side :
    knots <- c(x[1L:deg] + x[1L] - x[ord], x,
	       x[ndat + (1L:deg) - deg] + x[ndat] - x[ndat - deg])
    derivs <- c(2, integer(ndat), 2) # 2nd derivs coerced to 0 in solve() below
    x	   <- c(x[1L], x, x[ndat])
## Solving the system of equations for the spline coefficients can be
## simplified by using banded matrices but the required Linpack routines
## are not loaded as part of S.
##  z <- .C("spline_basis",
##	as.double(knots),
##	as.integer(length(knots) - 4),
##	as.integer(4),
##	as.double(x),
##	as.integer(derivs),
##	as.integer(ndat + 2),
##	design = array(0, c(4, ndat + 2)),
##	offsets = integer(ndat + 2))
##  abd <- array(0, c(7, ndat + 2))
##  abd[4:7, 2:ndat] <- z$design[, 2:ndat]
##  abd[5:7, 1] <- z$design[-4, 1]
##  abd[4:6, ndat + 1] <- z$design[-1, ndat + 1]
##  abd[3:5, ndat + 2] <- z$design[-1, ndat + 2]
##  z <- .Fortran("dgbfa",
##	abd = abd,
##	lda = as.integer(7),
##	n = as.integer(ndat + 2),
##	ml = 2L,
##	mu = 2L,
##	ipvt = integer(ndat + 2),
##	info = integer(1))
##  zz <- .Fortran("dgbsl",
##	abd = z$abd,
##	lda = z$lda,
##	n = z$n,
##	ml = z$ml,
##	mu = z$mu,
##	ipvt = z$ipvt,
##	b = c(0, y, 0),
##	job = 1L)
    des <- splineDesign(knots, x, ord, derivs)
    coeff <- solve(as.matrix(des), c(0, frm$y, 0))
    value <- list(knots = knots, coefficients = coeff, order = ord)
    attr(value, "formula") <-
	do.call("~", list(substitute(obj2), substitute(obj1)))
    class(value) <- c("nbSpline", "bSpline", "spline")
    if (bSpline) return(value)
    ## else convert from B- to poly-Spline:
    value <- polySpline(value)
    coeff <- coef(value)
    coeff[ , 1] <- frm$y
    coeff[1, deg] <- coeff[nrow(coeff), deg] <- 0
    value$coefficients <- coeff
    value
}

interpSpline.formula <-
    function(obj1, obj2, bSpline = FALSE, period = NULL, na.action = na.fail)
{
    form <- as.formula(obj1)
    if (length(form) != 3)
	stop("'formula' must be of the form \"y ~ x\"")
    local <- if (missing(obj2)) sys.parent(1) else as.data.frame(obj2)
    value <- interpSpline(as.numeric(eval(form[[3L]], local)),
			  as.numeric(eval(form[[2L]], local)),
			  bSpline = bSpline, period = period,
			  na.action = na.action)
    attr(value, "formula") <- form
    value
}

periodicSpline <-
    ## Determine the periodic interpolation spline.
    function(obj1, obj2, knots, period = 2 * pi, ord = 4)
    UseMethod("periodicSpline")

periodicSpline.default <-
    function(obj1, obj2, knots, period = 2 * pi, ord = 4)
{
    x <- as.numeric(obj1)
    y <- as.numeric(obj2)
    lenx <- length(x)
    if(lenx != length(y))
	stop("lengths of 'x' and 'y' must match")
    ind <- order(x)
    x <- x[ind]
    if(length(unique(x)) != lenx)
	stop("values of 'x' must be distinct")
    if(any((x[-1L] - x[ - lenx]) <= 0))
	stop("values of 'x' must be strictly increasing")
    if(ord < 2) stop("'ord' must be >= 2")
    o1 <- ord - 1
    if(!missing(knots)) {
	period <- knots[length(knots) - o1] - knots[1L]
    }
    else {
	knots <- c(x[(lenx - (ord - 2)):lenx] - period, x, x[1L:ord] + period)
    }
    if((x[lenx] - x[1L]) >= period)
	stop("the range of 'x' values exceeds one period")
    y <- y[ind]
    coeff.mat <- splineDesign(knots, x, ord)
    sys.mat <- coeff.mat[, (1L:lenx)]
    sys.mat[, seq_len(o1)] <- sys.mat[, seq_len(o1)] +
	coeff.mat[, lenx + seq_len(o1)]
    coeff <- qr.coef(qr(sys.mat), y)
    coeff <- c(coeff, coeff[seq_len(o1)])
    value <- list(knots = knots, coefficients = coeff, order = ord,
		  period = period)
    attr(value, "formula") <- do.call("~", as.list(sys.call())[3:2])
    class(value) <- c("pbSpline", "bSpline", "spline")
    value
}

periodicSpline.formula <- function(obj1, obj2, knots, period = 2 * pi, ord = 4)
{
    form <- as.formula(obj1)
    if (length(form) != 3)
	stop("'formula' must be of the form \"y ~ x\"")
    local <- if (missing(obj2)) sys.parent(1) else as.data.frame(obj2)
    ## 'missing(knots)' is transfered :
    value <-  periodicSpline(as.numeric(eval(form[[3L]], local)),
			     as.numeric(eval(form[[2L]], local)),
			     knots = knots, period = period, ord = ord)
    attr(value, "formula") <- form
    value
}

polySpline <-
    ## Constructor for polynomial representation of splines
    function(object, ...) UseMethod("polySpline")

polySpline.polySpline <- function(object, ...) object

as.polySpline <-
    ## Conversion of an object to a polynomial spline representation
    function(object, ...) polySpline(object, ...)

polySpline.bSpline <- function(object, ...)
{
    ord <- splineOrder(object)
    knots <- splineKnots(object)
    if(is.unsorted(knots))
	stop("knot positions must be non-decreasing")
    knots <- knots[ord:(length(knots) + 1L - ord)]
    coeff <- array(0, c(length(knots), ord))
    coeff[, 1] <- asVector(predict(object, knots))
    if(ord > 1) {
	for(i in 2:ord) {
	    coeff[, i] <- asVector(predict(object, knots, deriv = i - 1))/
		prod(1L:(i - 1))
	}
    }
    value <- list(knots = knots, coefficients = coeff)
    attr(value, "formula") <- attr(object, "formula")
    class(value) <- c("polySpline", "spline")
    value
}

polySpline.nbSpline <- function(object, ...)
{
    value <- NextMethod("polySpline")
    class(value) <- c("npolySpline", "polySpline", "spline")
    value
}

polySpline.pbSpline <- function(object, ...)
{
    value <- NextMethod("polySpline")
    value[["period"]] <- object$period
    class(value) <- c("ppolySpline", "polySpline", "spline")
    value
}

## A couple of accessor functions for the virtual class of splines.

splineKnots <- ## Extract the knot positions
    function(object) UseMethod("splineKnots")

splineKnots.spline <- function(object) object$knots

splineOrder <- ## Extract the order of the spline
    function(object) UseMethod("splineOrder")

splineOrder.bSpline <- function(object) object$order

splineOrder.polySpline <- function(object) ncol(coef(object))

## xyVector is a class of numeric vectors that represent responses and
## carry with them the corresponding inputs x.	Very similar in purpose
## to the "track" class in JMC's book.  All methods for predict
## applied to spline objects produce such objects as their value.

xyVector <- ## Constructor for the xyVector class
    function(x, y)
{
    x <- as.vector(x)
    y <- as.vector(y)
    if(length(x) != length(y)) {
	stop("lengths of 'x' and 'y' must be the same")
    }
    structure(list(x = x, y = y), class = "xyVector")
}

asVector <- ## coerce object to a vector.
    function(object) UseMethod("asVector")

asVector.xyVector <- function(object) object$y

as.data.frame.xyVector <- function(x, ...) data.frame(x = x$x, y = x$y)

plot.xyVector <- function(x, ...)
{
    plot(x = x$x, y = x$y, ...)
###  xyplot(y ~ x, as.data.frame(x), ...)
}

predict.polySpline <- function(object, x, nseg = 50, deriv = 0, ...)
{
    knots <- splineKnots(object)
    coeff <- coef(object)
    cdim <- dim(coeff)
    ord <- cdim[2L]
    if(missing(x))
	x <- seq.int(knots[1L], knots[cdim[1L]], length.out = nseg + 1)
    i <- as.numeric(cut(x, knots))
    i[x == knots[1L]] <- 1
    delx <- x - knots[i]
    deriv <- as.integer(deriv)[1L]
    if(deriv < 0 || deriv >= ord)
	stop(gettextf("'deriv' must be between 0 and %d", ord - 1),
             domain = NA)
    while(deriv > 0) {
	ord <- ord - 1
	coeff <- t(t(coeff[, -1]) * (1L:ord))
	deriv <- deriv - 1
    }
    y <- coeff[i, ord]
    if(ord > 1) {
	for(j in (ord - 1):1)
	    y <- y * delx + coeff[i, j]
    }
    xyVector(x = x, y = y)
}

predict.bSpline <- function(object, x, nseg = 50, deriv = 0, ...)
{
    knots <- splineKnots(object)
    if(is.unsorted(knots))
	stop("knot positions must be non-decreasing")
    ord <- splineOrder(object)
    if(deriv < 0 || deriv >= ord)
	stop(gettextf("'deriv' must be between 0 and %d", ord - 1),
             domain = NA)
    ncoeff <- length(coeff <- coef(object))
    if(missing(x)) {
	x <- seq.int(knots[ord], knots[ncoeff + 1], length.out = nseg + 1)
	accept <- TRUE
    } else accept <- knots[ord] <= x & x <= knots[ncoeff + 1]
    y <- x
    y[!accept] <- NA
    y[accept] <- .Call(C_spline_value, knots, coeff, ord, x[accept], deriv)
    xyVector(x = x, y = y)
}

predict.nbSpline <- function(object, x, nseg = 50, deriv = 0, ...)
{
    value <- NextMethod("predict")
    if(!any(is.na(value$y))) # when x were inside knots
	return(value)

    x <- value$x
    y <- value$y
    ## Compute y[] for x[] outside knots:
    knots <- splineKnots(object)
    ord <- splineOrder(object)
    ncoeff <- length(coef(object))
    bKn <- knots[c(ord,ncoeff + 1)]
    coeff <- array(0, c(2L, ord))
    ## Extrapolate using a + b*(x - boundary.knot) (ord=4 specific?)
    coeff[,1] <- asVector(predict(object, bKn))
    coeff[,2] <- asVector(predict(object, bKn, deriv = 1))
    deriv <- as.integer(deriv)## deriv < ord already tested in NextMethod()
    while(deriv) {
	ord <- ord - 1
	## could be simplified when coeff has <= 2 non-zero cols:
	coeff <- t(t(coeff[, -1]) * (1L:ord))# 1L:ord = the 'k' in k* x^{k-1}
	deriv <- deriv - 1
    }
    nc <- ncol(coeff)
    if(any(which <- (x < bKn[1L]) & is.na(y)))
	y[which] <- if(nc==0) 0 else if(nc==1) coeff[1, 1]
	else coeff[1, 1] + coeff[1, 2] * (x[which] - bKn[1L])
    if(any(which <- (x > bKn[2L]) & is.na(y)))
	y[which] <- if(nc==0) 0 else if(nc==1) coeff[1, 1]
	else coeff[2, 1] + coeff[2, 2] * (x[which] - bKn[2L])
    xyVector(x = x, y = y)
}

predict.pbSpline <- function(object, x, nseg = 50, deriv = 0, ...)
{
    knots <- splineKnots(object)
    ord <- splineOrder(object)
    period <- object$period
    ncoeff <- length(coef(object))
    if(missing(x))
	x <- seq.int(knots[ord], knots[ord] + period, length.out = nseg + 1)
    x.original <- x
    if(any(ind <- x < knots[ord]))
	x[ind] <- x[ind] + period * (1 + (knots[ord] - x[ind]) %/% period)
    if(any(ind <- x > knots[ncoeff + 1]))
	x[ind] <- x[ind] - period * (1 + (x[ind] - knots[ncoeff +1]) %/% period)
    xyVector(x = x.original,
	     y = .Call(C_spline_value, knots, coef(object), ord, x, deriv))
}

predict.npolySpline <- function(object, x, nseg = 50, deriv = 0, ...)
{
    value <- NextMethod()
    if(!any(is.na(value$y))) # when x were inside knots
	return(value)

    x <- value$x
    y <- value$y
    ## Compute y[] for x[] outside knots:
    nk <- length(knots <- splineKnots(object))
    coeff <- coef(object)[ - (2:(nk - 1)), ] # only need col 1L:2
    ord <- dim(coeff)[2L]
    if(ord >= 3) coeff[, 3:ord] <- 0
    deriv <- as.integer(deriv)
    while(deriv) {
	ord <- ord - 1
	## could be simplified when coeff has <= 2 non-zero cols:
	coeff <- t(t(coeff[, -1]) * (1L:ord))# 1L:ord = the 'k' in k* x^{k-1}
	deriv <- deriv - 1
    }
    nc <- ncol(coeff)
    if(any(which <- (x < knots[1L]) & is.na(y)))
	y[which] <- if(nc==0) 0 else if(nc==1) coeff[1, 1]
	else coeff[1, 1] + coeff[1, 2] * (x[which] - knots[1L])
    if(any(which <- (x > knots[nk]) & is.na(y)))
	y[which] <- if(nc==0) 0 else if(nc==1) coeff[1, 1]
	else coeff[2, 1] + coeff[2, 2] * (x[which] - knots[nk])

    xyVector(x = x, y = y)
}

predict.ppolySpline <- function(object, x, nseg = 50, deriv = 0, ...)
{
    knots <- splineKnots(object)
    nknot <- length(knots)
    period <- object$period
    if(missing(x))
	x <- seq.int(knots[1L], knots[1L] + period, length.out = nseg + 1)
    x.original <- x

    if(any(ind <- x < knots[1L]))
	x[ind] <- x[ind] + period * (1 + (knots[1L] - x[ind]) %/% period)
    if(any(ind <- x > knots[nknot]))
	x[ind] <- x[ind] - period * (1 + (x[ind] - knots[nknot]) %/% period)

    value <- NextMethod("predict")
    value$x <- x.original
    value
}

## The plot method for all spline objects

plot.spline <- function(x, ...)
{
###  args <- list(formula = y ~ x, data = as.data.frame(predict(x)), type = "l")
    args <- list(x = as.data.frame(predict(x)), type = "l")
    if(length(form <- attr(x, "formula")) == 3) {
	args <- c(args, list(xlab = deparse(form[[3L]]), ylab = deparse(form[[2L]])))
    }
    args <- c(list(...), args)
###  do.call("xyplot", args)
    do.call("plot", args[unique(names(args))])
}

print.polySpline <- function(x, ...)
{
    coeff <- coef(x)
    dnames <- dimnames(coeff)
    if (is.null(dnames[[2L]]))
	dimnames(coeff) <-
	    list(format(splineKnots(x)),
		 c("constant", "linear", "quadratic", "cubic",
		   paste0(4:29, "th"))[1L:(dim(coeff)[2L])])
    cat("polynomial representation of spline")
    if (!is.null(form <- attr(x, "formula")))
	cat(" for", deparse(as.vector(form)))
    cat("\n")
    print(coeff, ...)
    invisible(x)
}

print.ppolySpline <- function(x, ...)
{
    cat("periodic ")
    value <- NextMethod("print")
    cat("\nPeriod:", format(x[["period"]]), "\n")
    value
}

print.bSpline <- function(x, ...)
{
    value <- c(rep(NA, splineOrder(x)), coef(x))
    names(value) <- format(splineKnots(x), digits = 5)
    cat("bSpline representation of spline")
    if (!is.null(form <- attr(x, "formula")))
	cat(" for", deparse(as.vector(form)))
    cat("\n")
    print(value, ...)
    invisible(x)
}

## backSpline - a class of monotone inverses to an interpolating spline.
## Used mostly for the inverse of the signed square root profile function.

backSpline <- function(object) UseMethod("backSpline")

backSpline.npolySpline <- function(object)
{
    knots <- splineKnots(object)
    nk <- length(knots)
    nkm1 <- nk - 1
    kdiff <- diff(knots)
    if(any(kdiff <= 0))
	stop("knot positions must be strictly increasing")
    coeff <- coef(object)
    bknots <- coeff[, 1]
    adiff <- diff(bknots)
    if(all(adiff < 0))
	revKnots <- TRUE
    else if(all(adiff > 0))
	revKnots <- FALSE
    else
	stop("spline must be monotone")
    bcoeff <- array(0, dim(coeff))
    bcoeff[, 1] <- knots
    bcoeff[, 2] <- 1/coeff[, 2]
    a <- array(c(adiff^2, 2 * adiff, adiff^3, 3 * adiff^2),
	       c(nkm1, 2L, 2L))
    b <- array(c(kdiff - adiff * bcoeff[ - nk, 2L],
		 bcoeff[-1L, 2L] - bcoeff[ - nk, 2L]), c(nkm1, 2))
    for(i in 1L:(nkm1))
	bcoeff[i, 3L:4L] <- solve(a[i,, ], b[i,  ])
    bcoeff[nk, 2L:4L] <- NA
    if(nk > 2L) {
	bcoeff[1L, 4L] <- bcoeff[nkm1, 4L] <- 0
	bcoeff[1L, 2L:3L] <- solve(array(c(adiff[1L], 1, adiff[1L]^2,
                                           2 * adiff[1L]), c(2L, 2L)),
                                   c(kdiff[1L], 1/coeff[2L, 2L]))
	bcoeff[nkm1, 3L] <- (kdiff[nkm1] - adiff[nkm1] *
			    bcoeff[nkm1, 2L])/adiff[nkm1]^2
    }
    if(bcoeff[1L, 3L] > 0) {
	bcoeff[1L, 3L] <- 0
	bcoeff[1L, 2L] <- kdiff[1L]/adiff[1L]
    }
    if(bcoeff[nkm1, 3L] < 0) {
	bcoeff[nkm1, 3L] <- 0
	bcoeff[nkm1, 2L] <- kdiff[nkm1]/adiff[nkm1]
    }
    value <- if(!revKnots) list(knots = bknots, coefficients = bcoeff)
    else {
	ikn <- length(bknots):1L
	list(knots = bknots[ikn], coefficients = bcoeff[ikn,])
    }
    attr(value, "formula") <- do.call("~", as.list(attr(object, "formula"))[3L:2L])
    class(value) <- c("polySpline", "spline")
    value
}

backSpline.nbSpline <- function(object) backSpline(polySpline(object))