## ---- cache = FALSE, include=FALSE---------------------------------------
knitr::opts_chunk$set(collapse = T, comment = "#>", 
                      fig.width = 6, fig.height = 4, fig.align = "center")

required <- c("simmer.plot", "dplyr")

if (!all(unlist(lapply(required, function(pkg) requireNamespace(pkg, quietly = TRUE)))))
  knitr::opts_chunk$set(eval = FALSE)

## ---- message=FALSE------------------------------------------------------
library(simmer)
library(simmer.plot)
library(parallel)
set.seed(1234)

## ------------------------------------------------------------------------
lambda <- 2
mu <- 4
rho <- lambda/mu # = 2/4

mm1.trajectory <- trajectory() %>%
  seize("resource", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("resource", amount=1)

mm1.env <- simmer() %>%
  add_resource("resource", capacity=1, queue_size=Inf) %>%
  add_generator("arrival", mm1.trajectory, function() rexp(1, lambda)) %>% 
  run(until=2000)

## ------------------------------------------------------------------------
# Theoretical value
mm1.N <- rho/(1-rho)

# Evolution of the average number of customers in the system
plot(mm1.env, "resources", "usage", "resource", items="system") +
  geom_hline(yintercept=mm1.N)

## ------------------------------------------------------------------------
plot(mm1.env, "resources", "usage", "resource", items=c("queue", "server"), steps=TRUE) +
  xlim(0, 20) + ylim(0, 4)

## ------------------------------------------------------------------------
mm1.arrivals <- get_mon_arrivals(mm1.env)
mm1.t_system <- mm1.arrivals$end_time - mm1.arrivals$start_time

mm1.T <- mm1.N / lambda
mm1.T ; mean(mm1.t_system)

## ---- eval=FALSE---------------------------------------------------------
#  envs <- mclapply(1:100, function(i) {
#    simmer() %>%
#      add_resource("resource", capacity=1, queue_size=Inf) %>%
#      add_generator("arrival", mm1.trajectory, function() rexp(100, lambda)) %>%
#      run(1000/lambda) %>%
#      wrap()
#  }, mc.set.seed=FALSE)

## ---- eval=FALSE---------------------------------------------------------
#  t_system <- get_mon_arrivals(envs) %>%
#    dplyr::mutate(t_system = end_time - start_time) %>%
#    dplyr::group_by(replication) %>%
#    dplyr::summarise(mean = mean(t_system))
#  
#  t.test(t_system$mean)
#  #>
#  #> 	One Sample t-test
#  #>
#  #> data:  t_system$mean
#  #> t = 112.35, df = 99, p-value < 2.2e-16
#  #> alternative hypothesis: true mean is not equal to 0
#  #> 95 percent confidence interval:
#  #>  0.4883420 0.5059016
#  #> sample estimates:
#  #> mean of x
#  #> 0.4971218

## ------------------------------------------------------------------------
lambda; 1/mean(diff(subset(mm1.arrivals, finished==TRUE)$start_time))

## ------------------------------------------------------------------------
qqplot(mm1.t_system, rexp(1000, 1/mm1.T))
abline(0, 1, lty=2, col="red")

## ------------------------------------------------------------------------
lambda <- 2
mu <- 4

mm23.trajectory <- trajectory() %>%
  seize("server", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("server", amount=1)

mm23.env <- simmer() %>%
  add_resource("server", capacity=2, queue_size=1) %>%
  add_generator("arrival", mm23.trajectory, function() rexp(1, lambda)) %>%
  run(until=2000)

## ------------------------------------------------------------------------
mm23.arrivals <- get_mon_arrivals(mm23.env)
mm23.arrivals %>%
  dplyr::summarise(rejection_rate = sum(!finished)/length(finished))

## ------------------------------------------------------------------------
mm23.t_system <- mm23.arrivals$end_time - mm23.arrivals$start_time
# Comparison with M/M/1 times
qqplot(mm1.t_system, mm23.t_system)
abline(0, 1, lty=2, col="red")

## ------------------------------------------------------------------------
mean_pkt_size <- 100        # bytes
lambda1 <- 2                # pkts/s
lambda3 <- 0.5              # pkts/s
lambda4 <- 0.6              # pkts/s
rate <- 2.2 * mean_pkt_size # bytes/s

# set an exponential message size of mean mean_pkt_size
set_msg_size <- function(.)
  set_attribute(., "size", function() rexp(1, 1/mean_pkt_size))

# seize an M/D/1 queue by id; the timeout is function of the message size
md1 <- function(., id)
  seize(., paste0("md1_", id), 1) %>%
  timeout(function(attrs) attrs[["size"]] / rate) %>%
  release(paste0("md1_", id), 1)

## ------------------------------------------------------------------------
to_queue_1 <- trajectory() %>%
  set_msg_size() %>%
  md1(1) %>%
  leave(0.25) %>%
  md1(2) %>%
  branch(function() (runif(1) > 0.65) + 1, c(F, F),
         trajectory() %>%
           md1(3),
         trajectory() %>%
           md1(4))

to_queue_3 <- trajectory() %>%
  set_msg_size() %>%
  md1(3)

to_queue_4 <- trajectory() %>%
  set_msg_size() %>%
  md1(4)

## ------------------------------------------------------------------------
env <- simmer()
for (i in 1:4) env %>% 
  add_resource(paste0("md1_", i))
env %>%
  add_generator("arrival1_", to_queue_1, function() rexp(1, lambda1), mon=2) %>%
  add_generator("arrival3_", to_queue_3, function() rexp(1, lambda3), mon=2) %>%
  add_generator("arrival4_", to_queue_4, function() rexp(1, lambda4), mon=2) %>%
  run(4000)

## ------------------------------------------------------------------------
res <- get_mon_arrivals(env, per_resource = TRUE) %>%
  dplyr::select(name, resource) %>%
  dplyr::filter(resource %in% c("md1_3", "md1_4"))
arr <- get_mon_arrivals(env) %>%
  dplyr::mutate(waiting_time = end_time - (start_time + activity_time),
                generator = regmatches(name, regexpr("arrival[[:digit:]]", name))) %>%
  dplyr::left_join(res) %>%
  dplyr::group_by(generator, resource)

dplyr::summarise(arr, average = sum(waiting_time) / n())
get_n_generated(env, "arrival1_") + get_n_generated(env, "arrival4_")
dplyr::count(arr)

