# splsmu.s                         J F Monahan, June 2001
 # smoothing cubic splines demonstration using Eppright problem
 xy <- scan('eppright.dat')
 n <- length(xy)/2
 n
 xy <- matrix(xy,2,n)
 x <- xy[1,]
 y <- xy[2,]
# natural spline 
 splsmu <- function(xk,y,n,t) {
   h <- c( 0, (xk[2:n] - xk[1:(n-1)] ) )
   z <- c( 0, 1/h[2:n] )                      # temporarily
   M <- (y[2:n] - y[1:(n-1)] ) / h[2:n]
   M <- 6*( M[2:(n-1)] - M[1:(n-2)] )
   Dd <- 2*(h[2:(n-1)]+h[3:n]) +
        6*t*( z[2:(n-1)]*z[2:(n-1)] + z[3:n]*z[3:n] +
                 (z[2:(n-1)]+z[3:n])*(z[2:(n-1)]+z[3:n]) )
   D1 <- h[3:(n-1)] - 6*t*( (z[2:(n-2)] + z[3:(n-1)])*z[3:(n-1)]
                          + (z[3:(n-1)] + z[4:n])*z[3:(n-1)]  )
   D2 <- 6*t*z[3:(n-2)]*z[4:(n-1)]
 # Cholesky factorization of (D2, D1, Dd, D1, D2)
  # first few and start solving
   Dd[1] <- sqrt( Dd[1] )
   M[1] <- M[1]/Dd[1]
   D1[1] <- D1[1] / Dd[1]
   Dd[2] <- sqrt( Dd[2] - D1[1]*D1[1] )
   M[2] <- ( M[2] - D1[1]*M[1] ) / Dd[2]
   for( j in 3:(n-2) ) {
   D2[j-2] <- D2[j-2] / Dd[j-2]
   D1[j-1] <- ( D1[j-1] - D2[j-2]*D1[j-2] ) / Dd[j-1]
   Dd[j] <- sqrt( Dd[j] - D1[j-1]*D1[j-1] - D2[j-2]*D2[j-2] )
   M[j] <- ( M[j] - D2[j-2]*M[j-2] - D1[j-1]*M[j-1] )/Dd[j]
                   }   # end of loop on j
 # solve system backwards
   M[n-2] <- M[n-2] / Dd[n-2]
   M[n-3] <- (M[n-3] - D1[n-3]*M[n-2])/Dd[n-3]
   for( j in (n-4):1 ) M[j]<-(M[j]-D1[j]*M[j+1]-D2[j]*M[j+2])/Dd[j]
 # natural spline
   M <- c( 0, M, 0 )
 # fitted values
   z[1] <- y[1] - t*z[2]*M[2]
   z[2] <- y[2] + t*((z[2]+z[3])*M[2] - z[3]*M[3] )
   for( j in 3:(n-2) ) {
    z[j] <- y[j] - t*(z[j]*M[j-1]-(z[j]+z[j+1])*M[j]+z[j+1]*M[j+1])
                       }
   z[n-1] <- y[n-1] - t*( z[n-1]*M[n-2] - (z[n-1]+z[n])*M[n-1] )
   z[n] <- y[n] - t*z[n]*M[n-1]
   return( list(xk=xk,z=z,n=n,M=M,h=h) )
                            } # end of function splsmu
# evaluate spline at x vector
  splev <- function(x,splinfo) {
     j <- 1 + as.numeric( cut(x,splinfo$xk) )
     du <- splinfo$xk[j] - x
     dl <- x - splinfo$xk[j-1]
     Mu <- splinfo$M[j]
     Ml <- splinfo$M[j-1]
     hj <- splinfo$h[j]   
     ( (Ml*du*du*du + Mu*dl*dl*dl)/6
                   + (splinfo$z[j-1]-Ml*hj*hj/6 )*du
                   + (splinfo$z[j]-Mu*hj*hj/6 )*dl )/hj
                                }  # end of function splev 
   splepp <- splsmu(x,y,n,n*1.38)
   splepp$M
   splepp$n
 fspl <- splev(x,splepp)

 out <- paste( 'x', format(x),'y',format(y),
               'fit',format(fspl),'der',format(splepp$M) )
 cat(out,sep="\n")
 rm(list=ls())