[Rcpp-devel] Is there a way to pass references to R

[Søren Højsgaard]

Dear all,

To supplement Dirk's remarks about using RcppEigen in connection with sparse matrices: Here is an example where I use RcppEigen to moralized a directed acyclic graph represented by an adjacency matrix which is sparse (a dgCMatrix from R's Matrix package). The code can actually be made more efficient by replacing the for-loops by a different type of iterator which ensures that only non-zero entries are visited. The "triplet-trick" is the "recommended" way of building up a sparse matrix - as far as I can tell. There may be many other possible improvents which I am unaware of (as I don't know c++).

Best regards
Søren


RcppExport SEXP C_moralizeM ( SEXP XX_){
  using Eigen::Map;
  using namespace Rcpp;
  typedef Eigen::MappedSparseMatrix<double> MSpMat;
  typedef Eigen::SparseMatrix<double> SpMat;
  SpMat   X(as<MSpMat>(XX_));
  
  typedef Eigen::Triplet<double> T;
  std::vector<T> triplets;
  triplets.reserve(X.nonZeros() * 2);
  
  int nrX(X.rows());
  int kk, ll, vv;
  for (vv=0; vv<nrX; vv++){ /* consider vertex vv */
    for (kk=0; kk<nrX; kk++){
      if (X.coeff(kk, vv) != 0){     /* yes, kk->vv */
	for (ll=kk+1; ll<nrX; ll++){
	  if (X.coeff(ll, vv) != 0){ /* yes, ll->vv */
	    if ((X.coeff(kk, ll)==0) && (X.coeff(ll, kk)==0)){ /* kk not~ ll */
	      triplets.push_back(T(kk, ll, 1));
	      triplets.push_back(T(ll, kk, 1));
	    }
	  }
	}
      }
    }
  }
  
  SpMat ans(X.rows(), X.cols());
  ans.setFromTriplets(triplets.begin(), triplets.end());
  SpMat Xt(X.transpose());
  ans = ans + Xt + X;
  
  for (kk=0; kk<nrX; kk++){
    for (ll=kk+1; ll<nrX; ll++){
      if (ans.coeff(kk,ll)!=0){
	ans.coeffRef(kk,ll)=1;
	ans.coeffRef(ll,kk)=1;
      }
    }
  }
  ans.makeCompressed();
  return(wrap(ans));
}

-----Original Message-----
From: rcpp-devel-bounces at lists.r-forge.r-project.org [mailto:rcpp-devel-bounces at lists.r-forge.r-project.org] On Behalf Of Dirk Eddelbuettel
Sent: 22. december 2012 01:41
To: matzke at berkeley.edu
Cc: rcpp-devel at lists.r-forge.r-project.org
Subject: Re: [Rcpp-devel] Is there a way to pass references to R


Nick,

On 21 December 2012 at 16:16, Nick Matzke wrote:
| I'm working on a problem that involves calculating a very large, 
| sparse transition probability matrix, P.  It takes significant time to 
| run through the whole matrix, check all the transitions, and put in 
| the appropriate probabilities, based on multiplying parameter values 
| j, d, e, etc.
| 
| Each time I update j, or d, or e, I would rather not go re-process 
| which cells of P get j, j*d, d*e, etc.  It would be better to just do 
| this once, with j, d, and e each just referring to a memory address.  
| I could then update the values at that memory address and re-use the 
| matrix P without re-building it again and again.
| 
| Currently, my script builds P via an Rcpp function, then passes it to 
| R, where it gets used (with the values
| calculated) a few hundred times (via another Rcpp function); then j, 
| d, and e are updated and I have to re-calculate P.
| 
| So, my question is, is there a way to make P, currently a vector of 
| float vectors, into a vector of float pointers, which can be passed to 
| R and then back to Rcpp?  Or is this just ridiculous?

There will be lots of different ways to skin this cat. Sparse matrices may be one, I have no experience there. They now exist in Armadillo, though without [at current] glue code for as<> and wrap in RcppArmadillo; this does exist in Eigen and RcppEigen. Others may help you on sparse stuff.

You can also define you own structure to hold those vectors, and then simply pass one single pointer around using the Rcpp::XPtr wrapper class for R's external pointer. If data volume is an issue, external pointers are your friends. (Some of) the database class package uses this, and as does the bigmemory familiy of packages.  

But you should do some measuring and profiling before you go off and rearchitect your application based on a hunch.  Great slide deck seen today and retweeted:  

  The only good intuition: 'I should time time this' 
  -- Andrei Alexandrescu

  Via http://isocpp.org/blog/2012/12/three-optimization-tips-alexandrescu

and I think he has that exactly right.  And he has about a gazillion times as much street cred on C++ as I do ....

Dirk

PS You need to work on your signature. Not sure it's long and detailed enough.

--
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com _______________________________________________
Rcpp-devel mailing list
Rcpp-devel at lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel