0

When I try to multiply an arma::mat and a NumericVector in Rcpp using operator*, I get the following error:

no match for operator*

Here is an example of what I'm trying to multiply:

NumericVector exampleNumericVector(L-1);
arma::mat exampleMatrix(L-1,L-1);
exampleMatrix*(Rcpp::qnorm(exampleNumericVector,0,1,true,false));

I have tried looking through the excellent Armadillo documentation, and the closest I could find to a solution is the prod() function, but this has not corrected my issue. I tried overloading the operator*, as well. However, this was also unsuccessful.

Update: Reproducible example

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
arma::vec foo(Rcpp::NumericVector x) {
  // [[Rcpp::depends(RcppArmadillo)]]
  Rcpp::NumericVector y = Rcpp::qnorm(x, 0, 1, true, false);
  int n = y.size();
  arma::mat M(n,n);
  M.fill(1.0);                  // nonsense content, we don't care
  arma::vec v(y);               // instantiate Arma vec from Rcpp vec
  arma::vec res = M * v;
  return res;
}

// [[Rcpp::export]]
Rcpp::List foo(const int n, const int M, const int L, const double mu){
  
  // initialize variables
  arma::mat P = arma::zeros(M,L);
  arma::mat SigInvTmp = arma::zeros(L-1,L-1);
  double muTmp = 0;
  double s2Tmp = 0;
  arma::mat Sig = arma::zeros(L,L);
  arma::mat Sigee = arma::zeros(L-1,L-1);
  arma::mat Sigie = arma::zeros(1,L-1);
  arma::mat Sigei = arma::zeros(L-1,1);
  Rcpp::NumericVector Pie(L-1);
  
  for(int i = 0; i < M; i++){
    
    arma::uvec iIdx(1); // vector of index to take for i
    iIdx(0) = i;
    
    for(int l = 0; l < L; l++){
      
      srand(time(0)); 
      
      // We need to get the vector of indices excluding l
      arma::uvec idx(L-1); // vector of indices to subset by
      arma::uword ii = 0;  // the integer we'll add at each elem of idx
      for (arma::uword j = 0; j < (L-1); ++j) {
        if (ii == l) {   // if ii equals l, we need to skip l
          ii += 1;     
        }
        idx[j] = ii;     // then we store ii for this elem
        ii += 1;         // and increment ii
      }
      
      // Single index for l
      arma::uvec lIdx(1); // vector of index to take for l
      lIdx(0) = l;
      
      Sigee = Sig.submat(idx,idx);
      Sigie = Sig.submat(lIdx,idx);
      Sigei = Sig.submat(idx,lIdx);
      Pie = P.submat(iIdx,idx);
      
      SigInvTmp = inv(Sigee);
      arma::vec qnormVec = foo(Pie);
      muTmp = mu + as_scalar(Sigie*SigInvTmp*(qnormVec-mu));
      s2Tmp = sqrt(as_scalar(Sig(l,l)-Sigie*SigInvTmp*Sigei));
    }
  }
  return Rcpp::List::create(Rcpp::Named("mu") = muTmp, Rcpp::Named("s2") = s2Tmp);
}
Ron Snow
  • 239
  • 1
  • 7
  • 1
    Likely it is simplest to cast the matrix into a NumericMatrix or the vector into armadillo matrix. – Oliver Aug 28 '20 at 22:56
  • @Oliver I was considering both of those options with a preference for the latter. However, I can't find the way to convert the NumericVector to an Armadillo matrix – Ron Snow Aug 28 '20 at 23:04
  • @Oliver Is typecasting the same for going from numericVector to arma::mat? I can't find the function to do so – Ron Snow Aug 30 '20 at 18:40
  • 1
    Hi @RonSnow. Dirk's answer is a direct method your specific solution. In c++ typecasting is done using `(new-type)variable` possible with a reference pointer. It does require that a method has been defined that can cast the method from one to another. His answer is a much more general and usable. While Dirk can be very forward in his attitude (mainly due to the amount of answers he gives I believe), his answers are brilliant. For your specific problem i would move `// [[Rcpp::depends(RcppArmadillo)]]` down to the actual function, while I am not sure this is the answer. – Oliver Aug 30 '20 at 18:59
  • @Oliver Thank you for following up. I implemented Dirk's function, but I am still receiving the error "static assertion failed: cannot convert type to SEXP." I am thinking this is a result of my line of code "muTmp = mu + as_scalar(Sigie*SigInvTmp*(qnormVec-mu));" What do you think? – Ron Snow Aug 30 '20 at 19:09

1 Answers1

4

Your example is neither complete, nor minimal, nor reproducible.

So in an attempt to approve on this, a full example. You feed is vector of p values to obtain quantile values for, it creates a (nonsense) matrix of suitable dimension, creates an arma vector from the Rcpp vector (key: they all convert to/from SEXP) and multiplies. The result is returned.

Code
#include <RcppArmadillo.h>

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

// [[Rcpp::export]]
arma::vec foo(Rcpp::NumericVector x) {
  Rcpp::NumericVector y = Rcpp::qnorm(x, 0, 1, true, false);
  int n = y.size();
  arma::mat M(n,n);
  M.fill(1.0);                  // nonsense content, we don't care
  arma::vec v(y);               // instantiate Arma vec from Rcpp vec
  arma::vec res = M * v;
  return res;
}

/*** R
foo(c(0.95, 0.975))
*/
Demo
R> sourceCpp("~/git/stackoverflow/63641766/answer.cpp")

R> foo(c(0.95, 0.975))
        [,1]
[1,] 3.60482
[2,] 3.60482
R> 
Dirk Eddelbuettel
  • 360,940
  • 56
  • 644
  • 725
  • Thank you for your example. I tried applying your method to my problem, and I'm getting the following error: static assertion failed: cannot convert type to SEXP. – Ron Snow Aug 29 '20 at 16:21
  • 2
    Please read https://stackoverflow.com/help/minimal-reproducible-example -- our mindreader device is in the shop so we can't tell what code you wrote. – Dirk Eddelbuettel Aug 29 '20 at 16:24
  • My apologies- you have a good point. I updated the post to include a minimal-reproducible-example. – Ron Snow Aug 29 '20 at 16:46