[Lme4-commits] r1560 - in pkg/lme4Eigen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 2 23:35:07 CET 2012
Author: dmbates
Date: 2012-02-02 23:35:07 +0100 (Thu, 02 Feb 2012)
New Revision: 1560
Modified:
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/R/profile.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/predModule.cpp
pkg/lme4Eigen/src/predModule.h
Log:
rePos reference class and its use to calculate the postVar attribute for ranef
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2012-02-02 22:32:56 UTC (rev 1559)
+++ pkg/lme4Eigen/R/AllClass.R 2012-02-02 22:35:07 UTC (rev 1560)
@@ -33,16 +33,16 @@
##'
##'
##' @param X dense model matrix for the fixed-effects parameters, to be stored
-##' in the \code{X} field.
+##' in the \code{X} field.
##' @param Zt transpose of the sparse model matrix for the random effects. It
-##' is stored in the \code{Zt} field.
+##' is stored in the \code{Zt} field.
##' @param Lambdat transpose of the sparse lower triangular relative variance
-##' factor (stored in the \code{Lambdat} field).
+##' factor (stored in the \code{Lambdat} field).
##' @param Lind integer vector of the same length as the \code{"x"} slot in the
-##' \code{Lambdat} field. Its elements should be in the range 1 to the length
-##' of the \code{theta} field.
+##' \code{Lambdat} field. Its elements should be in the range 1 to the length
+##' of the \code{theta} field.
##' @param theta numeric vector of variance component parameters (stored in the
-##' \code{theta} field).
+##' \code{theta} field).
##' @section Methods:
##' \describe{
##' \item{new(X, Zt, Lambdat, Lind, theta):}{Create a new \code{\linkS4class{merPredD}} object}
@@ -876,3 +876,63 @@
##' @export
setClass("nlmerMod", representation(resp="nlsResp"), contains="merMod")
+
+##' Generator object for the rePos (random-effects positions) class
+##'
+##' The generator object for the \code{\linkS4class{rePos}} class used
+##' to determine the positions and orders of random effects associated
+##' with particular random-effects terms in the model.
+##' @param mer an object of class \code{"\linkS4class{merMod}"}
+##' @note Arguments to the \code{new} methods must be named arguments.
+##' @section Methods:
+##' \describe{
+##' \item{\code{new(mer=mer)}}{Create a new
+##' \code{\linkS4class{rePos}} object.}
+##' }
+##' @seealso \code{\linkS4class{rePos}}
+##' @keywords classes
+##' @export
+rePos <-
+ setRefClass("rePos",
+ fields =
+ list(
+ cnms = "list",
+ flist = "list",
+ ncols = "integer",
+ nctot = "integer",
+ nlevs = "integer",
+ offsets = "integer",
+ terms = "list"
+ ),
+ methods =
+ list(
+ initialize = function(mer, ...) {
+ stopifnot((ntrms <- length(Cnms <- mer at cnms)) > 0L,
+ (length(Flist <- mer at flist)) > 0L,
+ length(asgn <- as.integer(attr(Flist, "assign"))) == ntrms)
+ cnms <<- Cnms
+ flist <<- Flist
+ ncols <<- unname(vapply(cnms, length, 0L))
+ nctot <<- unname(as.vector(tapply(ncols, asgn, sum)))
+ nlevs <<- unname(vapply(flist, function(el) length(levels(el)), 0L))
+ offsets <<- c(0L, cumsum(sapply(seq_along(asgn),
+ function(i) ncols[i] * nlevs[asgn[i]])))
+ terms <<- lapply(seq_along(flist), function(i) which(asgn == i))
+ }
+ )
+ )
+##' Class \code{"rePos"}
+##'
+##' A reference class for determining the positions in the random-effects vector
+##' that correspond to particular random-effects terms in the model formula
+##'
+##' @name rePos-class
+##' @docType class
+##' @section Extends: All reference classes extend and inherit methods from
+##' \code{"\linkS4class{envRefClass}"}.
+##' @keywords classes
+##' @examples
+##'
+##' showClass("rePos")
+##'
+rePos$lock("cnms", "flist", "ncols", "nctot", "nlevs", "terms")
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-02-02 22:32:56 UTC (rev 1559)
+++ pkg/lme4Eigen/R/lmer.R 2012-02-02 22:35:07 UTC (rev 1560)
@@ -1157,9 +1157,10 @@
ans <- ans[whchL]
if (postVar) {
- vv <- .Call(merPredDcondVar, object at pp$ptr(), sigma(object), object at flist)
+ sigsqr <- sigma(object)^2
+ vv <- .Call(merPredDcondVar, object at pp$ptr(), as.environment(rePos$new(object)))
for (i in seq_along(ans))
- attr(ans[[i]], "postVar") <- vv[[i]]
+ attr(ans[[i]], "postVar") <- vv[[i]] * sigsqr
}
if (drop)
ans <- lapply(ans, function(el)
Modified: pkg/lme4Eigen/R/profile.R
===================================================================
--- pkg/lme4Eigen/R/profile.R 2012-02-02 22:32:56 UTC (rev 1559)
+++ pkg/lme4Eigen/R/profile.R 2012-02-02 22:35:07 UTC (rev 1560)
@@ -680,6 +680,7 @@
##' @title Transform to the variance scale
##' @param pr a mixed-effects model profile
##' @return a transformed mixed-effects model profile
+##' @export
varianceProf <- function(pr) {
stopifnot(inherits(pr, "thpr"))
spl <- attr(pr, "forward")
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2012-02-02 22:32:56 UTC (rev 1559)
+++ pkg/lme4Eigen/src/external.cpp 2012-02-02 22:35:07 UTC (rev 1560)
@@ -33,7 +33,6 @@
using lme4Eigen::lmerResp;
using lme4Eigen::merPredD;
using lme4Eigen::nlsResp;
- using lme4Eigen::reTrms;
using optimizer::Golden;
using optimizer::Nelder_Mead;
@@ -577,10 +576,9 @@
END_RCPP;
}
- SEXP merPredDcondVar(SEXP ptr, SEXP sigma, SEXP cnms, SEXP flist) {
+ SEXP merPredDcondVar(SEXP ptr, SEXP rho) {
BEGIN_RCPP;
- reTrms trms(cnms, flist);
- return wrap(XPtr<merPredD>(ptr)->condVar(::Rf_asReal(sigma), trms));
+ return wrap(XPtr<merPredD>(ptr)->condVar(Rcpp::Environment(rho)));
END_RCPP;
}
@@ -792,37 +790,6 @@
return ::Rf_ScalarReal(XPtr<nlsResp>(ptr_)->updateMu(as<MVec>(gamma)));
END_RCPP;
}
-
- SEXP reTrms_Create(SEXP cnms, SEXP flist) {
- BEGIN_RCPP;
- reTrms *ans = new reTrms(cnms, flist);
- return wrap(XPtr<reTrms>(ans, true));
- END_RCPP;
- }
-
- SEXP reTrms_ncols(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<reTrms>(ptr_)->ncols());
- END_RCPP;
- }
-
- SEXP reTrms_nctot(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<reTrms>(ptr_)->nctot());
- END_RCPP;
- }
-
- SEXP reTrms_nlevs(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<reTrms>(ptr_)->nlevs());
- END_RCPP;
- }
-
- SEXP reTrms_offsets(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<reTrms>(ptr_)->offsets());
- END_RCPP;
- }
}
#include <R_ext/Rdynload.h>
@@ -908,7 +875,7 @@
CALLDEF(merPredDb, 2), // methods
CALLDEF(merPredDbeta, 2),
- CALLDEF(merPredDcondVar, 4),
+ CALLDEF(merPredDcondVar, 2),
CALLDEF(merPredDlinPred, 2),
CALLDEF(merPredDinstallPars, 2),
CALLDEF(merPredDsolve, 1),
@@ -940,11 +907,6 @@
CALLDEF(nls_Laplace, 4), // methods
CALLDEF(nls_updateMu, 2),
- CALLDEF(reTrms_Create, 2),
- CALLDEF(reTrms_ncols, 1),
- CALLDEF(reTrms_nctot, 1),
- CALLDEF(reTrms_nlevs, 1),
- CALLDEF(reTrms_offsets, 1),
{NULL, NULL, 0}
};
Modified: pkg/lme4Eigen/src/predModule.cpp
===================================================================
--- pkg/lme4Eigen/src/predModule.cpp 2012-02-02 22:32:56 UTC (rev 1559)
+++ pkg/lme4Eigen/src/predModule.cpp 2012-02-02 22:35:07 UTC (rev 1560)
@@ -88,55 +88,34 @@
return d_X * beta(f) + d_Zt.adjoint() * b(f);
}
- Rcpp::List merPredD::condVar(const double& scale, const reTrms& trm) const {
- if (scale < 0 || !R_finite(scale))
- throw runtime_error("scale must be non-negative and finite");
- Eigen::VectorXi nc = trm.ncols(), nl = trm.nlevs(), nct = trm.nctot(),
- off = trm.offsets();
- const Rcpp::List ll(trm.flist());
- int nf(ll.size());
+ Rcpp::List merPredD::condVar(const Rcpp::Environment& rho) const {
+ const Rcpp::List ll(as<Rcpp::List>(rho["flist"])), trmlst(as<Rcpp::List>(rho["terms"]));
+ const int nf(ll.size());
+ const MiVec nc(as<MiVec>(rho["ncols"])), nl(as<MiVec>(rho["nlevs"])),
+ nct(as<MiVec>(rho["nctot"])), off(as<MiVec>(rho["offsets"]));
Rcpp::List ans(nf);
- Rcpp::CharacterVector nms = ll.names();
- ans.names() = clone(nms);
+ ans.names() = clone(as<Rcpp::CharacterVector>(ll.names()));
-#if 0
+ const SpMatrixd d_Lambda(d_Lambdat.adjoint());
for (int i = 0; i < nf; i++) {
- int ncti = nct[i], nli = nl[i];
- VectorXi trms(trm.terms(i));
- int *cset = new int[ncti], nct2 = ncti * ncti;
-
- NumericVector ansi(Dimension(ncti, ncti, nli));
+ int ncti(nct[i]), nli(nl[i]);
+ Rcpp::NumericVector ansi(Rcpp::Dimension(ncti, ncti, nli));
ans[i] = ansi;
- double *ai = ansi.begin();
-
- for (int j = 0; j < nli; j++) {
- int kk = 0;
- for (int jj = 0; jj < trms.size(); jj++) {
- int tjj = trms[jj];
- for (int k = 0; k < nc[tjj]; k++)
- cset[kk++] = off[tjj] + j * nc[tjj] + k;
+ const MiVec trms(as<MiVec>(trmlst(i)));
+ if (trms.size() == 1) { // simple case
+ int offset = off[trms[0] - 1];
+ for (int j = 0; j < nli; ++j) {
+ MatrixXd Lv(d_Lambda.innerVectors(offset + j * ncti, ncti));
+ d_L.solveInPlace(Lv, CHOLMOD_A);
+ MatrixXd rr(MatrixXd(ncti, ncti).setZero().
+ selfadjointView<Eigen::Lower>().rankUpdate(Lv.adjoint()));
+ std::copy(rr.data(), rr.data() + rr.size(), &ansi[j * ncti * ncti]);
}
- CHM_SP cols =
- M_cholmod_submatrix(&d_Lambda, (int*)NULL, -1, cset, ncti,
- 1/*values*/, 1/*sorted*/, &c);
- CHM_SP sol = d_L.spsolve(CHOLMOD_A, cols);
- CHM_SP tcols = M_cholmod_transpose(cols, 1/*values*/, &c);
- M_cholmod_free_sparse(&cols, &c);
- CHM_SP var = M_cholmod_ssmult(tcols, sol, 0/*stype*/,
- 1/*values*/, 1/*sorted*/, &c);
- M_cholmod_free_sparse(&sol, &c);
- M_cholmod_free_sparse(&tcols, &c);
- CHM_DN dvar = M_cholmod_sparse_to_dense(var, &c);
- M_cholmod_free_sparse(&var, &c);
- Memcpy(ai + j * nct2, (double*)dvar->x, nct2);
- M_cholmod_free_dense(&dvar, &c);
+ } else {
+ throw std::runtime_error("multiple terms per factor not yet written");
}
- delete[] cset;
- transform(ansi.begin(), ansi.end(), ansi.begin(),
- bind2nd(multiplies<double>(), scale * scale));
}
-#endif
return ans;
}
Modified: pkg/lme4Eigen/src/predModule.h
===================================================================
--- pkg/lme4Eigen/src/predModule.h 2012-02-02 22:32:56 UTC (rev 1559)
+++ pkg/lme4Eigen/src/predModule.h 2012-02-02 22:35:07 UTC (rev 1560)
@@ -10,7 +10,6 @@
#define LME4_PREDMODULE_H
#include <RcppEigen.h>
-#include "reTrms.h"
namespace lme4Eigen {
#if 0
@@ -140,7 +139,7 @@
VectorXd linPred(const Scalar& f) const;
VectorXd u(const Scalar& f) const;
- Rcpp::List condVar(const Scalar&, const reTrms&) const;
+ Rcpp::List condVar(const Rcpp::Environment&) const;
Scalar CcNumer() const {return d_CcNumer;}
Scalar ldL2() const {return d_ldL2;}
More information about the Lme4-commits
mailing list