---
title: "Basic Predator-Prey"
author: "Bill Foote"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(cmdstanr)
library(dplyr)
library(ggplot2)
library(posterior)
library(readr)
library(readsdr) # Development version
library(tidyr)
```


```{r}
model_path <- "model/stmx_files/SEIR.stmx"
mdl        <- read_xmile(model_path)
mdl$deSolve_components
#
constants_df <- sd_constants(mdl) |> 
  mutate(value = format(value, scientific = FALSE))
constants_df
#
meas_mdl <- list("y ~ neg_binomial_2(net_flow(C), 10)")

set.seed(789)

dat_df <- sd_measurements(n_meas       = 1,
                          ds_inputs    = mdl$deSolve_components,
                          meas_model   = meas_mdl,
                          start_time   = 0,
                          stop_time    = 70,
                          timestep     = 1/16,
                          integ_method = "rk4")
summary(dat_df)
p <- dat_df |>
ggplot(aes(time, measurement)) +
  geom_point(colour = "#7B92DC") +
  labs(x = "Day", y = "Incidence [People per day]") +
  theme_classic()
p
#
meas_mdl <- list("y ~ neg_binomial_2(net_flow(C), phi)")

est_params <- list(sd_prior(par_name  = "par_beta", 
                            dist      = "lognormal",
                            dist_pars =  c(0, 1)),
                   sd_prior("par_rho", "beta", c(2, 2)),
                   # This parameter only affects initial conditions 
                   sd_prior("I0", "lognormal", c(0, 1), type = "init"))

stan_file <- sd_Bayes(filepath         = model_path,
                      meas_mdl         = meas_mdl,
                      estimated_params = est_params)

cat(stan_file)
```

```{r}
#
cat(Lynx_Hare_model)
#
stan_path <- "./Stan_files/SEIR.stan"
#write_file(stan_file, file = stan_path)
stan_d <- list(n_obs = nrow(dat_df),
               y     = dat_df$measurement,
               ts    = 1:max(dat_df$time))

mod           <- cmdstan_model(stan_path)

fit <- mod$sample(data              = stan_d,
                  chains          = 4,
                  parallel_chains = 4,
                  iter_warmup     = 1000,
                  iter_sampling   = 1000,
                  refresh         = 100,
                  save_warmup     = FALSE)
SEIR_mod <- mod
path2Stan <-  "Stan_files/SEIR_mod"
save(SEIR_mod, file="Stan_files/SEIR_mod")
#
# fit$cmdstan_summary()
#
fit$cmdstan_diagnose()
#
posterior_df <- as_draws_df(fit$draws())
str(posterior_df)
# Useful expanded version of base summary() 
# Actually summarizes cmdstan_summary()
summarize_predicted_observations <- function(posterior_df, var_name) {
  
  y_tilde <- readsdr::extract_timeseries_var(var_name, posterior_df)
  
  y_tilde |> group_by(time) |>
  summarize(min   = min( value ),
            q5.5  = quantile(value, 0.055),
            q25   = quantile(value, 0.25),
            q50   = quantile(value, 0.50),
            q75   = quantile(value, 0.75),
            q94.5 = quantile(value, 0.945),
            max   = max( value ),
            mean  = mean( value ),
            stdev = sd( value ),
            mad( value )
            )
}
#
summary_y <- summarize_predicted_observations(posterior_df, "sim_y")
p <- summary_y |>
ggplot(aes(time, mean)) +
  geom_line(colour = "grey60") +
  geom_ribbon(aes(ymin = q25, ymax = q75), fill = "grey90", alpha = 0.5) +
  geom_ribbon(aes(ymin = q5.5, ymax = q94.5), fill = "grey90", alpha = 0.25) +
  geom_point(data = dat_df, aes(time, measurement), colour = "#7B92DC") +
  labs(y = "Incidence", x = "Day") +
  theme_classic() +
  theme(axis.line = element_line(colour = "grey70"),
        axis.text = element_text(colour = "grey70"),
        axis.ticks = element_line(colour = "grey70"))
p
#
# Parameter estimation
#
pars_df <- posterior_df[, c("par_beta", "I0", "par_rho", "inv_phi")] |> 
  mutate(par_phi = 1 / inv_phi) |> 
  dplyr::select(-inv_phi) |> 
  mutate(iter = row_number()) |> 
  pivot_longer(-iter, names_to = "par")
3.1.3 Parameter estimates (marginal distributions)
Vertical lines denote the true value.

actual_vals <- data.frame(par   = c("I0", "par_beta", "par_rho", "par_phi"),
                          value = c(1, 1, 0.75, 10))
ggplot(pars_df, aes(value)) +
  geom_histogram(colour = "white", fill = "grey60", alpha = 0.75) +
  facet_wrap(vars(par), scales = "free") +
  geom_vline(data = actual_vals, aes(xintercept = value), 
             colour = "#7B92DC", linewidth = 1.5) +
  theme_classic()
```



