[Rcpp-devel] Integrating RcppParallel and RcppEigen

Joshua Wiley jwiley.psych at gmail.com
Tue Apr 14 10:41:07 CEST 2015


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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150414/28af7b83/attachment-0001.html>


More information about the Rcpp-devel mailing list