[Lme4-commits] r1496 - in pkg/lme4Eigen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jan 7 01:36:31 CET 2012
Author: dmbates
Date: 2012-01-07 01:36:31 +0100 (Sat, 07 Jan 2012)
New Revision: 1496
Modified:
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/predModule.cpp
Log:
More work on nlmer.
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2012-01-06 15:00:07 UTC (rev 1495)
+++ pkg/lme4Eigen/R/AllClass.R 2012-01-07 00:36:31 UTC (rev 1496)
@@ -35,9 +35,10 @@
u0 = "numeric"),
methods =
list(
+### FIXME: probably don't need S as Xwts is required for nlmer
initialize = function(X, Zt, Lambdat, Lind, theta, S, ...) {
if (!nargs()) return
- X <<- X
+ X <<- as(X, "matrix")
Zt <<- as(Zt, "dgCMatrix")
Lambdat <<- as(Lambdat, "dgCMatrix")
Lind <<- as.integer(Lind)
@@ -57,10 +58,12 @@
V <<- array(0, c(n, p))
VtV <<- array(0, c(p, p))
Vtr <<- numeric(p)
- beta0 <<- numeric(p)
+ b0 <- list(...)$beta0
+ beta0 <<- if (is.null(b0)) numeric(p) else b0
delb <<- numeric(p)
delu <<- numeric(q)
- u0 <<- numeric(q)
+ uu <- list(...)$u0
+ u0 <<- if (is.null(uu)) numeric(q) else uu
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
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-01-06 15:00:07 UTC (rev 1495)
+++ pkg/lme4Eigen/R/lmer.R 2012-01-07 00:36:31 UTC (rev 1496)
@@ -97,6 +97,23 @@
.Call(glmerAGQ, pp$ptr(), resp$ptr(), fac, GQmat, pars[dpars],
u0, pars[-dpars], tolPwrss)
}
+ } else if (is(rho$resp, "nlsResp")) {
+ if (nAGQ < 2L) {
+ rho$nlmerLaplace <- nlmerLaplace
+ ff <- switch(nAGQ + 1L,
+ function(theta)
+ .Call(nlmerLaplace, pp$ptr(), resp$ptr(), theta,
+ u0, beta0, verbose, FALSE, tolPwrss),
+ function(pars)
+ .Call(nlmerLaplace, pp$ptr(), resp$ptr(), pars[dpars], u0,
+ pars[-dpars], verbose, TRUE, tolPwrss))
+ } else {
+ rho$nlmerAGQ <- nlmerAGQ
+ rho$GQmat <- GHrule(nAGQ)
+ ff <- function(pars)
+ .Call(nlmerAGQ, pp$ptr(), resp$ptr(), fac, GQmat, pars[dpars],
+ u0, pars[-dpars], tolPwrss)
+ }
}
if (is.null(ff)) stop("code not yet written")
environment(ff) <- rho
@@ -668,7 +685,8 @@
nlmer <- function(formula, data, family = gaussian, start = NULL,
verbose = 0L, nAGQ = 1L, doFit = TRUE,
subset, weights, na.action, offset,
- contrasts = NULL, devFunOnly=0L, ...)
+ contrasts = NULL, devFunOnly=0L,
+ tolPwrss = 1e-10, optimizer=c("bobyqa","NelderMead"), ...)
{
if (!missing(family)) stop("code not yet written")
mf <- mc <- match.call()
@@ -724,11 +742,18 @@
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, S=s, Xwts = rr$sqrtXwt,
- beta0=qr.coef(qrX, unlist(lapply(pnames, get, envir = nlenv))))
- ))
- return(list(pp=pp, resp=rr))
+ rho <- as.environment(list(verbose=verbose, tolPwrss=tolPwrss))
+ parent.env(rho) <- parent.frame()
+ rho$pp <- do.call(merPredD$new, c(reTrms[c("Zt","theta","Lambdat","Lind")],
+ list(X=X, S=s, Xwts = rr$sqrtXwt,
+ beta0=qr.coef(qrX, unlist(lapply(pnames, get, envir = nlenv))))
+ ))
+ rho$resp <- rr
+ rho$u0 <- rho$pp$u0
+ rho$beta0 <- rho$pp$beta0
+ rho$lower <- reTrms$lower
+ devfun <- mkdevfun(rho, 0L) # deviance as a function of theta only
+ devfun
}
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2012-01-06 15:00:07 UTC (rev 1495)
+++ pkg/lme4Eigen/src/external.cpp 2012-01-07 00:36:31 UTC (rev 1496)
@@ -309,6 +309,60 @@
END_RCPP;
}
+ static inline double prss(lmResp *rp, merPredD *pp, double fac) {
+ return rp->wrss() + (fac ? pp->sqrL(fac) : pp->u0().squaredNorm());
+ }
+
+ void nstepFac(nlsResp *rp, merPredD *pp, int verb) {
+ double prss0(prss(rp, pp, 0.));
+
+ for (double fac = 1.; fac > 0.001; fac /= 2.) {
+ double prss1 = rp->updateMu(pp->linPred(fac)) + pp->sqrL(fac);
+ if (verb > 3)
+ ::Rprintf("pwrss0=%10g, diff=%10g, fac=%6.4f\n",
+ prss0, prss0 - prss1, fac);
+ if (prss1 < prss0) {
+ pp->installPars(fac);
+ return;
+ }
+ }
+ throw runtime_error("step factor reduced below 0.001 without reducing pwrss");
+ }
+
+#define MAXITER 30
+ static void prssUpdate(nlsResp *rp, merPredD *pp, int verb, bool uOnly, double tol) {
+ bool cvgd(false);
+ for (int it=0; it < MAXITER; ++it) {
+ rp->updateMu(pp->linPred(0.));
+ pp->updateXwts(rp->sqrtXwt());
+ pp->updateDecomp();
+ pp->updateRes(rp->wtres());
+ if ((uOnly ? pp->solveU() : pp->solve())/pwrss(rp, pp, 0.) < tol) {
+ cvgd = true;
+ break;
+ }
+ nstepFac(rp, pp, verb);
+ }
+ if (!cvgd) throw runtime_error("pwrss failed to converge in 30 iterations");
+ }
+
+ SEXP nlmerLaplace(SEXP pp_, SEXP rp_, SEXP theta_, SEXP u0_, SEXP beta0_,
+ SEXP verbose_, SEXP uOnly_, SEXP tol_) {
+ BEGIN_RCPP;
+
+ XPtr<nlsResp> rp(rp_);
+ XPtr<merPredD> pp(pp_);
+
+ pp->setTheta(as<MVec>(theta_));
+ pp->setU0(as<MVec>(u0_));
+ pp->setBeta0(as<MVec>(beta0_));
+ prssUpdate(rp, pp, ::Rf_asInteger(verbose_), ::Rf_asLogical(uOnly_),
+ ::Rf_asReal(tol_));
+ return ::Rf_ScalarReal(rp->Laplace(pp->ldL2(), pp->ldRX2(), pp->sqrL(1.)));
+
+ END_RCPP;
+ }
+
SEXP golden_Create(SEXP lower_, SEXP upper_) {
BEGIN_RCPP;
Golden *ans = new Golden(::Rf_asReal(lower_), ::Rf_asReal(upper_));
@@ -803,6 +857,8 @@
CALLDEF(NelderMead_xeval, 1),
CALLDEF(NelderMead_xpos, 1),
+ CALLDEF(nlmerLaplace, 8),
+
CALLDEF(nls_Create, 11), // generate external pointer
CALLDEF(nls_Laplace, 4), // methods
Modified: pkg/lme4Eigen/src/predModule.cpp
===================================================================
--- pkg/lme4Eigen/src/predModule.cpp 2012-01-06 15:00:07 UTC (rev 1495)
+++ pkg/lme4Eigen/src/predModule.cpp 2012-01-07 00:36:31 UTC (rev 1496)
@@ -52,8 +52,10 @@
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
+// Rcpp::Rcout << "before setTheta" << std::endl;
setTheta(d_theta); // starting values into Lambda
d_L.cholmod().final_ll = 1; // force an LL' decomposition
+// Rcpp::Rcout << "before updateLamtUt" << std::endl;
updateLamtUt();
d_L.analyzePattern(d_LamtUt); // perform symbolic analysis
if (d_L.info() != Eigen::Success)
@@ -64,20 +66,21 @@
// 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());
+ MVec(d_LamtUt._valuePtr(), d_LamtUt.nonZeros()).setZero();
for (Index j = 0; j < d_Ut.cols(); ++j) {
for(MSpMatrixd::InnerIterator rhsIt(d_Ut, j); rhsIt; ++rhsIt) {
- Scalar y(rhsIt.value());
- Index k(rhsIt.index());
+ Scalar y(rhsIt.value());
+ Index k(rhsIt.index());
MSpMatrixd::InnerIterator prdIt(d_LamtUt, j);
for (MSpMatrixd::InnerIterator lhsIt(d_Lambdat, k); lhsIt; ++lhsIt) {
- Index i = lhsIt.index();
+ Index i = lhsIt.index();
while (prdIt && prdIt.index() != i) ++prdIt;
if (!prdIt) throw runtime_error("logic error in updateLamtUt");
- prdIt.valueRef() += lhsIt.value() * y;
+ prdIt.valueRef() += lhsIt.value() * y;
}
}
}
+// Rcpp::Rcout << "end of updateLamtUt" << std::endl;
}
VectorXd merPredD::b(const double& f) const {return d_Lambdat.adjoint() * u(f);}
@@ -94,6 +97,9 @@
void merPredD::updateL() {
updateLamtUt();
+ // More complicated code to handle the case of zeros in
+ // potentially nonzero positions. The factorize_p method is
+ // for a SparseMatrix<double>, not a MappedSparseMatrix<double>.
SpMatrixd m(d_LamtUt.rows(), d_LamtUt.cols());
m.resizeNonZeros(d_LamtUt.nonZeros());
std::copy(d_LamtUt._valuePtr(),
@@ -150,9 +156,11 @@
}
void merPredD::updateXwts(const VectorXd& sqrtXwt) {
+// Rcpp::Rcout << "in updateXwts" << std::endl;
if (d_Xwts.size() != sqrtXwt.size())
throw invalid_argument("updateXwts: dimension mismatch");
std::copy(sqrtXwt.data(), sqrtXwt.data() + sqrtXwt.size(), d_Xwts.data());
+// Rcpp::Rcout << "after copy, d_Xwts = " << d_Xwts.adjoint() << std::endl;
if (sqrtXwt.size() == d_V.rows()) {
d_V = sqrtXwt.asDiagonal() * d_X;
for (int j = 0; j < d_N; ++j)
@@ -160,24 +168,45 @@
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());
+ SpMatrixd W(d_V.rows(), 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.insertBack(j % d_V.rows(), j) = *pt;
}
W.finalize();
+ // Rcpp::Rcout << "in nlmer block: W(" << W.rows() << "," << W.cols()
+ // << "), outer: " << MiVec(W._outerIndexPtr(), W.outerSize()).adjoint()
+ // << std::endl;
+ // Rcpp::Rcout << "inner: " << MiVec(W._innerIndexPtr(), W.nonZeros()).adjoint()
+ // << std::endl;
+ // Rcpp::Rcout << "values: " << MVec(W._valuePtr(), W.nonZeros()).adjoint()
+ // << std::endl;
d_V = W * d_X;
-// 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];
+ SpMatrixd Ut(d_Zt * W.adjoint());
+ // Rcpp::Rcout << "Ut(" << Ut.rows() << "," << Ut.cols()
+ // << "), outer: " << MiVec(Ut._outerIndexPtr(), Ut.outerSize()).adjoint()
+ // << std::endl;
+ // Rcpp::Rcout << "inner: " << MiVec(Ut._innerIndexPtr(), Ut.nonZeros()).adjoint()
+ // << std::endl;
+ // Rcpp::Rcout << "values: " << MVec(Ut._valuePtr(), Ut.nonZeros()).adjoint()
+ // << std::endl;
+ if (Ut.cols() != d_Ut.cols())
+ throw std::runtime_error("Size mismatch in updateXwts");
+
+ // More complex code to handle the pruning of zeros
+ MVec(d_Ut._valuePtr(), d_Ut.nonZeros()).setZero();
+ for (int j = 0; j < d_Ut.cols(); ++j) {
+ MSpMatrixd::InnerIterator lhsIt(d_Ut, j);
+ SpMatrixd::InnerIterator rhsIt(Ut, j);
+ Index k(rhsIt.index());
+ while (lhsIt && lhsIt.index() != k) ++lhsIt;
+ if (lhsIt.index() != k)
+ throw std::runtime_error("Pattern mismatch in updateXwts");
+ lhsIt.valueRef() = rhsIt.value();
+ }
}
-// 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();
}
More information about the Lme4-commits
mailing list