Commit 93a4e91b authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

New upstream version 1.8-26

parent d6b7d412
......@@ -17,6 +17,16 @@ present, and reparameterization needs checking (also for bam).
* OpenBLAS 0.3.2/3 appears not to be thread safe at present - can
cause problems when opemmp multithreading.
1.8-26
* LINPACK dependency removed.
* Added service routine chol.down to down date a Cholesky factor on row/col
deletion.
* liu2 had a check messed up when vectorized. Fix to stop vector being
checked for equality to zero.
1.8-25
** bam(...,discrete=TRUE) methods improved. Cross products now usually faster
......@@ -28,7 +38,7 @@ present, and reparameterization needs checking (also for bam).
sarse matrices). Useful for other packages using constructors and
smooth2random, for 'fs' smooths.
* The mrf smooth constructor contained an obsolete hack in which the term
* The mrf smooth constructor contained a obsolete hack in which the term
dimension was set to 2 to avoid plotting when used as a te marginal. This
messed up side constraints for terms where a mrf smooth was a main effect
and te marginal. Fixed.
......
Package: mgcv
Version: 1.8-25
Version: 1.8-26
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: 2018-10-26 08:31:24 UTC; sw283
Packaged: 2018-11-20 20:40:48 UTC; sw283
Repository: CRAN
Date/Publication: 2018-10-26 14:50:06 UTC
Date/Publication: 2018-11-21 11:50:02 UTC
7471749e7e75a48aba7d56c7d0451464 *ChangeLog
28c24bc01e4b54e0368fc857108a6efb *DESCRIPTION
27ebfd695d348a41dde36452454d4f80 *ChangeLog
8c7c5ba5650690f9409b109c928718b1 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
baf2b1ee09bd0ba3b653ab4247b15ce4 *NAMESPACE
6fbddaad9c9fad126ce709cbd09e56d3 *NAMESPACE
459045bffc77934a4f287258f31f8aa6 *R/bam.r
de69a45fe00e28a0309193df283c5829 *R/coxph.r
6c79693fe31558339cd694883d0beeb1 *R/efam.r
......@@ -12,8 +12,8 @@ a00617b7da9f9de5fff27f13b5086c76 *R/gam.fit3.r
f2368e22ccf048cf3d1e210e970013dd *R/gamlss.r
ce69039c15d4db9cfda5d1309af8b242 *R/gamm.r
10facb791e4cfd123d183f05660119c6 *R/jagam.r
70308e442854adfe95980e3a6c05357f *R/mgcv.r
2cea3312b32a86d53e0304b1279da048 *R/misc.r
081f37414ae038f5608a1a15e34bb34a *R/mgcv.r
a186527989eee0303f966f80a5fd7b84 *R/misc.r
03772998ab05887f2eeae10dd6efe983 *R/mvam.r
eed19ceca4c4114023f12cda3f622659 *R/plots.r
530e432792d072e2587cbd0896b620cf *R/smooth.r
......@@ -42,11 +42,12 @@ fd0cfd64be579f9fbecdbb7f2b8ec1eb *man/Sl.initial.repara.Rd
60670020425f8749b81a8d8c3f168880 *man/Sl.setup.Rd
69ae63833438a3af2963e39482f1d72f *man/Tweedie.Rd
8087ab00d10b44c99c37f49bf90e19cd *man/anova.gam.Rd
d235547d7f960f11c6ca2cb0ed571ef6 *man/bam.Rd
ca103588880287accfbde4265cc56f54 *man/bam.Rd
ab5e37c3bf8803de63b63c3bdc5909cd *man/bam.update.Rd
cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd
745cbf31eb14fc1c5916fc634c74d998 *man/bug.reports.mgcv.Rd
530b8b4bacffa9353561f19ecfccfe19 *man/cSplineDes.Rd
05d2b5d4b20621a356a4dff2f88c3fb5 *man/chol.down.Rd
1e4a88fca144eda0f985b362ff5b3972 *man/choose.k.Rd
c03748964ef606621418e428ae49b103 *man/columb.Rd
9906a83ce29a3b17044bc2a3c9940cee *man/concurvity.Rd
......@@ -125,7 +126,7 @@ ee9352ba4c531a8def16deddcab9a9fd *man/pdIdnot.Rd
8280e8aff6d5a004449d9e19661e3920 *man/plot.gam.Rd
c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
1a9d83c9fc67e5f0fc85d66d3112f4ef *man/predict.bam.Rd
40056c77f064706e6933de405d2339be *man/predict.gam.Rd
650cd7383fe7d54920e2e3273745780e *man/predict.gam.Rd
41b1b0883f2419c8e2d2dd913fc7d677 *man/print.gam.Rd
6d0ce4e574fabceffdbedd46c91364cb *man/qq.gam.Rd
22b7dcbc8ff4096365fa98ce56b957c9 *man/rTweedie.Rd
......@@ -188,13 +189,13 @@ fdbb7d250e7baa7d73f26837246480f3 *po/pl.po
4801447cdbc0638d0aaf8ff723ba4b6b *src/discrete.c
0aa17445bd208f941d79e5dd70181297 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
4bd936321fe26c1c77cb810f7f627258 *src/init.c
818b7134dc51a935b07f1f2792757635 *src/init.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c
81a300809799304235bd032949bc047a *src/mat.c
6c0695593ed52375475c9bd15493d38a *src/mat.c
e4cef7f1753153fbab242d1c4d4f7e7f *src/matrix.c
de37b0972199b796654405efc007f25b *src/matrix.h
8df4b96961491d76989b50856237ee2d *src/mgcv.c
3b7f9c7eff3a15b27c4cc5bd475afe6d *src/mgcv.h
3cc1f17f3c2f04e0ba21639453959038 *src/mgcv.h
8e3efe7a0d6b9619671cf2f7471f2112 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
563938b7bb6504ab10df5376c4360220 *src/qp.c
......
useDynLib(mgcv, .registration = TRUE, .fixes = "C_")
export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
export(anova.gam, bam, bam.update,bandchol, betar,
choldrop,cox.ph,concurvity,
cSplineDes,dDeta,
exclude.too.far,extract.lme.cov, extract.lme.cov2,
formXtViX, full.score, formula.gam,fixDependence,fix.family.link,
......
......@@ -3208,14 +3208,16 @@ liu2 <- function(x, lambda, h = rep(1,length(lambda)),lower.tail=FALSE) {
lh <- lh*lambda
c3 <- sum(lh)
if (x<=0 || c2 <= 0) return(1)
xpos <- x > 0
res <- 1 + 0 * x
if (sum(xpos)==0 || c2 <= 0) return(res)
s1 <- c3/c2^1.5
s2 <- sum(lh*lambda)/c2^2
sigQ <- sqrt(2*c2)
t <- (x-muQ)/sigQ
t <- (x[xpos]-muQ)/sigQ
if (s1^2>s2) {
a <- 1/(s1-sqrt(s1^2-s2))
......@@ -3224,16 +3226,15 @@ liu2 <- function(x, lambda, h = rep(1,length(lambda)),lower.tail=FALSE) {
} else {
a <- 1/s1
delta <- 0
if (c3==0) return(1)
if (c3==0) return(res)
l <- c2^3/c3^2
}
muX <- l+delta
sigX <- sqrt(2)*a
return(pchisq(t*sigX+muX,df=l,ncp=delta,lower.tail=lower.tail))
}
res[xpos] <- pchisq(t*sigX+muX,df=l,ncp=delta,lower.tail=lower.tail)
res
} ## liu2
simf <- function(x,a,df,nq=50) {
## suppose T = sum(a_i \chi^2_1)/(chi^2_df/df). We need
......
......@@ -227,6 +227,21 @@ dchol <- function(dA,R) {
return(matrix(oo$dR,p,p))
} ## dchol
choldrop <- function(R,k) {
## routine to update Cholesky factor R of A on dropping row/col k of A.
## R can be upper triangular, in which case (R'R=A) or lower triangular in
## which case RR'=A...
n <- as.integer(ncol(R))
k1 <- as.integer(k-1)
ut <- as.integer(as.numeric(R[1,2]!=0))
if (k<1||k>n) return(R)
Rup <- matrix(0,n-1,n-1)
#oo <- .C(C_chol_down,R=as.double(R),Rup=as.double(Rup),n=as.integer(n),k=as.integer(k-1),ut=as.integer(ut))
.Call(C_mgcv_chol_down,R,Rup,n,k1,ut)
#matrix(oo$Rup,n-1,n-1)
Rup
} ## choldrop
vcorr <- function(dR,Vr,trans=TRUE) {
## Suppose b = sum_k op(dR[[k]])%*%z*r_k, z ~ N(0,Ip), r ~ N(0,Vr). vcorr returns cov(b).
## dR is a list of p by p matrices. 'op' is 't' if trans=TRUE and I() otherwise.
......
......@@ -200,9 +200,10 @@ The extended families given in \code{\link{family.mgcv}} can also be used. The e
\references{
Wood, S.N., Goude, Y. & Shaw S. (2015) Generalized additive models for large datasets. Journal of the Royal Statistical Society, Series C 64(1): 139-155.
\url{http://dx.doi.org/10.1111/rssc.12068}
Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association.
Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association. 112(519):1199-1210
\url{http://dx.doi.org/10.1080/01621459.2016.1195744}
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}
......
\name{choldrop}
\alias{choldrop}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Down date a Cholesky factor on dropping a row/col}
\description{Given a Cholesky factor, \code{R}, of a matrix, \code{A}, finds the Cholesky factor of \code{A[-k,-k]},
where \code{k} is an integer.
}
\usage{
choldrop(R,k)
}
%- maybe also `usage' for other objects documented here.
\arguments{
\item{R}{Cholesky factor of a matrix, \code{A}.}
\item{k}{row and column of \code{A} to drop.}
}
\details{If \code{R} is upper triangular then \code{t(R[,-k])\%*\%R[,-k] == A[-k,-k]}, but \code{R[,-k]} has elements on the first sub-diagonal, from its kth column onwards. To get from this to a triangular Cholesky factor of \code{A[-k,-k]} we can apply a sequence of Given rotations from the left to eliminate the sub-diagonal elements. The routine does this. If \code{R} is a lower triangular factor then Givens rotations from the right are needed to remove the extra elements. If \code{n} is the dimension of \code{R} then the update has O(n^2) computational cost.
Note that the update is vector oriented, and is hence not susceptible to speed up by use of an optimized BLAS. The update is set up to be relatively Cache friendly, in that in the upper triangular case successive Givens rotations are stored for sequential application columnwise, rather than being applied rowwise as soon as they are computed. Even so, the upper triangular update is slightly slower than the lower triangular update.
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}}
\examples{
require(mgcv)
set.seed(0)
n <- 6
A <- crossprod(matrix(runif(n*n),n,n))
R0 <- chol(A)
k <- 3
Rd <- choldrop(R0,k)
range(Rd-chol(A[-k,-k]))
Rd;chol(A[-k,-k])
## same but using lower triangular factor A = LL'
L <- t(R0)
Ld <- choldrop(L,k)
range(Ld-t(chol(A[-k,-k])))
Ld;t(chol(A[-k,-k]))
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..
......@@ -200,19 +200,6 @@ mean(res);var(res)
res <- colSums(log(abs(Xp \%*\% t(br))))
##################################################################
# illustration of unsafe scale dependent transforms in smooths....
##################################################################
b0 <- gam(y~s(x0)+s(x1)+s(x2)+x3,data=dat) ## safe
b1 <- gam(y~s(x0)+s(I(x1/2))+s(x2)+scale(x3),data=dat) ## safe
b2 <- gam(y~s(x0)+s(scale(x1))+s(x2)+scale(x3),data=dat) ## unsafe
pd <- dat; pd$x1 <- pd$x1/2; pd$x3 <- pd$x3/2
par(mfrow=c(1,2))
plot(predict(b0,pd),predict(b1,pd),main="b0 and b1 predictions match")
abline(0,1,col=2)
plot(predict(b0,pd),predict(b2,pd),main="b2 unsafe, doesn't match")
abline(0,1,col=2)
##################################################################
......@@ -245,6 +232,21 @@ se <- sqrt(x0\%*\%b$Vp\%*\%t(x0));se ## get standard error
predict(b,newdata=data.frame(x0=xn[1],x1=xn[2],
x2=xn[3],x3=xn[4]),se=TRUE)
##################################################################
# illustration of unsafe scale dependent transforms in smooths....
##################################################################
b0 <- gam(y~s(x0)+s(x1)+s(x2)+x3,data=dat) ## safe
b1 <- gam(y~s(x0)+s(I(x1/2))+s(x2)+scale(x3),data=dat) ## safe
b2 <- gam(y~s(x0)+s(scale(x1))+s(x2)+scale(x3),data=dat) ## unsafe
pd <- dat; pd$x1 <- pd$x1/2; pd$x3 <- pd$x3/2
par(mfrow=c(1,2))
plot(predict(b0,pd),predict(b1,pd),main="b0 and b1 predictions match")
abline(0,1,col=2)
plot(predict(b0,pd),predict(b2,pd),main="b2 unsafe, doesn't match")
abline(0,1,col=2)
####################################################################
## Differentiating the smooths in a model (with CIs for derivatives)
####################################################################
......
......@@ -10,7 +10,8 @@
R_CallMethodDef CallMethods[] = {
{"mgcv_pmmult2", (DL_FUNC) &mgcv_pmmult2,5},
{"mgcv_Rpiqr", (DL_FUNC) &mgcv_Rpiqr,5},
{ "mgcv_tmm",(DL_FUNC)&mgcv_tmm,5},
{ "mgcv_tmm",(DL_FUNC)&mgcv_tmm,5},
{ "mgcv_chol_down",(DL_FUNC)&mgcv_chol_down,5},
{ "mgcv_Rpbsi",(DL_FUNC)&mgcv_Rpbsi,2},
{ "mgcv_RPPt",(DL_FUNC)&mgcv_RPPt,3},
{ "mgcv_Rpchol",(DL_FUNC)&mgcv_Rpchol,4},
......@@ -34,6 +35,7 @@ R_CMethodDef CEntries[] = {
{"Xbd", (DL_FUNC) &Xbd,15},
{"vcorr", (DL_FUNC) &vcorr, 5},
{"dchol", (DL_FUNC) &dchol, 4},
{"chol_down", (DL_FUNC) &chol_down, 5},
{"mgcv_omp", (DL_FUNC) &mgcv_omp, 1},
{"coxpred", (DL_FUNC) &coxpred, 14},
{"coxpp", (DL_FUNC) &coxpp, 10},
......
......@@ -52,7 +52,7 @@
#include <stdio.h>
#include <math.h>
#include <R.h>
#include <R_ext/Linpack.h> /* only needed for pivoted chol - see note in mgcv_chol */
//#include <R_ext/Linpack.h> /* only needed for pivoted chol - see note in mgcv_chol */
#include <R_ext/Lapack.h>
#include <R_ext/BLAS.h>
#ifdef OPENMP_ON
......@@ -1710,6 +1710,86 @@ void dchol(double *dA, double *R, double *dR,int *p) {
} /* dchol */
void mgcv_chol_down(SEXP r,SEXP ru,SEXP N,SEXP K, SEXP UT) {
/* wrapper for calling chol_down using .Call */
double *R,*Rup;
int *n,*k,*ut;
R = REAL(r);Rup=REAL(ru);
n = INTEGER(N);
k = INTEGER(K);
ut = INTEGER(UT);
chol_down(R,Rup,n,k,ut);
}
void chol_down(double *R,double *Rup, int *n,int *k,int *ut) {
/* R is an n by n choleski factor of an n by n matrix A. We want the downdated
factor for A[-k,-k] returned in n-1 by n-1 matrix Rup.
If ut!=0 then R'R = A, with R and Rup upper triangular.
Otherwise RR'=A, with R and Rup lower triangular. The latter update
is more Cache friendly, since the update is then from the right and operates
columnwise. Calls from R should ideally be made from a wrapper called from
.Call, since otherwise copying can be the dominant cost.
Currently called directly with .C for accuracy testing. Assumes Rup zeroed on entry.
*/
int i,j,n1;
double x,*Ri1,*Ri,*Rj,c,s,*Re,*ca,*sa,*sp,*cp;
n1 = *n-1;
if (*ut) { /* upper trianglar col oriented computation */
ca = R + 2;sa = ca + *n; /* Givens storage */
for (i=0;i<n1;i++) { /* work across columns */
Ri = Rup + i * n1; /* destination column */
if (i < *k) {
Rj = R + i * *n;
Re = Rj + i;
} else {
Rj = R + (i+1) * *n; /* source column */
Re = Rj + *k;
}
for (;Rj<=Re;Ri++,Rj++) *Ri=*Rj; /* copy column starts */
/* now the updates */
if (i >= *k) { /* cols from beyond dropped col k need updating */
/* first the stored rotations */
Re = Rup + i * n1 + i;
Ri1=Ri;Ri--;
for (cp=ca,sp=sa;Ri<Re;cp++,sp++,Ri++,Ri1++,Rj++) {
*Ri1 = - *sp * *Ri + *cp * *Rj;
*Ri = *cp * *Ri + *sp * *Rj;
}
/* Now the final rotation */
x = sqrt(*Ri * *Ri + *Rj * *Rj);
c = *Ri / x; s = *Rj / x;*Ri = x;
if (i<n1-1) { *cp = c;*sp=s;} /* store rotation */
}
}
/* now clear Givens storage */
Re = R + *n;
for (;ca<Re;ca++,sa++) *ca = *sa = 0.0;
} else { /* lower triangular */
/* First copy the input rows before row k */
for (i=0;i< *k;i++) {
for (Ri=Rup+i*n1,Rj=R+i * *n,Re=Ri+*k;Ri<Re;Ri++,Rj++) *Ri = *Rj;
}
/* now copy cols 0:k from row k+1 */
for (i=0;i<=*k;i++) {
for (Ri=Rup+i*n1,Re=Ri+n1,Ri+= *k,Rj=R+i * *n+ *k+1;Ri<Re;Ri++,Rj++) *Ri = *Rj;
}
for (i = *k;i<n1;i++) { /* work down diagonal */
Ri = Rup + i + i * n1; /* rotate into here */
Rj = R + (i+1) + (i+1) * *n; /* zeroing this */
x = sqrt(*Ri * *Ri + *Rj * *Rj);
c = *Ri / x; s = *Rj / x;*Ri = x;
Rj++;Ri++;Ri1=Ri+n1;Re = Rup + (i+1) * n1;
for (;Ri<Re;Rj++,Ri++,Ri1++) {
*Ri1 = -s * *Ri + c * *Rj;
*Ri = c * *Ri + s * *Rj;
}
}
}
} /* chol_down */
void mgcv_chol(double *a,int *pivot,int *n,int *rank)
/* a stored in column order, this routine finds the pivoted choleski decomposition of matrix a
library(mgcv)
......@@ -1723,14 +1803,15 @@ void mgcv_chol(double *a,int *pivot,int *n,int *rank)
rD<-rD[,ind]
L<-mroot(D)
D;t(rD)%*%rD;L%*%t(L)
NOTE: This uses LINPACK - dpstf2.f is LAPACK version, but not in R headers yet!
*/
{ double *work,*p1,*p2,*p;
int piv=1;
work=(double *)CALLOC((size_t) *n,sizeof(double));
F77_CALL(dchdc)(a,n,n,work,pivot,&piv,rank);
{ double *work,*p1,*p2,*p,*p3,tol=-1.0;
int info=1;
char uplo='U';
work=(double *)CALLOC((size_t) *n * 2,sizeof(double));
//F77_CALL(dchdc)(a,n,n,work,pivot,&info,rank); /* LINPACK */
F77_CALL(dpstrf)(&uplo,n,a,n,pivot,rank,&tol,work,&info); /* LAPACK */
/* zero stuff below the leading diagonal */
for (p2=a+ *n,p1=a+1;p2<a+ *n * *n;p1+= *n+1,p2+= *n) for (p=p1;p<p2;p++) *p=0.0;
for (p2=a + *n,p1=a+1,p3=a+ *n * *n;p2<p3;p1+= *n+1,p2+= *n) for (p=p1;p<p2;p++) *p=0.0;
FREE(work);
}
......
......@@ -161,6 +161,8 @@ void mgcv_Rpbsi(SEXP A, SEXP NT);
void mgcv_RPPt(SEXP a,SEXP r, SEXP NT);
SEXP mgcv_Rpchol(SEXP Amat,SEXP PIV,SEXP NT,SEXP NB);
void dchol(double *dA, double *R, double *dR,int *p);
void chol_down(double *R,double *Rup,int *n,int *k,int *ut);
void mgcv_chol_down(SEXP r,SEXP ru,SEXP N,SEXP K, SEXP UT);
void vcorr(double *dR,double *Vr,double *Vb,int *p,int *M);
SEXP mgcv_Rpforwardsolve(SEXP R, SEXP B,SEXP NT);
SEXP mgcv_Rpbacksolve(SEXP R, SEXP B,SEXP NT);
......
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