[Lme4-commits] r1416 - in pkg/lme4Eigen: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 14 21:54:51 CEST 2011
Author: dmbates
Date: 2011-10-14 21:54:51 +0200 (Fri, 14 Oct 2011)
New Revision: 1416
Modified:
pkg/lme4Eigen/DESCRIPTION
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/predModule.cpp
pkg/lme4Eigen/src/predModule.h
pkg/lme4Eigen/src/respModule.cpp
pkg/lme4Eigen/src/respModule.h
Log:
nlmer added. Increment calculation for PWRSS works but not yet incorporated in R code.
Modified: pkg/lme4Eigen/DESCRIPTION
===================================================================
--- pkg/lme4Eigen/DESCRIPTION 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/DESCRIPTION 2011-10-14 19:54:51 UTC (rev 1416)
@@ -1,6 +1,6 @@
Package: lme4Eigen
-Version: 0.9996875-3
-Date: 2011-10-01
+Version: 0.9996875-4
+Date: 2011-10-14
Title: Linear mixed-effects models using Eigen and S4
Author: Douglas Bates <bates at stat.wisc.edu>,
Martin Maechler <maechler at R-project.org> and
@@ -10,7 +10,7 @@
The models and their components are represented using S4 classes and
methods. The core computational algorithms are implemented using the
Eigen C++ library for numerical linear algebra and RcppEigen "glue".
-Depends: methods, R(>= 2.14.0), Matrix(>= 0.9996875-1), lattice, stats, minqa(>= 1.1.15), Rcpp(>= 0.9.5), RcppEigen(>= 0.1.2)
+Depends: methods, R(>= 2.14.0), Matrix(>= 1.0), lattice, stats, minqa(>= 1.1.15)
LinkingTo: Rcpp, RcppEigen
Imports: graphics, grid, nlme, splines, MatrixModels(>= 0.2.0)
Suggests: mlmRev, MEMSS, boot, MASS, RUnit
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/NAMESPACE 2011-10-14 19:54:51 UTC (rev 1416)
@@ -7,11 +7,10 @@
import("grid")
import("lattice")
import("splines")
-import("Rcpp")
+#import("Rcpp")
importFrom("graphics", plot)
importFrom("nlme", fixef, ranef, VarCorr)
-#importFrom("stats4", AIC, BIC, logLik)# so S4 methods are used!
importFrom("stats"
, AIC
, BIC
@@ -83,7 +82,7 @@
## , update
)
-#exportPattern("^[^\\.]")
+# exportPattern("^[^\\.]")
# and the rest (S3 generics; regular functions):
export(
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/R/AllClass.R 2011-10-14 19:54:51 UTC (rev 1416)
@@ -154,7 +154,7 @@
},
updateXwts = function(wts) {
'update Ut and V from Zt and X using X weights'
- .Call(merPredDupdateXwts, ptr, as.numeric(wts))
+ .Call(merPredDupdateXwts, ptr, wts)
}
)
)
@@ -167,8 +167,7 @@
Ptr = "externalptr",
ptr = function (value) { # generates external pointer if necessary
if (!missing(value)) stop("ptr field cannot be assigned")
- try(
- {
+ try({
if (.Call(isNullExtPtr, Ptr)) {
Ptr <<- .Call(lm_Create, y)
if (length(Offset)) .Call(lm_setOffset, Ptr, Offset)
@@ -182,16 +181,14 @@
if (missing(value))
return(tryCatch(.Call(lm_offset, ptr), error=function(e)Offset))
## lm_setOffset checks value's length, etc.
- .Call(lm_setOffset, ptr, value <- as.numeric(value))
- Offset <<- value
+ Offset <<- .Call(lm_setOffset, ptr, value <- as.numeric(value))
},
Weights = "numeric",
weights = function(value) {
if (missing(value))
return(tryCatch(.Call(lm_weights, ptr), error=function(e)Weights))
## lm_setWeights checks value's length, non-negativity, etc.
- .Call(lm_setWeights, ptr, value <- as.numeric(value))
- Weights <<- value
+ Weights <<- .Call(lm_setWeights, ptr, value <- as.numeric(value))
}),
methods =
list(
@@ -227,28 +224,25 @@
setRefClass("lmerResp",
fields=
list(
- ptr = function (value) {
+ ptr=function (value) {
if (!missing(value)) stop("ptr field cannot be assigned")
- try(
+ try({
if (.Call(isNullExtPtr, Ptr)) {
if (!length(y)) stop("field y must have positive length")
Ptr <<- .Call(lmer_Create, y)
if (length(Offset)) .Call(lm_setOffset, Ptr, Offset)
if (length(Weights)) .Call(lm_setWeights, Ptr, Weights)
- }, silent=TRUE)
+ }
+ }, silent=TRUE)
Ptr
},
REML="integer",
reml=function(value) {
if (missing(value))
- return(tryCatch(.Call(lmer_REML, ptr), error=function(e)0L))
+ return(tryCatch(.Call(lmer_REML, ptr), error=function(e)REML))
stopifnot(length(value <- as.integer(value)) > 0L,
(value <- value[1]) >= 0L)
-### FIXME: Change the cls_set<whatever> C functions to return the
-### value that was set so these can be collapsed to one-liners like
-### REML <<- .Call(lmer_setREML, ptr, value) ?
- .Call(lmer_setREML, ptr, value)
- REML <<- value
+ REML <<- .Call(lmer_setREML, ptr, value)
}),
contains="lmResp",
methods=
@@ -270,16 +264,17 @@
family="family",
ptr = function (value) {
if (!missing(value)) stop("ptr field cannot be assigned")
- try (
- if (.Call(isNullExtPtr, Ptr)) {
- if (!length(y))
- stop("field y must have positive length")
- if (!length(family$family))
- stop("field family must be initialized")
- Ptr <<- .Call(glm_Create, family, y)
- if (length(Offset)) .Call(lm_setOffset, Ptr, Offset)
- if (length(Weights)) .Call(lm_setWeights, Ptr, Weights)
- }, silent=TRUE)
+ try({
+ if (.Call(isNullExtPtr, Ptr)) {
+ if (!length(y))
+ stop("field y must have positive length")
+ if (!length(family$family))
+ stop("field family must be initialized")
+ Ptr <<- .Call(glm_Create, family, y)
+ if (length(Offset)) .Call(lm_setOffset, Ptr, Offset)
+ if (length(Weights)) .Call(lm_setWeights, Ptr, Weights)
+ }
+ }, silent=TRUE)
Ptr
},
N="integer",
@@ -361,13 +356,14 @@
nlsResp <-
setRefClass("nlsResp",
fields=
- list(nlmod="language",
+ list(N="integer",
+ nlmod="formula",
nlenv="environment",
pnames="character",
ptr = function (value) {
if (!missing(value)) stop("ptr field cannot be assigned")
try (if (.Call(isNullExtPtr, Ptr))
- Ptr <<- .Call(nls_Create, nlmod, nlenv, y),
+ Ptr <<- .Call(nls_Create, y, nlmod[[2]], nlenv, pnames, N),
silent=TRUE)
Ptr
}),
@@ -384,6 +380,8 @@
})
)
+nlsResp$lock("nlmod", "nlenv", "pnames")
+
glmFamily <- # used in tests of family definitions
setRefClass("glmFamily",
fields=
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/R/lmer.R 2011-10-14 19:54:51 UTC (rev 1416)
@@ -9,7 +9,7 @@
if(length(l... <- list(...))) {
if (!is.null(l...$family)) { # call glmer if family specified
mc[[1]] <- as.name("glmer")
- return( eval(mc, parent.frame()) )
+ return(eval(mc, parent.frame()) )
}
## Check for method argument which is no longer used
if (!is.null(method <- l...$method)) {
@@ -45,8 +45,9 @@
form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
X <- model.Matrix(form, fr, contrasts, sparse = FALSE, row.names = FALSE) ## sparseX not yet
p <- ncol(X)
- pp <- new("merPredD", X=X, Zt=reTrms$Zt, Lambdat=reTrms$Lambdat, Lind=reTrms$Lind,
- theta=reTrms$theta)
+ if ((qrX <- qr(X))$rank < p)
+ stop(gettextf("rank of X = %d < ncol(X) = %d", qrX$rank, p))
+ pp <- do.call(merPredD$new, c(list(X=X), reTrms[c("Zt","theta","Lambdat","Lind")]))
resp <- mkRespMod2(fr)
if (REML) resp$REML <- p
@@ -83,7 +84,7 @@
mkdevfun <- function(pp, resp) {
if (is(resp, "lmerResp"))
- return (function(theta) .Call(lmer_Deviance, pp$ptr, resp$ptr, theta))
+ return(function(theta) .Call(lmer_Deviance, pp$ptr, resp$ptr, theta))
stop("unknown response type: ", class(resp))
}
@@ -142,7 +143,7 @@
form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
X <- model.Matrix(form, fr, contrasts, sparse = FALSE, row.names = FALSE) ## sparseX not yet
p <- ncol(X)
- pp <- merPredD$new(X=X, Zt=reTrms$Zt, Lambdat=reTrms$Lambdat, Lind=reTrms$Lind, theta=reTrms$theta)
+ pp <- do.call(merPredD$new, c(list(X=X), reTrms[c("Zt","theta","Lambdat","Lind")]))
# response module
resp <- mkRespMod2(fr, family=family)
# initial step from working response
@@ -233,14 +234,14 @@
lower=reTrms$lower, devcomp=list(cmp=cmp, dims=dims), pp=pp, resp=resp)
}## {glmer}
-##' Create an lmerResp, glmerResp or (later) nlmerResp instance
+##' Create an lmerResp, glmResp or nlsResp instance
##'
-##' @title Create a [ng]lmerResp instance
+##' @title Create an lmerResp, glmResp or nlsResp instance
##' @param fr a model frame
-##' @param family the optional glm family (glmRespMod only)
-##' @param nlenv the nonlinear model evaluation environment (nlsRespMod only)
-##' @param nlmod the nonlinear model function (nlsRespMod only)
-##' @return a lmerResp (or glmerResp) instance
+##' @param family the optional glm family (glmResp only)
+##' @param nlenv the nonlinear model evaluation environment (nlsResp only)
+##' @param nlmod the nonlinear model function (nlsResp only)
+##' @return an lmerResp or glmResp or nlsResp instance
mkRespMod2 <- function(fr, family = NULL, nlenv = NULL, nlmod = NULL) {
y <- model.response(fr)
if(length(dim(y)) == 1) {
@@ -251,15 +252,16 @@
}
rho <- new.env()
rho$y <- y
- n <- N <- nrow(fr)
+ n <- nrow(fr)
if (!is.null(nlenv)) {
- stopifnot(is.numeric(val <- eval(nlmod, nlenv)),
- length(val) == N,
+ stopifnot(is.language(nlmod),
+ is.environment(nlenv),
+ is.numeric(val <- eval(nlmod, nlenv)),
+ length(val) == n,
is.matrix(gr <- attr(val, "gradient")),
mode(gr) == "numeric",
- nrow(gr) == N,
+ nrow(gr) == n,
!is.null(pnames <- colnames(gr)))
- N <- length(gr)
}
if (!is.null(offset <- model.offset(fr))) {
if (length(offset) == 1L) offset <- rep.int(offset, N)
@@ -270,29 +272,30 @@
stopifnot(length(weights) == n, all(weights >= 0))
rho$weights <- unname(weights)
}
- if (is.null(family) && is.null(nlenv))
- return(do.call(new, c(list(Class="lmerResp"), as.list(rho))))
- if (!is.null(family)) {
- stopifnot(inherits(family, "family"))
- # need weights for initialize evaluation
- if (!exists("weights", rho, inherits=FALSE))
- rho$weights <- rep.int(1, n)
- rho$nobs <- n
- eval(family$initialize, rho)
- family$initialize <- NULL # remove clutter from str output
- ll <- as.list(rho)
- ans <- do.call(new, c(list(Class="glmResp", family=family),
- ll[setdiff(names(ll),
- c("m", "nobs", "mustart"))]))
- ans$updateMu(if (!is.null(es <- model.extract(fr, "etastart"))) es else
- family$linkfun(get("mustart", rho)))
- return(ans)
+ if (is.null(family)) {
+ if (is.null(nlenv)) return(do.call(lmerResp$new, as.list(rho)))
+ return(do.call(nlsResp$new,
+ c(list(nlenv=nlenv,
+ nlmod=substitute(~foo, list(foo=nlmod)),
+ pnames=pnames, N=length(gr)), as.list(rho))))
}
- stopifnot(is.language(nlmod), is.environment(nlenv))
- do.call(new, c(list(Class="nlsResp", nlenv=nlenv, nlmod=nlmod,
- pnames=pnames), as.list(rho)))
+ stopifnot(inherits(family, "family"))
+ # need weights for initialize evaluation
+ if (!exists("weights", rho, inherits=FALSE))
+ rho$weights <- rep.int(1, n)
+ rho$nobs <- n
+ eval(family$initialize, rho)
+ family$initialize <- NULL # remove clutter from str output
+ ll <- as.list(rho)
+ ans <- do.call(new, c(list(Class="glmResp", family=family),
+ ll[setdiff(names(ll),
+ c("m", "nobs", "mustart"))]))
+ ans$updateMu(if (!is.null(es <- model.extract(fr, "etastart"))) es else
+ family$linkfun(get("mustart", rho)))
+ ans
}
+
###' Create Z, Lambda, Lind, etc.
##'
##' @param bars a list of parsed random-effects terms
@@ -395,7 +398,7 @@
ll$flist <- fl
ll$cnms <- cnms
ll
-} ## {mkReTrms2}
+} ## {mkReTrms}
##' Determine a step factor that will reduce the pwrss
##'
@@ -652,8 +655,8 @@
##' @return an object of S4 class "merMod"
nlmer <- function(formula, data, family = gaussian, start = NULL,
verbose = 0L, nAGQ = 1L, doFit = TRUE,
- subset, weights, na.action, mustart, etastart,
- sparseX = FALSE, contrasts = NULL, control = list(), ...)
+ subset, weights, na.action, offset,
+ contrasts = NULL, devFunOnly=FALSE, ...)
{
if (!missing(family)) stop("code not yet written")
mf <- mc <- match.call()
@@ -690,6 +693,7 @@
nlenv <- new.env() # inherit from this environment (or environment(formula)?)
lapply(all.vars(nlmod),
function(nm) assign(nm, fr[[nm]], envir = nlenv))
+ rr <- mkRespMod2(fr, nlenv=nlenv, nlmod=nlmod)
## Second, extend the frame and convert the nlpar columns to indicators
n <- nrow(fr)
@@ -698,24 +702,19 @@
frE[[nm]] <- as.numeric(rep(nm == pnames, each = n))
# random-effects module
reTrms <- mkReTrms(findbars(formula[[3]]), frE, s = s)
- dcmp <- updateDcmp(reTrms, .dcmp())
- dcmp$dims["nAGQ"] <- as.integer(nAGQ)[1]
fe.form <- nlform
fe.form[[3]] <- formula[[3]]
- feMod <- mkFeModule(fe.form, frE, contrasts, reTrms, sparseX = sparseX)
- # should this check be in mkFeModule?
- p <- length(feMod at coef)
- if ((qrX <- qr(feMod at X))$rank < p)
+ X <- model.Matrix(nobars(fe.form), frE, contrasts,
+ sparse = FALSE, row.names = FALSE) ## sparseX not yet
+ p <- ncol(X)
+ if ((qrX <- qr(X))$rank < p)
stop(gettextf("rank of X = %d < ncol(X) = %d", qrX$rank, p))
- feMod at coef <- qr.coef(qrX, unlist(lapply(pnames, get, envir = nlenv)))
- respMod <- mkRespMod2(fr, nlenv = nlenv, nlmod = nlmod)
- ans <- new("merMod",
- call = mc,
- devcomp = updateDcmp(respMod, updateDcmp(feMod, dcmp)),
- frame = fr, re = reTrms, fe = feMod, resp = respMod)
- if (!doFit) return(ans)
- PIRLSest(ans, verbose, control, nAGQ)
+ pp <- do.call(merPredD$new, c(reTrms[c("Zt","theta","Lambdat","Lind")],
+ list(X=X,
+ beta0=qr.coef(qrX, unlist(lapply(pnames, get, envir = nlenv))))
+ ))
+ return(list(pp=pp, resp=rr))
}
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/src/external.cpp 2011-10-14 19:54:51 UTC (rev 1416)
@@ -8,9 +8,14 @@
#include "respModule.h"
extern "C" {
+ using Eigen::Map;
+ using Eigen::MatrixXd;
using Eigen::VectorXd;
+ using Rcpp::CharacterVector;
+ using Rcpp::Environment;
using Rcpp::IntegerVector;
+ using Rcpp::Language;
using Rcpp::List;
using Rcpp::NumericVector;
using Rcpp::S4;
@@ -24,6 +29,7 @@
using lme4Eigen::lmResp;
using lme4Eigen::lmerResp;
using lme4Eigen::merPredD;
+ using lme4Eigen::nlsResp;
using std::runtime_error;
@@ -130,7 +136,7 @@
SEXP glm_updateMu(SEXP ptr_, SEXP gamma) {
BEGIN_RCPP;
- return ::Rf_ScalarReal(XPtr<glmResp>(ptr_)->updateMu(as<VectorXd>(gamma)));
+ return ::Rf_ScalarReal(XPtr<glmResp>(ptr_)->updateMu(as<Map<VectorXd> >(gamma)));
END_RCPP;
}
@@ -145,27 +151,27 @@
SEXP glmFamily_link(SEXP ptr, SEXP mu) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->linkFun(as<VectorXd>(mu)));
+ return wrap(XPtr<glmFamily>(ptr)->linkFun(as<Map<VectorXd> >(mu)));
END_RCPP;
}
SEXP glmFamily_linkInv(SEXP ptr, SEXP eta) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->linkInv(as<VectorXd>(eta)));
+ return wrap(XPtr<glmFamily>(ptr)->linkInv(as<Map<VectorXd> >(eta)));
END_RCPP;
}
SEXP glmFamily_devResid(SEXP ptr, SEXP mu, SEXP weights, SEXP y) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->devResid(as<VectorXd>(mu),
- as<VectorXd>(weights),
- as<VectorXd>(y)));
+ return wrap(XPtr<glmFamily>(ptr)->devResid(as<Map<VectorXd> >(mu),
+ as<Map<VectorXd> >(weights),
+ as<Map<VectorXd> >(y)));
END_RCPP;
}
SEXP glmFamily_muEta(SEXP ptr, SEXP eta) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->muEta(as<VectorXd>(eta)));
+ return wrap(XPtr<glmFamily>(ptr)->muEta(as<Map<VectorXd> >(eta)));
END_RCPP;
}
@@ -268,12 +274,14 @@
SEXP lm_setOffset(SEXP ptr_, SEXP offset) {
BEGIN_RCPP;
XPtr<lmResp>(ptr_)->setOffset(as<VectorXd>(offset));
+ return offset;
END_RCPP;
}
SEXP lm_setWeights(SEXP ptr_, SEXP weights) {
BEGIN_RCPP;
XPtr<lmResp>(ptr_)->setWeights(as<VectorXd>(weights));
+ return weights;
END_RCPP;
}
@@ -327,7 +335,7 @@
SEXP lm_updateMu(SEXP ptr_, SEXP gamma) {
BEGIN_RCPP;
- return ::Rf_ScalarReal(XPtr<lmerResp>(ptr_)->updateMu(as<VectorXd>(gamma)));
+ return ::Rf_ScalarReal(XPtr<lmerResp>(ptr_)->updateMu(as<Map<VectorXd> >(gamma)));
END_RCPP;
}
@@ -342,7 +350,9 @@
SEXP lmer_setREML(SEXP ptr_, SEXP REML) {
BEGIN_RCPP;
- XPtr<lmerResp>(ptr_)->setReml(::Rf_asInteger(REML));
+ int reml = ::Rf_asInteger(REML);
+ XPtr<lmerResp>(ptr_)->setReml(reml);
+ return ::Rf_ScalarInteger(reml);
END_RCPP;
}
@@ -392,7 +402,7 @@
// setters
SEXP merPredDsetBeta0(SEXP ptr, SEXP beta0) {
BEGIN_RCPP;
- XPtr<merPredD>(ptr)->setBeta0(as<Eigen::VectorXd>(beta0));
+ XPtr<merPredD>(ptr)->setBeta0(as<Map<VectorXd> >(beta0));
END_RCPP;
}
@@ -404,7 +414,7 @@
SEXP merPredDsetU0(SEXP ptr, SEXP u0) {
BEGIN_RCPP;
- XPtr<merPredD>(ptr)->setU0(as<Eigen::VectorXd>(u0));
+ XPtr<merPredD>(ptr)->setU0(as<Map<VectorXd> >(u0));
END_RCPP;
}
@@ -606,16 +616,42 @@
SEXP merPredDupdateRes(SEXP ptr, SEXP wtres) {
BEGIN_RCPP;
- XPtr<merPredD>(ptr)->updateRes(as<VectorXd>(wtres));
+ XPtr<merPredD>(ptr)->updateRes(as<Map<VectorXd> >(wtres));
END_RCPP;
}
SEXP merPredDupdateXwts(SEXP ptr, SEXP wts) {
BEGIN_RCPP;
- XPtr<merPredD>(ptr)->updateXwts(as<VectorXd>(wts));
+ XPtr<merPredD>(ptr)->updateXwts(as<Map<MatrixXd> >(wts));
END_RCPP;
}
+ SEXP nls_Create(SEXP ys, SEXP mods, SEXP envs, SEXP pnms, SEXP Ns) {
+ BEGIN_RCPP;
+ nlsResp *ans =
+ new nlsResp(NumericVector(ys), Language(mods),
+ Environment(envs), CharacterVector(pnms),
+ ::Rf_asInteger(Ns));
+ return wrap(XPtr<nlsResp>(ans, true));
+ END_RCPP;
+ }
+
+ SEXP nls_Laplace(SEXP ptr_, SEXP ldL2, SEXP ldRX2, SEXP sqrL) {
+ BEGIN_RCPP;
+ return ::Rf_ScalarReal(XPtr<nlsResp>(ptr_)->
+ Laplace(::Rf_asReal(ldL2),
+ ::Rf_asReal(ldRX2),
+ ::Rf_asReal(sqrL)));
+ END_RCPP;
+ }
+
+ SEXP nls_updateMu(SEXP ptr_, SEXP gamma) {
+ BEGIN_RCPP;
+ return ::Rf_ScalarReal(XPtr<nlsResp>(ptr_)->updateMu(as<Map<VectorXd> >(gamma)));
+ END_RCPP;
+ }
+
+
}
#include <R_ext/Rdynload.h>
@@ -725,6 +761,11 @@
CALLDEF(merPredDupdateDecomp, 1),
CALLDEF(merPredDupdateRes, 2),
CALLDEF(merPredDupdateXwts, 2),
+
+ CALLDEF(nls_Create, 5), // generate external pointer
+
+ CALLDEF(nls_Laplace, 4), // methods
+ CALLDEF(nls_updateMu, 2),
{NULL, NULL, 0}
};
Modified: pkg/lme4Eigen/src/predModule.cpp
===================================================================
--- pkg/lme4Eigen/src/predModule.cpp 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/src/predModule.cpp 2011-10-14 19:54:51 UTC (rev 1416)
@@ -49,9 +49,9 @@
d_theta(clone(theta)),
d_Lind(Lind),
d_n(d_X.rows()),
+ d_nnz(-1),
d_p(d_X.cols()),
d_q(d_Zt.rows()),
- d_nnz(-1),
d_Lambdat(Lambdat),
d_RZX(d_q, d_p),
d_V(d_X),
@@ -63,7 +63,8 @@
d_beta0(d_p),
d_u0(d_q),
d_Ut(d_Zt),
- d_RX(d_p)
+ d_RX(d_p),
+ d_LamtUtRestructure(false)
{ // Check consistency of dimensions
if (d_n != d_Zt.cols())
throw invalid_argument("Z dimension mismatch");
@@ -76,7 +77,6 @@
setTheta(d_theta); // starting values into Lambda
d_L.cholmod().final_ll = 1; // force an LL' decomposition
-// FIXME: updating d_Ut should be done in a method to accomodate the nlmer case
d_LamtUt = d_Lambdat * d_Ut;
d_L.analyzePattern(d_LamtUt); // perform symbolic analysis
if (d_L.info() != Eigen::Success)
@@ -84,16 +84,23 @@
}
void merPredD::updateLamtUt() {
+ if (d_LamtUtRestructure) {
+ d_LamtUt = d_Lambdat * d_Ut;
+ return;
+ }
+ // This complicated code bypasses problems caused by Eigen's
+ // sparse/sparse matrix multiplication pruning zeros. The
+ // Cholesky decomposition croaks if the structure of d_LamtUt changes.
fill(d_LamtUt._valuePtr(), d_LamtUt._valuePtr() + d_LamtUt.nonZeros(), Scalar());
for (Index j = 0; j < d_Ut.cols(); ++j) {
- for (SpMatrixd::InnerIterator rhsIt(d_Ut, j); rhsIt; ++rhsIt) {
+ for(SpMatrixd::InnerIterator rhsIt(d_Ut, j); rhsIt; ++rhsIt) {
SpMatrixd::InnerIterator prdIt(d_LamtUt, j);
Scalar y = rhsIt.value();
for (SpMatrixd::InnerIterator lhsIt(d_Lambdat, rhsIt.index()); lhsIt; ++lhsIt) {
- Index i = lhsIt.index();
+ Index i = lhsIt.index();
while (prdIt.index() != i && prdIt) ++prdIt;
if (!prdIt) throw runtime_error("logic error in updateLamtUt");
- prdIt.valueRef() += lhsIt.value() * y;
+ prdIt.valueRef() += lhsIt.value() * y;
}
}
}
@@ -157,14 +164,28 @@
return d_CcNumer;
}
- void merPredD::updateXwts(const VectorXd& sqrtXwt) {
+ void merPredD::updateXwts(const MatrixXd& sqrtXwt) {
if (d_X.rows() != sqrtXwt.size())
throw invalid_argument("updateXwts: dimension mismatch");
- if (d_V.rows() == d_X.rows()) { //FIXME: Generalize this for nlmer
- DiagonalMatrix<double, Dynamic> W(sqrtXwt.asDiagonal());
- d_V = W * d_X;
- d_Ut = d_Zt * W;
- } else throw invalid_argument("updateXwts: no provision for nlmer yet");
+ if (sqrtXwt.cols() == 1) {
+ DiagonalMatrix<double, Dynamic> W(sqrtXwt.col(1).asDiagonal());
+ d_V = W * d_X;
+ d_Ut = d_Zt * W;
+ } else {
+ int n = sqrtXwt.rows();
+ SpMatrixd W(n, sqrtXwt.size());
+ const double *pt = sqrtXwt.data();
+ W.reserve(sqrtXwt.size());
+ for (Index j = 0; j < W.cols(); ++j, ++pt) {
+ W.startVec(j);
+ W.insertBack(j % n, j) = *pt;
+ }
+ W.finalize();
+ d_V = W * d_X;
+ d_Ut = d_Zt * W.adjoint();
+ }
+ if (d_LamtUt.rows() != d_Ut.rows() ||
+ d_LamtUt.cols() != d_Ut.cols()) d_LamtUtRestructure = true;
d_VtV.setZero().selfadjointView<Upper>().rankUpdate(d_V.adjoint());
updateL();
}
Modified: pkg/lme4Eigen/src/predModule.h
===================================================================
--- pkg/lme4Eigen/src/predModule.h 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/src/predModule.h 2011-10-14 19:54:51 UTC (rev 1416)
@@ -158,7 +158,7 @@
MdgCMatrix d_Zt;
NumericVector d_theta;
IntegerVector d_Lind;
- Index d_n, d_p, d_q, d_nnz;
+ Index d_n, d_nnz, d_p, d_q;
dgCMatrix d_Lambdat;
Scalar d_CcNumer, d_ldL2, d_ldRX2;
MatrixXd d_RZX, d_V, d_VtV;
@@ -166,6 +166,7 @@
SpMatrixd d_Ut, d_LamtUt;
ChmDecomp d_L;
LLT<MatrixXd> d_RX;
+ bool d_LamtUtRestructure;
public:
merPredD(S4, S4, S4, IntegerVector, NumericVector);
@@ -213,13 +214,14 @@
void installPars(const Scalar& f);
void setBeta0(const VectorXd&);
+ void setS(int);
void setTheta(const NumericVector&);
void setU0(const VectorXd&);
void updateDecomp();
void updateL();
void updateLamtUt();
void updateRes(const VectorXd&);
- void updateXwts(const VectorXd&);
+ void updateXwts(const MatrixXd&);
};
}
Modified: pkg/lme4Eigen/src/respModule.cpp
===================================================================
--- pkg/lme4Eigen/src/respModule.cpp 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/src/respModule.cpp 2011-10-14 19:54:51 UTC (rev 1416)
@@ -22,13 +22,13 @@
d_weights(VectorXd::Constant(y.size(), 1.0)),
d_offset( VectorXd::Zero(y.size())),
d_mu( y.size()),
- d_sqrtXwt(y.size()),
d_sqrtrwt(y.size()),
- d_wtres( y.size()) {
+ d_wtres( y.size()),
+ d_sqrtXwt(y.size(), 1) {
d_sqrtrwt = d_weights.cwiseSqrt();
if (d_mu.size() == d_offset.size()) {
d_mu = d_offset;
- copy(d_sqrtrwt.begin(), d_sqrtrwt.end(), d_sqrtXwt.begin());
+ copy(d_sqrtrwt.data(), d_sqrtrwt.data() + y.size(), d_sqrtXwt.data());
}
updateWrss();
}
@@ -140,19 +140,23 @@
d_n = n;
}
- nlmerResp::nlmerResp(NumericVector ys, Language mm, Environment ee, CharacterVector pp)
- : lmResp(ys), d_nlenv(ee), d_nlmod(mm), d_pnames(pp) {
+ nlsResp::nlsResp(NumericVector ys, Language mm, Environment ee, CharacterVector pp, int N)
+ : lmResp(ys), d_nlenv(ee), d_nlmod(mm), d_pnames(pp), d_N(N) {
+ if (d_N <= 0 || d_N % d_y.size())
+ throw invalid_argument("N must be a positive multiple of y.size()");
+ d_offset = VectorXd::Zero(N);
+ d_sqrtXwt = MatrixXd::Zero(d_y.size(), N / d_y.size());
}
- double nlmerResp::Laplace(double ldL2, double ldRX2, double sqrL) const {
+ double nlsResp::Laplace(double ldL2, double ldRX2, double sqrL) const {
double lnum = 2.* PI * (d_wrss + sqrL), n = d_y.size();
return ldL2 + n * (1 + log(lnum / n));
}
- double nlmerResp::updateMu(VectorXd const &gamma) {
+ double nlsResp::updateMu(VectorXd const &gamma) {
int n = d_y.size();
VectorXd gam = gamma + d_offset;
- const double *gg = gam.begin();
+ const double *gg = gam.data();
for (int p = 0; p < d_pnames.size(); p++) {
string pn(d_pnames[p]);
@@ -162,9 +166,10 @@
NumericVector rr = d_nlmod.eval(SEXP(d_nlenv));
if (rr.size() != n)
throw invalid_argument("dimension mismatch");
- copy(rr.begin(), rr.end(), d_mu.begin());
- NumericMatrix rrg = rr.attr("gradient");
- copy(rrg.begin(), rrg.end(), d_sqrtXwt.begin());
+ copy(rr.begin(), rr.end(), d_mu.data());
+ NumericMatrix gr = rr.attr("gradient");
+ d_sqrtXwt.resize(gr.rows(), gr.cols());
+ copy(gr.begin(), gr.end(), d_sqrtXwt.data());
return updateWrss();
}
Modified: pkg/lme4Eigen/src/respModule.h
===================================================================
--- pkg/lme4Eigen/src/respModule.h 2011-10-03 16:06:28 UTC (rev 1415)
+++ pkg/lme4Eigen/src/respModule.h 2011-10-14 19:54:51 UTC (rev 1416)
@@ -14,7 +14,9 @@
namespace lme4Eigen {
using Eigen::Map;
using Eigen::VectorXd;
+ using Eigen::MatrixXd;
+ using Rcpp::as;
using Rcpp::CharacterVector;
using Rcpp::Environment;
using Rcpp::Language;
@@ -27,13 +29,14 @@
double d_wrss;
const NumericVector d_yR;
const Map<VectorXd> d_y;
- VectorXd d_weights, d_offset, d_mu, d_sqrtXwt, d_sqrtrwt, d_wtres;
+ VectorXd d_weights, d_offset, d_mu, d_sqrtrwt, d_wtres;
+ MatrixXd d_sqrtXwt;
public:
lmResp(NumericVector);
+ const MatrixXd& sqrtXwt() const {return d_sqrtXwt;}
const VectorXd& mu() const {return d_mu;}
const VectorXd& offset() const {return d_offset;}
- const VectorXd& sqrtXwt() const {return d_sqrtXwt;}
const VectorXd& sqrtrwt() const {return d_sqrtrwt;}
const VectorXd& weights() const {return d_weights;}
const VectorXd& wtres() const {return d_wtres;}
@@ -85,13 +88,14 @@
void setN(const VectorXd&);
};
- class nlmerResp : public lmResp {
+ class nlsResp : public lmResp {
protected:
Environment d_nlenv;
Language d_nlmod;
CharacterVector d_pnames;
+ int d_N;
public:
- nlmerResp(NumericVector,Language,Environment,CharacterVector);
+ nlsResp(NumericVector,Language,Environment,CharacterVector,int);
double Laplace(double, double, double) const;
double updateMu(const VectorXd&);
More information about the Lme4-commits
mailing list