[Rcpp-devel] Rcpp Parallel and Rcpp Armadillo

[Romain François]

Some pointers. 

When you use an arma::mat passed by value in an Rcpp::export, this means copying all of the data of the underlying R object into armadillo. I’d suggest you use a reference to const to avoid that, i.e. 

mat contrib1(const mat& X1) { … }

Then in pQnorm, you do: 

NumericMatrix x_q = Rcpp::as<Rcpp::NumericMatrix>(wrap(xx_q));

That is yet again, copying all of the data from the arma::mat into an Rcpp matrix. 

You then return a arma::mat, which data is copied implicitly as the return of contrib1. 

I’d suggest you do all this without armadillo, which you don’t really use except for inducing a lot of extra copies of data. 

To anwser your last question, R uses a garbage collector, so the memory is not automatically reclaimed as soon as it is no longer needed. 

Hope this helps. 

Romain

> Le 10 déc. 2014 à 15:01, Maxime To <maxime.to at outlook.fr> a écrit :
> 
> Hi, 
> 
> I changed the function as indicated by Dirk and I modify the functions and the program does work now.
> However, I am still puzzled by the memory use of the program. when I run a loop of my function in R as in the code below, it seems that the program does not free the memory used in the previous iterations... which is annoying when I need to optimize on my final object.
> 
> So I was wondering whether it was a question of declaration of object in my code?
> 
> ------------------------------------------------------------------------------------------------------------------
> 
> sourceCpp("Rcpp/test.cpp") #
> qwe = matrix(runif(10000), nrow = 100)
> a = contrib1(qwe)
> b = qnorm(qwe)
> a - b
> 
> for (i in 1:20000) a = contrib1(qwe)
> ----------------------------------------------------------
> // test.cpp
> 
> #include <RcppArmadillo.h>
> #include <cmath>
> #include <algorithm>
> #include <RcppParallel.h>
> #include <boost/math/distributions/inverse_gaussian.hpp>
>  
> using namespace Rcpp;
> using namespace arma;
> using namespace std;
> using namespace RcppParallel;
>  
> // [[Rcpp::depends(RcppArmadillo, RcppParallel, BH)]]
>  
> double qnorm_f(const double& x_q) {
>     boost::math::normal s;
>     return boost::math::quantile(s, x_q);
> };
> 
>  
>  
> struct Qnorm : public Worker
> {
>    // source matrix
>    const RMatrix<double> input_q;
>     
>    // destination matrix
>    RMatrix<double> output_q;
>     
>    // initialize with source and destination
>    Qnorm(const NumericMatrix input_q, NumericMatrix output_q)
>       : input_q(input_q), output_q(output_q) {}
>     
>    // take the Pnorm of the range of elements requested
>    void operator()(std::size_t begin, std::size_t end) {
>       std::transform(input_q.begin() + begin,
>                      input_q.begin() + end,
>                      output_q.begin() + begin,
>                      ::qnorm_f);
>    }
> };
>  
> mat pQnorm(mat xx_q) {
>  
>     NumericMatrix x_q = Rcpp::as<Rcpp::NumericMatrix>(wrap(xx_q));
>    
>     // allocate the output matrix
>     const NumericMatrix output_q(x_q.nrow(), x_q.ncol());
>    
>     // Pnorm functor (pass input and output matrices)
>     Qnorm qnorm_temp(x_q, output_q);
>    
>     // call parallelFor to do the work
>     parallelFor(0, x_q.length(), qnorm_temp);
>    
>     // return the output matrix
>     mat outmat_q(output_q.begin(), output_q.nrow(),output_q.ncol());
>     return outmat_q;
>  
> }
>  
> // [[Rcpp::export]]
> mat contrib1(mat X1) {
>  
>     mat test    = pQnorm(X1);
>     mat results = test;
> 
>     return results;
> }
> 
> ----------------------------------------------------------
> 
> > Date: Tue, 9 Dec 2014 09:07:10 -0600
> > To: qkou at umail.iu.edu <mailto:qkou at umail.iu.edu>
> > CC: maxime.to at outlook.fr <mailto:maxime.to at outlook.fr>; rcpp-devel at lists.r-forge.r-project.org <mailto:rcpp-devel at lists.r-forge.r-project.org>
> > Subject: Re: [Rcpp-devel] Rcpp Parallel and Rcpp Armadillo
> > From: edd at debian.org <mailto:edd at debian.org>
> > 
> > 
> > On 9 December 2014 at 09:46, Qiang Kou wrote:
> > | What do you mean by "doesn't work" ? Compiling error or the result is not
> > | right?
> > | 
> > | I just tried the code, and it seems the code can compile and work.
> > 
> > I am generally very careful about calling back to anything related to R from
> > functions to be parallelized. So for
> > 
> > inline double f(double x) { return ::Rf_pnorm5(x, 0.0, 1.0, 1, 0); }
> > 
> > I think going with an equivalent pnorm() function from Boost / Bh may be better.
> > 
> > But I am shooting from my hip here as I have not had time to look at this,
> > having been out way too late at a nice concert :) 
> > 
> > Dirk
> > 
> > -- 
> > http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org
> _______________________________________________
> 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20141210/fbd549b5/attachment-0001.html>