## -------------------------------------------------------------------
library(SimInf)

## -------------------------------------------------------------------
transitions <- "S -> beta * S * I / (S + I) -> I"

## -------------------------------------------------------------------
transitions <- c(
  "S -> beta * S * I / (S + I + R) -> I",
  "I -> gamma * I -> R"
)

## -------------------------------------------------------------------
model <- mparse(
  transitions = transitions,
  compartments = c("S", "I", "R"),
  gdata = c(beta = 0.16, gamma = 0.077),
  u0 = data.frame(S = 99, I = 1, R = 0),
  tspan = 1:100
)

## ----fig.width=7, fig.height=4, fig.align="left", fig.cap="**Figure 1.** Classic SIR epidemic curve generated with mparse."----
plot(run(model, seed = 22))

## -------------------------------------------------------------------
transitions <- c(
  "S -> beta * S * I / N -> I",
  "I -> gamma * I -> R",
  "N <- S + I + R"
)

## -------------------------------------------------------------------
transitions <- c(
  "S -> beta * S * I / N -> I",
  "I -> gamma * I -> R",
  "(int)N <- S + I + R"
)

## -------------------------------------------------------------------
model <- mparse(
  transitions = transitions,
  compartments = c("S", "I", "R"),
  gdata = c(beta = 0.16, gamma = 0.077),
  u0 = data.frame(S = 99, I = 1, R = 0),
  tspan = 1:100
)

## ----fig.width=7, fig.height=4, fig.align="left", fig.cap="**Figure 2.** SIR epidemic curve using a defined variable for population size. The results are identical to Figure 1."----
plot(run(model, seed = 22))

## -------------------------------------------------------------------
transitions <- c(
  "S -> N > 0 ? beta * S * I / N : 0 -> I",
  "I -> gamma * I -> R",
  "(int)N <- S + I + R"
)

## ----fig.width=7, fig.height=4, fig.align="left", fig.cap="**Figure 3.** SIR model using a ternary operator to prevent division by zero. Although the curve is identical to previous examples (as the population did not reach zero), this syntax ensures the simulation continues safely if the node becomes empty."----
model <- mparse(
  transitions = transitions,
  compartments = c("S", "I", "R"),
  gdata = c(beta = 0.16, gamma = 0.077),
  u0 = data.frame(S = 99, I = 1, R = 0),
  tspan = 1:100
)

plot(run(model, seed = 22))

## -------------------------------------------------------------------
transitions <- c(
  "S -> N > 0 ? beta_farm * S * I / N : 0 -> I",
  "I -> gamma * I -> R",
  "(int)N <- S + I + R"
)

## -------------------------------------------------------------------
gdata <- c(gamma = 0.077)

ldata <- data.frame(
  beta_farm = c(0.1, 0.4)  # Farm 1: 0.1, Farm 2: 0.4
)

## -------------------------------------------------------------------
u0 <- data.frame(
  S = c(99, 95),  # Farm 1: 99 S, Farm 2: 95 S
  I = c(1, 5),    # Farm 1: 1 I, Farm 2: 5 I
  R = c(0, 0)     # Both start with 0 R
)

## -------------------------------------------------------------------
model <- mparse(
  transitions = transitions,
  compartments = c("S", "I", "R"),
  gdata = gdata,
  ldata = ldata,
  u0 = u0,
  tspan = 1:100
)

## ----fig.width=7, fig.height=4, fig.align="left", fig.cap="**Figure 4.** Epidemic curves for two farms with different transmission rates in a single stochastic realization. Farm 2 (higher `beta_farm`) shows a faster outbreak compared to Farm 1, reflecting the expected impact of the higher transmission rate, though exact outcome varies due to randomness."----
result <- run(model, seed = 22)
plot(result, range = FALSE)

## ----fig.width=7, fig.height=4, fig.align="left", fig.cap="**Figure 5.** Trajectory for Farm 2 only, showing the rapid spread due to the high transmission rate in this specific realization."----
plot(result, index = 2)

## -------------------------------------------------------------------
transitions <- c(
  "@ -> bR * R -> R",
  "R -> (dR + (bR - dR) * R / K) * R -> @",
  "R -> alpha / (1 + w * R) * R * F -> @",
  "@ -> bF * alpha / (1 + w * R) * R * F -> F",
  "F -> dF * F -> @"
)

## -------------------------------------------------------------------
parameters <- c(
  bR = 2, bF = 2, dR = 1, K = 1000,
  alpha = 0.007, w = 0.0035, dF = 2
)

u0 <- data.frame(R = 1000, F = 100)

## -------------------------------------------------------------------
model <- mparse(
  transitions = transitions,
  compartments = c("R", "F"),
  gdata = parameters,
  u0 = u0,
  tspan = 1:100
)

## ----fig.width=7, fig.height=6, fig.align="left", fig.cap="**Figure 6.** One stochastic realization of the Rosenzweig-MacArthur predator-prey model. The populations exhibit cyclic oscillations characteristic of the deterministic limit cycle, but stochastic fluctuations eventually lead to the extinction of the predator population in this specific run."----
result <- run(model, seed = 60)
plot(result)

## ----fig.width=7, fig.height=6, fig.align="left", fig.cap="**Figure 7.** Phase plane trajectory of the predator-prey model. The path spirals outward from the initial conditions, tracing the characteristic limit cycle of the Rosenzweig-MacArthur model before stochastic fluctuations drive the predator population to extinction."----
plot(F ~ R, data = trajectory(result), type = "l", col = "darkgreen",
     xlab = "Prey (R)", ylab = "Predators (F)", main = "")

## -------------------------------------------------------------------
transitions <- c(
  "S -> N > 0 ? beta * S * I / N : 0 -> I1",
  "I1 -> gamma * I1 -> I2",
  "I2 -> gamma * I2 -> I3",
  "I3 -> gamma * I3 -> I4",
  "I4 -> gamma * I4 -> R",
  "I <- I1 + I2 + I3 + I4",
  "N <- S + I + R"
)

## -------------------------------------------------------------------
model <- mparse(
  transitions = transitions,
  compartments = c("S", "I1", "I2", "I3", "I4", "R"),
  gdata = c(beta = 0.16, gamma = 0.077),
  u0 = data.frame(S = 100, I1 = 5, I2 = 0, I3 = 0, I4 = 0, R = 0),
  tspan = 1:100
)

## ----fig.width=7, fig.height=4, fig.align="left", fig.cap="**Figure 8.** One realization of a stochastic SIR model with an Erlang-distributed infectious period (shape=4). The trajectory shows the infection progressing sequentially through the four infectious stages (`I1 -> I2 -> I3 -> I4`) before individuals recover."----
result <- run(model, seed = 22)
plot(result)

## ----fig.width=7, fig.height=4, fig.align="left", fig.cap="**Figure 9.** Prevalence of infection over time, calculated as the sum of individuals in all infectious stages (`I1 + I2 + I3 + I4`). This demonstrates how to aggregate multiple compartments to track the total burden of disease in a stage-structured model."----
plot(result, I1 + I2 + I3 + I4 ~ ., col = "blue", lwd = 2, ylab = "Prevalence")