Commit 6e1e032b authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Debian changes 1.8-22-1

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

  * New upstream release
parents 960a6316 a1ae9b9e
...@@ -4,6 +4,28 @@ ...@@ -4,6 +4,28 @@
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-22
* Fix of bug whereby testing for OpenMP and nthreads>1 in bam, would fail if
OpenMP was missing.
1.8-21
* When functions were added to families within mgcv some very large
environments could end up attached to those functions, for no good reason.
The problem originated from the dispatch of the generic 'fix.family.link'
and then propagated via fix.family.var and fix.family.ls. This is now avoided,
resulting in smaller gam objects on disk and lower R memory usage.
Thanks to Niels Richard Hansen for uncovering this.
* Another bug fix for prediction from discrete fit bam models with an offset,
this time when there were more than 50000 data. Also fix to bam fitting when
the number of data was an integer multiple of the chunk size + 1.
* check.term was missing a 'stop' so that some unhandled nesting structures
in bam(...,discrete=TRUE) failed with an unhelpful error, instead of a
helpful one. Fixed.
1.8-20 1.8-20
* bam(,discrete=TRUE) could produce garbage with ti(x,z,k=c(6,5),mc=c(T,F)) * bam(,discrete=TRUE) could produce garbage with ti(x,z,k=c(6,5),mc=c(T,F))
......
Package: mgcv Package: mgcv
Version: 1.8-20 Version: 1.8-22
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-09-09 07:02:31 UTC; sw283 Packaged: 2017-09-18 10:38:41 UTC; sw283
Repository: CRAN Repository: CRAN
Date/Publication: 2017-09-09 16:26:09 UTC Date/Publication: 2017-09-19 00:27:56 UTC
c26bec397cf6ee2a88d4c295945b8cab *ChangeLog bc8bc8be9cec065f0c7ec0cf670516f2 *ChangeLog
ac713195702fe496609fec839dbdcc9f *DESCRIPTION ad2a852c227dd80d86b24784d12dba42 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2 eb723b61539feef013de476e68b5c50a *GPL-2
3d1bd78b1a6f876c1160b7e64ea79cc0 *NAMESPACE 9fa89dc9361930dca536d4e011e4e605 *NAMESPACE
043efd04d042fba89ddd6781c0220ab9 *R/bam.r b051474f30e8e779f2440e8ecd5bbd51 *R/bam.r
7b419683b0948cf6da009f614078fe90 *R/coxph.r 7b419683b0948cf6da009f614078fe90 *R/coxph.r
777a0d67a1f7fa14bf87bc312064061b *R/efam.r 777a0d67a1f7fa14bf87bc312064061b *R/efam.r
f610b2695c525c2d9d8f8174e6760548 *R/fast-REML.r dfdb821247da3e780de0d4599b88735d *R/fast-REML.r
4dc75c18c9feceb73d207b8ac01b2f5f *R/gam.fit3.r 0a9cc97524994d1563b93769284266c4 *R/gam.fit3.r
60aa872ec8901e4dc6273e653f77c83f *R/gam.fit4.r f0fd93b43b36abf05e6666b4a71c20fa *R/gam.fit4.r
1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r 1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r
b99d9de0c46a230e081b9bbbbbb3c617 *R/gamlss.r b99d9de0c46a230e081b9bbbbbb3c617 *R/gamlss.r
c934ba653cb1a904132b3f15a01e41c5 *R/gamm.r c934ba653cb1a904132b3f15a01e41c5 *R/gamm.r
10facb791e4cfd123d183f05660119c6 *R/jagam.r 10facb791e4cfd123d183f05660119c6 *R/jagam.r
48c3f7080c31084058aaf3de06015e06 *R/mgcv.r b9084037d058ef5d787e2818fbc29bae *R/mgcv.r
2feca0dc9d354de7bc707c67a5891364 *R/misc.r 2feca0dc9d354de7bc707c67a5891364 *R/misc.r
03772998ab05887f2eeae10dd6efe983 *R/mvam.r 03772998ab05887f2eeae10dd6efe983 *R/mvam.r
ca13f63d1112d31258eb3c68464fa653 *R/plots.r ca13f63d1112d31258eb3c68464fa653 *R/plots.r
...@@ -42,7 +42,7 @@ fd0cfd64be579f9fbecdbb7f2b8ec1eb *man/Sl.initial.repara.Rd ...@@ -42,7 +42,7 @@ fd0cfd64be579f9fbecdbb7f2b8ec1eb *man/Sl.initial.repara.Rd
60670020425f8749b81a8d8c3f168880 *man/Sl.setup.Rd 60670020425f8749b81a8d8c3f168880 *man/Sl.setup.Rd
69ae63833438a3af2963e39482f1d72f *man/Tweedie.Rd 69ae63833438a3af2963e39482f1d72f *man/Tweedie.Rd
8087ab00d10b44c99c37f49bf90e19cd *man/anova.gam.Rd 8087ab00d10b44c99c37f49bf90e19cd *man/anova.gam.Rd
a7e0ce83164f1e34d0a9d99d0f9853f3 *man/bam.Rd 6180a4e9ea206e1a350b993dade0a869 *man/bam.Rd
ab5e37c3bf8803de63b63c3bdc5909cd *man/bam.update.Rd ab5e37c3bf8803de63b63c3bdc5909cd *man/bam.update.Rd
cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd
745cbf31eb14fc1c5916fc634c74d998 *man/bug.reports.mgcv.Rd 745cbf31eb14fc1c5916fc634c74d998 *man/bug.reports.mgcv.Rd
...@@ -151,7 +151,7 @@ b45d8e71bda4ceca8203dffea577e441 *man/smooth.construct.re.smooth.spec.Rd ...@@ -151,7 +151,7 @@ b45d8e71bda4ceca8203dffea577e441 *man/smooth.construct.re.smooth.spec.Rd
0bfe981f2c3e6ea5b8d5372076ccde53 *man/smooth.construct.sos.smooth.spec.Rd 0bfe981f2c3e6ea5b8d5372076ccde53 *man/smooth.construct.sos.smooth.spec.Rd
3cb4e59f915c8d64b90754eaeeb5a86f *man/smooth.construct.t2.smooth.spec.Rd 3cb4e59f915c8d64b90754eaeeb5a86f *man/smooth.construct.t2.smooth.spec.Rd
8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd 8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd
a088e69cf148d07a78ce95a69759c95c *man/smooth.construct.tp.smooth.spec.Rd c522c270c217e5b83cf8f3e95220a03f *man/smooth.construct.tp.smooth.spec.Rd
ae5e27524e37d57505754639455f18a5 *man/smooth.terms.Rd ae5e27524e37d57505754639455f18a5 *man/smooth.terms.Rd
4c49358e5e6a70d1eca69a2ccaa77609 *man/smooth2random.Rd 4c49358e5e6a70d1eca69a2ccaa77609 *man/smooth2random.Rd
844f9653d74441293d05a24dd3e2876a *man/smoothCon.Rd 844f9653d74441293d05a24dd3e2876a *man/smoothCon.Rd
...@@ -185,19 +185,19 @@ ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po ...@@ -185,19 +185,19 @@ ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars 03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
342aa30c8f6f1e99ffa2576a6f29d7ce *src/coxph.c 342aa30c8f6f1e99ffa2576a6f29d7ce *src/coxph.c
0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c 0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c
dba99f7d7cc412dd9255f6307c8c7fa7 *src/gdi.c f6c4ba80ce1b71be7f4a44fac1b08c28 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h 2436f9b328e80370ce2203dbf1dd813c *src/general.h
6ea9eff0b4b804a4e3ebf830c0ecacb9 *src/init.c a287f92033406194afe0301ee076fda8 *src/init.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c a9151b5236852eef9e6947590bfcf88a *src/magic.c
654ff83187dc0f7ef4e085f3348f70d2 *src/mat.c 654ff83187dc0f7ef4e085f3348f70d2 *src/mat.c
0545dabf3524a110d616ea5e6373295d *src/matrix.c e4cef7f1753153fbab242d1c4d4f7e7f *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h de37b0972199b796654405efc007f25b *src/matrix.h
48cde0e19d5dd54b131ba66c777c0ec2 *src/mgcv.c 8df4b96961491d76989b50856237ee2d *src/mgcv.c
f2eacd39b434ddef7a3bf9d53fb6fe77 *src/mgcv.h b6278c9a4a32bdc16487a96e6a0045b6 *src/mgcv.h
97e3717e95a70b1470b4c3071e144d17 *src/misc.c 97e3717e95a70b1470b4c3071e144d17 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c 465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
8f480dc455f9ff011c3e8f059efec2c5 *src/qp.c 563938b7bb6504ab10df5376c4360220 *src/qp.c
cd563899be5b09897d1bf36a7889caa0 *src/qp.h 073a4b5b0bc6e869c5b35478c47facf1 *src/qp.h
d5673b88f6f3d85c62a1337f49abba24 *src/soap.c d5673b88f6f3d85c62a1337f49abba24 *src/soap.c
dcac8c02b5f9be28d13efc834fc88d55 *src/sparse-smooth.c dcac8c02b5f9be28d13efc834fc88d55 *src/sparse-smooth.c
fe0444bece322bc229e46b3d1c150779 *src/tprs.c fe0444bece322bc229e46b3d1c150779 *src/tprs.c
......
...@@ -85,7 +85,7 @@ influence,logLik,lm,mad, ...@@ -85,7 +85,7 @@ influence,logLik,lm,mad,
make.link,median,model.frame,model.offset,model.matrix,nlm, make.link,median,model.frame,model.offset,model.matrix,nlm,
na.pass,napredict,na.omit,naresid,optim,pchisq,pnorm,pt,pf, na.pass,napredict,na.omit,naresid,optim,pchisq,pnorm,pt,pf,
power,predict,printCoefmat,quantile, power,predict,printCoefmat,quantile,
qbeta,qbinom,qchisq,qnbinom,qgamma,qnorm,qpois,qqline,qqnorm,qqplot, qbeta,qbinom,qcauchy,qchisq,qnbinom,qgamma,qnorm,qpois,qqline,qqnorm,qqplot,
reformulate,residuals, reformulate,residuals,
rbeta,rbinom,rgamma,rnbinom,rnorm,rpois,runif,sd, rbeta,rbinom,rgamma,rnbinom,rnorm,rpois,runif,sd,
termplot,terms.formula,terms,uniroot,var,vcov,weights) termplot,terms.formula,terms,uniroot,var,vcov,weights)
......
...@@ -109,8 +109,8 @@ qr.up <- function(arg) { ...@@ -109,8 +109,8 @@ qr.up <- function(arg) {
wt <- c(wt,w) wt <- c(wt,w)
w <- sqrt(w) w <- sqrt(w)
## note assumption that nt=1 in following qr.update - i.e. each cluster node is strictly serial ## note assumption that nt=1 in following qr.update - i.e. each cluster node is strictly serial
if (b == 1) qrx <- qr.update(w*X[good,],w*z,use.chol=arg$use.chol) if (b == 1) qrx <- qr.update(w*X[good,,drop=FALSE],w*z,use.chol=arg$use.chol)
else qrx <- qr.update(w*X[good,],w*z,qrx$R,qrx$f,qrx$y.norm2,use.chol=arg$use.chol) else qrx <- qr.update(w*X[good,,drop=FALSE],w*z,qrx$R,qrx$f,qrx$y.norm2,use.chol=arg$use.chol)
rm(X);if(arg$gc.level>1) gc() ## X can be large: remove and reclaim rm(X);if(arg$gc.level>1) gc() ## X can be large: remove and reclaim
} }
qrx$dev <- dev;qrx$wt <- wt;qrx$eta <- eta qrx$dev <- dev;qrx$wt <- wt;qrx$eta <- eta
...@@ -191,8 +191,8 @@ check.term <- function(term,rec) { ...@@ -191,8 +191,8 @@ check.term <- function(term,rec) {
ii <- which(rec$vnames%in%term) ii <- which(rec$vnames%in%term)
if (length(ii)) { ## at least one variable already discretized if (length(ii)) { ## at least one variable already discretized
if (length(term)==rec$d[min(ii)]) { ## dimensions match previous discretization if (length(term)==rec$d[min(ii)]) { ## dimensions match previous discretization
if (sum(!(term%in%rec$vnames[ii]))) ("bam can not discretize with this nesting structure") if (sum(!(term%in%rec$vnames[ii]))) stop("bam can not discretize with this nesting structure")
else return(rec$ki[min(ii)]) ## all names match previous - retun index of previous else return(rec$ki[min(ii)]) ## all names match previous - return index of previous
} else stop("bam can not discretize with this nesting structure") } else stop("bam can not discretize with this nesting structure")
} else return(0) ## no match } else return(0) ## no match
} ## check.term } ## check.term
...@@ -953,8 +953,8 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas ...@@ -953,8 +953,8 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
wt[ind] <- w ## wt <- c(wt,w) wt[ind] <- w ## wt <- c(wt,w)
w <- sqrt(w) w <- sqrt(w)
## note that QR may be parallel using npt>1, even under serial accumulation... ## note that QR may be parallel using npt>1, even under serial accumulation...
if (b == 1) qrx <- qr.update(w*X[good,],w*z,use.chol=use.chol,nt=npt) if (b == 1) qrx <- qr.update(w*X[good,,drop=FALSE],w*z,use.chol=use.chol,nt=npt)
else qrx <- qr.update(w*X[good,],w*z,qrx$R,qrx$f,qrx$y.norm2,use.chol=use.chol,nt=npt) else qrx <- qr.update(w*X[good,,drop=FALSE],w*z,qrx$R,qrx$f,qrx$y.norm2,use.chol=use.chol,nt=npt)
rm(X);if(gc.level>1) gc() ## X can be large: remove and reclaim rm(X);if(gc.level>1) gc() ## X can be large: remove and reclaim
} }
if (use.chol) { ## post proc to get R and f... if (use.chol) { ## post proc to get R and f...
...@@ -1855,7 +1855,7 @@ AR.resid <- function(rsd,rho=0,AR.start=NULL) { ...@@ -1855,7 +1855,7 @@ AR.resid <- function(rsd,rho=0,AR.start=NULL) {
bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,na.action=na.omit, bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,na.action=na.omit,
offset=NULL,method="fREML",control=list(),select=FALSE,scale=0,gamma=1,knots=NULL,sp=NULL, offset=NULL,method="fREML",control=list(),select=FALSE,scale=0,gamma=1,knots=NULL,sp=NULL,
min.sp=NULL,paraPen=NULL,chunk.size=10000,rho=0,AR.start=NULL,discrete=FALSE, min.sp=NULL,paraPen=NULL,chunk.size=10000,rho=0,AR.start=NULL,discrete=FALSE,
cluster=NULL,nthreads=NA,gc.level=1,use.chol=FALSE,samfrac=1,coef=NULL, cluster=NULL,nthreads=1,gc.level=1,use.chol=FALSE,samfrac=1,coef=NULL,
drop.unused.levels=TRUE,G=NULL,fit=TRUE,drop.intercept=NULL,...) drop.unused.levels=TRUE,G=NULL,fit=TRUE,drop.intercept=NULL,...)
## Routine to fit an additive model to a large dataset. The model is stated in the formula, ## Routine to fit an additive model to a large dataset. The model is stated in the formula,
...@@ -1894,7 +1894,7 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n ...@@ -1894,7 +1894,7 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
warning("discretization only available with fREML") warning("discretization only available with fREML")
} else { } else {
if (!is.null(cluster)) warning("discrete method does not use parallel cluster - use nthreads instead") if (!is.null(cluster)) warning("discrete method does not use parallel cluster - use nthreads instead")
if (nthreads>1 && !mgcv.omp()) warning("openMP not available: single threaded computation only") if (is.finite(nthreads) && nthreads>1 && !mgcv.omp()) warning("openMP not available: single threaded computation only")
} }
} }
if (inherits(family,"extended.family")) { if (inherits(family,"extended.family")) {
......
...@@ -474,7 +474,7 @@ Sl.addS <- function(Sl,A,rho) { ...@@ -474,7 +474,7 @@ Sl.addS <- function(Sl,A,rho) {
ind <- (Sl[[b]]$start:Sl[[b]]$stop)[Sl[[b]]$ind] ind <- (Sl[[b]]$start:Sl[[b]]$stop)[Sl[[b]]$ind]
if (length(Sl[[b]]$S)==1) { ## singleton if (length(Sl[[b]]$S)==1) { ## singleton
B <- exp(rho[k]);diag <- -1 B <- exp(rho[k]);diag <- -1
er <- .Call(C_mgcv_madi,A,B,ind,diag) dummy <- .Call(C_mgcv_madi,A,B,ind,diag)
## diag(A)[ind] <- diag(A)[ind] + exp(rho[k]) ## penalty is identity times sp ## diag(A)[ind] <- diag(A)[ind] + exp(rho[k]) ## penalty is identity times sp
k <- k + 1 k <- k + 1
} else { } else {
......
This diff is collapsed.
...@@ -570,27 +570,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(), ...@@ -570,27 +570,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
rSncol=as.integer(rSncol),deriv=as.integer(deriv), rSncol=as.integer(rSncol),deriv=as.integer(deriv),
fixed.penalty = as.integer(rp$fixed.penalty),nt=as.integer(control$nthreads), fixed.penalty = as.integer(rp$fixed.penalty),nt=as.integer(control$nthreads),
type=as.integer(gdi.type),dVkk=as.double(rep(0,nSp^2))) type=as.integer(gdi.type),dVkk=as.double(rep(0,nSp^2)))
## test code used to ensure type 0 and type 1 produce identical results, when both should work.
# oot <- .C(C_gdi2,
# X=as.double(x[good,]),E=as.double(Sr),Es=as.double(Eb),rS=as.double(unlist(rS)),
# U1 = as.double(U1),sp=as.double(exp(sp)),theta=as.double(theta),
# z=as.double(z),w=as.double(w),wz=as.double(wz),wf=as.double(wf),Dth=as.double(dd$Dth),
# Det=as.double(dd$Deta),
# Det2=as.double(dd$Deta2),Dth2=as.double(dd$Dth2),Det.th=as.double(dd$Detath),
# Det2.th=as.double(dd$Deta2th),Det3=as.double(dd$Deta3),Det.th2 = as.double(dd$Detath2),
# Det4 = as.double(dd$Deta4),Det3.th=as.double(dd$Deta3th), Deta2.th2=as.double(dd$Deta2th2),
# beta=as.double(coef),b1=as.double(rep(0,ntot*ncol(x))),w1=rep(0,ntot*length(z)),
# D1=as.double(rep(0,ntot)),D2=as.double(rep(0,ntot^2)),
# P=as.double(0),P1=as.double(rep(0,ntot)),P2 = as.double(rep(0,ntot^2)),
# ldet=as.double(1-2*(scoreType=="ML")),ldet1 = as.double(rep(0,ntot)),
# ldet2 = as.double(rep(0,ntot^2)),
# rV=as.double(rep(0,ncol(x)^2)),
# rank.tol=as.double(.Machine$double.eps^.75),rank.est=as.integer(0),
# n=as.integer(sum(good)),q=as.integer(ncol(x)),M=as.integer(nSp),
# n.theta=as.integer(length(theta)), Mp=as.integer(Mp),Enrow=as.integer(rows.E),
# rSncol=as.integer(rSncol),deriv=as.integer(deriv),
# fixed.penalty = as.integer(rp$fixed.penalty),nt=as.integer(control$nthreads),
# type=as.integer(1))
rV <- matrix(oo$rV,ncol(x),ncol(x)) ## rV%*%t(rV)*scale gives covariance matrix rV <- matrix(oo$rV,ncol(x),ncol(x)) ## rV%*%t(rV)*scale gives covariance matrix
rV <- T %*% rV rV <- T %*% rV
## derivatives of coefs w.r.t. sps etc... ## derivatives of coefs w.r.t. sps etc...
......
...@@ -97,8 +97,7 @@ pcls <- function(M) ...@@ -97,8 +97,7 @@ pcls <- function(M)
# M$sp - array of theta_i's # M$sp - array of theta_i's
# Ain, bin and p are not in the object needed to call mgcv.... # Ain, bin and p are not in the object needed to call mgcv....
# #
{ nar<-c(length(M$y),length(M$p),dim(M$Ain)[1],dim(M$C)[1],0) { nar<-c(length(M$y),length(M$p),dim(M$Ain)[1],dim(M$C)[1])
H<-0
## sanity checking ... ## sanity checking ...
if (nrow(M$X)!=nar[1]) stop("nrow(M$X) != length(M$y)") if (nrow(M$X)!=nar[1]) stop("nrow(M$X) != length(M$y)")
if (ncol(M$X)!=nar[2]) stop("ncol(M$X) != length(M$p)") if (ncol(M$X)!=nar[2]) stop("ncol(M$X) != length(M$p)")
...@@ -140,13 +139,12 @@ pcls <- function(M) ...@@ -140,13 +139,12 @@ pcls <- function(M)
M$C <- matrix(0,0,0) M$C <- matrix(0,0,0)
M$p <- rep(0,ncol(M$X)) M$p <- rep(0,ncol(M$X))
nar[2] <- length(M$p) nar[2] <- length(M$p)
nar[4] <- 0
qra.exist <- TRUE qra.exist <- TRUE
if (ncol(M$X)>nrow(M$X)) stop("Model matrix not full column rank") if (ncol(M$X)>nrow(M$X)) stop("Model matrix not full column rank")
} }
} }
o<-.C(C_RPCLS,as.double(M$X),as.double(M$p),as.double(M$y),as.double(M$w),as.double(M$Ain),as.double(M$bin) o<-.C(C_RPCLS,as.double(M$X),as.double(M$p),as.double(M$y),as.double(M$w),as.double(M$Ain),as.double(M$bin)
,as.double(M$C),as.double(H),as.double(Sa),as.integer(M$off),as.integer(df),as.double(M$sp), ,as.double(M$C),as.double(Sa),as.integer(M$off),as.integer(df),as.double(M$sp),
as.integer(length(M$off)),as.integer(nar)) as.integer(length(M$off)),as.integer(nar))
p <- array(o[[2]],length(M$p)) p <- array(o[[2]],length(M$p))
if (qra.exist) p <- qr.qy(qra,c(rep(0,j),p)) if (qra.exist) p <- qr.qy(qra,c(rep(0,j),p))
...@@ -2679,6 +2677,8 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclu ...@@ -2679,6 +2677,8 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclu
nb <- length(object$coefficients) nb <- length(object$coefficients)
} }
if (type=="lpmatrix") block.size <- NULL ## nothing gained by blocking in this case - and offset handling easier this way
## split prediction into blocks, to avoid running out of memory ## split prediction into blocks, to avoid running out of memory
if (is.null(block.size)) { if (is.null(block.size)) {
## use one block as predicting using model frame ## use one block as predicting using model frame
......
mgcv (1.8-22-1) unstable; urgency=medium
* New upstream release
-- Dirk Eddelbuettel <edd@debian.org> Wed, 20 Sep 2017 06:22:38 -0500
mgcv (1.8-20-1) unstable; urgency=medium mgcv (1.8-20-1) unstable; urgency=medium
* New upstream release * New upstream release
......
...@@ -21,7 +21,7 @@ bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL, ...@@ -21,7 +21,7 @@ bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
na.action=na.omit, offset=NULL,method="fREML",control=list(), na.action=na.omit, offset=NULL,method="fREML",control=list(),
select=FALSE,scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL, select=FALSE,scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL,
paraPen=NULL,chunk.size=10000,rho=0,AR.start=NULL,discrete=FALSE, paraPen=NULL,chunk.size=10000,rho=0,AR.start=NULL,discrete=FALSE,
cluster=NULL,nthreads=NA,gc.level=1,use.chol=FALSE,samfrac=1, cluster=NULL,nthreads=1,gc.level=1,use.chol=FALSE,samfrac=1,
coef=NULL,drop.unused.levels=TRUE,G=NULL,fit=TRUE,drop.intercept=NULL,...) coef=NULL,drop.unused.levels=TRUE,G=NULL,fit=TRUE,drop.intercept=NULL,...)
} }
%- maybe also `usage' for other objects documented here. %- maybe also `usage' for other objects documented here.
...@@ -135,7 +135,7 @@ single machine). See details and example code. ...@@ -135,7 +135,7 @@ single machine). See details and example code.
} }
\item{nthreads}{Number of threads to use for non-cluster computation (e.g. combining results from cluster nodes). \item{nthreads}{Number of threads to use for non-cluster computation (e.g. combining results from cluster nodes).
if \code{NA} set to \code{max(1,length(cluster))}.} If \code{NA} set to \code{max(1,length(cluster))}.}
\item{gc.level}{to keep the memory footprint down, it helps to call the garbage collector often, but this takes \item{gc.level}{to keep the memory footprint down, it helps to call the garbage collector often, but this takes
a substatial amount of time. Setting this to zero means that garbage collection only happens when R decides it should. Setting to 2 gives frequent garbage collection. 1 is in between.} a substatial amount of time. Setting this to zero means that garbage collection only happens when R decides it should. Setting to 2 gives frequent garbage collection. 1 is in between.}
......
...@@ -7,9 +7,9 @@ ...@@ -7,9 +7,9 @@
\description{\code{\link{gam}} can use isotropic smooths of any number of variables, specified via terms like \description{\code{\link{gam}} can use isotropic smooths of any number of variables, specified via terms like
\code{s(x,z,bs="tp",m=3)} (or just \code{s(x,z)} as this is the default basis). These terms are based on thin plate \code{s(x,z,bs="tp",m=3)} (or just \code{s(x,z)} as this is the default basis). These terms are based on thin plate
regression splines. \code{m} specifies the order of the derivatives in the thin plate spline penalty. If \code{m} is regression splines. \code{m} specifies the order of the derivatives in the thin plate spline penalty.
a vector of length 2 and the second element is zero, then the penalty null space of the smooth is not included in
the smooth: this is useful if you need to test whether a smooth could be replaced by a linear term, for example. If \code{m} is a vector of length 2 and the second element is zero, then the penalty null space of the smooth is not included in the smooth: this is useful if you need to test whether a smooth could be replaced by a linear term, or construct models with odd nesting structures.
Thin plate regression splines are constructed by starting with the Thin plate regression splines are constructed by starting with the
basis and penalty for a full thin plate spline and then truncating this basis in basis and penalty for a full thin plate spline and then truncating this basis in
...@@ -117,11 +117,23 @@ Wood, S.N. (2003) Thin plate regression splines. J.R.Statist.Soc.B 65(1):95-114 ...@@ -117,11 +117,23 @@ Wood, S.N. (2003) Thin plate regression splines. J.R.Statist.Soc.B 65(1):95-114
\examples{ \examples{
require(mgcv); n <- 100; set.seed(2) require(mgcv); n <- 100; set.seed(2)
x <- runif(n); y <- x + x^2*.2 + rnorm(n) *.1 x <- runif(n); y <- x + x^2*.2 + rnorm(n) *.1
## is smooth significantly different from straight line? ## is smooth significantly different from straight line?
summary(gam(y~s(x,m=c(2,0))+x,method="REML")) ## not quite summary(gam(y~s(x,m=c(2,0))+x,method="REML")) ## not quite
## is smooth significatly different from zero? ## is smooth significatly different from zero?
summary(gam(y~s(x),method="REML")) ## yes! summary(gam(y~s(x),method="REML")) ## yes!
## Fool bam(...,discrete=TRUE) into (strange) nested
## model fit...
set.seed(2) ## simulate some data...
dat <- gamSim(1,n=400,dist="normal",scale=2)
dat$x1a <- dat$x1 ## copy x1 so bam allows 2 copies of x1
## Following removes identifiability problem, by removing
## linear terms from second smooth, and then re-inserting
## the one that was not a duplicate (x2)...
b <- bam(y~s(x0,x1)+s(x1a,x2,m=c(2,0))+x2,data=dat,discrete=TRUE)
## example of knot based tprs... ## example of knot based tprs...
k <- 10; m <- 2 k <- 10; m <- 2
y <- y[order(x)];x <- x[order(x)] y <- y[order(x)];x <- x[order(x)]
......
...@@ -95,91 +95,6 @@ double trAB(double *A,double *B,int *n, int *m) ...@@ -95,91 +95,6 @@ double trAB(double *A,double *B,int *n, int *m)
return(tr); return(tr);
} }
void get_bSb0(double *bSb,double *bSb1, double *bSb2,double *sp,double *E,
double *rS,int *rSncol,int *Enrow, int *q,int *M,double *beta,
double *b1, double *b2,int *deriv)
/* Routine to obtain beta'Sbeta and its derivatives w.r.t. the log smoothing
parameters, this is part of REML calculation...
S= E'E. NOTE: change!!
b1 and b2 contain first and second derivatives of q-vector beta w.r.t.
\pho_k. They are packed as follows....
* b1 will contain dbeta/d\rho_0, dbeta/d\rho_1 etc. So, for example, dbeta_i/d\rho_j
(indices starting at zero) is located in b1[q*j+i].
* b2 will contain d^2beta/d\rho_0d\rho_0, d^2beta/d\rho_1d\rho_0,... but rows will not be
stored if they duplicate an existing row (e.g. d^2beta/d\rho_0d\rho_1 would not be
stored as it already exists and can be accessed by interchanging the sp indices).
So to get d^2beta_k/d\rho_id\rho_j:
i) if i<j interchange the indices
ii) off = (j*m-(j+1)*j/2+i)*q (m is number of sp's)
iii) v2[off+k] is the required derivative.
*/
{ double *Sb,*Skb,*work,*work1,*p1,*p0,*p2,xx;
int i,j,bt,ct,one=1,m,k,rSoff,mk,km;
work = (double *)CALLOC((size_t)*q,sizeof(double));
Sb = (double *)CALLOC((size_t)*q,sizeof(double));
bt=0;ct=0;mgcv_mmult(work,E,beta,&bt,&ct,Enrow,&one,q);
bt=1;ct=0;mgcv_mmult(Sb,E,work,&bt,&ct,q,&one,Enrow); /* S \hat \beta */
for (*bSb=0.0,i=0;i<*q;i++) *bSb += beta[i] * Sb[i]; /* \hat \beta' S \hat \beta */
if (*deriv <=0) {FREE(work);FREE(Sb);return;}
work1 = (double *)CALLOC((size_t)*q,sizeof(double));
Skb = (double *)CALLOC((size_t)*M * *q,sizeof(double));
for (p1=Skb,rSoff=0,i=0;i<*M;i++) { /* first part of first derivatives */
/* form S_k \beta * sp[i]... */
bt=1;ct=0;mgcv_mmult(work,rS + rSoff ,beta,&bt,&ct,rSncol+i,&one,q);
for (j=0;j<rSncol[i];j++) work[j] *= sp[i];
bt=0;ct=0;mgcv_mmult(p1,rS + rSoff ,work,&bt,&ct,q,&one,rSncol+i);
rSoff += *q * rSncol[i];
/* now the first part of the first derivative */
for (xx=0.0,j=0;j<*q;j++,p1++) xx += beta[j] * *p1;
bSb1[i] = xx;
}
if (*deriv>1) for (m=0;m < *M;m++) { /* Hessian */
bt=0;ct=0;mgcv_mmult(work1,E,b1+m * *q,&bt,&ct,Enrow,&one,q);
bt=1;ct=0;mgcv_mmult(work,E,work1,&bt,&ct,q,&one,Enrow); /* S dbeta/drho_m */
for (k=m;k < *M;k++) {
km=k * *M + m;mk=m * *M + k; /* second derivatives needed */
/* d2beta'/drho_k drho_m S beta */
for (xx=0.0,p0=Sb,p1=Sb + *q;p0<p1;p0++,b2++) xx += *b2 * *p0;
bSb2[km] = 2*xx;
/* dbeta'/drho_k S d2beta/drho_m */
for (xx=0.0,p0=b1+k* *q,p1=p0 + *q,p2=work;p0<p1;p0++,p2++) xx += *p2 * *p0;
bSb2[km] += 2*xx;
/* dbeta'/drho_k S_m beta sp[m] */
for (xx=0.0,p0=Skb + k* *q,p1=p0 + *q,p2= b1+m* *q;p0<p1;p0++,p2++) xx += *p2 * *p0;
bSb2[km] += 2*xx;
/* dbeta'/drho_m S_k beta sp[k] */
for (xx=0.0,p0=Skb + m* *q,p1=p0 + *q,p2= b1+k* *q;p0<p1;p0++,p2++) xx += *p2 * *p0;
bSb2[km] += 2*xx;
if (k==m) bSb2[km] += bSb1[k]; else bSb2[mk] = bSb2[km];
}
} /* done Hessian */
/* Now finish off the first derivatives */
bt=1;ct=0;mgcv_mmult(work,b1,Sb,&bt,&ct,M,&one,q);
for (i=0;i<*M;i++) bSb1[i] += 2*work[i];
FREE(Sb);FREE(work);FREE(Skb);FREE(work1);
} /* end get_bSb0 */
void get_bSb(double *bSb,double *bSb1, double *bSb2,double *sp,double *E, void get_bSb(double *bSb,double *bSb1, double *bSb2,double *sp,double *E,
double *rS,int *rSncol,int *Enrow, int *q,int *M,int *M0, double *rS,int *rSncol,int *Enrow, int *q,int *M,int *M0,
...@@ -892,159 +807,6 @@ void get_stableS(double *S,double *Qf,double *sp,double *sqrtS, int *rSncol, int ...@@ -892,159 +807,6 @@ void get_stableS(double *S,double *Qf,double *sp,double *sqrtS, int *rSncol, int
}/* end get_stableS */ }/* end get_stableS */
void get_ddetXWXpS0(double *det1,double *det2,double *P,double *K,double *sp,
double *rS,int *rSncol,double *Tk,double *Tkm,int *n,int *q,int *r,int *M,
int *deriv,int nthreads)
/* obtains derivatives of |X'WX + S| wrt the log smoothing parameters, as required for REML.
The determinant itself has to be obtained during intial decompositions: see gdi().
* P is q by r
* K is n by r
* this routine assumes that sp contains smoothing parameters, rather than log smoothing parameters.
* Note that P and K are as in Wood (2008) JRSSB 70, 495-518.
* uses nthreads via openMP - assumes nthreads already set and nthreads already reset to 1
if openMP not present
*/
{ double *diagKKt,xx,*KtTK,*PtrSm,*PtSP,*trPtSP,*work,*pdKK,*p1,*pTkm;
int m,k,bt,ct,j,one=1,km,mk,*rSoff,deriv2,max_col;
int tid;
#ifdef OMP_REPORT
Rprintf("get_ddetXWXpS0...");
#endif
if (nthreads<1) nthreads = 1;
if (*deriv==2) deriv2=1; else deriv2=0;
/* obtain diag(KK') */
if (*deriv) {
diagKKt = (double *)CALLOC((size_t)*n,sizeof(double));
xx = diagABt(diagKKt,K,K,n,r);
} else { /* nothing to do */
return;
}
/* set up work space */
work = (double *)CALLOC((size_t)*n * nthreads,sizeof(double));
tid=0; /* thread identifier defaults to zero if openMP not available */
/* now loop through the smoothing parameters to create K'TkK */
if (deriv2) {
KtTK = (double *)CALLOC((size_t)(*r * *r * *M),sizeof(double));
#ifdef OPENMP_ON
#pragma omp parallel private(k,j,tid) num_threads(nthreads)
#endif
{ /* open parallel section */
#ifdef OPENMP_ON
#pragma omp for
#endif
for (k=0;k < *M;k++) {
#ifdef OPENMP_ON
tid = omp_get_thread_num(); /* thread running this bit */
#endif
j = k * *r * *r;
getXtWX(KtTK+ j,K,Tk + k * *n,n,r,work + *n * tid);
}
} /* end of parallel section */
} else { KtTK=(double *)NULL;} /* keep compiler happy */
/* start first derivative */
bt=1;ct=0;mgcv_mmult(det1,Tk,diagKKt,&bt,&ct,M,&one,n); /* tr(TkKK') */
/* Finish first derivative and create create P'SmP if second derivs needed */
max_col = *q;
for (j=0;j<*M;j++) if (max_col<rSncol[j]) max_col=rSncol[j]; /* under ML can have q < max(rSncol) */
PtrSm = (double *)CALLOC((size_t)(*r * max_col * nthreads),sizeof(double)); /* storage for P' rSm */
trPtSP = (double *)CALLOC((size_t) *M,sizeof(double));
if (deriv2) {
PtSP = (double *)CALLOC((size_t)(*M * *r * *r ),sizeof(double));
} else { PtSP = (double *) NULL;}
rSoff = (int *)CALLOC((size_t)*M,sizeof(int));
rSoff[0] = 0;for (m=0;m < *M-1;m++) rSoff[m+1] = rSoff[m] + rSncol[m];
tid = 0;
#ifdef OPENMP_ON
#pragma omp parallel private(m,bt,ct,tid) num_threads(nthreads)
#endif
{ /* parallel section start */
#ifdef OPENMP_ON
#pragma omp for
#endif
for (m=0;m < *M;m++) { /* loop through penalty matrices */
#ifdef OPENMP_ON
tid = omp_get_thread_num(); /* thread running this bit */
#endif
bt=1;ct=0;mgcv_mmult(PtrSm + tid * *r * max_col,P,rS+rSoff[m] * *q,&bt,&ct,r,rSncol+m,q);
/*rSoff += rSncol[m];*/
trPtSP[m] = sp[m] * diagABt(work + *n * tid,PtrSm + tid * *r * max_col,
PtrSm + tid * *r * max_col,r,rSncol+m); /* sp[m]*tr(P'S_mP) */
det1[m] += trPtSP[m]; /* completed first derivative */
if (deriv2) { /* get P'S_mP */
bt=0;ct=1;mgcv_mmult(PtSP+ m * *r * *r,PtrSm + tid * *r * max_col,
PtrSm+ tid * *r * max_col ,&bt,&ct,r,r,rSncol+m);
}
}
} /* end of parallel section */
FREE(rSoff);
/* Now accumulate the second derivatives */
// #ifdef OPENMP_ON
//#pragma omp parallel private(m,k,km,mk,xx,tid,pdKK,p1,pTkm) num_threads(nthreads)
//#endif
if (deriv2)
{ /* start of parallel section */
// if (deriv2)
#ifdef OPENMP_ON
#pragma omp parallel for private(m,k,km,mk,xx,tid,pdKK,p1,pTkm) num_threads(nthreads)
#endif
for (m=0;m < *M;m++) {
#ifdef OPENMP_ON
tid = omp_get_thread_num(); /* thread running this bit */
#endif
if (m==0) pTkm = Tkm; else pTkm = Tkm + (m * *M - (m*(m-1))/2) * *n;
for (k=m;k < *M;k++) {
km=k * *M + m;mk=m * *M + k;
/* tr(Tkm KK') */
/*for (xx=0.0,pdKK=diagKKt,p1=pdKK + *n;pdKK<p1;pdKK++,Tkm++) xx += *Tkm * *pdKK;*/
for (xx=0.0,pdKK=diagKKt,p1=pdKK + *n;pdKK<p1;pdKK++,pTkm++) xx += *pTkm * *pdKK;
det2[km] = xx;
/* - tr(KTkKK'TmK) */
det2[km] -= diagABt(work + *n * tid,KtTK + k * *r * *r,KtTK+ m * *r * *r,r,r);
/* sp[k]*tr(P'S_kP) */
if (k==m) det2[km] += trPtSP[m];
/* -sp[m]*tr(K'T_kKP'S_mP) */
det2[km] -= sp[m]*diagABt(work + *n * tid,KtTK + k * *r * *r,PtSP + m * *r * *r,r,r);
/* -sp[k]*tr(K'T_mKP'S_kP) */
det2[km] -= sp[k]*diagABt(work + *n * tid,KtTK + m * *r * *r,PtSP + k * *r * *r,r,r);
/* -sp[m]*sp[k]*tr(P'S_kPP'S_mP) */
det2[km] -= sp[m]*sp[k]*diagABt(work + *n * tid,PtSP + k * *r * *r,PtSP + m * *r * *r,r,r);
det2[mk] = det2[km];
}
}
} /* end of parallel section */
/* free up some memory */
if (deriv2) {FREE(PtSP);FREE(KtTK);}
FREE(diagKKt);FREE(work);
FREE(PtrSm);FREE(trPtSP);
#ifdef OMP_REPORT
Rprintf("done\n");
#endif
} /* end get_ddetXWXpS0 */
void get_ddetXWXpS(double *det1,double *det2,double *P,double *K,double *sp, void get_ddetXWXpS(double *det1,double *det2,double *P,double *K,double *sp,
double *rS,int *rSncol,double *Tk,double *Tkm,int *n,int *q,int *r,int *M,int *M0, double *rS,int *rSncol,double *Tk,double *Tkm,int *n,int *q,int *r,int *M,int *M0,
int *deriv,int nthreads) int *deriv,int nthreads)
......
...@@ -40,7 +40,7 @@ R_CMethodDef CEntries[] = { ...@@ -40,7 +40,7 @@ R_CMethodDef CEntries[] = {
{"mvn_ll", (DL_FUNC) &mvn_ll,15}, {"mvn_ll", (DL_FUNC) &mvn_ll,15},
{"RMonoCon", (DL_FUNC) &RMonoCon, 7}, {"RMonoCon", (DL_FUNC) &RMonoCon, 7},
{"RuniqueCombs", (DL_FUNC) &RuniqueCombs, 4}, {"RuniqueCombs", (DL_FUNC) &RuniqueCombs, 4},
{"RPCLS", (DL_FUNC) &RPCLS, 14}, {"RPCLS", (DL_FUNC) &RPCLS, 13},
{"construct_tprs", (DL_FUNC) &construct_tprs, 13}, {"construct_tprs", (DL_FUNC) &construct_tprs, 13},
{"crspl", (DL_FUNC) &crspl,8}, {"crspl", (DL_FUNC) &crspl,8},