Commit 1b84d0b8 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-20

parent e0667888
...@@ -4,6 +4,25 @@ ...@@ -4,6 +4,25 @@
Currently deprecated and liable to be removed: Currently deprecated and liable to be removed:
- gam performance iteration (1.8-19, Sept 2017) - gam performance iteration (1.8-19, Sept 2017)
1.8-20
* bam(,discrete=TRUE) could produce garbage with ti(x,z,k=c(6,5),mc=c(T,F))
because tensor re-ordering for efficiency failed to re-order mc (this is
a *very* specialist bug!). Thanks to Fabian Scheipl.
* plot(...,residuals=TRUE) weighted the working residuals by the sqrt working
weights divided by the mean sqrt working weight. The standardization by the
mean sqrt weight was non standard and has been removed.
* Fix to bad bug in bam(...,discrete=TRUE) offset handling, and predict.bamd
modified to avoid failure predicting with offset. Thanks to Paul Shearer.
* fix of typo in bgam.fit, which caused failure of extended families
when dataset larger than chunk size. Thanks Martijn Wieling.
* bam(...,discrete=TRUE)/bgam.fitd modified to use fisher weights with
extended families if rho!=0.
1.8-19 1.8-19
** bam() now accepts extended families (i.e. nb, tw, ocat etc) ** bam() now accepts extended families (i.e. nb, tw, ocat etc)
......
Package: mgcv Package: mgcv
Version: 1.8-19 Version: 1.8-20
Author: Simon Wood <simon.wood@r-project.org> Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org> Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with Automatic Smoothness Title: Mixed GAM Computation Vehicle with Automatic Smoothness
...@@ -18,6 +18,6 @@ LazyLoad: yes ...@@ -18,6 +18,6 @@ LazyLoad: yes
ByteCompile: yes ByteCompile: yes
License: GPL (>= 2) License: GPL (>= 2)
NeedsCompilation: yes NeedsCompilation: yes
Packaged: 2017-08-29 08:25:01 UTC; sw283 Packaged: 2017-09-09 07:02:31 UTC; sw283
Repository: CRAN Repository: CRAN
Date/Publication: 2017-08-29 11:55:09 UTC Date/Publication: 2017-09-09 16:26:09 UTC
c7596ce5c0e893fce457f6365ce48834 *ChangeLog c26bec397cf6ee2a88d4c295945b8cab *ChangeLog
51b6fb01bdd37b7fed320bb84a8d8e97 *DESCRIPTION ac713195702fe496609fec839dbdcc9f *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2 eb723b61539feef013de476e68b5c50a *GPL-2
3d1bd78b1a6f876c1160b7e64ea79cc0 *NAMESPACE 3d1bd78b1a6f876c1160b7e64ea79cc0 *NAMESPACE
d181df8bdf90150190e6defeffae1faf *R/bam.r 043efd04d042fba89ddd6781c0220ab9 *R/bam.r
7b419683b0948cf6da009f614078fe90 *R/coxph.r 7b419683b0948cf6da009f614078fe90 *R/coxph.r
777a0d67a1f7fa14bf87bc312064061b *R/efam.r 777a0d67a1f7fa14bf87bc312064061b *R/efam.r
f610b2695c525c2d9d8f8174e6760548 *R/fast-REML.r f610b2695c525c2d9d8f8174e6760548 *R/fast-REML.r
...@@ -15,10 +15,10 @@ c934ba653cb1a904132b3f15a01e41c5 *R/gamm.r ...@@ -15,10 +15,10 @@ c934ba653cb1a904132b3f15a01e41c5 *R/gamm.r
48c3f7080c31084058aaf3de06015e06 *R/mgcv.r 48c3f7080c31084058aaf3de06015e06 *R/mgcv.r
2feca0dc9d354de7bc707c67a5891364 *R/misc.r 2feca0dc9d354de7bc707c67a5891364 *R/misc.r
03772998ab05887f2eeae10dd6efe983 *R/mvam.r 03772998ab05887f2eeae10dd6efe983 *R/mvam.r
9e65a5249f1be974bcfd7a4942f92364 *R/plots.r ca13f63d1112d31258eb3c68464fa653 *R/plots.r
bdc529fe735f22e20aa0cde6708cde25 *R/smooth.r bdc529fe735f22e20aa0cde6708cde25 *R/smooth.r
666d7fd36fda68b928993d5388b0d7fc *R/soap.r 666d7fd36fda68b928993d5388b0d7fc *R/soap.r
76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r aeacfacc81c1726bdaee38e7f173c844 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda 40874e3ced720a596750f499ded8a60a *data/columb.rda
794d1e804e26577500008820b9196f24 *inst/CITATION 794d1e804e26577500008820b9196f24 *inst/CITATION
...@@ -121,7 +121,7 @@ ee9352ba4c531a8def16deddcab9a9fd *man/pdIdnot.Rd ...@@ -121,7 +121,7 @@ ee9352ba4c531a8def16deddcab9a9fd *man/pdIdnot.Rd
8bc429d92aa9f58c4c43f2852e1f8123 *man/pdTens.Rd 8bc429d92aa9f58c4c43f2852e1f8123 *man/pdTens.Rd
1721f1b266d9e14827e8226e2cb74a81 *man/pen.edf.Rd 1721f1b266d9e14827e8226e2cb74a81 *man/pen.edf.Rd
931c3aefb8b5e42aa230cfedd281bed1 *man/place.knots.Rd 931c3aefb8b5e42aa230cfedd281bed1 *man/place.knots.Rd
073b8d671e322ceaa2c6e04cddb41344 *man/plot.gam.Rd 74b9891dcf56eae77cbd11b04fa0623a *man/plot.gam.Rd
c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
1a9d83c9fc67e5f0fc85d66d3112f4ef *man/predict.bam.Rd 1a9d83c9fc67e5f0fc85d66d3112f4ef *man/predict.bam.Rd
6ea139c7a291c7dbf897a3f39afc97c0 *man/predict.gam.Rd 6ea139c7a291c7dbf897a3f39afc97c0 *man/predict.gam.Rd
...@@ -187,19 +187,19 @@ ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po ...@@ -187,19 +187,19 @@ ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po
0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c 0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c
dba99f7d7cc412dd9255f6307c8c7fa7 *src/gdi.c dba99f7d7cc412dd9255f6307c8c7fa7 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h 2436f9b328e80370ce2203dbf1dd813c *src/general.h
71a592d681b3de5e77c9e9dd4117efde *src/init.c 6ea9eff0b4b804a4e3ebf830c0ecacb9 *src/init.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c a9151b5236852eef9e6947590bfcf88a *src/magic.c
654ff83187dc0f7ef4e085f3348f70d2 *src/mat.c 654ff83187dc0f7ef4e085f3348f70d2 *src/mat.c
0545dabf3524a110d616ea5e6373295d *src/matrix.c 0545dabf3524a110d616ea5e6373295d *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h 6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
48cde0e19d5dd54b131ba66c777c0ec2 *src/mgcv.c 48cde0e19d5dd54b131ba66c777c0ec2 *src/mgcv.c
547f557167fca7b219a55d0535dea74a *src/mgcv.h f2eacd39b434ddef7a3bf9d53fb6fe77 *src/mgcv.h
97e3717e95a70b1470b4c3071e144d17 *src/misc.c 97e3717e95a70b1470b4c3071e144d17 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c 465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
8f480dc455f9ff011c3e8f059efec2c5 *src/qp.c 8f480dc455f9ff011c3e8f059efec2c5 *src/qp.c
cd563899be5b09897d1bf36a7889caa0 *src/qp.h cd563899be5b09897d1bf36a7889caa0 *src/qp.h
d5673b88f6f3d85c62a1337f49abba24 *src/soap.c d5673b88f6f3d85c62a1337f49abba24 *src/soap.c
50c9ed355507971d3d764fbbb90c8a25 *src/sparse-smooth.c dcac8c02b5f9be28d13efc834fc88d55 *src/sparse-smooth.c
fe0444bece322bc229e46b3d1c150779 *src/tprs.c fe0444bece322bc229e46b3d1c150779 *src/tprs.c
5bd85bf0319a7b7c755cf49c91a7cd94 *src/tprs.h 5bd85bf0319a7b7c755cf49c91a7cd94 *src/tprs.h
38e593a85a6fd0bb4fbed836f3361406 *tests/bam.R 38e593a85a6fd0bb4fbed836f3361406 *tests/bam.R
......
...@@ -94,7 +94,7 @@ qr.up <- function(arg) { ...@@ -94,7 +94,7 @@ qr.up <- function(arg) {
dd <- dDeta(y,mu,weights,theta=arg$theta,arg$family,0) dd <- dDeta(y,mu,weights,theta=arg$theta,arg$family,0)
## note: no handling of infinities and wz case yet ## note: no handling of infinities and wz case yet
w <- dd$EDeta2 * .5 w <- dd$EDeta2 * .5
w <- w #w <- w
z <- (eta1-arg$offset[ind]) - dd$Deta.EDeta2 z <- (eta1-arg$offset[ind]) - dd$Deta.EDeta2
good <- is.finite(z)&is.finite(w) good <- is.finite(z)&is.finite(w)
w[!good] <- 0 ## drop if !good w[!good] <- 0 ## drop if !good
...@@ -530,7 +530,7 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL, ...@@ -530,7 +530,7 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
if (iter>1) { if (iter>1) {
## form eta = X%*%beta ## form eta = X%*%beta
eta <- Xbd(G$Xd,coef,G$kd,G$ks,G$ts,G$dt,G$v,G$qc,G$drop) eta <- Xbd(G$Xd,coef,G$kd,G$ks,G$ts,G$dt,G$v,G$qc,G$drop) + offset
lsp.full <- G$lsp0 lsp.full <- G$lsp0
if (n.sp>0) lsp.full <- lsp.full + if (is.null(G$L)) lsp[1:n.sp] else G$L %*% lsp[1:n.sp] if (n.sp>0) lsp.full <- lsp.full + if (is.null(G$L)) lsp[1:n.sp] else G$L %*% lsp[1:n.sp]
Sb <- Sl.Sb(Sl,lsp.full,prop$beta) ## store S beta to allow rapid step halving Sb <- Sl.Sb(Sl,lsp.full,prop$beta) ## store S beta to allow rapid step halving
...@@ -580,9 +580,21 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL, ...@@ -580,9 +580,21 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
dd <- dDeta(y,mu,G$w,theta=theta,family,0) dd <- dDeta(y,mu,G$w,theta=theta,family,0)
## note: no handling of infinities and wz case yet ## note: no handling of infinities and wz case yet
if (rho==0) {
w <- dd$Deta2 * .5 w <- dd$Deta2 * .5
w <- w
z <- (eta-offset) - dd$Deta.Deta2 z <- (eta-offset) - dd$Deta.Deta2
} else { ## use fisher weights
w <- dd$EDeta2 * .5
z <- (eta-offset) - dd$Deta.EDeta2
}
#if (rho!=0) {
# ind <- which(w<0)
# if (length(ind)>0) { ## substitute Fisher weights
# w[ind] <- dd$EDeta2[ind] * .5
# z[ind] <- (eta[ind]-offset[ind]) - dd$Deta.EDeta2[ind]
# }
#}
good <- is.finite(z)&is.finite(w) good <- is.finite(z)&is.finite(w)
w[!good] <- 0 ## drop if !good w[!good] <- 0 ## drop if !good
z[!good] <- 0 ## irrelevant z[!good] <- 0 ## irrelevant
...@@ -926,7 +938,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas ...@@ -926,7 +938,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
## note: no handling of infinities and wz case yet ## note: no handling of infinities and wz case yet
w <- dd$EDeta2 * .5 w <- dd$EDeta2 * .5
w <- w #w <- w
z <- (eta1-offset[ind]) - dd$Deta.EDeta2 z <- (eta1-offset[ind]) - dd$Deta.EDeta2
good <- is.finite(z)&is.finite(w) good <- is.finite(z)&is.finite(w)
w[!good] <- 0 ## drop if !good w[!good] <- 0 ## drop if !good
...@@ -1061,7 +1073,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas ...@@ -1061,7 +1073,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
if (efam && iter>1) { ## estimate theta if (efam && iter>1) { ## estimate theta
scale1 <- if (!is.null(family$scale)) family$scale else scale scale1 <- if (!is.null(family$scale)) family$scale else scale
if (family$n.theta>0||scale<0) theta <- estimate.theta(theta,family,y,linkinv(eta),scale=scale1,wt=G$w,tol=1e-7) if (family$n.theta>0||scale<0) theta <- estimate.theta(theta,family,G$y,linkinv(eta),scale=scale1,wt=G$w,tol=1e-7)
if (!is.null(family$scale) && family$scale<0) { if (!is.null(family$scale) && family$scale<0) {
scale <- exp(theta[family$n.theta+1]) scale <- exp(theta[family$n.theta+1])
theta <- theta[1:family$n.theta] theta <- theta[1:family$n.theta]
...@@ -1177,7 +1189,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas ...@@ -1177,7 +1189,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
if (!is.null(posr$deviance)) object$deviance <- posr$deviance if (!is.null(posr$deviance)) object$deviance <- posr$deviance
if (!is.null(posr$null.deviance)) object$null.deviance <- posr$null.deviance if (!is.null(posr$null.deviance)) object$null.deviance <- posr$null.deviance
} }
if (is.null(object$null.deviance)) object$null.deviance <- sum(family$dev.resids(y,weighted.mean(y,G$w),G$w,theta)) if (is.null(object$null.deviance)) object$null.deviance <- sum(family$dev.resids(G$y,weighted.mean(G$y,G$w),G$w,theta))
} }
if (!conv) if (!conv)
...@@ -1647,10 +1659,17 @@ predict.bamd <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,excl ...@@ -1647,10 +1659,17 @@ predict.bamd <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,excl
if (!is.null(exclude)) warning("exclude ignored by discrete prediction at present") if (!is.null(exclude)) warning("exclude ignored by discrete prediction at present")
## newdata has to be processed first to avoid, e.g. dropping different subsets of data
## for parametric and smooth components....
newdata <- predict.gam(object,newdata=newdata,type="newdata",se.fit=se.fit,terms=terms,exclude=exclude, newdata <- predict.gam(object,newdata=newdata,type="newdata",se.fit=se.fit,terms=terms,exclude=exclude,
block.size=block.size,newdata.guaranteed=newdata.guaranteed, block.size=block.size,newdata.guaranteed=newdata.guaranteed,
na.action=na.action,...) na.action=na.action,...)
## Next line needed to avoid treating newdata as a model frame and then
## having incorrect labels for offset, for example....
attr(newdata,"terms") <- NULL
## Parametric terms have to be dealt with safely, but without forming all terms ## Parametric terms have to be dealt with safely, but without forming all terms
## or a full model matrix. Strategy here is to use predict.gam, having removed ## or a full model matrix. Strategy here is to use predict.gam, having removed
## key smooth related components from model object, so that it appears to be ## key smooth related components from model object, so that it appears to be
...@@ -1793,8 +1812,6 @@ predict.bamd <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,excl ...@@ -1793,8 +1812,6 @@ predict.bamd <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,excl
tero <- function(sm) { tero <- function(sm) {
## te smooth spec re-order so that largest marginal is last. ## te smooth spec re-order so that largest marginal is last.
maxd <- 0 maxd <- 0
...@@ -1806,6 +1823,7 @@ tero <- function(sm) { ...@@ -1806,6 +1823,7 @@ tero <- function(sm) {
ind <- 1:ns;ind[maxi] <- ns;ind[ns] <- maxi ind <- 1:ns;ind[maxi] <- ns;ind[ns] <- maxi
sm$margin <- sm$margin[ind] sm$margin <- sm$margin[ind]
sm$fix <- sm$fix[ind] sm$fix <- sm$fix[ind]
if (!is.null(sm$mc)) sm$mc <- sm$mc[ind]
sm$term <- rep("",0) sm$term <- rep("",0)
for (i in 1:ns) sm$term <- c(sm$term,sm$margin[[i]]$term) for (i in 1:ns) sm$term <- c(sm$term,sm$margin[[i]]$term)
sm$label <- paste0(substr(sm$label,1,3),paste0(sm$term,collapse=","),")",collapse="") sm$label <- paste0(substr(sm$label,1,3),paste0(sm$term,collapse=","),")",collapse="")
......
...@@ -1071,7 +1071,7 @@ plot.gam <- function(x,residuals=FALSE,rug=NULL,se=TRUE,pages=0,select=NULL,scal ...@@ -1071,7 +1071,7 @@ plot.gam <- function(x,residuals=FALSE,rug=NULL,se=TRUE,pages=0,select=NULL,scal
if (is.null(w.resid)) { ## produce working resids if info available if (is.null(w.resid)) { ## produce working resids if info available
if (is.null(x$residuals)||is.null(x$weights)) partial.resids <- FALSE else { if (is.null(x$residuals)||is.null(x$weights)) partial.resids <- FALSE else {
wr <- sqrt(x$weights) wr <- sqrt(x$weights)
w.resid <- x$residuals*wr/mean(wr) # weighted working residuals w.resid <- x$residuals*wr#/mean(wr) # weighted working residuals
} }
} }
if (partial.resids) fv.terms <- predict(x,type="terms") ## get individual smooth effects if (partial.resids) fv.terms <- predict(x,type="terms") ## get individual smooth effects
......
...@@ -100,6 +100,10 @@ apply.spline <- function(spl,y) { ...@@ -100,6 +100,10 @@ apply.spline <- function(spl,y) {
kd.vis <- function(kd,X,cex=.5) { kd.vis <- function(kd,X,cex=.5) {
## code visualizes a kd tree for points in rows of X ## code visualizes a kd tree for points in rows of X
## kd <- kd.tree(X) produces correct tree. ## 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") if (ncol(X)!=2) stop("only deals with 2D case")
##n <- nrow(X) ##n <- nrow(X)
d <- ncol(X) d <- ncol(X)
...@@ -156,53 +160,82 @@ nearest <- function(k,X,gt.zero = FALSE,get.a=FALSE) { ...@@ -156,53 +160,82 @@ nearest <- function(k,X,gt.zero = FALSE,get.a=FALSE) {
list(ni=ni,dist=dist,a=a) list(ni=ni,dist=dist,a=a)
} # nearest } # 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) { kd.tree <- function(X) {
## function to obtain kd tree for points in rows of X ## function to obtain kd tree for points in rows of X
n <- nrow(X) ## number of points kd <- .Call(C_Rkdtree,X)
d <- ncol(X) ## dimension of points kd
## 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.nearest <- function(kd,X,x,k) { 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.
attr(X,"kd_ptr") <- attr(kd,"kd_ptr")
nei <- .Call(C_Rkdnearest,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)
attr(X,"kd_ptr") <- attr(kd,"kd_ptr")
off <- rep(as.integer(0),nrow(x)+1)
ni <- .Call(C_Rkradius,X,x,as.double(r),off) + 1
list(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 ## 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 ## (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. ## k nearest neighbours in X, to the points in the rows of x.
## outputs: ni[i,] lists k nearest neighbours of X[i,]. ## outputs: ni[i,] lists k nearest neighbours of X[i,].
## dost[i,] is distance to those neighbours. ## dost[i,] is distance to those neighbours.
## note R indexing of output ## note R indexing of output
n <- nrow(X) # n <- nrow(X)
m <- nrow(x) # m <- nrow(x)
ni <- matrix(0,m,k) # 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), # 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)) # 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)) # list(ni=matrix(oo$ni+1,m,k),dist=matrix(oo$dist,m,k))
} #}
kd.radius <- function(kd,X,x,r) { #kd.radius <- function(kd,X,x,r) {
## find all points in kd tree (kd,X) in radius r of points in x. ## find all points in kd tree (kd,X) in radius r of points in x.
## kd should come from kd.tree(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)] ## neighbours of x[i,] in X are the rows given by ni[off[i]:(off[i+1]-1)]
m <- nrow(x); # m <- nrow(x);
off <- rep(0,m+1) # off <- rep(0,m+1)
## do the work... ## 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)), # 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)) # as.integer(m),off=as.integer(off),ni=as.integer(0),op=as.integer(0))
off <- oo$off # off <- oo$off
ni <- rep(0,off[m+1]) # ni <- rep(0,off[m+1])
## extract to R and clean up... ## 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)), # 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)) # 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. # list(off=off+1,ni=oo$ni+1) ## note R indexing here.
} ## kd.radius #} ## kd.radius
tieMatrix <- function(x) { tieMatrix <- function(x) {
## takes matrix x, and produces sparse matrix P that maps list of unique ## takes matrix x, and produces sparse matrix P that maps list of unique
......
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
...@@ -22,8 +22,8 @@ ...@@ -22,8 +22,8 @@
\item{residuals}{If \code{TRUE} then partial residuals are added to plots of 1-D smooths. If \code{FALSE} \item{residuals}{If \code{TRUE} then partial residuals are added to plots of 1-D smooths. If \code{FALSE}
then no residuals are added. If this is an array of the correct length then it is used as the array of then no residuals are added. If this is an array of the correct length then it is used as the array of
residuals to be used for producing partial residuals. If \code{TRUE} then the residuals to be used for producing partial residuals. If \code{TRUE} then the
residuals are the working residuals from the IRLS iteration weighted by the residuals are the working residuals from the IRLS iteration weighted by the (square root)
IRLS weights. Partial residuals for a smooth term are the IRLS weights, in order that they have constant variance if the model is correct. Partial residuals for a smooth term are the
residuals that would be obtained by dropping the term concerned from the model, while leaving all other residuals that would be obtained by dropping the term concerned from the model, while leaving all other
estimates fixed (i.e. the estimates for the term plus the residuals).} estimates fixed (i.e. the estimates for the term plus the residuals).}
......
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
...@@ -18,6 +18,9 @@ R_CallMethodDef CallMethods[] = { ...@@ -18,6 +18,9 @@ R_CallMethodDef CallMethods[] = {
{ "mgcv_Rpbacksolve",(DL_FUNC)&mgcv_Rpbacksolve,3}, { "mgcv_Rpbacksolve",(DL_FUNC)&mgcv_Rpbacksolve,3},
{ "mgcv_Rpcross",(DL_FUNC)&mgcv_Rpcross,3}, { "mgcv_Rpcross",(DL_FUNC)&mgcv_Rpcross,3},
{ "mgcv_madi",(DL_FUNC)&mgcv_madi,4}, { "mgcv_madi",(DL_FUNC)&mgcv_madi,4},
{ "Rkdtree",(DL_FUNC)&Rkdtree,1},
{"Rkdnearest",(DL_FUNC)&Rkdnearest,3},
{"Rkradius",(DL_FUNC)&Rkradius,4},
{NULL, NULL, 0} {NULL, NULL, 0}
}; };
...@@ -64,9 +67,9 @@ R_CMethodDef CEntries[] = { ...@@ -64,9 +67,9 @@ R_CMethodDef CEntries[] = {
{"rksos",(DL_FUNC)&rksos,3}, {"rksos",(DL_FUNC)&rksos,3},
{"gen_tps_poly_powers",(DL_FUNC)&gen_tps_poly_powers,4}, {"gen_tps_poly_powers",(DL_FUNC)&gen_tps_poly_powers,4},
{"k_nn",(DL_FUNC)&k_nn,8}, {"k_nn",(DL_FUNC)&k_nn,8},
{"Rkdtree",(DL_FUNC)&Rkdtree,5}, // {"Rkdtree",(DL_FUNC)&Rkdtree,5},
{"Rkdnearest",(DL_FUNC)&Rkdnearest,9}, //{"Rkdnearest",(DL_FUNC)&Rkdnearest,9},
{"Rkradius",(DL_FUNC)&Rkradius,9}, //{"Rkradius",(DL_FUNC)&Rkradius,9},
{"sspl_construct",(DL_FUNC)&sspl_construct,9}, {"sspl_construct",(DL_FUNC)&sspl_construct,9},
{"sspl_mapply",(DL_FUNC)&sspl_mapply,9}, {"sspl_mapply",(DL_FUNC)&sspl_mapply,9},
{"tri2nei",(DL_FUNC)&tri2nei,5}, {"tri2nei",(DL_FUNC)&tri2nei,5},
......
...@@ -190,9 +190,12 @@ typedef struct { ...@@ -190,9 +190,12 @@ typedef struct {
void k_newn_work(double *Xm,kdtree_type kd,double *X,double *dist,int *ni,int*m,int *n,int *d,int *k); void k_newn_work(double *Xm,kdtree_type kd,double *X,double *dist,int *ni,int*m,int *n,int *d,int *k);
void k_nn(double *X,double *dist,double *a,int *ni,int *n,int *d,int *k,int *get_a); void k_nn(double *X,double *dist,double *a,int *ni,int *n,int *d,int *k,int *get_a);
void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat); //void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat);
void Rkdnearest(double *X,int *idat,double *ddat,int *n,double *x, int *m, int *ni, double *dist,int *k); SEXP Rkdtree(SEXP x);
void Rkradius(double *r,int *idat,double *ddat,double *X,double *x,int *m,int *off,int *ni,int *op); //void Rkdnearest(double *X,int *idat,double *ddat,int *n,double *x, int *m, int *ni, double *dist,int *k);
SEXP Rkdnearest(SEXP Xr, SEXP xr,SEXP k);
//void Rkradius(double *r,int *idat,double *ddat,double *X,double *x,int *m,int *off,int *ni,int *op);
SEXP Rkradius(SEXP Xr, SEXP xr,SEXP rr,SEXP offr);
double xidist(double *x,double *X,int i,int d, int n); double xidist(double *x,double *X,int i,int d, int n);
int closest(kdtree_type *kd, double *X,double *x,int n,int *ex,int nex); int closest(kdtree_type *kd, double *X,double *x,int n,int *ex,int nex);
void kd_tree(double *X,int *n, int *d,kdtree_type *kd); void kd_tree(double *X,int *n, int *d,kdtree_type *kd);
......
File mode changed from 100755 to 100644
...@@ -348,7 +348,7 @@ void kd_tree(double *X,int *n, int *d,kdtree_type *kd) { ...@@ -348,7 +348,7 @@ void kd_tree(double *X,int *n, int *d,kdtree_type *kd) {
} /* end of kd_tree */ } /* end of kd_tree */
void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat) { void Rkdtree0(double *X,int *n, int *d,int *idat,double *ddat) {
/* Routine to export kdtree data to R /* Routine to export kdtree data to R
m <- 2; m <- 2;
while (m<n) m <- m*2 while (m<n) m <- m*2
...@@ -360,7 +360,8 @@ void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat) { ...@@ -360,7 +360,8 @@ void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat) {
ddat is an nb*2*d vector ddat is an nb*2*d vector
idat is a 2 + 7*nb vector idat is a 2 + 7*nb vector
- together they encode the kd tree - together they encode the kd tree
OLD routine based on .C calls and copying kd tree to
R vectors.
*/ */
kdtree_type kd; kdtree_type kd;
kd_tree(X,n,d,&kd); /* create kd tree */ kd_tree(X,n,d,&kd); /* create kd tree */
...@@ -368,6 +369,57 @@ void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat) { ...@@ -368,6 +369,57 @@ void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat) {
free_kdtree(kd); /* free structure */ free_kdtree(kd); /* free structure */
} }
static void kdFinalizer(SEXP ptr) {
kdtree_type *kd;
kd = (kdtree_type *) R_ExternalPtrAddr(ptr);
//Rprintf("%p %p\n",kd,ptr);
free_kdtree(*kd); /* free tree structure */
//Rprintf("*!");
FREE(kd); /* free kd itself */
}
SEXP Rkdtree(SEXP x) {
/* Create a kd tree and return a handle for it to R as attribute
of (dummy) integer returned here. Idea is that further routines
requiring the tree then don't need to waste time copying it.
library(mgcv)
X <- matrix(runif(3000000),1000000,3)
kd <- mgcv:::kd.tree(X)
x <- matrix(runif(300),100,3)
nei <- mgcv:::kd.nearest(kd,X,x,1)
bb <- mgcv:::kd.radius(kd,X,x,.01)
where X is a matrix each row of which is a d-point.
*/
kdtree_type *kd;
double *X;
int *dim,n,d;
SEXP DIM, ans, ptr;
static SEXP kd_symb = NULL;
if (!kd_symb) kd_symb = install("kd_ptr"); /* register symbol for attribute */
X = REAL(x);
DIM = getAttrib(x, install("dim"));
dim=INTEGER(DIM);
n=dim[0];d=dim[1];
kd = (kdtree_type *) CALLOC((size_t)1,sizeof(kdtree_type));
kd_tree(X,&n,&d,kd); /* create kd tree */
/* something to attach the pointer to as an attribute.... */
ans = PROTECT(allocVector(INTSXP, 1));
INTEGER(ans)[0] =0;
/* make pointer to kd (i.e. a pointer to a pointer - a handle) ...*/
ptr = R_MakeExternalPtr(kd,R_NilValue, R_NilValue);
PROTECT(ptr);
/* Register the routine to call when R object to
which ptr belongs is destroyed... */
R_RegisterCFinalizerEx(ptr, kdFinalizer, TRUE);
/* attach ptr as attibute to 'ans' ... */
setAttrib(ans, kd_symb, ptr);
//Rprintf("%p %p\n",kd,ptr);
UNPROTECT(2);
return ans;
} /* Rkdtree */
void update_heap(double *h,int *ind,int n) { void update_heap(double *h,int *ind,int n) {
/* h contains n numbers, such that h[i] > h[2*i+1] and /* h contains n numbers, such that h[i] > h[2*i+1] and
...@@ -677,9 +729,59 @@ void k_radius(double r, kdtree_type kd, double *X,double *x,int *list,int *nlist ...@@ -677,9 +729,59 @@ void k_radius(double r, kdtree_type kd, double *X,double *x,int *list,int *nlist
} }
} /* k_radius */ } /* k_radius */
SEXP Rkradius(SEXP Xr,SEXP xr,SEXP rr,SEXP offr) {
/* Xr is matrix of points with attribute "kd_ptr" which is a handle to a kd tree.
xr is a matrix of m rows. The routine finds all elements of Xr within r of each
row of xr. off is an m+1 vector. Returns a vector ni such that ni[off[i]:(off[i+1]-1)]
contains the indices (rows) in Xr of the neighbours of the ith row of xr.
void Rkradius(double *r,int *idat,double *ddat,double *X,double *x,int *m,int *off,int *ni,int *op) { */
/* Given kd tree defined by idat, ddat and X, from R, this routine finds all points in double *X,*x,*dis,*r,*xx;
kdtree_type *kd;
int *dim,m,d,*off,*nei,*list,nn,i,j,n_buff=0,nlist,*ni;
SEXP DIM,ptr,neir;
static SEXP kd_symb = NULL, dim_sym = NULL;
if (!dim_sym) dim_sym = install("dim");
if (!kd_symb) kd_symb = install("kd_ptr"); /* register symbol for attribute */
//DIM = getAttrib(Xr, dim_sym);
//dim = INTEGER(DIM); n = dim[0];
DIM = getAttrib(xr, dim_sym);
dim = INTEGER(DIM); m = dim[0];
Rprintf("0 ");
X = REAL(Xr);x = REAL(xr);
r = REAL(rr);
ptr = getAttrib(Xr, kd_symb);
kd = (kdtree_type *) R_ExternalPtrAddr(ptr);
d = kd->d; /* dimension */
off = INTEGER(offr);
Rprintf("1 ");
/* get the r-radius neighbour information... */
list = (int *)CALLOC((size_t)kd->n,sizeof(int)); /* list of neighbours of ith point */
n_buff = kd->n*10;
nei = (int *)CALLOC((size_t)n_buff,sizeof(int)); /* global list of neighbours */
xx=x;nn=0;off[0]=0;
for (i=0;i<m;i++) { /* work through points in x */
k_radius(*r, *kd, X,xx,list,&nlist);
if (nn+nlist>n_buff) { /* expand nei */
n_buff *= 2;
nei = (int *)R_chk_realloc(nei,(size_t)n_buff*sizeof(int));
}
for (j=nn;j<nn+nlist;j++) nei[j] = list[j-nn];
nn += nlist;
off[i+1] = nn;
xx += d; /* next point */
}
neir = PROTECT(allocVector(INTSXP, nn));
ni = INTEGER(neir);Rprintf("2 ");
for (dim=nei;dim<nei+nn;dim++,ni++) *ni = *dim;
FREE(list);FREE(nei);
UNPROTECT(1);
return neir;
} /* Rkradius */
void Rkradius0(double *r,int *idat,double *ddat,double *X,double *x,int *m,int *off,int *ni,int *op) {
/* This is old version based on copying tree to R structures, rather than retunring tree as external pointer.
Given kd tree defined by idat, ddat and X, from R, this routine finds all points in
the tree less than distance r from each point in x. x contains the points stored end-to-end. the tree less than distance r from each point in x. x contains the points stored end-to-end.
Routine must be called twice. First with op==0, which does the work, but only returns the Routine must be called twice. First with op==0, which does the work, but only returns the
length required for ni, in off[m+1]. length required for ni, in off[m+1].
...@@ -818,7 +920,7 @@ void k_newn_work(double *Xm,kdtree_type kd,double *X,double *dist,int *ni,int*m, ...@@ -818,7 +920,7 @@ void k_newn_work(double *Xm,kdtree_type kd,double *X,double *dist,int *ni,int*m,
*n = pcount; *n = pcount;
} /* k_newn_work */ } /* k_newn_work */
void Rkdnearest(double *X,int *idat,double *ddat,int *n,double *x, int *m, int *ni, double *dist,int *k) { void Rkdnearest0(double *X,int *idat,double *ddat,int *n,double *x, int *m, int *ni, double *dist,int *k) {
/* given points in n rows of X and a kd tree stored in idat, ddat in R, find the /* given points in n rows of X and a kd tree stored in idat, ddat in R, find the
k neares neighbours to each row of x m by d matrix x. k neares neighbours to each row of x m by d matrix x.
* outputs * outputs
...@@ -834,6 +936,36 @@ void Rkdnearest(double *X,int *idat,double *ddat,int *n,double *x, int *m, int * ...@@ -834,6 +936,36 @@ void Rkdnearest(double *X,int *idat,double *ddat,int *n,double *x, int *m, int *