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
