[Lme4-commits] r1473 - in pkg/lme4Eigen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 8 00:26:00 CET 2011
Author: dmbates
Date: 2011-12-08 00:25:59 +0100 (Thu, 08 Dec 2011)
New Revision: 1473
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:
Allow regeneration of pointers in merPredD objects.
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2011-12-06 23:25:08 UTC (rev 1472)
+++ pkg/lme4Eigen/R/AllClass.R 2011-12-07 23:25:59 UTC (rev 1473)
@@ -25,7 +25,8 @@
V = "matrix",
VtV = "matrix",
Vtr = "numeric",
- X = "ddenseModelMatrix",
+ X = "matrix",
+ Xwts = "numeric",
Zt = "dgCMatrix",
beta0 = "numeric",
delb = "numeric",
@@ -36,7 +37,7 @@
list(
initialize = function(X, Zt, Lambdat, Lind, theta, S, ...) {
if (!nargs()) return
- X <<- as(X, "ddenseModelMatrix")
+ X <<- X
Zt <<- as(Zt, "dgCMatrix")
Lambdat <<- as(Lambdat, "dgCMatrix")
Lind <<- as.integer(Lind)
@@ -51,21 +52,20 @@
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
+ RZX <<- array(0, c(q, p))
+ Utr <<- numeric(q)
+ 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)
+ u0 <<- numeric(q)
+ 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))
+ LamtUt <<- Lambdat %*% Ut
+ updateXwts(Xwts)
},
CcNumer = function() {
'returns the numerator of the orthogonality convergence criterion'
@@ -75,10 +75,6 @@
'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())
- ## },
P = function() {
'returns the permutation vector for the sparse Cholesky factor'
.Call(merPredDPvec, ptr())
@@ -87,34 +83,14 @@
'returns the dense downdated Cholesky factor for the fixed-effects parameters'
.Call(merPredDRX, ptr())
},
+ RXi = function() {
+ 'returns the inverse of the dense downdated Cholesky factor for the fixed-effects parameters'
+ .Call(merPredDRXi, ptr())
+ },
RXdiag = function() {
'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())
- ## },
b = function(fac) {
'random effects on original scale for step factor fac'
.Call(merPredDb, ptr(), as.numeric(fac))
@@ -123,14 +99,6 @@
'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())
- ## },
ldL2 = function() {
'twice the log determinant of the sparse Cholesky factor'
.Call(merPredDldL2, ptr())
@@ -154,10 +122,14 @@
ptr = function() {
'returns the external pointer, regenerating if necessary'
if (length(theta)) {
- if (.Call(isNullExtPtr, Ptr))
- Ptr <<- .Call(merPredDCreate, X, Lambdat, LamtUt, Lind,
- RZX, Ut, Utr, V, VtV, Vtr, Zt,
- beta0, delb, delu, theta, u0)
+ if (.Call(isNullExtPtr, Ptr)) {
+ Ptr <<- .Call(merPredDCreate, as(X, "matrix"), Lambdat,
+ LamtUt, Lind, RZX, Ut, Utr, V, VtV, Vtr,
+ Xwts, Zt, beta0, delb, delu, theta, u0)
+ .Call(merPredDsetTheta, Ptr, theta)
+ .Call(merPredDupdateXwts, Ptr, Xwts)
+ .Call(merPredDupdateDecomp, Ptr)
+ }
}
Ptr
},
@@ -203,8 +175,8 @@
}
)
)
-merPredD$lock("Lambdat", "LamtUt", "Lind", "Ptr", "RZX", "Ut", "Utr", "V", "VtV",
- "Vtr", "X", "Zt", "beta0", "delb", "delu", "theta", "u0")
+merPredD$lock("Lambdat", "LamtUt", "Lind", "RZX", "Ut", "Utr", "V", "VtV", "Vtr",
+ "X", "Xwts", "Zt", "beta0", "delb", "delu", "theta", "u0")
lmResp <- # base class for response modules
setRefClass("lmResp",
@@ -219,6 +191,11 @@
y = "numeric"),
methods =
list(
+ allInfo = function() {
+ 'return all the information available on the object'
+ data.frame(y=y, offset=offset, weights=weights, mu=mu,
+ rwt=sqrtrwt, wres=wtres, Xwt=sqrtXwt)
+ },
initialize = function(...) {
if (!nargs()) return()
ll <- list(...)
@@ -306,11 +283,10 @@
},
allInfo = function() {
'return all the information available on the object'
- data.frame(y=y, n=n, offset=offset, weights=weights,
- eta=eta, mu=mu, muEta=muEta(), variance=variance(),
- sqrtrwt=sqrtrwt, wtres=wtres, sqrtXwt=sqrtXwt,
- sqrtWrkWt=sqrtWrkWt(), wrkResids=wrkResids(),
- wrkResp=wrkResp(), devRes=devResid())
+ cbind(callSuper(),
+ data.frame(eta=eta, muEta=muEta(), var=variance(),
+ WrkWt=sqrtWrkWt(), wrkRes=wrkResids(),
+ wrkResp=wrkResp(), devRes=devResid()))
},
devResid = function() {
'returns the vector of deviance residuals'
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2011-12-06 23:25:08 UTC (rev 1472)
+++ pkg/lme4Eigen/R/lmer.R 2011-12-07 23:25:59 UTC (rev 1473)
@@ -44,7 +44,7 @@
## 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
+ 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))
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2011-12-06 23:25:08 UTC (rev 1472)
+++ pkg/lme4Eigen/src/external.cpp 2011-12-07 23:25:59 UTC (rev 1473)
@@ -8,23 +8,22 @@
#include "respModule.h"
extern "C" {
- using Eigen::Map;
- using Eigen::MatrixXd;
- using Eigen::VectorXd;
+ using Eigen::MatrixXd;
+ using Eigen::VectorXd;
+ typedef Eigen::Map<VectorXd> MVec;
- using Rcpp::CharacterVector;
- using Rcpp::Environment;
- using Rcpp::IntegerVector;
- using Rcpp::Language;
- using Rcpp::List;
- using Rcpp::Named;
- using Rcpp::NumericVector;
- using Rcpp::S4;
- using Rcpp::XPtr;
- using Rcpp::as;
- using Rcpp::wrap;
+ using Rcpp::CharacterVector;
+ using Rcpp::Environment;
+ using Rcpp::IntegerVector;
+ using Rcpp::Language;
+ using Rcpp::List;
+ using Rcpp::Named;
+ using Rcpp::NumericVector;
+ using Rcpp::XPtr;
+ using Rcpp::as;
+ using Rcpp::wrap;
- using glm::glmFamily;
+ using glm::glmFamily;
using lme4Eigen::glmResp;
using lme4Eigen::lmResp;
@@ -32,6 +31,8 @@
using lme4Eigen::merPredD;
using lme4Eigen::nlsResp;
+ using std::runtime_error;
+
SEXP Eigen_SSE() {
BEGIN_RCPP;
return wrap(Eigen::SimdInstructionSetsInUse());
@@ -61,12 +62,6 @@
END_RCPP;
}
- SEXP glm_eta(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<glmResp>(ptr_)->eta());
- END_RCPP;
- }
-
SEXP glm_family(SEXP ptr_) {
BEGIN_RCPP;
return wrap(XPtr<glmResp>(ptr_)->family());
@@ -85,12 +80,6 @@
END_RCPP;
}
- SEXP glm_n(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<glmResp>(ptr_)->n());
- END_RCPP;
- }
-
SEXP glm_resDev(SEXP ptr_) {
BEGIN_RCPP;
return ::Rf_ScalarReal(XPtr<glmResp>(ptr_)->resDev());
@@ -137,7 +126,7 @@
SEXP glm_updateMu(SEXP ptr_, SEXP gamma) {
BEGIN_RCPP;
- return ::Rf_ScalarReal(XPtr<glmResp>(ptr_)->updateMu(as<Map<VectorXd> >(gamma)));
+ return ::Rf_ScalarReal(XPtr<glmResp>(ptr_)->updateMu(as<MVec>(gamma)));
END_RCPP;
}
@@ -152,27 +141,27 @@
SEXP glmFamily_link(SEXP ptr, SEXP mu) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->linkFun(as<Map<VectorXd> >(mu)));
+ return wrap(XPtr<glmFamily>(ptr)->linkFun(as<MVec>(mu)));
END_RCPP;
}
SEXP glmFamily_linkInv(SEXP ptr, SEXP eta) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->linkInv(as<Map<VectorXd> >(eta)));
+ return wrap(XPtr<glmFamily>(ptr)->linkInv(as<MVec>(eta)));
END_RCPP;
}
SEXP glmFamily_devResid(SEXP ptr, SEXP mu, SEXP weights, SEXP y) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->devResid(as<Map<VectorXd> >(mu),
- as<Map<VectorXd> >(weights),
- as<Map<VectorXd> >(y)));
+ return wrap(XPtr<glmFamily>(ptr)->devResid(as<MVec>(mu),
+ as<MVec>(weights),
+ as<MVec>(y)));
END_RCPP;
}
SEXP glmFamily_muEta(SEXP ptr, SEXP eta) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->muEta(as<Map<VectorXd> >(eta)));
+ return wrap(XPtr<glmFamily>(ptr)->muEta(as<MVec>(eta)));
END_RCPP;
}
@@ -199,20 +188,25 @@
return;
}
}
- throw std::runtime_error("step factor reduced below 0.001 without reducing pwrss");
+ throw runtime_error("step factor reduced below 0.001 without reducing pwrss");
}
+#define MAXITER 30
void pwrssUpdate(glmResp *rp, merPredD *pp, int verb, bool uOnly, double tol) {
- // FIXME: limit the number of iterations?
- do {
+ bool cvgd(false);
+ for (int it=0; it < MAXITER; ++it) {
rp->updateMu(pp->linPred(0.));
rp->updateWts();
pp->updateXwts(rp->sqrtXwt());
pp->updateDecomp();
pp->updateRes(rp->wtres());
- if ((uOnly ? pp->solveU() : pp->solve())/pwrss(rp, pp, 0.) < tol) break;
+ if ((uOnly ? pp->solveU() : pp->solve())/pwrss(rp, pp, 0.) < tol) {
+ cvgd = true;
+ break;
+ }
stepFac(rp, pp, verb);
- } while (true);
+ }
+ if (!cvgd) throw runtime_error("pwrss failed to converge in 30 iterations");
}
SEXP glmerPwrssUpdate(SEXP pp_, SEXP rp_, SEXP verb_, SEXP uOnly_, SEXP tol) {
@@ -247,7 +241,7 @@
XPtr<glmResp> rp(rp_);
XPtr<merPredD> pp(pp_);
- pp->setTheta(as<Map<VectorXd> >(theta_));
+ pp->setTheta(as<MVec>(theta_));
pp->setU0(as<VectorXd>(u0_));
pp->setBeta0(as<VectorXd>(beta0_));
pwrssUpdate(rp, pp, ::Rf_asInteger(verbose_), ::Rf_asLogical(uOnly_),
@@ -276,68 +270,24 @@
SEXP lm_setOffset(SEXP ptr_, SEXP offset) {
BEGIN_RCPP;
XPtr<lmResp>(ptr_)->setOffset(as<VectorXd>(offset));
- return offset;
END_RCPP;
}
SEXP lm_setWeights(SEXP ptr_, SEXP weights) {
BEGIN_RCPP;
XPtr<lmResp>(ptr_)->setWeights(as<VectorXd>(weights));
- return weights;
END_RCPP;
}
- SEXP lm_mu(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<lmResp>(ptr_)->mu());
- END_RCPP;
- }
-
- SEXP lm_offset(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<lmResp>(ptr_)->offset());
- END_RCPP;
- }
-
- SEXP lm_sqrtXwt(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<lmResp>(ptr_)->sqrtXwt());
- END_RCPP;
- }
-
- SEXP lm_sqrtrwt(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<lmResp>(ptr_)->sqrtrwt());
- END_RCPP;
- }
-
- SEXP lm_weights(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<lmResp>(ptr_)->weights());
- END_RCPP;
- }
-
SEXP lm_wrss(SEXP ptr_) {
BEGIN_RCPP;
return ::Rf_ScalarReal(XPtr<lmResp>(ptr_)->wrss());
END_RCPP;
}
- SEXP lm_wtres(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<lmResp>(ptr_)->wtres());
- END_RCPP;
- }
-
- SEXP lm_y(SEXP ptr_) {
- BEGIN_RCPP;
- return wrap(XPtr<lmResp>(ptr_)->y());
- END_RCPP;
- }
-
SEXP lm_updateMu(SEXP ptr_, SEXP gamma) {
BEGIN_RCPP;
- return ::Rf_ScalarReal(XPtr<lmerResp>(ptr_)->updateMu(as<Map<VectorXd> >(gamma)));
+ return ::Rf_ScalarReal(XPtr<lmerResp>(ptr_)->updateMu(as<MVec>(gamma)));
END_RCPP;
}
@@ -359,12 +309,6 @@
END_RCPP;
}
- SEXP lmer_REML(SEXP ptr_) {
- BEGIN_RCPP;
- return ::Rf_ScalarInteger(XPtr<lmerResp>(ptr_)->REML());
- END_RCPP;
- }
-
SEXP lmer_Laplace(SEXP ptr_, SEXP ldL2, SEXP ldRX2, SEXP sqrL) {
BEGIN_RCPP;
return ::Rf_ScalarReal(XPtr<lmerResp>(ptr_)->Laplace(::Rf_asReal(ldL2),
@@ -377,22 +321,13 @@
BEGIN_RCPP;
XPtr<lmerResp> rpt(rptr_);
XPtr<merPredD> ppt(pptr_);
-//Rcpp::Rcout << "Extracted pointers" << std::endl;
- ppt->setTheta(as<Map<VectorXd> >(theta_));
-//Rcpp::Rcout << "Set theta" << std::endl;
+ ppt->setTheta(as<MVec>(theta_));
ppt->updateDecomp();
-//Rcpp::Rcout << "Updated decomp" << std::endl;
rpt->updateMu(ppt->linPred(0.));
-//Rcpp::Rcout << "Updated mu for fac 0" << std::endl;
ppt->updateRes(rpt->wtres());
-//Rcpp::Rcout << "Updated residual" << std::endl;
ppt->solve();
-//Rcpp::Rcout << "Solve" << std::endl;
rpt->updateMu(ppt->linPred(1.));
-//Rcpp::Rcout << "Updated mu for fac 1" << std::endl;
- double ans(rpt->Laplace(ppt->ldL2(), ppt->ldRX2(), ppt->sqrL(1.)));
-//Rcpp::Rcout << "ans = " << ans << std::endl;
- return ::Rf_ScalarReal(ans);
+ return ::Rf_ScalarReal(rpt->Laplace(ppt->ldL2(), ppt->ldRX2(), ppt->sqrL(1.)));
END_RCPP;
}
@@ -400,38 +335,22 @@
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) {
+ SEXP Vtr, SEXP Xwts, SEXP Zt, SEXP beta0,
+ SEXP delb, SEXP delu, SEXP theta, SEXP u0) {
BEGIN_RCPP;
- S4 X(Xs);
- merPredD *ans = new merPredD(X, Lambdat, LamtUt, Lind, RZX, Ut, Utr,
- V, VtV, Vtr, Zt, beta0, delb, delu, theta, u0);
+ merPredD *ans = new merPredD(Xs, Lambdat, LamtUt, Lind, RZX, Ut, Utr, V, VtV,
+ Vtr, Xwts, Zt, beta0, delb, delu, theta, u0);
return wrap(XPtr<merPredD>(ans, true));
END_RCPP;
}
// setters
- SEXP merPredDsetBeta0(SEXP ptr, SEXP beta0) {
- BEGIN_RCPP;
- XPtr<merPredD>(ptr)->setBeta0(as<Map<VectorXd> >(beta0));
- return beta0;
- END_RCPP;
- }
-
SEXP merPredDsetTheta(SEXP ptr, SEXP theta) {
BEGIN_RCPP;
- XPtr<merPredD>(ptr)->setTheta(as<Map<VectorXd> >(theta));
+ XPtr<merPredD>(ptr)->setTheta(as<MVec>(theta));
return theta;
END_RCPP;
}
-
- SEXP merPredDsetU0(SEXP ptr, SEXP u0) {
- BEGIN_RCPP;
- XPtr<merPredD>(ptr)->setU0(as<Map<VectorXd> >(u0));
- return u0;
- END_RCPP;
- }
-
// getters
SEXP merPredDCcNumer(SEXP ptr) {
BEGIN_RCPP;
@@ -445,12 +364,6 @@
END_RCPP;
}
- SEXP merPredDLamtUt(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->LamtUt());
- END_RCPP;
- }
-
SEXP merPredDPvec(SEXP ptr) {
BEGIN_RCPP;
return wrap(XPtr<merPredD>(ptr)->Pvec());
@@ -475,67 +388,6 @@
END_RCPP;
}
- SEXP merPredDRZX(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->RZX());
- END_RCPP;
- }
-
- SEXP merPredDUt(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->Ut());
- END_RCPP;
- }
-
- SEXP merPredDUtr(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->Utr());
- END_RCPP;
- }
-
- SEXP merPredDV(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->V());
- END_RCPP;
- }
-
- SEXP merPredDVtV(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->VtV());
- END_RCPP;
- }
-
- SEXP merPredDVtr(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->Vtr());
- END_RCPP;
- }
-
-#if 0
- SEXP merPredDZt(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->Zt());
- END_RCPP;
- }
-
- SEXP merPredDbeta0(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->beta0());
- END_RCPP;
- }
-#endif
- SEXP merPredDdelb(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->delb());
- END_RCPP;
- }
-
- SEXP merPredDdelu(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->delu());
- END_RCPP;
- }
-
SEXP merPredDldL2(SEXP ptr) {
BEGIN_RCPP;
return ::Rf_ScalarReal(XPtr<merPredD>(ptr)->ldL2());
@@ -548,19 +400,6 @@
END_RCPP;
}
-#if 0
- SEXP merPredDtheta(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->theta());
- END_RCPP;
- }
-
- SEXP merPredDu0(SEXP ptr) {
- BEGIN_RCPP;
- return wrap(XPtr<merPredD>(ptr)->u0());
- END_RCPP;
- }
-#endif
SEXP merPredDunsc(SEXP ptr) {
BEGIN_RCPP;
return wrap(XPtr<merPredD>(ptr)->unsc());
@@ -617,7 +456,6 @@
END_RCPP;
}
- // methods
SEXP merPredDupdateDecomp(SEXP ptr) {
BEGIN_RCPP;
XPtr<merPredD>(ptr)->updateDecomp();
@@ -630,7 +468,6 @@
END_RCPP;
}
-
SEXP merPredDupdateLamtUt(SEXP ptr) {
BEGIN_RCPP;
XPtr<merPredD>(ptr)->updateLamtUt();
@@ -639,16 +476,18 @@
SEXP merPredDupdateRes(SEXP ptr, SEXP wtres) {
BEGIN_RCPP;
- XPtr<merPredD>(ptr)->updateRes(as<Map<VectorXd> >(wtres));
+ XPtr<merPredD>(ptr)->updateRes(as<MVec>(wtres));
END_RCPP;
}
SEXP merPredDupdateXwts(SEXP ptr, SEXP wts) {
BEGIN_RCPP;
- XPtr<merPredD>(ptr)->updateXwts(as<Map<VectorXd> >(wts));
+ XPtr<merPredD>(ptr)->updateXwts(as<MVec>(wts));
END_RCPP;
}
+ // 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) {
BEGIN_RCPP;
@@ -671,7 +510,7 @@
SEXP nls_updateMu(SEXP ptr_, SEXP gamma) {
BEGIN_RCPP;
- return ::Rf_ScalarReal(XPtr<nlsResp>(ptr_)->updateMu(as<Map<VectorXd> >(gamma)));
+ return ::Rf_ScalarReal(XPtr<nlsResp>(ptr_)->updateMu(as<MVec>(gamma)));
END_RCPP;
}
@@ -685,29 +524,27 @@
CALLDEF(Eigen_SSE, 0),
- CALLDEF(glm_Create, 10), // generate external pointer
+ CALLDEF(glm_Create, 10), // generate external pointer
- CALLDEF(glm_setN, 2), // setters
+ CALLDEF(glm_setN, 2), // setters
- CALLDEF(glm_devResid, 1), // getters
- CALLDEF(glm_eta, 1),
+ CALLDEF(glm_devResid, 1), // getters
CALLDEF(glm_family, 1),
CALLDEF(glm_link, 1),
CALLDEF(glm_muEta, 1),
- CALLDEF(glm_n, 1),
CALLDEF(glm_resDev, 1),
CALLDEF(glm_sqrtWrkWt, 1),
CALLDEF(glm_variance, 1),
CALLDEF(glm_wrkResids, 1),
CALLDEF(glm_wrkResp, 1),
- CALLDEF(glm_Laplace, 4), // methods
+ CALLDEF(glm_Laplace, 4), // methods
CALLDEF(glm_updateMu, 2),
CALLDEF(glm_updateWts, 1),
CALLDEF(glmFamily_Create, 1), // generate external pointer
- CALLDEF(glmFamily_link, 2), // methods
+ CALLDEF(glmFamily_link, 2), // methods
CALLDEF(glmFamily_linkInv, 2),
CALLDEF(glmFamily_devResid, 4),
CALLDEF(glmFamily_muEta, 2),
@@ -719,51 +556,32 @@
CALLDEF(isNullExtPtr, 1),
- CALLDEF(lm_Create, 7), // generate external pointer
+ CALLDEF(lm_Create, 7), // generate external pointer
- CALLDEF(lm_setOffset, 2), // setters
+ CALLDEF(lm_setOffset, 2), // setters
CALLDEF(lm_setWeights, 2),
- CALLDEF(lm_mu, 1), // getters
- CALLDEF(lm_offset, 1),
- CALLDEF(lm_sqrtXwt, 1),
- CALLDEF(lm_sqrtrwt, 1),
- CALLDEF(lm_weights, 1),
- CALLDEF(lm_wrss, 1),
- CALLDEF(lm_wtres, 1),
- CALLDEF(lm_y, 1),
+ CALLDEF(lm_wrss, 1), // getter
- CALLDEF(lm_updateMu, 2), // method
+ CALLDEF(lm_updateMu, 2), // method
- CALLDEF(lmer_Create, 7), // generate external pointer
+ CALLDEF(lmer_Create, 7), // generate external pointer
- CALLDEF(lmer_setREML, 2), // setter
+ CALLDEF(lmer_setREML, 2), // setter
- CALLDEF(lmer_REML, 1), // getter
-
- CALLDEF(lmer_Deviance, 3), // method
+ CALLDEF(lmer_Deviance, 3), // methods
CALLDEF(lmer_Laplace, 4),
- CALLDEF(merPredDCreate, 16), // generate external pointer
+ CALLDEF(merPredDCreate, 17), // generate external pointer
CALLDEF(merPredDsetTheta, 2), // setters
- CALLDEF(merPredDsetBeta0, 2),
- CALLDEF(merPredDsetU0, 2),
- CALLDEF(merPredDCcNumer, 1), // getters
+ CALLDEF(merPredDCcNumer, 1), // getters
CALLDEF(merPredDL, 1),
- CALLDEF(merPredDLamtUt, 1),
CALLDEF(merPredDPvec, 1),
CALLDEF(merPredDRX, 1),
+ CALLDEF(merPredDRXi, 1),
CALLDEF(merPredDRXdiag, 1),
- CALLDEF(merPredDRZX, 1),
- CALLDEF(merPredDUt, 1),
- CALLDEF(merPredDUtr, 1),
- CALLDEF(merPredDV, 1),
- CALLDEF(merPredDVtV, 1),
- CALLDEF(merPredDVtr, 1),
- CALLDEF(merPredDdelb, 1),
- CALLDEF(merPredDdelu, 1),
CALLDEF(merPredDldL2, 1),
CALLDEF(merPredDldRX2, 1),
CALLDEF(merPredDunsc, 1),
@@ -782,9 +600,9 @@
CALLDEF(merPredDupdateRes, 2),
CALLDEF(merPredDupdateXwts, 2),
- CALLDEF(nls_Create, 10), // generate external pointer
+ CALLDEF(nls_Create, 10), // generate external pointer
- CALLDEF(nls_Laplace, 4), // methods
+ CALLDEF(nls_Laplace, 4), // methods
CALLDEF(nls_updateMu, 2),
{NULL, NULL, 0}
Modified: pkg/lme4Eigen/src/predModule.cpp
===================================================================
--- pkg/lme4Eigen/src/predModule.cpp 2011-12-06 23:25:08 UTC (rev 1472)
+++ pkg/lme4Eigen/src/predModule.cpp 2011-12-07 23:25:59 UTC (rev 1473)
@@ -8,47 +8,48 @@
#include "predModule.h"
namespace lme4Eigen {
- using std::invalid_argument;
+ using Rcpp::as;
- typedef Map<MatrixXd> MMat;
- typedef Map<VectorXd> MVec;
- typedef Map<VectorXi> MiVec;
+ using std::invalid_argument;
+ using std::runtime_error;
- merPredD::merPredD(S4 X, SEXP Lambdat, SEXP LamtUt, SEXP Lind,
+ typedef Eigen::Map<MatrixXd> MMat;
+ typedef Eigen::Map<VectorXd> MVec;
+ typedef Eigen::Map<VectorXi> MiVec;
+
+ merPredD::merPredD(SEXP 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)),
- d_Lind( as<MiVec>(Lind)),
- d_n( d_X.rows()),
- d_nnz( -1),
- d_p( d_X.cols()),
- d_q( d_Zt.rows()),
- d_Lambdat( as<MSpMatrixd>(Lambdat)),
+ SEXP Vtr, SEXP Xwts, SEXP Zt, SEXP beta0,
+ SEXP delb, SEXP delu, SEXP theta, SEXP u0)
+ : d_X( as<MMat>(X)),
d_RZX( as<MMat>(RZX)),
d_V( as<MMat>(V)),
d_VtV( as<MMat>(VtV)),
+ d_Zt( as<MSpMatrixd>(Zt)),
+ d_Ut( as<MSpMatrixd>(Ut)),
+ d_LamtUt( as<MSpMatrixd>(LamtUt)),
+ d_Lambdat( as<MSpMatrixd>(Lambdat)),
+ d_theta( as<MVec>(theta)),
d_Vtr( as<MVec>(Vtr)),
d_Utr( as<MVec>(Utr)),
+ d_Xwts( as<MVec>(Xwts)),
d_beta0( as<MVec>(beta0)),
d_delb( as<MVec>(delb)),
d_delu( as<MVec>(delu)),
d_u0( as<MVec>(u0)),
- d_Ut( as<MSpMatrixd>(Ut)),
- d_LamtUt( as<MSpMatrixd>(LamtUt)),
+ d_Lind( as<MiVec>(Lind)),
+ d_N( d_X.rows()),
+ d_p( d_X.cols()),
+ d_q( d_Zt.rows()),
d_RX( d_p)
// d_LamtUtRestructure(false)
{ // Check consistency of dimensions
- if (d_n != d_Zt.cols())
+ if (d_N != d_Zt.cols())
throw invalid_argument("Z dimension mismatch");
if (d_Lind.size() != d_Lambdat.nonZeros())
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_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
@@ -148,15 +149,12 @@
}
void merPredD::updateXwts(const VectorXd& sqrtXwt) {
-// Rcpp::Rcout << "X[" << d_X.rows() << ", " << d_X.cols() << "]"
-// << "V[" << d_V.rows() << ", " << d_V.cols() << "]"
-// << ", sqrtXwt.size() = " << sqrtXwt.size() << std::endl;
-
- if (d_X.rows() != sqrtXwt.size())
+ if (d_Xwts.size() != sqrtXwt.size())
throw invalid_argument("updateXwts: dimension mismatch");
+ std::copy(sqrtXwt.data(), sqrtXwt.data() + sqrtXwt.size(), d_Xwts.data());
if (sqrtXwt.size() == d_V.rows()) {
d_V = sqrtXwt.asDiagonal() * d_X;
- for (int j = 0; j < d_n; ++j)
+ 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];
@@ -219,10 +217,11 @@
std::copy(newU0.data(), newU0.data() + d_q, d_u0.data());
}
- IntegerVector merPredD::Pvec() const {
- const cholmod_factor* cf = d_L.factor();
- int* ppt = (int*)cf->Perm;
- return IntegerVector(ppt, ppt + cf->n);
+ VectorXi merPredD::Pvec() const {
+ int* ppt((int*)d_L.factor()->Perm);
+ VectorXi ans(d_q);
+ std::copy(ppt, ppt + d_q, ans.data());
+ return ans;
}
MatrixXd merPredD::RX() const {
@@ -235,7 +234,8 @@
MatrixXd merPredD::unsc() const {
return MatrixXd(MatrixXd(d_p, d_p).setZero().
- selfadjointView<Lower>().rankUpdate(RXi()));
+ selfadjointView<Eigen::Lower>().
+ rankUpdate(RXi()));
}
VectorXd merPredD::RXdiag() const {
Modified: pkg/lme4Eigen/src/predModule.h
===================================================================
--- pkg/lme4Eigen/src/predModule.h 2011-12-06 23:25:08 UTC (rev 1472)
+++ pkg/lme4Eigen/src/predModule.h 2011-12-07 23:25:59 UTC (rev 1473)
@@ -12,24 +12,7 @@
#include <RcppEigen.h>
namespace lme4Eigen {
- using Eigen::CholmodDecomposition;
- using Eigen::LLT;
- using Eigen::Lower;
- using Eigen::Map;
- using Eigen::MatrixXd;
- using Eigen::SelfAdjointView;
- using Eigen::SparseMatrix;
- using Eigen::VectorXd;
- using Eigen::VectorXi;
-
- using Rcpp::IntegerVector;
- using Rcpp::List;
- using Rcpp::NumericVector;
- using Rcpp::S4;
- using Rcpp::as;
-
- using std::runtime_error;
-
+#if 0
class dgCMatrix : public SparseMatrix<double> { // sparse Matrix, copies contents of R object
public:
dgCMatrix(const S4& xp)
@@ -73,7 +56,6 @@
modelMatrix(xp), d_xp(xp) {}
};
-#if 0
class dsparseModelMatrix : public MappedSparseMatrix<double>, public modelMatrix {
protected:
S4 d_xp;
@@ -117,84 +99,83 @@
}
#endif
+ using Eigen::LLT;
+ using Eigen::MatrixXd;
+ using Eigen::VectorXd;
+ using Eigen::VectorXi;
+
class merPredD {
public:
- typedef ddenseModelMatrix XType;
- typedef XType::Scalar Scalar;
- typedef XType::Index Index;
- typedef CholmodDecomposition<SparseMatrix<double> > ChmDecomp;
- typedef SparseMatrix<double> SpMatrixd;
- typedef Eigen::MappedSparseMatrix<double> MSpMatrixd;
+ typedef Eigen::Map<MatrixXd> MMap;
+ typedef Eigen::Map<VectorXd> MVec;
+ typedef Eigen::Map<VectorXi> MiVec;
+ typedef MatrixXd::Scalar Scalar;
+ typedef MatrixXd::Index Index;
+ typedef Eigen::SparseMatrix<double> SpMatrixd;
+ typedef Eigen::CholmodDecomposition<SpMatrixd> ChmDecomp;
+ typedef Eigen::MappedSparseMatrix<double> MSpMatrixd;
protected:
- XType d_X;
- MSpMatrixd d_Zt;
- Map<VectorXd> d_theta;
- Map<VectorXi> d_Lind;
- Index d_n, d_nnz, d_p, d_q;
- MSpMatrixd d_Lambdat;
- Scalar d_CcNumer, d_ldL2, d_ldRX2;
- 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;
+ MMap d_X, d_RZX, d_V, d_VtV;
+ MSpMatrixd d_Zt, d_Ut, d_LamtUt, d_Lambdat;
+ MVec d_theta, d_Vtr, d_Utr, d_Xwts, d_beta0, d_delb, d_delu, d_u0;
+ MiVec d_Lind;
+ Index d_N, d_p, d_q;
+ Scalar d_CcNumer, d_ldL2, d_ldRX2;
+ ChmDecomp d_L;
+ LLT<MatrixXd> d_RX;
public:
- merPredD(S4, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP,
- SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
+ merPredD(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP,
+ SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
- IntegerVector Pvec() const;
+ VectorXi Pvec() const;
- MatrixXd RX() const;
- MatrixXd RXi() const;
- MatrixXd unsc() const;
+ MatrixXd RX() const;
+ MatrixXd RXi() const;
+ MatrixXd unsc() const;
- VectorXd RXdiag() const;
- VectorXd b(const Scalar& f) const;
- VectorXd beta(const Scalar& f) const;
- VectorXd linPred(const Scalar& f) const;
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/lme4 -r 1473
More information about the Lme4-commits
mailing list