Commit a1ae9b9e authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-22

parent 1b84d0b8
......@@ -4,6 +4,28 @@
Currently deprecated and liable to be removed:
- 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
* bam(,discrete=TRUE) could produce garbage with ti(x,z,k=c(6,5),mc=c(T,F))
......
Package: mgcv
Version: 1.8-20
Version: 1.8-22
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-09-09 07:02:31 UTC; sw283
Packaged: 2017-09-18 10:38:41 UTC; sw283
Repository: CRAN
Date/Publication: 2017-09-09 16:26:09 UTC
Date/Publication: 2017-09-19 00:27:56 UTC
c26bec397cf6ee2a88d4c295945b8cab *ChangeLog
ac713195702fe496609fec839dbdcc9f *DESCRIPTION
bc8bc8be9cec065f0c7ec0cf670516f2 *ChangeLog
ad2a852c227dd80d86b24784d12dba42 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
3d1bd78b1a6f876c1160b7e64ea79cc0 *NAMESPACE
043efd04d042fba89ddd6781c0220ab9 *R/bam.r
9fa89dc9361930dca536d4e011e4e605 *NAMESPACE
b051474f30e8e779f2440e8ecd5bbd51 *R/bam.r
7b419683b0948cf6da009f614078fe90 *R/coxph.r
777a0d67a1f7fa14bf87bc312064061b *R/efam.r
f610b2695c525c2d9d8f8174e6760548 *R/fast-REML.r
4dc75c18c9feceb73d207b8ac01b2f5f *R/gam.fit3.r
60aa872ec8901e4dc6273e653f77c83f *R/gam.fit4.r
dfdb821247da3e780de0d4599b88735d *R/fast-REML.r
0a9cc97524994d1563b93769284266c4 *R/gam.fit3.r
f0fd93b43b36abf05e6666b4a71c20fa *R/gam.fit4.r
1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r
b99d9de0c46a230e081b9bbbbbb3c617 *R/gamlss.r
c934ba653cb1a904132b3f15a01e41c5 *R/gamm.r
10facb791e4cfd123d183f05660119c6 *R/jagam.r
48c3f7080c31084058aaf3de06015e06 *R/mgcv.r
b9084037d058ef5d787e2818fbc29bae *R/mgcv.r
2feca0dc9d354de7bc707c67a5891364 *R/misc.r
03772998ab05887f2eeae10dd6efe983 *R/mvam.r
ca13f63d1112d31258eb3c68464fa653 *R/plots.r
......@@ -42,7 +42,7 @@ fd0cfd64be579f9fbecdbb7f2b8ec1eb *man/Sl.initial.repara.Rd
60670020425f8749b81a8d8c3f168880 *man/Sl.setup.Rd
69ae63833438a3af2963e39482f1d72f *man/Tweedie.Rd
8087ab00d10b44c99c37f49bf90e19cd *man/anova.gam.Rd
a7e0ce83164f1e34d0a9d99d0f9853f3 *man/bam.Rd
6180a4e9ea206e1a350b993dade0a869 *man/bam.Rd
ab5e37c3bf8803de63b63c3bdc5909cd *man/bam.update.Rd
cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd
745cbf31eb14fc1c5916fc634c74d998 *man/bug.reports.mgcv.Rd
......@@ -151,7 +151,7 @@ b45d8e71bda4ceca8203dffea577e441 *man/smooth.construct.re.smooth.spec.Rd
0bfe981f2c3e6ea5b8d5372076ccde53 *man/smooth.construct.sos.smooth.spec.Rd
3cb4e59f915c8d64b90754eaeeb5a86f *man/smooth.construct.t2.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
4c49358e5e6a70d1eca69a2ccaa77609 *man/smooth2random.Rd
844f9653d74441293d05a24dd3e2876a *man/smoothCon.Rd
......@@ -185,19 +185,19 @@ ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
342aa30c8f6f1e99ffa2576a6f29d7ce *src/coxph.c
0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c
dba99f7d7cc412dd9255f6307c8c7fa7 *src/gdi.c
f6c4ba80ce1b71be7f4a44fac1b08c28 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
6ea9eff0b4b804a4e3ebf830c0ecacb9 *src/init.c
a287f92033406194afe0301ee076fda8 *src/init.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c
654ff83187dc0f7ef4e085f3348f70d2 *src/mat.c
0545dabf3524a110d616ea5e6373295d *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
48cde0e19d5dd54b131ba66c777c0ec2 *src/mgcv.c
f2eacd39b434ddef7a3bf9d53fb6fe77 *src/mgcv.h
e4cef7f1753153fbab242d1c4d4f7e7f *src/matrix.c
de37b0972199b796654405efc007f25b *src/matrix.h
8df4b96961491d76989b50856237ee2d *src/mgcv.c
b6278c9a4a32bdc16487a96e6a0045b6 *src/mgcv.h
97e3717e95a70b1470b4c3071e144d17 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
8f480dc455f9ff011c3e8f059efec2c5 *src/qp.c
cd563899be5b09897d1bf36a7889caa0 *src/qp.h
563938b7bb6504ab10df5376c4360220 *src/qp.c
073a4b5b0bc6e869c5b35478c47facf1 *src/qp.h
d5673b88f6f3d85c62a1337f49abba24 *src/soap.c
dcac8c02b5f9be28d13efc834fc88d55 *src/sparse-smooth.c
fe0444bece322bc229e46b3d1c150779 *src/tprs.c
......
......@@ -85,7 +85,7 @@ influence,logLik,lm,mad,
make.link,median,model.frame,model.offset,model.matrix,nlm,
na.pass,napredict,na.omit,naresid,optim,pchisq,pnorm,pt,pf,
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,
rbeta,rbinom,rgamma,rnbinom,rnorm,rpois,runif,sd,
termplot,terms.formula,terms,uniroot,var,vcov,weights)
......
......@@ -109,8 +109,8 @@ qr.up <- function(arg) {
wt <- c(wt,w)
w <- sqrt(w)
## 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)
else qrx <- qr.update(w*X[good,],w*z,qrx$R,qrx$f,qrx$y.norm2,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,,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
}
qrx$dev <- dev;qrx$wt <- wt;qrx$eta <- eta
......@@ -191,8 +191,8 @@ check.term <- function(term,rec) {
ii <- which(rec$vnames%in%term)
if (length(ii)) { ## at least one variable already discretized
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")
else return(rec$ki[min(ii)]) ## all names match previous - retun index of previous
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 - return index of previous
} else stop("bam can not discretize with this nesting structure")
} else return(0) ## no match
} ## check.term
......@@ -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)
w <- sqrt(w)
## 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)
else qrx <- qr.update(w*X[good,],w*z,qrx$R,qrx$f,qrx$y.norm2,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,,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
}
if (use.chol) { ## post proc to get R and f...
......@@ -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,
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,
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,...)
## 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
warning("discretization only available with fREML")
} else {
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")) {
......
......@@ -474,7 +474,7 @@ Sl.addS <- function(Sl,A,rho) {
ind <- (Sl[[b]]$start:Sl[[b]]$stop)[Sl[[b]]$ind]
if (length(Sl[[b]]$S)==1) { ## singleton
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
k <- k + 1
} else {
......
This diff is collapsed.
......@@ -570,27 +570,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
rSncol=as.integer(rSncol),deriv=as.integer(deriv),
fixed.penalty = as.integer(rp$fixed.penalty),nt=as.integer(control$nthreads),
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 <- T %*% rV
## derivatives of coefs w.r.t. sps etc...
......
......@@ -97,8 +97,7 @@ pcls <- function(M)
# M$sp - array of theta_i's
# 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)
H<-0
{ nar<-c(length(M$y),length(M$p),dim(M$Ain)[1],dim(M$C)[1])
## sanity checking ...
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)")
......@@ -140,13 +139,12 @@ pcls <- function(M)
M$C <- matrix(0,0,0)
M$p <- rep(0,ncol(M$X))
nar[2] <- length(M$p)
nar[4] <- 0
qra.exist <- TRUE
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)
,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))
p <- array(o[[2]],length(M$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
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
if (is.null(block.size)) {
## use one block as predicting using model frame
......
......@@ -21,7 +21,7 @@ bam(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,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,
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,...)
}
%- maybe also `usage' for other objects documented here.
......@@ -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).
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
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 @@
\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
regression splines. \code{m} specifies the order of the derivatives in the thin plate spline penalty. 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, for example.
regression splines. \code{m} specifies the order of the derivatives in the thin plate spline penalty.
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
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
\examples{
require(mgcv); n <- 100; set.seed(2)
x <- runif(n); y <- x + x^2*.2 + rnorm(n) *.1
## is smooth significantly different from straight line?
summary(gam(y~s(x,m=c(2,0))+x,method="REML")) ## not quite
## is smooth significatly different from zero?
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...
k <- 10; m <- 2
y <- y[order(x)];x <- x[order(x)]
......
......@@ -95,91 +95,6 @@ double trAB(double *A,double *B,int *n, int *m)
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,
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
}/* 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,
double *rS,int *rSncol,double *Tk,double *Tkm,int *n,int *q,int *r,int *M,int *M0,
int *deriv,int nthreads)
......
......@@ -40,7 +40,7 @@ R_CMethodDef CEntries[] = {
{"mvn_ll", (DL_FUNC) &mvn_ll,15},
{"RMonoCon", (DL_FUNC) &RMonoCon, 7},
{"RuniqueCombs", (DL_FUNC) &RuniqueCombs, 4},
{"RPCLS", (DL_FUNC) &RPCLS, 14},
{"RPCLS", (DL_FUNC) &RPCLS, 13},
{"construct_tprs", (DL_FUNC) &construct_tprs, 13},
{"crspl", (DL_FUNC) &crspl,8},
{"predict_tprs", (DL_FUNC) &predict_tprs, 12},
......
This diff is collapsed.
......@@ -15,16 +15,12 @@ extern long matrallocd;
/* The user routines */
void mtest(void);
void mcopy(matrix *A,matrix *B);
void sort(matrix);
matrix initmat(int rows,int cols);
void freemat(matrix A);
void vmult(matrix *A,matrix *b,matrix *c,int t);
void matmult(matrix C,matrix A,matrix B,int tA,int tB);
void multi(int n,matrix C, ...);
void invert(matrix *a);
void tricholeski(matrix *T,matrix *l0,matrix *l1);
double dot(matrix a,matrix b);
double enorm(matrix d);
void householder(matrix *u,matrix a,matrix b,int t1);
......@@ -33,20 +29,11 @@ void HQmult(matrix C,matrix U,int p,int t);
void QT(matrix Q,matrix A,int Qfull);
void Rsolv(matrix *R,matrix *p,matrix *y, int transpose);
int QR(matrix *Q,matrix *R);
void UTU(matrix *T,matrix *U);
void bidiag(matrix *A,matrix *wl,matrix *ws,matrix *V);
void OrthoMult(matrix *Q,matrix *A,int off,int rows,int t,int pre,int o_pre);
void root(matrix *M,matrix *C,double tol);
void matrixintegritycheck(void);
double mean(matrix a);
double pythag(double a,double b);
void svd(matrix *a,matrix *w,matrix *v);
matrix svdroot(matrix A,double reltol);
void svd_bidiag(matrix *U, matrix *w, matrix *ws,matrix *V);
void msort(matrix a);
void RArrayFromMatrix(double *a,int r,matrix *M);
matrix Rmatrix(double *A,int r,int c);
matrix initvec(int rows);
void interchange(matrix *M,int i,int j,int col);
#endif
......@@ -410,7 +410,7 @@ void RMonoCon(double *Ad,double *bd,double *xd,int *control,double *lower,double
void RPCLS(double *Xd,double *pd,double *yd, double *wd,double *Aind,double *bd,
double *Afd,double *Hd,double *Sd,
double *Afd,double *Sd,
int *off,int *dim,double *theta, int *m,int *nar)
/* Interface routine for PCLS the constrained penalized weighted least squares solver.
......@@ -419,7 +419,7 @@ void RPCLS(double *Xd,double *pd,double *yd, double *wd,double *Aind,double *bd
np=nar[1] - number of parameters
nai=nar[2] - number of inequality constraints
naf=nar[3] - number of fixed constraints
getH=nar[4] - 0 for no hat matrix, 1 to produce one.
Problem to be solved is:
......@@ -438,7 +438,7 @@ void RPCLS(double *Xd,double *pd,double *yd, double *wd,double *Aind,double *bd
on exit p contains the best fit parameter vector.
*/
{ matrix y,X,p,w,Ain,Af,b,H,*S;
{ matrix y,X,p,w,Ain,Af,b,*S;
int n,np,i,*active;
np=nar[1];n=nar[0];
......@@ -452,20 +452,20 @@ void RPCLS(double *Xd,double *pd,double *yd, double *wd,double *Aind,double *bd
if (nar[2]>0) b=Rmatrix(bd,(long)nar[2],1L);else b.r=0L;
if (*m) S=(matrix *)CALLOC((size_t) *m,sizeof(matrix));
else S=&H; /* avoid spurious compiler warning */
else S=NULL; /* avoid spurious compiler warning */
for (i=0;i< *m;i++) S[i]=initmat((long)dim[i],(long)dim[i]);
RUnpackSarray(*m,S,Sd);
if (nar[4]) H=initmat(y.r,y.r); else H.r=H.c=0L;
//if (nar[4]) H=initmat(y.r,y.r); else H.r=H.c=0L;
active=(int *)CALLOC((size_t)(p.r+1),sizeof(int)); /* array for active constraints at best fit active[0] will be number of them */
/* call routine that actually does the work */
PCLS(&X,&p,&y,&w,&Ain,&b,&Af,&H,S,off,theta,*m,active);
PCLS(&X,&p,&y,&w,&Ain,&b,&Af,S,off,theta,*m,active);
/* copy results back into R arrays */
for (i=0;i<p.r;i++) pd[i]=p.V[i];
if (H.r) RArrayFromMatrix(Hd,H.r,&H);
//if (H.r) RArrayFromMatrix(Hd,H.r,&H);
/* clear up .... */
FREE(active);
......@@ -473,7 +473,7 @@ void RPCLS(double *Xd,double *pd,double *yd, double *wd,double *Aind,double *bd
if (*m) FREE(S);
freemat(X);freemat(p);freemat(y);freemat(w);
if (H.r) freemat(H);
//if (H.r) freemat(H);
if (Ain.r) freemat(Ain);
if (Af.r) freemat(Af);
if (b.r) freemat(b);
......
......@@ -102,7 +102,7 @@ void rwMatrix(int *stop,int *row,double *w,double *X,int *n,int *p,int *trans,do
void in_out(double *bx, double *by, double *break_code, double *x,double *y,int *in, int *nb, int *n);
void Rlanczos(double *A,double *U,double *D,int *n, int *m, int *lm,double *tol,int *nt);
void RuniqueCombs(double *X,int *ind,int *r, int *c);
void RPCLS(double *Xd,double *pd,double *yd, double *wd,double *Aind,double *bd,double *Afd,double *Hd,double *Sd,int *off,int *dim,double *theta, int *m,int *nar);
void RPCLS(double *Xd,double *pd,double *yd, double *wd,double *Aind,double *bd,double *Afd,double *Sd,int *off,int *dim,double *theta, int *m,int *nar);
void RMonoCon(double *Ad,double *bd,double *xd,int *control,double *lower,double *upper,int *n);
/*void MinimumSeparation(double *gx,double *gy,int *gn,double *dx,double *dy, int *dn,double *dist);*/
void MinimumSeparation(double *x,int *n, int *d,double *t,int *m,double *dist);
......
......@@ -516,7 +516,7 @@ void QPCLS(matrix *Z,matrix *X, matrix *p, matrix *y,matrix *Ain,matrix *b,matri
void PCLS(matrix *X,matrix *p,matrix *y,matrix *w,matrix *Ain,matrix *b,
matrix *Af,matrix *H,matrix *S,int *off,double *theta,int m,int *active)
matrix *Af,matrix *S,int *off,double *theta,int m,int *active)
/* Routine for Penalized Constrained Least Squares problems.
PCLS() is an interface routine for QPCLS for solving the general problem class:
......@@ -538,9 +538,6 @@ void PCLS(matrix *X,matrix *p,matrix *y,matrix *w,matrix *Ain,matrix *b,
... where F = [ X'W^0.5, B^0.5']' and z = [y'W^0.5, 0]'. This rewrite is
performed and then QPCLS is called to obtain the solution.
If H->r==y->r on entry, then an influence (or "hat") matrix is returned in H.
At present the calculation of H is inefficient and none too stable.
On exit active[] contains a list of the active inequlity constraints in elements
1->active[0]. This array should be initialized to length p.r+1 on entry.
......@@ -548,9 +545,9 @@ void PCLS(matrix *X,matrix *p,matrix *y,matrix *w,matrix *Ain,matrix *b,