Nothing
R/misc.r
In mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation
Defines functions
Ztb
Zb
pdev
neicov
minres
AddBVB
isa
mat.rowsum
temp.seed
eigen.approx
lanczos
treig
pmmult
pqr.qy
pqr.R
pqr
block.reorder
pRRt
pcrossprod
pforwardsolve
pchol
pbsi
pqr2
pinv
vcorr
cholup
choldrop
dchol
diagXVXd
Xbd
XWyd
XWXd
mgcv.omp
trichol
bandchol
sdiag
rmvn
blas.thread.test
dpnorm
Documented in
bandchol
blas.thread.test
choldrop
cholup
diagXVXd
dpnorm
rmvn
sdiag
trichol
Xbd
XWXd
XWyd
## (c) Simon N. Wood 2011-2023
## Many of the following are simple wrappers for C functions
dpnorm <- function(x0,x1) {
## Cancellation avoiding evaluation of pnorm(x1)-pnorm(x0)
## first avoid 1-1 problems by exchanging and changing sign of double +ve
ii <- x1>0&x0>0
d <- x0[ii];x0[ii] <- -x1[ii];x1[ii] <- -d
## now deal with points that are so close that cancellation error
## too large - might as well use density times interval width
ii <- abs(x1-x0) < sqrt(.Machine$double.eps)*dnorm((x1+x0)/2)
p <- x0; d <- x1[ii]-x0[ii]; m <- (x1[ii]+x0[ii])/2
p[ii] <- dnorm(m)*d
p[!ii] <- pnorm(x1[!ii]) - pnorm(x0[!ii])
p
} ## dpnorm
"%.%" <- function(a,b) {
if (inherits(a,"dgCMatrix")||inherits(b,"dgCMatrix"))
tensor.prod.model.matrix(list( ## following is coercion to double, general (no special structure), compressed column
as(as(as(a, "dMatrix"), "generalMatrix"), "CsparseMatrix"),
as(as(as(b, "dMatrix"), "generalMatrix"), "CsparseMatrix")
)) else tensor.prod.model.matrix(list(as.matrix(a),as.matrix(b)))
# tensor.prod.model.matrix(list(as(a,"dgCMatrix"),as(b,"dgCMatrix"))) else - deprecated
# tensor.prod.model.matrix(list(as.matrix(a),as.matrix(b)))
}
blas.thread.test <- function(n=1000,nt=4) {
inter <- interactive()
if (inter) prg <- txtProgressBar(min = 0, max = n, initial = 0,
char = "=",width = NA, title="Progress", style = 3)
for (i in 1:n) {
m <- sample(213:654,1)
p <- sample(3:153,1)
X <- matrix(runif(m*p),m,p)
er <- pqr(X,nt=nt)
rr <- range(pqr.qy(er,pqr.R(er))-X[,er$pivot])
if (rr[1] < -1e-4||rr[2] > 1e-4) {
break;
}
if (inter) setTxtProgressBar(prg, i)
}
if (inter) close(prg)
if (rr[1] < -1e-4||rr[2] > 1e-4) {
cat("BLAS thread safety problem at iteration",i,"\n")
} else cat("No problem encountered in",i,"iterations\n")
} ## blas.thread.test
rmvn <- function(n,mu,V) {
## generate multivariate normal deviates. e.g.
## V <- matrix(c(2,1,1,2),2,2); mu <- c(1,1);n <- 1000;z <- rmvn(n,mu,V);crossprod(sweep(z,2,colMeans(z)))/n
p <- ncol(V)
R <- mroot(V,rank=ncol(V)) ## RR' = V
if (is.matrix(mu)) {
if (ncol(mu)!=p||nrow(mu)!=n) stop("mu dimensions wrong")
z <- matrix(rnorm(p*n),n,p)%*%t(R) + mu
} else {
if (length(mu)!=p) stop("mu dimensions wrong")
z <- t(R%*% matrix(rnorm(p*n),p,n) + mu)
if (n==1) z <- as.numeric(z)
}
z
} ## rmvn
sdiag <- function(A,k=0) {
## extract sub or super diagonal of matrix (k=0 is leading)
p <- ncol(A)
n <- nrow(A)
if (k>p-1||-k > n-1) return()
if (k >= 0) {
i <- 1:n
j <- (k+1):p
} else {
i <- (-k+1):n
j <- 1:p
}
if (length(i)>length(j)) i <- i[1:length(j)] else j <- j[1:length(i)]
ii <- i + (j-1) * n
A[ii]
} ## sdiag
"sdiag<-" <- function(A,k=0,value) {
p <- ncol(A)
n <- nrow(A)
if (k>p-1||-k > n-1) return()
if (k >= 0) {
i <- 1:n
j <- (k+1):p
} else {
i <- (-k+1):n
j <- 1:p
}
if (length(i)>length(j)) i <- i[1:length(j)] else j <- j[1:length(i)]
ii <- i + (j-1) * n
A[ii] <- value
A
} ## "sdiag<-"
bandchol <- function(B) {
## obtain R such that R'R = A. Where A is banded matrix contained in R.
n <- ncol(B)
k <- 0
if (n==nrow(B)) { ## square matrix. Extract the diagonals
A <- B*0
for (i in 1:n) {
b <- sdiag(B,i-1)
if (sum(b!=0)!=0) {
k <- i ## largest index of a non-zero band
A[i,1:length(b)] <- b
}
}
B <- A[1:k,]
}
oo <- .C(C_band_chol,B=as.double(B),n=as.integer(n),k=as.integer(nrow(B)),info=as.integer(0))
if (oo$info<0) stop("something wrong with inputs to LAPACK routine")
if (oo$info>0) stop("not positive definite")
B <- matrix(oo$B,nrow(B),n)
if (k>0) { ## was square on entry, so also on exit...
A <- A * 0
for (i in 1:k) sdiag(A,i-1) <- B[i,1:(n-i+1)]
B <- A
}
B
} ## bandchol
trichol <- function(ld,sd) {
## obtain chol factor R of symm tridiag matrix, A, with leading diag
## ld and sub/super diags sd. R'R = A. On exit ld is diag of R and
## sd its super diagonal.
n <- length(ld)
if (n<2) stop("don't be silly")
if (n!=length(sd)+1) stop("sd should have exactly one less entry than ld")
oo <- .C(C_tri_chol,ld=as.double(ld),sd=as.double(sd),n=as.integer(n),info=as.integer(0))
if (oo$info<0) stop("something wrong with inputs to LAPACK routine")
if (oo$info>0) stop("not positive definite")
ld <- sqrt(oo$ld)
sd <- oo$sd*ld[1:(n-1)]
list(ld=ld,sd=sd)
}
mgcv.omp <- function() {
## does open MP appear to be available?
oo <- .C(C_mgcv_omp,a=as.integer(-1))
if (oo$a==1) TRUE else FALSE
}
## discretized covariate routines...
XWXd <- function(X,w,k,ks,ts,dt,v,qc,nthreads=1,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1,lt=NULL,rt=NULL) {
## Form X'WX given weights in w and X in compressed form in list X.
## each element of X is a (marginal) model submatrix. Full version
## is given by X[[i]][k[,i],] (see below for summation convention).
## list X relates to length(ts) separate
## terms. ith term starts at matrix ts[i] and has dt[i] marginal matrices.
## For summation convention, k[,ks[j,1]:ks[j,2]] gives index columns
## for matrix j, thereby allowing summation over matrix covariates....
## i.e. for q in ks[j,1]:ks[j,2] sum up X[[j]][k[,q],]
## Terms with several marginals are tensor products and may have
## constraints (if qc[i]>1), stored as a householder vector in v[[i]].
## check ts and k index start (assumed 1 here)
## if drop is non-NULL it contains index of rows/cols to drop from result
## * lt is array of terms to include in left matrix (assumed in ascending coef index order)
## * rt is array of terms to include in right matrix (assumed in ascending coef index order)
## * if both NULL all terms are included, if only one is NULL then used for left and right.
m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
n <- length(w);ptfull <- pt <- 0;
for (i in 1:nt) {
fullsize <- prod(p[ts[i]:(ts[i]+dt[i]-1)])
ptfull <- ptfull + fullsize
pt <- pt + fullsize - if (qc[i]>0) 1 else if (qc[i]<0) v[[i]][v[[i]][1]+2] else 0
}
if (inherits(X[[1]],"dgCMatrix")) { ## the marginals are sparse
if (length(ar.stop)>1||ar.stop!=-1) warning("AR not available with sparse marginals")
## create list for passing to C
if (any(qc<0)) stop("sparse method for Kronecker product contrasts not implemented")
m <- list(Xd=X,kd=k,ks=ks,v=v,ts=ts,dt=dt,qc=qc)
m$off <- attr(X,"off"); m$r <- attr(X,"r")
if (is.null(m$off)||is.null(m$r)) stop("reverse indices missing from sparse discrete marginals")
### code that could create the marginal reverse indices...
#for (j in 1:nrow(m$ks)) {
#nr <- nrow(m$Xd[[j]]) ## make sure we always tab to final stored row
#for (i in m$ks[j,1]:(m$ks[j,2]-1)) {
# m$r[,i] <- (1:length(m$kd[,i]))[order(m$kd[,i])]
# m$off[[i]] <- cumsum(c(1,tabulate(m$kd[,i],nbins=nr)))-1
#}
m$offstart <- cumsum(c(0,lapply(m$off,length)))
m$off <- unlist(m$off)
## Now C base all indices...
m$ks <- m$ks - 1; m$kd <- m$kd - 1; m$r <- m$r - 1; m$ts <- m$ts-1
nthreads <- as.integer(nthreads)
w <- as.double(w)
if ((!is.null(lt)||!is.null(rt))&&!is.null(drop)) {
lpip <- attr(X,"lpip") ## list of coefs for each term
rpi <- unlist(lpip[rt])
lpi <- unlist(lpip[lt])
if (is.null(lt)) lpi <- rpi
else if (is.null(rt)) rpi <- lpi
ldrop <- which(lpi %in% drop)
rdrop <- which(lpi %in% drop)
} else rdrop <- ldrop <- drop
if (is.null(lt)&&is.null(rt)) {
lt <- rt <- 1:nt
} else if (is.null(lt)) {
lt <- rt
} else if (is.null(rt)) rt <- lt;
lt <- as.integer(lt-1)
rt <- as.integer(rt-1)
XWX <- .Call(C_sXWXd,m,w,lt,rt,nthreads)
if (!is.null(drop)) {
Dl <- Diagonal(ncol(XWX),1)
XWX <- Dl[-ldrop,] %*% XWX %*% t(Dl[-rdrop,])
}
return(XWX) ## note that this is sparse
}
## block oriented code...
if (is.null(lt)&&is.null(lt)) {
# old .C code - can't handle long vector k
#oo <- .C(C_XWXd0,XWX =as.double(rep(0,ptfull^2)),X= as.double(unlist(X)),w=as.double(w),
# k=as.integer(k-1),ks=as.integer(ks-1),m=as.integer(m),p=as.integer(p), n=as.integer(n),
# ns=as.integer(nx), ts=as.integer(ts-1), as.integer(dt), nt=as.integer(nt),
# v = as.double(unlist(v)),qc=as.integer(qc),nthreads=as.integer(nthreads),
# ar.stop=as.integer(ar.stop-1),ar.weights=as.double(ar.w))
#XWX <- if (is.null(drop)) matrix(oo$XWX[1:pt^2],pt,pt) else matrix(oo$XWX[1:pt^2],pt,pt)[-drop,-drop]
XWX <- numeric(ptfull^2)
.Call(C_CXWXd0,XWX,as.double(unlist(X)),w,k-1L,as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(ts-1L), as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(nthreads),
as.integer(ar.stop-1L),as.double(ar.w))
XWX <- if (is.null(drop)) matrix(XWX[1:pt^2],pt,pt) else matrix(XWX[1:pt^2],pt,pt)[-drop,-drop]
} else {
lpip <- attr(X,"lpip") ## list of coefs for each term
rpi <- unlist(lpip[rt])
lpi <- unlist(lpip[lt])
if (is.null(lt)) {
lpi <- rpi
nrs <- lt <- 0
ncs <- length(rt)
} else {
nrs <- length(lt)
if (is.null(rt)) {
rt <- ncs <- 0
rpi <- lpi
} else ncs <- length(rt)
}
#oo <- .C(C_XWXd1,XWX =as.double(rep(0,ptfull^2)),X= as.double(unlist(X)),w=as.double(w),
# k=as.integer(k-1),ks=as.integer(ks-1),m=as.integer(m),p=as.integer(p), n=as.integer(n),
# ns=as.integer(nx), ts=as.integer(ts-1), dt=as.integer(dt), nt=as.integer(nt),
# v = as.double(unlist(v)),qc=as.integer(qc),nthreads=as.integer(nthreads),
# ar.stop=as.integer(ar.stop-1),ar.weights=as.double(ar.w),rs=as.integer(lt-1),
# cs=as.integer(rt-1),nrs=as.integer(nrs),ncs=as.integer(ncs))#)
#XWX <- matrix(oo$XWX[1:(length(lpi)*length(rpi))],length(lpi),length(rpi))
XWX <- numeric(ptfull^2)
.Call(C_CXWXd1,XWX,as.double(unlist(X)),w,k-1L,as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(ts-1L), as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(nthreads),
as.integer(ar.stop-1L),as.double(ar.w),as.integer(lt-1L),as.integer(rt-1L))
XWX <- matrix(XWX[1:(length(lpi)*length(rpi))],length(lpi),length(rpi))
if (!is.null(drop)) {
ldrop <- which(lpi %in% drop)
rdrop <- which(lpi %in% drop)
if (length(ldrop)>0||length(rdrop)>0) XWX <-
if (length(ldrop==0)) XWX[,-rdrop] else if (length(rdrop)==0) XWX[-ldrop,] else XWX[-ldrop,-rdrop]
}
}
XWX
} ## XWXd
XWyd <- function(X,w,y,k,ks,ts,dt,v,qc,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1,lt=NULL) {
## X'Wy...
## if lt if not NULL then it lists the discrete terms to include (from X)
## returned vector/matrix only includes rows for selected terms
m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
n <- length(w);
if (is.null(lt)) {
pt <- 0
for (i in 1:nt) pt <- pt + prod(p[ts[i]:(ts[i]+dt[i]-1)]) - if (qc[i]>0) 1 else if (qc[i]<0) v[[i]][v[[i]][1]+2] else 0
lt <- 1:nt
} else {
lpip <- attr(X,"lpip") ## list of coefs for each term
lpi <- unlist(lpip[lt]) ## coefs corresponding to terms selected by lt
if (!is.null(drop)) drop <- which(lpi %in% drop) ## rebase drop
pt <- length(lpi)
}
cy <- if (is.matrix(y)) ncol(y) else 1
if (inherits(X[[1]],"dgCMatrix")) { ## the marginals are sparse
## create list for passing to C
m <- list(Xd=X,kd=k,ks=ks,v=v,ts=ts,dt=dt,qc=qc)
m$off <- attr(X,"off"); m$r <- attr(X,"r")
if (is.null(m$off)||is.null(m$r)) stop("reverse indices missing from sparse discrete marginals")
m$offstart <- cumsum(c(0,lapply(m$off,length)))
m$off <- unlist(m$off)
## Now C base all indices...
m$ks <- m$ks - 1; m$kd <- m$kd - 1; m$r <- m$r - 1; m$ts <- m$ts-1
Wy <- as.double(w*y);lt <- as.integer(lt-1)
XWy <- .Call(C_sXyd,m,Wy,lt)
if (cy>1) XWy <- matrix(XWy,ncol=cy)
if (!is.null(drop)) XWy <- if (cy>1) XWy[-drop,] else XWy[-drop]
} else { ## dense marginals case
## old .C code - can't handle long vector k
#oo <- .C(C_XWyd,XWy=rep(0,pt*cy),y=as.double(y),X=as.double(unlist(X)),w=as.double(w),k=as.integer(k-1),
# ks=as.integer(ks-1),
# m=as.integer(m),p=as.integer(p),n=as.integer(n),cy=as.integer(cy), nx=as.integer(nx), ts=as.integer(ts-1),
# dt=as.integer(dt),nt=as.integer(nt),v=as.double(unlist(v)),qc=as.integer(qc),
# ar.stop=as.integer(ar.stop-1),ar.row=as.integer(ar.row-1),ar.weights=as.double(ar.w),
# cs=as.integer(lt-1),ncs=as.integer(length(lt)))
#if (cy>1) XWy <- if (is.null(drop)) matrix(oo$XWy,pt,cy) else matrix(oo$XWy,pt,cy)[-drop,] else
#XWy <- if (is.null(drop)) oo$XWy else oo$XWy[-drop]
XWy <- numeric(pt*cy)
.Call(C_CXWyd,XWy,as.double(y),as.double(unlist(X)),as.double(w),k-1L,as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(cy),as.integer(ts-1L), as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(ar.stop-1L),
as.integer(ar.row-1L),as.double(ar.w),as.integer(lt-1L))
if (cy>1) { XWy <- if (is.null(drop)) matrix(XWy,pt,cy) else matrix(XWy,pt,cy)[-drop,] } else {
XWy <- if (is.null(drop)) XWy else XWy[-drop]
}
}
XWy
} ## XWyd
Xbd <- function(X,beta,k,ks,ts,dt,v,qc,drop=NULL,lt=NULL) {
## note that drop may contain the index of columns of X to drop before multiplying by beta.
## equivalently we can insert zero elements into beta in the appropriate places.
## if lt if not NULL then it lists the discrete terms to include (from X)
n <- if (is.matrix(k)) nrow(k) else length(k) ## number of data
m <- unlist(lapply(X,nrow)) ## number of rows in each discrete model matrix
p <- unlist(lapply(X,ncol)) ## number of cols in each discrete model matrix
nx <- length(X) ## number of model matrices
if (length(p)!=nx) stop("something wrong with matrix list - not all matrices?")
nt <- length(ts) ## number of terms
if (!is.null(drop)) {
b <- if (is.matrix(beta)) matrix(0,nrow(beta)+length(drop),ncol(beta)) else rep(0,length(beta)+length(drop))
if (is.matrix(beta)) b[-drop,] <- beta else b[-drop] <- beta
beta <- b
}
if (is.null(lt)) lt <- 1:nt
bc <- if (is.matrix(beta)) ncol(beta) else 1 ## number of columns in beta
## The C code mechanism for dealing with lt is very basic, and requires that beta is re-ordered and
## truncated to relate only to the selected terms, in the order they are selected.
lpip <- attr(X,"lpip")
if (!is.null(lpip)) { ## then X list may not be in coef order...
lpip <- unlist(lpip[lt])
beta <- if (is.matrix(beta)) beta[lpip,] else beta[lpip] ## select params required in correct order
}
if (inherits(X[[1]],"dgCMatrix")) { ## the marginals are sparse
## create list for passing to C
m <- list(Xd=X,kd=k,ks=ks,v=v,ts=ts,dt=dt,qc=qc)
m$off <- attr(X,"off"); m$r <- attr(X,"r")
if (is.null(m$off)||is.null(m$r)) stop("reverse indices missing from sparse discrete marginals")
m$offstart <- cumsum(c(0,lapply(m$off,length)))
m$off <- unlist(m$off)
## Now C base all indices...
m$ks <- m$ks - 1; m$kd <- m$kd - 1; m$r <- m$r - 1; m$ts <- m$ts-1
beta <- as.matrix(beta);storage.mode(beta) <- "double"
lt <- as.integer(lt-1)
Xb <- .Call(C_sXbd,m,beta,lt)
if (bc>1) Xb <- matrix(Xb,ncol=bc)
} else { ## dense marginals case
#oo <- .C(C_Xbd,f=as.double(rep(0,n*bc)),beta=as.double(beta),X=as.double(unlist(X)),k=as.integer(k-1),
# ks = as.integer(ks-1),
# m=as.integer(m),p=as.integer(p), n=as.integer(n), nx=as.integer(nx), ts=as.integer(ts-1),
# as.integer(dt), as.integer(nt),as.double(unlist(v)),as.integer(qc),as.integer(bc),as.integer(lt-1),as.integer(length(lt)))
#Xb <- if (is.matrix(beta)) matrix(oo$f,n,bc) else oo$f
f <- numeric(n*bc)
.Call(C_CXbd,f,beta,as.double(unlist(X)),k-1L, as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(ts-1L),as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(bc), as.integer(lt-1L))
Xb <- if (is.matrix(beta)) matrix(f,n,bc) else f
}
return(Xb)
} ## Xbd
diagXVXd <- function(X,V,k,ks,ts,dt,v,qc,drop=NULL,nthreads=1,lt=NULL,rt=NULL) {
## discrete computation of diag(XVX')
n <- if (is.matrix(k)) nrow(k) else length(k)
m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
if (is.null(lt)&&is.null(rt)) rt <- lt <- 1:nt else {
if (is.null(rt)) rt <- lt else if (is.null(lt)) lt <- rt
}
if (is.null(rt)) rt <- 1:nt
if (inherits(X[[1]],"dgCMatrix")) { ## the marginals are sparse
## create list for passing to C
m <- list(Xd=X,kd=k,ks=ks,v=v,ts=ts,dt=dt,qc=qc)
m$off <- attr(X,"off"); m$r <- attr(X,"r")
if (is.null(m$off)||is.null(m$r)) stop("reverse indices missing from sparse discrete marginals")
m$offstart <- cumsum(c(0,lapply(m$off,length)))
m$off <- unlist(m$off)
## Now C base all indices...
m$ks <- m$ks - 1; m$kd <- m$kd - 1; m$r <- m$r - 1; m$ts <- m$ts-1
if (!is.null(drop)) {
D <- Diagonal(ncol(V)+length(drop),1)[-drop,]
V <- t(D) %*% V %*% D
}
## The C code mechanism for dealing with rt and lt is very basic, and requires that V is
## re-ordered and truncated to relate only to the selected terms, in the order they are selected.
lpip <- attr(X,"lpip")
if (!is.null(lpip)) { ## then X list may not be in coef order...
lpi <- unlist(lpip[lt])
rpi <- unlist(lpip[rt])
V <- V[lpi,rpi,drop=FALSE] ## select part of V required in correct order
}
lt <- as.integer(lt-1);rt <- as.integer(rt-1)
D <- .Call(C_sdiagXVXt,m , V, lt, rt)
} else { ## dense marginals
if (!is.null(drop)) {
pv <- ncol(V)+length(drop)
V0 <- matrix(0,pv,pv)
V0[-drop,-drop] <- V
V <- V0;rm(V0)
} else pv <- ncol(V)
## The C code mechanism for dealing with rt and lt is very basic, and requires that V is
## re-ordered and truncated to relate only to the selected terms, in the order they are selected.
lpip <- attr(X,"lpip")
if (!is.null(lpip)) { ## then X list may not be in coef order...
lpi <- unlist(lpip[lt])
rpi <- unlist(lpip[rt])
V <- V[lpi,rpi,drop=FALSE] ## select part of V required in correct order
}
# oo <- .C(C_diagXVXt,diag=as.double(rep(0,n)),V=as.double(V),X=as.double(unlist(X)),k=as.integer(k-1),
# ks=as.integer(ks-1),m=as.integer(m),p=as.integer(p), n=as.integer(n), nx=as.integer(nx),
# ts=as.integer(ts-1), as.integer(dt), as.integer(nt),as.double(unlist(v)),as.integer(qc),as.integer(nrow(V)),
# as.integer(ncol(V)),as.integer(nthreads),as.integer(lt-1),as.integer(length(lt)),as.integer(rt-1),as.integer(length(rt)))
# D <- oo$diag
D <- numeric(n)
.Call(C_CdiagXVXt,D,V,as.double(unlist(X)),k-1L,as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(ts-1L), as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(nthreads),
as.integer(lt-1L),as.integer(rt-1L))
}
D
} ## diagXVXd
dchol <- function(dA,R) {
## if dA contains matrix dA/dx where R is chol factor s.t. R'R = A
## then this routine returns dR/dx...
p <- ncol(R)
oo <- .C(C_dchol,dA=as.double(dA),R=as.double(R),dR=as.double(R*0),p=as.integer(ncol(R)))
return(matrix(oo$dR,p,p))
} ## dchol
choldrop <- function(R,k) {
## routine to update Cholesky factor R of A on dropping row/col k of A.
## R can be upper triangular, in which case (R'R=A) or lower triangular in
## which case RR'=A...
n <- as.integer(ncol(R))
k1 <- as.integer(k-1)
ut <- as.integer(as.numeric(R[1,2]!=0))
if (k<1||k>n) return(R)
Rup <- matrix(0,n-1,n-1)
#oo <- .C(C_chol_down,R=as.double(R),Rup=as.double(Rup),n=as.integer(n),k=as.integer(k-1),ut=as.integer(ut))
.Call(C_mgcv_chol_down,R,Rup,n,k1,ut)
#matrix(oo$Rup,n-1,n-1)
Rup
} ## choldrop
cholup <- function(R,u,up=TRUE) {
## routine to update Cholesky factor R to the factor of R'R + uu' (up == TRUE)
## or R'R - uu' (up=FALSE).
n <- as.integer(ncol(R))
up <- as.integer(up)
eps <- as.double(.Machine$double.eps)
R1 <- R * 1.0
.Call(C_mgcv_chol_up,R1,u,n,up,eps)
if (up==0) if ((n>1 && R1[2,1] < -1)||(n==1&&u[1]>R[1])) stop("update not positive definite")
R1
} ## cholup
vcorr <- function(dR,Vr,trans=TRUE) {
## Suppose b = sum_k op(dR[[k]])%*%z*r_k, z ~ N(0,Ip), r ~ N(0,Vr). vcorr returns cov(b).
## dR is a list of p by p matrices. 'op' is 't' if trans=TRUE and I() otherwise.
p <- ncol(dR[[1]])
M <- if (trans) ncol(Vr) else -ncol(Vr) ## sign signals transpose or not to C code
if (abs(M)!=length(dR)) stop("internal error in vcorr, please report to simon.wood@r-project.org")
oo <- .C(C_vcorr,dR=as.double(unlist(dR)),Vr=as.double(Vr),Vb=as.double(rep(0,p*p)),
p=as.integer(p),M=as.integer(M))
return(matrix(oo$Vb,p,p))
} ## vcorr
pinv <- function(X,svd=FALSE) {
## a pseudoinverse for n by p, n>p matrices
qrx <- qr(X,tol=0,LAPACK=TRUE)
R <- qr.R(qrx);Q <- qr.Q(qrx)
rr <- Rrank(R)
if (svd&&rr<ncol(R)) {
piv <- 1:ncol(X); piv[qrx$pivot] <- 1:ncol(X)
er <- svd(R[,piv])
d <- er$d*0;d[1:rr] <- 1/er$d[1:rr]
X <- Q%*%er$u%*%(d*t(er$v))
} else {
Ri <- R*0
Ri[1:rr,1:rr] <- backsolve(R[1:rr,1:rr],diag(rr))
X[,qrx$pivot] <- Q%*%t(Ri)
}
X
} ## end pinv
pqr2 <- function(x,nt=1,nb=30) {
## Function for parallel pivoted qr decomposition of a matrix using LAPACK
## householder routines. Currently uses a block algorithm.
## library(mgcv); n <- 4000;p<-3000;x <- matrix(runif(n*p),n,p)
## system.time(qrx <- qr(x,LAPACK=TRUE))
## system.time(qrx2 <- mgcv:::pqr2(x,2))
## system.time(qrx3 <- mgcv:::pqr(x,2))
## range(qrx2$qr-qrx$qr)
p <- ncol(x)
beta <- rep(0.0,p)
piv <- as.integer(rep(0,p))
## need to force a copy of x, otherwise x will be over-written
## by .Call *in environment from which function is called*
x <- x*1
rank <- .Call(C_mgcv_Rpiqr,x,beta,piv,nt,nb)
ret <- list(qr=x,rank=rank,qraux=beta,pivot=piv+1)
attr(ret,"useLAPACK") <- TRUE
class(ret) <- "qr"
ret
} ## pqr2
pbsi <- function(R,nt=1,copy=TRUE) {
## parallel back substitution inversion of upper triangular R
## library(mgcv); n <- 500;p<-400;x <- matrix(runif(n*p),n,p)
## qrx <- qr(x);R <- qr.R(qrx)
## system.time(Ri <- mgcv:::pbsi(R,2))
## system.time(Ri2 <- backsolve(R,diag(p)));range(Ri-Ri2)
if (copy) R <- R * 1 ## ensure that R modified only within pbsi
.Call(C_mgcv_Rpbsi,R,nt)
R
} ## pbsi
pchol <- function(A,nt=1,nb=40) {
## parallel Choleski factorization.
## library(mgcv);
## set.seed(2);n <- 200;r <- 190;A <- tcrossprod(matrix(runif(n*r),n,r))
## system.time(R <- chol(A,pivot=TRUE));system.time(L <- mgcv:::pchol(A));range(R[1:r,]-L[1:r,])
## k <- 30;range(R[1:k,1:k]-L[1:k,1:k])
## system.time(L <- mgcv:::pchol(A,nt=2,nb=30))
## piv <- attr(L,"pivot");attr(L,"rank");range(crossprod(L)-A[piv,piv])
## should nb be obtained from 'ILAENV' as page 23 of Lucas 2004??
piv <- as.integer(rep(0,ncol(A)))
A <- A*1 ## otherwise over-write in calling env!
rank <- .Call(C_mgcv_Rpchol,A,piv,nt,nb)
attr(A,"pivot") <- piv+1;attr(A,"rank") <- rank
A
}
pforwardsolve <- function(R,B,nt=1) {
## parallel forward solve via simple col splitting...
if (!is.matrix(B)) B <- as.matrix(B)
.Call(C_mgcv_Rpforwardsolve,R,B,nt)
}
pcrossprod <- function(A,trans=FALSE,nt=1,nb=30) {
## parallel cross prod A'A or AA' if trans==TRUE...
if (!is.matrix(A)) A <- as.matrix(A)
if (trans) A <- t(A)
.Call(C_mgcv_Rpcross,A,nt,nb)
}
pRRt <- function(R,nt=1) {
## parallel RR' for upper triangular R
## following creates index of lower triangular elements...
## n <- 4000;a <- rep(1:n,n);b <- rep(1:n,each=n);which(a>=b) -> ii;a[ii]+(b[ii]-1)*n->ii ## lower
## n <- 4000;a <- rep(1:n,n);b <- rep(1:n,each=n);which(a<=b) -> ii;a[ii]+(b[ii]-1)*n->ii ## upper
## library(mgcv);R <- matrix(0,n,n);R[ii] <- runif(n*(n+1)/2)
## Note: A[a-b<=0] <- 0 zeroes upper triangle
## system.time(A <- mgcv:::pRRt(R,2))
## system.time(A2 <- tcrossprod(R));range(A-A2)
n <- nrow(R)
A <- matrix(0,n,n)
.Call(C_mgcv_RPPt,A,R,nt)
A
}
block.reorder <- function(x,n.blocks=1,reverse=FALSE) {
## takes a matrix x divides it into n.blocks row-wise blocks, and re-orders
## so that the blocks are stored one after the other.
## e.g. library(mgcv); x <- matrix(1:18,6,3);xb <- mgcv:::block.reorder(x,2)
## x;xb;mgcv:::block.reorder(xb,2,TRUE)
r = nrow(x);cols = ncol(x);
if (n.blocks <= 1) return(x);
if (r%%n.blocks) {
nb = ceiling(r/n.blocks)
} else nb = r/n.blocks;
oo <- .C(C_row_block_reorder,x=as.double(x),as.integer(r),as.integer(cols),
as.integer(nb),as.integer(reverse));
matrix(oo$x,r,cols)
} ## block.reorder
pqr <- function(x,nt=1) {
## parallel QR decomposition, using openMP in C, and up to nt threads (only if worthwhile)
## library(mgcv);n <- 20;p<-4;X <- matrix(runif(n*p),n,p);er <- mgcv:::pqr(X,nt=2)
## range(mgcv:::pqr.qy(er,mgcv:::pqr.R(er))-X[,er$pivot])
x.c <- ncol(x);r <- nrow(x)
oo <- .C(C_mgcv_pqr,x=as.double(c(x,rep(0,nt*x.c^2))),as.integer(r),as.integer(x.c),
pivot=as.integer(rep(0,x.c)), tau=as.double(rep(0,(nt+1)*x.c)),as.integer(nt))
list(x=oo$x,r=r,c=x.c,tau=oo$tau,pivot=oo$pivot+1,nt=nt)
}
pqr.R <- function(x) {
## x is an object returned by pqr. This extracts the R factor...
## e.g. as pqr then...
## R <- mgcv:::pqr.R(er); R0 <- qr.R(qr(X,tol=0))
## svd(R)$d;svd(R0)$d
oo <- .C(C_getRpqr,R=as.double(rep(0,x$c^2)),as.double(x$x),as.integer(x$r),as.integer(x$c),
as.integer(x$c),as.integer(x$nt))
matrix(oo$R,x$c,x$c)
}
pqr.qy <- function(x,a,tr=FALSE) {
## x contains a parallel QR decomp as computed by pqr. a is a matrix. computes
## Qa or Q'a depending on tr.
## e.g. as above, then...
## a <- diag(p);Q <- mgcv:::pqr.qy(er,a);crossprod(Q)
## X[,er$pivot+1];Q%*%R
## Qt <- mgcv:::pqr.qy(er,diag(n),TRUE);Qt%*%t(Qt);range(Q-t(Qt))
## Q <- qr.Q(qr(X,tol=0));z <- runif(n);y0<-t(Q)%*%z
## mgcv:::pqr.qy(er,z,TRUE)->y
## z <- runif(p);y0<-Q%*%z;mgcv:::pqr.qy(er,z)->y
if (is.matrix(a)) a.c <- ncol(a) else a.c <- 1
if (tr) {
if (is.matrix(a)) { if (nrow(a) != x$r) stop("a has wrong number of rows") }
else if (length(a) != x$r) stop("a has wrong number of rows")
} else {
if (is.matrix(a)) { if (nrow(a) != x$c) stop("a has wrong number of rows") }
else if (length(a) != x$c) stop("a has wrong number of rows")
a <- c(a,rep(0,a.c*(x$r-x$c)))
}
oo <- .C(C_mgcv_pqrqy,a=as.double(a),as.double(x$x),as.double(x$tau),as.integer(x$r),
as.integer(x$c),as.integer(a.c),as.integer(tr),as.integer(x$nt))
if (tr) return(matrix(oo$a[1:(a.c*x$c)],x$c,a.c)) else
return(matrix(oo$a,x$r,a.c))
}
pmmult <- function(A,B,tA=FALSE,tB=FALSE,nt=1) {
## parallel matrix multiplication (not for use on vectors or thin matrices)
## library(mgcv);r <- 10;c <- 5;n <- 8
## A <- matrix(runif(r*n),r,n);B <- matrix(runif(n*c),n,c);range(A%*%B-mgcv:::pmmult(A,B,nt=1))
## A <- matrix(runif(r*n),n,r);B <- matrix(runif(n*c),n,c);range(t(A)%*%B-mgcv:::pmmult(A,B,TRUE,FALSE,nt=1))
## A <- matrix(runif(r*n),n,r);B <- matrix(runif(n*c),c,n);range(t(A)%*%t(B)-mgcv:::pmmult(A,B,TRUE,TRUE,nt=1))
## A <- matrix(runif(r*n),r,n);B <- matrix(runif(n*c),c,n);range(A%*%t(B)-mgcv:::pmmult(A,B,FALSE,TRUE,nt=1))
if (tA) { n = nrow(A);r = ncol(A)} else {n = ncol(A);r = nrow(A)}
if (tB) { c = nrow(B)} else {c = ncol(B)}
C <- rep(0,r * c)
oo <- .C(C_mgcv_pmmult,C=as.double(C),as.double(A),as.double(B),as.integer(tA),as.integer(tB),as.integer(r),
as.integer(c),as.integer(n),as.integer(nt));
matrix(oo$C,r,c)
}
treig <- function(ld,sd,vec=FALSE,descend=FALSE) {
## eigen decomposition of tri-diagonal matrix with leading diagonal ld and
## sub-diagonal sd...
n <- length(ld)
v <- if (vec) rep(0,n*n) else 0
oo <- .C(C_mgcv_trisymeig,d=as.double(ld),g=as.double(sd),v=as.double(v),n=as.integer(n),
get.vec=as.integer(vec),descending=as.integer(descend))
v <- if (vec) matrix(oo$v,n,n) else NA
list(values=oo$d,vectors=v)
} ## treig
lanczos <- function(A,v0,M,Av = function(A,v) A%*%v,n=ncol(A)) {
## Apply M steps of Lanczos starting at v0 for n by n +ve semi definite matrix A.
## Av is a function forming the product of matrix A with vector v.
## A can be a matrix, in which case the default Av applies, or
## it could be a list of arguments used to define the multiplication
## in some other way (e.g. as a sequence of matrix products) defined
## in a custom Av...
## The function is used to find an approximate eigen value CDF for A
## as described in Lin, Saad and Yang (2016) SIAM Review 58(1), 34-65
## section 3.2.1 in particular.
v0.norm <- sqrt(sum(v0^2))
gamma <- epsilon <- rep(0,M)
q <- matrix(v0/v0.norm,n,M)
for (j in 1:M) {
c <- Av(A,q[,j])
if (j==1) vAv <- sum(v0*c)*v0.norm ## v0'Av0 - useful for estimating tr(A)
gamma[j] <- sum(c*q[,j])
c <- c - gamma[j]*q[,j]
if (j>1) {
c <- c - epsilon[j-1]*q[,j-1]
cq <- drop(t(c) %*% q[,1:j,drop=FALSE])
c <- c - colSums(cq*t(q[,1:j]))
cq <- drop(t(c) %*% q[,1:j,drop=FALSE])
c <- c - colSums(cq*t(q[,1:j]))
}
epsilon[j] <- sqrt(sum(c^2))
if (j<M) q[,j+1] <- c/epsilon[j]
}
et <- treig(gamma,epsilon,descend=FALSE,vec=TRUE)
## compute error bounds on the eigenvalues of
## A using the method described in section 3.2 of
## Parlett, BN (1998) The Symmetric Eigenvalue Problem, SIAM
err <- abs(et$vectors[M,])*epsilon[M]
theta <- c(0,et$values)
tau <- et$vectors[1,]
eta <- c(0,cumsum(tau^2))
## theta is vector of eigenvalues at which CDF jumps
## tau^2 is jump size, eta is CDF at theta
## lam.ub is upper bound on larget eigenvalue
list(theta=theta,tau=tau,eta=eta,err=err,vAv=vAv)
} ## lanczos
eigen.approx <- function(A,Av = function(A,v) A%*%v,M=20,n.rep=20,n=ncol(A),seed=1) {
## get the approximate eigenvalues of n by n +ve semi def matrix A. Av is
## the function for multiplying a vector by the matrix defined by A. If A
## is simply a matrix then the default Av is sufficient. M is the number of
## Lanczos steps to use, and n.rep the number of random replicates to average
## over. Routine restores RNG to pre-call state on exit.
## Based on Lin, Saad and Yang (2016) SIAM Review 58(1), 34-65
## section 3.2.1, but extended to only use this approximation for the
## eigen-values that have yet to converge.
## The CDF approximation is based on their Appendix C proposal, rather than
## using a Gausssian kernel approximation to the pdf and cdf and then
## inverting by tabulation. This is because the latter tends to oversmooth the
## CDF in a way that is unhelpful for rank deficient matrices.
## The Gaussian kernel approach would probably be prefereable for full rank matrices,
## since it then benefits from the extra stability of kernel smoothing.
if (is.finite(seed)) a <- temp.seed(seed) ## seed RNG and store state
eva <- rep(0,n)
trA <- rep(0,n.rep)
tol <- .Machine$double.eps^.5
for (r in 1:n.rep) {
v0 <- rnorm(n)
lz <- lanczos(A,v0,M=M,Av=Av,n=n)
trA[r] <- lz$vAv
## following is suggested in Appendix C of LSY, and is quite important
## to avoid slight downward bias...
eta1 <- c(0,(lz$eta[1:M]*0.5+0.5*lz$eta[1:M+1]))
eta1[2] <- lz$eta[2] ## correction to avoid over-estimation in lower tail if rank def
conv <- lz$err<lz$theta[M+1]*tol ## these eigenvalues are converged
upper.uconv <- if (any(!conv)) max(which(!conv)) else 0 ## last uncoverged
n.conv <- M - upper.uconv ## number converged
lz$theta[lz$theta<0] <- 0
theta.conv <- if (n.conv) lz$theta[(upper.uconv+1):M+1] else rep(0,0)
if (upper.uconv) {
eta <- eta1[1:(upper.uconv+1)]
theta <- lz$theta[1:(upper.uconv+1)]
eta <- eta/max(eta)
nri <- c(diff(eta)!=0,TRUE) ## strip out duplicates
eva <- eva + c(approx(eta[nri],theta[nri],seq(0,1,length=n-n.conv),method="linear",rule=2)$y,theta.conv)
} else {
eva <- eva + c(rep(0,n-n.conv),theta.conv)
}
}
if (is.finite(seed)) temp.seed(a) ## restore RNG state
eva <- eva/n.rep
trA.sd <- sd(trA)/sqrt(n.rep);trA <- mean(trA)
if (abs(sum(eva)-trA)>2.5*trA.sd) { ## evidence for bias in eigen-spectrum
eva <- eva*trA/sum(eva) ## correction
}
eva
} ## eigen.approx
temp.seed <- function(x) {
## when called with a numeric x stores the state of the RNG and sets its seed to
## x. Returns an object of class "rng.state". When called with an object of this
## class created on a previus call to this function, resets the RNG to the state
## on entry to that previous call.
if (inherits(x,"rng.state")) {
RNGkind(x$kind[1],x$kind[2])
assign(".Random.seed",x$seed,envir=.GlobalEnv) ## RNG behaves as if it had not been used
} else {
seed <- try(get(".Random.seed",envir=.GlobalEnv),silent=TRUE) ## store RNG seed
if (inherits(seed,"try-error")) {
runif(1)
seed <- get(".Random.seed",envir=.GlobalEnv)
}
kind <- RNGkind(NULL)
RNGkind("default","default")
set.seed(x) ## ensure repeatability
x <- list(seed=seed,kind=kind)
class(x) <- "rng.state"
return(x)
}
} ## temp.seed
mat.rowsum <- function(X,m,k) {
## Let X be n by p and m of length M. Produces an M by p matrix B
## where the ith row of B is the sum of the rows X[k[j],] where
## j = (m[i-1]+1):m[i]. m[0] is taken as 0.
## n <- 10;p <- 5;X <- matrix(runif(n*p),n,p)
## m <- c(3,5,8,11);k <- c(1,4,3,6,1,5,7,10,9,5,6)
## mgcv:::mat.rowsum(X,m,k)
if (max(k)>nrow(X)||min(k)<1) stop("index vector has invalid entries")
k <- k - 1 ## R to C index conversion
.Call(C_mrow_sum,X,m,k)
} ## mat.rowsum
isa <- function(R,nt=1) {
## Finds the elements of (R'R)^{-1} on NZP(R+R').
if (!inherits(R,c("dgCMatrix","dtCMatrix"))) stop("isa requires a dg/tCMatrix")
nt <- round(nt)
if (nt<1) nt = 1
Hpi <- R + t(R)
.Call(C_isa1p,t(R),Hpi,nt)
Hpi
} ## isa
AddBVB <- function(A,Bt,VBt) {
## Add B %*% V %*% t(B) to calss 'dgCMatrix' A returning result on NZP(A) only
## (i.e. discarding elements of BVB' not in NZP(A)), Bt is the transpose
## of B. B and VBt are class 'matrix'
A@x <- A@x * 1.0 ## force copy, otherwise A and return value modified
.Call(C_AddBVB,A,Bt,VBt)
A
} ## AddBVB
minres <- function(R,u,b) {
## routine to solve (R'R-uu')x = b using minres algorithm.
## set.seed(0);n <- 100;p <- 20;X <- matrix(runif(n*p)-.5,n,p);R <- chol(crossprod(X));b <- runif(p);k <- 1;
## solve(crossprod(X[-k,]),b);mgcv:::minres(R,t(X[k,]),b)
x <- b; p <- length(b);
m <- if (is.matrix(u)) ncol(u) else 1
work <- rep(0,p*(m+7)+m)
oo <- .C(C_minres,R=as.double(R), u=as.double(u),b=as.double(b), x=as.double(x), p=as.integer(p),m=as.integer(m),work=as.double(work))
cat("\n niter : ",oo$m,"\n")
oo$x
}
neicov <- function(Dd,nei) {
## wrapper for nei_cov. Dd is n by p matrix of leave one out perturbations to
## coef vectors. nei is neighbourhood structure.
p <- ncol(Dd)
V <- matrix(0,p,p)
k <- nei$k-1
.Call(C_nei_cov,V,Dd,nei$m,k)
(V+t(V))/2
} ## neicov
## Routines to force matrix to pdef
pdev <- function(A) {
## Force A to meet necessary conditions for +ve def from Thm 4.2.8. of Golub and van Load 4th ed.
## Non-positive diagonals are set to diagonal dominance. Off diagonals are then set to just meet
## necessary conditions.
## NOTE: A is modified directly in situ, so care needed to force a copy to be kept if it's needed
## A0 <- A will not work, as R only actually copies when R modifies!
## A0 <- A; A0[1,1] <- A0[1,1] + 1;A0[1,1] <- A0[1,1] - 1 works.
if (inherits(A,"Matrix")) {
if (!inherits(A,"dgCMatrix")) A <- as(as(as(A, "dMatrix"), "generalMatrix"), "CsparseMatrix")
da <- diag(A); ii <- which(da==0)
if (length(ii)) diag(A)[ii] <- -1e-30 ## Avoid C code having to insert extra non-zeroes
mod <- .Call(C_spdev,A)
} else { ## dense matrix
mod <- .Call(C_dpdev,A)
}
if (mod>0) attr(A,"modified") <- mod
return(A)
}
## following are wrappers for KP STZ constraints - intended for testing only
Zb <- function(b0,v,qc,p,w) {
b1 <- rep(0,p)
oo <- .C(C_Zb,b1=as.double(b1),as.double(b0),as.double(v),as.integer(qc),as.integer(p),as.double(w))
oo$b1
}
Ztb <- function(b0,v,qc,di,p,w) {
## p is length(b0)/di
w <- rep(0,2*p)
M <- v[1]
pp <- p
for (i in 1:M) pp <- pp/v[i+1];
p0 <- prod(v[1+1:M]-1)*pp
b1 <- rep(0,p0*di)
oo <- .C(C_Ztb,b1=as.double(b1),as.double(b0),as.double(v),as.integer(qc),as.integer(di),as.integer(p),as.double(w))
oo$b1
}
Try the mgcv package in your browser
Any scripts or data that you put into this service are public.
mgcv documentation built on May 29, 2024, 4:34 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Ztb Zb pdev neicov minres AddBVB isa mat.rowsum temp.seed eigen.approx lanczos treig pmmult pqr.qy pqr.R pqr block.reorder pRRt pcrossprod pforwardsolve pchol pbsi pqr2 pinv vcorr cholup choldrop dchol diagXVXd Xbd XWyd XWXd mgcv.omp trichol bandchol sdiag rmvn blas.thread.test dpnorm
Documented in bandchol blas.thread.test choldrop cholup diagXVXd dpnorm rmvn sdiag trichol Xbd XWXd XWyd
## (c) Simon N. Wood 2011-2023
## Many of the following are simple wrappers for C functions
dpnorm <- function(x0,x1) {
## Cancellation avoiding evaluation of pnorm(x1)-pnorm(x0)
## first avoid 1-1 problems by exchanging and changing sign of double +ve
ii <- x1>0&x0>0
d <- x0[ii];x0[ii] <- -x1[ii];x1[ii] <- -d
## now deal with points that are so close that cancellation error
## too large - might as well use density times interval width
ii <- abs(x1-x0) < sqrt(.Machine$double.eps)*dnorm((x1+x0)/2)
p <- x0; d <- x1[ii]-x0[ii]; m <- (x1[ii]+x0[ii])/2
p[ii] <- dnorm(m)*d
p[!ii] <- pnorm(x1[!ii]) - pnorm(x0[!ii])
p
} ## dpnorm
"%.%" <- function(a,b) {
if (inherits(a,"dgCMatrix")||inherits(b,"dgCMatrix"))
tensor.prod.model.matrix(list( ## following is coercion to double, general (no special structure), compressed column
as(as(as(a, "dMatrix"), "generalMatrix"), "CsparseMatrix"),
as(as(as(b, "dMatrix"), "generalMatrix"), "CsparseMatrix")
)) else tensor.prod.model.matrix(list(as.matrix(a),as.matrix(b)))
# tensor.prod.model.matrix(list(as(a,"dgCMatrix"),as(b,"dgCMatrix"))) else - deprecated
# tensor.prod.model.matrix(list(as.matrix(a),as.matrix(b)))
}
blas.thread.test <- function(n=1000,nt=4) {
inter <- interactive()
if (inter) prg <- txtProgressBar(min = 0, max = n, initial = 0,
char = "=",width = NA, title="Progress", style = 3)
for (i in 1:n) {
m <- sample(213:654,1)
p <- sample(3:153,1)
X <- matrix(runif(m*p),m,p)
er <- pqr(X,nt=nt)
rr <- range(pqr.qy(er,pqr.R(er))-X[,er$pivot])
if (rr[1] < -1e-4||rr[2] > 1e-4) {
break;
}
if (inter) setTxtProgressBar(prg, i)
}
if (inter) close(prg)
if (rr[1] < -1e-4||rr[2] > 1e-4) {
cat("BLAS thread safety problem at iteration",i,"\n")
} else cat("No problem encountered in",i,"iterations\n")
} ## blas.thread.test
rmvn <- function(n,mu,V) {
## generate multivariate normal deviates. e.g.
## V <- matrix(c(2,1,1,2),2,2); mu <- c(1,1);n <- 1000;z <- rmvn(n,mu,V);crossprod(sweep(z,2,colMeans(z)))/n
p <- ncol(V)
R <- mroot(V,rank=ncol(V)) ## RR' = V
if (is.matrix(mu)) {
if (ncol(mu)!=p||nrow(mu)!=n) stop("mu dimensions wrong")
z <- matrix(rnorm(p*n),n,p)%*%t(R) + mu
} else {
if (length(mu)!=p) stop("mu dimensions wrong")
z <- t(R%*% matrix(rnorm(p*n),p,n) + mu)
if (n==1) z <- as.numeric(z)
}
z
} ## rmvn
sdiag <- function(A,k=0) {
## extract sub or super diagonal of matrix (k=0 is leading)
p <- ncol(A)
n <- nrow(A)
if (k>p-1||-k > n-1) return()
if (k >= 0) {
i <- 1:n
j <- (k+1):p
} else {
i <- (-k+1):n
j <- 1:p
}
if (length(i)>length(j)) i <- i[1:length(j)] else j <- j[1:length(i)]
ii <- i + (j-1) * n
A[ii]
} ## sdiag
"sdiag<-" <- function(A,k=0,value) {
p <- ncol(A)
n <- nrow(A)
if (k>p-1||-k > n-1) return()
if (k >= 0) {
i <- 1:n
j <- (k+1):p
} else {
i <- (-k+1):n
j <- 1:p
}
if (length(i)>length(j)) i <- i[1:length(j)] else j <- j[1:length(i)]
ii <- i + (j-1) * n
A[ii] <- value
A
} ## "sdiag<-"
bandchol <- function(B) {
## obtain R such that R'R = A. Where A is banded matrix contained in R.
n <- ncol(B)
k <- 0
if (n==nrow(B)) { ## square matrix. Extract the diagonals
A <- B*0
for (i in 1:n) {
b <- sdiag(B,i-1)
if (sum(b!=0)!=0) {
k <- i ## largest index of a non-zero band
A[i,1:length(b)] <- b
}
}
B <- A[1:k,]
}
oo <- .C(C_band_chol,B=as.double(B),n=as.integer(n),k=as.integer(nrow(B)),info=as.integer(0))
if (oo$info<0) stop("something wrong with inputs to LAPACK routine")
if (oo$info>0) stop("not positive definite")
B <- matrix(oo$B,nrow(B),n)
if (k>0) { ## was square on entry, so also on exit...
A <- A * 0
for (i in 1:k) sdiag(A,i-1) <- B[i,1:(n-i+1)]
B <- A
}
B
} ## bandchol
trichol <- function(ld,sd) {
## obtain chol factor R of symm tridiag matrix, A, with leading diag
## ld and sub/super diags sd. R'R = A. On exit ld is diag of R and
## sd its super diagonal.
n <- length(ld)
if (n<2) stop("don't be silly")
if (n!=length(sd)+1) stop("sd should have exactly one less entry than ld")
oo <- .C(C_tri_chol,ld=as.double(ld),sd=as.double(sd),n=as.integer(n),info=as.integer(0))
if (oo$info<0) stop("something wrong with inputs to LAPACK routine")
if (oo$info>0) stop("not positive definite")
ld <- sqrt(oo$ld)
sd <- oo$sd*ld[1:(n-1)]
list(ld=ld,sd=sd)
}
mgcv.omp <- function() {
## does open MP appear to be available?
oo <- .C(C_mgcv_omp,a=as.integer(-1))
if (oo$a==1) TRUE else FALSE
}
## discretized covariate routines...
XWXd <- function(X,w,k,ks,ts,dt,v,qc,nthreads=1,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1,lt=NULL,rt=NULL) {
## Form X'WX given weights in w and X in compressed form in list X.
## each element of X is a (marginal) model submatrix. Full version
## is given by X[[i]][k[,i],] (see below for summation convention).
## list X relates to length(ts) separate
## terms. ith term starts at matrix ts[i] and has dt[i] marginal matrices.
## For summation convention, k[,ks[j,1]:ks[j,2]] gives index columns
## for matrix j, thereby allowing summation over matrix covariates....
## i.e. for q in ks[j,1]:ks[j,2] sum up X[[j]][k[,q],]
## Terms with several marginals are tensor products and may have
## constraints (if qc[i]>1), stored as a householder vector in v[[i]].
## check ts and k index start (assumed 1 here)
## if drop is non-NULL it contains index of rows/cols to drop from result
## * lt is array of terms to include in left matrix (assumed in ascending coef index order)
## * rt is array of terms to include in right matrix (assumed in ascending coef index order)
## * if both NULL all terms are included, if only one is NULL then used for left and right.
m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
n <- length(w);ptfull <- pt <- 0;
for (i in 1:nt) {
fullsize <- prod(p[ts[i]:(ts[i]+dt[i]-1)])
ptfull <- ptfull + fullsize
pt <- pt + fullsize - if (qc[i]>0) 1 else if (qc[i]<0) v[[i]][v[[i]][1]+2] else 0
}
if (inherits(X[[1]],"dgCMatrix")) { ## the marginals are sparse
if (length(ar.stop)>1||ar.stop!=-1) warning("AR not available with sparse marginals")
## create list for passing to C
if (any(qc<0)) stop("sparse method for Kronecker product contrasts not implemented")
m <- list(Xd=X,kd=k,ks=ks,v=v,ts=ts,dt=dt,qc=qc)
m$off <- attr(X,"off"); m$r <- attr(X,"r")
if (is.null(m$off)||is.null(m$r)) stop("reverse indices missing from sparse discrete marginals")
### code that could create the marginal reverse indices...
#for (j in 1:nrow(m$ks)) {
#nr <- nrow(m$Xd[[j]]) ## make sure we always tab to final stored row
#for (i in m$ks[j,1]:(m$ks[j,2]-1)) {
# m$r[,i] <- (1:length(m$kd[,i]))[order(m$kd[,i])]
# m$off[[i]] <- cumsum(c(1,tabulate(m$kd[,i],nbins=nr)))-1
#}
m$offstart <- cumsum(c(0,lapply(m$off,length)))
m$off <- unlist(m$off)
## Now C base all indices...
m$ks <- m$ks - 1; m$kd <- m$kd - 1; m$r <- m$r - 1; m$ts <- m$ts-1
nthreads <- as.integer(nthreads)
w <- as.double(w)
if ((!is.null(lt)||!is.null(rt))&&!is.null(drop)) {
lpip <- attr(X,"lpip") ## list of coefs for each term
rpi <- unlist(lpip[rt])
lpi <- unlist(lpip[lt])
if (is.null(lt)) lpi <- rpi
else if (is.null(rt)) rpi <- lpi
ldrop <- which(lpi %in% drop)
rdrop <- which(lpi %in% drop)
} else rdrop <- ldrop <- drop
if (is.null(lt)&&is.null(rt)) {
lt <- rt <- 1:nt
} else if (is.null(lt)) {
lt <- rt
} else if (is.null(rt)) rt <- lt;
lt <- as.integer(lt-1)
rt <- as.integer(rt-1)
XWX <- .Call(C_sXWXd,m,w,lt,rt,nthreads)
if (!is.null(drop)) {
Dl <- Diagonal(ncol(XWX),1)
XWX <- Dl[-ldrop,] %*% XWX %*% t(Dl[-rdrop,])
}
return(XWX) ## note that this is sparse
}
## block oriented code...
if (is.null(lt)&&is.null(lt)) {
# old .C code - can't handle long vector k
#oo <- .C(C_XWXd0,XWX =as.double(rep(0,ptfull^2)),X= as.double(unlist(X)),w=as.double(w),
# k=as.integer(k-1),ks=as.integer(ks-1),m=as.integer(m),p=as.integer(p), n=as.integer(n),
# ns=as.integer(nx), ts=as.integer(ts-1), as.integer(dt), nt=as.integer(nt),
# v = as.double(unlist(v)),qc=as.integer(qc),nthreads=as.integer(nthreads),
# ar.stop=as.integer(ar.stop-1),ar.weights=as.double(ar.w))
#XWX <- if (is.null(drop)) matrix(oo$XWX[1:pt^2],pt,pt) else matrix(oo$XWX[1:pt^2],pt,pt)[-drop,-drop]
XWX <- numeric(ptfull^2)
.Call(C_CXWXd0,XWX,as.double(unlist(X)),w,k-1L,as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(ts-1L), as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(nthreads),
as.integer(ar.stop-1L),as.double(ar.w))
XWX <- if (is.null(drop)) matrix(XWX[1:pt^2],pt,pt) else matrix(XWX[1:pt^2],pt,pt)[-drop,-drop]
} else {
lpip <- attr(X,"lpip") ## list of coefs for each term
rpi <- unlist(lpip[rt])
lpi <- unlist(lpip[lt])
if (is.null(lt)) {
lpi <- rpi
nrs <- lt <- 0
ncs <- length(rt)
} else {
nrs <- length(lt)
if (is.null(rt)) {
rt <- ncs <- 0
rpi <- lpi
} else ncs <- length(rt)
}
#oo <- .C(C_XWXd1,XWX =as.double(rep(0,ptfull^2)),X= as.double(unlist(X)),w=as.double(w),
# k=as.integer(k-1),ks=as.integer(ks-1),m=as.integer(m),p=as.integer(p), n=as.integer(n),
# ns=as.integer(nx), ts=as.integer(ts-1), dt=as.integer(dt), nt=as.integer(nt),
# v = as.double(unlist(v)),qc=as.integer(qc),nthreads=as.integer(nthreads),
# ar.stop=as.integer(ar.stop-1),ar.weights=as.double(ar.w),rs=as.integer(lt-1),
# cs=as.integer(rt-1),nrs=as.integer(nrs),ncs=as.integer(ncs))#)
#XWX <- matrix(oo$XWX[1:(length(lpi)*length(rpi))],length(lpi),length(rpi))
XWX <- numeric(ptfull^2)
.Call(C_CXWXd1,XWX,as.double(unlist(X)),w,k-1L,as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(ts-1L), as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(nthreads),
as.integer(ar.stop-1L),as.double(ar.w),as.integer(lt-1L),as.integer(rt-1L))
XWX <- matrix(XWX[1:(length(lpi)*length(rpi))],length(lpi),length(rpi))
if (!is.null(drop)) {
ldrop <- which(lpi %in% drop)
rdrop <- which(lpi %in% drop)
if (length(ldrop)>0||length(rdrop)>0) XWX <-
if (length(ldrop==0)) XWX[,-rdrop] else if (length(rdrop)==0) XWX[-ldrop,] else XWX[-ldrop,-rdrop]
}
}
XWX
} ## XWXd
XWyd <- function(X,w,y,k,ks,ts,dt,v,qc,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1,lt=NULL) {
## X'Wy...
## if lt if not NULL then it lists the discrete terms to include (from X)
## returned vector/matrix only includes rows for selected terms
m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
n <- length(w);
if (is.null(lt)) {
pt <- 0
for (i in 1:nt) pt <- pt + prod(p[ts[i]:(ts[i]+dt[i]-1)]) - if (qc[i]>0) 1 else if (qc[i]<0) v[[i]][v[[i]][1]+2] else 0
lt <- 1:nt
} else {
lpip <- attr(X,"lpip") ## list of coefs for each term
lpi <- unlist(lpip[lt]) ## coefs corresponding to terms selected by lt
if (!is.null(drop)) drop <- which(lpi %in% drop) ## rebase drop
pt <- length(lpi)
}
cy <- if (is.matrix(y)) ncol(y) else 1
if (inherits(X[[1]],"dgCMatrix")) { ## the marginals are sparse
## create list for passing to C
m <- list(Xd=X,kd=k,ks=ks,v=v,ts=ts,dt=dt,qc=qc)
m$off <- attr(X,"off"); m$r <- attr(X,"r")
if (is.null(m$off)||is.null(m$r)) stop("reverse indices missing from sparse discrete marginals")
m$offstart <- cumsum(c(0,lapply(m$off,length)))
m$off <- unlist(m$off)
## Now C base all indices...
m$ks <- m$ks - 1; m$kd <- m$kd - 1; m$r <- m$r - 1; m$ts <- m$ts-1
Wy <- as.double(w*y);lt <- as.integer(lt-1)
XWy <- .Call(C_sXyd,m,Wy,lt)
if (cy>1) XWy <- matrix(XWy,ncol=cy)
if (!is.null(drop)) XWy <- if (cy>1) XWy[-drop,] else XWy[-drop]
} else { ## dense marginals case
## old .C code - can't handle long vector k
#oo <- .C(C_XWyd,XWy=rep(0,pt*cy),y=as.double(y),X=as.double(unlist(X)),w=as.double(w),k=as.integer(k-1),
# ks=as.integer(ks-1),
# m=as.integer(m),p=as.integer(p),n=as.integer(n),cy=as.integer(cy), nx=as.integer(nx), ts=as.integer(ts-1),
# dt=as.integer(dt),nt=as.integer(nt),v=as.double(unlist(v)),qc=as.integer(qc),
# ar.stop=as.integer(ar.stop-1),ar.row=as.integer(ar.row-1),ar.weights=as.double(ar.w),
# cs=as.integer(lt-1),ncs=as.integer(length(lt)))
#if (cy>1) XWy <- if (is.null(drop)) matrix(oo$XWy,pt,cy) else matrix(oo$XWy,pt,cy)[-drop,] else
#XWy <- if (is.null(drop)) oo$XWy else oo$XWy[-drop]
XWy <- numeric(pt*cy)
.Call(C_CXWyd,XWy,as.double(y),as.double(unlist(X)),as.double(w),k-1L,as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(cy),as.integer(ts-1L), as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(ar.stop-1L),
as.integer(ar.row-1L),as.double(ar.w),as.integer(lt-1L))
if (cy>1) { XWy <- if (is.null(drop)) matrix(XWy,pt,cy) else matrix(XWy,pt,cy)[-drop,] } else {
XWy <- if (is.null(drop)) XWy else XWy[-drop]
}
}
XWy
} ## XWyd
Xbd <- function(X,beta,k,ks,ts,dt,v,qc,drop=NULL,lt=NULL) {
## note that drop may contain the index of columns of X to drop before multiplying by beta.
## equivalently we can insert zero elements into beta in the appropriate places.
## if lt if not NULL then it lists the discrete terms to include (from X)
n <- if (is.matrix(k)) nrow(k) else length(k) ## number of data
m <- unlist(lapply(X,nrow)) ## number of rows in each discrete model matrix
p <- unlist(lapply(X,ncol)) ## number of cols in each discrete model matrix
nx <- length(X) ## number of model matrices
if (length(p)!=nx) stop("something wrong with matrix list - not all matrices?")
nt <- length(ts) ## number of terms
if (!is.null(drop)) {
b <- if (is.matrix(beta)) matrix(0,nrow(beta)+length(drop),ncol(beta)) else rep(0,length(beta)+length(drop))
if (is.matrix(beta)) b[-drop,] <- beta else b[-drop] <- beta
beta <- b
}
if (is.null(lt)) lt <- 1:nt
bc <- if (is.matrix(beta)) ncol(beta) else 1 ## number of columns in beta
## The C code mechanism for dealing with lt is very basic, and requires that beta is re-ordered and
## truncated to relate only to the selected terms, in the order they are selected.
lpip <- attr(X,"lpip")
if (!is.null(lpip)) { ## then X list may not be in coef order...
lpip <- unlist(lpip[lt])
beta <- if (is.matrix(beta)) beta[lpip,] else beta[lpip] ## select params required in correct order
}
if (inherits(X[[1]],"dgCMatrix")) { ## the marginals are sparse
## create list for passing to C
m <- list(Xd=X,kd=k,ks=ks,v=v,ts=ts,dt=dt,qc=qc)
m$off <- attr(X,"off"); m$r <- attr(X,"r")
if (is.null(m$off)||is.null(m$r)) stop("reverse indices missing from sparse discrete marginals")
m$offstart <- cumsum(c(0,lapply(m$off,length)))
m$off <- unlist(m$off)
## Now C base all indices...
m$ks <- m$ks - 1; m$kd <- m$kd - 1; m$r <- m$r - 1; m$ts <- m$ts-1
beta <- as.matrix(beta);storage.mode(beta) <- "double"
lt <- as.integer(lt-1)
Xb <- .Call(C_sXbd,m,beta,lt)
if (bc>1) Xb <- matrix(Xb,ncol=bc)
} else { ## dense marginals case
#oo <- .C(C_Xbd,f=as.double(rep(0,n*bc)),beta=as.double(beta),X=as.double(unlist(X)),k=as.integer(k-1),
# ks = as.integer(ks-1),
# m=as.integer(m),p=as.integer(p), n=as.integer(n), nx=as.integer(nx), ts=as.integer(ts-1),
# as.integer(dt), as.integer(nt),as.double(unlist(v)),as.integer(qc),as.integer(bc),as.integer(lt-1),as.integer(length(lt)))
#Xb <- if (is.matrix(beta)) matrix(oo$f,n,bc) else oo$f
f <- numeric(n*bc)
.Call(C_CXbd,f,beta,as.double(unlist(X)),k-1L, as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(ts-1L),as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(bc), as.integer(lt-1L))
Xb <- if (is.matrix(beta)) matrix(f,n,bc) else f
}
return(Xb)
} ## Xbd
diagXVXd <- function(X,V,k,ks,ts,dt,v,qc,drop=NULL,nthreads=1,lt=NULL,rt=NULL) {
## discrete computation of diag(XVX')
n <- if (is.matrix(k)) nrow(k) else length(k)
m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
if (is.null(lt)&&is.null(rt)) rt <- lt <- 1:nt else {
if (is.null(rt)) rt <- lt else if (is.null(lt)) lt <- rt
}
if (is.null(rt)) rt <- 1:nt
if (inherits(X[[1]],"dgCMatrix")) { ## the marginals are sparse
## create list for passing to C
m <- list(Xd=X,kd=k,ks=ks,v=v,ts=ts,dt=dt,qc=qc)
m$off <- attr(X,"off"); m$r <- attr(X,"r")
if (is.null(m$off)||is.null(m$r)) stop("reverse indices missing from sparse discrete marginals")
m$offstart <- cumsum(c(0,lapply(m$off,length)))
m$off <- unlist(m$off)
## Now C base all indices...
m$ks <- m$ks - 1; m$kd <- m$kd - 1; m$r <- m$r - 1; m$ts <- m$ts-1
if (!is.null(drop)) {
D <- Diagonal(ncol(V)+length(drop),1)[-drop,]
V <- t(D) %*% V %*% D
}
## The C code mechanism for dealing with rt and lt is very basic, and requires that V is
## re-ordered and truncated to relate only to the selected terms, in the order they are selected.
lpip <- attr(X,"lpip")
if (!is.null(lpip)) { ## then X list may not be in coef order...
lpi <- unlist(lpip[lt])
rpi <- unlist(lpip[rt])
V <- V[lpi,rpi,drop=FALSE] ## select part of V required in correct order
}
lt <- as.integer(lt-1);rt <- as.integer(rt-1)
D <- .Call(C_sdiagXVXt,m , V, lt, rt)
} else { ## dense marginals
if (!is.null(drop)) {
pv <- ncol(V)+length(drop)
V0 <- matrix(0,pv,pv)
V0[-drop,-drop] <- V
V <- V0;rm(V0)
} else pv <- ncol(V)
## The C code mechanism for dealing with rt and lt is very basic, and requires that V is
## re-ordered and truncated to relate only to the selected terms, in the order they are selected.
lpip <- attr(X,"lpip")
if (!is.null(lpip)) { ## then X list may not be in coef order...
lpi <- unlist(lpip[lt])
rpi <- unlist(lpip[rt])
V <- V[lpi,rpi,drop=FALSE] ## select part of V required in correct order
}
# oo <- .C(C_diagXVXt,diag=as.double(rep(0,n)),V=as.double(V),X=as.double(unlist(X)),k=as.integer(k-1),
# ks=as.integer(ks-1),m=as.integer(m),p=as.integer(p), n=as.integer(n), nx=as.integer(nx),
# ts=as.integer(ts-1), as.integer(dt), as.integer(nt),as.double(unlist(v)),as.integer(qc),as.integer(nrow(V)),
# as.integer(ncol(V)),as.integer(nthreads),as.integer(lt-1),as.integer(length(lt)),as.integer(rt-1),as.integer(length(rt)))
# D <- oo$diag
D <- numeric(n)
.Call(C_CdiagXVXt,D,V,as.double(unlist(X)),k-1L,as.integer(ks-1L),as.integer(m),as.integer(p),
as.integer(ts-1L), as.integer(dt),as.double(unlist(v)),as.integer(qc),as.integer(nthreads),
as.integer(lt-1L),as.integer(rt-1L))
}
D
} ## diagXVXd
dchol <- function(dA,R) {
## if dA contains matrix dA/dx where R is chol factor s.t. R'R = A
## then this routine returns dR/dx...
p <- ncol(R)
oo <- .C(C_dchol,dA=as.double(dA),R=as.double(R),dR=as.double(R*0),p=as.integer(ncol(R)))
return(matrix(oo$dR,p,p))
} ## dchol
choldrop <- function(R,k) {
## routine to update Cholesky factor R of A on dropping row/col k of A.
## R can be upper triangular, in which case (R'R=A) or lower triangular in
## which case RR'=A...
n <- as.integer(ncol(R))
k1 <- as.integer(k-1)
ut <- as.integer(as.numeric(R[1,2]!=0))
if (k<1||k>n) return(R)
Rup <- matrix(0,n-1,n-1)
#oo <- .C(C_chol_down,R=as.double(R),Rup=as.double(Rup),n=as.integer(n),k=as.integer(k-1),ut=as.integer(ut))
.Call(C_mgcv_chol_down,R,Rup,n,k1,ut)
#matrix(oo$Rup,n-1,n-1)
Rup
} ## choldrop
cholup <- function(R,u,up=TRUE) {
## routine to update Cholesky factor R to the factor of R'R + uu' (up == TRUE)
## or R'R - uu' (up=FALSE).
n <- as.integer(ncol(R))
up <- as.integer(up)
eps <- as.double(.Machine$double.eps)
R1 <- R * 1.0
.Call(C_mgcv_chol_up,R1,u,n,up,eps)
if (up==0) if ((n>1 && R1[2,1] < -1)||(n==1&&u[1]>R[1])) stop("update not positive definite")
R1
} ## cholup
vcorr <- function(dR,Vr,trans=TRUE) {
## Suppose b = sum_k op(dR[[k]])%*%z*r_k, z ~ N(0,Ip), r ~ N(0,Vr). vcorr returns cov(b).
## dR is a list of p by p matrices. 'op' is 't' if trans=TRUE and I() otherwise.
p <- ncol(dR[[1]])
M <- if (trans) ncol(Vr) else -ncol(Vr) ## sign signals transpose or not to C code
if (abs(M)!=length(dR)) stop("internal error in vcorr, please report to simon.wood@r-project.org")
oo <- .C(C_vcorr,dR=as.double(unlist(dR)),Vr=as.double(Vr),Vb=as.double(rep(0,p*p)),
p=as.integer(p),M=as.integer(M))
return(matrix(oo$Vb,p,p))
} ## vcorr
pinv <- function(X,svd=FALSE) {
## a pseudoinverse for n by p, n>p matrices
qrx <- qr(X,tol=0,LAPACK=TRUE)
R <- qr.R(qrx);Q <- qr.Q(qrx)
rr <- Rrank(R)
if (svd&&rr<ncol(R)) {
piv <- 1:ncol(X); piv[qrx$pivot] <- 1:ncol(X)
er <- svd(R[,piv])
d <- er$d*0;d[1:rr] <- 1/er$d[1:rr]
X <- Q%*%er$u%*%(d*t(er$v))
} else {
Ri <- R*0
Ri[1:rr,1:rr] <- backsolve(R[1:rr,1:rr],diag(rr))
X[,qrx$pivot] <- Q%*%t(Ri)
}
X
} ## end pinv
pqr2 <- function(x,nt=1,nb=30) {
## Function for parallel pivoted qr decomposition of a matrix using LAPACK
## householder routines. Currently uses a block algorithm.
## library(mgcv); n <- 4000;p<-3000;x <- matrix(runif(n*p),n,p)
## system.time(qrx <- qr(x,LAPACK=TRUE))
## system.time(qrx2 <- mgcv:::pqr2(x,2))
## system.time(qrx3 <- mgcv:::pqr(x,2))
## range(qrx2$qr-qrx$qr)
p <- ncol(x)
beta <- rep(0.0,p)
piv <- as.integer(rep(0,p))
## need to force a copy of x, otherwise x will be over-written
## by .Call *in environment from which function is called*
x <- x*1
rank <- .Call(C_mgcv_Rpiqr,x,beta,piv,nt,nb)
ret <- list(qr=x,rank=rank,qraux=beta,pivot=piv+1)
attr(ret,"useLAPACK") <- TRUE
class(ret) <- "qr"
ret
} ## pqr2
pbsi <- function(R,nt=1,copy=TRUE) {
## parallel back substitution inversion of upper triangular R
## library(mgcv); n <- 500;p<-400;x <- matrix(runif(n*p),n,p)
## qrx <- qr(x);R <- qr.R(qrx)
## system.time(Ri <- mgcv:::pbsi(R,2))
## system.time(Ri2 <- backsolve(R,diag(p)));range(Ri-Ri2)
if (copy) R <- R * 1 ## ensure that R modified only within pbsi
.Call(C_mgcv_Rpbsi,R,nt)
R
} ## pbsi
pchol <- function(A,nt=1,nb=40) {
## parallel Choleski factorization.
## library(mgcv);
## set.seed(2);n <- 200;r <- 190;A <- tcrossprod(matrix(runif(n*r),n,r))
## system.time(R <- chol(A,pivot=TRUE));system.time(L <- mgcv:::pchol(A));range(R[1:r,]-L[1:r,])
## k <- 30;range(R[1:k,1:k]-L[1:k,1:k])
## system.time(L <- mgcv:::pchol(A,nt=2,nb=30))
## piv <- attr(L,"pivot");attr(L,"rank");range(crossprod(L)-A[piv,piv])
## should nb be obtained from 'ILAENV' as page 23 of Lucas 2004??
piv <- as.integer(rep(0,ncol(A)))
A <- A*1 ## otherwise over-write in calling env!
rank <- .Call(C_mgcv_Rpchol,A,piv,nt,nb)
attr(A,"pivot") <- piv+1;attr(A,"rank") <- rank
A
}
pforwardsolve <- function(R,B,nt=1) {
## parallel forward solve via simple col splitting...
if (!is.matrix(B)) B <- as.matrix(B)
.Call(C_mgcv_Rpforwardsolve,R,B,nt)
}
pcrossprod <- function(A,trans=FALSE,nt=1,nb=30) {
## parallel cross prod A'A or AA' if trans==TRUE...
if (!is.matrix(A)) A <- as.matrix(A)
if (trans) A <- t(A)
.Call(C_mgcv_Rpcross,A,nt,nb)
}
pRRt <- function(R,nt=1) {
## parallel RR' for upper triangular R
## following creates index of lower triangular elements...
## n <- 4000;a <- rep(1:n,n);b <- rep(1:n,each=n);which(a>=b) -> ii;a[ii]+(b[ii]-1)*n->ii ## lower
## n <- 4000;a <- rep(1:n,n);b <- rep(1:n,each=n);which(a<=b) -> ii;a[ii]+(b[ii]-1)*n->ii ## upper
## library(mgcv);R <- matrix(0,n,n);R[ii] <- runif(n*(n+1)/2)
## Note: A[a-b<=0] <- 0 zeroes upper triangle
## system.time(A <- mgcv:::pRRt(R,2))
## system.time(A2 <- tcrossprod(R));range(A-A2)
n <- nrow(R)
A <- matrix(0,n,n)
.Call(C_mgcv_RPPt,A,R,nt)
A
}
block.reorder <- function(x,n.blocks=1,reverse=FALSE) {
## takes a matrix x divides it into n.blocks row-wise blocks, and re-orders
## so that the blocks are stored one after the other.
## e.g. library(mgcv); x <- matrix(1:18,6,3);xb <- mgcv:::block.reorder(x,2)
## x;xb;mgcv:::block.reorder(xb,2,TRUE)
r = nrow(x);cols = ncol(x);
if (n.blocks <= 1) return(x);
if (r%%n.blocks) {
nb = ceiling(r/n.blocks)
} else nb = r/n.blocks;
oo <- .C(C_row_block_reorder,x=as.double(x),as.integer(r),as.integer(cols),
as.integer(nb),as.integer(reverse));
matrix(oo$x,r,cols)
} ## block.reorder
pqr <- function(x,nt=1) {
## parallel QR decomposition, using openMP in C, and up to nt threads (only if worthwhile)
## library(mgcv);n <- 20;p<-4;X <- matrix(runif(n*p),n,p);er <- mgcv:::pqr(X,nt=2)
## range(mgcv:::pqr.qy(er,mgcv:::pqr.R(er))-X[,er$pivot])
x.c <- ncol(x);r <- nrow(x)
oo <- .C(C_mgcv_pqr,x=as.double(c(x,rep(0,nt*x.c^2))),as.integer(r),as.integer(x.c),
pivot=as.integer(rep(0,x.c)), tau=as.double(rep(0,(nt+1)*x.c)),as.integer(nt))
list(x=oo$x,r=r,c=x.c,tau=oo$tau,pivot=oo$pivot+1,nt=nt)
}
pqr.R <- function(x) {
## x is an object returned by pqr. This extracts the R factor...
## e.g. as pqr then...
## R <- mgcv:::pqr.R(er); R0 <- qr.R(qr(X,tol=0))
## svd(R)$d;svd(R0)$d
oo <- .C(C_getRpqr,R=as.double(rep(0,x$c^2)),as.double(x$x),as.integer(x$r),as.integer(x$c),
as.integer(x$c),as.integer(x$nt))
matrix(oo$R,x$c,x$c)
}
pqr.qy <- function(x,a,tr=FALSE) {
## x contains a parallel QR decomp as computed by pqr. a is a matrix. computes
## Qa or Q'a depending on tr.
## e.g. as above, then...
## a <- diag(p);Q <- mgcv:::pqr.qy(er,a);crossprod(Q)
## X[,er$pivot+1];Q%*%R
## Qt <- mgcv:::pqr.qy(er,diag(n),TRUE);Qt%*%t(Qt);range(Q-t(Qt))
## Q <- qr.Q(qr(X,tol=0));z <- runif(n);y0<-t(Q)%*%z
## mgcv:::pqr.qy(er,z,TRUE)->y
## z <- runif(p);y0<-Q%*%z;mgcv:::pqr.qy(er,z)->y
if (is.matrix(a)) a.c <- ncol(a) else a.c <- 1
if (tr) {
if (is.matrix(a)) { if (nrow(a) != x$r) stop("a has wrong number of rows") }
else if (length(a) != x$r) stop("a has wrong number of rows")
} else {
if (is.matrix(a)) { if (nrow(a) != x$c) stop("a has wrong number of rows") }
else if (length(a) != x$c) stop("a has wrong number of rows")
a <- c(a,rep(0,a.c*(x$r-x$c)))
}
oo <- .C(C_mgcv_pqrqy,a=as.double(a),as.double(x$x),as.double(x$tau),as.integer(x$r),
as.integer(x$c),as.integer(a.c),as.integer(tr),as.integer(x$nt))
if (tr) return(matrix(oo$a[1:(a.c*x$c)],x$c,a.c)) else
return(matrix(oo$a,x$r,a.c))
}
pmmult <- function(A,B,tA=FALSE,tB=FALSE,nt=1) {
## parallel matrix multiplication (not for use on vectors or thin matrices)
## library(mgcv);r <- 10;c <- 5;n <- 8
## A <- matrix(runif(r*n),r,n);B <- matrix(runif(n*c),n,c);range(A%*%B-mgcv:::pmmult(A,B,nt=1))
## A <- matrix(runif(r*n),n,r);B <- matrix(runif(n*c),n,c);range(t(A)%*%B-mgcv:::pmmult(A,B,TRUE,FALSE,nt=1))
## A <- matrix(runif(r*n),n,r);B <- matrix(runif(n*c),c,n);range(t(A)%*%t(B)-mgcv:::pmmult(A,B,TRUE,TRUE,nt=1))
## A <- matrix(runif(r*n),r,n);B <- matrix(runif(n*c),c,n);range(A%*%t(B)-mgcv:::pmmult(A,B,FALSE,TRUE,nt=1))
if (tA) { n = nrow(A);r = ncol(A)} else {n = ncol(A);r = nrow(A)}
if (tB) { c = nrow(B)} else {c = ncol(B)}
C <- rep(0,r * c)
oo <- .C(C_mgcv_pmmult,C=as.double(C),as.double(A),as.double(B),as.integer(tA),as.integer(tB),as.integer(r),
as.integer(c),as.integer(n),as.integer(nt));
matrix(oo$C,r,c)
}
treig <- function(ld,sd,vec=FALSE,descend=FALSE) {
## eigen decomposition of tri-diagonal matrix with leading diagonal ld and
## sub-diagonal sd...
n <- length(ld)
v <- if (vec) rep(0,n*n) else 0
oo <- .C(C_mgcv_trisymeig,d=as.double(ld),g=as.double(sd),v=as.double(v),n=as.integer(n),
get.vec=as.integer(vec),descending=as.integer(descend))
v <- if (vec) matrix(oo$v,n,n) else NA
list(values=oo$d,vectors=v)
} ## treig
lanczos <- function(A,v0,M,Av = function(A,v) A%*%v,n=ncol(A)) {
## Apply M steps of Lanczos starting at v0 for n by n +ve semi definite matrix A.
## Av is a function forming the product of matrix A with vector v.
## A can be a matrix, in which case the default Av applies, or
## it could be a list of arguments used to define the multiplication
## in some other way (e.g. as a sequence of matrix products) defined
## in a custom Av...
## The function is used to find an approximate eigen value CDF for A
## as described in Lin, Saad and Yang (2016) SIAM Review 58(1), 34-65
## section 3.2.1 in particular.
v0.norm <- sqrt(sum(v0^2))
gamma <- epsilon <- rep(0,M)
q <- matrix(v0/v0.norm,n,M)
for (j in 1:M) {
c <- Av(A,q[,j])
if (j==1) vAv <- sum(v0*c)*v0.norm ## v0'Av0 - useful for estimating tr(A)
gamma[j] <- sum(c*q[,j])
c <- c - gamma[j]*q[,j]
if (j>1) {
c <- c - epsilon[j-1]*q[,j-1]
cq <- drop(t(c) %*% q[,1:j,drop=FALSE])
c <- c - colSums(cq*t(q[,1:j]))
cq <- drop(t(c) %*% q[,1:j,drop=FALSE])
c <- c - colSums(cq*t(q[,1:j]))
}
epsilon[j] <- sqrt(sum(c^2))
if (j<M) q[,j+1] <- c/epsilon[j]
}
et <- treig(gamma,epsilon,descend=FALSE,vec=TRUE)
## compute error bounds on the eigenvalues of
## A using the method described in section 3.2 of
## Parlett, BN (1998) The Symmetric Eigenvalue Problem, SIAM
err <- abs(et$vectors[M,])*epsilon[M]
theta <- c(0,et$values)
tau <- et$vectors[1,]
eta <- c(0,cumsum(tau^2))
## theta is vector of eigenvalues at which CDF jumps
## tau^2 is jump size, eta is CDF at theta
## lam.ub is upper bound on larget eigenvalue
list(theta=theta,tau=tau,eta=eta,err=err,vAv=vAv)
} ## lanczos
eigen.approx <- function(A,Av = function(A,v) A%*%v,M=20,n.rep=20,n=ncol(A),seed=1) {
## get the approximate eigenvalues of n by n +ve semi def matrix A. Av is
## the function for multiplying a vector by the matrix defined by A. If A
## is simply a matrix then the default Av is sufficient. M is the number of
## Lanczos steps to use, and n.rep the number of random replicates to average
## over. Routine restores RNG to pre-call state on exit.
## Based on Lin, Saad and Yang (2016) SIAM Review 58(1), 34-65
## section 3.2.1, but extended to only use this approximation for the
## eigen-values that have yet to converge.
## The CDF approximation is based on their Appendix C proposal, rather than
## using a Gausssian kernel approximation to the pdf and cdf and then
## inverting by tabulation. This is because the latter tends to oversmooth the
## CDF in a way that is unhelpful for rank deficient matrices.
## The Gaussian kernel approach would probably be prefereable for full rank matrices,
## since it then benefits from the extra stability of kernel smoothing.
if (is.finite(seed)) a <- temp.seed(seed) ## seed RNG and store state
eva <- rep(0,n)
trA <- rep(0,n.rep)
tol <- .Machine$double.eps^.5
for (r in 1:n.rep) {
v0 <- rnorm(n)
lz <- lanczos(A,v0,M=M,Av=Av,n=n)
trA[r] <- lz$vAv
## following is suggested in Appendix C of LSY, and is quite important
## to avoid slight downward bias...
eta1 <- c(0,(lz$eta[1:M]*0.5+0.5*lz$eta[1:M+1]))
eta1[2] <- lz$eta[2] ## correction to avoid over-estimation in lower tail if rank def
conv <- lz$err<lz$theta[M+1]*tol ## these eigenvalues are converged
upper.uconv <- if (any(!conv)) max(which(!conv)) else 0 ## last uncoverged
n.conv <- M - upper.uconv ## number converged
lz$theta[lz$theta<0] <- 0
theta.conv <- if (n.conv) lz$theta[(upper.uconv+1):M+1] else rep(0,0)
if (upper.uconv) {
eta <- eta1[1:(upper.uconv+1)]
theta <- lz$theta[1:(upper.uconv+1)]
eta <- eta/max(eta)
nri <- c(diff(eta)!=0,TRUE) ## strip out duplicates
eva <- eva + c(approx(eta[nri],theta[nri],seq(0,1,length=n-n.conv),method="linear",rule=2)$y,theta.conv)
} else {
eva <- eva + c(rep(0,n-n.conv),theta.conv)
}
}
if (is.finite(seed)) temp.seed(a) ## restore RNG state
eva <- eva/n.rep
trA.sd <- sd(trA)/sqrt(n.rep);trA <- mean(trA)
if (abs(sum(eva)-trA)>2.5*trA.sd) { ## evidence for bias in eigen-spectrum
eva <- eva*trA/sum(eva) ## correction
}
eva
} ## eigen.approx
temp.seed <- function(x) {
## when called with a numeric x stores the state of the RNG and sets its seed to
## x. Returns an object of class "rng.state". When called with an object of this
## class created on a previus call to this function, resets the RNG to the state
## on entry to that previous call.
if (inherits(x,"rng.state")) {
RNGkind(x$kind[1],x$kind[2])
assign(".Random.seed",x$seed,envir=.GlobalEnv) ## RNG behaves as if it had not been used
} else {
seed <- try(get(".Random.seed",envir=.GlobalEnv),silent=TRUE) ## store RNG seed
if (inherits(seed,"try-error")) {
runif(1)
seed <- get(".Random.seed",envir=.GlobalEnv)
}
kind <- RNGkind(NULL)
RNGkind("default","default")
set.seed(x) ## ensure repeatability
x <- list(seed=seed,kind=kind)
class(x) <- "rng.state"
return(x)
}
} ## temp.seed
mat.rowsum <- function(X,m,k) {
## Let X be n by p and m of length M. Produces an M by p matrix B
## where the ith row of B is the sum of the rows X[k[j],] where
## j = (m[i-1]+1):m[i]. m[0] is taken as 0.
## n <- 10;p <- 5;X <- matrix(runif(n*p),n,p)
## m <- c(3,5,8,11);k <- c(1,4,3,6,1,5,7,10,9,5,6)
## mgcv:::mat.rowsum(X,m,k)
if (max(k)>nrow(X)||min(k)<1) stop("index vector has invalid entries")
k <- k - 1 ## R to C index conversion
.Call(C_mrow_sum,X,m,k)
} ## mat.rowsum
isa <- function(R,nt=1) {
## Finds the elements of (R'R)^{-1} on NZP(R+R').
if (!inherits(R,c("dgCMatrix","dtCMatrix"))) stop("isa requires a dg/tCMatrix")
nt <- round(nt)
if (nt<1) nt = 1
Hpi <- R + t(R)
.Call(C_isa1p,t(R),Hpi,nt)
Hpi
} ## isa
AddBVB <- function(A,Bt,VBt) {
## Add B %*% V %*% t(B) to calss 'dgCMatrix' A returning result on NZP(A) only
## (i.e. discarding elements of BVB' not in NZP(A)), Bt is the transpose
## of B. B and VBt are class 'matrix'
A@x <- A@x * 1.0 ## force copy, otherwise A and return value modified
.Call(C_AddBVB,A,Bt,VBt)
A
} ## AddBVB
minres <- function(R,u,b) {
## routine to solve (R'R-uu')x = b using minres algorithm.
## set.seed(0);n <- 100;p <- 20;X <- matrix(runif(n*p)-.5,n,p);R <- chol(crossprod(X));b <- runif(p);k <- 1;
## solve(crossprod(X[-k,]),b);mgcv:::minres(R,t(X[k,]),b)
x <- b; p <- length(b);
m <- if (is.matrix(u)) ncol(u) else 1
work <- rep(0,p*(m+7)+m)
oo <- .C(C_minres,R=as.double(R), u=as.double(u),b=as.double(b), x=as.double(x), p=as.integer(p),m=as.integer(m),work=as.double(work))
cat("\n niter : ",oo$m,"\n")
oo$x
}
neicov <- function(Dd,nei) {
## wrapper for nei_cov. Dd is n by p matrix of leave one out perturbations to
## coef vectors. nei is neighbourhood structure.
p <- ncol(Dd)
V <- matrix(0,p,p)
k <- nei$k-1
.Call(C_nei_cov,V,Dd,nei$m,k)
(V+t(V))/2
} ## neicov
## Routines to force matrix to pdef
pdev <- function(A) {
## Force A to meet necessary conditions for +ve def from Thm 4.2.8. of Golub and van Load 4th ed.
## Non-positive diagonals are set to diagonal dominance. Off diagonals are then set to just meet
## necessary conditions.
## NOTE: A is modified directly in situ, so care needed to force a copy to be kept if it's needed
## A0 <- A will not work, as R only actually copies when R modifies!
## A0 <- A; A0[1,1] <- A0[1,1] + 1;A0[1,1] <- A0[1,1] - 1 works.
if (inherits(A,"Matrix")) {
if (!inherits(A,"dgCMatrix")) A <- as(as(as(A, "dMatrix"), "generalMatrix"), "CsparseMatrix")
da <- diag(A); ii <- which(da==0)
if (length(ii)) diag(A)[ii] <- -1e-30 ## Avoid C code having to insert extra non-zeroes
mod <- .Call(C_spdev,A)
} else { ## dense matrix
mod <- .Call(C_dpdev,A)
}
if (mod>0) attr(A,"modified") <- mod
return(A)
}
## following are wrappers for KP STZ constraints - intended for testing only
Zb <- function(b0,v,qc,p,w) {
b1 <- rep(0,p)
oo <- .C(C_Zb,b1=as.double(b1),as.double(b0),as.double(v),as.integer(qc),as.integer(p),as.double(w))
oo$b1
}
Ztb <- function(b0,v,qc,di,p,w) {
## p is length(b0)/di
w <- rep(0,2*p)
M <- v[1]
pp <- p
for (i in 1:M) pp <- pp/v[i+1];
p0 <- prod(v[1+1:M]-1)*pp
b1 <- rep(0,p0*di)
oo <- .C(C_Ztb,b1=as.double(b1),as.double(b0),as.double(v),as.integer(qc),as.integer(di),as.integer(p),as.double(w))
oo$b1
}
Try the mgcv package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.