Hello!
I am very new to Rcpp and C++, although I have some experience with the C/R
interface, so please bear with me...
In the context of a Gibbs sampler for a multivariate normal mixture model, I
need to center all the observations by the mean of each component of the
mixture. In pure R I have a function similar to the one below.
> ### 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
+ }
I am trying to transport the code to Rcpp. However, since K, M, and N may be
large and I am running a possibly large number of iterations of a Gibbs
sampler, I want to use a list in the global environment as working space so as
to avoid creating and destroying large matrices at every Gibbs sampler
iteration. So I do the following
> ### For the Rcpp part, set aside working space in the global environment
> workspace <- lapply(1:K, function(k) matrix(0.0, N, M))
Then I try to access the matrices in the 'workspace' list and modify them, but
so far I have had no success. The first attemp, bolow, does not compile...
> 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);
+ }
+
+ return R_NilValue;
+ '
> ## the following does not compile
> library(Rcpp); library(inline)
> f_cpp <- cxxfunction(signature(x = "numeric", mu = "numeric"), code_f, plugin
> = "Rcpp")
file70363a0a2d30.cpp: In function ‘SEXPREC* file70363a0a2d30(SEXP, SEXP)’:
file70363a0a2d30.cpp:42:26: error: invalid conversion from ‘SEXP’ to ‘size_t
{aka unsigned int}’ [-fpermissive]
/usr/local/lib/R/site-library/Rcpp/include/Rcpp/vector/Vector.h:403:18: error:
initializing argument 1 of ‘Rcpp::Vector<RTYPE>::Proxy
Rcpp::Vector<RTYPE>::operator()(const size_t&, const size_t&) [with int RTYPE =
19, Rcpp::Vector<RTYPE>::Proxy = Rcpp::internal::generic_proxy<19>, size_t =
unsigned int]’ [-fpermissive]
file70363a0a2d30.cpp:42:43: error: invalid conversion from ‘SEXP’ to ‘size_t
{aka unsigned int}’ [-fpermissive]
/usr/local/lib/R/site-library/Rcpp/include/Rcpp/vector/Vector.h:403:18: error:
initializing argument 1 of ‘Rcpp::Vector<RTYPE>::Proxy
Rcpp::Vector<RTYPE>::operator()(const size_t&, const size_t&) [with int RTYPE =
19, Rcpp::Vector<RTYPE>::Proxy = Rcpp::internal::generic_proxy<19>, size_t =
unsigned int]’ [-fpermissive]
make: *** [file70363a0a2d30.o] Error 1
ERROR(s) during compilation: source code errors or compiler configuration
errors!
Program source:
1:
2: // includes from the plugin
3:
4: #include <Rcpp.h>
5:
6:
7: #ifndef BEGIN_RCPP
8: #define BEGIN_RCPP
9: #endif
10:
11: #ifndef END_RCPP
12: #define END_RCPP
13: #endif
14:
15: using namespace Rcpp;
16:
17:
18: // user includes
19:
20:
21: // declarations
22: extern "C" {
23: SEXP file70363a0a2d30( SEXP x, SEXP mu) ;
24: }
25:
26: // definition
27:
28: SEXP file70363a0a2d30( SEXP x, SEXP mu ){
29: BEGIN_RCPP
30:
31: Rcpp::NumericMatrix cpp_x(x);
32: Rcpp::NumericMatrix cpp_mu(mu);
33: int K = cpp_mu.ncol();
34: int M = cpp_mu.nrow();
35: int N = cpp_x.nrow();
36: Environment glob = Environment::global_env();
37: List cpp_work(glob.get("workspace"));
38:
39: for (int k=0; k < K; k++) {
40: cpp_work[k] = clone(cpp_x);
41: for (int i=0; i < M; i++)
42: cpp_work(_, i) = cpp_work(_, i) - cpp_mu(i, k);
43: }
44:
45: return R_NilValue;
46:
47: END_RCPP
48: }
49:
50:
Error in compileCode(f, code, language = language, verbose = verbose) :
Compilation ERROR, function(s)/method(s) not created! file70363a0a2d30.cpp:
In function ‘SEXPREC* file70363a0a2d30(SEXP, SEXP)’:
file70363a0a2d30.cpp:42:26: error: invalid conversion from ‘SEXP’ to ‘size_t
{aka unsigned int}’ [-fpermissive]
/usr/local/lib/R/site-library/Rcpp/include/Rcpp/vector/Vector.h:403:18: error:
initializing argument 1 of ‘Rcpp::Vector<RTYPE>::Proxy
Rcpp::Vector<RTYPE>::operator()(const size_t&, const size_t&) [with int RTYPE =
19, Rcpp::Vector<RTYPE>::Proxy = Rcpp::internal::generic_proxy<19>, size_t =
unsigned int]’ [-fpermissive]
file70363a0a2d30.cpp:42:43: error: invalid conversion from ‘SEXP’ to ‘size_t
{aka unsigned int}’ [-fpermissive]
/usr/local/lib/R/site-library/Rcpp/include/Rcpp/vector/Vector.h:403:18: error:
initializing argument 1 of ‘Rcpp::Vector<RTYPE>::Proxy
Rcpp::Vector<RTYPE>::operator()(const size_t&, const size_t&) [with int RTYPE =
19, Rcpp::Vector<RTYPE>::Proxy = Rcpp::internal::generic_proxy<19>, size
In addition: Warning message:
running command '/usr/lib/R/bin/R CMD SHLIB file70363a0a2d30.cpp 2>
file70363a0a2d30.cpp.err.txt' had status 1
Then I tried to modify my code, which now compiles just fine, but does not give
me what I wanted...
> ## try again
> 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++) {
+ NumericMatrix tmpMatrix(as<SEXP>(cpp_work[k]));
+ tmpMatrix = clone(cpp_x);
+ for (int i=0; i < M; i++)
+ tmpMatrix(_, i) = tmpMatrix(_, i) - cpp_mu(i, k);
+ }
+
+ return R_NilValue;
+ '
> ## this time it compiles
> f_cpp <- cxxfunction(signature(x = "numeric", mu = "numeric"), code_f, plugin
> = "Rcpp")
> ## but it does not have the efect I was hoping for, i.e., modifying
> ## the object 'workspace' in the global environment
> f_cpp(x, mu)
NULL
> all.equal(workspace, f_R())
[1] "Component 1: Mean absolute difference: 0.72524" "Component 2: Mean
absolute difference: 1.7288"
Can anybody help me?
Thank you in advance,
Giovanni
Giovanni Petris
Associate Professor
Department of Mathematical Sciences
University of Arkansas - Fayetteville, AR 72701
Ph: (479) 575-6324, 575-8630 (fax)
http://definetti.uark.edu/~gpetris/
_______________________________________________
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
