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

Pierre.Gloaguen at ifremer.fr Pierre.Gloaguen at ifremer.fr
Wed Oct 8 01:20:01 CEST 2014


Hello Dirk, and everyone,

I sadly was too optimistic thinking that a gc() would solve the problem.
The main_function_C I wrote was putted in a more general function  
named simu_cond_ea_C

List simu_cond_ea_C(arma::vec X0,arma::vec XF,
                     double t0, double tF, NumericVector times,
                     arma::mat Gamma, NumericVector piks,
                     arma::mat muks,List Cks,
                     int max_iter = 500, bool do_Lamp=true){
   arma::vec X0_tmp = X0;
   arma::vec XF_tmp = XF;
   if(do_Lamp){
     arma::mat Gi = arma::inv(Gamma);
     X0_tmp = Gi*X0_tmp;
     XF_tmp = Gi*XF_tmp;
   }
   List tmp = main_function_C(X0,XF,t0,tF,Gamma,piks,muks,Cks,max_iter);
   arma::mat skeleton = tmp["skel"];
   bool accept = as<bool>(tmp["accepted"]);
   arma::vec tmp_vec = skeleton.col(2);
   NumericVector skel_times(tmp_vec.begin(),tmp_vec.end());
   IntegerVector segmentation = findInterval_C(times,skel_times);
   arma::mat res(skeleton);
   IntegerVector idx_seg= unique(segmentation);
   IntegerVector::iterator end_loop = idx_seg.end();
   for(IntegerVector::iterator it=idx_seg.begin();it!=end_loop;++it){
     NumericVector times_tmp = times[segmentation==*it];
     arma::mat tmp_bb = rubb_it_C(arma::trans(skeleton(*it-1,arma::span(0,1))),
                                  arma::trans(skeleton(*it,arma::span(0,1))),
                                  skel_times[*it-1],
                                  skel_times[*it],
                                  times_tmp);
     int n_tmp = times_tmp.size();
     arma::mat tmp_res(n_tmp,3);
     tmp_res.col(2) = as<arma::vec>(times_tmp);
     tmp_res.submat(0,0,n_tmp-1,1) = tmp_bb.submat(1,0,n_tmp,1);
     res = rbind_C(res,tmp_res);
   }
   arma::mat out(res);
   arma::uvec indx_time = sort_index(res.col(2));
   int nobs = res.n_rows;
   for(int i=0;i < nobs;++i){
     out.row(i) = res.row(indx_time[i]);
   }
   return List::create(Named("skel") = out,
                         Named("accepted") = accept);

}

(again all the functions needed for its compilation are attached).
The loop
for(i in 1:10000){
  test <- simu_cond_ea_C(myparameters)
  gc()
}
seems to work, but it somehow leads to an instable R session.
And, moreover, when I embed it in my Rprogram in a function (where  
there is a loop also):
simu_mc_xus_R <- function(X,theta0,Gamma,N,do_Lamp=F){
     t0.glob <- X[1,3]
     tF.glob <- X[nrow(X),3]
     us <- stl_sort(runif(N,t0.glob,tF.glob))
     which.seg <- findInterval_C(us,X[,3])
     xus_tmp <- matrix(nrow=N,ncol=3)
     xus_tmp[,3] <- us
     for(i in 1:N){
       X0 <- X[which.seg[i],1:2]
       XF <- X[which.seg[i]+1,1:2]
       t0 <- X[which.seg[i],3]
       tF <- X[which.seg[i]+1,3]
       tmp <- simu_cond_ea_C(X0,XF,t0=t0,tF=tF,times=us[i],
                             Gamma=Gamma,
                             Cks=theta0$Cks,
                             muks=theta0$muks,
                             piks=theta0$piks,do_Lamp=do_Lamp)$skel
       xus_tmp[i,1:2] <- tmp[tmp[,3]==us[i],1:2]
       gc()
     }
     xus_tmp
   }

then I encounter a problem for a big N (10000 for instance), it either
1) crash the R session with no message but "r encountered a fatal error"
2) or print this message
Error in gc() :
   GC encountered a node (0x000000000ce11a10) with an unknown SEXP  
type: <unknown> at memory.c:1636

I don't have much experience in dealing with memory, so I don't know  
where the problem could come from. My cpp code? R intern memory problem?
It's worth noting the object that I stock are not big, I wrote this  
program before in R, it was slower but it never failed for memory  
reasons.

I run on R 3.1.1, platform  x86_64-w64-mingw32/x64 (64-bit)
I work on windows 7.
I attach the cpp_functions, the R_programm that embeds the last  
function, and a R file debug_test wich gives R code to test the  
functions and lead to my problems.

