model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    d <- f*r/(1-f)
    N <- S + I + R
    dS <- w*R - beta*S*I/N
    dI <- beta*S*I/N - (d+r)*I
    dR <- r*I - w*R
    return(list(c(dS,dI,dR)))  
  }) 
} 

Tweak <- function(state, parms, nsol) {
  with(as.list(c(state,parms)), {
    d <- f*r/(1-f)
    nsol$D <- d*nsol$I   # D defines excess deaths
    return(nsol)
  })
}
tweak <- "nsol<-Tweak(state,parms,nsol)"

p <- c(beta=0.25,r=0.1,f=0.002,w=0.006)  
# Recovery takes 10d, 0.2% dies, R0=beta/r=2.5 -> beta=rR0=0.25. 
s <- c(S=2.2e6,I=1e4,R=0)
# Daily excess deaths from 1 April 2020 to 1 April 2021
data <- read.csv("data/excessDeaths.csv",header=TRUE)  
timePlot(data,draw=points,main="Excess deaths 1 April 2020 to 1 April 2021")
run(200,1,tweak=tweak,ymax=150,show="D")
free <- c("beta","r","f","w","I")
# Fit the first peak (the first 200 days)
fit1 <- fit(data[1:200,],free=free,tweak=tweak,lower=0,upper=c(1,0.5,0.1,0.1,1e5),ymax=150,show="D",method="Pseudo")
summary(fit1)
p[free[1:4]] <- fit1$par[1:4]
s["I"] <- fit1$par["I"]

d <- with(as.list(p),f*r/(1-f))
R0 <- with(as.list(p),beta/(r+d))
rho0 <- with(as.list(p),(R0-1)*(d+r))
print(c(d=d,R0=R0,p)); print(s); 
print(c(rho0=rho0,tStart=log(fit1$par["I"])/rho0))

# Find confidence ranges of the parameters by bootstrapping the data (100x)
# fit2 <- fit(data[1:200,],free=free,tweak=tweak,lower=0,ymax=150,show="D",bootstrap=100)

# Try to fit both peaks (all data)
fit3 <- fit(data,free=free,tweak=tweak,lower=0,upper=c(1,0.5,0.1,0.1,1e5),ymax=150,show="D",method="Pseudo")




# Start with one infected individual and shift time 45 days.
# s["I"] <- 1
# run(250,1,ymax=150,tweak=tweak,show="D")
# points(data$time+45,data$D)

# In case you prefer to download the data yourself:
# library("readxl")
# alldata <- read_excel("excess_deaths_manaus_01apr2020_01may2021_compact.xlsx")
# daily <- alldata[!is.na(alldata$Covid_period_excess_deaths),]
# ndays <- nrow(daily)
# data <- data.frame(cbind(seq(0,ndays-1),daily$Covid_period_excess_deaths))
# colnames(data) <- c("time","D")
# timePlot(data,draw="points")