Commit 960a6316 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Debian changes 1.8-20-1

mgcv (1.8-20-1) unstable; urgency=medium

  * New upstream release
parents f2009880 1b84d0b8
......@@ -4,6 +4,25 @@
Currently deprecated and liable to be removed:
- 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
** bam() now accepts extended families (i.e. nb, tw, ocat etc)
......
Package: mgcv
Version: 1.8-19
Version: 1.8-20
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with Automatic Smoothness
......@@ -18,6 +18,6 @@ LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
NeedsCompilation: yes
Packaged: 2017-08-29 08:25:01 UTC; sw283
Packaged: 2017-09-09 07:02:31 UTC; sw283
Repository: CRAN
Date/Publication: 2017-08-29 11:55:09 UTC
Date/Publication: 2017-09-09 16:26:09 UTC
c7596ce5c0e893fce457f6365ce48834 *ChangeLog
51b6fb01bdd37b7fed320bb84a8d8e97 *DESCRIPTION
c26bec397cf6ee2a88d4c295945b8cab *ChangeLog
ac713195702fe496609fec839dbdcc9f *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
3d1bd78b1a6f876c1160b7e64ea79cc0 *NAMESPACE
d181df8bdf90150190e6defeffae1faf *R/bam.r
043efd04d042fba89ddd6781c0220ab9 *R/bam.r
7b419683b0948cf6da009f614078fe90 *R/coxph.r
777a0d67a1f7fa14bf87bc312064061b *R/efam.r
f610b2695c525c2d9d8f8174e6760548 *R/fast-REML.r
......@@ -15,10 +15,10 @@ c934ba653cb1a904132b3f15a01e41c5 *R/gamm.r
48c3f7080c31084058aaf3de06015e06 *R/mgcv.r
2feca0dc9d354de7bc707c67a5891364 *R/misc.r
03772998ab05887f2eeae10dd6efe983 *R/mvam.r
9e65a5249f1be974bcfd7a4942f92364 *R/plots.r
ca13f63d1112d31258eb3c68464fa653 *R/plots.r
bdc529fe735f22e20aa0cde6708cde25 *R/smooth.r
666d7fd36fda68b928993d5388b0d7fc *R/soap.r
76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r
aeacfacc81c1726bdaee38e7f173c844 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
794d1e804e26577500008820b9196f24 *inst/CITATION
......@@ -121,7 +121,7 @@ ee9352ba4c531a8def16deddcab9a9fd *man/pdIdnot.Rd
8bc429d92aa9f58c4c43f2852e1f8123 *man/pdTens.Rd
1721f1b266d9e14827e8226e2cb74a81 *man/pen.edf.Rd
931c3aefb8b5e42aa230cfedd281bed1 *man/place.knots.Rd
073b8d671e322ceaa2c6e04cddb41344 *man/plot.gam.Rd
74b9891dcf56eae77cbd11b04fa0623a *man/plot.gam.Rd
c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
1a9d83c9fc67e5f0fc85d66d3112f4ef *man/predict.bam.Rd
6ea139c7a291c7dbf897a3f39afc97c0 *man/predict.gam.Rd
......@@ -187,19 +187,19 @@ ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po
0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c
dba99f7d7cc412dd9255f6307c8c7fa7 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
71a592d681b3de5e77c9e9dd4117efde *src/init.c
6ea9eff0b4b804a4e3ebf830c0ecacb9 *src/init.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c
654ff83187dc0f7ef4e085f3348f70d2 *src/mat.c
0545dabf3524a110d616ea5e6373295d *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
48cde0e19d5dd54b131ba66c777c0ec2 *src/mgcv.c
547f557167fca7b219a55d0535dea74a *src/mgcv.h
f2eacd39b434ddef7a3bf9d53fb6fe77 *src/mgcv.h
97e3717e95a70b1470b4c3071e144d17 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
8f480dc455f9ff011c3e8f059efec2c5 *src/qp.c
cd563899be5b09897d1bf36a7889caa0 *src/qp.h
d5673b88f6f3d85c62a1337f49abba24 *src/soap.c
50c9ed355507971d3d764fbbb90c8a25 *src/sparse-smooth.c
dcac8c02b5f9be28d13efc834fc88d55 *src/sparse-smooth.c
fe0444bece322bc229e46b3d1c150779 *src/tprs.c
5bd85bf0319a7b7c755cf49c91a7cd94 *src/tprs.h
38e593a85a6fd0bb4fbed836f3361406 *tests/bam.R
......
......@@ -94,7 +94,7 @@ qr.up <- function(arg) {
dd <- dDeta(y,mu,weights,theta=arg$theta,arg$family,0)
## note: no handling of infinities and wz case yet
w <- dd$EDeta2 * .5
w <- w
#w <- w
z <- (eta1-arg$offset[ind]) - dd$Deta.EDeta2
good <- is.finite(z)&is.finite(w)
w[!good] <- 0 ## drop if !good
......@@ -530,7 +530,7 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
if (iter>1) {
## 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
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
......@@ -580,10 +580,22 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
dd <- dDeta(y,mu,G$w,theta=theta,family,0)
## note: no handling of infinities and wz case yet
w <- dd$Deta2 * .5
w <- w
z <- (eta-offset) - dd$Deta.Deta2
good <- is.finite(z)&is.finite(w)
if (rho==0) {
w <- dd$Deta2 * .5
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)
w[!good] <- 0 ## drop if !good
z[!good] <- 0 ## irrelevant
#dev <- sum(family$dev.resids(G$y,mu,G$w,theta))
......@@ -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
w <- dd$EDeta2 * .5
w <- w
#w <- w
z <- (eta1-offset[ind]) - dd$Deta.EDeta2
good <- is.finite(z)&is.finite(w)
w[!good] <- 0 ## drop if !good
......@@ -1061,7 +1073,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
if (efam && iter>1) { ## estimate theta
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) {
scale <- exp(theta[family$n.theta+1])
theta <- theta[1:family$n.theta]
......@@ -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$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)
......@@ -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")
## 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,
block.size=block.size,newdata.guaranteed=newdata.guaranteed,
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
## 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
......@@ -1793,8 +1812,6 @@ predict.bamd <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,excl
tero <- function(sm) {
## te smooth spec re-order so that largest marginal is last.
maxd <- 0
......@@ -1806,6 +1823,7 @@ tero <- function(sm) {
ind <- 1:ns;ind[maxi] <- ns;ind[ns] <- maxi
sm$margin <- sm$margin[ind]
sm$fix <- sm$fix[ind]
if (!is.null(sm$mc)) sm$mc <- sm$mc[ind]
sm$term <- rep("",0)
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="")
......
......@@ -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(x$residuals)||is.null(x$weights)) partial.resids <- FALSE else {
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
......
......@@ -100,6 +100,10 @@ apply.spline <- function(spl,y) {
kd.vis <- function(kd,X,cex=.5) {
## code visualizes a kd tree for points in rows of X
## kd <- kd.tree(X) produces correct tree.
## this worked with the structures used when
## kd tree was written out to R vecotrs and the read
## back in. Does not work with revised approach (would
## need an explicit helper function to do the write out)
if (ncol(X)!=2) stop("only deals with 2D case")
##n <- nrow(X)
d <- ncol(X)
......@@ -156,53 +160,82 @@ nearest <- function(k,X,gt.zero = FALSE,get.a=FALSE) {
list(ni=ni,dist=dist,a=a)
} # nearest
#kd.tree <- function(X) {
## function to obtain kd tree for points in rows of X
## old version based on writing out tree to R
# n <- nrow(X) ## number of points
# d <- ncol(X) ## dimension of points
# ## compute the number of boxes in the kd tree, nb
# m <- 2;while (m < n) m <- m* 2;
# nb = n * 2 - m %/% 2 - 1;
# if (nb > m-1) nb = m - 1;
# ## compute the storage requirements for the tree
# nd = 1 + d * nb * 2 ## number of doubles
# ni = 3 + 5 * nb + 2*n ## number of integers
# oo <- .C(C_Rkdtree,as.double(X),as.integer(n),as.integer(d),idat = as.integer(rep(0,ni)),
# ddat = as.double(rep(0,nd)))
# list(idat=oo$idat,ddat=oo$ddat)
#}
kd.tree <- function(X) {
## function to obtain kd tree for points in rows of X
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 <- .Call(C_Rkdtree,X)
kd
}
kd.nearest <- function(kd,X,x,k) {
## Finds k nearest neigbours of each row of x within X. X has
## corresponding kd tree kd, stored as an external pointer:
## attribute "kd_ptr" of kd. Returns array of indices to rows in
## X, along with corresponding distances.
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
## (produced by a call to kd.tree(X)), then this routine finds the
## k nearest neighbours in X, to the points in the rows of x.
## outputs: ni[i,] lists k nearest neighbours of X[i,].
## dost[i,] is distance to those neighbours.
## note R indexing of output
n <- nrow(X)
m <- nrow(x)
ni <- matrix(0,m,k)
oo <- .C(C_Rkdnearest,as.double(X),as.integer(kd$idat),as.double(kd$ddat),as.integer(n),as.double(x),
as.integer(m), ni=as.integer(ni), dist=as.double(ni),as.integer(k))
list(ni=matrix(oo$ni+1,m,k),dist=matrix(oo$dist,m,k))
}
kd.radius <- function(kd,X,x,r) {
# n <- nrow(X)
# m <- nrow(x)
# ni <- matrix(0,m,k)
# oo <- .C(C_Rkdnearest,as.double(X),as.integer(kd$idat),as.double(kd$ddat),as.integer(n),as.double(x),
# as.integer(m), ni=as.integer(ni), dist=as.double(ni),as.integer(k))
# list(ni=matrix(oo$ni+1,m,k),dist=matrix(oo$dist,m,k))
#}
#kd.radius <- function(kd,X,x,r) {
## find all points in kd tree (kd,X) in radius r of points in x.
## kd should come from kd.tree(X).
## neighbours of x[i,] in X are the rows given by ni[off[i]:(off[i+1]-1)]
m <- nrow(x);
off <- rep(0,m+1)
# m <- nrow(x);
# off <- rep(0,m+1)
## do the work...
oo <- .C(C_Rkradius,as.double(r),as.integer(kd$idat),as.double(kd$ddat),as.double(X),as.double(t(x)),
as.integer(m),off=as.integer(off),ni=as.integer(0),op=as.integer(0))
off <- oo$off
ni <- rep(0,off[m+1])
# oo <- .C(C_Rkradius,as.double(r),as.integer(kd$idat),as.double(kd$ddat),as.double(X),as.double(t(x)),
# as.integer(m),off=as.integer(off),ni=as.integer(0),op=as.integer(0))
# off <- oo$off
# ni <- rep(0,off[m+1])
## extract to R and clean up...
oo <- .C(C_Rkradius,as.double(r),as.integer(kd$idat),as.double(kd$ddat),as.double(X),as.double(t(x)),
as.integer(m),off=as.integer(off),ni=as.integer(ni),op=as.integer(1))
list(off=off+1,ni=oo$ni+1) ## note R indexing here.
} ## kd.radius
# oo <- .C(C_Rkradius,as.double(r),as.integer(kd$idat),as.double(kd$ddat),as.double(X),as.double(t(x)),
# as.integer(m),off=as.integer(off),ni=as.integer(ni),op=as.integer(1))
# list(off=off+1,ni=oo$ni+1) ## note R indexing here.
#} ## kd.radius
tieMatrix <- function(x) {
## takes matrix x, and produces sparse matrix P that maps list of unique
......
Explanation for binary files inside source package according to
http://lists.debian.org/debian-devel/2013/09/msg00332.html
Files: inst/external/CAex_slots.rda
Documentation: man/CAex.Rd
-> "Difficult" Eigen factorization matrix from contributed email
covering 180 days starting with Jan 2, 2007.
Files: inst/external/KNex_slots.rda
Documentation: man/KNex.Rd
-> Koenker and Ng exampple from SparseM package
Files: inst/external/USCounties_slots.rda
Documentation: man/KNex.Rd
-> US Counties data from Ord (1975)
-- Dirk Eddelbuettel <edd@debian.org>, Thu, 17 Aug 2017 20:59:21 -0500
mgcv (1.8-20-1) unstable; urgency=medium
* New upstream release
-- Dirk Eddelbuettel <edd@debian.org> Mon, 11 Sep 2017 20:53:07 -0500
mgcv (1.8-19-1) unstable; urgency=medium
* New upstream release
......
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 @@
\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
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
IRLS weights. Partial residuals for a smooth term are the
residuals are the working residuals from the IRLS iteration weighted by the (square root)
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
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[] = {
{ "mgcv_Rpbacksolve",(DL_FUNC)&mgcv_Rpbacksolve,3},
{ "mgcv_Rpcross",(DL_FUNC)&mgcv_Rpcross,3},
{ "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}
};
......@@ -64,9 +67,9 @@ R_CMethodDef CEntries[] = {
{"rksos",(DL_FUNC)&rksos,3},
{"gen_tps_poly_powers",(DL_FUNC)&gen_tps_poly_powers,4},
{"k_nn",(DL_FUNC)&k_nn,8},
{"Rkdtree",(DL_FUNC)&Rkdtree,5},
{"Rkdnearest",(DL_FUNC)&Rkdnearest,9},
{"Rkradius",(DL_FUNC)&Rkradius,9},
// {"Rkdtree",(DL_FUNC)&Rkdtree,5},
//{"Rkdnearest",(DL_FUNC)&Rkdnearest,9},
//{"Rkradius",(DL_FUNC)&Rkradius,9},
{"sspl_construct",(DL_FUNC)&sspl_construct,9},
{"sspl_mapply",(DL_FUNC)&sspl_mapply,9},
{"tri2nei",(DL_FUNC)&tri2nei,5},
......
......@@ -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_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 Rkdnearest(double *X,int *idat,double *ddat,int *n,double *x, int *m, int *ni, double *dist,int *k);
void Rkradius(double *r,int *idat,double *ddat,double *X,double *x,int *m,int *off,int *ni,int *op);
//void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat);
SEXP Rkdtree(SEXP x);
//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);
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);
......
File mode changed from 100755 to 100644
......@@ -348,7 +348,7 @@ void kd_tree(double *X,int *n, int *d,kdtree_type *kd) {
} /* 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
m <- 2;
while (m<n) m <- m*2
......@@ -360,7 +360,8 @@ void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat) {
ddat is an nb*2*d vector
idat is a 2 + 7*nb vector
- together they encode the kd tree
OLD routine based on .C calls and copying kd tree to
R vectors.
*/
kdtree_type kd;
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) {
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) {
/* 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
}
} /* 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.
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].
......@@ -818,7 +920,7 @@ void k_newn_work(double *Xm,kdtree_type kd,double *X,double *dist,int *ni,int*m,
*n = pcount;
} /* 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
k neares neighbours to each row of x m by d matrix x.
* outputs
......@@ -834,6 +936,36 @@ void Rkdnearest(double *X,int *idat,double *ddat,int *n,double *x, int *m, int *
FREE(kd.box); /* free storage created by kd_read */
}
SEXP Rkdnearest(SEXP Xr,SEXP xr,SEXP kr) {
/* Takes n by d point matrix Xr, with kd_ptr attribute pointing to the corresponding
kd tree and m by p point matrix xr. Finds the k nearest neighbours in Xr to each
point in xr returning the matrix of nearest neighbours to each point, with the
corresponding distances as an attribute */
double *X,*x,*dis;
kdtree_type *kd;
int *dim,n,m,d,*k,*nei;
SEXP DIM,ptr,neir,dir;
static SEXP kd_symb = NULL, dim_sym = NULL,dist_sym = NULL;
if (!dim_sym) dim_sym = install("dim");if (!dist_sym) dist_sym = install("dist");
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];
X = REAL(Xr);x = REAL(xr);
k = INTEGER(kr);
ptr = getAttrib(Xr, kd_symb);
kd = (kdtree_type *) R_ExternalPtrAddr(ptr);
d = kd->d; /* dimension */
neir = PROTECT(allocMatrix(INTSXP,m,*k));
nei = INTEGER(neir);
dir = PROTECT(allocMatrix(REALSXP,m,*k));
dis = REAL(dir);
k_newn_work(x,*kd,X,dis,nei,&m,&n,&d,k);
setAttrib(neir, dist_sym,dir);
UNPROTECT(2);
return neir;
}
void k_nn_work(kdtree_type kd,double *X,double *dist,int *ni,int *n,int *d,int *k) {
/* Given a kd tree, this routine does the actual work of finding the nearest neighbours.
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment