[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