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

Import Upstream version 1.3-19

parent 98c7e43b
Package: mgcv
Version: 1.3-18
Version: 1.3-19
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: GAMs with GCV smoothness estimation and GAMMs by REML/PQL
......@@ -12,4 +12,4 @@ Imports: graphics, stats
Suggests: nlme (>= 3.1-64), MASS (>= 7.2-2)
LazyLoad: yes
License: GPL version 2 or later
Packaged: Thu Jul 20 14:08:34 2006; simon
Packaged: Tue Sep 12 14:57:25 2006; simon
......@@ -6,7 +6,7 @@ anova.gam,coef.pdIdnot,coef.pdTens,
exclude.too.far,extract.lme.cov, extract.lme.cov2,
formXtViX, full.score, formula.gam,fixDependence,fix.family.link,
fix.family.var, gam, gam2derivative, gam2objective,
gam3objective, gamm, gam.check, gam.control, gam.fit2,
gam3objective, gamm, gam.check, gam.control, gam.fit2,gam.fit3,
gam.fit, gam.outer, gamm.setup , influence.gam, interpret.gam,
gam.setup,
gam.side,gam.method,
......@@ -41,9 +41,13 @@ anova.gam,coef.pdIdnot,coef.pdTens,
tensor.prod.model.matrix,tensor.prod.penalties,
uniquecombs, vcov.gam, vis.gam)
importFrom(graphics, plot)
importFrom(stats, predict, residuals, influence, formula, logLik, anova,
cooks.distance,vcov,coef)
importFrom(grDevices,cm.colors,gray,heat.colors,terrain.colors,topo.colors)
importFrom(graphics,axis,box,contour,hist,lines,mtext, par, persp,plot,points,polygon,strheight,strwidth,text)
importFrom(stats,.checkMFClasses,.getXlevels, as.formula, anova,anova.glmlist,coef,cooks.distance,cor,delete.response,family,fitted,formula,
gaussian,glm,influence,lm,logLik,median,model.frame,model.matrix,model.offset,na.pass,napredict,naresid,nlm,optim,
pchisq,pf,pnorm,pt,
predict,printCoefmat ,quantile,qqnorm, reformulate,residuals,rnorm,runif,termplot,terms,terms.formula, uniroot,
var,vcov)
S3method(anova, gam)
S3method(influence, gam)
......
......@@ -6,7 +6,7 @@ rep(1, nobs), start = NULL, etastart = NULL,
mustart = NULL, offset = rep(0, nobs), family = gaussian(),
control = gam.control(), intercept =
TRUE,deriv=TRUE,gamma=1,scale=1,pearson=FALSE,
printWarn=TRUE)
printWarn=TRUE,...)
## deriv, sp, S, H added to arg list.
## need to modify family before call.
{ x <- as.matrix(x)
......@@ -15,7 +15,7 @@ rep(1, nobs), start = NULL, etastart = NULL,
rownames(y)
else names(y)
conv <- FALSE
nobs <- NROW(y)
n <- nobs <- NROW(y) ## n is just to keep codetools happy
nvars <- ncol(x)
EMPTY <- nvars == 0
if (is.null(weights))
......@@ -347,7 +347,7 @@ rep(1, nobs), start = NULL, etastart = NULL,
gam2derivative <- function(lsp,args)
gam2derivative <- function(lsp,args,...)
## Performs IRLS GAM fitting for smoothing parameters given in lsp
## and returns the derivatives of the GCV or UBRE score w.r.t the
## smoothing parameters for the model.
......@@ -356,13 +356,13 @@ gam2derivative <- function(lsp,args)
{ b<-gam.fit2(x=args$X, y=args$y, sp=lsp, S=args$S,rS=args$rS,off=args$off, H=args$H,
offset = args$offset,family = args$family,weights=args$w,deriv=TRUE,
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson,
printWarn=FALSE)
printWarn=FALSE,...)
if (args$scoreType == "GCV") ret <- b$GCV1 else ret <- b$UBRE1
ret
}
gam2objective <- function(lsp,args,printWarn=FALSE)
gam2objective <- function(lsp,args,printWarn=FALSE,...)
## Performs IRLS GAM fitting for smoothing parameters given in lsp
## and returns the GCV or UBRE score for the model.
## args is a list containing the arguments for gam.fit2
......@@ -371,13 +371,13 @@ gam2objective <- function(lsp,args,printWarn=FALSE)
b<-gam.fit2(x=args$X, y=args$y, sp=lsp, S=args$S,rS=args$rS,off=args$off, H=args$H,
offset = args$offset,family = args$family,weights=args$w,deriv=FALSE,
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson,
printWarn=printWarn)
printWarn=printWarn,...)
if (args$scoreType == "GCV") ret <- b$GCV else ret <- b$UBRE
attr(ret,"full.fit") <- b
ret
}
gam3objective <- function(lsp,args)
gam3objective <- function(lsp,args,...)
## Performs IRLS GAM fitting for smoothing parameters given in lsp
## and returns the GCV or UBRE score for the model.
## args is a list containing the arguments for gam.fit2
......@@ -386,7 +386,7 @@ gam3objective <- function(lsp,args)
b<-gam.fit2(x=args$X, y=args$y, sp=lsp, S=args$S,rS=args$rS,off=args$off, H=args$H,
offset = args$offset,family = args$family,weights=args$w,deriv=TRUE,
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson,
printWarn=FALSE)
printWarn=FALSE,...)
if (args$scoreType == "GCV") ret <- b$GCV else ret <- b$UBRE
attr(ret,"full.fit") <- b
if (args$scoreType == "GCV") at <- b$GCV1 else at <- b$UBRE1
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
1.3-19
* numerous changes to NAMESPACE and gamm related functions to pass
codetools checks.
* gam.method() modified to allow GACV as an option for outer GCV
model selection.
* Improved outer iteration added via gdi.c coupled with gam.fit3.
Exact first and second derivatives of GACV/GCV/UBRE/AIC are now
available via new iteration methods. These improve the
speed and reliability of fitting in the generalized case.
* magic.c::mgcv_mmult modified so that all inner loop calculations are
optimal (i.e. inner loop pointers increments are all 1).
* `smooth.construct' functions for "cc" and "cr" smooths now increase `k'
to the minimum possible value (and warn), if it's too low.
* `gam' modified to allow passing of `mustart' etc to gam.fit and
gam.fit2, properly
* `gam' modified to fix a bug whereby fitting in two steps using argument
`G' could fail when some sp's are to be estimated and some fixed.
* an argument `in.out' added to `gam' to allow user initialization of
smoothing parameters when using `outer' iteration in the generalized
case. This can speed up analyses that rely on several refits of the same
model.
1.3-18
* gamm modifed so that weights dealt with properly if lme type varFunc
......
......@@ -12,7 +12,7 @@ outer iteration (see \code{\link{gam.outer}}) when the the outer method is
\code{"nlm.fd"} (see \code{\link{gam.method}}).
}
\usage{
full.score(sp,G,family,control,gamma,pearson)
full.score(sp,G,family,control,gamma,pearson,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......@@ -28,9 +28,12 @@ full.score(sp,G,family,control,gamma,pearson)
\item{pearson}{whether the GCV/UBRE score should be based on the Pearson
statistic or the deviance. The deviance is usually preferable.}
\item{...}{other arguments, typically for passing on to \code{gam.fit}.}
}
\value{ The value of the GCV/UBRE score, with attribute \code{"full.gam.object"}
which is the full object returned by \code{\link{gam.fit}}.
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}}
......
......@@ -38,7 +38,7 @@ variables, and performance iteration for smoothing parameter estimation (see
gam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
na.action,offset=NULL,control=gam.control(),method=gam.method(),
scale=0,knots=NULL,sp=NULL,min.sp=NULL,H=NULL,gamma=1,
fit=TRUE,G=NULL,...)
fit=TRUE,G=NULL,in.out,...)
}
%- maybe also `usage' for other objects documented here.
\details{
......@@ -289,10 +289,18 @@ such a multiplier to be supplied if fit method is \code{"magic"}.}
\item{fit}{If this argument is \code{TRUE} then \code{gam} sets up the model and fits it, but if it is
\code{FALSE} then the model is set up and an object \code{G} is returned which is the output from
\code{\link{gam.setup}} plus some extra items required to complete the GAM fitting process.}
\item{G}{Usually \code{NULL}, but may contain the object returned by a previous call to \code{gam} with
\code{fit=FALSE}, in which case all other arguments are ignored except for \code{gamma}, \code{family}, \code{control} and \code{fit}.}
\code{fit=FALSE}, in which case all other arguments are ignored except for
\code{gamma}, \code{in.out}, \code{control}, \code{method} and \code{fit}.}
\item{in.out}{optional list for initializing outer iteration. If supplied
then this must contain two elements: \code{sp} should be an array of
initialization values for all smoothing parameters (there must be a value for
all smoothing parameters, whether fixed or to be estimated, but those for
fixed s.p.s are not used); \code{scale} is the typical scale of the GCV/UBRE function,
for passing to the outer optimizer.}
\item{...}{further arguments for
passing on e.g. to \code{gam.fit}} }
passing on e.g. to \code{gam.fit}} (such as \code{mustart}). }
\value{
If \code{fit == FALSE} the function returns a list \code{G} of items needed to
......
......@@ -13,8 +13,8 @@ choise of fitting method, see \code{\link{gam.method}}.
gam.control(irls.reg=0.0,epsilon = 1e-06, maxit = 100,globit = 20,
mgcv.tol=1e-7,mgcv.half=15,nb.theta.mult=10000, trace = FALSE,
rank.tol=.Machine$double.eps^0.5,absorb.cons=TRUE,
max.tprs.knots=5000,nlm=list(),optim=list(),
outerPIsteps=4)
max.tprs.knots=5000,nlm=list(),optim=list(),newton=list(),
outerPIsteps=1)
}
\arguments{
\item{irls.reg}{For most models this should be 0. The iteratively re-weighted least squares method
......@@ -77,7 +77,10 @@ examples in \code{\link{gam}}, rather than simply increasing this default.}
used for outer estimation of smoothing parameters. See details.}
\item{optim}{list of control parameters to pass to \code{\link{optim}} if this
is used for outer estimation of smoothin parameters. See details.}
is used for outer estimation of smoothing parameters. See details.}
\item{newton}{list of control parameters to pass to default Newton optimizer
used for outer estimation of log smoothing parameters. See details.}
\item{outerPIsteps}{The number of performance interation steps used to
initialize outer iteration. Less than 1 means
......@@ -100,6 +103,16 @@ calculations should be checked numerically - defaults to \code{FALSE}. Any of
these which are not supplied and named in the list are set to their default
values.
Outer iteration using \code{newton} is controlled by the list \code{newton}
with the following elements : \code{conv.tol} (default
1e-6) is the relative convergence tolerance; \code{maxNstep} is the maximum
length allowed for an element of the Newton search direction (default 5);
\code{maxSstep} is the maximum length allowed for an element of the steepest
descent direction (only used if Newton fails - default 2); \code{maxHalf} is
the maximum number of step halvings to permit before giving up (default 30);
\code{use.svd} determines whether to use an SVD step or QR step for the part
of fitting that determines problem rank, etc. (defaults to \code{FALSE}).
Outer iteration using \code{\link{optim}} is controlled using list
\code{optim}, which currently has one element: \code{factr} which takes
default value 1e7.
......
\name{gam.fit2}
\alias{gam.fit2}
\alias{gam.fit3}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{P-IRLS GAM estimation with GCV \& UBRE derivative calculation}
\description{Estimation of GAM smoothing parameters is most stable if
......@@ -7,20 +8,30 @@ optimization of the UBRE or GCV score is outer to the penalized iteratively
re-weighted least squares scheme used to estimate the model given smoothing
parameters.
This routine estimates a GAM given log smoothing paramaters, and evaluates the
first derivative of the GCV and UBRE scores of the model with respect to the
These routines estimates a GAM given log smoothing paramaters, and evaluate
derivatives of the GCV and UBRE scores of the model with respect to the
log smoothing parameters. Calculation of exact derivatives is generally faster
than approximating them by finite differencing, as well as generally improving
the reliability of GCV/UBRE score minimization.
Not normally called directly, but rather a service routine for \code{\link{gam}}.
\code{gam.fit2} evaluates first derivatives, by accumulating them as P-IRLS
progresses. \code{gam.fit3} uses a more efficient approach in which the P-IRLS
is first run to convergence, and only then are the derivatives evaluated by a
separate iteration. \code{gam.fit3} can evaluate second as well as first derivatives.
Not normally called directly, but rather service routines for \code{\link{gam}}.
}
\usage{
gam.fit2(x, y, sp, S=list(),rS=list(),off, H=NULL,
weights = rep(1, nobs), start = NULL, etastart = NULL,
mustart = NULL, offset = rep(0, nobs), family = gaussian(),
control = gam.control(), intercept = TRUE,deriv=TRUE,
gamma=1,scale=1,pearson=FALSE,printWarn=TRUE)
gamma=1,scale=1,pearson=FALSE,printWarn=TRUE,...)
gam.fit3(x, y, sp, S=list(),rS=list(),off, H=NULL,
weights = rep(1, nobs), start = NULL, etastart = NULL,
mustart = NULL, offset = rep(0, nobs), family = gaussian(),
control = gam.control(), intercept = TRUE,deriv=2,use.svd=TRUE,
gamma=1,scale=1,printWarn=TRUE,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......@@ -60,7 +71,13 @@ as possible, but as many rows as there are parameters.}
\code{FALSE}}
\item{deriv}{ Should derivatives of the GCV and UBRE scores be calculated?
\code{TRUE} or \code{FALSE}}
\code{TRUE} or \code{FALSE} for \code{gam.fit2}, but 0, 1 or 2,
indicating the maximum order of differentiation to apply, for
\code{gam.fit3}.}
\item{use.svd}{Should the algorithm use SVD (\code{TRUE}) or the cheaper QR
(\code{FALSE}) as the second matrix decomposition of the final
derivative/coefficient evaluation method? Only used by \code{gam.fit3}. }
\item{gamma}{The weight given to each degree of freedom in the GCV and UBRE
scores can be varied (usually increased) using this parameter.}
......@@ -69,11 +86,13 @@ scores can be varied (usually increased) using this parameter.}
\item{pearson}{The GCV/UBRE score can be based either on the Pearson statistic
or the deviance. The latter is generally to be preferred, as it is less prone
to severe undersmoothing.}
to severe undersmoothing. Only used by \code{gam.fit2}.}
\item{printWarn}{Set to \code{FALSE} to suppress some warnings. Useful in
order to ensure that some warnings are only printed if they apply to the final
fitted model, rather than an intermediate used in optimization.}
\item{...}{Other arguments: ignored.}
}
\details{ This routine is basically \code{\link{glm.fit}} with some
modifications to allow (i) for quadratic penalties on the log likelihood;
......
......@@ -9,7 +9,8 @@ estimation criterion for a gam.
It is used to set argument \code{method} of \code{\link{gam}}.
}
\usage{
gam.method(am="magic",gam="outer",outer="nlm",pearson=FALSE,family=NULL)
gam.method(am="magic",gam="outer",outer="newton",gcv="deviance",
family=NULL)
}
\arguments{
\item{am}{ Which method to use for a pure additive model (i.e. identity link,
......@@ -26,17 +27,28 @@ engine. \code{"perf.outer"} for \code{\link{magic}} based performance
iteration followed by outer iteration (see details). \code{"outer"} for pure
outer iteration.}
\item{outer}{ The optimization approach to use for outer
iteration. \code{"nlm"} to use \code{\link{nlm}} with exact first derivatives
\item{outer}{ The optimization approach to use to optimize log smoothing
parameters by outer
iteration.\code{"newton"} for modified Newton method backed up by steepest
descent, based on exact first and second derivatives. \code{"nlm"} to use
\code{\link{nlm}} with exact first derivatives
to optimize the smoothness selection criterion. \code{"nlm.fd"} to use
\code{nlm} with finite differenced first dericatives (slower and less
\code{nlm} with finite differenced first derivatives (slower and less
reliable). \code{"optim"} to use the \code{"L-BFGS-B"} quasi-Newton method
option of routine \code{optim}, with exact first derivatives.}
\item{pearson}{When using any sort of outer iteration, there is a choise
between using GCV/UBRE scores based on the Pearson statistic, or based on the
deviance. The Pearson versions result in smoother models, but this is often
oversmoothing, particularly for `low count' data and binary regression.}
\item{gcv}{One of \code{"deviance"}, \code{"GACV"} or \code{"pearson"},
specifying the flavour of GCV to use with outer iteration. \code{"deviance"}
simply replaces the residual sum of squares term in a GCV score with the
deviance, following Hastie and Tibshirani (1990, section 6.9). \code{"GACV"}
(only available with outer method \code{"newton"})
uses a varient of Gu and Xiang's (2001) generalized approximate cross
validation, modified to deal with arbitrary link-error
combinations. \code{"pearson"} uses a suggestion from the literature of
replacing the RSS in a GCV score with the Pearson statistic: this is a very
poor choise (see Wood, 2006, section 4.5.4). This latter option is not
available with outer methods \code{"newton"} and \code{"nlm"}.
}
\item{family}{The routine is called by \code{\link{gam}} to check the supplied
method argument. In this circumstance the family argument is passed, to check
......@@ -88,28 +100,23 @@ particularly with low n binomial data, or low mean counts. Hence the default
is to use (i).
Several alternative optimisation methods can be used for outer
optimization. \code{nlm} can be used with finite differenced first
derivatives. This is not ideal theoretically, since it is possible for the
finite difference estimates of derivatives to be very badly in error on rare
optimization.Usually the fastest and most
reliable approach is to use a modified Newton optimizer with exact first and
second derivatives, and this is the default. \code{nlm} can be used with
finite differenced first derivatives. This is not ideal theoretically, since
it is possible for the finite difference estimates of derivatives to be very
badly in error on rare
occasions when the P-IRLS convergence tolerance is close to being matched
exactly, so that two components of a finite differenced derivative require
different numbers of iterations of P-IRLS in their evaluation. An alternative
is provided in which \code{nlm} uses numerically exact first derivatives, this
is faster and less problematic than the other scheme. Finally, a quasi-Newton
scheme with exact derivtives can be used instead, based on \code{optim}. In practice this usually
seems to be slower than the \code{nlm} method.
It is possible to iterate the performance iteration to convergence and then
improve the smoothing parameter estimates using outer iteration: only a few steps are
usually required in the outer iteration in this case, so it may be quite
efficient, but it is not recommended if the performance iteration itself is
non-convergent. When using `pure' outer iteration, a single step of the
performance iteration is in fact taken first, to obtain estimates of the scale
of the GCV/UBRE objective: starting values for the smoothing parameters are
obtained using \code{\link{initial.sp}}.
is faster and less problematic than the other scheme. An further alternative is to use a quasi-Newton
scheme with exact derivtives, based on \code{optim}. In practice this usually
seems to be slower than the \code{nlm} method.
In summary: performance iteration is fast, but can fail to converge. Outer
iteration is slower, but more reliable. At present only performance iteration
iteration is a little slower, but more reliable. At present only performance iteration
is available for negative binomial families.
}
......
......@@ -21,7 +21,7 @@ the Pearson statistic: see \code{\link{gam.method}}.
Not normally called directly, but rather a service routine for \code{\link{gam}}.
}
\usage{
gam.outer(lsp,fscale,family,control,method,gamma,G)
gam.outer(lsp,fscale,family,control,method,gamma,G,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......@@ -42,6 +42,7 @@ the optimization method to use.}
\item{G}{List produced by \code{\link{gam.setup}}, containing most of what's
needed to actually fit GAM.}
\item{...}{other arguments, typically for passing on to \code{gam.fit2} (ultimately).}
}
\details{
Estimation of smoothing parameters by optimizing GCV scores obtained at
......
......@@ -12,9 +12,9 @@ smoothing parameters, in a manner sutable for use by \code{\link{optim}} or \cod
Not normally called directly, but rather service routines for \code{\link{gam.outer}}.
}
\usage{
gam2objective(lsp,args,printWarn=FALSE)
gam2derivative(lsp,args)
gam3objective(lsp,args)
gam2objective(lsp,args,printWarn=FALSE,...)
gam2derivative(lsp,args,...)
gam3objective(lsp,args,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......@@ -25,6 +25,8 @@ gam3objective(lsp,args)
\item{printWarn}{Should \code{\link{gam.fit2}} print some warnings? Used to
suppress warnings that are only of interest for the final fitted model, until
that model is reached.}
\item{...}{Other arguments for passing to \code{gam.fit2}.}
}
\details{ \code{gam2objective} and \code{gam2derivative} are functions suitable
for calling by \code{\link{optim}}, to evaluate the GCV/UBRE score and it's
......
......@@ -29,8 +29,8 @@ To use this function effectively it helps to be quite familiar with the use of
\usage{
gamm(formula,random=NULL,correlation=NULL,family=gaussian(),
data=list(),weights=NULL,subset=NULL,na.action,knots=NULL,
control=lmeControl(niterEM=0,optimMethod="L-BFGS-B"),niterPQL=20,
verbosePQL=TRUE,method="ML",...)
control=nlme::lmeControl(niterEM=0,optimMethod="L-BFGS-B"),
niterPQL=20,verbosePQL=TRUE,method="ML",...)
}
\arguments{
......
......@@ -9,8 +9,9 @@
\usage{
plot.gam(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scale=-1,
n=100,n2=40,pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,
ylab=NULL,main=NULL,ylim=NULL,xlim=NULL,too.far=0.1,all.terms=FALSE,
shade=FALSE,shade.col="gray80",shift=0,trans=I,...)
ylab=NULL,main=NULL,ylim=NULL,xlim=NULL,too.far=0.1,
all.terms=FALSE,shade=FALSE,shade.col="gray80",
shift=0,trans=I,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......
......@@ -503,3 +503,13 @@ msgstr ""
msgid "weights must be like glm weights for generalized case"
msgstr ""
msgid "Unkwown flavour of GCV"
msgstr ""
msgid "GACV only supported with newton optimization, GCV type reset"
msgstr ""
msgid "Pearson based GCV is unsupported for newton or nlm outer methods, reset"
msgstr ""
This diff is collapsed.
/* Symbol registration initialization: original provided by Brian Ripley.
Anything called from R should be registered here. (See also NAMESPACE:1)
Anything called from R should be registered here (and declared in mgcv.h).
(See also NAMESPACE:1)
*/
#include <R.h>
#include <Rinternals.h>
......@@ -18,6 +19,10 @@ R_CMethodDef CEntries[] = {
{"predict_tprs", (DL_FUNC) &predict_tprs, 12},
{"MinimumSeparation", (DL_FUNC) &MinimumSeparation, 7},
{"magic", (DL_FUNC) &magic, 16},
{"mgcv_mmult", (DL_FUNC) &mgcv_mmult,8},
{"gdi",(DL_FUNC) &gdi,36},
{"R_cond",(DL_FUNC) &R_cond,5} ,
{"pls_fit",(DL_FUNC)&pls_fit,10},
{NULL, NULL, 0}
};
......
......@@ -64,11 +64,11 @@ void mgcv_AtA(double *AA,double *A,int *q,int *n)
matrix(.C("mgcv_AtA",as.double(rep(0,q*q)),as.double(A),as.integer(q),
as.integer(n),PACKAGE="mgcv")[[1]],q,q)
*/
{ double xx,*p,*p1,*p2,*p3;
{ double xx,*p,*p1,*p2,*p3,*p4;
int nq,i,j;
nq= *n * *q;
for (i=0,p=A;i < *q;p+= *n,i++) for (j=i,p1=p;j< *q;p1+= *n,j++)
{ for (xx=0.0,p2=p,p3=p1;p2<p + *n;p2++,p3++) xx += *p2 * *p3;
{ for (xx=0.0,p2=p,p3=p1,p4=p+ *n;p2<p4;p2++,p3++) xx += *p2 * *p3;
AA[i * *q + j] = AA[ j * *q + i]=xx;
}
}
......@@ -77,54 +77,77 @@ void mgcv_AtA(double *AA,double *A,int *q,int *n)
void mgcv_mmult(double *A,double *B,double *C,int *bt,int *ct,int *r,int *c,int *n)
/* Forms r by c product of B and C, transposing each according to bt and ct.
n is the common dimension of the two matrices, which are stored in R
default column order form:
r<-3;c<-4;n<-5
default column order form. Algorithm is inner loop optimized in each case
(i.e. inner loops only update pointers by 1, rather than jumping).
r<-1000;c<-1000;n<-1000
A<-matrix(0,r,c);B<-matrix(rnorm(r*n),n,r);C<-matrix(rnorm(n*c),c,n)
system.time(
A<-matrix(.C("mgcv_mmult",as.double(A),as.double(B),as.double(C),as.integer(1),as.integer(1),
as.integer(r),as.integer(c),as.integer(n))[[1]],r,c)
t(B)%*%t(C)-A
as.integer(r),as.integer(c),as.integer(n))[[1]],r,c))
range(as.numeric(t(B)%*%t(C)-A))
A<-matrix(0,r,c);B<-matrix(rnorm(r*n),n,r);C<-matrix(rnorm(n*c),n,c)
system.time(
A<-matrix(.C("mgcv_mmult",as.double(A),as.double(B),as.double(C),as.integer(1),as.integer(0),
as.integer(r),as.integer(c),as.integer(n))[[1]],r,c)
t(B)%*%C-A
as.integer(r),as.integer(c),as.integer(n))[[1]],r,c))
range(as.numeric(t(B)%*%C-A))
A<-matrix(0,r,c);B<-matrix(rnorm(r*n),r,n);C<-matrix(rnorm(n*c),c,n)
system.time(
A<-matrix(.C("mgcv_mmult",as.double(A),as.double(B),as.double(C),as.integer(0),as.integer(1),
as.integer(r),as.integer(c),as.integer(n))[[1]],r,c)
B%*%t(C)-A
as.integer(r),as.integer(c),as.integer(n))[[1]],r,c))
range(as.numeric(B%*%t(C)-A))
A<-matrix(0,r,c);B<-matrix(rnorm(r*n),r,n);C<-matrix(rnorm(n*c),n,c)
system.time(
A<-matrix(.C("mgcv_mmult",as.double(A),as.double(B),as.double(C),as.integer(0),as.integer(0),
as.integer(r),as.integer(c),as.integer(n))[[1]],r,c)
B%*%C-A
as.integer(r),as.integer(c),as.integer(n))[[1]],r,c))
range(as.numeric(B%*%C-A))
*/
{ double xx,*bp,*cp,*bp1,*cp1,*ap;
int br,cr,i;
{ double xx,*bp,*cp,*bp1,*cp1,*cp2,*cp3,*ap,*ap1;
int br,cr,i,j;
if (*bt)
{ if (*ct) /* A=B'C' */
{ br= *n;cr= *c;
for (ap=A,cp1=C;cp1<C+ *c;cp1++) for (bp1=B;bp1<B+ *r*br;bp1+=br,ap++)
{ for (xx=0.0,bp=bp1,cp=cp1;bp<bp1 + *n;bp++,cp+=cr) xx += *bp * *cp;/* B[k+br*i]*C[j+cr*k];*/
*ap=xx;
}
} else /* A=B'C */
{ br= *n;cr= *n;
for (ap=A,cp1=C;cp1< C+ *c*cr;cp1+=cr) for (bp=B,i=0;i< *r;i++,ap++)
{ for (xx=0.0,cp=cp1;cp< cp1+ *n;cp++,bp++) xx += *cp * *bp; /* B[k+br*i]*C[k+cr*j];*/
{ /* this one is really awkward: have to use first row of C' as working storage
for current A row; move most slowly through B */
for (i=0;i<*r;i++) /* update A *row-wise*, one full sweep of C per i */
{ cp1 = C + *c;
/* back up row 1 of C' (in A), and initialize current row of A (stored in row 1 of C') */
xx = *B;
for (ap=A,cp=C;cp<cp1;cp++,ap+= *r) { *ap = *cp; *cp *= xx;}
B++;cp2=cp1;
for (j=1;j< *n;j++,B++)
for (xx= *B,cp=C;cp<cp1;cp++,cp2++) *cp += xx * *cp2;
/* row i of A is now in row 1 of C', need to swap back */
for (ap=A,cp=C;cp<cp1;cp++,ap+= *r) { xx = *ap; *ap = *cp; *cp = xx;}
A++;
}
} else /* A=B'C - easiest case: move most slowly through A*/
{ br= *n;cr= *n;cp2 = C + *c * cr;
for (ap=A,cp1=C;cp1< cp2;cp1+=cr) for (bp=B,i=0;i< *r;i++,ap++)
{ for (xx=0.0,cp=cp1,cp3=cp1+ *n;cp< cp3;cp++,bp++) xx += *cp * *bp; /* B[k+br*i]*C[k+cr*j];*/
*ap=xx;
}
}
} else
{ if (*ct) /* A=BC' */
{ br= *r;cr= *c;
for (ap=A,cp1=C;cp1< C + *c;cp1++) for (bp1=B;bp1<B+ *r;bp1++,ap++)
{ for (xx=0.0,bp=bp1,cp=cp1;bp< bp1+ *n*br;bp+=br,cp+=cr) xx += *bp * *cp;/* B[i+br*k]*C[j+cr*k];*/
*ap=xx;
{ if (*ct) /* A=BC' - update A column-wise, move most slowly through C (but in big steps)*/
{ cp = C;
for (j=0;j< *c;j++) /* go through columns of A, one full B sweep per j */
{ ap1 = A + *r;C=cp;xx = *C;bp=B;
for (ap=A;ap<ap1;ap++,bp++) *ap = xx * *bp ;
C += *c;
for (i=1;i< *n;i++,C+= *c)
for (xx=*C,ap=A;ap<ap1;ap++,bp++) *ap += xx * *bp;
A=ap1;cp++;
}
} else /* A=BC */
{ br= *r;cr= *n;
for (ap=A,cp1=C;cp1< C + *c*cr;cp1+=cr) for (bp1=B;bp1<B+ *r;ap++,bp1++)
{ for (xx=0.0,cp=cp1,bp=bp1;cp<cp1+ *n;cp++,bp+=br) xx += *bp * *cp; /* B[i+br*k]*C[k+cr*j];*/
*ap=xx;
} else /* A=BC - update A column-wise, moving most slowly through C */
{ for (j=0;j< *c;j++) /* go through columns of A, one full B sweep per j */
{ ap1 = A + *r;xx = *C;bp=B;
for (ap=A;ap<ap1;ap++,bp++) *ap = xx * *bp ;
C++;
for (i=1;i< *n;i++,C++)
for (xx=*C,ap=A;ap<ap1;ap++,bp++) *ap += xx * *bp;
A=ap1;
}
}
}
......
......@@ -9,10 +9,22 @@ void update_beta(double *X,double *Sr,double *rS,double *theta,double *w,
int *m, int *n,int *q, int *get_trA,int *deriv,
double *rank_tol,double *beta, double *trA, double *beta1,
double *trA1,double *rV,int *rank_est);
void magic(double *y,double *X,double *sp,double *def_sp,double *S,double *H,
double *gamma,double *scale, int *control,int *cS,double *rank_tol,
double *tol,double *b,double *rV,double *norm_const,int *n_score);
void gdi(double *X,double *E,double *rS,
double *sp,double *z,double *w,double *mu, double *eta, double *y,
double *p_weights,double *g1,double *g2,double *g3,double *V0,
double *V1,double *V2,double *beta,double *D1,double *D2,
double *P0,double *P1,double *P2,double *trA,
double *trA1,double *trA2,double *rV,double *rank_tol,double *conv_tol, int *rank_est,
int *n,int *q, int *M,int *Encol,int *rSncol,int *deriv,int *use_svd);
void pls_fit(double *y,double *X,double *w,double *E,int *n,int *q,int *cE,double *eta,
double *penalty,double *rank_tol);
/* various service routines */
void RQT(double *A,int *r,int *c);
void RuniqueCombs(double *X,int *r, int *c);
......@@ -30,7 +42,7 @@ void mgcv_mmult(double *A,double *B,double *C,int *bt,int *ct,int *r,int *c,int
void mgcv_svd_full(double *x,double *vt,double *d,int *r,int *c);
void mgcv_symeig(double *A,double *ev,int *n,int *use_dsyevd);
void mroot(double *A,int *rank,int *n);
void R_cond(double *R,int *r,int *c,double *work,double *Rcondition);
/* basis constructor/prediction routines*/
......
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