rm(list = ls())
require(deSolve)
## Loading required package: deSolve
require(lhs)
## Loading required package: lhs
## Warning: package 'lhs' was built under R version 4.0.2
legrand <- function (t, y, params) {
  
S.h <- y[1] #create local variable S, first element of y
E.h <- y[2] 
I.h <- y[3]
S.l <- y[4]
E.l <- y[5]
I.l <- y[6]
Q <- y[7]
R <- y[8]
V <- y[9]

with(
  as.list(params),
  {
dS.h<--q*beta*(I.h+I.l)*S.h/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)-c*(1-exp((-1/K)*V))*S.h
dE.h<-q*beta*(I.h+I.l)*S.h/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)+c*(1-exp((-1/K)*V))*S.h-lambda*E.h
dI.h<-lambda*E.h-b*g*I.h-aH*h*(1-g)*I.h-gammah*(1-h)*(1-g)*I.h
dS.l<--(1-p)*q*beta*(I.h+I.l)*S.l/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)-(1-p)*c*(1-exp((-1/K)*V))*S.l
dE.l<-(1-p)*q*beta*(I.h+I.l)/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)+(1-p)*c*(1-exp((-1/K)*V))*S.l-lambda*E.l
dI.l<-lambda*E.l-b*I.l
dQ<-b*I.l+g*b*I.h-aQ*h*Q-gammaQ*(1-h)*Q
dR<-gammah*(1-h)*(1-g)*I.h+gammaQ*(1-h)*Q
dV<-omega*I.h+(1-p)*omega*I.l-delta*V
  dy<-c(dS.h,dE.h,dI.h,dS.l,dE.l,dI.l,dQ,dR,dV) #combine results into one vector dy
