R/spbase.R

In JOPS: Practical Smoothing with P-Splines

#' Compute a sparse B-spline basis on evenly spaced knots
#'
#' @description Constructs a sparse B-spline basis on evenly spaced knots.
#'
#' @import spam
#'
#' @param x a vector of argument values, at which the B-spline basis functions
#' are to be evaluated.
#' @param xl the lower limit of the domain of \code{x} (default \code{min(x)}).
#' @param xr the upper limit of the domain of \code{x} (default \code{max(x)}) .
#' @param nseg the number of evenly spaced segments between \code{xl} and \code{xr} (default 10).
#' @param bdeg the degree of the basis, usually 1, 2, or 3 (default).
#' @return A sparse matrix (in \code{spam} format) with \code{length(x)}  of rows= and \code{nseg + bdeg} columns.
#'
#' @author Paul Eilers
#'
#' @references Eilers, P.H.C. and Marx, B.D. (1996). Flexible smoothing with
#' B-splines and penalties (with comments and rejoinder), \emph{Statistical Science}, 11: 89-121.
#' @references Eilers, P.H.C. and Marx, B.D. (2021). \emph{Practical Smoothing, The Joys of
#' P-splines.} Cambridge University Press.
#' @examples
#' library(JOPS)
#' # Basis  on grid
#' x = seq(0, 4, length = 1000)
#' B = spbase(x, 0, 4, nseg = 50, bdeg = 3)
#' nb1 = ncol(B)
#' matplot(x, B, type = 'l', lty = 1, lwd = 1, xlab = 'x', ylab = '')
#' cat('Dimensions of B:', nrow(B), 'by', ncol(B), 'with', length(B@entries), 'non-zero elements' )
#'
#' @export

spbase = function(x, xl = min(x), xr = max(x), nseg = 10, bdeg = 3) {

  # Check domain and adjust it if necessary
  if (xl > min(x)) {
    xl = min(x)
    warning("Left boundary adjusted to min(x) = ", xl)
  }
  if (xr <  max(x)) {
    xr = max(x)
    warning("Right boundary adjusted to max(x) = ", xr)
  }

  # Reduce x to first interval between knots
  m = length(x)
  dx = (xr - xl) / nseg
  ix <- floor((x - xl) / (1.0000001 * dx))
  xx = (x - xl) - ix * dx

  # Full basis for reduced x
  Br = bbase(xx, xl = 0, xr = dx, nseg = 1, bdeg = bdeg)

  # Compute proper rows, columns
  nr = ncol(Br)
  rw = rep(1:m, each = nr)
  cl = rep(1:nr, m) + rep(ix, each = nr)

  # Make the sparse matrix
  b = as.vector(t(Br))
  Bs = spam(list(i = rw, j = cl, b), nrow = m, ncol = (nseg + bdeg))
  att1 = attributes(Bs)
  att2 = list(x = x, xl = xl, xr = xr, nseg = nseg, bdeg = bdeg, type = 'spbase')
  attributes(Bs) <- append(att1, att2)
  return(Bs)
}