[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