[Rcpp-commits] r3348 - pkg/RcppEigen/vignettes

[noreply at r-forge.r-project.org]

Author: dmbates
Date: 2011-11-14 20:30:12 +0100 (Mon, 14 Nov 2011)
New Revision: 3348

Modified:
   pkg/RcppEigen/vignettes/RcppEigen-intro-jss.tex
   pkg/RcppEigen/vignettes/code.R
Log:
Commit changes prior to moving back into inst/doc


Modified: pkg/RcppEigen/vignettes/RcppEigen-intro-jss.tex
===================================================================
--- pkg/RcppEigen/vignettes/RcppEigen-intro-jss.tex	2011-11-14 19:28:06 UTC (rev 3347)
+++ pkg/RcppEigen/vignettes/RcppEigen-intro-jss.tex	2011-11-14 19:30:12 UTC (rev 3348)
@@ -1,5 +1,5 @@
 
-\documentclass[shortnames,article]{jss}
+\documentclass[shortnames,article,nojss]{jss}
 %\VignetteIndexEntry{RcppEigen-intro}
 %\VignetteKeywords{linear algebra, template programming, C++, R, Rcpp}
 %\VignettePackage{RcppEigen}

Modified: pkg/RcppEigen/vignettes/code.R
===================================================================
--- pkg/RcppEigen/vignettes/code.R	2011-11-14 19:28:06 UTC (rev 3347)
+++ pkg/RcppEigen/vignettes/code.R	2011-11-14 19:30:12 UTC (rev 3348)
@@ -102,8 +102,7 @@
 
 cholDetCpp <- '
 const MatrixXd      ata(AtA(as<MapMatd>(AA)));
-const MatrixXd     Lmat(ata.llt().matrixL());
-const double       detL(Lmat.diagonal().prod());
+const double       detL(MatrixXd(ata.llt().matrixL()).diagonal().prod());
 const VectorXd     Dvec(ata.ldlt().vectorD());
 return List::create(Named("d1") = detL * detL,
                     Named("d2") = Dvec.prod(),
@@ -176,13 +175,10 @@
 all.equal(unname(residuals(fm1)), residuals(fmSVD))
 
 
-print(summary(fmVLV <- fastLm(y ~ f1 * f2, dd, method=5L)), signif.stars=FALSE)
+fmVLV <- fastLm(y ~ f1 * f2, dd, method=5L)
 all.equal(coef(fmSVD), coef(fmVLV))
-all.equal(unname(fitted(fm1)), fitted(fmSVD))
-all.equal(unname(residuals(fm1)), residuals(fmSVD))
 
 
-
 ## section 5
 
 badtransCpp <- '
@@ -199,20 +195,18 @@
 
 ## section 6
 sparseProdCpp <- '
-using Eigen::Map;
 using Eigen::MappedSparseMatrix;
 using Eigen::SparseMatrix;
-using Eigen::VectorXd;
 
 const MappedSparseMatrix<double>  A(as<MappedSparseMatrix<double> >(AA));
-const Map<VectorXd>               y(as<Map<VectorXd> >(yy));
+const MapVecd                     y(as<MapVecd>(yy));
 const SparseMatrix<double>       At(A.adjoint());
 return List::create(Named("At")  = At,
                     Named("Aty") = At * y);
 '
 
 sparse1 <- cxxfunction(signature(AA = "dgCMatrix", yy = "numeric"),
-                       sparseProdCpp, "RcppEigen")
+                       sparseProdCpp, "RcppEigen", incl)
 data(KNex, package="Matrix")
 rr <- sparse1(KNex$mm, KNex$y)
 stopifnot(all.equal(rr$At, t(KNex$mm)),
@@ -220,35 +214,29 @@
 
 
 sparseLSCpp <- '
-using   Eigen::Lower;
-using   Eigen::VectorXd;
-typedef Eigen::Map<VectorXd>               MapVec;
 typedef Eigen::MappedSparseMatrix<double>  MSpMat;
 typedef Eigen::SparseMatrix<double>         SpMat;
 typedef Eigen::SimplicialLDLt<SpMat>       SpChol;
 typedef Eigen::CholmodDecomposition<SpMat> CholMD;
 
 const SpMat      At(as<MSpMat>(AA).adjoint());
-const VectorXd  Aty(At * as<MapVec>(yy));
+const VectorXd  Aty(At * as<MapVecd>(yy));
 const SpChol     Ch(At * At.adjoint());
-if (Ch.info() != Eigen::Success)
-   return R_NilValue;
+if (Ch.info() != Eigen::Success) return R_NilValue;
 const CholMD      L(At);
-if (L.info() != Eigen::Success)
-   return R_NilValue;
+if (L.info() != Eigen::Success) return R_NilValue;
 return List::create(Named("L")        = wrap(L),
                     Named("betahatS") = Ch.solve(Aty),
                     Named("betahatC") = L.solve(Aty),
                     Named("perm")     = Ch.permutationP().indices());
 '
 
-
 sparse2 <- cxxfunction(signature(AA = "dgCMatrix", yy = "numeric"),
-                       sparseLSCpp, "RcppEigen")
+                       sparseLSCpp, "RcppEigen", incl)
 str(rr <-  sparse2(KNex$mm, KNex$y))
 res <- as.vector(solve(Ch <- Cholesky(crossprod(KNex$mm)),
                        crossprod(KNex$mm, KNex$y)))
 stopifnot(all.equal(rr$betahatS, res), all.equal(rr$betahatC, res))
-all.equal(rr$L, Ch)   # not sure yet why these are different.  The one from Eigen is smaller
-## It's because the Ch was created from mm'mm and L was created directly from mm'
+                                        # factors are different sizes
+c(nnzL=length(rr$L at x), nnzCh=length(Ch at x))
 all(rr$perm == Ch at perm) # fill-reducing permutations are different