###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.005)
###x <- sort(runif(100,0,1)) ###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 

imageS <- function(lambda=1,...) {
  OmegaN <- U %*% diag(1/(1+lambda*d)) %*% t(U)
  K <- dim(OmegaN)[1]
  image(t(OmegaN[K-(1:K),]),...)
}

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)
  qplot(xy$x,xy$y,geom="line",ylab="Weight") + geom_point(aes(x=x,y=y),shape=20)
}


visU <- function(k=1,n=1000) {
  y <-  U[,k]
  xy <- spline(x,y,n=n)
  qplot(xy$x,xy$y,geom="line",ylab="Eigenfunction value") + geom_point(aes(x=x,y=y),shape=20)
}

##Examples, try playing with lambda and k

imageS(lambda=0.001,col=terrain.colors(40))
visS(k=100,lambda=0.001)
visU(k=20)