Any clue about the origin of this problem would be welcomed!

Pierre Gloaguen

Pierre.Gloaguen at ifremer.fr a écrit :

> Dirk Eddelbuettel <edd at debian.org> a écrit :
>
>>
>> On 6 October 2014 at 23:36, Pierre.Gloaguen at ifremer.fr wrote:
>> | Hello,
>> |
>> | I have found out that the problem was in the R loop, the garbage
>> | collection of R wasn't perform efficiently. Indeed, when I force the
>> | garbage collection to be done, using R function gc(), the Rsession
>> | won't crash, although the execution of the loop will be slower.
>> | This leads me to another question. Is there anyway to force the
>> | garbage collection inside a Rcpp function?
>>
>> Yes, the Rcpp::Function() class allows you to call R functions from C++.
>>
> Ok, i will try this, thank you for the help!
>> Part of me thinks, though, that once you are in C++ you may be able to just
>> control your memory, not call back, compute your result and then report it
>> back to R.  But such an idealised situation may not work in your case -- and
>> sorry that I did not have time to work through your code in any detail.
>
> No worries! I'm glad I found out by myself after all this times
> Have a nice day,
>
> Pierre
>>
>> Dirk
>>
>>
>> | like this
>> | NumericVector myfunction(NumericVector x){
>> |     for(int= i = 0; i < x.size; ++i){
>> | //do my stuff;
>> | gc();
>> | }
>> | }
>> |
>> | Thanks for any help!
>> |
>> | Pierre Gloaguen
>> |
>> | pgloague at ifremer.fr a écrit :
>> |
>> | > 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
>> | >
>> |
>> |
>> |
>> | _______________________________________________
>> | Rcpp-devel mailing list
>> | Rcpp-devel at lists.r-forge.r-project.org
>> | https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>>
>> --
>> http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org
>>
>
>
>
>

-------------- 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(!*i) return false;
  }
  return true;
}

// [[Rcpp::export]]
arma::mat rbind_C(arma::mat x, arma::mat y){// rbind two matrices
  int nx = x.n_rows;
  int ny = y.n_rows;
  int ncol = x.n_cols;
  arma::mat out(nx+ny,ncol);
  out.submat(0,0,nx-1,ncol-1)=x;
  out.submat(nx,0,nx+ny-1,ncol-1)=y;
  return out;
}

// [[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]]
IntegerVector findInterval_C(NumericVector x, NumericVector breaks) {// equivalent of findIntervalfunction
  IntegerVector out(x.size());
  NumericVector::iterator it, pos;
  IntegerVector::iterator out_it;
  for(it = x.begin(), out_it = out.begin(); it != x.end(); 
      ++it, ++out_it) {
    pos = std::upper_bound(breaks.begin(), breaks.end(), *it);
    *out_it = std::distance(breaks.begin(), pos);
  }

  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);
}
// [[Rcpp::export]]
List simu_cond_ea_C(arma::vec X0,arma::vec XF,
                    double t0, double tF, NumericVector times,
                    arma::mat Gamma, NumericVector piks,
                    arma::mat muks,List Cks,
                    int max_iter = 500, bool do_Lamp=true){
  arma::vec X0_tmp = X0;
  arma::vec XF_tmp = XF;
  if(do_Lamp){
    arma::mat Gi = arma::inv(Gamma);
    X0_tmp = Gi*X0_tmp;
    XF_tmp = Gi*XF_tmp;
  }
  List tmp = main_function_C(X0,XF,t0,tF,Gamma,piks,muks,Cks,max_iter);
  arma::mat skeleton = tmp["skel"];
  bool accept = as<bool>(tmp["accepted"]);
  arma::vec tmp_vec = skeleton.col(2);
  NumericVector skel_times(tmp_vec.begin(),tmp_vec.end());
  IntegerVector segmentation = findInterval_C(times,skel_times);
  arma::mat res(skeleton);
  IntegerVector idx_seg= unique(segmentation);
  IntegerVector::iterator end_loop = idx_seg.end();
  for(IntegerVector::iterator it=idx_seg.begin();it!=end_loop;++it){
    NumericVector times_tmp = times[segmentation==*it]; 
    arma::mat tmp_bb = rubb_it_C(arma::trans(skeleton(*it-1,arma::span(0,1))),
                                 arma::trans(skeleton(*it,arma::span(0,1))),
                                 skel_times[*it-1],
                                 skel_times[*it],
                                 times_tmp);
    int n_tmp = times_tmp.size();
    arma::mat tmp_res(n_tmp,3);
    tmp_res.col(2) = as<arma::vec>(times_tmp);
    tmp_res.submat(0,0,n_tmp-1,1) = tmp_bb.submat(1,0,n_tmp,1);
    res = rbind_C(res,tmp_res);
  }
  arma::mat out(res);
  arma::uvec indx_time = sort_index(res.col(2));
  int nobs = res.n_rows;
  for(int i=0;i < nobs;++i){
    out.row(i) = res.row(indx_time[i]);
  }
  return List::create(Named("skel") = out,
                        Named("accepted") = accept);
  
}

