[Rcpp-devel] R session crashes when largely using a Rcpp sourced function

Pierre.Gloaguen at ifremer.fr Pierre.Gloaguen at ifremer.fr
Sat Oct 4 01:57:01 CEST 2014


Hello everyone,

I am a new user of Rcpp, I use it to rewrite a whole R program.
at the end of the program, my "final function is the following"

#include <RcppArmadillo.h>
#include <Rcpp.h>
#define _USE_MATH_DEFINES
#include <math.h>
using namespace Rcpp;

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

List main_function_C(arma::vec X0,arma::vec XF,double t0,
                    double tF, arma::mat Gamma,
                    NumericVector piks,arma::mat muks, List Cks,
                    int max_iter = 500){
     List bounds = bounds_fun_C(piks,Cks,Gamma);
     double sup_bound = 0.5 *(as<double>(bounds["alpha_bar"])
                              + as<double>(bounds["delta_max"])
                              - as<double>(bounds["delta_min"]));
     bool accept = false;
     int kappa = 0;
     arma::mat omegas(2,2);//arbitrary size, will change in loop
     NumericVector Psi(5);//arbitrary size, will change in loop
     for(int i = 0; i < max_iter; ++i){
       kappa = rpois(1,(tF-t0)*sup_bound)[0];
       if(kappa > 0){
         Psi = stl_sort(runif(kappa,t0,tF));
         NumericVector Upsilon = runif(kappa, 0, sup_bound);
         omegas = rubb_it_C(X0, XF, t0, tF, Psi);
         NumericVector phi_omega(kappa);
         for(int i = 1; i < kappa+2; i++){
           phi_omega[i-1] = phi_C(arma::trans(omegas.row(i)),
                                  piks,muks,Cks,Gamma);
           }
         accept = all_C((phi_omega < Upsilon));
       }//end if kappa >0
       else{
         accept = true;
       }
       if(accept){
         break;
       }
     }//end for
     arma::mat res(kappa+2,3);//result matrix,
     if(kappa > 0){
       res.submat(0,0,kappa+1,1) = omegas;
       res(0,2) = t0;
       for(int i =1; i < kappa+1;++i){
        res(i,2) = Psi[i-1];
        //Psi is of size kappa the index then goes till kappa-1
       }
       res(kappa+1,2) = tF;
     }
     else{
       res(0,arma::span(0,1)) = trans(X0);
       res(1,arma::span(0,1)) = trans(XF);
       res(0,2) = t0;
       res(kappa+1,2) = tF;
     }
     return List::create(Named("skel") = res,
                         Named("accepted") = accept);
}

This function calls other funtion that I have written, all of them are  
attached on the cpp_funcs.cpp file.
The sourceCpp performs well and I can, from R obtain the following result

    sourceCpp("cpp_funcs.cpp")
    X0 <- c(0.5,0.5)
    XF <- c(1,1)
    t0 <- 0
    tF <- 1
    Gam <- diag(0.1,2)
    ncomp=2
    pis <- c(0.5,0.5)
    mu <- matrix(c(-1,1,
                    1,1),ncol=2,nrow=ncomp)
    cov <- list(diag(0.5,2),diag(1,2))
    set.seed(123)
  test <- main_function_C(X0,XF,t0,tF,Gam,pis,mu,cov)
  test
#$skel
            [,1]     [,2]      [,3]
#[1,]  0.5000000 0.500000 0.0000000
#[2,] -0.1261731 1.313880 0.4089769
#[3,]  0.5564577 1.069211 0.7883051
#[4,]  1.0000000 1.000000 1.0000000

#$accepted
#[1] TRUE

PROBLEM:

When I run the exact same code as above a large number of times, as this

    for(i in 1:1000){
    ## same code as above
    }

The r session aborts all the time, giving no error message but "R  
encountered a fatal error".
I do not encounter this problem with the intermediate functions in the  
attached file (and, of course, I'm not asking for a debugging of  
those), I only have it with this one

Is this a problem with my code? with how the loop is written?The List object?
Is this a problem of the memory?

It does it either on Rstudio or Rgui
I work on windows 7, with R version 3.0.2
Rcpp 0.11.2 and RcppArmadillo 0.4.450.1.0

I also attach a R debugging file giving parameters for all  
intermediate functions.

