[Rcpp-devel] Rcpp Parallel and Rcpp Armadillo

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.

-----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;
}

----------------------------------------------------------

[Rcpp-devel] Rcpp Parallel and Rcpp Armadillo

Thread (2 messages)