diff -Nru mgcv-1.8-27/ChangeLog mgcv-1.8-28/ChangeLog --- mgcv-1.8-27/ChangeLog 2019-02-02 18:46:31.000000000 +0000 +++ mgcv-1.8-28/ChangeLog 2019-03-21 08:29:07.000000000 +0000 @@ -17,6 +17,18 @@ * OpenBLAS 0.3.2/3 appears not to be thread safe at present - can cause problems when opemmp multithreading. +1.8-28 + +* fix of obscure sp naming bug. + +* changed some contour default colours from green to blue (they overlay + heatmaps, so green was not clever). + +* Tweedie likelihood evaluation code made slightly more robust - for a model + with machine zero scale parameter estimate it could segfault, as series + maximum location could then overflow integer storage. Fixed + upper limit + imposed on series length (warning if it's not enough). + 1.8-27 ** Added routine 'ginla' for fully Bayesian inference based on an integrated @@ -25,7 +37,7 @@ * Tweedie location scale family added: 'twlss'. * gam.fit5 modified to distinguish more carefully between +ve semi definite - and +ve definite. Previously could fail, claiming indefinteness when it + and +ve definite. Previously could fail, claiming indefiniteness when it should not have. Affects general families. * bam was ignoring supplied scale parameter in extended family cases - fixed. @@ -71,7 +83,7 @@ sarse matrices). Useful for other packages using constructors and smooth2random, for 'fs' smooths. -* The mrf smooth constructor contained a obsolete hack in which the term +* The mrf smooth constructor contained an 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. diff -Nru mgcv-1.8-27/debian/changelog mgcv-1.8-28/debian/changelog --- mgcv-1.8-27/debian/changelog 2019-03-27 03:12:19.000000000 +0000 +++ mgcv-1.8-28/debian/changelog 2019-03-27 03:12:19.000000000 +0000 @@ -1,14 +1,14 @@ -mgcv (1.8-27-1cran1xenial0) xenial; urgency=medium +mgcv (1.8-28-1cran1xenial0) xenial; urgency=medium - * Compilation for Ubuntu 16.04.5 LTS + * Compilation for Ubuntu 16.04.6 LTS - -- Michael Rutter Sun, 10 Feb 2019 20:50:24 +0000 + -- Michael Rutter Wed, 27 Mar 2019 03:10:31 +0000 -mgcv (1.8-27-1cran1) testing; urgency=low +mgcv (1.8-28-1cran1) testing; urgency=low * cran2deb svn: 362M with DB version 1. - -- cran2deb4ubuntu Sat, 09 Feb 2019 09:11:07 -0500 + -- cran2deb4ubuntu Sat, 23 Mar 2019 07:47:41 -0400 mgcv (1.8-9-1cran1) testing; urgency=low @@ -60,6 +60,13 @@ -- cran2deb4ubuntu Sat, 30 Aug 2014 08:40:20 -0400 +mgcv (1.8-27-1cran1) testing; urgency=low + + * cran2deb svn: 362M with DB version 1. + + -- cran2deb4ubuntu Sat, 09 Feb 2019 09:11:15 -0500 + + mgcv (1.8-26-1cran1) testing; urgency=low * cran2deb svn: 362M with DB version 1. diff -Nru mgcv-1.8-27/DESCRIPTION mgcv-1.8-28/DESCRIPTION --- mgcv-1.8-27/DESCRIPTION 2019-02-06 15:00:03.000000000 +0000 +++ mgcv-1.8-28/DESCRIPTION 2019-03-21 11:40:07.000000000 +0000 @@ -1,5 +1,5 @@ Package: mgcv -Version: 1.8-27 +Version: 1.8-28 Author: Simon Wood Maintainer: Simon Wood Title: Mixed GAM Computation Vehicle with Automatic Smoothness @@ -20,6 +20,6 @@ ByteCompile: yes License: GPL (>= 2) NeedsCompilation: yes -Packaged: 2019-02-06 10:57:24 UTC; sw283 +Packaged: 2019-03-21 08:56:32 UTC; sw283 Repository: CRAN -Date/Publication: 2019-02-06 15:00:03 UTC +Date/Publication: 2019-03-21 11:40:07 UTC diff -Nru mgcv-1.8-27/man/bam.Rd mgcv-1.8-28/man/bam.Rd --- mgcv-1.8-27/man/bam.Rd 2018-12-11 20:35:40.000000000 +0000 +++ mgcv-1.8-28/man/bam.Rd 2019-03-18 12:42:15.000000000 +0000 @@ -13,7 +13,7 @@ of \code{bam} is much lower memory footprint than \code{\link{gam}}, but it can also be much faster, for large datasets. \code{bam} can also compute on a cluster set up by the \link[parallel]{parallel} package. -An alternative fitting approach (Wood et al. 2017) is provided by the \code{discrete==TRUE} method. In this case a method based on discretization of covariate values and C code level parallelization (controlled by the \code{nthreads} argument instead of the \code{cluster} argument) is used. This extends both the data set and model size that are practical. +An alternative fitting approach (Wood et al. 2017, Li and Wood, 2019) is provided by the \code{discrete==TRUE} method. In this case a method based on discretization of covariate values and C code level parallelization (controlled by the \code{nthreads} argument instead of the \code{cluster} argument) is used. This extends both the data set and model size that are practical. } \usage{ @@ -187,7 +187,7 @@ If \code{cluster} is provided as a cluster set up using \code{\link[parallel]{makeCluster}} (or \code{\link[parallel]{makeForkCluster}}) from the \code{parallel} package, then the rate limiting QR decomposition of the model matrix is performed in parallel using this cluster. Note that the speed ups are often not that great. On a multi-core machine it is usually best to set the cluster size to the number of physical cores, which is often less than what is reported by \code{\link[parallel]{detectCores}}. Using more than the number of physical cores can result in no speed up at all (or even a slow down). Note that a highly parallel BLAS may negate all advantage from using a cluster of cores. Computing in parallel of course requires more memory than computing in series. See examples. -When \code{discrete=TRUE} the covariate data are first discretized. Discretization takes place on a smooth by smooth basis, or in the case of tensor product smooths (or any smooth that can be represented as such, such as random effects), separately for each marginal smooth. The required spline bases are then evaluated at the discrete values, and stored, along with index vectors indicating which original observation they relate to. Fitting is by a version of performance oriented iteration/PQL using REML smoothing parameter selection on each iterative working model (as for the default method). The iteration is based on the derivatives of the REML score, without computing the score itself, allowing the expensive computations to be reduced to one parallel block Cholesky decomposition per iteration (plus two basic operations of equal cost, but easily parallelized). Unlike standard POI/PQL, only one step of the smoothing parameter update for the working model is taken at each step (rather than iterating to the optimal set of smoothing parameters for each working model). At each step a weighted model matrix crossproduct of the model matrix is required - this is efficiently computed from the pre-computed basis functions evaluated at the discretized covariate values. Efficient computation with tensor product terms means that some terms within a tensor product may be re-ordered for maximum efficiency. +When \code{discrete=TRUE} the covariate data are first discretized. Discretization takes place on a smooth by smooth basis, or in the case of tensor product smooths (or any smooth that can be represented as such, such as random effects), separately for each marginal smooth. The required spline bases are then evaluated at the discrete values, and stored, along with index vectors indicating which original observation they relate to. Fitting is by a version of performance oriented iteration/PQL using REML smoothing parameter selection on each iterative working model (as for the default method). The iteration is based on the derivatives of the REML score, without computing the score itself, allowing the expensive computations to be reduced to one parallel block Cholesky decomposition per iteration (plus two basic operations of equal cost, but easily parallelized). Unlike standard POI/PQL, only one step of the smoothing parameter update for the working model is taken at each step (rather than iterating to the optimal set of smoothing parameters for each working model). At each step a weighted model matrix crossproduct of the model matrix is required - this is efficiently computed from the pre-computed basis functions evaluated at the discretized covariate values. Efficient computation with tensor product terms means that some terms within a tensor product may be re-ordered for maximum efficiency. See Wood et al (2017) and Li and Wood (2019) for full details. When \code{discrete=TRUE} parallel computation is controlled using the \code{nthreads} argument. For this method no cluster computation is used, and the \code{parallel} package is not required. Note that actual speed up from parallelization depends on the BLAS installed and your hardware. With the (R default) reference BLAS using several threads can make a substantial difference, but with a single threaded tuned BLAS, such as openblas, the effect is less marked (since cache use is typically optimized for one thread, and is then sub optimal for several). However the tuned BLAS is usually much faster than using the reference BLAS, however many threads you use. If you have a multi-threaded BLAS installed then you should leave \code{nthreads} at 1, since calling a multi-threaded BLAS from multiple threads usually slows things down: the only exception to this is that you might choose to form discrete matrix cross products (the main cost in the fitting routine) in a multi-threaded way, but use single threaded code for other computations: this can be achieved by e.g. \code{nthreads=c(2,1)}, which would use 2 threads for discrete inner products, and 1 for most code calling BLAS. Not that the basic reason that multi-threaded performance is often disappointing is that most computers are heavily memory bandwidth limited, not flop rate limited. It is hard to get data to one core fast enough, let alone trying to get data simultaneously to several cores. @@ -204,6 +204,10 @@ 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} + +Li, Z & S.N. Wood (2019) Faster model matrix crossproducts for large generalized linear models with discretized covariates. Statistics and Computing. +\url{https://doi.org/10.1007/s11222-019-09864-2} + } \author{ Simon N. Wood \email{simon.wood@r-project.org} diff -Nru mgcv-1.8-27/man/gamSim.Rd mgcv-1.8-28/man/gamSim.Rd --- mgcv-1.8-27/man/gamSim.Rd 2017-04-11 14:08:24.000000000 +0000 +++ mgcv-1.8-28/man/gamSim.Rd 2019-03-18 12:42:15.000000000 +0000 @@ -16,7 +16,7 @@ \item{n}{ number of data to simulate.} -\item{dist}{character string which may be used to spcify the distribution of +\item{dist}{character string which may be used to specify the distribution of the response.} \item{scale}{Used to set noise level.} diff -Nru mgcv-1.8-27/man/ginla.Rd mgcv-1.8-28/man/ginla.Rd --- mgcv-1.8-27/man/ginla.Rd 2018-12-15 18:53:23.000000000 +0000 +++ mgcv-1.8-28/man/ginla.Rd 2019-03-18 12:42:15.000000000 +0000 @@ -2,7 +2,7 @@ \alias{ginla} %- Also NEED an `\alias' for EACH other topic documented here. \title{GAM Integrated Nested Laplace Approximation Newton Enhanced} -\description{Apply Integrated Nested Laplace Approximation (INLA, Rue et al. 2009) to models estimable by \code{\link{gam}} or \code{\link{bam}}, using the INLA variant described in Wood (2018). Produces marginal posterior densities for each coefficient, selected coefficients or linear transformations of the coefficient vector. +\description{Apply Integrated Nested Laplace Approximation (INLA, Rue et al. 2009) to models estimable by \code{\link{gam}} or \code{\link{bam}}, using the INLA variant described in Wood (2019). Produces marginal posterior densities for each coefficient, selected coefficients or linear transformations of the coefficient vector. } \usage{ ginla(G,A=NULL,nk=16,nb=100,J=1,interactive=FALSE,int=0,approx=0) @@ -36,7 +36,7 @@ \references{ Rue, H, Martino, S. & Chopin, N. (2009) Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society, Series B. 71: 319-392. -Wood (2018, submitted) Simplified Integrated Laplace Approximation. +Wood (2019) Simplified Integrated Laplace Approximation. In press Biometrika. } \section{WARNINGS}{ diff -Nru mgcv-1.8-27/man/ldetS.Rd mgcv-1.8-28/man/ldetS.Rd --- mgcv-1.8-27/man/ldetS.Rd 2017-07-29 12:21:00.000000000 +0000 +++ mgcv-1.8-28/man/ldetS.Rd 2019-03-21 08:56:25.000000000 +0000 @@ -4,7 +4,7 @@ \alias{ldetS} \title{Getting log generalized determinant of penalty matrices} \usage{ -ldetS(Sl, rho, fixed, np, root = FALSE, repara = TRUE, nt = 1) +ldetS(Sl, rho, fixed, np, root = FALSE, repara = TRUE, nt = 1,deriv=2) } \arguments{ \item{Sl}{the output of \code{Sl.setup}.} @@ -21,6 +21,8 @@ a re-parameterization object supplied in the returned object.} \item{nt}{number of parallel threads to use.} + +\item{deriv}{maximum order of derivatives to compute.} } \value{ A list containing: \itemize{ diff -Nru mgcv-1.8-27/man/smooth.construct.bs.smooth.spec.Rd mgcv-1.8-28/man/smooth.construct.bs.smooth.spec.Rd --- mgcv-1.8-27/man/smooth.construct.bs.smooth.spec.Rd 2018-05-30 16:18:39.000000000 +0000 +++ mgcv-1.8-28/man/smooth.construct.bs.smooth.spec.Rd 2019-03-18 12:42:15.000000000 +0000 @@ -84,7 +84,7 @@ fv <- predict(b,pd,se=TRUE) ul <- fv$fit + fv$se.fit*2; ll <- fv$fit - fv$se.fit*2 plot(x,y,xlim=c(-.7,1.7),ylim=range(c(y,ll,ul)),main= - "1st order penalties: red tps; blue bs over (0,1); black bs over (-.5,1.5)") + "Order 1 penalties: red tps; blue bs on (0,1); black bs on (-.5,1.5)") ## penalty defined on (-.5,1.5) gives plausible predictions and intervals ## over this range... lines(pd$x,fv$fit);lines(pd$x,ul,lty=2);lines(pd$x,ll,lty=2) diff -Nru mgcv-1.8-27/man/smooth.construct.Rd mgcv-1.8-28/man/smooth.construct.Rd --- mgcv-1.8-27/man/smooth.construct.Rd 2017-04-11 14:08:22.000000000 +0000 +++ mgcv-1.8-28/man/smooth.construct.Rd 2019-03-18 12:42:15.000000000 +0000 @@ -211,10 +211,10 @@ ## using this, since mgcv can happily handle non-identity ## penalties.) -smooth.construct.tr.smooth.spec<-function(object,data,knots) +smooth.construct.tr.smooth.spec<-function(object,data,knots) { ## a truncated power spline constructor method function ## object$p.order = null space dimension -{ m <- object$p.order[1] + m <- object$p.order[1] if (is.na(m)) m <- 2 ## default if (m<1) stop("silly m supplied") if (object$bs.dim<0) object$bs.dim <- 10 ## default @@ -249,9 +249,9 @@ object } -Predict.matrix.tr.smooth<-function(object,data) +Predict.matrix.tr.smooth<-function(object,data) { ## prediction method function for the `tr' smooth class -{ x <- data[[object$term]] + x <- data[[object$term]] x <- x - object$x.shift # stabilizing shift m <- object$m; # spline order (3=cubic) k<-object$knots # knot locations diff -Nru mgcv-1.8-27/man/uniquecombs.Rd mgcv-1.8-28/man/uniquecombs.Rd --- mgcv-1.8-27/man/uniquecombs.Rd 2017-04-11 14:08:22.000000000 +0000 +++ mgcv-1.8-28/man/uniquecombs.Rd 2019-03-18 12:42:15.000000000 +0000 @@ -54,7 +54,7 @@ either have no \code{as.character} method, or whose \code{as.character} method is a many to one mapping, then the routine is likely to fail. -If the character representation of a dataframe variable (other than of class factor of character) contains \code{*} then in priniciple the method could fail (but with a warning). +If the character representation of a dataframe variable (other than of class factor of character) contains \code{*} then in principle the method could fail (but with a warning). } \examples{ diff -Nru mgcv-1.8-27/MD5 mgcv-1.8-28/MD5 --- mgcv-1.8-27/MD5 2019-02-06 15:00:03.000000000 +0000 +++ mgcv-1.8-28/MD5 2019-03-21 11:40:07.000000000 +0000 @@ -1,22 +1,22 @@ -fd791b31993ade7c507f85afff3242e2 *ChangeLog -136bd8b6da7fe438d10fa69479e60c31 *DESCRIPTION +5fa8085364f960a26ac7690edfd6e423 *ChangeLog +1ca0566710cd19035e444ca30d94eda8 *DESCRIPTION eb723b61539feef013de476e68b5c50a *GPL-2 a7edea9aaf9e16a95b7844f9161515d6 *NAMESPACE -d35108707f4ccc3b5237626aa6dbe546 *R/bam.r +fac2c9f6811b7e760da5d55bac9b6270 *R/bam.r 2c132dc64eedf0c17d34f7061346f2d0 *R/coxph.r -5e3e5d8699554ed9a3593195055d99ba *R/efam.r -efcb37e8ceeba5da36cb1e3ac36a2bbe *R/fast-REML.r -4eaaa4960d578cd66a4201e6d0dc1d3a *R/gam.fit3.r -9ba6e588b9bbd50a3165fdfb8c16d23a *R/gam.fit4.r +36b00131cb4d9d98b45d8d4acb4e5953 *R/efam.r +b6953debd773ebbaab28afe97047e34c *R/fast-REML.r +81a226e30a5e69e0620cfbb6029a283c *R/gam.fit3.r +922b0786419a82a2127a94b062befc4d *R/gam.fit4.r 1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r -021b3c1566fa81b6260c17b199b011c8 *R/gamlss.r +8e457764d72567f936b70d5811473f9e *R/gamlss.r ce69039c15d4db9cfda5d1309af8b242 *R/gamm.r f4c20e36b4f1d42c5bf902270ba3515a *R/inla.r 10facb791e4cfd123d183f05660119c6 *R/jagam.r -1a82805404b2027b46a124f020da30e2 *R/mgcv.r +05722ab28faea920a01ca48b07ff30ba *R/mgcv.r 97e439c8d1f8a0562d1383d4e4877fe2 *R/misc.r b73809727c6eb8f7c7b0a20e9625fddc *R/mvam.r -eed19ceca4c4114023f12cda3f622659 *R/plots.r +1d3c2d5dae39bdb462ed75411829957d *R/plots.r ae5ad4af1b30efb65b5deb05f6b8e778 *R/smooth.r 7398607a9ba7e85276bcf09df0d9344b *R/soap.r bde1774ce7903cabc57b3196f8872ea8 *R/sparse.r @@ -44,7 +44,7 @@ 60670020425f8749b81a8d8c3f168880 *man/Sl.setup.Rd 69ae63833438a3af2963e39482f1d72f *man/Tweedie.Rd 8087ab00d10b44c99c37f49bf90e19cd *man/anova.gam.Rd -4114840d3ed1b238b0430eed07fbd44c *man/bam.Rd +9e06549c2c3ca2e4945dfb31693512a4 *man/bam.Rd ab5e37c3bf8803de63b63c3bdc5909cd *man/bam.update.Rd cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd 745cbf31eb14fc1c5916fc634c74d998 *man/bug.reports.mgcv.Rd @@ -81,14 +81,14 @@ b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd eb8648cc6b3c9374b899abd2f2b12f7b *man/gam2objective.Rd ed77ce6e1b941625485706d7e307b816 *man/gamObject.Rd -a2593990d6a87f7b783d0435062dec02 *man/gamSim.Rd +89148f2dc12caff5073ac70c8873b33f *man/gamSim.Rd e5d2541f32dab56972f58b0773eba50c *man/gamlss.etamu.Rd c7f140d128d1d1d76909499900faf49e *man/gamlss.gH.Rd 623851672072e852e0c0ff03a96de0d6 *man/gamm.Rd 222535dd19201cfd929efd3305b13f43 *man/gaulss.Rd 398a5c12285401c1d37a8edb58780bc3 *man/get.var.Rd a62dd487f34f476f7f830ed9d1bc58dc *man/gevlss.Rd -3158643d44410c71f15075fab0e516e0 *man/ginla.Rd +09dba82ee7459800056dbd0511788ebf *man/ginla.Rd 39b47f30a7ea45382d62ca1753d876a8 *man/identifiability.Rd 4f96476abbf9692f52030d3859580a31 *man/in.out.Rd 6c33026ebb458483d34a04548c05d664 *man/inSide.Rd @@ -98,7 +98,7 @@ 4bc4d96708fc6454ad81207299608d28 *man/jagam.Rd d37e837db3089e3c0eb105669c04f0c8 *man/k.check.Rd 87d942b17bee05bb662270b894b183e6 *man/ldTweedie.Rd -c7073b72336b0385080f2f4a40694416 *man/ldetS.Rd +adc52f0e82ccbbe115bf15045b0a15a4 *man/ldetS.Rd 1b314907562e9236747ec1e344ebe491 *man/linear.functional.terms.Rd ab35592c0a29d1574579811ea1c6ec39 *man/logLik.gam.Rd b1c95a20afd6eb0825c00b46b8c3cbfa *man/ls.size.Rd @@ -142,9 +142,9 @@ 898e7cc2def2ee234475e68d0b904b29 *man/sdiag.Rd d54f4042e212fca7704cf8428bdaea38 *man/single.index.Rd 6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd -8020154bd5c709d11f0e7cf043df886d *man/smooth.construct.Rd +40a5e35173c808dde102e2bfc492ac9d *man/smooth.construct.Rd 3f3e0cd76b77e207ee5a6ff89e5f7a9f *man/smooth.construct.ad.smooth.spec.Rd -fc7fba34b89fdf29f66744d1fdadda79 *man/smooth.construct.bs.smooth.spec.Rd +d2f3fb49a85066ef08d395a244d08a4d *man/smooth.construct.bs.smooth.spec.Rd 2f0463d1aca0b8798da6e681bd4c6e54 *man/smooth.construct.cr.smooth.spec.Rd f5e6d0f5122f61c336827b3615482157 *man/smooth.construct.ds.smooth.spec.Rd bd58515156b5e07c006e69f29aa830d1 *man/smooth.construct.fs.smooth.spec.Rd @@ -171,7 +171,7 @@ f22f1cee0ff2b70628846d1d0f8e9a66 *man/trichol.Rd 87e6b4437d00fab4fc814f4cefa3795c *man/trind.generator.Rd 6e975eef6b1214e0be93fc641f982f67 *man/twlss.Rd -96c48dd705710f639d76a0d0cc3fb128 *man/uniquecombs.Rd +bc350bfd3f4f8316d3b29b247292a16d *man/uniquecombs.Rd a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd 281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd cbb69f16706da27b44fc62c92dcae161 *man/ziP.Rd @@ -195,12 +195,12 @@ 2436f9b328e80370ce2203dbf1dd813c *src/general.h 87f4add06d78bdaebf588eadbdbffbbf *src/init.c a9151b5236852eef9e6947590bfcf88a *src/magic.c -613c7bccb3d307824cf278b829190087 *src/mat.c +cd7abd2d8c808315207aeb566065efbb *src/mat.c e4cef7f1753153fbab242d1c4d4f7e7f *src/matrix.c de37b0972199b796654405efc007f25b *src/matrix.h 8df4b96961491d76989b50856237ee2d *src/mgcv.c c7b7d5e90758bb042f206c2e9e93b129 *src/mgcv.h -8e3efe7a0d6b9619671cf2f7471f2112 *src/misc.c +0e9733504235e94635c6688ce6bb6bb1 *src/misc.c 465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c 563938b7bb6504ab10df5376c4360220 *src/qp.c 073a4b5b0bc6e869c5b35478c47facf1 *src/qp.h diff -Nru mgcv-1.8-27/R/bam.r mgcv-1.8-28/R/bam.r --- mgcv-1.8-27/R/bam.r 2018-12-12 19:58:07.000000000 +0000 +++ mgcv-1.8-28/R/bam.r 2019-03-18 12:42:14.000000000 +0000 @@ -1,5 +1,5 @@ ## routines for very large dataset generalized additive modelling. -## (c) Simon N. Wood 2009-2017 +## (c) Simon N. Wood 2009-2019 ls.size <- function(x) { @@ -249,7 +249,6 @@ nmarg <- if (inherits(gp$smooth.spec[[i]],"tensor.smooth.spec")) length(gp$smooth.spec[[i]]$margin) else 1 maxj <- if (gp$smooth.spec[[i]]$by=="NA") nmarg else nmarg + 1 mi <- if (is.null(m)||length(m)==1) m else m[i] -# if (!is.null(gp$smooth.spec[[i]]$id)) stop("discretization can not handle smooth ids") j <- 0 for (jj in 1:maxj) { ## loop through marginals if (jj==1&&maxj!=nmarg) termi <- gp$smooth.spec[[i]]$by else { @@ -599,24 +598,15 @@ w <- dd$EDeta2 * .5 z <- (eta-offset) - dd$Deta.EDeta2 } - #if (rho!=0) { - # ind <- which(w<0) - # if (length(ind)>0) { ## substitute Fisher weights - # w[ind] <- dd$EDeta2[ind] * .5 - # z[ind] <- (eta[ind]-offset[ind]) - dd$Deta.EDeta2[ind] - # } - #} good <- is.finite(z)&is.finite(w) w[!good] <- 0 ## drop if !good z[!good] <- 0 ## irrelevant - #dev <- sum(family$dev.resids(G$y,mu,G$w,theta)) } else { ## exponential family mu.eta.val <- mu.eta(eta) good <- mu.eta.val != 0 mu.eta.val[!good] <- .1 ## irrelvant as weight is zero z <- (eta - offset) + (G$y - mu)/mu.eta.val w <- (G$w * mu.eta.val^2)/variance(mu) - #dev <- sum(family$dev.resids(G$y,mu,G$w)) } @@ -655,7 +645,7 @@ ## block diagonal penalty object, Sl, set up before loop if (iter==1) { ## need to get initial smoothing parameters - lambda.0 <- initial.sp(qrx$R,G$S,G$off,XX=TRUE) ## note that this uses the untrasformed X'X in qrx$R + lambda.0 <- initial.sp(qrx$R,G$S,G$off,XX=TRUE) ## note that this uses the untransformed X'X in qrx$R ## convert intial s.p.s to account for L lsp0 <- log(lambda.0) ## initial s.p. if (!is.null(G$L)) lsp0 <- diff -Nru mgcv-1.8-27/R/efam.r mgcv-1.8-28/R/efam.r --- mgcv-1.8-27/R/efam.r 2019-01-24 09:40:07.000000000 +0000 +++ mgcv-1.8-28/R/efam.r 2019-03-18 12:42:14.000000000 +0000 @@ -935,7 +935,7 @@ #if (vec) lsth1 <- LS[,c(4,2)] LS <- colSums(LS) #if (!vec) lsth1 <- c(LS[4],LS[2]) - lsth1 <- c(LS[4],LS[2]) + lsth1 <- c(LS[4],LS[2]) ## deriv w.r.t. p then log scale lsth2 <- matrix(c(LS[5],LS[6],LS[6],LS[3]),2,2) list(ls=LS[1],lsth1=lsth1,lsth2=lsth2) } diff -Nru mgcv-1.8-27/R/fast-REML.r mgcv-1.8-28/R/fast-REML.r --- mgcv-1.8-27/R/fast-REML.r 2019-02-01 20:00:20.000000000 +0000 +++ mgcv-1.8-28/R/fast-REML.r 2019-03-20 13:18:38.000000000 +0000 @@ -401,7 +401,7 @@ list(det = 2*sum(log(diag(R))+log(d[piv])),det1=dS1,det2=dS2,E=E) } ## ldetSblock -ldetS <- function(Sl,rho,fixed,np,root=FALSE,repara=TRUE,nt=1) { +ldetS <- function(Sl,rho,fixed,np,root=FALSE,repara=TRUE,nt=1,deriv=2) { ## Get log generalized determinant of S stored blockwise in an Sl list. ## If repara=TRUE multi-term blocks will be re-parameterized using gam.reparam, and ## a re-parameterization object supplied in the returned object. @@ -445,13 +445,13 @@ ind <- k.sp:(k.sp+m-1) ## index for smoothing parameters ## call gam.reparam to deal with this block ## in a stable way... - grp <- if (repara) gam.reparam(Sl[[b]]$rS,lsp=rho[ind],deriv=2) else - ldetSblock(Sl[[b]]$rS,rho[ind],deriv=2,root=root,nt=nt) + grp <- if (repara) gam.reparam(Sl[[b]]$rS,lsp=rho[ind],deriv=deriv) else + ldetSblock(Sl[[b]]$rS,rho[ind],deriv=deriv,root=root,nt=nt) Sl[[b]]$lambda <- exp(rho[ind]) ldS <- ldS + grp$det ## next deal with the derivatives... grp$det1 <- grp$det1[!fixed[ind]] ## discard derivatives for fixed components - grp$det2 <- grp$det2[!fixed[ind],!fixed[ind]] + grp$det2 <- if (deriv>1) grp$det2[!fixed[ind],!fixed[ind]] else 0 ##NULL nd <- length(grp$det1) if (nd>0) { ## then not all sp's are fixed dind <- k.deriv:(k.deriv+nd-1) @@ -735,7 +735,7 @@ SA } ## end Sl.termMult -d.detXXS <- function(Sl,PP,nt=1) { +d.detXXS <- function(Sl,PP,nt=1,deriv=2) { ## function to obtain derivatives of log |X'X+S| given unpivoted PP' where ## P is inverse of R from the QR of the augmented model matrix. SPP <- Sl.termMult(Sl,PP,full=FALSE,nt=nt) ## SPP[[k]] is S_k PP' @@ -744,11 +744,13 @@ for (i in 1:nd) { indi <- attr(SPP[[i]],"ind") d1[i] <- sum(diag(SPP[[i]][,indi,drop=FALSE])) - for (j in i:nd) { - indj <- attr(SPP[[j]],"ind") - d2[i,j] <- d2[j,i] <- -sum(t(SPP[[i]][,indj,drop=FALSE])*SPP[[j]][,indi,drop=FALSE]) - } - d2[i,i] <- d2[i,i] + d1[i] + if (deriv==2) { + for (j in i:nd) { + indj <- attr(SPP[[j]],"ind") + d2[i,j] <- d2[j,i] <- -sum(t(SPP[[i]][,indj,drop=FALSE])*SPP[[j]][,indi,drop=FALSE]) + } + d2[i,i] <- d2[i,i] + d1[i] + } } list(d1=d1,d2=d2) } ## end d.detXXS @@ -803,10 +805,12 @@ ## alternative all in one code - matches loop results, but ## timing close to identical - modified for parallel exec - D <- matrix(unlist(Skb),nrow(XX),nd) + D <- matrix(unlist(Skb),length(beta),nd) bSb1 <- colSums(beta*D) D1 <- .Call(C_mgcv_Rpforwardsolve,R,D[piv,]/d[piv],nt) ## note R transposed internally unlike forwardsolve db[piv,] <- -.Call(C_mgcv_Rpbacksolve,R,D1,nt)/d[piv] + + if (is.null(XX)) return(list(bSb1=bSb1,db=db)) ## return early ## XX.db <- XX%*%db XX.db <- .Call(C_mgcv_pmmult2,XX,db,0,0,nt) @@ -825,7 +829,7 @@ nobs=0,Mp=0,nt=c(1,1),tol=0,gamma=1) { ## given X'WX in XX and f=X'Wy solves the penalized least squares problem ## with penalty defined by Sl and rho, and evaluates a REML Newton step, the REML -## gradiant and the the estimated coefs bhat. If phi.fixed=FALSE then we need +## gradient and the the estimated coefs bhat. If phi.fixed=FALSE then we need ## yy = y'Wy in order to get derivsatives w.r.t. phi. ## NOTE: with an optimized BLAS nt==1 is likely to be much faster than ## nt > 1 @@ -1198,7 +1202,7 @@ } rank <- 0 for (b in 1:length(Sl)) rank <- rank + Sl[[b]]$rank ## the total penalty rank - ## Also add Mp, the total mull space dimension to return list. + ## Also add Mp, the total null space dimension to return list. list(X=X,Sl=Sl,undrop=id$undrop,rank=rank,Mp=ncol(X)-rank) } ## end Sl.Xprep diff -Nru mgcv-1.8-27/R/gam.fit3.r mgcv-1.8-28/R/gam.fit3.r --- mgcv-1.8-27/R/gam.fit3.r 2019-02-04 16:21:47.000000000 +0000 +++ mgcv-1.8-28/R/gam.fit3.r 2019-03-18 13:29:26.000000000 +0000 @@ -28,7 +28,7 @@ rSncol <- unlist(lapply(rS,ncol)) M <- length(lsp) if (length(rS)>M) fixed.penalty <- TRUE else fixed.penalty <- FALSE - + d.tol <- .Machine$double.eps^.3 ## group `similar sized' penalties, to save work r.tol <- .Machine$double.eps^.75 ## This is a bit delicate -- too large and penalty range space can be supressed. @@ -1672,8 +1672,6 @@ } ## newton - - bfgs <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights, control,gamma,scale,conv.tol=1e-6,maxNstep=3,maxSstep=2, maxHalf=30,printWarn=FALSE,scoreType="GCV",start=NULL, @@ -1810,6 +1808,8 @@ initial$dVkk <- diag(t(L0) %*% b$dVkk %*% L0) initial$score <- score;initial$grad <- grad; initial$scale.est <- b$scale.est + start0 <- coef(b) + mustart0 <- fitted(b) rm(b) B <- diag(length(initial$grad)) ## initial Hessian @@ -1819,7 +1819,7 @@ b <- gam.fit3(x=X, y=y, sp=L%*%ilsp+lsp0,Eb=Eb,UrS=UrS, offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=1, control=control,gamma=gamma,scale=scale,printWarn=FALSE, - start=start,mustart=mustart, + start=start0,mustart=mustart0, scoreType=scoreType,null.coef=null.coef, pearson.extra=pearson.extra,dev.extra=dev.extra,n.true=n.true,Sl=Sl,...) if (reml) { @@ -2497,7 +2497,7 @@ } ## fix.family.var -fix.family.ls<-function(fam) +fix.family.ls <- function(fam) # adds ls the log saturated likelihood and its derivatives # w.r.t. the scale parameter to the family object. { if (!inherits(fam,"family")) stop("fam not a family object") @@ -2951,6 +2951,12 @@ w2pp=as.double(y*0),y=as.double(y),eps=as.double(.Machine$double.eps^2),n=as.integer(length(y)), th=as.double(theta[ind]),rho=as.double(rho[ind]),a=as.double(a),b=as.double(b)) } + if (oo$eps < -.5) { + if (oo$eps < -1.5) { ## failure of series in C code + oo$w2 <- oo$w1 <- oo$w2p <- oo$w1p <- oo$w2pp <- rep(NA,length(y)) + } + else warning("Tweedie density may be unreliable - series not fully converged") + } phii <- phi[ind] if (!work.param) { ## transform working param derivatives to p/phi derivs... if (length(dthp1)!=n) dthp1 <- array(dthp1,dim=n) @@ -2960,7 +2966,7 @@ oo$w1 <- oo$w1/phii oo$w2p <- oo$w2p*dthp1i^2 + dthp2[ind] * oo$w1p oo$w1p <- oo$w1p*dthp1i - oo$w2pp <- oo$w2pp*dthp1i/phii ## this appears to be wrong + oo$w2pp <- oo$w2pp*dthp1i/phii } @@ -3014,7 +3020,7 @@ ld[ind,6] <- ld[ind,6] + oo$w2pp ## d2 logf / dphi dp } -if (FALSE) { ## DEBUG disconnetion of density terms +if (FALSE) { ## DEBUG disconnection of density terms ld[ind,1] <- oo$w ## log density ld[ind,2] <- oo$w1 ## d log f / dphi ld[ind,3] <- oo$w2 ## d2 logf / dphi2 diff -Nru mgcv-1.8-27/R/gam.fit4.r mgcv-1.8-28/R/gam.fit4.r --- mgcv-1.8-27/R/gam.fit4.r 2019-02-02 18:06:18.000000000 +0000 +++ mgcv-1.8-28/R/gam.fit4.r 2019-03-20 13:18:38.000000000 +0000 @@ -326,13 +326,6 @@ ## if (!is.null(start)) null.coef <- start - can be on edge of feasible - not good null.eta <- as.numeric(x%*%null.coef + as.numeric(offset)) - #old.pdev <- sum(dev.resids(y, linkinv(null.eta), weights,theta)) + t(null.coef)%*%St%*%null.coef - - #if (!is.null(start)) { ## check it's at least better than null.coef - # pdev <- sum(dev.resids(y, linkinv(x%*%start+as.numeric(offset)), weights,theta)) + t(start)%*%St%*%start - # if (pdev>old.pdev) start <- mustart <- etastart <- NULL - #} - ## call the families initialization code... if (is.null(mustart)) { @@ -372,7 +365,7 @@ mu <- linkinv(eta);etaold <- eta conv <- boundary <- FALSE dd <- dDeta(y,mu,weights,theta,family,0) ## derivatives of deviance w.r.t. eta - w <- dd$Deta2 * .5; + w <- dd$Deta2 * .5 wz <- w*(eta-offset) - .5*dd$Deta z <- (eta-offset) - dd$Deta.Deta2 good <- is.finite(z)&is.finite(w) @@ -907,7 +900,6 @@ if (is.null(weights)) weights <- rep.int(1, nobs) if (is.null(offset)) offset <- rep.int(0, nobs) - ## get log likelihood, grad and Hessian (w.r.t. coefs - not s.p.s) ... llf <- family$ll ll <- llf(y,x,coef,weights,family,offset=offset,deriv=1) @@ -919,6 +911,7 @@ check.deriv <- FALSE; eps <- 1e-5 drop <- NULL;bdrop <- rep(FALSE,q) ## by default nothing dropped perturbed <- 0 ## counter for number of times perturbation tried on possible saddle + for (iter in 1:(2*control$maxit)) { ## main iteration ## get Newton step... if (check.deriv) { ## code for checking derivatives when debugging diff -Nru mgcv-1.8-27/R/gamlss.r mgcv-1.8-28/R/gamlss.r --- mgcv-1.8-27/R/gamlss.r 2019-02-02 19:03:20.000000000 +0000 +++ mgcv-1.8-28/R/gamlss.r 2019-03-20 13:18:38.000000000 +0000 @@ -195,144 +195,6 @@ } # gamlss.etamu -gamlss.gH0 <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=NULL,D=NULL) { -## X[,jj[[i]]] is the ith model matrix. -## lj contains jth derivatives of the likelihood for each datum, -## columns are w.r.t. different combinations of parameters. -## ij is the symmetric array indexer for the jth order derivs... -## e.g. l4[,i4[i,j,l,m]] contains derivatives with -## respect to parameters indexed by i,j,l,m -## d1b and d2b are first and second derivatives of beta w.r.t. sps. -## fh is a factorization of the penalized hessian, while D contains the corresponding -## Diagonal pre-conditioning weights. -## deriv: 0 - just grad and Hess -## 1 - diagonal of first deriv of Hess wrt sps -## 2 - first deriv of Hess wrt sps -## 3 - everything. - K <- length(jj) - p <- ncol(X);n <- nrow(X) - trHid2H <- d1H <- d2H <- NULL ## defaults - - ## the gradient... - lb <- rep(0,p) - for (i in 1:K) { ## first derivative loop - lb[jj[[i]]] <- colSums(l1[,i]*X[,jj[[i]],drop=FALSE]) - } - - ## the Hessian... - lbb <- matrix(0,p,p) - for (i in 1:K) for (j in i:K) { - lbb[jj[[i]],jj[[j]]] <- t(X[,jj[[i]],drop=FALSE])%*%(l2[,i2[i,j]]*X[,jj[[j]],drop=FALSE]) - lbb[jj[[j]],jj[[i]]] <- t(lbb[jj[[i]],jj[[j]]]) - } - - if (deriv>0) { - ## the first derivative of the Hessian, using d1b - ## the first derivates of the coefficients wrt the sps - m <- ncol(d1b) ## number of smoothing parameters - ## stack the derivatives of the various linear predictors on top - ## of each other... - d1eta <- matrix(0,n*K,m) - ind <- 1:n - for (i in 1:K) { - d1eta[ind,] <- X[,jj[[i]],drop=FALSE]%*%d1b[jj[[i]],] - ind <- ind + n - } - } - - if (deriv==1) { - d1H <- matrix(0,p,m) ## only store diagonals of d1H - for (l in 1:m) { - for (i in 1:K) { - v <- rep(0,n);ind <- 1:n - for (q in 1:K) { - v <- v + l3[,i3[i,i,q]] * d1eta[ind,l] - ind <- ind + n - } - d1H[jj[[i]],l] <- colSums(X[,jj[[i]],drop=FALSE]*(v*X[,jj[[i]],drop=FALSE])) - } - } - } ## if deriv==1 - - if (deriv>1) { - d1H <- list() - for (l in 1:m) { - d1H[[l]] <- matrix(0,p,p) - for (i in 1:K) for (j in i:K) { - v <- rep(0,n);ind <- 1:n - for (q in 1:K) { - v <- v + l3[,i3[i,j,q]] * d1eta[ind,l] - ind <- ind + n - } - ## d1H[[l]][jj[[j]],jj[[i]]] <- - d1H[[l]][jj[[i]],jj[[j]]] <- t(X[,jj[[i]],drop=FALSE])%*%(v*X[,jj[[j]],drop=FALSE]) - d1H[[l]][jj[[j]],jj[[i]]] <- t(d1H[[l]][jj[[i]],jj[[j]]]) - } - } - } ## if deriv>1 - - if (deriv>2) { - ## need tr(Hp^{-1} d^2H/drho_k drho_j) - ## First form the expanded model matrix... - VX <- Xe <- matrix(0,K*n,ncol(X)) - ind <- 1:n - for (i in 1:K) { - Xe[ind,jj[[i]]] <- X[,jj[[i]]] - ind <- ind + n - } - ## Now form Hp^{-1} Xe'... - if (is.list(fh)) { ## then the supplied factor is an eigen-decomposition - d <- fh$values;d[d>0] <- 1/d[d>0];d[d<=0] <- 0 - Xe <- t(D*((fh$vectors%*%(d*t(fh$vectors)))%*%(D*t(Xe)))) - } else { ## the supplied factor is a choleski factor - ipiv <- piv <- attr(fh,"pivot");ipiv[piv] <- 1:p - Xe <- t(D*(backsolve(fh,forwardsolve(t(fh),(D*t(Xe))[piv,]))[ipiv,])) - } - ## now compute the required trace terms - d2eta <- matrix(0,n*K,ncol(d2b)) - ind <- 1:n - for (i in 1:K) { - d2eta[ind,] <- X[,jj[[i]],drop=FALSE]%*%d2b[jj[[i]],] - ind <- ind + n - } - trHid2H <- rep(0,ncol(d2b)) - kk <- 0 ## counter for second derivatives - for (k in 1:m) for (l in k:m) { ## looping over smoothing parameters... - kk <- kk + 1 - for (i in 1:K) for (j in 1:K) { - v <- rep(0,n);ind <- 1:n - for (q in 1:K) { ## accumulate the diagonal matrix for X_i'diag(v)X_j - v <- v + d2eta[ind,kk]*l3[,i3[i,j,q]] - ins <- 1:n - for (s in 1:K) { - v <- v + d1eta[ind,k]*d1eta[ins,l]*l4[,i4[i,j,q,s]] - ins <- ins + n - } - ind <- ind + n - } - if (i==j) { - rind <- 1:n + (i-1)*n - VX[rind,jj[[i]]] <- v * X[,jj[[i]]] - } else { - rind1 <- 1:n + (i-1)*n - rind2 <- 1:n + (j-1)*n - VX[rind2,jj[[i]]] <- v * X[,jj[[i]]] - VX[rind1,jj[[j]]] <- v * X[,jj[[j]]] - } - } - trHid2H[kk] <- sum(Xe*VX) - } - } ## if deriv>2 - - list(lb=lb, ## grad w.r.t. coefs - lbb=lbb, ## hess w.r.t. coefs, H - d1H=d1H, ## grad of H wrt sps - ## d2H=d2H, - trHid2H=trHid2H) ## tr(H^{-1}d^2H/dspsp) -} ## end of gamlss.gH0 - - - gamlss.gH <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=NULL,D=NULL) { ## X[,jj[[i]]] is the ith model matrix. ## lj contains jth derivatives of the likelihood for each datum, diff -Nru mgcv-1.8-27/R/mgcv.r mgcv-1.8-28/R/mgcv.r --- mgcv-1.8-27/R/mgcv.r 2019-02-01 22:55:51.000000000 +0000 +++ mgcv-1.8-28/R/mgcv.r 2019-03-18 12:42:14.000000000 +0000 @@ -1080,7 +1080,7 @@ # idea here is that terms are set up in accordance with information given in split$smooth.spec # appropriate basis constructor is called depending on the class of the smooth # constructor returns penalty matrices model matrix and basis specific information - ## sm[[i]] <- smoothCon(split$smooth.spec[[i]],data,knots,absorb.cons,scale.penalty=scale.penalty,sparse.cons=sparse.cons) ## old code + id <- split$smooth.spec[[i]]$id if (is.null(id)||!idLinksBases) { ## regular evaluation sml <- smoothCon(split$smooth.spec[[i]],data,knots,absorb.cons,scale.penalty=scale.penalty, @@ -1126,11 +1126,12 @@ if (is.null(sm[[i]]$L)) Li <- diag(length.S) else Li <- sm[[i]]$L if (length.S > 0) { ## there are smoothing parameters to name - if (length.S == 1) spn <- sm[[i]]$label else { + if (length.S == 1) lspn <- sm[[i]]$label else { Sname <- names(sm[[i]]$S) - if (is.null(Sname)) spn <- paste(sm[[i]]$label,1:length.S,sep="") else - spn <- paste(sm[[i]]$label,Sname,sep="") + lspn <- if (is.null(Sname)) paste(sm[[i]]$label,1:length.S,sep="") else + paste(sm[[i]]$label,Sname,sep="") ## names for all sp's } + spn <- lspn[1:ncol(Li)] ## names for actual working sps } ## extend the global L matrix... @@ -1143,7 +1144,7 @@ cbind(matrix(0,nrow(Li),ncol(L)),Li)) if (length.S > 0) { ## there are smoothing parameters to name sp.names <- c(sp.names,spn) ## extend the sp name vector - lsp.names <- c(lsp.names,spn) ## extend full.sp name vector + lsp.names <- c(lsp.names,lspn) ## extend full.sp name vector } } else { ## it's a repeat id => shares existing sp's L0 <- matrix(0,nrow(Li),ncol(L)) @@ -1153,7 +1154,7 @@ L0[,idx[[id]]$c:(idx[[id]]$c+ncol(Li)-1)] <- Li L <- rbind(L,L0) if (length.S > 0) { ## there are smoothing parameters to name - lsp.names <- c(lsp.names,spn) ## extend full.sp name vector + lsp.names <- c(lsp.names,lspn) ## extend full.sp name vector } } } diff -Nru mgcv-1.8-27/R/plots.r mgcv-1.8-28/R/plots.r --- mgcv-1.8-27/R/plots.r 2018-07-25 20:00:13.000000000 +0000 +++ mgcv-1.8-28/R/plots.r 2019-03-18 12:42:14.000000000 +0000 @@ -32,6 +32,10 @@ } } else if (family=="binomial") { fam$qf <- function(p,mu,wt,scale) { + if (all.equal(wt,ceiling(wt))!=TRUE) { + wt <- ceiling(wt) + warning("non-integer binomial denominator: quantiles incorrect") + } qbinom(p,wt,mu)/(wt + as.numeric(wt==0)) } } else if (family=="Gamma") { @@ -440,7 +444,7 @@ pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL, ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80", shift=0,trans=I,by.resids=FALSE,scheme=0,hcolors=heat.colors(100), - contour.col=3,...) { + contour.col=4,...) { ## plot method function for sos.smooth terms if (scheme>1) return(plot.mgcv.smooth(x,P=P,data=data,label=label,se1.mult=se1.mult,se2.mult=se2.mult, partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2, @@ -702,7 +706,7 @@ pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL, ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80", shift=0,trans=I,by.resids=FALSE,scheme=0,hcolors=heat.colors(50), - contour.col=3,...) { + contour.col=4,...) { ## default plot method for smooth objects `x' inheriting from "mgcv.smooth" ## `x' is a smooth object, usually part of a `gam' fit. It has an attribute ## 'coefficients' containing the coefs for the smooth, but usually these @@ -1551,7 +1555,7 @@ surf.col<-surf.col/(max.z-min.z) surf.col<-round(surf.col*nCol) con.col <-1 - if (color=="heat") { pal<-heat.colors(nCol);con.col<-3;} + if (color=="heat") { pal<-heat.colors(nCol);con.col<-4;} else if (color=="topo") { pal<-topo.colors(nCol);con.col<-2;} else if (color=="cm") { pal<-cm.colors(nCol);con.col<-1;} else if (color=="terrain") { pal<-terrain.colors(nCol);con.col<-2;} diff -Nru mgcv-1.8-27/src/mat.c mgcv-1.8-28/src/mat.c --- mgcv-1.8-27/src/mat.c 2018-11-26 17:07:06.000000000 +0000 +++ mgcv-1.8-28/src/mat.c 2019-03-21 08:56:32.000000000 +0000 @@ -1736,9 +1736,12 @@ 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. + is more Cache friendly with less effort, since the update is then from + the right and operates columnwise. However, code below is column oriented in both + cases, storing the givens rotations as they are computed to facilitate + column orientation when R is upper triangular. + Calls from R should ideally be made from a wrapper called from .Call, + since otherwise copying can be the dominant cost. */ int i,j,n1; @@ -1864,7 +1867,7 @@ /* now construct the next hyperbolic rotation u[j] <-> R[j,j] */ z0 = z / *x; /* sqrt(z^2+R[j,j]^2) */ - if (fabs(z0>=1)) { /* downdate not +ve def */ + if (fabs(z0)>=1) { /* downdate not +ve def */ //Rprintf("j = %d d = %g ",j,z0); if (*n>1) R[1] = -2.0;return; /* signals error */ } diff -Nru mgcv-1.8-27/src/misc.c mgcv-1.8-28/src/misc.c --- mgcv-1.8-27/src/misc.c 2018-09-22 07:21:46.000000000 +0000 +++ mgcv-1.8-28/src/misc.c 2019-03-21 08:56:32.000000000 +0000 @@ -184,6 +184,10 @@ rho=log(phi), and th defines p = (a + b * exp(th))/(exp(th)+1). note, 1 .5||j_max<1) j_max++; + if (x - j_max > .5||j_max<1) j_max++; + if (fabs(j_max-x)>1) { /* index has integer overflowed */ + failed = 1; + break; + } j_max -= j0; /* converted to buffer index */ j = j_max+j0; jalogy = j*alogy[i]; wdW2d2W= wdlogwdp=dWpp=0.0; wi=w1i=w2i=0.0; // 1.0; + /* j_max could be > jal_1 the currently allocated, or outside [j_lo,j_hi] + the currently initialized. In either case we need fill in all the buffer + values from the initialized set to j_max, so we might as well reset j_max + to the appropriate buffer edge. + */ + if (j_max>j_hi) j_max = j_hi; if (j_max jal-1) j_hi = jal-1; /* set j_hi to last element filled */ - if (!ok) { /* need to expand buffer storage*/ + if (!ok) if (jal .5||j_max<1) j_max++; + if (fabs(j_max-x)>1) { /* index has integer overflowed */ + failed = 1; + break; + } j = j_max; onep = 1 - p;onep2 = onep * onep; @@ -573,6 +593,7 @@ //for (j=j_max+j0,jb=j_max;jb<=j_hi;jb++,j++) { // note initially wi etc initialized to 1 and summation starts 1 later incr = 1; lgammaj1 = lgamma((double)j+1); // lgamma(j+1) to be computed by recursion + k=0; while (!ok) { wbj = j * w_base - lgammaj1 - lgamma(-j * alpha); wb1j = -j/onep; @@ -621,6 +642,8 @@ lgammaj1 += -log(j+1); if (wj < wmin||j<1) ok=1; // finished } + k++; + if (k>=jal_lim) ok = series_too_long = 1; /* avoid going on for ever */ } /* end of upsweep/downsweep */ //Rprintf("wdlogwdp = %g\n",wdlogwdp); /* Summation now complete: need to do final transformations */ @@ -632,6 +655,8 @@ w1p[i] = wdlogwdp/wi; } /* end of looping through y */ + if (series_too_long) *eps = -1.0; + if (failed) *eps = -2.0; } /* tweedious2 */