## (nombre de grains de pollen / cm2) x (2 cm2 ; surface de la portion de feuille) x (40 répétitions)
tox.dose=c(0,23,51,104,205,288)
tox=tox.dose*2*40

## proportion de mortalité
mortality=c(42, 70, 45, 101, 269, 340)/400
names(mortality)=paste("dose",tox.dose,sep="_")

## display data
plot(tox.dose,mortality,ylim=c(0,1), 
     xlab="Toxic dose (# toxic pollen grains / cm2)",
     ylab="Death proportion")

library(MCMCpack) ## load this package for drawing data from the Dirichlet distribution
library(abc) ## load this package for running functions related to ABC
library(psych) ## load this package for generating plots with functions pairs.panels() and violinBy()

par(mfrow=c(2,2))
for(b in c(0.01,0.1,1,4)){
    hist(rdirichlet(n=40,alpha=rep(1/b,20)),xlab="Proportion",xlim=c(0,1),main=paste("beta =",b))
    abline(v=1/20,col=2)
}

a2=0.1## try both values 0.1 and 0.02
tox.seq=0:100
par(mfrow=c(2,2))
for(a3 in c(0.5,1,2,10)){
    plot(tox.seq,exp(-(a2*tox.seq)^a3),type="l",ylim=c(0,1),
         xlab="Individual internal toxic level",
         ylab="Relative daily survival probability",
         main=paste("alpha2 =",a2,"; alpha3 =",a3))
    points(tox/400,exp(-(a2*tox/400)^a3),pch=19,col=2)
}

## simulator
mortality.pop=function(n,nj,X,kout,alpha1,alpha2,alpha3,beta, ngroup=1){
    death.prop=rep(0,length(X))
    size=n/ngroup
    for(i in 1:length(X)){
        prop.feed=rdirichlet(n=ngroup,alpha=rep(1/beta,size*2))
        Q=cbind(as.numeric(prop.feed[,1:size]),as.numeric(prop.feed[,size+1:size]))
        rho=cbind(X[i]*Q[,1])%*%((1-kout)^((1:nj)-1)) +
            cbind(X[i]*Q[,2])%*%c(0,(1-kout)^((2:nj)-2))
        event=ceiling(alpha1*exp(-(alpha2*rho)^alpha3)-matrix(runif(nj*n),n,nj))
        surv=apply(event,1,min)
        death.prop[i]=mean(1-surv)
    }
    return(death.prop)
}

## run and show 1 simulation ; play with alpha2 (0.1; 0.01; 0.001)
mortality.sim=mortality.pop(n=400,nj=7,tox, kout=0.5, alpha1= 0.99, alpha2=0.01, 
    alpha3=4, beta=0.25, ngroup=40)
cat("Doses toxiques:\n")
print(tox.dose)
cat("Taux de mortalité observés aux différentes doses toxiques:\n")
print(mortality.sim)
plot(tox.dose,mortality.sim,ylim=c(0,1), 
     xlab="Toxic dose (# toxic pollen grains / cm2)",
     ylab="Death proportion")
points(tox.dose,mortality,col=2)  ## death proportions observed for Inachis io

## plot priors
par(mfrow=c(2,3))
plot(seq(-5,10,0.01),dunif(seq(-5,10,0.01),0,4),type="l",xlab=expression(beta),ylab="Density")
plot(seq(-2,3,0.01),dunif(seq(-2,3,0.01),0,1),type="l",xlab=expression(k[out]),ylab="Density")
plot(seq(0,1.1,0.01),dbeta(seq(0,1.1,0.01),shape1=50,shape2=1),type="l",xlab=expression(alpha[1]),ylab="Density")
plot(seq(0,0.01,0.0001),dexp(seq(0,0.01,0.0001),rate=1000),type="l",xlab=expression(alpha[2]),ylab="Density")
plot(seq(-5,15,0.01),seq(-5,15,0.01)==1,type="l",xlab=expression(alpha[3]),ylab="Density",
    main="Fixed alpha3")
plot(seq(-5,15,0.01),dunif(seq(-5,15,0.01),0,10),type="l",xlab=expression(alpha[3]),ylab="Density",
    main="Variable alpha3")

## run simulations for ABC (try first with alpha3=1)
time0=proc.time()
nsimul=10^4
param=cbind(runif(nsimul,0,1), ## kout
    rbeta(nsimul,shape1=50,shape2=1), ## alpha1
    rexp(nsimul,rate=1000), ## alpha2
    rep(1,nsimul), ## alpha3=1
    runif(nsimul,0,4)) ## beta
