# pharmacist example
# opens 9 A.M.
# expects 32 prescriptions to serve daily before 5 P.M
# working on a prescription takes on average 10 mins, std. dev. 4 mins
# no new prescriptions after 5 P.M.
# stays in the shop for finishing prescriptions
# 1) what is the average time that he will depart his store at night?
# 2) what proportion of days will he still be working at 5:30 P.M.?
# 3) what is the average time it will take him to fill a prescription?
# 4) what proportion of prescriptions will be filled within 30 minutes?
# ** Let us assume to change acceptance policy: only accept new ones when less than 5 are on the stack
# 5) How many prescriptions are lost on average?
# 6) How answers to 1-4 change?

# working time -> 510 mins
# opening time -> 480 mins

# STEP 1 --> build a *probability model*
# Assumptions:
# A) Probabilistic model for the arrival of prescriptions? Constant interarrival? Changing according to the time of the day?
# Arrivals -> Poisson Process with rate (exp. number of events per unit of time) Lambda = 32 (per day)
# Interarrival times -> Exp. distribution with same rate (i.e. expected value 1 / Lambda)
get_next_delay <- function() {
  r = runif(1)
  return( - log(r) * 15)
}

# B) Probabilistic model for the service time? Does it depend also on how many prescriptions are in queue (or time of the day)?
# Service times -> Normal distribution mean 10, stdev 4
get_service_time <- function() {
  return( rnorm(1,m=10,sd=4) )
}
# C) Do these models change from one day to another (e.g. Monday, Tuesday ...) ?

# STEP 2 --> search for numerical solutions
# Analytical solution might be *too complex* --> SIMULATE!
# i.e. DRAW RANDOM NUMBERS to simulate possible outcomes on a specific (hypotetical) day

simulate <- function (max_time, accept_time) {
  
  # evaluation measures
  cum_prescription_time = 0;
  total_prescriptions = 0;
  fast_prescriptions = 0;
  
  # system status
  queue_length = 0;
  busy = FALSE;
  events_time = c();
  events_type = c();
  events_start = c();
  
  # current event status
  current_time = get_next_delay();
  current_type = 'A';
  current_start = 0;
  
  # iteratively, put events to a queue, retrieve the earliest one and handle it
  while (current_time < max_time) {
    
    # What if acceptance policy is "no queue limits"?
    if (current_type == 'A' && current_time < accept_time) {
    # What if acceptance policy is "up to 5 in queue"?
    # if (current_type == 'A' && current_time < accept_time && queue_length < 5) {

      total_prescriptions = total_prescriptions + 1
      
      if (busy) {
        queue_length = queue_length + 1;
      } else {
        serviced_event_time = get_service_time() + current_time
        events_time = c(events_time, serviced_event_time)
        events_type = c(events_type, 'S')
        events_start = c(events_start, current_time)
        cum_prescription_time = cum_prescription_time + (serviced_event_time - current_time)
        if (serviced_event_time - current_time <= 30) fast_prescriptions = fast_prescriptions + 1;
        busy = TRUE
      }
      
      arrival_event_time = get_next_delay()
      events_time = c(events_time,  current_time + arrival_event_time)
      events_type = c(events_type, 'A')
      # just for symmetry
      events_start = c(events_start, current_time)
      
    }
    if (current_type == 'S') {
      
      if (queue_length == 0) {
        busy = FALSE
      } else {
        queue_length = queue_length - 1
        serviced_event_time = get_service_time() + current_time
        events_time = c(events_time, serviced_event_time)
        events_type = c(events_type, 'S')
        events_start = c(events_start, current_time)
        cum_prescription_time = cum_prescription_time + (serviced_event_time - current_start)
        if (serviced_event_time - current_start <= 30) fast_prescriptions = fast_prescriptions + 1;
        busy = TRUE
      }
      
    }

    if (length(events_time) == 0) 
      # 1) what is the average time that he will depart his store at night?
      return(max(current_time, accept_time));
      # return(current_time);
      # 2) what proportion of days will he still be working at 5:30 P.M.?
      # return(FALSE);
      # 3) what is the average time it will take him to fill a prescription?
      # return (cum_prescription_time / total_prescriptions)
      # 4) what proportion of prescriptions will be filled within 30 minutes?
      # return (fast_prescriptions / total_prescriptions)
    
    

    k = which.min(events_time);
    current_time = events_time[k]
    current_type = events_type[k]
    current_start = events_start[k]
    events_time = events_time[-k]
    events_type = events_type[-k]
    events_start = events_start[-k]
  }
  # 1) what is the average time that he will depart his store at night?
  return (current_time)
  # 2) what proportion of days will he still be working at 5:30 P.M.?
  # return (TRUE)
  # 3) what is the average time it will take him to fill a prescription?
  # return (cum_prescription_time / total_prescriptions);
  # 4) what proportion of prescriptions will be filled within 30 minutes?
  # return (fast_prescriptions / total_prescriptions)
  
}

# run the test and observe how the results are distributed
r = c(); for (i in 1:10000) r = c(r, simulate(510, 480))
# hist(r, breaks = 100)
# compute average results on different tests, and observe how the average results are distributed
v = c(); for (k in 0:199) v = c(v,mean(r[(k*50):((k+1)*50)]));
hist(v, breaks = 15)