---
title: "deSolve to denim"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{deSolve to denim}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r}
library(denim)
library(deSolve)
```

## Migrate deSolve code to denim

### Original code in deSolve

The model used for demonstrating the process of migrating code from `deSolve` to `denim` is as followed

```{r}
# --- Model definition in deSolve
transition_func <- function(t, state, param){
  with(as.list( c(state, param) ), {
      gamma_rate = 1/scale
      
      dS = -beta*S*I/N
      dI1 = beta*S*I/N - gamma_rate*I1
      dI2 = gamma_rate*I1 - gamma_rate*I2
      dI =  dI1 + dI2
      dR = gamma_rate*I2
      list(c(dS, dI, dI1, dI2, dR))
  })
}

# ---- Model configuration 
parameters <- c(beta = 0.3, scale = 3, N = 1000) 
initialValues <- c(S = 999, I = 1, I1 = 1, I2=0, R=0)

# ---- Run simulation
times <- seq(0, 100) # simulation duration
ode_mod <- ode(y = initialValues, times = times, parms = parameters, func = transition_func)

# --- show output
ode_mod <- as.data.frame(ode_mod)
head(ode_mod[-1, c("time", "S", "I", "R")])
```

### Model definition

Unlike `deSolve` where transitions between compartments are defined by a system of ODEs, transitions in `denim` must be defined by (i) a distribution of dwell-time, (ii) a mathematical expression, (iii) a fixed number or proportion.

User must first identify the transitions that best describe the ones in their `deSolve` model.

```{=html}
<details>
  <summary>Identify distribution</summary>
```
```{r, eval=FALSE}
# --- Model definition in deSolve
transition_func <- function(t, state, param){
  with(as.list( c(state, param) ), {
      gamma_rate = 1/scale
      
      # For S -> I transition, since it involves parameters (beta, N), 
      # the best transition to describe this is using a mathematical formula
      dS = -beta*S*I/N
      
      # For I -> R transition, linear chain trick is applied --> implies Erlang distributed dwell time
      # Hence, we can use d_gamma from denim
      dI1 = beta*S*I/N - gamma_rate*I1
      dI2 = gamma_rate*I1 - gamma_rate*I2
      dI =  dI1 + dI2
      dR = gamma_rate*I2
      list(c(dS, dI, dI1, dI2, dR))
  })
}
```

</details>

With the transitions identified, user can then define the model in `denim`.

When using `denim`, the model structure is given as a `list` of *key*-*value* pairs where

-   *key* is a string showing the transition direction between compartments

-   *value* is the built-in distribution function that describe the transition

```{r}
# --- Transition def for denim
transitions <- list(
  "S -> I" = "beta * S * I/N * timestepDur",
  "I -> R" = d_gamma(3, 2) # shape is 2 from number of I sub compartments
)
```

Note that when converting an ODE to a math expression, we must multiply the original expression with time step duration (variable `timestepDur` in the model definition code). This is because we are essentially using Euler's method to estimate the solution for the ODE.

### Model configurations

Similar to `deSolve`, `denim` also ask users to provide the [initial values]{.underline} and any [additional parameters]{.underline} in the form of named vectors.

For the example `deSolve` code, while users can use the `initalValues` from the `deSolve` code as is (`denim` will ignore unused `I1`, `I2` compartments as these sub-compartments will be automatically computed internally), it is recommended to remove redundant compartments (in this example, `I1` and `I2`).

For parameters, since `rate` and `scale` are already defined in the distribution functions, users only need to keep `beta` and `N` from the initial parameters vector. We also need to specify the value for `timestepDur` here.

```{r}
# remove I1, I2 compartments
denim_initialValues <- c(S = 999, I = 1, R=0)
denim_parameters <- c(beta = 0.3, N = 1000, timestepDur=0.01) 
```

***Initialization of sub-compartments*****:** when there are multiple sub-compartments (e.g., compartment `I` consist of `I1` and `I2` sub-compartments), the initial population is always assigned to the first sub-compartment. In our example, since `I = 1`, denim will assign `I1 = 1` and `I2 = 0`.

### Simulation

Lastly, users need to define the simulation duration and time step for `denim` to run. Unlike `deSolve` which takes a time sequence, `denim` only require the [simulation duration]{.underline} and [time step.]{.underline}

Since `denim` is a discrete time model, time step must be set to a small value for the result to closely follow that of `deSolve` (in this example, 0.01).

```{r}
mod <- sim(transitions = transitions,
             initialValues = denim_initialValues, 
             parameters = denim_parameters,
             simulationDuration = 100,
             timeStep = 0.01)
```

## Compare output

The following plots show output from `denim` and `deSolve`

```{r, echo=FALSE, fig.width=8, fig.height=5}
# ---- Plot S compartment
plot(x = mod$Time, y = mod$S,xlab = "Time", ylab = "Count", main="S compartment",
     col = "#4876ff", type="l", lwd=3)
lines(ode_mod$time, ode_mod$S, lwd=3, lty=3)
legend(x = 15, y = 4e5,legend=c("denim", "deSolve"), col = c("#4876ff", "black"), lty=c(1,3))

# ---- Plot I compartment
plot(x = mod$Time, y = mod$I, xlab = "Time", ylab = "Count", main="I compartment",
      col = "#4876ff", type="l", lwd=2)
lines(ode_mod$time, ode_mod$I, lwd=3, lty=3)
legend(x = 15, y = 1e5,legend=c("denim", "deSolve"), col = c("#4876ff", "black"), lty=c(1,3))

# ---- Plot R compartment
plot(x = mod$Time, y = mod$R, xlab = "Time", ylab = "Count", main="R compartment",
     col = "#4876ff", type="l", lwd=2)
lines(ode_mod$time, ode_mod$R, lwd=3, lty=3)
legend(x = 15, y = 4e5,legend=c("denim", "deSolve"), col = c("#4876ff", "black"), lty=c(1,3))
```