## ----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-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)

## ----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, 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
  })
}

## ----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-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)

## ----sobol-results, cache=TRUE, eval=LOCAL------------------------------------
print(sob_mm1)

## ----sobol-plot, message=FALSE, warning=FALSE, cache=TRUE, eval=LOCAL---------
Sobol4R::autoplot(sob_mm1)

## ----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)