## B-spline basis functions.
## Single basis functions.
bsbasis <- function(z, knots, j, degree) {
  if(degree == 0)
    B <- 1 * (knots[j] <= z & z < knots[j + 1])
  if(degree > 0) {
    b1 <- (z - knots[j]) / (knots[j + degree] - knots[j])
    b2 <- (knots[j + degree + 1] - z) /
      (knots[j + degree + 1] - knots[j + 1])
    B <- b1 * bsbasis(z, knots, j, degree - 1) +
      b2 * bsbasis(z, knots, j + 1, degree - 1)
  }
  B[is.na(B)] <- 0
  return(B)
}

## And the complete design matrix.
bs <- function(z, degree = 3, knots = NULL) {
  ## Compute knots.
  if(is.null(knots))
    knots <- 40
  if(length(knots) < 2) {
    step <- (max(z) - min(z)) / (knots - 1)
    knots <- seq(min(z) - degree * step,
      max(z) + degree * step, by = step)
  }

  ## Evaluate each basis function
  ## and return the full design matrix B.
  B <- NULL
  for(j in 1:(length(knots) - degree - 1))
    B <- cbind(B, bsbasis(z, knots, j, degree))
  return(B)
}

## Penalty matrix based on differences.
P <- function(order = 2, k = 7) {
  D <- diag(k)
  for(i in 1:order)
    D <- diff(D)
  K <- crossprod(D, D)
  return(K)
}