Environment-based lme4

intro.md

Background

The current version of the (lme4)[https://github.com/lme4/lme4] package for (R)[http://www.r-project.org) is based on a rather questionable practice of storing external pointers from C++ objects in an R reference class object. The idea was that when the reference class object was created the corresponding C++ class instance would be created. All changes to this object would be made through the external pointer in the C++ code.

This is risky because any changes to the R objects in the reference class could cause changes in the location of its values, after which all bets are off because the R object and the C++ object are referring to different memory locations.

It turns out that this does happen in the development version of R, which will become R-3.1.0, if the LAZY_DUPLICATE_OK flag is set. The symptom is that deepcopy methods are not copying and I suspect this is because the lazy duplication doesn't duplicate before the external pointer is constructed. We were living outside the law and we got caught.

Suggested changes

I think we should abandon the reference classes and go back to explicitly using an environment to hold the various arrays representing the model. The point is that C++ code using Rcpp templates can access R objects in an environment and can assign new objects that then become visible in R. For example, to update the sparse Cholesky decomposition for a new value of theta, the beginning of a method is

update_Lambdat.cpp

// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-

#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

typedef Eigen::MappedSparseMatrix<double>         MSpMat;
typedef Eigen::SparseMatrix<double>                SpMat;
typedef Eigen::SimplicialLDLT<SpMat>              SpChol;
typedef Eigen::MatrixXd                              Mat;
typedef Eigen::VectorXd                              Vec;
typedef Eigen::Map<Mat>                             MMat;
typedef Eigen::Map<Vec>                             MVec;
typedef Eigen::VectorXi                             iVec;
typedef Eigen::Map<iVec>                           MiVec;
typedef Eigen::CholmodDecomposition<SpMat>     ChmDecomp;

using Rcpp::as;
using Rcpp::wrap;
using Rcpp::Environment;

// [[Rcpp::export]]
ChmDecomp update_predD(Environment rho) {
  // Calls to rho.get can be preceded by calls to rho.exists to make sure the name is available
  MSpMat        Zt(as<MSpMat>(rho.get("Zt")));
  MSpMat   Lambdat(as<MSpMat>(rho.get("Lambdat")));
  MiVec       Lind(as<MiVec>(rho.get("Lind")));
  MVec       theta(as<MVec>(rho.get("theta")));
  if (Lind.size() != Lambdat.nonZeros()) throw std::invalid_argument("dimension mismatch");
  int    *lipt = Lind.data();
	double *LamX = Lambdat.valuePtr(), *thpt = theta.data();
	for (int i = 0; i < Lind.size(); ++i) {
	    LamX[i] = thpt[lipt[i] - 1];
	}
  SpMat  LamtZt(Lambdat * Zt);
  SpMat  LamtZtZLam(LamtZt * LamtZt.adjoint());
  rho.assign("ZtZ",wrap(LamtZtZLam));
  return ChmDecomp(LamtZtZLam);//  ChmDecomp      L = as<ChmDecomp>(rho.get("L"));
}