On 18 October 2012 at 21:54, Davor Cubranic wrote:
| OK, 'sweep' is a fair bit (15x) slower than the Rcpp code Giovanni
| wrote, but you still should be using the return value of the C++
| function to get the result back to R, rather than writing directly into
| function arguments.
|
| Here is a simple benchmark, comparing Giovanni's 'f_cpp' function with
| my solution using 'sweep', with N=10000, M=500, K=10:
|
| > XX <- array(x, dim=c(N, M, K))
| > benchmark(f_cpp=f_cpp(x, mu), sweep=sweep(XX, 2:3, mu),
| columns=c('test', 'relative', 'elapsed'),
| order='relative')
| test relative elapsed
| 1 f_cpp 1.000 69.306
| 2 sweep 15.269 1058.241
|
| I would have also used f_R, but it runs out of memory with these dimensions.
Good idea to benchmark. My "cleaner" approach doesn't shine, it is slower as
Giovanni feared but I still think his (faster!!) approach is distasteful.
R> benchmark(f_cpp=f_cpp(x, mu), sweep=sweep(XX, 2:3, mu), frcpp=f(x, mu,
workspace),
+ columns=c('test', 'relative', 'elapsed'),
+ order='relative')
test relative elapsed
1 f_cpp 1.000 8.947
3 frcpp 2.372 21.222
2 sweep 33.478 299.525
R>
What I would really do here is to create a simple struct or rather class that
upon initialization creates the workspace and holds it.
Doesn't fit as easily in the inline paradigm though.
Dirk
| Davor
|
|
| On 12-10-18 08:02 PM, Davor Cubranic wrote:
| > Giovanni,
| >
| > Anything you pass from R into your C++ routine should be considered
| > read-only. It's typically a bad idea to do in-place updates of the R
| > object you passed in as an argument.
| >
| > In addition, I doubt that the code you wrote (addition/subtraction of
| > matrix columns) is going to be that much faster in C++ over R simply
| > because R already does those operations with optimized C-level code.
| > Maybe not as you wrote it, but you can use available R's functionality a
| > lot more effectively with something like this:
| >
| > sweep(array(x, dim=c(N, M, K)), 2:3, mu)
| >
| > This should do exactly what you're trying to do with your C++ function,
| > except storing it as a 3D array instead of a list of 2D matrices. (It's
| > more efficient and works better with everything else in R.)
| >
| > If it's the large number of iterations in your Gibbs sampler that you're
| > worried about, then write the entire algorithm in C++ and return the
| > final result.
| >
| > Davor
| >
| >
| > On 12-10-18 07:28 PM, Giovanni Petris wrote:
| >> Thanks Dirk!
| >>
| >> By the way, I just figured out a working solution for what I wanted to
| >> achieve - 'code_f' below.
| >>
| >> Best,
| >> Giovanni
| >> =====
| >>
| >> ### Small reproducible example
| >> K <- 2; N <- 5; M <- 3
| >> x <- matrix(rnorm(N * M), N, M)
| >> mu <- matrix(rnorm(M * K, mean = 1 : K), M, K, TRUE)
| >>
| >> ### What I would do in R
| >> f_R <- function() {
| >> ## center observations in 'x' using columns of 'mu'
| >> w <- lapply(1L:K, function(i) x - rep(mu[, i], each = N))
| >> ## eventually more stuff here...
| >> w
| >> }
| >>
| >> ### For the Rcpp part, set aside working space in the global environment
| >> workspace <- lapply(1:K, function(k) matrix(0.0, N, M))
| >>
| >> code_f <- '
| >> Rcpp::NumericMatrix cpp_x(x);
| >> Rcpp::NumericMatrix cpp_mu(mu);
| >> int K = cpp_mu.ncol();
| >> int M = cpp_mu.nrow();
| >> int N = cpp_x.nrow();
| >> Environment glob = Environment::global_env();
| >> List cpp_work(glob.get("workspace"));
| >> NumericMatrix tmpMatrix;
| >>
| >> for (int k=0; k < K; k++) {
| >> tmpMatrix = as<SEXP>(cpp_work[k]);
| >> for (int i=0; i < M; i++)
| >> tmpMatrix(_, i) = cpp_x(_, i) - cpp_mu(i, k);
| >> }
| >>
| >> return R_NilValue;
| >> '
| >> ## library(Rcpp); library(inline)
| >> f_cpp <- cxxfunction(signature(x = "numeric", mu = "numeric"), code_f,
| >> plugin = "Rcpp")
| >>
| >> ## try it
| >> f_cpp(x, mu)
| >> all.equal(workspace, f_R())
| >>
| >>
| >> ________________________________________
| >> From: Dirk Eddelbuettel [e...@debian.org]
| >> Sent: Thursday, October 18, 2012 9:15 PM
| >> To: Giovanni Petris
| >> Cc: Dirk Eddelbuettel; rcpp-devel@lists.r-forge.r-project.org
| >> Subject: RE: [Rcpp-devel] How to modifying the elements of a list in
| >> the global environment?
| >>
| >> On 19 October 2012 at 01:32, Giovanni Petris wrote:
| >> | Hi Dirk,
| >> |
| >> | Thank you for the quick reply. I will look for more examples on the
| >> net.
| >> |
| >> | About your suggestion of allocating scrap space inside the C++
| >> routine, am I wrong to think that when the matrices are large and the
| >> function is called repeatedly within a Gibbs sampler loop, this is not
| >> a very efficient approach?
| >>
| >> Yes, if you want to use it several times you can pass it around for
| >> reuse.
| >> As I mentioned, a reference is just a pointer and hence unrelated to
| >> the size
| >> of the object.
| >>
| >> Dirk
| >>
| >> | Thanks again!
| >> |
| >> | Best,
| >> | Giovanni
| >> |
| >> | ________________________________________
| >> | From: Dirk Eddelbuettel [e...@debian.org]
| >> | Sent: Thursday, October 18, 2012 8:06 PM
| >> | To: Dirk Eddelbuettel
| >> | Cc: Giovanni Petris; rcpp-devel@lists.r-forge.r-project.org
| >> | Subject: Re: [Rcpp-devel] How to modifying the elements of a list in
| >> the global environment?
| >> |
| >> | Giovanni,
| >> |
| >> | On 18 October 2012 at 17:46, Dirk Eddelbuettel wrote:
| >> | | | > code_f <- '
| >> | | | + Rcpp::NumericMatrix cpp_x(x);
| >> | | | + Rcpp::NumericMatrix cpp_mu(mu);
| >> | | | + int K = cpp_mu.ncol();
| >> | | | + int M = cpp_mu.nrow();
| >> | | | + int N = cpp_x.nrow();
| >> | | | + Environment glob = Environment::global_env();
| >> | | | + List cpp_work(glob.get("workspace"));
| >> | | | +
| >> | | | + for (int k=0; k < K; k++) {
| >> | | | + cpp_work[k] = clone(cpp_x);
| >> | | | + for (int i=0; i < M; i++)
| >> | | | + cpp_work(_, i) = cpp_work(_, i) - cpp_mu(i, k);
| >> |
| >> | That can't work. You have the types confused.
| >> |
| >> | You could just pass the list object down, and then pick from it.
| >> Working
| >> | examples for that are in the list archives, in the RcppExamples
| >> package, and
| >> | in other places on the net.
| >> |
| >> | Or, as I suggested, just allocate the scrap space inside the C++
| >> routine.
| >> |
| >> | Hth, Dirk
| >> |
| >> | --
| >> | Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
| >>
| >> --
| >> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
| >> _______________________________________________
| >> 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
| >>
| >
| > _______________________________________________
| > 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
|
| _______________________________________________
| 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
--
Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
_______________________________________________
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
