---
title: "Sobol sensitivity analysis of an M/M/1 queue in simmer"
shorttitle: "Sobol sensitivity analysis of an M/M/1 queue in simmer"
author:
- name: "Frédéric Bertrand"
  affiliation:
  - Cedric, Cnam, Paris
  email: frederic.bertrand@lecnam.net
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Sobol sensitivity analysis of an M/M/1 queue in simmer}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "figures/sobol-stochastic-",
  fig.width = 7,
  fig.height = 5,
  dpi = 150,
  message = FALSE,
  warning = FALSE,
  eval=FALSE
)

LOCAL <- identical(Sys.getenv("LOCAL"), "TRUE")

library(sensitivity)
library(Sobol4R)
set.seed(4669)
```

# Model description

We consider a simple M/M/1 queue simulated with the `simmer` package.  
Two input parameters are uncertain

- `lambda` arrival rate (interarrival times distributed as Exp(lambda))  
- `mu` service rate (service times distributed as Exp(mu))  

We study the impact of these two parameters on the mean waiting time in the
queue at stationarity using Sobol indices.

# Simulation model
```{r model-once, cache=TRUE, eval=LOCAL}
# M/M/1 simmer model for Sobol analysis using time in system as QoI

library(simmer)
# One replication: mean time in system (sojourn time)
simulate_mm1_once_ts <- function(lambda, mu,
                                 horizon = 1000,
                                 warmup  = 200) {
  traj <- trajectory("customer") %>%
    seize("server") %>%
    timeout(function() rexp(1, rate = mu)) %>%
    release("server")
  
  env <- simmer("MM1_queue") %>%
    add_resource("server", capacity = 1, queue_size = Inf) %>%
    add_generator(
      name_prefix  = "cust",
      trajectory   = traj,
      distribution = function() rexp(1, rate = lambda)
    )
  
  env <- run(env, until = horizon)
  arr <- get_mon_arrivals(env)
  
  # Keep only post warmup
  arr <- arr[arr$end_time > warmup, , drop = FALSE]
  
  if (nrow(arr) == 0L) {
    return(0)
  }
  
  # Time in system = end_time - start_time
  ts <- arr$end_time - arr$start_time
  
  m <- mean(ts, na.rm = TRUE)
  if (!is.finite(m)) m <- 0
  m
}

# Quick sanity check
simulate_mm1_once_ts(lambda = 1/5, mu = 1/4)
```

To reduce Monte Carlo noise, we define a quantity of interest (QoI) equal to the
mean waiting time over several independent replications for the same parameter
pair `(lambda, mu)`.

```{r model-qoi, cache=TRUE, eval=LOCAL}

# QoI: average mean time in system over replications
simulate_mm1_qoi_ts <- function(lambda, mu,
                                horizon = 1000,
                                warmup  = 200,
                                nrep   = 20L) {
  nrep <- as.integer(nrep)
  if (nrep < 1L) stop("'nrep' must be at least 1")
  vals <- replicate(
    nrep,
    simulate_mm1_once_ts(lambda = lambda, mu = mu,
                         horizon = horizon, warmup = warmup)
  )
  mean(vals)
}

simulate_mm1_qoi_ts(lambda = 1/5, mu = 1/4)
```

# Sobol model wrapper

The Sobol routine expects a model that takes a design matrix `X` and returns a
numeric vector of length `nrow(X)`. Here `X` has two columns

1. `lambda`  
2. `mu`

```{r sobol-model, cache=TRUE, eval=LOCAL}
# Model wrapper for sensitivity::sobol
mm1_sobol4r_model_ts <- function(X,
                               horizon = 1000,
                               warmup  = 200,
                               nrep   = 20L) {
  X <- as.data.frame(X)
  apply(X, 1L, function(par) {
    lambda <- par[1]
    mu     <- par[2]
    val <- simulate_mm1_qoi_ts(lambda = lambda, mu = mu,
                               horizon = horizon,
                               warmup  = warmup,
                               nrep   = nrep)
    if (!is.finite(val)) val <- 0
    val
  })
}
```

# Sobol design for lambda and mu

We assume independent uniform priors

- `lambda` in `[0.10, 0.30]`  
- `mu` in `[0.20, 0.60]`  

```{r design, cache=TRUE, eval=LOCAL}
n <- 200

X1 <- data.frame(
  lambda = runif(n, min = 0.10, max = 0.30),
  mu     = runif(n, min = 0.20, max = 0.60)
)

X2 <- data.frame(
  lambda = runif(n, min = 0.10, max = 0.30),
  mu     = runif(n, min = 0.20, max = 0.60)
)

# Helper to build a design region with non trivial queues but stable system
# choose lambda < mu, with rho close to 1
example_design_mm1_ts <- function(n = 100) {
  lambda <- runif(n, min = 0.18, max = 0.22)
  mu     <- runif(n, min = 0.23, max = 0.27)
  data.frame(lambda = lambda, mu = mu)
}
```

# Sobol indices

```{r sobol-run, message=FALSE, cache=TRUE, eval=LOCAL}
library(simmer)
library(sensitivity)

set.seed(4669)

n  <- 150
X1 <- example_design_mm1_ts(n)
X2 <- example_design_mm1_ts(n)

sob_mm1 <- sobol(
  model = NULL,
  X1    = X1,
  X2    = X2,
  order = 2,
  nboot = 50
)

Y <- mm1_sobol4r_model_ts(
  sob_mm1$X,
  horizon = 1000,
  warmup  = 200,
  nrep   = 20
)

if (any(is.na(Y))) stop("Model returned NA values")
if (any(is.infinite(Y))) stop("Model returned infinite values")
summary(Y)             # should now show positive, variable values
sob_mm1 <- tell(sob_mm1, Y)
```

```{r sobol-results, cache=TRUE, eval=LOCAL}
print(sob_mm1)
```

If `ggplot2` is installed, we can visualise the indices with
the autoplot method provided by the package.

```{r sobol-plot, message=FALSE, warning=FALSE, cache=TRUE, eval=LOCAL}
Sobol4R::autoplot(sob_mm1)
```

# Summary

We have shown how to

- build a simple M/M/1 queue in `simmer`  
- define a scalar summary statistic (mean waiting time) as a QoI  
- wrap the simulator into a function compatible with `sensitivity::sobol`  
- compute Sobol indices for the two key parameters `lambda` and `mu`  

The same pattern can be reused for more complex discrete event simulations
implemented with `simmer`.


```{r, cache=TRUE, eval=LOCAL}
n <- 200
X1 <- data.frame(
  lambda = runif(n, 0.05, 0.20),
  mu_reg = runif(n, 0.30, 0.80),
  mu_exam = runif(n, 0.20, 0.60)
)
X2 <- data.frame(
  lambda = runif(n, 0.05, 0.20),
  mu_reg = runif(n, 0.30, 0.80),
  mu_exam = runif(n, 0.20, 0.60)
)

sob_clinic <- sobol(model = NULL,
                    X1    = X1,
                    X2    = X2,
                    order = 2,
                    nboot = 50)

Yc <- sobol4r_clinic_model(sob_clinic$X,
                           cap_reg = 2, 
                           cap_exam = 3,
                           horizon = 2000,
                           warmup_prob = 0.2,
                           nrep = 10)
sob_clinic <- tell(sob_clinic, Yc)
print(sob_clinic)
Sobol4R::autoplot(sob_clinic)
```

