### Functions for computing the basis functions for
### natural cubic splines and the penalty matrix Omega_N.
### Also a function for computing the second derivatives.

natCubicSplines <- function(x, knots = x){
  K <- length(knots)
  evalList <- list()
  evalList[[1]] <- rep(1, length.out = length(x))
  evalList[[2]] <- x
  for(i in 1:(K-2)) {
    evalList[[i+2]] <- (pmax(x - knots[i], 0)^3 - pmax(x - knots[K], 0)^3)/(knots[K] - knots[i]) -
       (pmax(x - knots[K-1], 0)^3 - pmax(x - knots[K], 0)^3)/(knots[K] - knots[K-1])
  }
  result <- as.data.frame(evalList)
  names(result) <- x
  return(as.matrix(result))
}

OmegaNij <- function(i, j, x) {
  K <- length(x)
  mx <- pmax(x[i], x[j])
  mi <- pmin(x[i], x[j])
  ij <-   ((18*(x[K-1] - mx)*(x[K-1] - mi)^2 + 6*(mx - mi)^3 - 6*(x[K-1] - mi)^3) +
              12*(x[K-1] - x[i])*(x[K-1] - x[j])*(x[K] - x[K-1]))/((x[K] - x[i])*(x[K] - x[j]))
  return(ij)
}

natSplinesDeriv <- function(x, knots = x){
  K <- length(knots)
  evalList <- list()
  evalList[[1]] <- rep(0, length.out = length(x))
  evalList[[2]] <- rep(0, length.out = length(x))
  for(i in 1:(K-2)) {
    evalList[[i+2]] <- 6*((pmax(x - knots[i], 0) - pmax(x - knots[K], 0))/(knots[K] - knots[i]) -
       (pmax(x - knots[K-1], 0) - pmax(x - knots[K], 0))/(knots[K] - knots[K-1]))
  }
  result <- as.data.frame(evalList)
  names(result) <- knots
  return(as.matrix(result))
}

### Example on [1, 2] with equidistant knots

x <- seq(1, 2, 0.01)
### x <- sort(runif(100, 0, 1)) ## Try this for non-equidistant distribution of the x values.  
### x <- sort(rexp(100)) ## Try this for non-equidistant distribution of the x values.  
N <-  natCubicSplines(x)
K <- length(x)
OmegaN <- matrix(0, ncol = K, nrow = K)
OmegaN[3:K, 3:K] <- outer(1:(K-2), 1:(K-2), OmegaNij, x = x)
S <- N %*% solve(t(N) %*% N + OmegaN, t(N)) ## This is not numerically stable if there are too many knots!
d <- svd(S)
U <- d$u
d <- 1/d$d-1

### The following three functions can be used to visualize the smoother S
### The functions imageS and visS are comparable to the plots in Figure 5.8 in the book.
### The function visU shows the Demmler-Reinsch basis functions 
require(ggplot2)

imageS <- function(lambda = 1, cuts = 50, exc = 2) {
  S <- U %*% diag(1/(1 + lambda*d)) %*% t(U)
  K <- dim(S)[1]
  inc <- seq(exc + 1, K - exc - 1)
  image(Matrix(S)[inc, inc],
        cuts = cuts,
        col.regions = colorRampPalette(c("green", "red"))(cuts+10))
}

visS <- function(k = 1, lambda = 1, n = 1000, ...) {
  y <- as.vector(U[k, ] %*% diag(1/(1 + lambda*d)) %*% t(U))
  xy <- spline(x, y, n = n)
  plot(xy$x, xy$y, type = "l", ylab = "Weight", ...)
  points(x, y, pch = 20)
  lines(range(x), c(0, 0), lty = 2)
}

visU <- function(k = 1, n = 1000, ...) {
  y <-  U[, k]
  xy <- spline(x, y, n = n)
  plot(xy$x, xy$y, type = "l", ylab = "Eigenfunction value", ...)
  points(x, y, pch = 20)
  lines(range(x), c(0, 0), lty = 2)  
}

##Examples, try playing with lambda and k

imageS(lambda = 0.001, exc = 20)
visS(k = 4, lambda = 0.001)
visU(k = 20)

### A small "dynamic" visualization

lapply(seq_along(x),
       function(i) {
         visS(k = i, lambda = 0.001, ylim = c(-0.01, 0.1))
         Sys.sleep(0.2)
       }
       )