# splsmu.r          smoothing cubic splines
 #
  all <- read.table("eppright.dat")
  x <- all[,1]
  y <- all[,2]
 # direct calculation
  n <- length(x)
  nm1 <- n-1
  h <- diff(x)                # length n-1
  t <- 72*1.38                # smoothing parameter
  Estar <- diag(2*(h[-1]+h[-nm1]),n-2)            # diagonal
  diag(Estar[2:(n-2),1:(n-3)]) <- h[2:(n-2)]      # subdiagonal
  diag(Estar[1:(n-3),2:(n-2)]) <- h[2:(n-2)]      # superdiagonal
  Estar[1:6,1:6]              # write out part of this matrix
  C <- matrix(0,n-2,n)        
  diag(C[,1:(n-2)]) <- 1/h[-nm1]                  # left
  diag(C[,3:n]) <- 1/h[-1]                        # right
  diag(C[,2:nm1]) <- -(1/h[-1]+1/h[-nm1])         # middle or diagonal
  C[1:6,1:8]                  # let's see part of this
  Lt <- chol(Estar+6*t*C%*%t(C))                  # Cholesky factor
  Lt[1:6,1:6]                 # let's see part of this
  mstar <- backsolve(Lt,forwardsolve(t(Lt),6*C%*%y)) # slopes
  z <- y - t*t(C)%*%mstar                         # fitted values
 #
 # compare to native function      v -- note smoothing parameter
  smuspl <- smooth.spline(x,y,spar=0)
  smuspl$lambda
  smuspl$lev
  smuspl <- smooth.spline(x,y,spar=.654278)       # close to lam=1.38
  smuspl$lambda
  smuspl$lev
 # 
 # compare output
  cbind(x,y,smuspl$y,z,c(0,mstar,0))              
 # clean up & go
  rm(list=ls())
  q()