list(dy)
  }
)
}
times <- seq(0,180,by=1) 
params <- c(q=0.2,beta=13,p=.2,c=0,lambda=1/4.43,b=1/0.77,g=0,gammaQ=0.1,gammah=1/2.7,aQ=1/1.93,aH=1/2.7,h=0.00082,omega=0,delta=1,K=10000)
ystart <- c(S.h=(.2)*100000,E.h=0,I.h=1,S.l=(1-0.2)*100000,E.l=0,I.l=0,Q=0,R=0,V=0)
out <- as.data.frame(lsoda(ystart,times,legrand,params))
plot.legrand<-function(out)
  {
  par(mfrow=c(3,2))
  plot(out$time, out$S.h, xlab='Time', ylab='Susceptible(High Risk)',type='l',col='steelblue4', lwd=3)
  plot(out$time, out$E.h, xlab='Time', ylab='Exposed(High Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$I, xlab='Time', ylab='Infectious(High Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$H, xlab='Time', ylab='Susceptible(Low Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$F, xlab='Time', ylab='Exposed(Low Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$I, xlab='Time', ylab='Infectious(Low Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$R, xlab='Time', ylab='Quarantined', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$R, xlab='Time', ylab='Recovered', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$R, xlab='Time', ylab='V', type='l', col='steelblue4', lwd=3)
}
plot.legrand(out)

get.size <- function(out) tail(out$R,1)
output0 <- rep(NA, 180)
times <- seq(1,180)
for(T in times){
  params$T <- T
  out2 <- as.data.frame(ode(ystart, times, legrand, params))
  output0[T] <- get.size(out2)  
}
## Warning in params$T <- T: Coercing LHS to a list
plot(times, output0, log='y', xlab='Intervention time', ylab='Epidemic size', las=1)

h1<-1000
set.seed(6242015)
lhs<-maximinLHS(h1,16)
#using g,q,p,delta
g.min<-0
g.max<-1
q.min<-0
q.max<-1
p.min<-0
p.max<-1
delta.min<-1/7
delta.max<-24
b.max <- 1/0.1
b.min <- 1/5
K.max <- 1000000
K.min <- 1
beta.max <- 20
beta.min <- 1
c.max <- 20
c.min <- 1
S0.max <- 100000
S0.min <- 1
omega.max <- 1
omega.min <- 0
aH.max <- 1/0.1
aH.min <- 1/10
aQ.max <- 1/0.1
aQ.min <- 1/10
h.max <- 1
h.min <- 0
gammah.max <- 1/0.1
gammah.min <- 1/10
gammaQ.max <- 1/0.1
gammaQ.min <- 1/15
lambda.max <- 1/0.1
lambda.min <- 1/10
params.set<-cbind(
g=lhs[,1]*(g.max-g.min)+g.min,
q=lhs[,2]*(q.max-q.min)+q.min,
p=lhs[,3]*(p.max-p.min)+p.min,
delta=lhs[,4]*(delta.max-delta.min)+delta.min,
b = lhs[,5]*(b.max-b.min)+b.min,
K = lhs[,6]*(K.max-K.min)+K.min,
beta = lhs[,7]*(beta.max-beta.min)+beta.min,
c = lhs[,8]*(c.max-c.min)+c.min,
S0 = lhs[,9]*(S0.max-S0.min)+S0.min,
omega = lhs[,10]*(omega.max-omega.min)+omega.min,
aH = lhs[,11]*(aH.max-aH.min)+aH.min,
aQ = lhs[,12]*(aQ.max-aQ.min)+aQ.min,
h = lhs[,13]*(h.max-h.min)+h.min,
gammah = lhs[,14]*(gammah.max-gammah.min)+gammah.min,
gammaQ = lhs[,15]*(gammaQ.max-gammaQ.min)+gammaQ.min,
lambda = lhs[,16]*(lambda.max-lambda.min)+lambda.min
)
levels<-15
h2<-25
j <- 1
 data <- data.frame(matrix(rep(NA,levels*h2*18),nrow=levels*h2))
 for(i in 1:h2){
 for(T in times){

data[j,1:17] <- params <- as.list(c(params.set[i,], T=T))
out <- (as.data.frame(ode(ystart,times, legrand, params)))
data[j,18] <- (get.size(out))
j <- j+1
}
}
names(data) <- c(names(params),'outbreak.size')
save(data, file='data.Rdata')
load('data.Rdata')
plot(times, output0, type='l', lwd=3, ylim=c(10,2e6), log='y',
xlab='Intervention times',
ylab='Epidemic size')
points(data$T, data$outbreak.size, pch=19, cex=0.3, col='blue')

boxplot(data$outbreak.size~data$T, ylim=c(10,2e6), border='blue', log='y', pch='.')
par(new=TRUE)
plot(times, output0, type='l', lwd=3, ylim=c(10,2e6), log='y', xlab='Intervention times', ylab='Epidemic size', axes=FALSE)

library(sensitivity)
## Warning: package 'sensitivity' was built under R version 4.0.2
bonferroni.alpha <- 0.05/16
prcc <- pcc(data[,1:16], data[,18], nboot = 1000, rank=TRUE, conf=1-bonferroni.alpha)
save(prcc, file='prcc.Rdata')
library(sensitivity)
load('prcc.Rdata')
summary <- print(prcc)
## 
## Call:
## pcc(X = data[, 1:16], y = data[, 18], rank = TRUE, nboot = 1000,     conf = 1 - bonferroni.alpha)
## 
## Partial Rank Correlation Coefficients (PRCC):
##            original          bias  std. error   min. c.i.   max. c.i.
## g      -0.096849049  1.461844e-04 0.015951652 -0.14235246 -0.05279976
## q       0.609990848  5.829482e-05 0.006114992  0.59280322  0.62903725
## p      -0.437167920 -1.670473e-04 0.012335664 -0.47368763 -0.40539258
## delta   0.531313700 -6.659116e-05 0.009913353  0.50215313  0.56198914
## b      -0.217691867  4.402021e-04 0.012923684 -0.26018748 -0.17784593
## K       0.498031882  1.882083e-04 0.010119482  0.47210707  0.53045533
## beta    0.626631802  1.118212e-04 0.008501931  0.60375692  0.65668660
## c      -0.378626882 -4.423300e-05 0.010244162 -0.41168993 -0.35093238
## S0     -0.358721842 -6.513590e-04 0.012432332 -0.39791391 -0.32435660
## omega  -0.023413221  3.266295e-05 0.012667098 -0.05944101  0.01283078
## aH      0.107691837  9.027792e-06 0.013550347  0.07005795  0.15053890
## aQ     -0.306344440  4.892748e-04 0.011649492 -0.34653075 -0.27069634
## h      -0.910581358  1.360331e-05 0.001962806 -0.91723244 -0.90452916
## gammah  0.229512479 -3.996812e-05 0.012037543  0.18947740  0.26247263
## gammaQ  0.470185048 -5.020093e-04 0.011790250  0.43524259  0.50365929
## lambda -0.003476009  3.975086e-04 0.014986865 -0.05080559  0.04491250
par(mar=c(7,4,4,2)+0.1)
plot(summary$original, main='Partial rank correlation coefficients', ylim=c(-1,1),
     xlab='', ylab='Coefficient',
     axes=FALSE)
axis(2)
axis(1, at=seq(1:16), labels=row.names(summary), las=2)
mtext(text='Parameter', side=1, line=4.5)
box()
for(i in 1:16) lines(c(i,i),c(summary[i,4], summary[i,5]))
abline(h=0)

.