[Rcpp-devel] Integrating RcppParallel and RcppEigen

[Joshua Wiley]

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>