## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    cache = TRUE,
    collapse = TRUE,
    comment = "#>"
)

## ----eval=FALSE---------------------------------------------------------------
# install.packages("EasyABC")

## ----setup, include=FALSE-----------------------------------------------------
library(EasyABC)

## ----eval=FALSE---------------------------------------------------------------
# help(package="EasyABC")

## ----eval=FALSE---------------------------------------------------------------
# help(ABC_sequential)

## ----simple_prior-------------------------------------------------------------
my_prior=list(c("unif",0,1),c("normal",1,2))

## ----custom_prior-------------------------------------------------------------
my_prior=list(list(c("runif",1,0,1), c("dunif",0,1)))

## ----eval=FALSE---------------------------------------------------------------
# prior = list(c("unif",0,1),c("unif",0,10))
# ABC_rejection(model=a_model,prior=prior,nb_simul=3, prior_test="X2 > X1")

## ----eval=FALSE---------------------------------------------------------------
# extern "C" {
#     void trait_model(double *input,double *stat_to_return){
#         // compute output and fill the array stat_to_return
#     }
# }

## ----eval=FALSE---------------------------------------------------------------
# trait_model <- function(input=c(1,1,1,1,1,1)) {
#     .C("trait_model",input=input,stat_to_return=array(0,4))$stat_to_return
# }

## ----eval=FALSE---------------------------------------------------------------
# trait_model <- function(input=c(1,1,1,1,1,1)) {
#     .C("trait_model",input=c(input[1], 500, input[2:3], 1, input[4:5]),
#     stat_to_return=array(0,4))$stat_to_return
# }

## ----eval=FALSE---------------------------------------------------------------
# r<-read.table('input',head=F)
# sink('mod.input')
# cat(paste('1','p1','unif',round(r[1,],0),round(r[1,],0),sep='\t'))
# cat('\n')
# cat(paste('1','p2','unif',round(r[2,],0),round(r[2,],0),sep='\t'))
# cat('\n')
# cat(paste('1','p3','unif',round(r[3,],0),round(r[3,],0),sep='\t'))
# sink()

## ----eval=FALSE---------------------------------------------------------------
# prior=list(c("unif",500,1000),c("unif",100,500),c("unif",50,200))
# ABC_sim<-ABC_rejection(model=binary_model('./run_sim.sh'),prior=prior,nb_simul=3)

## ----eval=FALSE---------------------------------------------------------------
# mymodel <- function(x) {
#     library("rJava")
#     .jinit(classpath=".")
#     result = .jcall(J("Model"),"[D","run",.jarray(x))
#     result
# }

## ----eval=FALSE---------------------------------------------------------------
# prior=list(c("unif",0,1),c("normal",1,2))
# ABC_sim<-ABC_rejection(model=mymodel,prior=prior,nb_simul=3)

## ----toy_model----------------------------------------------------------------
toy_model<-function(x){
    c( x[1] + x[2] + rnorm(1,0,0.1) , x[1] * x[2] + rnorm(1,0,0.1) ) 
}

## ----toy_prior----------------------------------------------------------------
toy_prior=list(c("unif",0,1),c("normal",1,2))

## ----sum_stat_obs-------------------------------------------------------------
sum_stat_obs=c(1.5,0.5)

## ----ABC_rejection------------------------------------------------------------
set.seed(1)
n=10
p=0.2
ABC_rej<-ABC_rejection(model=toy_model, prior=toy_prior, nb_simul=n,
summary_stat_target=sum_stat_obs, tol=p)
ABC_rej

## ----ABC_rejection2-----------------------------------------------------------
# Run the ABC rejection on the model
set.seed(1)
n=10
ABC_rej<-ABC_rejection(model=toy_model, prior=toy_prior, nb_simul=n)
ABC_rej

## ----abcinstall,eval=FALSE----------------------------------------------------
# # Install if needed the "abc" package
# install.packages("abc")

## ----abc----------------------------------------------------------------------
# Post-process the simulations outputs
library(abc)
rej<-abc(sum_stat_obs, ABC_rej$param, ABC_rej$stats, tol=0.2, method="rejection")
# simulations selected:
rej$unadj.values
# their associated summary statistics:
rej$ss
# their normalized euclidean distance to the data summary statistics:
rej$dist

