rloopsrcpp

Converting code from R to C++ using Rcpp not returning any results


I've had success in using vectors, but I understand that I could power up my simulations using C++ to obtain a higher level.

Here follows the original for loop in R which runs pretty well:

library(plyr)
library(dplyr)
library(tidyverse)
library(DT)
library(Rcpp)
library(scales)

Mandela <- rpois(2000000, 0.01)
mean <- mean(Mandela)
sigma <- sd(Mandela)

#colocar distribui??o para adquirir par?metros
paths <- 100000
count<-12*30
interval<-12/count
sample <- matrix(0,nrow=(count+1),ncol=paths)
sample[1,]<- 0.01
rnormm <- as.matrix(sapply(1:(count+1), function(x) rnorm(1, mean = 0, sd = 1)))

t0 <- Sys.time()

for(i in 1:paths) {
  for(j in 2:(count+1))
  {
    sample[j,i]<-sample[j-1,i]*exp(interval*(mean-((sigma)^2)/2)+((interval)^.5)*rnormm[j-1,1]*sigma) #Expression for Geometric Brownian Motion
  }
}

t1 <- Sys.time()

tempoatual <- t1-t0

This chunk returns me the matrix I actually need to proceed with my calculus.

But when I convert it into C++ using Rcpp, my attempt results in failure, not returning any values.

cppFunction(
  
  "IntegerMatrix proc(IntegerMatrix x, IntegerMatrix random, int media, int variancia, int intervalo) {

  const int n = x.nrow();
  const int m = x.ncol();
  const int mean = media;
  const int sigma = variancia;
  const int interval = intervalo;
  IntegerMatrix y(n,m);
  IntegerMatrix z(random.nrow(),1);

  y = x;
  z = random;

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

    for (int i=1; i < n; i++) {

      y(i,j) = y(i-1,j) * exp(interval*(mean-((sigma)^2)/2)+((interval)^(1/2))*z(i-1,1)*sigma);

    }

  }

  return y;

}")

tempoatual <- t1-t0

t0 <- Sys.time()

teste <- proc(x = sample,random = rnormm,media = mean,variancia = sigma,intervalo = interval)

t1 <- Sys.time()

tempoRCPP <- t1-t0

Could you advise me where it all went wrong?


Solution

  • Your sample contains double values, which are coerced to integer when you pass them to the C++ function. as.integer(0.01) returns 0L. Replace IntegerMatrix with NumericMatrix and unique(c(teste)) will return [1] 0.01 NaN. I assume that is still not the intended result but leave further debugging to you.