//// [[Rcpp::export]]
//arma::mat simu_mc_xusC(arma::mat X, List theta0,arma::mat Gamma,
//                        int Nsamp, bool do_Lamp=false){
//  int nobs = X.n_rows;
//  double t0_glob = X(0,2);
//  double tF_glob = X(nobs-1,2);
//  NumericVector piks = theta0["piks"];
//  arma::mat muks = theta0["muks"];
//  List Cks = theta0["Cks"];
//  NumericVector us = stl_sort(runif(Nsamp,t0_glob,tF_glob));
//  arma::vec tmp = X.col(2);
//  NumericVector times_obs(tmp.begin(),tmp.end());
//  IntegerVector which_seg = findInterval_C(us,times_obs);
//  arma::mat res(Nsamp,3);
//  arma::mat Xtrans = arma::trans(X);
//  for(int i = 0; i <  Nsamp; ++i){
//    arma::vec X0 = Xtrans(arma::span(0,1),which_seg[i]-1).col(0);
//    arma::vec XF = Xtrans(arma::span(0,1),which_seg[i]).col(0);
//    double t0 =  X(which_seg[i]-1,2);
//    double tF = X(which_seg[i],2);
//    NumericVector times_tp(1);
//    times_tp.fill(us[i]);
//    List simu = main_function_C(X0,XF,t0,tF, times_tp,
//                                Gamma, piks, muks, Cks,
//                                500, do_Lamp);
//    bool accept = simu["accepted"];
//    arma::mat res_tmp = simu["skel"];
//    if(!accept){
//      std::cout << "Conditionnal simu failed for point" <<i+1<<", increase max_iter" << std::endl;
//    }
//    res.row(i) = res_tmp.row(1);
//  }
//  return res;
//}
-------------- 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")
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))
for(i in 1:1000){
  set.seed(123)
  test <- main_function_C(X0,XF,t0,tF,Gam,pis,mu,cov)
  gc()
}

# test simu_cond_ea ------------------------------------------------------------

#uses bounds_fun_C, stl_sort,phi_C,all_C,rubb_it_C
#adding simu_cond_ea
sourceCpp("cpp_funcs.cpp")
X0 <- c(0.5,0.5)
XF <- c(1,1)
t0 <- 0
tF <- 1
tim <- seq(t0+0.1,tF-0.1,0.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))
for(i in 1:10000){
  test <- simu_cond_ea_C(X0,XF,t0,tF,tim,Gam,pis,mu,cov) 
  gc()
  print(i)
}


# test simu_mc_xus.R ------------------------------------------------------
sourceCpp("cpp_funcs.cpp")
source("R_function_embedding_rcpp.R")
ncomp=2
X <- cbind(rnorm(10),rnorm(10),seq(0,1,length.out=10))
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))
theta0 <- list(piks=pis,muks=mu,Cks=cov)
Gam <- diag(0.1,2)
N <- 10000
test <- simu_mc_xus_R(X,theta0,Gam,N,FALSE)
shine(test)
-------------- next part --------------
simu_mc_xus_R <- function(X,theta0,Gamma,N,do_Lamp=F){
    t0.glob <- X[1,3]
    tF.glob <- X[nrow(X),3]
    us <- stl_sort(runif(N,t0.glob,tF.glob))
    which.seg <- findInterval_C(us,X[,3])
    xus_tmp <- matrix(nrow=N,ncol=3)
    xus_tmp[,3] <- us
    for(i in 1:N){
      X0 <- X[which.seg[i],1:2]
      XF <- X[which.seg[i]+1,1:2]
      t0 <- X[which.seg[i],3]
      tF <- X[which.seg[i]+1,3]
      tmp <- simu_cond_ea_C(X0,XF,t0=t0,tF=tF,times=us[i],
                            Gamma=Gamma,
                            Cks=theta0$Cks,
                            muks=theta0$muks,
                            piks=theta0$piks,do_Lamp=do_Lamp)$skel
      xus_tmp[i,1:2] <- tmp[tmp[,3]==us[i],1:2]
      gc()
    }
    xus_tmp
  }


More information about the Rcpp-devel mailing list