## ----ABC_Beaumont-------------------------------------------------------------
n=10
tolerance=c(1.25,0.75)
ABC_Beaumont<-ABC_sequential(method="Beaumont", model=toy_model,
prior=toy_prior, nb_simul=n, summary_stat_target=sum_stat_obs,
tolerance_tab=tolerance)
ABC_Beaumont

## ----ABC_Drovandi-------------------------------------------------------------
n=10
tolerance=0.75
c_drov=0.7
ABC_Drovandi<-ABC_sequential(method="Drovandi", model=toy_model,
prior=toy_prior, nb_simul=n, summary_stat_target=sum_stat_obs,
tolerance_tab=tolerance, c=c_drov)
ABC_Drovandi

## ----ABC_Delmoral-------------------------------------------------------------
n=10
alpha_delmo=0.5
tolerance=0.75
ABC_Delmoral<-ABC_sequential(method="Delmoral", model=toy_model,
prior=toy_prior, nb_simul=n, summary_stat_target=sum_stat_obs,
alpha=alpha_delmo, tolerance_target=tolerance)
ABC_Delmoral

## ----toy_prior2---------------------------------------------------------------
toy_prior2=list(c("unif",0,1),c("unif",0.5,1.5))

## ----ABC_Lenormand------------------------------------------------------------
n=10
pacc=0.4
ABC_Lenormand<-ABC_sequential(method="Lenormand", model=toy_model,
prior=toy_prior2, nb_simul=10, summary_stat_target=sum_stat_obs,
p_acc_min=pacc)
ABC_Lenormand

## ----ABC_Marjoram_original----------------------------------------------------
n=10
ABC_Marjoram_original<-ABC_mcmc(method="Marjoram_original", model=toy_model,
prior=toy_prior, summary_stat_target=sum_stat_obs, n_rec=n)
ABC_Marjoram_original

## ----ABC_Marjoram-------------------------------------------------------------
n=10
ABC_Marjoram<-ABC_mcmc(method="Marjoram", model=toy_model,
prior=toy_prior, summary_stat_target=sum_stat_obs, n_rec=n)
ABC_Marjoram

## ----ABC_Wegmann--------------------------------------------------------------
n=10
ABC_Wegmann<-ABC_mcmc(method="Wegmann", model=toy_model,
prior=toy_prior, summary_stat_target=sum_stat_obs, n_rec=n)
ABC_Wegmann

## ----SABCPrior----------------------------------------------------------------
r.prior <- function()   c(runif(1,0,1),rnorm(1,1,2))
d.prior <- function(x)  dunif(x[1],0,1)*dnorm(x[2],1,2)

## ----SABCParam----------------------------------------------------------------
n.sample <- 300
iter.max <- n.sample * 30
eps.init <- 2

## ----SABC, include=FALSE------------------------------------------------------
ABC_Albert <-SABC(r.model = toy_model,
r.prior = r.prior,
d.prior = d.prior,
n.sample = n.sample,
eps.init = eps.init,
iter.max = iter.max,
method = "informative",
y = sum_stat_obs
)

## ----SABCPlot1----------------------------------------------------------------
plot(ABC_Albert$E[,1:2])

## ----toy_model_parallel-------------------------------------------------------
toy_model_parallel<-function(x){
set.seed(x[1]) # so that each core is initialized with a different seed value.
c( x[2] + x[3] + rnorm(1,0,0.1) , x[2] * x[3] + rnorm(1,0,0.1) ) 
}

## ----ABC_rejection3-----------------------------------------------------------
set.seed(1)
n=10
p=0.2
ABC_rej<-ABC_rejection(model=toy_model_parallel, prior=toy_prior,
nb_simul=n, summary_stat_target=sum_stat_obs, tol=p, n_cluster=2,
use_seed=TRUE)
ABC_rej

## ----trait_prior--------------------------------------------------------------
trait_prior=list(c("unif",3,5),c("unif",-2.3,1.6),
c("unif",-25,125), c("unif",-0.7,3.2))

## ----sum_stat_obs2------------------------------------------------------------
sum_stat_obs=c(100,2.5,20,30000)

## ----ABC_rejection4-----------------------------------------------------------
set.seed(1)
n=10
p=0.2
ABC_rej<-ABC_rejection(model=trait_model, prior=trait_prior, nb_simul=n,
summary_stat_target=sum_stat_obs, tol=p, use_seed=TRUE)
ABC_rej

