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
