bsplinesetup.fun <- function(y, n=min(50, dim(y)[2]))
{
#
#   Given a matrix y of data, find a B-spline basis based on at most (n-2)  internal knots. 
#   In this version, the number of knots is reduced if there are not at least n time points 
#     observed data
#
#  On input:  The rows of y represent observations and the columns time points.
#
#  Output:ycoef	Matrix of spline coefficients, columns represent observations
#	bmat	Matrix giving values of bsplines at observation points
#	jmat	Matrix of integrals of crossproducts of B-splines
#	kmat	Matrix of integrals of crossproducts of second derivatives of B-splines
#	meancoef  Mean B-spline coefficients
#	funmean	Average function values
#
#     First find the knots
#
	if(!is.matrix(y))  y <- t(as.matrix(y))
	m <-  dim(y)[2]
	nd <- dim(y)[1]
	n <- min(n, m)
	intkn <- seq(from = 0, to = 1, length = n)
	kn <- c(0, 0, 0, intkn, 1, 1, 1)	
#
#   Now find the crossintegral matrix jmat 
#
	bmat1 <- spline.des(knots = kn, x = seq(from = 0, to = 1, 
		length = 101))[[4]]
	jmat <- 0.01 * t(bmat1) %*% bmat1	
#
#  find the matrix kmat of integrals of products of second derivatives
#
	zs <- spline.des(kn, intkn, deriv = rep(2,n))[[4]]
	z1 <- zs[-1,  ]
	z2 <- zs[-n,  ]
	kmat <- 2*t(z1)%*%z1 + 2*t(z2)%*%z2 + t(z1)%*%z2 + t(z2)%*%z1
	kmat <- kmat/(6*(n-1))	

#
#   Now find the matrix of values at the evaluation points, 
#     and the coefficients of the data.  
#   If there are more bsplines than evaluation points, 
#     find the fit with the least roughness.
#
	bmat <- spline.des(knots=kn, x=seq(from=0, to=1, length=m))[[4]]
	if((n + 2) <= m) {
		ycoef <- solve(bmat, t(y))
	}
	else {
		bsvd <- svd(bmat, nv = n + 2)
		baug <- rbind(bmat, t(bsvd$v[, (m + 1):(n + 2)]) %*% kmat)
		zmat <- matrix(0, nrow = n + 2 - m, ncol = nd)
		yaug <- rbind(t(y), zmat)
		ycoef <- solve(baug, yaug)
	}
#
#  Now find the mean coefficients and the values of the mean of the original data
#
	meancoef <- apply(ycoef, 1, mean)
	funmean <- bmat %*% meancoef
	return(ycoef, bmat, jmat, kmat, meancoef, funmean)
}
bvals.fun <- function(yc, nvals = 101)
{
#
#   Given a vector of B-spline coefficients yc, 
#
	kn <- c(0,0,0, seq(from=0, to=1, length=length(yc)-2), 1,1,1)
	xx <- seq(from = 0, to = 1, length = nx)
	bmat1 <- spline.des(knots = kn, x = xx)[[4]]
	yfun <- bmat1 %*% yc
	return(x = xx, y = yfun)
}
funsmooth.fun <-  function(zs, spar)
{
#  	carry out roughness penalty smoothing of the mean curve, following preprocessing
#   	by the function bsplinesetup.fun.          On input 
#	zs is the output of bsplinesetup.fun
#    	spar is the smoothing parameter
#	On output
#     the values g of the smoothed mean function at the evaluation points
#    	and the coefficients gcoef of the smoothed mean function are returned
#    	
#
	jlk <- zs$jmat + spar * zs$kmat
	gcoef <- solve(jlk, zs$jmat %*% zs$meancoef)
	g <- zs$bmat %*% gcoef
	return(g, gcoef)
}
funsmoothcv.fun <- function(zs, spar)
{
#  	calculates crossvalidation score for roughness penalty smoothing of the mean curve
# 	following preprocessing by the function bsplinesetup.fun.  On input 
#		zs is the output of bsplinesetup.fun
#    		spar is a vector of smoothing parameters
#  	The vector of cross-validation scores is returned.
#
	nsp <- length(spar)
	n <- dim(zs$ycoef)[2]
	cvscores <- rep(NA, nsp)
	for(j in (1:nsp)) {
		jlk <- zs$jmat + spar[j] * zs$kmat
		gcoef <- solve(jlk, zs$jmat %*% zs$meancoef)
		smat <- solve(jlk, zs$jmat %*% zs$ycoef)
		delresmat <- n * sweep(zs$ycoef, 1, gcoef) + smat -zs$ycoef
		delresmatj <- t(delresmat) %*% chol(zs$jmat)
		cvscores[j] <- sum(delresmatj^2)
	}
	return(cvscores)
}
funpca.fun <-  function(zz, spar = 1, zmean = zz$meancoef)
{
#
#  Takes the output of bsplinesetup.fun and carries out a smoothed functional PCA
#    with smoothing parameter equal to spar. This is taken straight from page 118
#    of Ramsay and Silverman, 1997, but with a vital correction.  
#       If centered=T then it is assumed that the 
#    mean or estimated mean has been subtracted from the coefficients zz$ycoef
#
	zm <- zz$jmat + spar * zz$kmat
	zm <- 0.5 * (zm + t(zm))
	umat <- chol(zm)
	lmat <- t(umat)
	ycmat <- zz$jmat %*% sweep(zz$ycoef, 1, zmean)
	dmat <- solve(lmat, ycmat)
	cmat <- dmat %*% t(dmat)/dim(ycmat)[2]
	eig <- eigen(cmat)
	emat <- eig$vectors[, 1:10]
	pcasmat <- solve(umat, emat)
	pcabmat <- zz$bmat %*% pcasmat
	vn <- apply(pcabmat, 2, vecnorm)
	pcabmat <- sweep(pcabmat, 2, vn, "/")
	return(pcabmat)
}