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.gloag...@ifremer.fr a écrit :
Dirk Eddelbuettel <e...@debian.org> a écrit :
On 6 October 2014 at 23:36, pierre.gloag...@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
|
| pgloa...@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@lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
--
http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org
#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;
//}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)
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
}_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