colnames(param)=c("kout","alpha1","alpha2","alpha3","beta")
stats=t(apply(param,1,function(u) mortality.pop(n=400, nj=7, tox, u[1],u[2],u[3],u[4],u[5],ngroup=40)))
colnames(stats)=names(mortality)
proc.time()-time0

param[1:10,]
stats[1:10,]

## show simulations
nsimul0=min(10^4,nsimul)
pairs.panels(param[1:nsimul0,],pch=".")
pairs.panels(rbind(stats[1:nsimul0,],mortality),pch=c(rep(19,nsimul0),21),bg="green")

## run simulations for ABC (try with alpha3 a priori uniformly distributed in [0,10])
time0=proc.time()
nsimul=10^5
param=cbind(runif(nsimul,0,1), ## kout
    rbeta(nsimul,shape1=50,shape2=1), ## alpha1
    rexp(nsimul,rate=1000), ## alpha2
    runif(nsimul,0,10), ## alpha3
    runif(nsimul,0,4)) ## beta
colnames(param)=c("kout","alpha1","alpha2","alpha3","beta")
stats=t(apply(param,1,function(u) mortality.pop(n=400, nj=7, tox, u[1],u[2],u[3],u[4],u[5],ngroup=40)))
colnames(stats)=names(mortality)
proc.time()-time0
              


## save param and stats as rds files                        
# saveRDS("save/param.rds")
# saveRDS("save/stats.rds") 
              
## load param and stats from existing files              
# param=readRDS("save/param.rds")
# stats=readRDS("save/stats.rds") 
# nsimul=nrow(stats)

param[1:10,]
stats[1:10,]

## show simulations
nsimul0=min(10^4,nsimul)
pairs.panels(param[1:nsimul0,],pch=".")
pairs.panels(rbind(stats[1:nsimul0,],mortality),pch=c(rep(19,nsimul0),21),bg="green")

## apply ABC (only print output)
rej0=abc(target=mortality, param=param, sumstat=stats, tol=0.001, method="rejection")
rej0
summary(rej0)
## param[rej0$region,] ## table with all retained parameters

## but what tolerance value?
time0=proc.time()
tol.seq=c(0.0001,0.0005,0.001,0.005,0.01)
rej.cv=cv4abc(param=param, sumstat=stats, abc.out=rej0, nval=200, tols=tol.seq)
proc.time()-time0
summary(rej.cv)

## save rej.cv as rds files   
# saveRDS(tol.seq,"save/tol.seq.rds")
# saveRDS(rej.cv,"save/rej.cv.rds")
              
## load param and stats from existing files              
# tol.seq=readRDS("save/tol.seq.rds")
# tol.seq
# rej.cv0=readRDS("save/rej.cv.rds")
# summary(rej.cv0)

## display mean error assessed by cross-validation over each parameter 
summ.rej.cv=summary(rej.cv)
par(mfrow=c(2,3))
apply(summ.rej.cv,2,function(u) plot(tol.seq,u,type="b",log="x"))
mean.error=rowMeans(summ.rej.cv)
plot(tol.seq,mean.error,type="b",log="x")
tol.opt=tol.seq[mean.error==min(mean.error)]

cat("Tolérance sélectionnée par validation croisée:\n")
as.numeric(tol.opt)

plot(rej.cv)

## apply ABC 
rej=abc(target=mortality, param=param, sumstat=stats, tol=tol.opt, method="rejection")
rej
summary(rej)
par(mfrow=c(2,3))
hist(rej)

pairs.panels(rbind(param[1:nsimul0,],param[rej$region,]),pch=c(rep(19,nsimul0),rep(21,sum(rej$region))),
             bg="green")

pairs.panels(param[rej$region,],pch=21,bg="green")

pairs.panels(rbind(stats[1:nsimul0,],mortality),pch=c(rep(19,nsimul0),21),
             bg="green")

pairs.panels(rbind(stats[rej$region,],mortality),pch=c(rep(19,sum(rej$region)),21),
             bg="green")

## Marginal posterior distributions of mortality proportions
median.mortality=apply(stats,2,median)
median.mortality
cols.index=round((max(median.mortality)-median.mortality)/(max(median.mortality)-min(median.mortality))*100)
violinBy(stats,col=heat.colors(100)[pmax(1,cols.index)],
    xlab="Toxic dose (# toxic pollen grains / cm2)",
         ylab="Death proportion at the 7th day")
points(mortality,pch=19)