1

I have a section of code in that is running quite slowly, and so I'm hoping to rewrite the function in , however, this is my first attempt using Rcpp, and I cannot get the code to compile or run.

In R, the section I am trying to rewrite is

  drawsamples<-lapply(1:numstudies,function(v){
      temp<-t(mapply(rmvn,n=1,mu=finalmu_b2_l[[v]],sigma=finalcov_b2_l[[v]]))
      rownames(temp)<-id[[v]]
      temp
    })

This code should return a nested ist of matrices. The list should be of length numstudies and each element of it should be a of dimension n[v] rows by ql columns. Each row of these matrices should be draws from a multivariate normal distribution, however the draws for each row will be controlled by a different mean vector finalmu_b2_l and a different covariance matrix finalcov_b2_l. Both finalmu_b2_l and finalcov_b2_l are nested lists, top level of length numstudies, each element of which is a list of length n[v].

So far, I have tried to write a piece of code that will just draw one of the matrices needed, but want to extend this to return a list of matrices as soon as possible. The code is below:

#include <RcppDist.h>
#include <RcppArmadillo.h>

using namespace Rcpp;

// [[Rcpp::depends(RcppDist,RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat draw_b(const int n,
                 const int ql,
                 const Rcpp::List mu,
                 const Rcpp::List cov) {
  arma::mat draws = arma::zeros(n,ql);
  for (int iter = 1; iter <= n; iter++) {
     draws.row(iter) = rmvnorm(1, mu[iter], cov[iter]);
  }
  return draws;
}

Any time I try to compile the code (by sourcing the .cpp file the code is held in, using function sourceCpp), I get an error reading:

invalid initialization of reference of type 'const::mat&'

, which I understand as an issue either with creating the draws matrix or with filling it?

Any advice, pointers or guidance on -

  1. why this code is not compiling, what I should be doing in Rcpp instead, and/or
  2. how to expand this to return a list of matrices rather than just a single matrix would be greatly appreciated.

EDIT

The original r code contains a mapply function nested within a lapply function. I wanted to have an output that is a list length numstudies, each element of which is a matrix. Each row of each matrix within this list is a relisation from a different multivariate normal distribution (so within each list element, for each row of the matrix, there is a unique mean vector and a unique covariance matrix controlling the multivariate normal I want to draw from). I'd written as a mapply nested within a lapply to automatically give the output format I wanted, whilst allowing draws from different distributions for each matrix row.

EDIT2

After changing the iteration to run from 0 instead of 1, the following code compiles:

#include <RcppDist.h>
#include <RcppArmadillo.h>

using namespace Rcpp;

// [[Rcpp::depends(RcppDist,RcppArmadillo)]]


//' @keywords internal
// [[Rcpp::export]]
arma::mat draw_b(const int & n,
                 const int & ql,
                 const Rcpp::List & mu,
                 const Rcpp::List & cov) {
  arma::mat draws = arma::zeros(n,ql);
  for (int iter = 0; iter <= n; iter++ ) {
     draws.row(iter) = rmvnorm(1, mu_temp, cov_temp);
  }
  return draws;
}

EDIT3

The code currently compiles, but does not draw samples. Instead I get the following error message:

error: Mat::init(): requested size is not compatible with column vector layout
Error in draw_b(n = n, ql = ql, mu = mu_example, cov = cov_example) : 
  Mat::init(): requested size is not compatible with column vector layout

I've prepared an example of what I would like this basic function to do (that just samples a matrix of realisations from different multivariate normal distributions.

Data:

list(n = 10, ql = 2, mu_example = list(c(0.342909965003207, -0.788070875792469
), c(-0.00499810116271365, 0.0114865660452949), c(-0.149753928200309, 
0.344162379034614), c(0.335176829763227, -0.770298692761465), 
    c(0.254206123984596, -0.584212951520601), c(0.379893097582703, 
    -0.873064992779416), c(0.137231089484867, -0.315382566602526
    ), c(0.405123380985852, -0.931048876501857), c(-0.00505917031396947, 
    0.0116269143128456), c(-0.0743318653279181, 0.170828451158346
    )), cov_example = list(structure(c(0.199912910315971, -0.459437048770092, 
-0.459437048770092, 4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971, 
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L, 
2L)), structure(c(0.199912910315971, -0.459437048770092, -0.459437048770092, 
4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971, 
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L, 
2L)), structure(c(0.199912910315971, -0.459437048770092, -0.459437048770092, 
4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971, 
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L, 
2L)), structure(c(0.199912910315971, -0.459437048770092, -0.459437048770092, 
4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971, 
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L, 
2L)), structure(c(0.199912910315971, -0.459437048770092, -0.459437048770092, 
4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971, 
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L, 
2L))))

