src/library/splines/R/splineClasses.R - R - Git at Google

 #  File src/library/splines/R/splineClasses.R
 #  Part of the R package, https://www.R-project.org
 #  Copyright (C) 2000-2018 The R Core Team
 #  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
 #  https://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 = 4L, derivs = 0L, outer.ok = FALSE,
              sparse = FALSE)
 {
     if((nk <- length(knots <- as.numeric(knots))) <= 0)
         stop("must have at least 'ord' knots")
     if(is.unsorted(knots)) knots <- sort.int(knots)
     x <- as.numeric(x)
     nx <- length(x)
     ## derivs is re-cycled to length(x) in C
     if(length(derivs) > nx)
 	stop("length of 'derivs' is larger than length of 'x'")
     if(length(derivs) < 1L) stop("empty 'derivs'")
     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)

     degree <- 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(x < knots[ord] | knots[nk - degree] < x)) {
         if(outer.ok) { ## x[] is allowed to be 'anywhere'
 	    in.x <- knots[1L] <= x & x <= knots[nk]
 	    if((x.out <- !all(in.x))) {
 		x <- x[in.x]
 		nnx <- length(x)
 	    }
 	    ## extend knots set "temporarily": the boundary knots must be repeated >= 'ord' times.
             ## NB: If these are already repeated originally, then, on the *right* only, we need
             ##    to make sure not to add more than needed
             dkn <- diff(knots)[(nk-1L):1] # >= 0, since they are sorted
 	    knots <- knots[c(rep.int(1L, degree),
                              seq_len(nk),
                              rep.int(nk, max(0L, ord - match(TRUE, dkn > 0))))]
 	} else
 	    stop(gettextf("the 'x' data must be in the range %g to %g unless you set '%s'",
 			  knots[ord], knots[nk - degree], "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 - degree - 1L
 	    ok <- 0 <= jj & jj < ncoef
 	    methods::as(methods::new("dgTMatrix", i = ii[ok] - 1L, j = jj[ok],
 				     x = as.double(temp[ok]), # vector, not matrix
 				     Dim = c(nx, ncoef)), "CsparseMatrix")
 	}
 	else
 	    methods::as(methods::new("dgTMatrix", i = ii - 1L, j = jj - 1L,
 				     x = as.double(temp), # vector
 				     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 - degree
 	    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, ord = 4L,
              na.action = na.fail, sparse = FALSE)
     UseMethod("interpSpline")

 ##>>> FIXME: any  ord != 4 needs adaption in splineDesign(), i.e.,
 ##>>> =====       --------  probably in spline_basis() in ../src/splines.c
 interpSpline.default <-
     function(obj1, obj2, bSpline = FALSE, period = NULL,
              ord = 4L, na.action = na.fail, sparse = FALSE)
 {
     ## spline order 'ord' == 'degree' + 1
     stopifnot(exprs = {
         (degree <- ord - 1L) >= 0
         length(degree) == 1L
         degree == as.integer(degree)
     })
     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")
     if(length(x) < ord)
         stop(gettextf("must have at least 'ord'=%d points", ord), domain=NA)
     ## 'degree' extra knots (shifted) out on each side :
     iDeg <- seq_len(degree)
     knots <- c(x[iDeg] + x[1L] - x[ord], x,
 	       x[ndat + iDeg - degree] + x[ndat] - x[ndat - degree])
     if(even <- (ord %% 2 == 0)) { ## natural boundary conditions:
 	nu <- ord %/% 2L          ## nu-th derivs coerced to 0  [in solve() below]
 	derivs <- c(nu, integer(ndat),  nu)
 	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, sparse=sparse)
     y <- c(0, frm$y, 0)
     coeff <- if(sparse) Matrix::solve(des, Matrix::..2dge(y), sparse=TRUE)
              else solve(des, y)
     value <- structure(
         list(knots = knots, coefficients = coeff, order = ord),
              formula = do.call("~", list(substitute(obj2), substitute(obj1))),
         class = c("nbSpline", "bSpline", "spline"))
     if (bSpline) return(value)
     ## else convert from B- to poly-Spline:
     value <- polySpline(value)
     coeff <- coef(value)
     coeff[ , 1L] <- frm$y
     coeff[1L, degree] <- coeff[nrow(coeff), degree] <- 0
     value$coefficients <- coeff
     value
 }

 interpSpline.formula <-
     function(obj1, obj2, bSpline = FALSE, period = NULL,
              ord = 4L, na.action = na.fail, sparse = FALSE)
 {
     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, ord=ord,
 			  na.action=na.action, sparse=sparse)
     attr(value, "formula") <- form
     value
 }

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

 periodicSpline.default <-
     function(obj1, obj2, knots, period = 2 * pi, ord = 4L)
 {
     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")
     degree <- ord - 1L
     if(!missing(knots)) {
 	period <- knots[length(knots) - degree] - 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)
     i1 <- seq_len(degree)
     sys.mat <- coeff.mat[, (1L:lenx)]
     sys.mat[, i1] <- sys.mat[, i1] + coeff.mat[, lenx + i1]
     coeff <- qr.coef(qr(sys.mat), y)
     coeff <- c(coeff, coeff[i1])
     structure(list(knots = knots, coefficients = coeff,
 		   order = ord, period = period),
 	      formula = do.call("~", as.list(sys.call())[3:2]),
 	      class = c("pbSpline", "bSpline", "spline"))
 }

 periodicSpline.formula <- function(obj1, obj2, knots, period = 2 * pi, ord = 4L)
 {
     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 :
     structure(periodicSpline.default(as.numeric(eval(form[[3L]], local)),
 				     as.numeric(eval(form[[2L]], local)),
 				     knots = knots, period = period, ord = ord),
 	      formula = form)
 }

 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 - 1L))/
 		prod(seq_len(i - 1L))
 	}
     }
     structure(list(knots = knots, coefficients = coeff),
 	      formula = attr(object, "formula"),
 	      class = c("polySpline", "spline"))
 }

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

 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 - 1L
 	coeff <- t(t(coeff[, -1]) * seq_len(ord))
 	deriv <- deriv - 1L
     }
     y <- coeff[i, ord]
     if(ord > 1) {
 	for(j in (ord - 1L):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 # (C_spline_value's FIXME)
     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") # predict.bSpline() -> NaN outside knots
     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)
     ## Because of C_spline_value's FIXME, we move the outside-knots x[] values:
     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() # typically predict.polySpline()
     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)
     if((ord <- ncol(coeff)) != 4)
 	stop("currently implemented only for cubic splines")
     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))
	# File src/library/splines/R/splineClasses.R
	# Part of the R package, https://www.R-project.org
	# Copyright (C) 2000-2018 The R Core Team
	# 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
	# https://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 = 4L, derivs = 0L, outer.ok = FALSE,
	sparse = FALSE)
	{
	if((nk <- length(knots <- as.numeric(knots))) <= 0)
	stop("must have at least 'ord' knots")
	if(is.unsorted(knots)) knots <- sort.int(knots)
	x <- as.numeric(x)
	nx <- length(x)
	## derivs is re-cycled to length(x) in C
	if(length(derivs) > nx)
	stop("length of 'derivs' is larger than length of 'x'")
	if(length(derivs) < 1L) stop("empty 'derivs'")
	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",
	"2ord -1", 2ord -1),
	domain = NA)

	degree <- 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(x < knots[ord] \| knots[nk - degree] < x)) {
	if(outer.ok) { ## x[] is allowed to be 'anywhere'
	in.x <- knots[1L] <= x & x <= knots[nk]
	if((x.out <- !all(in.x))) {
	x <- x[in.x]
	nnx <- length(x)
	}
	## extend knots set "temporarily": the boundary knots must be repeated >= 'ord' times.
	## NB: If these are already repeated originally, then, on the right only, we need
	## to make sure not to add more than needed
	dkn <- diff(knots)[(nk-1L):1] # >= 0, since they are sorted
	knots <- knots[c(rep.int(1L, degree),
	seq_len(nk),
	rep.int(nk, max(0L, ord - match(TRUE, dkn > 0))))]
	} else
	stop(gettextf("the 'x' data must be in the range %g to %g unless you set '%s'",
	knots[ord], knots[nk - degree], "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 - degree - 1L
	ok <- 0 <= jj & jj < ncoef
	methods::as(methods::new("dgTMatrix", i = ii[ok] - 1L, j = jj[ok],
	x = as.double(temp[ok]), # vector, not matrix
	Dim = c(nx, ncoef)), "CsparseMatrix")
	}
	else
	methods::as(methods::new("dgTMatrix", i = ii - 1L, j = jj - 1L,
	x = as.double(temp), # vector
	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 - degree
	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, ord = 4L,
	na.action = na.fail, sparse = FALSE)
	UseMethod("interpSpline")

	##>>> FIXME: any ord != 4 needs adaption in splineDesign(), i.e.,
	##>>> ===== -------- probably in spline_basis() in ../src/splines.c
	interpSpline.default <-
	function(obj1, obj2, bSpline = FALSE, period = NULL,
	ord = 4L, na.action = na.fail, sparse = FALSE)
	{
	## spline order 'ord' == 'degree' + 1
	stopifnot(exprs = {
	(degree <- ord - 1L) >= 0
	length(degree) == 1L
	degree == as.integer(degree)
	})
	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")
	if(length(x) < ord)
	stop(gettextf("must have at least 'ord'=%d points", ord), domain=NA)
	## 'degree' extra knots (shifted) out on each side :
	iDeg <- seq_len(degree)
	knots <- c(x[iDeg] + x[1L] - x[ord], x,
	x[ndat + iDeg - degree] + x[ndat] - x[ndat - degree])
	if(even <- (ord %% 2 == 0)) { ## natural boundary conditions:
	nu <- ord %/% 2L ## nu-th derivs coerced to 0 [in solve() below]
	derivs <- c(nu, integer(ndat), nu)
	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, sparse=sparse)
	y <- c(0, frm$y, 0)
	coeff <- if(sparse) Matrix::solve(des, Matrix::..2dge(y), sparse=TRUE)
	else solve(des, y)
	value <- structure(
	list(knots = knots, coefficients = coeff, order = ord),
	formula = do.call("~", list(substitute(obj2), substitute(obj1))),
	class = c("nbSpline", "bSpline", "spline"))
	if (bSpline) return(value)
	## else convert from B- to poly-Spline:
	value <- polySpline(value)
	coeff <- coef(value)
	coeff[ , 1L] <- frm$y
	coeff[1L, degree] <- coeff[nrow(coeff), degree] <- 0
	value$coefficients <- coeff
	value
	}

	interpSpline.formula <-
	function(obj1, obj2, bSpline = FALSE, period = NULL,
	ord = 4L, na.action = na.fail, sparse = FALSE)
	{
	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, ord=ord,
	na.action=na.action, sparse=sparse)
	attr(value, "formula") <- form
	value
	}

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

	periodicSpline.default <-
	function(obj1, obj2, knots, period = 2 * pi, ord = 4L)
	{
	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")
	degree <- ord - 1L
	if(!missing(knots)) {
	period <- knots[length(knots) - degree] - 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)
	i1 <- seq_len(degree)
	sys.mat <- coeff.mat[, (1L:lenx)]
	sys.mat[, i1] <- sys.mat[, i1] + coeff.mat[, lenx + i1]
	coeff <- qr.coef(qr(sys.mat), y)
	coeff <- c(coeff, coeff[i1])
	structure(list(knots = knots, coefficients = coeff,
	order = ord, period = period),
	formula = do.call("~", as.list(sys.call())[3:2]),
	class = c("pbSpline", "bSpline", "spline"))
	}

	periodicSpline.formula <- function(obj1, obj2, knots, period = 2 * pi, ord = 4L)
	{
	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 :
	structure(periodicSpline.default(as.numeric(eval(form[[3L]], local)),
	as.numeric(eval(form[[2L]], local)),
	knots = knots, period = period, ord = ord),
	formula = form)
	}

	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 - 1L))/
	prod(seq_len(i - 1L))
	}
	}
	structure(list(knots = knots, coefficients = coeff),
	formula = attr(object, "formula"),
	class = c("polySpline", "spline"))
	}

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

	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 - 1L
	coeff <- t(t(coeff[, -1]) * seq_len(ord))
	deriv <- deriv - 1L
	}
	y <- coeff[i, ord]
	if(ord > 1) {
	for(j in (ord - 1L):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 # (C_spline_value's FIXME)
	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") # predict.bSpline() -> NaN outside knots
	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)
	## Because of C_spline_value's FIXME, we move the outside-knots x[] values:
	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() # typically predict.polySpline()
	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)
	if((ord <- ncol(coeff)) != 4)
	stop("currently implemented only for cubic splines")
	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))