[Lme4-commits] r1495 - in pkg/lme4Eigen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jan 6 16:00:07 CET 2012
Author: dmbates
Date: 2012-01-06 16:00:07 +0100 (Fri, 06 Jan 2012)
New Revision: 1495
Modified:
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/predModule.cpp
pkg/lme4Eigen/src/respModule.cpp
pkg/lme4Eigen/src/respModule.h
Log:
Modifications for nlmer - not yet finished.
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2011-12-31 20:38:15 UTC (rev 1494)
+++ pkg/lme4Eigen/R/AllClass.R 2012-01-06 15:00:07 UTC (rev 1495)
@@ -57,7 +57,6 @@
V <<- array(0, c(n, p))
VtV <<- array(0, c(p, p))
Vtr <<- numeric(p)
- Xwts <<- rep.int(1, N)
beta0 <<- numeric(p)
delb <<- numeric(p)
delu <<- numeric(q)
@@ -65,6 +64,8 @@
Ut <<- if (S == 1) Zt + 0 else
Zt %*% sparseMatrix(i=seq_len(N), j=as.integer(gl(n, 1, N)), x=rep.int(1,N))
LamtUt <<- Lambdat %*% Ut
+ Xw <- list(...)$Xwts
+ Xwts <<- if (is.null(Xw)) rep.int(1, N) else as.numeric(Xw)
updateXwts(Xwts)
},
CcNumer = function() {
@@ -202,11 +203,16 @@
if (is.null(ll$y)) stop("y must be specified")
y <<- as.numeric(ll$y)
n <- length(y)
- mu <<- if (!is.null(ll$mu)) as.numeric(ll$mu) else numeric(n)
- offset <<- if (!is.null(ll$offset)) as.numeric(ll$offset) else numeric(n)
- weights <<- if (!is.null(ll$weights)) as.numeric(ll$weights) else rep.int(1,n)
- sqrtXwt <<- sqrt(weights)
- sqrtrwt <<- sqrt(weights)
+ mu <<- if (!is.null(ll$mu))
+ as.numeric(ll$mu) else numeric(n)
+ offset <<- if (!is.null(ll$offset))
+ as.numeric(ll$offset) else numeric(n)
+ weights <<- if (!is.null(ll$weights))
+ as.numeric(ll$weights) else rep.int(1,n)
+ sqrtXwt <<- if (!is.null(ll$sqrtXwt))
+ as.numeric(ll$sqrtXwt) else sqrt(weights)
+ sqrtrwt <<- if (!is.null(ll$sqrtrwt))
+ as.numeric(ll$sqrtrwt) else sqrt(weights)
wtres <<- sqrtrwt * (y - mu)
},
ptr = function() {
@@ -356,24 +362,41 @@
nlsResp <-
setRefClass("nlsResp",
fields=
- list(N="integer",
+ list(gam="numeric",
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, y, nlmod[[2]], nlenv, pnames, N),
- silent=TRUE)
- Ptr
- }),
+ pnames="character"
+ ),
contains="lmResp",
methods=
- list(
- Laplace=function(ldL2, ldRX2, sqrL) {
+ list(initialize = function(...) {
+ callSuper(...)
+ ll <- list(...)
+ if (is.null(ll$nlmod)) stop("nlmod must be specified")
+ nlmod <<- ll$nlmod
+ if (is.null(ll$nlenv)) stop("nlenv must be specified")
+ nlenv <<- ll$nlenv
+ if (is.null(ll$pnames)) stop("pnames must be specified")
+ pnames <<- ll$pnames
+ if (is.null(ll$gam)) stop("gam must be specified")
+ stopifnot(length(ll$gam) == length(offset))
+ gam <<- ll$gam
+ },
+ Laplace =function(ldL2, ldRX2, sqrL) {
'returns the profiled deviance or REML criterion'
.Call(nls_Laplace, ptr(), ldL2, ldRX2, sqrL)
},
+ ptr = function() {
+ 'returns the external pointer, regenerating if necessary'
+ if (length(y)) {
+ if (.Call(isNullExtPtr, Ptr)) {
+ Ptr <<- .Call(nls_Create, y, weights, offset, mu, sqrtXwt,
+ sqrtrwt, wtres, gam, nlmod[[2]], nlenv, pnames)
+ .Call(nls_updateMu, Ptr, gam)
+ }
+ }
+ Ptr
+ },
updateMu=function(gamma) {
'update mu, residuals, gradient, etc. given the linear predictor matrix'
.Call(nls_updateMu, ptr(), as.numeric(gamma))
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2011-12-31 20:38:15 UTC (rev 1494)
+++ pkg/lme4Eigen/R/lmer.R 2012-01-06 15:00:07 UTC (rev 1495)
@@ -273,27 +273,30 @@
nrow(gr) == n,
!is.null(pnames <- colnames(gr)))
N <- length(gr)
+ rho$mu <- as.vector(val)
+ rho$sqrtXwt <- as.vector(gr)
+ rho$gam <-
+ unname(unlist(lapply(pnames,
+ function(nm) get(nm, envir=nlenv))))
}
if (!is.null(offset <- model.offset(fr))) {
if (length(offset) == 1L) offset <- rep.int(offset, N)
stopifnot(length(offset) == N)
rho$offset <- unname(offset)
- }
+ } else rho$offset <- rep.int(0, N)
if (!is.null(weights <- model.weights(fr))) {
stopifnot(length(weights) == n, all(weights >= 0))
rho$weights <- unname(weights)
- }
+ } else rho$weights <- rep.int(1, n)
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))))
+ 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)
+ # need weights for initialize evaluation
rho$nobs <- n
eval(family$initialize, rho)
family$initialize <- NULL # remove clutter from str output
@@ -332,7 +335,7 @@
nseq <- seq_len(nc)
sm <- as(ff, "sparseMatrix")
if (nc > 1)
- sm <- do.call(rBind, lapply(nseq, function(i) sm))
+ sm <- do.call(Matrix::rBind, lapply(nseq, function(i) sm))
sm at x[] <- t(mm[])
## When nc > 1 switch the order of the rows of sm
## so the random effects for the same level of the
@@ -351,7 +354,7 @@
blist <- blist[ord]
nl <- nl[ord]
}
- Zt <- do.call(rBind, lapply(blist, "[[", "sm"))
+ Zt <- do.call(Matrix::rBind, lapply(blist, "[[", "sm"))
q <- nrow(Zt)
## Create and install Lambdat, Lind, etc. This must be done after
@@ -368,7 +371,7 @@
### instead of generating Lambda first then transposing?
Lambdat <-
t(do.call(sparseMatrix,
- do.call(rBind,
+ do.call(Matrix::rBind,
lapply(seq_along(blist), function(i)
{
mm <- matrix(seq_len(nb[i]), nc = nc[i],
@@ -703,7 +706,8 @@
lapply(all.vars(nlmod),
function(nm) assign(nm, fr[[nm]], envir = nlenv))
rr <- mkRespMod2(fr, nlenv=nlenv, nlmod=nlmod)
-
+ stopifnot(all(pnames %in% rr$pnames), all(rr$pnames %in% pnames))
+
## Second, extend the frame and convert the nlpar columns to indicators
n <- nrow(fr)
frE <- do.call(rbind, lapply(1:s, function(i) fr)) # rbind s copies of the frame
@@ -714,13 +718,14 @@
fe.form <- nlform
fe.form[[3]] <- formula[[3]]
- X <- model.Matrix(nobars(fe.form), frE, contrasts,
- sparse = FALSE, row.names = FALSE) ## sparseX not yet
+ X <- model.matrix(nobars(fe.form), frE, contrasts)
+## sparse = FALSE, row.names = FALSE) ## sparseX not yet
+ rownames(X) <- NULL
p <- ncol(X)
if ((qrX <- qr(X))$rank < p)
stop(gettextf("rank of X = %d < ncol(X) = %d", qrX$rank, p))
pp <- do.call(merPredD$new, c(reTrms[c("Zt","theta","Lambdat","Lind")],
- list(X=X,
+ list(X=X, S=s, Xwts = rr$sqrtXwt,
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-12-31 20:38:15 UTC (rev 1494)
+++ pkg/lme4Eigen/src/external.cpp 2012-01-06 15:00:07 UTC (rev 1495)
@@ -676,12 +676,10 @@
// nonlinear model response (also the base class for other response classes)
SEXP nls_Create(SEXP y, SEXP weights, SEXP offset, SEXP mu, SEXP sqrtXwt,
- SEXP sqrtrwt, SEXP wtres, SEXP mods, SEXP envs, SEXP pnms) {
+ SEXP sqrtrwt, SEXP wtres, SEXP gamma, SEXP mod, SEXP env, SEXP pnms) {
BEGIN_RCPP;
nlsResp *ans =
- new nlsResp(y, weights, offset, mu, sqrtXwt, sqrtrwt, wtres,
- Language(mods), Environment(envs),
- CharacterVector(pnms));
+ new nlsResp(y, weights, offset, mu, sqrtXwt, sqrtrwt, wtres, gamma, mod, env, pnms);
return wrap(XPtr<nlsResp>(ans, true));
END_RCPP;
}
@@ -805,7 +803,7 @@
CALLDEF(NelderMead_xeval, 1),
CALLDEF(NelderMead_xpos, 1),
- CALLDEF(nls_Create, 10), // generate external pointer
+ CALLDEF(nls_Create, 11), // generate external pointer
CALLDEF(nls_Laplace, 4), // methods
CALLDEF(nls_updateMu, 2),
Modified: pkg/lme4Eigen/src/predModule.cpp
===================================================================
--- pkg/lme4Eigen/src/predModule.cpp 2011-12-31 20:38:15 UTC (rev 1494)
+++ pkg/lme4Eigen/src/predModule.cpp 2012-01-06 15:00:07 UTC (rev 1495)
@@ -170,8 +170,11 @@
}
W.finalize();
d_V = W * d_X;
-//FIXME: work out the corresponding code for Ut and Zt
-// d_Ut = d_Zt * W.adjoint();
+// FIXME: Need to work out the equivalent code here for d_Ut = d_Zt * W.adjoint()
+ // for (int j = 0; j < d_N; ++j)
+ // for (MSpMatrixd::InnerIterator Uit(d_Ut, j), Zit(d_Zt, j);
+ // Uit && Zit; ++Uit, ++Zit)
+ // Uit.valueRef() = Zit.value() * sqrtXwt.data()[j];
}
// if (d_LamtUt.rows() != d_Ut.rows() ||
// d_LamtUt.cols() != d_Ut.cols()) d_LamtUtRestructure = true;
Modified: pkg/lme4Eigen/src/respModule.cpp
===================================================================
--- pkg/lme4Eigen/src/respModule.cpp 2011-12-31 20:38:15 UTC (rev 1494)
+++ pkg/lme4Eigen/src/respModule.cpp 2012-01-06 15:00:07 UTC (rev 1495)
@@ -1,6 +1,6 @@
// respModule.cpp: response modules using Eigen
//
-// Copyright (C) 2011 Douglas Bates, Martin Maechler and Ben Bolker
+// Copyright (C) 2011-2012 Douglas Bates, Martin Maechler and Ben Bolker
//
// This file is part of lme4.
@@ -142,10 +142,13 @@
}
nlsResp::nlsResp(SEXP y, SEXP weights, SEXP offset, SEXP mu, SEXP sqrtXwt,
- SEXP sqrtrwt, SEXP wtres, Language mm, Environment ee,
- CharacterVector pp)
+ SEXP sqrtrwt, SEXP wtres, SEXP gamma, SEXP mm, SEXP ee,
+ SEXP pp)
: lmResp(y, weights, offset, mu, sqrtXwt, sqrtrwt, wtres),
- d_nlenv(ee), d_nlmod(mm), d_pnames(pp) {
+ d_gamma(as<MVec>(gamma)),
+ d_nlenv(as<Environment>(ee)),
+ d_nlmod(as<Language>(mm)),
+ d_pnames(as<CharacterVector>(pp)) {
}
double nlsResp::Laplace(double ldL2, double ldRX2, double sqrL) const {
@@ -153,12 +156,15 @@
return ldL2 + n * (1 + log(lnum / n));
}
- double nlsResp::updateMu(VectorXd const &gamma) {
+ double nlsResp::updateMu(const VectorXd& gamma) {
int n = d_y.size();
- VectorXd gam = gamma + d_offset;
- const double *gg = gam.data();
+ if (gamma.size() != d_gamma.size())
+ throw invalid_argument("size mismatch in updateMu");
+ std::copy(gamma.data(), gamma.data() + gamma.size(), d_gamma.data());
+ const VectorXd lp(d_gamma + d_offset); // linear predictor
+ const double *gg = lp.data();
- for (int p = 0; p < d_pnames.size(); p++) {
+ for (int p = 0; p < d_pnames.size(); ++p) {
string pn(d_pnames[p]);
NumericVector pp = d_nlenv.get(pn);
copy(gg + n * p, gg + n * (p + 1), pp.begin());
Modified: pkg/lme4Eigen/src/respModule.h
===================================================================
--- pkg/lme4Eigen/src/respModule.h 2011-12-31 20:38:15 UTC (rev 1494)
+++ pkg/lme4Eigen/src/respModule.h 2012-01-06 15:00:07 UTC (rev 1495)
@@ -2,7 +2,7 @@
//
// respModule.h: response modules using Eigen
//
-// Copyright (C) 2011 Douglas Bates, Martin Maechler and Ben Bolker
+// Copyright (C) 2011-2012 Douglas Bates, Martin Maechler and Ben Bolker
//
// This file is part of lme4.
@@ -85,11 +85,12 @@
class nlsResp : public lmResp {
protected:
+ MVec d_gamma;
Environment d_nlenv;
Language d_nlmod;
CharacterVector d_pnames;
public:
- nlsResp(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,Language,Environment,CharacterVector);
+ nlsResp(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
double Laplace(double, double, double) const;
double updateMu(const Eigen::VectorXd&);
More information about the Lme4-commits
mailing list