## ----abcinstall2,eval=FALSE---------------------------------------------------
# install.packages("abc")

## ----abc2---------------------------------------------------------------------
library(abc)
set.seed(1)
n=10
p=0.2
ABC_rej<-ABC_rejection(model=trait_model, prior=trait_prior, nb_simul=n, use_seed=TRUE)
rej<-abc(sum_stat_obs, ABC_rej$param, ABC_rej$stats,
tol=0.2, method="rejection")
# simulations selected:
rej$unadj.values
# their associated summary statistics:
rej$ss
# their normalized euclidean distance to the data summary statistics:
rej$dist

## ----ABC_Beaumont2------------------------------------------------------------
n=10
tolerance=c(8,5)
ABC_Beaumont<-ABC_sequential(method="Beaumont", model=trait_model,
prior=trait_prior, nb_simul=n, summary_stat_target=sum_stat_obs,
tolerance_tab=tolerance, use_seed=TRUE)
ABC_Beaumont

## ----ABC_Drovandi2------------------------------------------------------------
n=10
tolerance=3
c_drov=0.7
ABC_Drovandi<-ABC_sequential(method="Drovandi", model=trait_model,
prior=trait_prior, nb_simul=n, summary_stat_target=sum_stat_obs,
tolerance_tab=tolerance, c=c_drov, use_seed=TRUE)
ABC_Drovandi

## ----ABC_Delmoral2------------------------------------------------------------
n=10
alpha_delmo=0.5
tolerance=3
ABC_Delmoral<-ABC_sequential(method="Delmoral", model=trait_model,
prior=trait_prior, nb_simul=n, summary_stat_target=sum_stat_obs,
alpha=alpha_delmo, tolerance_target=tolerance, use_seed=TRUE)
ABC_Delmoral

## ----ABC_Lenormand2-----------------------------------------------------------
n=10
pacc=0.4
ABC_Lenormand<-ABC_sequential(method="Lenormand", model=trait_model,
prior=trait_prior, nb_simul=n, summary_stat_target=sum_stat_obs,
p_acc_min=pacc, use_seed=TRUE)
ABC_Lenormand

## ----ABC_Marjoram_original2---------------------------------------------------
n=10
ABC_Marjoram_original<-ABC_mcmc(method="Marjoram_original", model=trait_model,
prior=trait_prior, summary_stat_target=sum_stat_obs, n_rec=n, use_seed=TRUE)
ABC_Marjoram_original

## ----ABC_Marjoram2------------------------------------------------------------
n=10
n_calib=10
tol_quant=0.2 
ABC_Marjoram<-ABC_mcmc(method="Marjoram", model=trait_model, prior=trait_prior,
summary_stat_target=sum_stat_obs,
n_rec=n, n_calibration=n_calib, tolerance_quantile=tol_quant, use_seed=TRUE)
ABC_Marjoram

## ----ABC_Wegmann2-------------------------------------------------------------
n=10
n_calib=10
tol_quant=0.2 
ABC_Wegmann<-ABC_mcmc(method="Wegmann", model=trait_model, prior=trait_prior,
summary_stat_target=sum_stat_obs,
n_rec=n, n_calibration=n_calib, tolerance_quantile=tol_quant, use_seed=TRUE)
ABC_Wegmann

## ----SABCPrior2---------------------------------------------------------------
r.prior <- function() c(runif(1,3,5), runif(1,-2.3,1.6),
    runif(1,-25,125),runif(1,-0.7,3.2),1)
d.prior <- function(x) dunif(x[1],3,5) * dunif(x[2],-2.3,1.6) *
    dunif(x[3],-25,125) * dunif(x[4],-0.7,3.2)

## ----SABCParam2---------------------------------------------------------------
n.sample <- 30
iter.max <- n.sample * 30
eps.init <- 20

## ----SABC2, include=FALSE-----------------------------------------------------
ABC_Albert <-SABC(r.model = trait_model,
    r.prior = r.prior,
    d.prior = d.prior,
    n.sample = n.sample,
    eps.init = eps.init,
    iter.max = iter.max,
    method = "noninformative",
    y = sum_stat_obs
    )

## ----SABCPlot2----------------------------------------------------------------
hist(ABC_Albert$E[,3],breaks=100)