<div dir="ltr">Hi JJ,<div><br></div><div>Thank you, that worked wonderfully. </div><div><br></div><div>Thanks to your help, I've been able to add a simple way of bootstrapping with some success, that seems about 25x faster than direct lm() calls using 4 threads.</div><div><br></div><div>I've put it online in case anyone else runs across this and is interested:</div><div><br></div><div><div><a href="https://gist.github.com/JWiley/d9cba55603471f75d438" target="_blank">https://gist.github.com/JWiley/d9cba55603471f75d438</a><br></div><div><a href="https://gist.github.com/JWiley/8fe93ae105b1fb244f93" target="_blank">https://gist.github.com/JWiley/8fe93ae105b1fb244f93</a><br></div></div><div><br></div><div>There are some residual issues (the results are not an identical match to lm(), and occasionally it can make R crash, so I am probably overwriting something in memory), but once I figure things out, I will fix the gist files too.<br></div><div><br></div><div>Thanks again for your help and the great package!</div><div><br></div><div>Cheers,</div><div><br></div><div>Josh</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Apr 14, 2015 at 7:54 PM, JJ Allaire <span dir="ltr"><<a href="mailto:jj.allaire@gmail.com" target="_blank">jj.allaire@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">The initialization of your const fields needs to happen within the<br>
initializer list (rather than inline where the fields are declared).<br>
If you substitute this code for the code in your example everything<br>
will compile:<br>
<br>
// QR Decomposition<br>
const CPivQR PQR;<br>
const Permutation Pmat;<br>
const int r;<br>
<span class=""><br>
// output<br>
RMatrix<double> betamat;<br>
Eigen::VectorXd betahat;<br>
<br>
// initialize with input and output<br>
CVLM(const Eigen::VectorXd y, const Eigen::MatrixXd X,<br>
Rcpp::NumericMatrix betamat)<br>
</span> : y(y), X(X), PQR(X), Pmat(PQR.colsPermutation()), r(PQR.rank()),<br>
betamat(betamat) {}<br>
<div><div class="h5"><br>
<br>
<br>
On Tue, Apr 14, 2015 at 4:41 AM, Joshua Wiley <<a href="mailto:jwiley.psych@gmail.com">jwiley.psych@gmail.com</a>> wrote:<br>
> Hi,<br>
><br>
> I'm interested in combining faster linear model estimation via QR<br>
> decompositions from RcppEigen with JJs new(er) RcppParallel package to do<br>
> cross validation and bootstrapping on linear models.<br>
><br>
> As a first step, I've been trying to merge the JSS paper on RcppEigen (<br>
> <a href="http://www.jstatsoft.org/v52/i05/paper" target="_blank">http://www.jstatsoft.org/v52/i05/paper</a> ) with about the documentation for<br>
> RcppParallel ( <a href="http://rcppcore.github.io/RcppParallel/#examples" target="_blank">http://rcppcore.github.io/RcppParallel/#examples</a> )<br>
><br>
> ignoring bootstrapping and cross validation for the time being and just<br>
> getting a linear model to run in the parallel framework. I've toyed with<br>
> each separately successfuly, and think I am somewhat close to getting them<br>
> together, but I am running into an error:<br>
><br>
> partest.cpp:28:20: error: 'X' is not a type<br>
> partest.cpp:29:26: error: 'PQR' is not a type<br>
> partest.cpp:29:29: error: expected ',' or '...' before '.' token<br>
> partest.cpp:30:15: error: 'PQR' is not a type<br>
> partest.cpp:30:18: error: expected ',' or '...' before '.' token<br>
> partest.cpp: In member function 'virtual void CVLM::operator()(std::size_t,<br>
> std::size_t)':<br>
> partest.cpp:44:28: warning: comparison between signed and unsigned integer<br>
> expressions [-Wsign-compare]<br>
> make: *** [partest.o] Error 1<br>
><br>
> My code is below (run via sourceCpp() ), and I can see exactly where the<br>
> first error occurs, but I don't understand why this is OK when run outside<br>
> of RcppParallel, but so very wrong in conjunction with RcppParallel.<br>
><br>
> Any pointers on either what's wrong or where I should be reading would be<br>
> deeply appreciated.<br>
><br>
> Josh (code below)<br>
><br>
><br>
> // [[Rcpp::depends(RcppParallel)]]<br>
> // [[Rcpp::depends(RcppEigen)]]<br>
> #include <Rcpp.h><br>
> #include <RcppEigen.h><br>
> #include <RcppParallel.h><br>
><br>
> using namespace std;<br>
> using namespace Rcpp;<br>
> using namespace RcppParallel;<br>
><br>
> using Eigen::Lower;<br>
> using Eigen::Map;<br>
> using Eigen::MatrixXd;<br>
> using Eigen::Upper;<br>
> using Eigen::VectorXd;<br>
> typedef Map<MatrixXd> MapMatd;<br>
> typedef Map<VectorXd> MapVecd;<br>
> typedef Eigen::ColPivHouseholderQR<MatrixXd> CPivQR;<br>
> typedef CPivQR::PermutationType Permutation;<br>
><br>
> struct CVLM : public Worker<br>
> {<br>
> // design matrix and outcome<br>
> const Eigen::VectorXd y;<br>
> const Eigen::MatrixXd X;<br>
><br>
> // QR Decomposition<br>
> const CPivQR PQR(X); // ERROR HAPPENS HERE<br>
> const Permutation Pmat(PQR.colsPermutation());<br>
> const int r(PQR.rank());<br>
><br>
><br>
> // output<br>
> RMatrix<double> betamat;<br>
> Eigen::VectorXd betahat;<br>
><br>
> // initialize with input and output<br>
> CVLM(const Eigen::VectorXd y, const Eigen::MatrixXd X, Rcpp::NumericMatrix<br>
> betamat)<br>
> : y(y), X(X), betamat(betamat) {}<br>
><br>
> void operator()(std::size_t begin, std::size_t end) {<br>
><br>
> //betahat = PQR.solve(y);<br>
> for(int j = begin; j < end; j++) {<br>
> betamat(j, j) = X(j, j) + y[j];<br>
> }<br>
> }<br>
> };<br>
><br>
> // [[Rcpp::export]]<br>
> NumericMatrix parallelFit(Eigen::VectorXd dv, Eigen::MatrixXd x) {<br>
><br>
> // allocate the output matrix<br>
> NumericMatrix betamat(x.rows(), x.cols());<br>
><br>
> // pass input and output<br>
> CVLM cvLM(dv, x, betamat);<br>
><br>
> // parallelFor to do it<br>
> parallelFor(0, x.rows(), cvLM);<br>
><br>
> return betamat;<br>
> }<br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
> --<br>
> Joshua F. Wiley<br>
> <a href="http://joshuawiley.com/" target="_blank">http://joshuawiley.com/</a><br>
> ---<br>
> Postdoctoral Research Fellow<br>
> Mary MacKillop Institute for Health Research<br>
> Australian Catholic University<br>
> ---<br>
> Senior Analyst, Elkhart Group Ltd.<br>
> <a href="http://elkhartgroup.com" target="_blank">http://elkhartgroup.com</a><br>
> Office: <a href="tel:260.673.5518" value="+12606735518">260.673.5518</a><br>
><br>
</div></div>> _______________________________________________<br>
> Rcpp-devel mailing list<br>
> <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
> <a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</blockquote></div><br>
</div></div>