[Lme4-commits] r1472 - in pkg/lme4Eigen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 7 00:25:08 CET 2011
Author: dmbates
Date: 2011-12-07 00:25:08 +0100 (Wed, 07 Dec 2011)
New Revision: 1472
Modified:
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
Log:
Add many more fields to the merPredD class and allow for serialize/unserialize. Lower the default tolPwrss as per testing by Ben.
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2011-12-06 17:42:13 UTC (rev 1471)
+++ pkg/lme4Eigen/R/AllClass.R 2011-12-06 23:25:08 UTC (rev 1472)
@@ -15,27 +15,57 @@
merPredD <-
setRefClass("merPredD", # Predictor class for mixed-effects models with dense X
fields =
- list(Ptr = "externalptr",
+ list(Lambdat = "dgCMatrix",
+ LamtUt = "dgCMatrix",
+ Lind = "integer",
+ Ptr = "externalptr",
+ RZX = "matrix",
+ Ut = "dgCMatrix",
+ Utr = "numeric",
+ V = "matrix",
+ VtV = "matrix",
+ Vtr = "numeric",
X = "ddenseModelMatrix",
Zt = "dgCMatrix",
- Lambdat = "dgCMatrix",
- Lind = "integer",
beta0 = "numeric",
- u0 = "numeric",
- theta = "numeric"),
+ delb = "numeric",
+ delu = "numeric",
+ theta = "numeric",
+ u0 = "numeric"),
methods =
list(
- initialize = function(X, Zt, Lambdat, Lind, theta, ...) {
+ initialize = function(X, Zt, Lambdat, Lind, theta, S, ...) {
if (!nargs()) return
X <<- as(X, "ddenseModelMatrix")
Zt <<- as(Zt, "dgCMatrix")
Lambdat <<- as(Lambdat, "dgCMatrix")
Lind <<- as.integer(Lind)
- theta <<- as.numeric(theta) # force a copy
- stopifnot(length(theta) > 0L, length(Lind) > 0L,
- all(sort(unique(Lind)) == seq_along(theta)))
- if (!length(u0)) u0 <<- numeric(nrow(Zt))
- if (!length(beta0)) beta0 <<- numeric(ncol(X))
+ theta <<- as.numeric(theta)
+ N <- nrow(X)
+ S <- as.integer(S)[1]
+ n <- N %/% S
+ p <- ncol(X)
+ q <- nrow(Zt)
+ stopifnot(length(theta) > 0L,
+ length(Lind) > 0L,
+ all(sort(unique(Lind)) == seq_along(theta)),
+ S > 0L,
+ n * S == N)
+ if (!length(RZX)) RZX <<- array(0, c(q, p))
+ if (!length(Utr)) Utr <<- numeric(q)
+ if (!length(V)) V <<- array(0, c(n, p))
+ if (!length(VtV)) VtV <<- array(0, c(p, p))
+ if (!length(Vtr)) Vtr <<- numeric(p)
+ if (!length(beta0)) beta0 <<- numeric(p)
+ if (!length(delb)) delb <<- numeric(p)
+ if (!length(delu)) delu <<- numeric(q)
+ if (!length(u0)) u0 <<- numeric(q)
+ if (!length(Ut at x)) {
+ 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))
+ }
+ if (!length(LamtUt at x)) LamtUt <<- Lambdat %*% Ut
+ updateXwts(rep.int(1, N))
},
CcNumer = function() {
'returns the numerator of the orthogonality convergence criterion'
@@ -45,10 +75,10 @@
'returns the current value of the sparse Cholesky factor'
.Call(merPredDL, ptr())
},
- LamtUt = function() {
- 'returns the current value of the product Lambdat %*% Ut'
- .Call(merPredDLamtUt, ptr())
- },
+ ## LamtUt = function() {
+ ## 'returns the current value of the product Lambdat %*% Ut'
+ ## .Call(merPredDLamtUt, ptr())
+ ## },
P = function() {
'returns the permutation vector for the sparse Cholesky factor'
.Call(merPredDPvec, ptr())
@@ -61,30 +91,30 @@
'returns the diagonal of the dense downdated Cholesky factor'
.Call(merPredDRXdiag, ptr())
},
- RZX = function() {
- 'returns the cross term in Cholesky decomposition for all coefficients'
- .Call(merPredDRZX, ptr())
- },
- Ut = function() {
- 'returns the transposed weighted random-effects model matrix'
- .Call(merPredDUt, ptr())
- },
- Utr = function() {
- 'returns the cross-product of the weighted random-effects model matrix\nand the weighted residuals'
- .Call(merPredUtr, ptr())
- },
- V = function() {
- 'returns the weighted fixed-effects model matrix'
- .Call(merPredDV, ptr())
- },
- VtV = function() {
- 'returns the weighted cross-product of the fixed-effects model matrix'
- .Call(merPredDVtV, ptr())
- },
- Vtr = function() {
- 'returns the weighted cross-product of the fixed-effects model matrix\nand the residuals'
- .Call(merPredVtr, ptr())
- },
+ ## RZX = function() {
+ ## 'returns the cross term in Cholesky decomposition for all coefficients'
+ ## .Call(merPredDRZX, ptr())
+ ## },
+ ## Ut = function() {
+ ## 'returns the transposed weighted random-effects model matrix'
+ ## .Call(merPredDUt, ptr())
+ ## },
+ ## Utr = function() {
+ ## 'returns the cross-product of the weighted random-effects model matrix\nand the weighted residuals'
+ ## .Call(merPredUtr, ptr())
+ ## },
+ ## V = function() {
+ ## 'returns the weighted fixed-effects model matrix'
+ ## .Call(merPredDV, ptr())
+ ## },
+ ## VtV = function() {
+ ## 'returns the weighted cross-product of the fixed-effects model matrix'
+ ## .Call(merPredDVtV, ptr())
+ ## },
+ ## Vtr = function() {
+ ## 'returns the weighted cross-product of the fixed-effects model matrix\nand the residuals'
+ ## .Call(merPredVtr, ptr())
+ ## },
b = function(fac) {
'random effects on original scale for step factor fac'
.Call(merPredDb, ptr(), as.numeric(fac))
@@ -93,14 +123,14 @@
'fixed-effects coefficients for step factor fac'
.Call(merPredDbeta, ptr(), as.numeric(fac))
},
- delb = function() {
- 'increment for the fixed-effects coefficients'
- .Call(merPredDdelb, ptr())
- },
- delu = function() {
- 'increment for the orthogonal random-effects coefficients'
- .Call(merPredDdelu, ptr())
- },
+ ## delb = function() {
+ ## 'increment for the fixed-effects coefficients'
+ ## .Call(merPredDdelb, ptr())
+ ## },
+ ## delu = function() {
+ ## 'increment for the orthogonal random-effects coefficients'
+ ## .Call(merPredDdelu, ptr())
+ ## },
ldL2 = function() {
'twice the log determinant of the sparse Cholesky factor'
.Call(merPredDldL2, ptr())
@@ -125,7 +155,9 @@
'returns the external pointer, regenerating if necessary'
if (length(theta)) {
if (.Call(isNullExtPtr, Ptr))
- Ptr <<- .Call(merPredDCreate, X, Zt, Lambdat, Lind, theta, u0, beta0)
+ Ptr <<- .Call(merPredDCreate, X, Lambdat, LamtUt, Lind,
+ RZX, Ut, Utr, V, VtV, Vtr, Zt,
+ beta0, delb, delu, theta, u0)
}
Ptr
},
@@ -171,7 +203,8 @@
}
)
)
-merPredD$lock("X", "Zt", "Lind", "Lambdat", "u0", "beta0", "theta")
+merPredD$lock("Lambdat", "LamtUt", "Lind", "Ptr", "RZX", "Ut", "Utr", "V", "VtV",
+ "Vtr", "X", "Zt", "beta0", "delb", "delu", "theta", "u0")
lmResp <- # base class for response modules
setRefClass("lmResp",
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2011-12-06 17:42:13 UTC (rev 1471)
+++ pkg/lme4Eigen/R/lmer.R 2011-12-06 23:25:08 UTC (rev 1472)
@@ -41,17 +41,15 @@
reTrms <- mkReTrms(findbars(formula[[3]]), fr)
if (any(unlist(lapply(reTrms$flist, nlevels)) >= nrow(fr)))
stop("number of levels of each grouping factor must be less than number of obs")
- ## fixed-effects model matrix X - remove random parts from formula:
+ ## fixed-effects model matrix X - remove random effects from formula:
form <- formula
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)
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")]))
-# pp <- do.call(merPredD$new, c(list(X=X, Theta=reTrms$theta), reTrms[c("Zt","Lambdat","Lind")]))
+ pp <- do.call(merPredD$new, c(list(X=X, S=1L), reTrms[c("Zt","theta","Lambdat","Lind")]))
resp <- mkRespMod2(fr, if(REML) p else 0L)
-# if (REML) resp$reml <- p
devfun <- mkdevfun(pp, resp, emptyenv())
if (devFunOnly) return(devfun)
@@ -80,7 +78,7 @@
sigmaML=sqrt(pwrss/n), sigmaREML=sqrt(pwrss/(n-p)))
new("lmerMod", call=mc, frame=fr, flist=reTrms$flist, cnms=reTrms$cnms,
- theta=pp$theta, beta=pp$delb(), u=pp$delu(), lower=reTrms$lower, Gp=reTrms$Gp,
+ theta=pp$theta, beta=pp$delb, u=pp$delu, lower=reTrms$lower, Gp=reTrms$Gp,
devcomp=list(cmp=cmp, dims=dims), pp=pp, resp=resp)
}## { lmer }
@@ -112,7 +110,7 @@
control = list(), start = NULL, verbose = 0L, nAGQ = 1L,
doFit = TRUE, compDev=TRUE, subset, weights, na.action, offset,
contrasts = NULL, mustart, etastart, devFunOnly = 0L,
- tolPwrss = 0.000001, ...)
+ tolPwrss = 1e-10, ...)
{
verbose <- as.integer(verbose)
mf <- mc <- match.call()
@@ -1045,7 +1043,7 @@
### FIXME: Should modify the call slot to set REML=FALSE. It is
### tricky to do so without causing the call to be evaluated
new("lmerMod", call=x at call, frame=x at frame, flist=x at flist,
- cnms=x at cnms, theta=pp$theta, beta=pp$delb(), u=pp$delu(),
+ cnms=x at cnms, theta=pp$theta, beta=pp$delb, u=pp$delu,
lower=x at lower, devcomp=list(cmp=cmp, dims=dims), pp=pp, resp=resp)
}
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2011-12-06 17:42:13 UTC (rev 1471)
+++ pkg/lme4Eigen/src/external.cpp 2011-12-06 23:25:08 UTC (rev 1472)
@@ -398,11 +398,14 @@
// dense predictor module for mixed-effects models
- SEXP merPredDCreate(SEXP Xs, SEXP Zt, SEXP Lambdat, SEXP Lind, SEXP theta,
- SEXP u0, SEXP beta0) {
+ SEXP merPredDCreate(SEXP Xs, SEXP Lambdat, SEXP LamtUt, SEXP Lind,
+ SEXP RZX, SEXP Ut, SEXP Utr, SEXP V, SEXP VtV,
+ SEXP Vtr, SEXP Zt, SEXP beta0, SEXP delb, SEXP delu,
+ SEXP theta, SEXP u0) {
BEGIN_RCPP;
S4 X(Xs);
- merPredD *ans = new merPredD(X, Zt, Lambdat, Lind, theta, u0, beta0);
+ merPredD *ans = new merPredD(X, Lambdat, LamtUt, Lind, RZX, Ut, Utr,
+ V, VtV, Vtr, Zt, beta0, delb, delu, theta, u0);
return wrap(XPtr<merPredD>(ans, true));
END_RCPP;
}
@@ -741,7 +744,7 @@
CALLDEF(lmer_Deviance, 3), // method
CALLDEF(lmer_Laplace, 4),
- CALLDEF(merPredDCreate, 7), // generate external pointer
+ CALLDEF(merPredDCreate, 16), // generate external pointer
CALLDEF(merPredDsetTheta, 2), // setters
CALLDEF(merPredDsetBeta0, 2),
Modified: pkg/lme4Eigen/src/predModule.cpp
===================================================================
--- pkg/lme4Eigen/src/predModule.cpp 2011-12-06 17:42:13 UTC (rev 1471)
+++ pkg/lme4Eigen/src/predModule.cpp 2011-12-06 23:25:08 UTC (rev 1472)
@@ -10,11 +10,14 @@
namespace lme4Eigen {
using std::invalid_argument;
+ typedef Map<MatrixXd> MMat;
typedef Map<VectorXd> MVec;
typedef Map<VectorXi> MiVec;
- merPredD::merPredD(S4 X, SEXP Zt, SEXP Lambdat, SEXP Lind,
- SEXP theta, SEXP u0, SEXP beta0)
+ merPredD::merPredD(S4 X, SEXP Lambdat, SEXP LamtUt, SEXP Lind,
+ SEXP RZX, SEXP Ut, SEXP Utr, SEXP V, SEXP VtV,
+ SEXP Vtr, SEXP Zt, SEXP beta0, SEXP delb, SEXP delu,
+ SEXP theta, SEXP u0)
: d_X( X),
d_Zt( as<MSpMatrixd>(Zt)),
d_theta( as<MVec>(theta)),
@@ -24,18 +27,19 @@
d_p( d_X.cols()),
d_q( d_Zt.rows()),
d_Lambdat( as<MSpMatrixd>(Lambdat)),
- d_RZX( d_q, d_p),
- d_V( d_X.rows(), d_X.cols()),
- d_VtV( d_p,d_p),
- d_Vtr( d_p),
- d_Utr( d_q),
- d_delb( d_p),
- d_delu( d_q),
+ d_RZX( as<MMat>(RZX)),
+ d_V( as<MMat>(V)),
+ d_VtV( as<MMat>(VtV)),
+ d_Vtr( as<MVec>(Vtr)),
+ d_Utr( as<MVec>(Utr)),
d_beta0( as<MVec>(beta0)),
+ d_delb( as<MVec>(delb)),
+ d_delu( as<MVec>(delu)),
d_u0( as<MVec>(u0)),
- d_Ut( d_Zt),
- d_RX( d_p),
- d_LamtUtRestructure(false)
+ d_Ut( as<MSpMatrixd>(Ut)),
+ d_LamtUt( as<MSpMatrixd>(LamtUt)),
+ d_RX( d_p)
+ // d_LamtUtRestructure(false)
{ // Check consistency of dimensions
if (d_n != d_Zt.cols())
throw invalid_argument("Z dimension mismatch");
@@ -43,42 +47,32 @@
throw invalid_argument("size of Lind does not match nonzeros in Lambda");
// checking of the range of Lind is now done in R code for reference class
// initialize beta0, u0, delb, delu and VtV
- d_beta0.setZero(); d_u0.setZero(); d_delu.setZero(); d_delb.setZero();
- d_V = d_X;
+// d_beta0.setZero(); d_u0.setZero(); d_delu.setZero(); d_delb.setZero();
+// d_V = d_X;
d_VtV.setZero().selfadjointView<Eigen::Upper>().rankUpdate(d_V.adjoint());
d_RX.compute(d_VtV); // ensure d_RX is initialized even in the 0-column X case
setTheta(d_theta); // starting values into Lambda
d_L.cholmod().final_ll = 1; // force an LL' decomposition
- d_LamtUt = d_Lambdat * d_Ut;
+ updateLamtUt();
d_L.analyzePattern(d_LamtUt); // perform symbolic analysis
if (d_L.info() != Eigen::Success)
throw runtime_error("CholeskyDecomposition.analyzePattern failed");
}
void merPredD::updateLamtUt() {
- if (d_LamtUtRestructure) {
- d_LamtUt = d_Lambdat * d_Ut;
- return;
- }
- // This code fails if Lambdat is a MappedSparseMatrix<double>,
- // which is why it is stored as a SparseMatrix<double>.
// 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.
std::fill(d_LamtUt._valuePtr(), d_LamtUt._valuePtr() + d_LamtUt.nonZeros(), Scalar());
-// Rcpp::Rcout << "Lambdat[" << d_Lambdat.rows() << ", " << d_Lambdat.cols() << "]: "
-// << "innerSize() = " << d_Lambdat.innerSize()
-// << ", outerSize() = " << d_Lambdat.outerSize() << std::endl;
for (Index j = 0; j < d_Ut.cols(); ++j) {
-// Rcpp::Rcout << "d_Ut column index: " << j << std::endl;
- for(SpMatrixd::InnerIterator rhsIt(d_Ut, j); rhsIt; ++rhsIt) {
+ for(MSpMatrixd::InnerIterator rhsIt(d_Ut, j); rhsIt; ++rhsIt) {
Scalar y(rhsIt.value());
Index k(rhsIt.index());
// Rcpp::Rcout << "d_Ut row k = " << k << ", v = " << y;
// Rcpp::Rcout << ", Lambdat: col has " << d_Lambdat.innerNonZeros(k)
// << " nonzeros" << std::endl;
- SpMatrixd::InnerIterator prdIt(d_LamtUt, j);
+ MSpMatrixd::InnerIterator prdIt(d_LamtUt, j);
for (MSpMatrixd::InnerIterator lhsIt(d_Lambdat, k); lhsIt; ++lhsIt) {
Index i = lhsIt.index();
// Rcpp::Rcout << "lhsIt: i = " << i << ", v = " << lhsIt.value() << std::endl;
@@ -109,7 +103,7 @@
void merPredD::updateL() {
updateLamtUt();
- d_L.factorize_p(d_LamtUt, ArrayXi(), 1.);
+ d_L.factorize_p(d_LamtUt, Eigen::ArrayXi(), 1.);
d_ldL2 = ::M_chm_factor_ldetL2(d_L.factor());
}
@@ -162,7 +156,10 @@
throw invalid_argument("updateXwts: dimension mismatch");
if (sqrtXwt.size() == d_V.rows()) {
d_V = sqrtXwt.asDiagonal() * d_X;
- d_Ut = d_Zt * sqrtXwt.asDiagonal();
+ 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];
} else {
int n = d_X.rows()/d_V.rows();
SpMatrixd W(n, sqrtXwt.size());
@@ -174,11 +171,12 @@
}
W.finalize();
d_V = W * d_X;
- d_Ut = d_Zt * W.adjoint();
+//FIXME: work out the corresponding code for Ut and Zt
+// 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());
+// if (d_LamtUt.rows() != d_Ut.rows() ||
+// d_LamtUt.cols() != d_Ut.cols()) d_LamtUtRestructure = true;
+ d_VtV.setZero().selfadjointView<Eigen::Upper>().rankUpdate(d_V.adjoint());
updateL();
}
Modified: pkg/lme4Eigen/src/predModule.h
===================================================================
--- pkg/lme4Eigen/src/predModule.h 2011-12-06 17:42:13 UTC (rev 1471)
+++ pkg/lme4Eigen/src/predModule.h 2011-12-06 23:25:08 UTC (rev 1472)
@@ -12,17 +12,13 @@
#include <RcppEigen.h>
namespace lme4Eigen {
- using Eigen::ArrayXi;
using Eigen::CholmodDecomposition;
- using Eigen::Dynamic;
using Eigen::LLT;
using Eigen::Lower;
using Eigen::Map;
- using Eigen::Matrix;
using Eigen::MatrixXd;
using Eigen::SelfAdjointView;
using Eigen::SparseMatrix;
- using Eigen::Upper;
using Eigen::VectorXd;
using Eigen::VectorXi;
@@ -137,15 +133,15 @@
Index d_n, d_nnz, d_p, d_q;
MSpMatrixd d_Lambdat;
Scalar d_CcNumer, d_ldL2, d_ldRX2;
- MatrixXd d_RZX, d_V, d_VtV;
- VectorXd d_Vtr, d_Utr, d_delb, d_delu;
- Map<VectorXd> d_beta0, d_u0;
- SpMatrixd d_Ut, d_LamtUt;
+ Map<MatrixXd> d_RZX, d_V, d_VtV;
+ Map<VectorXd> d_Vtr, d_Utr, d_beta0, d_delb, d_delu, d_u0;
+ MSpMatrixd d_Ut, d_LamtUt;
ChmDecomp d_L;
LLT<MatrixXd> d_RX;
- bool d_LamtUtRestructure;
+// bool d_LamtUtRestructure;
public:
- merPredD(S4, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
+ merPredD(S4, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP,
+ SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
IntegerVector Pvec() const;
@@ -168,22 +164,22 @@
const ChmDecomp& L() const {return d_L;}
- const MatrixXd& V() const {return d_V;}
- const MatrixXd& VtV() const {return d_VtV;}
- const MatrixXd& RZX() const {return d_RZX;}
+ const Map<MatrixXd>& V() const {return d_V;}
+ const Map<MatrixXd>& VtV() const {return d_VtV;}
+ const Map<MatrixXd>& RZX() const {return d_RZX;}
const Map<VectorXd>& theta() const {return d_theta;}
const MSpMatrixd& Lambdat() const {return d_Lambdat;}
- const SpMatrixd& LamtUt() const {return d_LamtUt;}
- const SpMatrixd& Ut() const {return d_Ut;}
+ const MSpMatrixd& LamtUt() const {return d_LamtUt;}
+ const MSpMatrixd& Ut() const {return d_Ut;}
const MSpMatrixd& Zt() const {return d_Zt;}
- const VectorXd& Utr() const {return d_Utr;}
- const VectorXd& Vtr() const {return d_Vtr;}
- const VectorXd& delb() const {return d_delb;}
- const VectorXd& delu() const {return d_delu;}
+ const Map<VectorXd>& Utr() const {return d_Utr;}
+ const Map<VectorXd>& Vtr() const {return d_Vtr;}
+ const Map<VectorXd>& delb() const {return d_delb;}
+ const Map<VectorXd>& delu() const {return d_delu;}
const Map<VectorXd>& beta0() const {return d_beta0;}
const Map<VectorXd>& u0() const {return d_u0;}
More information about the Lme4-commits
mailing list