Search code examples
c++rrcpparmadillo

Getting RcppArmadillo working in Rstudio using .cpp file


I'm a newbie with this, so sorry for the basic question. I have managed to get Rcpp working but need a function in Armadillo, and I am just not getting it right. I know the code is bad, like I said I am a newbie in C++, but I would be very grateful if someone could help me get this example going. My code from the .cpp file in Rstudio is as follows. I am running latest R, and Rstudio. Thanks.

#include <iostream>
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

/*** R
require(Rcpp)
require(RcppArmadillo)
#Population of size N
N = 350
#Set size
m = 3
#Number of cycles
C = 20
#Final number of analysed units
n = C * m
#make the r vector of ranks
r <- rep(1:m, length.out = n)

#initialise alpha matrix nxn
alpha <- matrix(0, ncol = N, nrow = n)
  */

// [[Rcpp::export]]
unsigned nChoosek( unsigned n, unsigned k )
{
  if (k > n) return 0;
  if (k * 2 > n) k = n-k;
  if (k == 0) return 1;

  int result = n;
  for( int i = 2; i <= k; ++i ) {
    result *= (n-i+1);
    result /= i;
  }
  return result;
}

// [[Rcpp::export]]
NumericMatrix alphaC(NumericVector r, int N, int m, int n, int C, 
NumericMatrix alpha){
  int S = r.size();
  int a = 0;
  int b = 0;
  double c = 0;

  for(int i = 0; i < N; ++i){

    for (int j = 0; j < S; ++j){

      a = nChoosek(i, (r[j] - 1));
      b = nChoosek((N - (i + 1)),( m - r[j]));
      c = nChoosek(N, m);

      alpha(j,i) = a*b/c;  

      //std::cout << "a*b/c: " << a*b/c << std::endl;
    }
  }
  return(alpha);
}

/*** R
#Create first order probabilities
pi <- 1 - apply((1 - alphaC(r,N,m,n,C,alpha)), 2, prod)
  plot(pi)
  */

/*** R
#Second order inclusion properties.
pii <- matrix(0, ncol = N, nrow = N)
  */

// [[Rcpp::export]]

int times(NumericVector vec){
  int sum = 1;
  sum = arma::prod(vec);
  return sum;
}

/*** R
vec = 1:5
times(vec)
*/

Solution

  • Your example code contains many things that are not relevant here. Please try to provide a minimal example next time, together with the error message you got. Anyway, you will need an Armadillo data structure to work with arma::prod. You can either create this directly as shown by @mt1022 in the comments, or have it done implicitly by defining a function with suitable arguments:

    #include <RcppArmadillo.h>
    // [[Rcpp::depends(RcppArmadillo)]]
    
    // [[Rcpp::export]]
    double times(arma::vec vec){
      double result = arma::prod(vec);
      return result;
    }
    
    /*** R
    vec <- 1:5 + 0.5
    times(vec)
    */
    

    Note that I have changed the return type to double, since the product of a numeric vector is in general not representable as a int.