## -----------------------------------------------------------------------------
library(bioinactivation)

## -----------------------------------------------------------------------------
data(isothermal_inactivation)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
library(ggplot2)
library(purrr)
library(tidyr)
library(dplyr)
ggplot(data = isothermal_inactivation) +
    geom_point(aes(x = time, y = log_diff, col = as.factor(temp))) +
    ggtitle("Example dataset: isothermal_inactivation")


## -----------------------------------------------------------------------------
data(dynamic_inactivation)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
ggplot(data = dynamic_inactivation) +
    geom_point(aes(x = time, y = log10(N))) +
    ggtitle("Example dataset: dynamic_inactivation. Observations")

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
ggplot(data = dynamic_inactivation) +
    geom_point(aes(x = time, y = temperature))  +
    ggtitle("Example dataset: dynamic_inactivation. Temperature profile")

## -----------------------------------------------------------------------------
data(laterosporus_iso)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
ggplot(data = laterosporus_iso) +
    geom_point(aes(x = time, y = log_diff, col = as.factor(temp))) +
    ggtitle("Example dataset: laterosporus_iso")


## -----------------------------------------------------------------------------
data(laterosporus_dyna)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
ggplot(data = laterosporus_dyna) +
    geom_point(aes(x = time, y = logN)) +
    ggtitle("Example dataset: laterosporus_dyna. Observations")

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
ggplot(data = laterosporus_dyna) +
    geom_point(aes(x = time, y = temp))  +
    ggtitle("Example dataset: laterosporus_dyna. Temperature profile")

## ---- echo=FALSE, fig.width=6, fig.height=4, fig.align='center'---------------

tibble(
  Temperature = seq(30, 80, length = 100),
  `k = 0.3` = log( 1 + exp(0.3*(Temperature - 60))),
  `k = 0.2` = log( 1 + exp(0.2*(Temperature - 60)))
  ) %>%
  pivot_longer(-Temperature, names_to = "k", values_to = "b") %>%
  ggplot() +
  geom_line(aes(x = Temperature, y = b, colour = k))

## -----------------------------------------------------------------------------
get_model_data()

## -----------------------------------------------------------------------------
example_model <- "Geeraerd"

## -----------------------------------------------------------------------------
times <- seq(0, 4.5, length=100)

## -----------------------------------------------------------------------------
model_data <- get_model_data(example_model)
print(model_data$parameters)
print(model_data$variables)

## -----------------------------------------------------------------------------
model_parms <- c(D_R = 3,
                 z = 10,
                 N_min = 1e2,
                 temp_ref = 100,
                 N0 = 1e5,
                 C_c0 = 1e1
                 )

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
temperature_profile <- data.frame(time = c(0, 4.5),
                                  temperature = c(70, 120))

ggplot(temperature_profile) +
  geom_line(aes(time, temperature))

## -----------------------------------------------------------------------------
prediction_results <- predict_inactivation(example_model, times, model_parms, temperature_profile)

## -----------------------------------------------------------------------------
head(prediction_results$simulation)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
p <- plot(prediction_results)
print(p)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
library(ggplot2)
p + theme_light() + xlab("time (min)")

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
plot(prediction_results, make_gg = FALSE,
     xlab = "Time (min)", ylab = "logN")

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
parms_no_shoulder <- c(D_R = 3,
                       z = 10,
                       N_min = 100,
                       temp_ref = 100,
                       N0 = 100000,
                       C_c0 = 0
                       )


prediction_no_shoulder <- predict_inactivation(example_model, times, parms_no_shoulder, temperature_profile)

plot(prediction_no_shoulder)


## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
parms_no_tail <- c(D_R = 3,
                   z = 10,
                   N_min = 0,
                   temp_ref = 100,
                   N0 = 100000,
                   C_c0 = 100
                   )


prediction_no_tail <- predict_inactivation(example_model, times, parms_no_tail, temperature_profile)

