[Rcpp-devel] Rcpp Parallel and Rcpp Armadillo

[Dirk Eddelbuettel]

On 10 December 2014 at 14:54, Maxime To wrote:
| Ok, thanks, I'll try to build on it. In the example I tried to isolate the problem, but in my real program I have lot of other matrix step using armadillo, that's why I put it in that way... I'd like to avoid armadillo, but it makes matrix calculus messier...
| 
| About the memory, it does not seem to work, when I run a long loop the program just crash because of full memory.

The Writing R Extensions manual has some material on how to use tools like valgrind.

Dirk

| 
| -----Message d'origine-----
| De : "Romain François" <romain at r-enthusiasts.com>
| Envoyé : ‎10/‎12/‎2014 14:32
| À : "Maxime To" <maxime.to at outlook.fr>
| Cc : "Dirk Eddelbuettel" <edd at debian.org>; "rcpp-devel at lists.r-forge.r-project.org" <rcpp-devel at lists.r-forge.r-project.org>
| Objet : Re: [Rcpp-devel] Rcpp Parallel and Rcpp Armadillo
| 
| 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
| > CC: maxime.to at outlook.fr; rcpp-devel at lists.r-forge.r-project.org
| > Subject: Re: [Rcpp-devel] Rcpp Parallel and Rcpp Armadillo
| > From: 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
| _______________________________________________
| 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
-- 
http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org