The R Project SVN
R-packages
(root)/trunk/mgcv/R/sparse.r – Rev 8207
Rev 7433 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
## (c) Simon N. Wood 2011-2013## functions for sparse smoothing.tri2nei <- function(T) {## Each row of matrix T gives the indices of the points making up## one triangle in d dimensions. T has d+1 columns. This routine## finds the neighbours of each point in that triangulation.## indices start at 1.n <- max(T)oo <- .C(C_tri2nei,T=as.integer(cbind(T-1,T*0)),as.integer(nrow(T)),as.integer(n),as.integer(ncol(T)-1),off=as.integer(rep(0,n)));## ni[1:off[1]] gives neighbours of point 1.## ni[(off[i-1]+1):off[i]] give neighbours of point i>1return(list(ni = oo$T[1:oo$off[n]]+1,off=oo$off))}tri.pen <- function(X,T) {## finds a sparse approximate TPS penalty, based on the points in X,## with triangulation T. Rows of X are points. Rows of T index vertices## of triangles in X.nn <- tri2nei(T) ## get neighbour list from T## now obtain generalized FD penalty...n <- nrow(X);d <- ncol(X);D <- rep(0,3*(nn$off[n]+n)) ## storage foroo <- .C(C_nei_penalty,as.double(X),as.integer(n),as.integer(d),D=as.double(D),ni=as.integer(nn$ni-1),ii=as.integer(nn$ni*0),off=as.integer(nn$off),as.integer(2),as.integer(0),kappa=as.double(rep(0,n)));## unpack into sparse matrices...ni <- oo$off[n]ii <- c(1:n,oo$ii[1:ni]+1) ## row indexjj <- c(1:n,oo$ni[1:ni]+1) ## col indexni <- length(ii)Kx <- sparseMatrix(i=ii,j=jj,x=oo$D[1:ni],dims=c(n,n))Kz <- sparseMatrix(i=ii,j=jj,x=oo$D[1:ni+ni],dims=c(n,n))Kxz <- sparseMatrix(i=ii,j=jj,x=oo$D[1:ni+2*ni],dims=c(n,n))list(Kx=Kx,Kz=Kz,Kxz=Kxz)}## Efficient stable full rank cubic spline routines, based on## deHoog and Hutchinson, 1987 and Hutchinson and deHoog,## 1985....setup.spline <- function(x,w=rep(1,length(x)),lambda=1,tol=1e-9) {## setup a cubic smoothing spline given data locations in x.## ties will be treated by removing duplicate x's, and averaging corresponding## y's. Averaging is \sum_i y_i w_i^2 / \sum_i w_i^2, and the weight## assigned to this average is then w_a^2 = \sum_i w_i^2...## spline object has to record duplication information, as well as## rotations defining spline.n <- length(x)ind <- order(x)x <- x[ind] ## sort xw <- w[ind]U <- V <- rep(0,4*n)diagA <- rep(0,n)lb <- rep(0,2*n)oo <- .C(C_sspl_construct,as.double(lambda),x=as.double(x),w=as.double(w),U=as.double(U),V=as.double(V),diagA=as.double(diagA),lb=as.double(lb),n=as.integer(n),tol=as.double(tol))list(trA=sum(oo$diagA), ## trace of influence matrixU=oo$U,V=oo$V, ## spline defining Givens rotationslb=oo$lb, ## final lower bandx=x, ## original x sequence, orderedind=ind, ## x0 <- x; x0[ind] <- x, puts original ordering in x0w=w, ## original weightsns=oo$n, ## number of unique x values (maximum spline rank)tol=tol) ## tolerance used to judge tied x values}apply.spline <- function(spl,y) {## Use cubic spline object spl, from setup.spline, to smooth data in y.if (is.matrix(y)) {m <- ncol(y)y <- y[spl$ind,] ## order as x} else {m <- 1y <- y[spl$ind]}n <- length(spl$x)oo <- .C(C_sspl_mapply,f = as.double(y),x=as.double(spl$x),as.double(spl$w),U=as.double(spl$U),as.double(spl$V),n=as.integer(spl$ns),nf=as.integer(n),tol=as.double(spl$tol),m=as.integer(m))if (is.matrix(y)) {y <- matrix(oo$f,n,m)y[spl$ind,] <- y ## original order} else {y[spl$ind] <- oo$f}y}## kd tree/k nearest neighbout related routines....kd.vis <- function(kd,X,cex=.5) {## code visualizes a kd tree for points in rows of X## kd <- kd.tree(X) produces correct tree.## this worked with the structures used when## kd tree was written out to R vecotrs and the read## back in. Does not work with revised approach (would## need an explicit helper function to do the write out)if (ncol(X)!=2) stop("only deals with 2D case")##n <- nrow(X)d <- ncol(X)nb <- kd$idat[1]dd <- matrix(kd$ddat[-1],nb,2*d,byrow=TRUE)lo <- dd[,1:d];hi <- dd[,1:d+d]rm(dd)ll <- min(X[,1]); ul<- max(X[,1])w <- ul-llind <- lo[,1] < ll-w/10;lo[ind,1] <- ll-w/10ind <- hi[,1] > ul+w/10;hi[ind,1] <- ul+w/10ll <- min(X[,2]);ul <- max(X[,2])w <- ul-llind <- lo[,2] < ll-w/10;lo[ind,2] <- ll-w/10ind <- hi[,2] > ul+w/10;hi[ind,2] <- ul+w/10plot(X[,1],X[,2],pch=19,cex=cex,col=2)for (i in 1:nb) {rect(lo[i,1],lo[i,2],hi[i,1],hi[i,2])}#points(X[,1],X[,2],pch=19,cex=cex,col=2)}nearest <- function(k,X,gt.zero = FALSE,get.a=FALSE) {## The rows of X contain coordinates of points.## For each point, this routine finds its k nearest## neighbours, returning a list of 2, n by k matrices:## ni - ith row indexes the rows of X containing## the k nearest neighbours of X[i,]## dist - ith row is the distances to the k nearest## neighbours.## a - area associated with each point, if get.a is TRUE## ties are broken arbitrarily.## gt.zero indicates that neighbours must have distances greater## than zero...if (gt.zero) {Xu <- uniquecombs(X);ind <- attr(Xu,"index") ## Xu[ind,] == X} else { Xu <- X; ind <- 1:nrow(X)}if (k>nrow(Xu)) stop("not enough unique values to find k nearest")nobs <- length(ind)n <- nrow(Xu)d <- ncol(Xu)dist <- matrix(0,n,k)if (get.a) a <- 1:n else a=1oo <- .C(C_k_nn,Xu=as.double(Xu),dist=as.double(dist),a=as.double(a),ni=as.integer(dist),n=as.integer(n),d=as.integer(d),k=as.integer(k),get.a=as.integer(get.a))dist <- matrix(oo$dist,n,k)[ind,]rind <- 1:nobsrind[ind] <- 1:nobsni <- matrix(rind[oo$ni+1],n,k)[ind,]if (get.a) a=oo$a[ind] else a <- NULLlist(ni=ni,dist=dist,a=a)} # nearest#kd.tree <- function(X) {## function to obtain kd tree for points in rows of X## old version based on writing out tree to R# n <- nrow(X) ## number of points# d <- ncol(X) ## dimension of points# ## compute the number of boxes in the kd tree, nb# m <- 2;while (m < n) m <- m* 2;# nb = n * 2 - m %/% 2 - 1;# if (nb > m-1) nb = m - 1;# ## compute the storage requirements for the tree# nd = 1 + d * nb * 2 ## number of doubles# ni = 3 + 5 * nb + 2*n ## number of integers# oo <- .C(C_Rkdtree,as.double(X),as.integer(n),as.integer(d),idat = as.integer(rep(0,ni)),# ddat = as.double(rep(0,nd)))# list(idat=oo$idat,ddat=oo$ddat)#}kd.tree <- function(X) {## Function to obtain kd tree for points in rows of X.## Contains tree dumped as a vector of doubles with an attribute that is a vector## of integers. Another attribute is the internal pointer to the tree.## Redundant structure allows tree to be saved to disk and re-loaded. Pointer is## set to NULL by this, but can tree can then be restored by kd.nearest and kd.radius## and pointer reset. This works because documented behaviour is never to copy## such external pointers within R - they are effectively global - so resetting resets## it for every copy of the tree.kd <- .Call(C_Rkdtree,X)kd}kd.nearest <- function(kd,X,x,k) {## Finds k nearest neigbours of each row of x within X. X has## corresponding kd tree kd, stored as an external pointer:## attribute "kd_ptr" of kd. Returns array of indices to rows in## X, along with corresponding distances.nei <- .Call(C_Rkdnearest,kd,X,x,as.integer(k)) + 1 ## C to R}kd.radius <- function(kd,X,x,r){## find all points in kd tree (kd,X) in radius r of points in x.## kd should come from kd.tree(X).## neighbours of x[i,] in X are the rows given by ni[off[i]:(off[i+1]-1)]# m <- nrow(x);# off <- rep(0,m+1)off <- rep(as.integer(0),nrow(x)+1)xt <- t(x) ## required transposed in Rkradiusni <- .Call(C_Rkradius,kd,X,xt,as.double(r),off) + 1list(ni=ni,off=off+1)}#kd.nearest <- function(kd,X,x,k) {## given a set of points in rows of X, and corresponding kd tree, kd## (produced by a call to kd.tree(X)), then this routine finds the## k nearest neighbours in X, to the points in the rows of x.## outputs: ni[i,] lists k nearest neighbours of X[i,].## dost[i,] is distance to those neighbours.## note R indexing of output# n <- nrow(X)# m <- nrow(x)# ni <- matrix(0,m,k)# oo <- .C(C_Rkdnearest,as.double(X),as.integer(kd$idat),as.double(kd$ddat),as.integer(n),as.double(x),# as.integer(m), ni=as.integer(ni), dist=as.double(ni),as.integer(k))# list(ni=matrix(oo$ni+1,m,k),dist=matrix(oo$dist,m,k))#}#kd.radius <- function(kd,X,x,r) {## find all points in kd tree (kd,X) in radius r of points in x.## kd should come from kd.tree(X).## neighbours of x[i,] in X are the rows given by ni[off[i]:(off[i+1]-1)]# m <- nrow(x);# off <- rep(0,m+1)## do the work...# oo <- .C(C_Rkradius,as.double(r),as.integer(kd$idat),as.double(kd$ddat),as.double(X),as.double(t(x)),# as.integer(m),off=as.integer(off),ni=as.integer(0),op=as.integer(0))# off <- oo$off# ni <- rep(0,off[m+1])## extract to R and clean up...# oo <- .C(C_Rkradius,as.double(r),as.integer(kd$idat),as.double(kd$ddat),as.double(X),as.double(t(x)),# as.integer(m),off=as.integer(off),ni=as.integer(ni),op=as.integer(1))# list(off=off+1,ni=oo$ni+1) ## note R indexing here.#} ## kd.radiustieMatrix <- function(x) {## takes matrix x, and produces sparse matrix P that maps list of unique## rows to full set. Matrix of unique rows is returned in xu.## If a smoothing penalty matrix, S, is set up based on rows of xu,## then P%*%solve(t(P)%*%P + S,t(P)) is hat matrix.x <- as.matrix(x)n <- nrow(x)xu <- uniquecombs(x)if (nrow(xu)==nrow(x)) return(NULL)ind <- attr(xu,"index")x <- as.matrix(x)n <- nrow(x)P <- sparseMatrix(i=1:n,j=ind,x=rep(1,n),dims=c(n,nrow(xu)))return(list(P=P,xu=xu))}## sparse smooths must be initialized with...## 1. a set of variable names, a blocking factor and a type.########################################################## routines for full rank cubic spline smoothers, based on# deHoog and Hutchinson, 1987.#########################################################spasm.construct.cus <- function(object,data) {## entry object inherits from "cus" & contains:## * terms, the name of the argument of the smooth## * block, the name of a blocking factor. Can be NULL.## return object also has...## * nobs - number of observations in total## * nblock - number of blocks.## * ind, list where ind[[i]] gives rows to which block i applies.## * spl, and empty list which will contain intialised cubic## spline smoothers for each block, once a smoothing parameter## has been supplied...##dat <- list()d <- length(object$terms)if (d != 1) stop("cubic spline only deals with 1D data")object$x <- get.var(object$term[1],data)object$nobs <- length(object$x)ind <- list()n <- length(object$x)## if there is a blocking factor then set up indexes## indicating which data go with which block...if (!is.null(object$block)) {block <- as.factor(get.var(object$block,data))nb <- length(levels(block))edf1 <- 0for (i in 1:nb) {ind[[i]] <- (1:n)[block==levels(block)[i]]edf1 <- edf1 + length(unique(object$x[ind[[i]]])) ## max edf for this block}} else { ## all one blocknb <- 1ind[[1]] <- 1:nedf1 <- length(unique(object$x))}object$nblock <- nbobject$ind <- ind## so ind[[i]] indexes the elements operated on by the ith smoother.object$spl <- list()object$edf0 <- 2*nb;object$edf1 <- edf1class(object) <- "cus"object}spasm.sp.cus <- function(object,sp,w=rep(1,object$nobs),get.trH=FALSE,block=0,centre=FALSE) {## Set up full cubic spline smooth, given new smoothing parameter and weights.## In particular, construct the Givens rotations defining the## smooth and compute the trace of the influence matrix.## If block is non-zero, then it specifies which block to set up, otherwise## all are set up. In either case w is assumed to be for the whole smoother,## although only the values for the specified block(s) are used.## If centre == TRUE then the spline is set up for centred smoothing, i.e.## the results sum to zero.## Note: w propto 1/std.dev(response)if (is.null(object$spl)) stop("object not fully initialized")trH <- 0if (block==0) block <- 1:object$nblockfor (i in block) {##n <- length(object$ind[[i]])object$spl[[i]] <- setup.spline(object$x[object$ind[[i]]],w=w[object$ind[[i]]],lambda=sp)trH <- trH + object$spl[[i]]$trA}object$sp=spif (get.trH) {if (centre) { ## require correction for DoF lost by centring...for (i in block) {one <- rep(1,length(object$ind[[i]]))object$centre <- FALSEtrH <- trH - mean(spasm.smooth(object,one,block=i))}}object$trH <- trH}object$centre <- centreobject}spasm.smooth.cus <- function(object,X,residual=FALSE,block=0) {## apply smooth, or its residual operation to X.## if block == 0 then apply whole thing, otherwise X must have the correct## number of rows for the smooth block.if (block>0) {## n <- length(object$ind[[block]])if (object$centre) {X0 <- apply.spline(object$spl[[block]],X)if (is.matrix(X0)) {x0 <- colMeans(X0)X0 <- sweep(X0,2,x0)} else X0 <- X0 - mean(X0)if (residual) X <- X - X0 else X <- X0} else {if (residual) X <- X - apply.spline(object$spl[[block]],X)else X <- apply.spline(object$spl[[block]],X)}} else for (i in 1:object$nblock) { ## work through all blocksind <- object$ind[[i]]if (is.matrix(X)) {X0 <- apply.spline(object$spl[[i]],X[ind,])if (object$centre) X0 <- sweep(X0,2,colMeans(X0))if (residual) X[ind,] <- X[ind,] - X0else X[ind,] <- X0} else {X0 <- apply.spline(object$spl[[i]],X[ind])if (object$centre) X0 <- X0 - mean(X0)if (residual) X[ind] <- X[ind] - X0else X[ind] <- X0}}X}########################################################### The default sparse smooth class, which does nothing...#########################################################spasm.construct.default <- function(object,data) {## This smooth simply returns 0, under all circumstances.## object might contain....## * block, the name of a blocking factor. Can be NULL.## return object also has...## * nblock - number of blocks.## * ind, list where ind[[i]] gives rows to which block i applies.n <- nrow(data)if (!is.null(object$block)) {block <- as.factor(get.var(object$block,data))nb <- length(levels(block))for (i in 1:nb) {ind[[i]] <- (1:n)[block==levels(block)[i]]}} else { ## all one blocknb <- 1ind[[1]] <- 1:n}object$nblock <- nbobject$ind <- ind## so ind[[i]] indexes the elements operated on by the ith smoother.class(object) <- "default"object}spasm.sp.default <- function(object,sp,w=0,get.trH=FALSE,block=0,centre=FALSE) {## Set up default null smoother. i.e. set trH=0,trH <- 0if (get.trH) object$trH <- trHobject$ldetH <- NAobject}spasm.smooth.default <- function(object,X,residual=FALSE,block=0) {## apply smooth, or its residual operation to X.## if block == 0 then apply whole thing, otherwise X must have the correct## number of rows for the smooth block.if (residual) return(X) else return(X*0)X}## generics for sparse smooth classes...spasm.construct <- function(object,data) UseMethod("spasm.construct")spasm.sp <- function(object,sp,w=rep(1,object$nobs),get.trH=TRUE,block=0,centre=FALSE) UseMethod("spasm.sp")## Note that w is assumed proportional to 1/std.dev(response)spasm.smooth <- function(object,X,residual=FALSE,block=0) UseMethod("spasm.smooth")spasm.range <- function(object,upper.prop=.5,centre=TRUE) {## get reasonable smoothing parameter range for sparse smooth in objectsp <- 1edf <- spasm.sp(object,sp,get.trH=TRUE,centre=centre)$trHwhile (edf < object$edf0*1.01+.5) {sp <- sp /100edf <- spasm.sp(object,sp,get.trH=TRUE,centre=centre)$trH}sp1 <- sp ## store smallest known goodwhile (edf > object$edf0*1.01+.5) {sp <- sp * 100edf <- spasm.sp(object,sp,get.trH=TRUE,centre=centre)$trH}sp0 <- spwhile (edf < object$edf1*upper.prop) {sp1 <- sp1 / 100edf <- spasm.sp(object,sp1,get.trH=TRUE,centre=centre)$trH}while (edf > object$edf1*upper.prop) {sp1 <- sp1 * 4edf <- spasm.sp(object,sp1,get.trH=TRUE,centre=centre)$trH}c(sp1,sp0) ## small, large}