Yes, we'd LOVE to see this turned into an Rcpp Gallery article once you've got everything ironed out!
Details on creating Gallery articles here: https://github.com/jjallaire/rcpp-gallery/wiki/Contributing-to-the-Rcpp-Gallery On Tue, Apr 14, 2015 at 5:36 PM, Joshua Wiley <jwiley.ps...@gmail.com> wrote: > Hi JJ, > > Thank you, that worked wonderfully. > > 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. > > I've put it online in case anyone else runs across this and is interested: > > https://gist.github.com/JWiley/d9cba55603471f75d438 > https://gist.github.com/JWiley/8fe93ae105b1fb244f93 > > 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. > > Thanks again for your help and the great package! > > Cheers, > > Josh > > > > > > > > On Tue, Apr 14, 2015 at 7:54 PM, JJ Allaire <jj.alla...@gmail.com> wrote: >> >> The initialization of your const fields needs to happen within the >> initializer list (rather than inline where the fields are declared). >> If you substitute this code for the code in your example everything >> will compile: >> >> // QR Decomposition >> const CPivQR PQR; >> const Permutation Pmat; >> const int r; >> >> // output >> RMatrix<double> betamat; >> Eigen::VectorXd betahat; >> >> // initialize with input and output >> CVLM(const Eigen::VectorXd y, const Eigen::MatrixXd X, >> Rcpp::NumericMatrix betamat) >> : y(y), X(X), PQR(X), Pmat(PQR.colsPermutation()), r(PQR.rank()), >> betamat(betamat) {} >> >> >> >> On Tue, Apr 14, 2015 at 4:41 AM, Joshua Wiley <jwiley.ps...@gmail.com> >> wrote: >> > Hi, >> > >> > I'm interested in combining faster linear model estimation via QR >> > decompositions from RcppEigen with JJs new(er) RcppParallel package to >> > do >> > cross validation and bootstrapping on linear models. >> > >> > As a first step, I've been trying to merge the JSS paper on RcppEigen ( >> > http://www.jstatsoft.org/v52/i05/paper ) with about the documentation >> > for >> > RcppParallel ( http://rcppcore.github.io/RcppParallel/#examples ) >> > >> > ignoring bootstrapping and cross validation for the time being and just >> > getting a linear model to run in the parallel framework. I've toyed >> > with >> > each separately successfuly, and think I am somewhat close to getting >> > them >> > together, but I am running into an error: >> > >> > partest.cpp:28:20: error: 'X' is not a type >> > partest.cpp:29:26: error: 'PQR' is not a type >> > partest.cpp:29:29: error: expected ',' or '...' before '.' token >> > partest.cpp:30:15: error: 'PQR' is not a type >> > partest.cpp:30:18: error: expected ',' or '...' before '.' token >> > partest.cpp: In member function 'virtual void >> > CVLM::operator()(std::size_t, >> > std::size_t)': >> > partest.cpp:44:28: warning: comparison between signed and unsigned >> > integer >> > expressions [-Wsign-compare] >> > make: *** [partest.o] Error 1 >> > >> > My code is below (run via sourceCpp() ), and I can see exactly where the >> > first error occurs, but I don't understand why this is OK when run >> > outside >> > of RcppParallel, but so very wrong in conjunction with RcppParallel. >> > >> > Any pointers on either what's wrong or where I should be reading would >> > be >> > deeply appreciated. >> > >> > Josh (code below) >> > >> > >> > // [[Rcpp::depends(RcppParallel)]] >> > // [[Rcpp::depends(RcppEigen)]] >> > #include <Rcpp.h> >> > #include <RcppEigen.h> >> > #include <RcppParallel.h> >> > >> > using namespace std; >> > using namespace Rcpp; >> > using namespace RcppParallel; >> > >> > using Eigen::Lower; >> > using Eigen::Map; >> > using Eigen::MatrixXd; >> > using Eigen::Upper; >> > using Eigen::VectorXd; >> > typedef Map<MatrixXd> MapMatd; >> > typedef Map<VectorXd> MapVecd; >> > typedef Eigen::ColPivHouseholderQR<MatrixXd> CPivQR; >> > typedef CPivQR::PermutationType Permutation; >> > >> > struct CVLM : public Worker >> > { >> > // design matrix and outcome >> > const Eigen::VectorXd y; >> > const Eigen::MatrixXd X; >> > >> > // QR Decomposition >> > const CPivQR PQR(X); // ERROR HAPPENS HERE >> > const Permutation Pmat(PQR.colsPermutation()); >> > const int r(PQR.rank()); >> > >> > >> > // output >> > RMatrix<double> betamat; >> > Eigen::VectorXd betahat; >> > >> > // initialize with input and output >> > CVLM(const Eigen::VectorXd y, const Eigen::MatrixXd X, >> > Rcpp::NumericMatrix >> > betamat) >> > : y(y), X(X), betamat(betamat) {} >> > >> > void operator()(std::size_t begin, std::size_t end) { >> > >> > //betahat = PQR.solve(y); >> > for(int j = begin; j < end; j++) { >> > betamat(j, j) = X(j, j) + y[j]; >> > } >> > } >> > }; >> > >> > // [[Rcpp::export]] >> > NumericMatrix parallelFit(Eigen::VectorXd dv, Eigen::MatrixXd x) { >> > >> > // allocate the output matrix >> > NumericMatrix betamat(x.rows(), x.cols()); >> > >> > // pass input and output >> > CVLM cvLM(dv, x, betamat); >> > >> > // parallelFor to do it >> > parallelFor(0, x.rows(), cvLM); >> > >> > return betamat; >> > } >> > >> > >> > >> > >> > >> > >> > >> > >> > >> > >> > >> > >> > -- >> > Joshua F. Wiley >> > http://joshuawiley.com/ >> > --- >> > Postdoctoral Research Fellow >> > Mary MacKillop Institute for Health Research >> > Australian Catholic University >> > --- >> > Senior Analyst, Elkhart Group Ltd. >> > http://elkhartgroup.com >> > Office: 260.673.5518 >> > >> > _______________________________________________ >> > Rcpp-devel mailing list >> > Rcpp-devel@lists.r-forge.r-project.org >> > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel > > _______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
