Commit 59a47911 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Debian changes 1.8-9-1

mgcv (1.8-9-1) unstable; urgency=low

  * New upstream release
parents 0159f649 890eb842
** denotes quite substantial/important changes
*** denotes really big changes
1.8-9
* C level fix in bam(...,discrete=TRUE) code. Some memory was mistakenly
allocated via 'calloc' rather than 'R_chk_calloc', but was then freed
via 'R_chk_free'. This could cause R to halt on some platforms.
1.8-8
** New "gp" smooth class (see ?gp.smooth) implemeting the Matern
......
Package: mgcv
Version: 1.8-8
Version: 1.8-9
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with GCV/AIC/REML Smoothness
......@@ -16,6 +16,6 @@ LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
NeedsCompilation: yes
Packaged: 2015-10-24 06:15:25 UTC; sw283
Packaged: 2015-10-30 09:42:04 UTC; sw283
Repository: CRAN
Date/Publication: 2015-10-24 09:17:33
Date/Publication: 2015-10-30 10:53:27
bcb7f62d1e1297b1d685a095b00220aa *ChangeLog
f133d9c24128ba8a977d1674d59f2a23 *DESCRIPTION
2c86df7b3fcdf158d1bb11de61f2a6be *ChangeLog
54d39cfcc0af404fd8bb91703f033f2e *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
b1b86938209b64b673e096da331769eb *NAMESPACE
9cf2059cc5ab5fd2f8a9f4439fec53e5 *R/bam.r
......@@ -70,7 +70,7 @@ b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd
a66a814cc4c6f806e824751fda519ae0 *man/gam2objective.Rd
717401fd7efa3b39d90418a5d1d0c216 *man/gamObject.Rd
a2593990d6a87f7b783d0435062dec02 *man/gamSim.Rd
1a69eb8822eac2d99f92fee8dfe0a818 *man/gamm.Rd
e5e279fdabb71a81f03a16fed5022d8e *man/gamm.Rd
df4a6db696749986fd5e20586fc9b718 *man/gaulss.Rd
39ae109127110152e97dc9f79e08bb14 *man/get.var.Rd
a2ea1233d43fac89e0cacbc09a8d31e2 *man/in.out.Rd
......@@ -87,7 +87,7 @@ b440c31b77c368b484d9665b2b3e89fb *man/jagam.Rd
5169af4be5fccf9fa79b6de08e9ea035 *man/magic.post.proc.Rd
59672d8599211ce9cf44252f92110920 *man/mgcv-FAQ.Rd
38c75ee7bb718c277da5bc649e6886e7 *man/mgcv-package.Rd
9742c9bd1d42aaebe27e2b4218d12672 *man/mgcv-parallel.Rd
f58c9c3919f6971e7f9f7b985a696a08 *man/mgcv-parallel.Rd
00ccf213c31910cd14f1df65a300eb33 *man/model.matrix.gam.Rd
bc9b89db7e7ff246749551c16f5f1f07 *man/mono.con.Rd
d33914a328f645af13f5a42914ca0f35 *man/mroot.Rd
......@@ -123,7 +123,7 @@ c438eb3e41cb1f51e72283380145d210 *man/scat.Rd
76013feaf70d00976bba0154b6f2c946 *man/smooth.construct.cr.smooth.spec.Rd
f5e6d0f5122f61c336827b3615482157 *man/smooth.construct.ds.smooth.spec.Rd
db75c958cbfb561914a3291ab58b9786 *man/smooth.construct.fs.smooth.spec.Rd
3a0545da6f80f00e0b7347f8083a10b3 *man/smooth.construct.gp.smooth.spec.Rd
e38f6ea7f89cb068976b73a177878906 *man/smooth.construct.gp.smooth.spec.Rd
4aaa84b520992fbc32b0c37f7f63c1dd *man/smooth.construct.mrf.smooth.spec.Rd
abe15377f471a2d8957a59c19eeef0bb *man/smooth.construct.ps.smooth.spec.Rd
6aaff8575f68775ed21930893ea9e03d *man/smooth.construct.re.smooth.spec.Rd
......@@ -160,12 +160,12 @@ dc1ef92ff4454734c3a24876e299b760 *po/ko.po
dfd4eec9edc7d1ab6354d47b6e2bd42f *po/pl.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
e71c1a1624b431fbab0a4c8f151d2a97 *src/coxph.c
48f0c9de0da60d48cf4b306c2fdd039a *src/discrete.c
5a8e93ae5f4fae722c1304cecc9c94b3 *src/gdi.c
c1365ce13e8571bdd8dafe90ddb4e6b8 *src/discrete.c
5b0e3daee87ad81384e40ddd2f017dfc *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
d95aa7d646833d1e56f2414f32637bb3 *src/init.c
9a5d7cb3cf93cbdfc08353fbd20a270e *src/magic.c
53cde45313485d2d15e6377461dd3f46 *src/mat.c
52745e0f5fe2fcaf1d3bd9b0297e78f0 *src/mat.c
aae0f298e384952bb8b6b923d40520a8 *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
ec2ae157a9e3bedcccc88b053d0a4e0c *src/mgcv.c
......
mgcv (1.8-9-1) unstable; urgency=low
* New upstream release
-- Dirk Eddelbuettel <edd@debian.org> Sat, 31 Oct 2015 09:29:51 -0500
mgcv (1.8-8-1) unstable; urgency=low
* New upstream release
......
......@@ -127,7 +127,7 @@ Some details of how GAMs are represented as mixed models and estimated using
\code{lme} or \code{gammPQL} in \code{gamm} can be found in Wood (2004 ,2006a,b). In addition \code{gamm} obtains a posterior covariance matrix for the parameters of all the fixed effects and the smooth terms. The approach is similar to that described in Lin & Zhang (1999) - the covariance matrix of the data (or pseudodata in the generalized case) implied by the weights, correlation and random effects structure is obtained, based on the estimates of the parameters of these terms and this is used to obtain the posterior covariance matrix of the fixed and smooth effects.
The bases used to represent smooth terms are the same as those used in \code{\link{gam}}, although adaptive
smoothing bases are not available.
smoothing bases are not available. Prediction from the returned \code{gam} object is straightforward using \code{\link{predict.gam}}, but this will set the random effects to zero. If you want to predict with random effects set to their predicted values then you can adapt the prediction code given in the examples below.
In the event of \code{lme} convergence failures, consider
modifying \code{option(mgcv.vc.logrange)}: reducing it helps to remove
......@@ -329,6 +329,22 @@ b <- gamm(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+
plot(b$gam,pages=1)
summary(b$gam)
vis.gam(b$gam)
## Prediction from gam object, optionally adding
## in random effects.
## Extract random effects and make names more convenient...
refa <- ranef(b$lme,level=5)
rownames(refa) <- substr(rownames(refa),start=9,stop=20)
refb <- ranef(b$lme,level=6)
rownames(refb) <- substr(rownames(refb),start=9,stop=20)
## make a prediction, with random effects zero...
p0 <- predict(b$gam,data.frame(x0=.3,x1=.6,x2=.98,x3=.77))
## add in effect for fa = "2" and fb="2/4"...
p <- p0 + refa["2",1] + refb["2/4",1]
## and a "spatial" example...
library(nlme);set.seed(1);n <- 100
dat <- gamSim(2,n=n,scale=0) ## standard example
......
......@@ -6,6 +6,9 @@
\code{mgcv} can make some use of multiple cores or a cluster. Function \code{\link{bam}} uses the
facilities provided in the \link[parallel]{parallel} package for this purpose. See examples below. Note that most multi-core machines are memory bandwidth limited, so parallel speed up tends to be rather variable.
\code{\link{bam}} can also use an alternative openMP based parallelization approach alongside discretisation of covariates to
achieve substantial speed ups. This is selected using the \code{discrete=TRUE} option to \code{bam}, with the number of threads controlled via the \code{nthreads} argument. See example below.
Function \code{\link{gam}} can use parallel threads on a (shared memory) multi-core
machine via \code{openMP} (where this is supported). To do this, set the desired number of threads by setting \code{nthreads} to the number of cores to use, in the \code{control} argument of \code{\link{gam}}. Note that, for the most part, only the dominant \eqn{O(np^2)}{O(np^2)} steps are parallelized (n is number of data, p number of parameters). For additive Gaussian models estimated by GCV, the speed up can be disappointing as these employ an \eqn{O(p^3)}{O(p^3)} SVD step that can also have substantial cost in practice.
......@@ -17,7 +20,7 @@ Note that on Intel and similar processors the maximum performance is usually ach
mgcv's work splitting always makes the simple assumption that all your cores are equal, and you are not sharing them with other floating point intensive threads.
Note that in addition to hyper-threading several features may lead to apparently poor scaling. The first is that many CPUs have a Turbo mode, whereby a few cores can be run at higher frequency, provided the overall power used by the CPU does not exceed design limits, however it is not possible for all cores on the CPU to run at this frequency. So as you add threads eventually the CPU frequency has to be reduced below the Turbo frequency, with the result that you don't get the expected speed up from adding cores. Secondly, most modern CPUs have their frequency set dynamically according to load. You may need to set the system power management policy to favour high performance in order to maximize the chance that all threads run at the speed you were hoping for (you can turn off dynamic power control in BIOS, but then you turn off the possibility of Turbo also).
In addition to hyper-threading several features may lead to apparently poor scaling. The first is that many CPUs have a Turbo mode, whereby a few cores can be run at higher frequency, provided the overall power used by the CPU does not exceed design limits, however it is not possible for all cores on the CPU to run at this frequency. So as you add threads eventually the CPU frequency has to be reduced below the Turbo frequency, with the result that you don't get the expected speed up from adding cores. Secondly, most modern CPUs have their frequency set dynamically according to load. You may need to set the system power management policy to favour high performance in order to maximize the chance that all threads run at the speed you were hoping for (you can turn off dynamic power control in BIOS, but then you turn off the possibility of Turbo also).
Because the computational burden in \code{mgcv} is all in the linear algebra, then parallel computation may provide reduced or no benefit with a tuned BLAS. This is particularly the case if you are using a multi threaded BLAS, but a BLAS that is tuned to make efficient use of a particular cache size may also experience loss of performance if threads have to share the cache.
......@@ -74,4 +77,9 @@ fv <- predict(b3,cluster=cl) ## parallel prediction
if (!is.null(cl)) stopCluster(cl)
b3
## Alternative using the discrete option with bam...
system.time(b4 <- bam(y ~ s(x0,bs=bs,k=7)+s(x1,bs=bs,k=7)+s(x2,bs=bs,k=k)
,data=dat,family=poisson(),discrete=TRUE,nthreads=2))
}
......@@ -5,9 +5,9 @@
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Low rank Gaussian process smooths}
\description{Gaussian process/kriging models based on a simple covariance functions can be written in a very similar form to thin plate and Duchon spline models (e.g. Handcock, Meier, Nychka, 1994), and low rank versions produced by the eigen approximation method of Wood (2003). Kammann and Wand (2003) suggest a particularly simple form of the Matern covariance function with only a single smoothing parameter to estimate, and this class implements this and other similar models.
\description{Gaussian process/kriging models based on simple covariance functions can be written in a very similar form to thin plate and Duchon spline models (e.g. Handcock, Meier, Nychka, 1994), and low rank versions produced by the eigen approximation method of Wood (2003). Kammann and Wand (2003) suggest a particularly simple form of the Matern covariance function with only a single smoothing parameter to estimate, and this class implements this and other similar models.
Usually invoked by an \code{s(...,bs="ms")} term in a \code{gam} formula. Argument \code{m} of selects the covariance function, sets the range parameter and any power parameter. If \code{m} is not supplied then it defaults to \code{NA} and the covariance function suggested by Kammann and Wand (2003) along with their suggested range parameter is used. Otherwise \code{m[1]} between 1 and 5 selects the correlation function from respectively, spherical, power exponential, and Matern with kappa = 1.5, 2.5 or 3.5. \code{m[2]} if present specifies the range parameter, with non-positive or absent indicating that the Kammann and Wand estimate should be used. \code{m[3]} can be used to specify the power for the power exponential which otherwise defaults to 1.
Usually invoked by an \code{s(...,bs="gp")} term in a \code{gam} formula. Argument \code{m} selects the covariance function, sets the range parameter and any power parameter. If \code{m} is not supplied then it defaults to \code{NA} and the covariance function suggested by Kammann and Wand (2003) along with their suggested range parameter is used. Otherwise \code{m[1]} between 1 and 5 selects the correlation function from respectively, spherical, power exponential, and Matern with kappa = 1.5, 2.5 or 3.5. \code{m[2]} if present specifies the range parameter, with non-positive or absent indicating that the Kammann and Wand estimate should be used. \code{m[3]} can be used to specify the power for the power exponential which otherwise defaults to 1.
}
\usage{
......@@ -33,7 +33,7 @@ smooth class documented under \code{\link{smooth.construct}}, this object will c
avoid any co-linearity problems that might otehrwise occur in the penalty null space basis of the term. }
\item{Xu}{A matrix of the unique covariate combinations for this smooth (the basis is constructed by first stripping
out duplicate locations).}
\item{UZ}{The matrix mapping the smoother parameters back to the parameters of a full Matern spline.}
\item{UZ}{The matrix mapping the smoother parameters back to the parameters of a full GP smooth.}
\item{null.space.dimension}{The dimension of the space of functions that have zero wiggliness according to the wiggliness penalty for this term.}
\item{gp.defn}{the type, range parameter and power parameter defining the correlation function. }
}
......
File mode changed from 100644 to 100755
......@@ -135,7 +135,7 @@ void tensorXb(double *f,double *X, double *C,double *work, double *beta,
no constraint is applied.
*/
char trans='N';
int pb=1,md,*kp,*kd,pd,i,j,one=1;
int pb=1,md,*kp,*kd,pd,i,j;
double *M,done=1.0,dzero=0.0,*p0,*p1,*p2,*p3,*pf,*pc,x;
M = X;
for (i=0;i<*dt-1;i++) {
......@@ -182,7 +182,7 @@ void Xbd(double *f,double *beta,double *X,int *k, int *m,int *p, int *n,
int *nx, int *ts, int *dt, int *nt,double *v,int *qc) {
/* Forms f = X beta for X stored in the packed form described in function XWX */
int i,j,q,*pt,*off,*voff,*tps,dC=0,c1;
double *f0,*pf,*p0,*p1,*p2,*C,*work,maxp=0;
double *f0,*pf,*p0,*p1,*p2,*C=NULL,*work,maxp=0;
/* obtain various indices */
pt = (int *) R_chk_calloc((size_t)*nt,sizeof(int)); /* the term dimensions */
off = (int *) R_chk_calloc((size_t)*nx+1,sizeof(int)); /* offsets for X submatrix starts */
......@@ -196,7 +196,6 @@ void Xbd(double *f,double *beta,double *X,int *k, int *m,int *p, int *n,
if (c1>dC) dC = c1; /* dimension of working matrix C */
}
if (j==0) pt[i] = p[q]; else pt[i] *= p[q]; /* term dimension */
//if (maxm<m[q]) maxm=m[q];
}
if (qc[i]>0) voff[i+1] = voff[i] + pt[i]; else voff[i+1] = voff[i]; /* start of ith v matrix */
if (maxp < pt[i]) maxp = pt[i];
......@@ -204,10 +203,10 @@ void Xbd(double *f,double *beta,double *X,int *k, int *m,int *p, int *n,
else tps[i+1] = tps[i] + pt[i] - 1; /* there is a tensor constraint to apply - reducing param count*/
}
/* now form the product term by term... */
pf=f0 = (double *)calloc((size_t)*n,sizeof(double));
pf=f0 = (double *)R_chk_calloc((size_t)*n,sizeof(double));
i = *n; if (i<maxp) i=maxp;
work = (double *)calloc((size_t)i,sizeof(double));
if (dC) C = (double *)calloc((size_t)dC,sizeof(double));
work = (double *)R_chk_calloc((size_t)i,sizeof(double));
if (dC) C = (double *)R_chk_calloc((size_t)dC,sizeof(double));
for (i=0;i < *nt;i++) { /* work through terms */
if (i==0) f0 = f; /* result written straight to f for i==0 */
if (dt[i]==1) singleXb(f0,work,X+off[ts[i]],beta+tps[i],k + *n *ts[i],m+ts[i], p+ts[i],n);
......@@ -226,9 +225,8 @@ void Xbd(double *f,double *beta,double *X,int *k, int *m,int *p, int *n,
void XWyd(double *XWy,double *y,double *X,double *w,int *k, int *m,int *p, int *n,
int *nx, int *ts, int *dt, int *nt,double *v,int *qc,
int *ar_stop,int *ar_row,double *ar_weights) {
double *Wy,*p0,*p1,*p2,*p3,done=1.0,dzero=0.0,*Xy0,*work,*work1,x;
double *Wy,*p0,*p1,*p2,*p3,*Xy0,*work,*work1,x;
int q,i,j,*pt,*off,*voff,*tps,maxm=0,maxp=0,one=1,zero=0;
char trans='T';
if (*ar_stop>=0) { /* model has AR component, requiring sqrt(weights) */
for (p0 = w,p1 = w + *n;p0<p1;p0++) *p0 = sqrt(*p0);
}
......@@ -268,7 +266,6 @@ void XWyd(double *XWy,double *y,double *X,double *w,int *k, int *m,int *p, int *
for (x=0.0,p0=Xy0,p1=p0 + pt[i],p2=v+voff[i];p0<p1;p0++,p2++) x += *p0 * *p2; /* x = v'Xy0 */
p0=XWy + tps[i];p1 = p0 + pt[i]-1;p2 = v+voff[i] + 1;p3=Xy0+1;
for (;p0<p1;p0++,p2++,p3++) *p0 = *p3 - x * *p2; /* (I-vv')Xy0 less first element */
//F77_CALL(dgemv)(&trans, pt+i, qc+i,&done,Q+Qoff[i],pt+i,Xy0,&one,&dzero,XWy + tps[i],&one);
} else { /* straight copy */
for (p0=Xy0,p1=p0+pt[i],p2=XWy+tps[i];p0<p1;p0++,p2++) *p2 = *p0;
}
......@@ -304,10 +301,9 @@ void XWXd(double *XWX,double *X,double *w,int *k, int *m,int *p, int *n, int *nx
AR models are handled via the 3 ar_* arrays. ar_stop[0] < 0 signals no AR.
*/
int r,c,i,j,q,*pt,*pd,*off,a,b,*tps,ptot,maxp=0,maxm=0,*voff,pa,pb,kk,dk,rk,*start,one=1,zero=0;
double done=1.0,dzero=0.0,*p0,*p1,*p2, *Xi,*temp,*tempn,*xwx,*xwx0,
int r,c,i,j,q,*pt,*pd,*off,a,b,*tps,ptot,maxp=0,maxm=0,*voff,pa,pb,kk,dk,*start,one=1,zero=0;
double *p0,*p1,*p2, *Xi,*temp,*tempn,*xwx,*xwx0,
*XiB,*tempB,*tempnB,*x0,*x1,x;
char trans='T',not_trans='N';
#ifndef SUPPORT_OPENMP
*nthreads = 1;
#endif
......@@ -346,7 +342,7 @@ void XWXd(double *XWX,double *X,double *w,int *k, int *m,int *p, int *n, int *nx
a=c;b=r;
}
/* split cols between threads... */
dk = pt[b] / *nthreads; rk = pt[b] % *nthreads;
dk = pt[b] / *nthreads; //rk = pt[b] % *nthreads;
if (dk * *nthreads < pt[b]) dk++;start[0]=0;
for (i=0;i<*nthreads;i++) {
start[i+1] = start[i] + dk;
......@@ -394,12 +390,6 @@ void XWXd(double *XWX,double *X,double *w,int *k, int *m,int *p, int *n, int *nx
/* if Xb is tensor, may need to apply constraint */
if (dt[a]>1&&qc[a]>0) { /* first term is a tensor with a constraint */
///* copy xwx to xwx0 */
//for (p0=xwx,p1=p0 + pt[b]*pt[a],p2=xwx0;p0<p1;p0++,p2++) *p2 = *p0;
/* apply Qa' */
//F77_CALL(dgemm)(&trans,&not_trans,qc+a,pt+b,pt+a, &done,
// Q+Qoff[a],pt+a,xwx0,pt+a,&dzero,xwx,qc+a);
//pa = qc[a]; /* rows of xwx */
x0=x1=xwx; /* pointers to columns of xwx */
/* col by col form (I-vv')xwx, dropping first row... */
for (j=0;j<pt[b];j++) {
......@@ -418,10 +408,6 @@ void XWXd(double *XWX,double *X,double *w,int *k, int *m,int *p, int *n, int *nx
for (p0=xwx0+j+pa,p1=v+voff[b],p2=p1+pt[b],p1++;p1<p2;x1+=pa,p1++,p0+=pa) *x1 = *p0 - *p1 *x;
}
pb=pt[b]-1;
///* apply Qa from right */
//F77_CALL(dgemm)(&not_trans,&not_trans,&pa,qc+b,pt+b, &done,
// xwx0,&pa, Q+Qoff[b],pt+b,&dzero,xwx,&pa);
//pb = qc[b]; /* rows of xwx */
} else pb=pt[b];
/* copy result into overall XWX*/
......@@ -436,6 +422,5 @@ void XWXd(double *XWX,double *X,double *w,int *k, int *m,int *p, int *n, int *nx
R_chk_free(start);
R_chk_free(XiB); R_chk_free(tempnB); R_chk_free(pt); R_chk_free(pd); R_chk_free(off);
R_chk_free(voff); R_chk_free(tps); R_chk_free(tempB); R_chk_free(xwx); R_chk_free(xwx0);
//R_chk_free(); R_chk_free(); R_chk_free(); R_chk_free(); R_chk_free();
} /* XWX */
......@@ -1947,7 +1947,7 @@ void gdiPK(double *work,double *X,double *E,double *Es,double *rS,double *U1,dou
z must contain Wz on entry.
*/
{ int i,j,k,*pivot,nt1,nr,left,tp,bt,ct,TRUE=1,FALSE=0,one=1;
double *zz,*WX,*tau,*R1,Rnorm,Enorm,Rcond,*Q,*tau1,*Ri,ldetI2D,*IQ,*d,*p0,*p1,*p2,*p3,*p4,xx,norm1,norm2;
double *zz=NULL,*WX,*tau,*R1,Rnorm,Enorm,Rcond,*Q=NULL,*tau1,*Ri,ldetI2D,*IQ,*d,*p0,*p1,*p2,*p3,*p4,xx,norm1,norm2;
nt1 = *nt;
if (*type==0) {
......
......@@ -847,7 +847,7 @@ int mgcv_piqr(double *x,int n, int p, double *beta, int *piv, int nt) {
#pragma omp parallel for private(i,j,p1,v,z,z1,zz,ii) num_threads(nt)
#endif
for (i=0;i<nth;i++) {
if (i == nth-1) j = cpf;
if (i == nth-1) j = cpf; else j = cpt;
p1 = p0 + n + n * cpt *i;
z1=p1+nh;
for (ii =0;ii<j;ii++,p1+=n,z1+=n) {
......
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