[Rcpp-commits] r3275 - in pkg/RcppEigen: inst/include vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 3 23:17:20 CET 2011
Author: dmbates
Date: 2011-11-03 23:17:20 +0100 (Thu, 03 Nov 2011)
New Revision: 3275
Modified:
pkg/RcppEigen/inst/include/RcppEigenWrap.h
pkg/RcppEigen/vignettes/RcppEigen-Intro.Rnw
Log:
Got the wrap method for CholmodDecomposition working and incorporated it in the vignette
Modified: pkg/RcppEigen/inst/include/RcppEigenWrap.h
===================================================================
--- pkg/RcppEigen/inst/include/RcppEigenWrap.h 2011-11-03 20:52:17 UTC (rev 3274)
+++ pkg/RcppEigen/inst/include/RcppEigenWrap.h 2011-11-03 22:17:20 UTC (rev 3275)
@@ -24,39 +24,48 @@
namespace Rcpp{
-#if 0 // Have to implement this in some other way
- template<>
- SEXP wrap(const Eigen::CholmodDecomposition<Eigen::SparseMatrix<double> >& obj) {
- const cholmod_factor* f = obj.factor();
- if (f->minor < f->n)
- throw std::runtime_error("CHOLMOD factorization was unsuccessful");
+ namespace RcppEigen{
- S4 ans(std::string(f->is_super ? "dCHMsuper" : "dCHMsimpl"));
- ans.slot("Dim") = Dimension(f->n, f->n);
- ans.slot("perm") = IntegerVector((int*)f->Perm, (int*)f->Perm + f->n);
- ans.slot("colcount") = IntegerVector((int*)f->ColCount, (int*)f->ColCount + f->n);
- IntegerVector tt(f->is_super ? 6 : 4);
- tt[0] = f->ordering; tt[1] = f->is_ll;
- tt[2] = f->is_super; tt[3] = f->is_monotonic;
- ans.slot("type") = tt;
- if (f->is_super) {
- tt[4] = f->maxcsize; tt[5] = f->maxesize;
- ans.slot("super") = IntegerVector((int*)f->super, ((int*)f->super) + f->nsuper + 1);
- ans.slot("pi") = IntegerVector((int*)f->pi, ((int*)f->pi) + f->nsuper + 1);
- ans.slot("px") = IntegerVector((int*)f->px, ((int*)f->px) + f->nsuper + 1);
- ans.slot("s") = IntegerVector((int*)f->s, ((int*)f->s) + f->ssize);
- ans.slot("x") = NumericVector((double*)f->x, ((double*)f->x) + f->xsize);
- } else {
- ans.slot("i") = IntegerVector((int*)f->i, ((int*)f->i) + f->nzmax);
- ans.slot("p") = IntegerVector((int*)f->p, ((int*)f->p) + f->n + 1);
- ans.slot("x") = NumericVector((double*)f->x, ((double*)f->x) + f->nzmax);
- ans.slot("nz") = IntegerVector((int*)f->nz, ((int*)f->nz) + f->n);
- ans.slot("nxt") = IntegerVector((int*)f->next, ((int*)f->next) + f->n + 2);
- ans.slot("prv") = IntegerVector((int*)f->prev, ((int*)f->prev) + f->n + 2);
+ template<typename T>
+ SEXP Eigen_cholmod_wrap(const Eigen::CholmodDecomposition<Eigen::SparseMatrix<T> >& obj) {
+ const cholmod_factor* f = obj.factor();
+ if (f->minor < f->n)
+ throw std::runtime_error("CHOLMOD factorization was unsuccessful");
+
+//FIXME: Should extend this selection according to T
+ S4 ans(std::string(f->is_super ? "dCHMsuper" : "dCHMsimpl"));
+ ans.slot("Dim") = Dimension(f->n, f->n);
+ ans.slot("perm") = ::Rcpp::wrap((int*)f->Perm, (int*)f->Perm + f->n);
+ ans.slot("colcount") = ::Rcpp::wrap((int*)f->ColCount, (int*)f->ColCount + f->n);
+ IntegerVector tt(f->is_super ? 6 : 4);
+ tt[0] = f->ordering; tt[1] = f->is_ll;
+ tt[2] = f->is_super; tt[3] = f->is_monotonic;
+ ans.slot("type") = tt;
+ if (f->is_super) {
+ tt[4] = f->maxcsize; tt[5] = f->maxesize;
+ ans.slot("super") = ::Rcpp::wrap((int*)f->super, ((int*)f->super) + f->nsuper + 1);
+ ans.slot("pi") = ::Rcpp::wrap((int*)f->pi, ((int*)f->pi) + f->nsuper + 1);
+ ans.slot("px") = ::Rcpp::wrap((int*)f->px, ((int*)f->px) + f->nsuper + 1);
+ ans.slot("s") = ::Rcpp::wrap((int*)f->s, ((int*)f->s) + f->ssize);
+ ans.slot("x") = ::Rcpp::wrap((T*)f->x, ((T*)f->x) + f->xsize);
+ } else {
+ ans.slot("i") = ::Rcpp::wrap((int*)f->i, ((int*)f->i) + f->nzmax);
+ ans.slot("p") = ::Rcpp::wrap((int*)f->p, ((int*)f->p) + f->n + 1);
+ ans.slot("x") = ::Rcpp::wrap((T*)f->x, ((T*)f->x) + f->nzmax);
+ ans.slot("nz") = ::Rcpp::wrap((int*)f->nz, ((int*)f->nz) + f->n);
+ ans.slot("nxt") = ::Rcpp::wrap((int*)f->next, ((int*)f->next) + f->n + 2);
+ ans.slot("prv") = ::Rcpp::wrap((int*)f->prev, ((int*)f->prev) + f->n + 2);
+ }
+ return ::Rcpp::wrap(ans);
}
- return wrap(ans);
- }
-#endif
+
+ } /* namespace RcppEigen */
+
+ template<typename T>
+ SEXP wrap(const Eigen::CholmodDecomposition<Eigen::SparseMatrix<T> >& obj) {
+ return RcppEigen::Eigen_cholmod_wrap(obj);
+ }
+
namespace RcppEigen{
// helper trait to identify if T is a plain object type
@@ -122,15 +131,7 @@
typename is_plain<T>::type()
) ;
}
-
-#if 0 // I don't think this is used now
- template <typename T>
- SEXP Eigen_wrap( const T& object, const ::Rcpp::Dimension& dim){
- ::Rcpp::RObject x = ::Rcpp::wrap(object.data(), object.data() + object.size());
- x.attr( "dim" ) = dim ;
- return x;
- }
-#endif
+
} /* namespace RcppEigen */
Modified: pkg/RcppEigen/vignettes/RcppEigen-Intro.Rnw
===================================================================
--- pkg/RcppEigen/vignettes/RcppEigen-Intro.Rnw 2011-11-03 20:52:17 UTC (rev 3274)
+++ pkg/RcppEigen/vignettes/RcppEigen-Intro.Rnw 2011-11-03 22:17:20 UTC (rev 3275)
@@ -1426,12 +1426,13 @@
const SpChol Ch(At * At.adjoint());
if (Ch.info() != Eigen::Success)
return R_NilValue;
-const CholMD Cm(At);
-if (Cm.info() != Eigen::Success)
+const CholMD L(At);
+if (L.info() != Eigen::Success)
return R_NilValue;
-return List::create(_["betahatS"] = Ch.solve(Aty),
- _["betahatC"] = Cm.solve(Aty),
- _["perm"] = Ch.permutationP().indices());
+return List::create(_["L"] = wrap(L),
+ _["betahatS"] = Ch.solve(Aty),
+ _["betahatC"] = L.solve(Aty),
+ _["perm"] = Ch.permutationP().indices());
'
@
\begin{figure}[htb]
@@ -1451,12 +1452,13 @@
\hlstd{}\hlkwb{const\ }\hlstd{SpChol}\hlstd{\ \ \ \ \ }\hlstd{}\hlkwd{Ch}\hlstd{}\hlopt{(}\hlstd{At\ }\hlopt{{*}\ }\hlstd{At}\hlopt{.}\hlstd{}\hlkwd{adjoint}\hlstd{}\hlopt{());}\hspace*{\fill}\\
\hlstd{}\hlkwa{if\ }\hlstd{}\hlopt{(}\hlstd{Ch}\hlopt{.}\hlstd{}\hlkwd{info}\hlstd{}\hlopt{()\ !=\ }\hlstd{Eigen}\hlopt{::}\hlstd{Success}\hlopt{)}\hspace*{\fill}\\
\hlstd{}\hlstd{\ \ \ }\hlstd{}\hlkwa{return\ }\hlstd{R\textunderscore NilValue}\hlopt{;}\hspace*{\fill}\\
- \hlstd{}\hlkwb{const\ }\hlstd{CholMD}\hlstd{\ \ \ \ \ }\hlstd{}\hlkwd{Cm}\hlstd{}\hlopt{(}\hlstd{At}\hlopt{);}\hspace*{\fill}\\
- \hlstd{}\hlkwa{if\ }\hlstd{}\hlopt{(}\hlstd{Cm}\hlopt{.}\hlstd{}\hlkwd{info}\hlstd{}\hlopt{()\ !=\ }\hlstd{Eigen}\hlopt{::}\hlstd{Success}\hlopt{)}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{CholMD}\hlstd{\ \ \ \ \ \ }\hlstd{}\hlkwd{L}\hlstd{}\hlopt{(}\hlstd{At}\hlopt{);}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{if\ }\hlstd{}\hlopt{(}\hlstd{L}\hlopt{.}\hlstd{}\hlkwd{info}\hlstd{}\hlopt{()\ !=\ }\hlstd{Eigen}\hlopt{::}\hlstd{Success}\hlopt{)}\hspace*{\fill}\\
\hlstd{}\hlstd{\ \ \ }\hlstd{}\hlkwa{return\ }\hlstd{R\textunderscore NilValue}\hlopt{;}\hspace*{\fill}\\
- \hlstd{}\hlkwa{return\ }\hlstd{List}\hlopt{::}\hlstd{}\hlkwd{create}\hlstd{}\hlopt{(}\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"betahatS"}\hlstd{}\hlopt{{]}\ =\ }\hlstd{Ch}\hlopt{.}\hlstd{}\hlkwd{solve}\hlstd{}\hlopt{(}\hlstd{Aty}\hlopt{),}\hspace*{\fill}\\
- \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"betahatC"}\hlstd{}\hlopt{{]}\ =\ }\hlstd{Cm}\hlopt{.}\hlstd{}\hlkwd{solve}\hlstd{}\hlopt{(}\hlstd{Aty}\hlopt{),}\hspace*{\fill}\\
- \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"perm"}\hlstd{}\hlopt{{]}}\hlstd{\ \ \ \ }\hlopt{=\ }\hlstd{Ch}\hlopt{.}\hlstd{}\hlkwd{permutationP}\hlstd{}\hlopt{().}\hlstd{}\hlkwd{indices}\hlstd{}\hlopt{());}\hlstd{}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{return\ }\hlstd{List}\hlopt{::}\hlstd{}\hlkwd{create}\hlstd{}\hlopt{(}\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"L"}\hlstd{}\hlopt{{]}}\hlstd{\ \ \ \ \ \ \ \ }\hlopt{=\ }\hlstd{}\hlkwd{wrap}\hlstd{}\hlopt{(}\hlstd{L}\hlopt{),}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"betahatS"}\hlstd{}\hlopt{{]}\ =\ }\hlstd{Ch}\hlopt{.}\hlstd{}\hlkwd{solve}\hlstd{}\hlopt{(}\hlstd{Aty}\hlopt{),}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"betahatC"}\hlstd{}\hlopt{{]}\ =\ }\hlstd{L}\hlopt{.}\hlstd{}\hlkwd{solve}\hlstd{}\hlopt{(}\hlstd{Aty}\hlopt{),}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"perm"}\hlstd{}\hlopt{{]}}\hlstd{\ \ \ \ \ }\hlopt{=\ }\hlstd{Ch}\hlopt{.}\hlstd{}\hlkwd{permutationP}\hlstd{}\hlopt{().}\hlstd{}\hlkwd{indices}\hlstd{}\hlopt{());}\hlstd{}\hspace*{\fill}\\
\mbox{}
\normalfont
\normalsize
@@ -1471,6 +1473,8 @@
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'
all(rr$perm == Ch at perm) # fill-reducing permutations are different
@
More information about the Rcpp-commits
mailing list