Sorry for the length of the function, I'm pretty sure the problem  
comes from the main function, but it sadly depends on others (which,  
as far as I know, don't pose problems)

Thanks in advance for any help,

Pierre Gloaguen
-------------- next part --------------
library(Rcpp)
library(RcppArmadillo)
library(gtools)
#At the beginning, all cpp_funcs.cpp functions should be commented
#they are tested equentially

# test all_C ---------------------------------------------------------------

#adding allC function
sourceCpp("cpp_funcs.cpp")
for(i in 1:1000){
  x <- sample(c(TRUE,FALSE),size=4,replace=T)
  test <- all_C(x)
}

# SEEMS OK

# test trace_C -------------------------------------------------------------

#adding traceC
sourceCpp("cpp_funcs.cpp")
for(i in 1:1000){
  x <- matrix(sample(100),nrow=10,ncol=10)
  test <- traceC(x)
}
#what if x is non square? If ncol > nrow, noprob
x <- matrix(sample(20),nrow=5,ncol=4)
test <- traceC(x)
#Normal error message, so that would no be this
#SEEMS OK

# test stl_sort ------------------------------------------------------------

#adding stl_sort
sourceCpp("cpp_funcs.cpp")
for(i in 1:1000){
  x <- runif(1000)
  test <- stl_sort(x)
}
#seems OK

# test Mahalanobis ------------------------------------------------------------

#adding Mahalanobis
sourceCpp("cpp_funcs.cpp")
for(i in 1:1000){
  x <- matrix(rnorm(100),ncol=2,nrow=50)
  y <- rnorm(2)
  cov <- rWishart(1,2,diag(1,2))[,,1]
  test <- Mahalanobis(x,y,cov)
}
#SEEMS OK 

# test dmvnorm_C ------------------------------------------------------------

#This function depends on Mahalanobis
#adding dmvnorm_C
sourceCpp("cpp_funcs.cpp")
for(i in 1:1000){
  x <- matrix(rnorm(100),ncol=2,nrow=50)
  y <- rnorm(2)
  cov <- rWishart(1,2,diag(1,2))[,,1]
  test <- dmvnorm_C(x,y,cov,log=FALSE)
}
#SEEMS OK

# test bounds_fun_C ------------------------------------------------------------

#This function depends on trace_C
#adding bounds_fun_C
sourceCpp("cpp_funcs.cpp")
for(i in 1:10000){
  ncomp<- sample(1:5,size=1)
  pis <- rdirichlet(1,rep(1,ncomp))[1,]
  tmp <- rWishart(ncomp,2,diag(1,2))
  cov <- lapply(1:ncomp,function(i) tmp[,,i]) 
  Gam <- rWishart(ncomp,2,diag(1,2))[,,1]
  test <- bounds_fun_C(pis,cov,Gam)
}
#SEEMS OK 

# test phi_C ------------------------------------------------------------

#This function depends on trace_C, dmvnorm_C and bounds_fun_C
#adding bounds_fun_C
sourceCpp("cpp_funcs.cpp")
for(i in 1:1000){
  ncomp<- sample(1:5,size=1)
  x <- rnorm(2)
  pis <- rdirichlet(1,rep(1,ncomp))[1,]
  tmp <- rWishart(ncomp,2,diag(1,2))
  mu <- matrix(rnorm(2*ncomp),ncol=2,nrow=ncomp)
  cov <- lapply(1:ncomp,function(i) tmp[,,i]) 
  Gam <- rWishart(ncomp,2,diag(1,2))[,,1]
  test <- phi_C(x,pis,mu,cov,Gam)
}
#SEEMS OK BUT there was an unexplained abort at some point

# test rubb_it_C ------------------------------------------------------------

#uses bounds_fun_C, stl_sort,phi_C,all_C
#adding rubb_it_C
sourceCpp("cpp_funcs.cpp")
for(i in 1:10000){
  ntime <- sample(1000,size=1)
  t0= 0
  tF = 1
  tmp_times <- sort(runif(ntime,t0,tF))
  X0 <- rnorm(2)
  XF <- rnorm(2)
  test <- rubb_it_C(X0,XF,t0,tF,tmp_times)
}
#SEEMS OK 

# test simu_ea_seg_C ------------------------------------------------------------

#uses bounds_fun_C, stl_sort,phi_C,all_C,rubb_it_C
#adding simu_ea_segC
sourceCpp("cpp_funcs.cpp")
for(i in 1:1000){
  X0 <- c(0.5,0.5)
  XF <- c(1,1)
  t0 <- 0
  tF <- 1
  Gam <- diag(0.1,2)
  ncomp=2
  pis <- c(0.5,0.5)
  mu <- matrix(c(-1,1,
                 1,1),ncol=2,nrow=ncomp)
  cov <- list(diag(0.5,2),diag(1,2))
  set.seed(123)
  test <- main_function_C(X0,XF,t0,tF,Gam,pis,mu,cov)
}
-------------- next part --------------
#include <RcppArmadillo.h>
#include <Rcpp.h>
#define _USE_MATH_DEFINES
#include <math.h>
using namespace Rcpp;

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

// [[Rcpp::export]]
bool all_C(LogicalVector x){// Equivalent to all function in R
  LogicalVector::iterator n = x.end();
  for(LogicalVector::iterator i = x.begin(); i != n; ++i){
    if(!x[*i]) return false;
  }
  return true;
}

// [[Rcpp::export]]
double trace_C(arma::mat x){// Computing the trace of a squared matrix
  int n = x.n_rows;
  double out=0.0;
  for(int i=0;i < n;++i){
    out += x(i,i);
  }
  return out;
}

// [[Rcpp::export]]
arma::vec Mahalanobis(arma::mat x, arma::rowvec center, arma::mat cov){
    int n = x.n_rows;
    arma::mat x_cen;
    x_cen.copy_size(x);
    for (int i=0; i < n; ++i) {
        x_cen.row(i) = x.row(i) - center;
    }
    return sum((x_cen * cov.i())%x_cen, 1);    
}

// [[Rcpp::export]]
arma::vec dmvnorm_C (arma::mat x,arma::rowvec mean,
                     arma::mat sigma, bool log){ 
 
    arma::vec distval = Mahalanobis(x,  mean, sigma);
 
    double logdet = sum(arma::log(arma::eig_sym(sigma)));
    double log2pi = 1.8378770664093454835606594728112352797227949472755668;
    arma::vec logretval = -( (x.n_cols * log2pi + logdet + distval)/2  ) ;
       if(log){
         return logretval;
       }
       else {
         return exp(logretval);}
}


// [[Rcpp::export]]
NumericVector stl_sort(NumericVector x) {//function to sort a Numeric vector
   NumericVector y = clone(x);
   std::sort(y.begin(), y.end());
   return y;
}

// [[Rcpp::export]]
List bounds_fun_C(NumericVector piks, List Cks, arma::mat Gamma){//intermediate function
  IntegerVector is = (seq_along(piks)-1);
  IntegerVector::iterator ncomp_it = is.end();
  int ncomp = Cks.size();
  NumericVector dets(ncomp);
  NumericVector norms_Cinv(ncomp);
  NumericVector traces(ncomp);
  NumericVector maxeigen(ncomp);
  NumericVector norms_GiCinv(ncomp);
  arma::mat Gi = arma::inv(Gamma);
  for(IntegerVector::iterator i=is.begin(); i!=ncomp_it;++i){
    arma::mat Ck = Cks[*i];
    arma::mat Cki = arma::inv(Ck);
    dets[*i] = arma::det(2*M_PI*Ck);
    norms_Cinv[*i] = arma::norm(Cki,2);
    maxeigen[*i] = arma::norm(Ck,2);
    traces[*i] = trace_C(Cki);
    norms_GiCinv[*i] = arma::norm(Gi*Cki,2);
  }
  double delta_min = -sum(piks * pow(dets,-0.5) * traces);
  double delta_max = 2*exp(-1)*sum(piks * pow(norms_Cinv,2) * pow(dets,-0.5) * maxeigen);
  double alpha_bar = exp(-1)*sum(piks * pow(norms_GiCinv,2) * maxeigen/dets);
  return List::create(Named("delta_min") = delta_min,
                      Named("delta_max") = delta_max,
                      Named("alpha_bar") = alpha_bar);
}

// [[Rcpp::export]]
double phi_C(arma::vec X, NumericVector piks,arma::mat muks,
             List Cks, arma::mat Gamma){//intermediate function
  int ncomp = piks.size();
  List tmp=bounds_fun_C(piks,Cks,Gamma);
  double delta_min = as<double>(tmp["delta_min"]);
  arma::vec GX = Gamma*X;
  arma::mat Xtmp(1,2);//to go in dmvnorm function
  Xtmp(0,0) = GX[0];
  Xtmp(0,1) = GX[1];
  arma::vec pi_gaus_dens(ncomp);
  arma::vec traces(ncomp);
  arma::vec norms_C_Ix(ncomp);
  arma::mat Gi = arma::inv(Gamma);
  arma::vec alpha(2);
  alpha.fill(0);
  double lapl_H = 0.0;
  for(int i = 0; i < ncomp;++i){
    arma::rowvec mu = muks.row(i);
    arma::mat Ck = Cks[i];
    arma::mat Cki = arma::inv(Ck);
    pi_gaus_dens[i] = piks[i]*dmvnorm_C(Xtmp,mu,Ck,false)[0];
    traces[i] = trace_C(Cki);
    norms_C_Ix[i] = pow(arma::norm(Cki*(GX-arma::trans(mu)),2),2);
    lapl_H += pi_gaus_dens[i]*(norms_C_Ix[i]-traces[i]);
    alpha += pi_gaus_dens[i]*Gi*Cki*(GX-arma::trans(mu));
  }
  double norm_alpha = pow(arma::norm(alpha,2),2);
  double res = 0.5*(norm_alpha+lapl_H-delta_min);
  return res;
}

arma::mat rubb_it_C(arma::vec X0,arma::vec XF,//Simulating a brownian bridge
                    double t0, double tF,NumericVector times) {
   double XF1 = XF[0];
   double XF2 = XF[1];
   int n = times.size();
   arma::vec X1(n+2);
   arma::vec X2(n+2);
   NumericVector ts(n+1);//times series for the loop
   ts[0] = t0;
   X1[0] = X0[0];
   X2[0] = X0[1];
   X1[n+1] = XF1;
   X2[n+1] = XF2;
   for(int i=1; i < (n+1); ++i){
     ts[i] = times[i-1];
     double tar = tF-ts[i];//time before the arrival
     double tdeb = ts[i]-ts[i-1];//time since the departure
     double ttot = tF-ts[i-1];//overall time
     X1[i] = (tar*X1[i-1]+tdeb*XF1)/ttot + rnorm(1,0,pow((tar*tdeb)/ttot,0.5))[0];
     X2[i] = (tar*X2[i-1]+tdeb*XF2)/ttot + rnorm(1,0,pow((tar*tdeb)/ttot,0.5))[0];
     }
  arma::mat res(n+2,2);
  res.col(0) = X1;
  res.col(1) = X2;
  return res;
}

// [[Rcpp::export]]
List main_function_C(arma::vec X0,arma::vec XF,double t0,
                   double tF, arma::mat Gamma,
                   NumericVector piks,arma::mat muks, List Cks,
                   int max_iter = 500){
    List bounds = bounds_fun_C(piks,Cks,Gamma);
    double sup_bound = 0.5 *(as<double>(bounds["alpha_bar"])
                             + as<double>(bounds["delta_max"])
                             - as<double>(bounds["delta_min"]));
    bool accept = false;
    int kappa = 0;
    arma::mat omegas(2,2);//arbitrary size, will change in loop
    NumericVector Psi(5);//arbitrary size, will change in loop
    for(int i = 0; i < max_iter; ++i){
      kappa = rpois(1,(tF-t0)*sup_bound)[0];
      if(kappa > 0){
        Psi = stl_sort(runif(kappa,t0,tF));
        NumericVector Upsilon = runif(kappa, 0, sup_bound);
        omegas = rubb_it_C(X0, XF, t0, tF, Psi);
        NumericVector phi_omega(kappa);
        for(int i = 1; i < kappa+2; i++){
          phi_omega[i-1] = phi_C(arma::trans(omegas.row(i)),
                                 piks,muks,Cks,Gamma);
          }
        accept = all_C((phi_omega < Upsilon));
      }
      else{
        accept = true;
      }
      if(accept){
        break;
      }
    }
    arma::mat res(kappa+2,3);//result matrix, 
    if(kappa > 0){
      res.submat(0,0,kappa+1,1) = omegas;
      res(0,2) = t0;
      for(int i =1; i < kappa+1;++i){
       res(i,2) = Psi[i-1];
       //Psi is of size kappa the index then goes till kappa-1
      }
      res(kappa+1,2) = tF;
    }
    else{
      res(0,arma::span(0,1)) = trans(X0);
      res(1,arma::span(0,1)) = trans(XF);
      res(0,2) = t0;
      res(kappa+1,2) = tF;
    }
    return List::create(Named("skel") = res,
                        Named("accepted") = accept);
}


More information about the Rcpp-devel mailing list