library(crone)
library(deSolve)

average_death <- read.table("death_rate_avg_0A_2.csv", header=TRUE, sep=",")

t_max=24
initial_param = c(3,0.25,0.02)


###########################
#  least squares function #
###########################

least.squares.fn <- function(par,data){
  
  delta_xphi <- par[1]
  A <- par[2]
  B <- par[3]
  
  run <- run_phys_disruption_eqn(delta_xphi,A,B)
  
  run_data_points <- run$ex[c(1,2,4,6,8,25)]
  
  
  if (delta_xphi < 0){
    error <- 1e8 
  }
  else{
    error <- 0
    for (i in 1:length(data)){
      timepoint <- i
      add_to_error <- ((data[i]-run_data_points[timepoint])^2)
      error <- error + add_to_error
    }
    error <- error/length(data)
  }
  
  return(error)
}




##############
#  ODE model #
##############


run_phys_disruption_eqn <- function(delta_xphi,A,B){
  #times <- seq(from=1,to=end.time,by=.1) 
  time <- 0:24
  ex <- rep(0,25)
  for (i in time){
    t <- i
    ex[i+1] <- A*exp(-delta_xphi*t)+B
  }
  output <- as.data.frame(cbind(time, ex))
  return(output)
}


###########################
#  OPTIMISING PARAMETERS  #
###########################

#fitted value: delta_xphi=1.11

model_fitting_death_rate <- function(start){
  fn <- function(par) least.squares.fn(par,data=average_death$live.bacteria.fraction)
  fit <- optim(start,fn)
  while(identical(signif(fit$par,digits=5),signif(start,digits=5))==FALSE) {
    start <- fit$par
    fit <- optim(start,fn)
  }
  delta_xphi_fit <- fit$par[1]
  A_fit <- fit$par[2]
  B_fit <- fit$par[3]
  runfit <- run_phys_disruption_eqn(delta_xphi_fit,A_fit,B_fit)
  plot(average_death$time, average_death$live.bacteria.fraction, col = "red", pch=18, ylim = c(-0.01,.3),
       xlab = "time (hours)",
       ylab = "ex")
  arrows(x0=average_death$time, y0=average_death$live.bacteria.fraction-average_death$sd, x1=average_death$time, y1=average_death$live.bacteria.fraction+average_death$sd, code=3, angle=90, length=0.05)
  lines(runfit$time, runfit$ex, col="blue")
  legend("topright", legend=c("data","model"),
         col=c("red", "blue"), pch= c(18,NA), lty=c(NA,1))
  outputs <- list(delta_xphi_fit, runfit)
  return(fit$par)
}
model_fitting_death_rate(initial_param)






###############################################
#  PARAMETER ESTIMATION USING THE ABC METHOD  #
###############################################
N <- 500
epsilon <-0.000085
accepted_values <- matrix(0,500,1)
i <- 0
while(i<N){
  prior.val <- c(runif(1,0,100),0.24,0.03)
  least.squares.run <- least.squares.fn(prior.val, data=average_death$live.bacteria.fraction)
  if(least.squares.run<epsilon){
    i <- i + 1
    accepted_values[i] <- prior.val[1]
  }
  cat("\r",round(i/N*100),'% done')
}
par(mfrow=c(1,1))
param_means <- c()
parameter <- c(expression(delta[x*phi]))
 for (i in 1:1) {
  hist(accepted_values[,i], main = expression(paste("Histogram of accepted values of ", delta[x*phi])), xlab = parameter[i])
  param_means <- c(param_means, mean(accepted_values[,i]))
}
print(param_means)