R code

  drawsampletest<-draw_b(n=n,
                           ql=ql,
                           mu=mu_example,
                           cov=cov_example)

Rcpp code

#include <RcppDist.h>
#include <RcppArmadillo.h>

using namespace Rcpp;

// [[Rcpp::depends(RcppDist,RcppArmadillo)]]


//' @keywords internal
// [[Rcpp::export]]
arma::mat draw_b(const int & n,
                 const int & ql,
                 const Rcpp::List & mu,
                 const Rcpp::List & cov) {
  arma::mat draws = arma::zeros(n,ql);
  for (int iter = 0; iter <= n; iter++ ) {
    arma::rowvec mu_temp = mu[iter];
    arma::mat cov_temp = cov[iter];
    draws.row(iter) = rmvnorm(1, mu_temp, cov_temp);
  }
  return draws;
}

Of course once this is working- I still need to extend this to draw a list of matrices rather than a single matrix

8
  • Have you tried filling columns instead of rows? Commented Apr 2, 2020 at 11:00
  • 1
    Also, make sure to use indices starting at 0. Commented Apr 2, 2020 at 11:02
  • 1
    I'm not sure your R code is doing what you think it is (or else I'm misunderstanding the description of your problem). Why do you have a mapply() call within the lapply() call? I think this question could benefit from a full minimal reproducible example. For tips on making one in the context of R, see How to make a great R reproducible example. In particular, adding your data (i.e. numstudies, finalmu_b2_l, and finalcov_b2_l) by editing your question to include the output of the R command dput() and adding some expected output would probably help. Commented Apr 2, 2020 at 11:36
  • 1
    I say this as general advice, but also b/c it's possible this is an XY problem; what you're really worried about is speed of execution, not that it occurs in C++, right? It could be that tweaking your R code improves almost as good as rewriting in C++ would. So, even if we could fix w/ your C++ code & advise you how to extend it to return a list w/o seeing your data, your data (or a subset of it, or made up data with the same structure) can help produce a higher quality answer by helping us understand the problem. Commented Apr 2, 2020 at 11:43
  • Also, by the way, your code compiles just fine for me? After making the change prompted by @F.Privé 's astute observation re: indexing, it seemed to function as expected as well. You might want to add your operating system, and R, Rcpp, and RcppDist versions (such as via sessionInfo()) Commented Apr 2, 2020 at 11:47

1 Answer 1

2

Here's a basic setup that should do what you want:

#include <RcppDist.h>
#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppDist,RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat draw_b(const int ql,
                 const Rcpp::List& mu,
                 const Rcpp::List& cov) {
    int n = mu.size();
    arma::mat draws = arma::zeros(n, ql);
    for ( arma::uword iter = 0; iter < n; iter++ ) {
        draws.row(iter) = rmvnorm(1, mu[iter], cov[iter]);
    }
    return draws;
}

// [[Rcpp::export]]
Rcpp::List get_list_of_draws(Rcpp::List mu, Rcpp::List Sigma, int ql) {
    int numstudies = mu.size();
    Rcpp::List res(numstudies);
    for ( int iter = 0; iter < numstudies; ++iter ) {
        Rcpp::List mu_temp = mu[iter];
        Rcpp::List cov_temp = Sigma[iter];
        res[iter] = draw_b(ql, mu_temp, cov_temp);
    }
    return res;
}

It seems to work as expected:

res <- get_list_of_draws(mu_example, cov_example, ql)
str(res)
# List of 1
#  $ : num [1:10, 1:2] -0.0156 -0.4717 -0.8134 0.5489 0.1215 ...

(Note though that when I set up mu_example and cov_example I wrapped them in list()s as you stated they should be lists of lists.)

Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.