library.dynam("splines.so")

"bs"<- function(x, df = NULL, knots = NULL, degree = 3, intercept = F, Boundary.knots
	 = range(x))
{
	nx <- names(x)
	x <- as.vector(x)
	nax <- is.na(x)
	if(nas <- any(nax))
		x <- x[!nax]
	if(!missing(Boundary.knots)) {
		Boundary.knots <- sort(Boundary.knots)
		outside <- (ol <- x < Boundary.knots[1]) | (or <- x > 
			Boundary.knots[2])
	}
	else outside <- rep(F, length = length(x))
	if(!missing(df) && missing(knots)) {
		nIknots <- df - (degree + 1) + (1 - intercept)
		if(nIknots < 0) {
			nIknots <- 0
			warning(paste("df was too small; have used ", (degree + 
				1) - (1 - intercept)))
		}
		if(nIknots > 0) {
			knots <- seq(from = 0, to = 1, length = nIknots + 2)[ - 
				c(1, nIknots + 2)]
			knots <- quantile(x[!outside], knots)
		}
		else knots <- NULL
	}
	Aknots <- sort(c(rep(Boundary.knots, degree + 1), knots))
	if(any(outside)) {
		warning("Some x values beyond boundary knots may cause ill-conditioned bases"
			)
		derivs <- seq(from = 0, to = degree)	
	# tricky way to do factorial
		scalef <- c(1, exp(cumsum(log(seq(degree)))))
		basis <- array(0, c(length(x), length(Aknots) - degree - 1))
		if(any(ol)) {
			k.pivot <- Boundary.knots[1]
			xl <- cbind(1, outer(x[ol] - k.pivot, 1:degree, "^"))
			tt <- spline.des(Aknots, rep(k.pivot, degree + 1), 
				degree + 1, derivs)$design
			basis[ol,  ] <- xl %*% (tt/scalef)
		}
		if(any(or)) {
			k.pivot <- Boundary.knots[2]
			xr <- cbind(1, outer(x[or] - k.pivot, 1:degree, "^"))
			tt <- spline.des(Aknots, rep(k.pivot, degree + 1), 
				degree + 1, derivs)$design
			basis[or,  ] <- xr %*% (tt/scalef)
		}
		if(any(inside <- !outside))
			basis[inside,  ] <- spline.des(Aknots, x[inside], 
				degree + 1)$design
	}
	else basis <- spline.des(Aknots, x, degree + 1)$design
	if(!intercept)
		basis <- basis[, -1]
	n.col <- ncol(basis)
	if(nas) {
		nmat <- matrix(NA, length(nax), n.col)
		nmat[!nax,  ] <- basis
		basis <- nmat
	}
	dimnames(basis) <- list(nx, 1:n.col)
	a <- list(degree = degree, knots = knots, Boundary.knots = 
		Boundary.knots, intercept = intercept, class = c("bs", "basis")
		)
	attributes(basis) <- c(attributes(basis), a)
	basis
}
"ns"<-function(x, df = NULL, knots = NULL, intercept = F, Boundary.knots = range(x))
{
	nx <- names(x)
	x <- as.vector(x)
	nax <- is.na(x)
	if(nas <- any(nax))
		x <- x[!nax]
	if(!missing(Boundary.knots)) {
		Boundary.knots <- sort(Boundary.knots)
		outside <- (ol <- x < Boundary.knots[1]) | (or <- x > 
			Boundary.knots[2])
	}
	else outside <- rep(F, length = length(x))
	if(!missing(df) && missing(knots)) {
# df = number(interior knots) + 1 + intercept
		nIknots <- df - 1 - intercept
		if(nIknots < 0) {
			nIknots <- 0
			warning(paste("df was too small; have used ", 1 + 
				intercept))
		}
		if(nIknots > 0) {
			knots <- seq(from = 0, to = 1, length = nIknots + 2)[ - 
				c(1, nIknots + 2)]
			knots <- quantile(x[!outside], knots)
		}
		else knots <- NULL
	}
	Aknots <- sort(c(rep(Boundary.knots, 4), knots))
	if(any(outside)) {
		basis <- array(0, c(length(x), nIknots + 4))
		if(any(ol)) {
			k.pivot <- Boundary.knots[1]
			xl <- cbind(1, x[ol] - k.pivot)
			tt <- spline.des(Aknots, rep(k.pivot, 2), 4, c(0, 1))$
				design
			basis[ol,  ] <- xl %*% tt
		}
		if(any(or)) {
			k.pivot <- Boundary.knots[2]
			xr <- cbind(1, x[or] - k.pivot)
			tt <- spline.des(Aknots, rep(k.pivot, 2), 4, c(0, 1))$
				design
			basis[or,  ] <- xr %*% tt
		}
		if(any(inside <- !outside))
			basis[inside,  ] <- spline.des(Aknots, x[inside], 4)$
				design
	}
	else basis <- spline.des(Aknots, x, 4)$design
	const <- spline.des(Aknots, Boundary.knots, 4, c(2, 2))$design
	if(!intercept) {
		const <- const[, -1]
		basis <- basis[, -1]
	}
	qr.const <- qr(t(const))
	basis <- as.matrix((t(qr.qty(qr.const, t(basis))))[,  - (1:2)])
	n.col <- ncol(basis)
	if(nas) {
		nmat <- matrix(NA, length(nax), n.col)
		nmat[!nax,  ] <- basis
		basis <- nmat
	}
	dimnames(basis) <- list(nx, 1:n.col)
	a <- list(degree = 4, knots = knots, Boundary.knots = Boundary.knots, 
		intercept = intercept, class = c("ns", "basis"))
	attributes(basis) <- c(attributes(basis), a)
	basis
}
"predict.bs"<-function(object, newx, ...)
{
	if(missing(newx))
		return(object)
	a <- c(list(x = newx), attributes(object)[c("degree", "knots", 
		"Boundary.knots", "intercept")])
	do.call("bs", a)
}
"predict.ns"<-function(object, newx, ...)
{
	if(missing(newx))
		return(object)
	a <- c(list(x = newx), attributes(object)[c("degree", "knots", 
		"Boundary.knots", "intercept")])
	do.call("ns", a)
}
"interp.spline" <- function(x, y, ord = 4)
{
	x <- as.numeric(x)
	lenx <- length(x)
	y <- as.numeric(y)
	leny <- length(y)
	if(leny!=lenx)
		stop("Lengths of x and y must be equal")
	ind <- order(x)
	x <- x[ind]
	y <- y[ind]
	ordm1 <- ord - 1
	knots <- c(x[1:ordm1] + x[1] - x[ord], x, x[lenx + (1:ordm1) - ordm1] +
		x[lenx] - x[lenx - ordm1])
	derivs <- c(2, integer(lenx), 2)
	x <- c(x[1], x, x[lenx])
	sys.mat <- spline.des(knots, x, ord = ord, derivs = derivs)$design
	coeff <- solve(sys.mat, c(0, y, 0))
	list(knots = knots, coeff = coeff, order = ord)
}
"periodic.spline" <- function(x, y, knots, period = 2 * pi, ord = 4)
{
	x <- as.vector(x, "double")
	y <- as.vector(y, "double")
	lenx <- length(x)
	if(lenx!=length(y))
		stop("Lengths of x and y must match")
	ind <- order(x)
	x <- x[ind]
	y <- y[ind]
	if(any((x[-1] - x[ - lenx]) <= 0))
		stop("Values of x must be strictly increasing")
	if(!missing(knots))
		period <- knots[length(knots) + 1 - ord] - knots[1]
	else knots <- c(x[(lenx - (ord - 2)):lenx] - period, x, x[1:ord] + 
			period)
	if((x[lenx] - x[1]) >= period)
		stop("The range of x values exceeds one period")
	coeff.mat <- spline.des(knots, x, ord)$design
	sys.mat <- coeff.mat[, (1:lenx)]
	sys.mat[, 1:(ord - 1)] <- sys.mat[, 1:(ord - 1)] + coeff.mat[, lenx +
		(1:(ord - 1))]
	qrstr <- qr(sys.mat)
	coeff <- qr.coef(qrstr, y)
	coeff <- c(coeff, coeff[1:(ord - 1)])
	list(knots = knots, coeff = coeff, order = ord, period = period, m = 
		coeff.mat)
}
"spline.value" <- function(bspstr, x = seq(knots[ord], knots[ncoeff + 1], len = lenx), deriv = 0,
	lenx = length(x))
{
##	library.dynam("splines", "splines.o")
	if(is.null(bspstr$knots))
		stop("bspstr must be a B-spline structure")
	if(missing(x) && missing(lenx))
		stop("either x or lenx must be given to spline.value")
	knots <- bspstr$knots
	ord <- bspstr$order
	ncoeff <- length(bspstr$coeff)
	x.orig <- x
	if(!is.null(period <- bspstr$period)) {
		ind <- x < knots[ord]
		if(any(ind))
			x[ind] <- x[ind] + period * (1 + (knots[ord] - x[ind]) %/%
				period)
		ind <- x > knots[ncoeff + 1]
		if(any(ind))
			x[ind] <- x[ind] - period * (1 + (x[ind] - knots[ncoeff +
				1]) %/% period)
	}
	xind <- order(x)
	x <- x[xind]
	rind <- order(xind)
	z <- .C("spline_value",
		as.double(bspstr$knots),
		as.double(bspstr$coeff),
		as.integer(length(bspstr$coeff)),
		as.integer(bspstr$order),
		as.double(x),
		as.integer(lenx),
		as.integer(deriv),
		y = double(lenx))
	list(x = x.orig, y = z$y[rind])
}
"spline.des" <- function(knots, x, ord = 4, derivs = integer(nx))
{
##	library.dynam("splines", "splines.o")
	knots <- sort(as.vector(knots))
	x <- as.vector(x)
	nk <- length(knots)
	nx <- length(x)
	ind <- order(x)
	sortx <- x[ind]
	ind <- order(ind)
	if(sortx[1] < knots[ord] || sortx[nx] > knots[nk + 1 - ord])
		stop(paste("The x data must be in the range", knots[ord], "to",
			knots[nk + 1 - ord]))
	if(length(derivs)!=nx)
		stop("length of derivs must match length of x")
	ncoef <- nk - ord
	temp <- .C("spline_basis",
		as.double(knots),
		as.integer(ncoef),
		as.integer(ord),
		as.double(sortx),
		as.integer(derivs),
		as.integer(nx),
		design = array(0, c(ord, nx)),
		offsets = integer(nx))
	design <- array(0, c(nx, ncoef))
	d.ind <- array(c(rep(1:nx, rep(ord, nx)), outer(1:ord, temp$offsets,
		"+")), c(nx * ord, 2))
	design[d.ind] <- temp$design
	list(knots = knots, order = ord, derivs = derivs, design = design[ind,  ])
}
"linear.interp" <- function(x, y, x0 = seq(x[1], x[lenx], len = lenx0), lenx0 = length(x0))
{
##	library.dynam("splines","splines.o")
	x <- as.numeric(x)
	y <- as.numeric(y)
	lenx <- length(x)
	if(length(y)!=lenx)
		stop("lengths of x and y must be the same")
	xind <- order(x)
	x <- x[xind]
	y <- y[xind]
	if(missing(x0) && missing(lenx0))
		stop("either x0 or lenx0 must be given")
	x0 <- as.numeric(x0)
	x0ind <- order(x0)
	x0ord <- x0[x0ind]
	nvals <- length(x0)
	if(x0ord[1] < x[1] || x0ord[nvals] > x[lenx])
		stop(paste("x0 values must be in the range", x[1], "to", x[
			lenx]))
	y <- .C("lin_interp",
		x,
		y,
		x0,
		y = double(nvals),
		as.integer(nvals))$y
	list(x = x0, y = y[order(x0ind)])
}
