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: q...@umail.iu.edu
> CC: maxime...@outlook.fr; rcpp-devel@lists.r-forge.r-project.org
> Subject: Re: [Rcpp-devel] Rcpp Parallel and Rcpp Armadillo
> From: e...@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 | e...@debian.org
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