plot(prediction_no_tail)


## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
plot(prediction_results, plot_temp = TRUE)

## -----------------------------------------------------------------------------
data(isothermal_inactivation)

## -----------------------------------------------------------------------------
get_isothermal_model_data()

## -----------------------------------------------------------------------------
model_name <- "Bigelow"

## -----------------------------------------------------------------------------
model_data <- get_isothermal_model_data(model_name)
model_data$params

## -----------------------------------------------------------------------------
known_params = list(temp_ref = 100)

## -----------------------------------------------------------------------------
starting_point <- list(z = 10,
                       D_R = 1.5
                       )

## -----------------------------------------------------------------------------
iso_fit <- fit_isothermal_inactivation(model_name, isothermal_inactivation,  
                                       starting_point, known_params)

## -----------------------------------------------------------------------------
summary(iso_fit$nls)
vcov(iso_fit$nls)  # Calculates variance-covariance matrix {stats}
confint(iso_fit$nls)  # Calculates confidence intervals {stats}

## -----------------------------------------------------------------------------
iso_fit$parameters

## -----------------------------------------------------------------------------
iso_fit$model

## -----------------------------------------------------------------------------
head(iso_fit$data)

## -----------------------------------------------------------------------------
plot(iso_fit, ylab = "Number of logarithmic reductions",
     xlab = "Time (min)")

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
plot(iso_fit, make_gg = TRUE) +
    ylab("Number of logarithmic reductions") +
    xlab("Time (min)")

## -----------------------------------------------------------------------------
data(dynamic_inactivation)

## -----------------------------------------------------------------------------
get_model_data()

## -----------------------------------------------------------------------------
simulation_model <- "Mafart"

## -----------------------------------------------------------------------------
dummy_temp <- data.frame(time = dynamic_inactivation$time,
                         temperature = dynamic_inactivation$temperature) %>%
  distinct()

## -----------------------------------------------------------------------------
model_data <- get_model_data(simulation_model)
model_data$parameters
model_data$variables

## -----------------------------------------------------------------------------
known_params = c(temp_ref = 90)

## -----------------------------------------------------------------------------
starting_points <- c(delta_ref = 10,
                     p = 1,
                     z = 10,
                     logN0 = 5)
upper_bounds <- c(delta_ref = 20,
                  p = 2,
                  z = 20,
                  logN0 = Inf)

lower_bounds <- c(delta_ref = 5,
                  p = .5,
                  z = 5,
                  logN0 = 4)

## -----------------------------------------------------------------------------
dynamic_fit <- fit_dynamic_inactivation(dynamic_inactivation, simulation_model, dummy_temp,  
                                        starting_points, upper_bounds, lower_bounds,  
                                        known_params)

## -----------------------------------------------------------------------------
summary(dynamic_fit$fit_results)

## -----------------------------------------------------------------------------
head(dynamic_fit$data)

## -----------------------------------------------------------------------------
dynamic_fit$best_prediction$model_parameters

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
p <- plot(dynamic_fit)
p + theme_light()

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
plot(dynamic_fit, plot_temp = TRUE)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
plot(dynamic_fit, make_gg = FALSE,
     xlab = "Time (min)", ylab = "logN")

## -----------------------------------------------------------------------------
set.seed(82619)

MCMC_fit <- fit_inactivation_MCMC(dynamic_inactivation, simulation_model, dummy_temp,  
                                        starting_points, upper_bounds, lower_bounds,  
                                        known_params, niter = 50)

## -----------------------------------------------------------------------------
summary(MCMC_fit$modMCMC)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
p <- plot(MCMC_fit)
p + xlab("Time (min)")

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
plot(MCMC_fit, plot_temp = TRUE)

## ---- fig.width=6, fig.height=4, fig.align='center'---------------------------
plot(MCMC_fit, make_gg = FALSE,
     xlab = "Time (min)", ylab = "logN")