# Three phase simulation of an accident and emergency clinic
rm(list=ls())
set.seed(600617)

# Inputs
NP.num <- 1  # number of nurse-practicioners
N.num <- 2   # number of nurses
D.num <- 2   # number of doctors
t.end <- 5000  # simulation end time in hours
fname <- sprintf("emergency-%d-%d-%d.txt",NP.num, N.num, D.num) # file output
reporting <- TRUE  # flag for diagnostic reporting

# Functions giving activity durations
Arr.time <- function() rexp(1, rate=4)  # mean 1/4
NP.time <- function() rgamma(1, shape=2, rate=12)  # mean 1/6
NC.time <- function() rgamma(1, shape=2, rate=12)  # mean 1/6
D.time <- function() rgamma(1, shape=3, rate=9)  # mean 1/3
NT.time <- function() rgamma(1, shape=3, rate=12) # mean 1/4


# Event functions

Arr <- function(t, entities, resources, e) {
  # patient arrival
  # put them in the nurse-practitioner queue
  resources[4] <- resources[4] + 1
  entities$available[e] <- TRUE
  entities$eventtime[e] <- NA
  entities$nextevent[e] <- NA
  entities$activity[e] <- "NPQ"
  entities$queue_position[e] <- resources[4]
  entities$arrival_time[e] <- t
  # create next patient
  entities <- rbind(list(FALSE, t + Arr.time(), "Arr", "Arr", NA, NA),
                    entities)
  return(list(entities=entities, resources=resources))
}

NPstart <- function(t, entities, resources, e=NULL) {
  # start nurse-practitioner classification
  if (resources[1] > 0 && resources[4] > 0) {
    # update nurse-practitioner queue positions
    entities$queue_position[entities$activity == "NPQ"] <-
      entities$queue_position[entities$activity == "NPQ"] - 1
    # j is head of nurse-practitioner queue
    j <- which(entities$activity == "NPQ" & entities$queue_position == 0)
    # update patient entity
    entities$available[j] <- FALSE
    entities$eventtime[j] <- t + NP.time()
    entities$nextevent[j] <- "NPend"
    entities$activity[j] <- "NP"
    entities$queue_position[j] <- NA
    # update resource numbers
    resources[1] <- resources[1] - 1   # one less NP
    resources[4] <- resources[4] - 1   # one less NPQ
    changed <- TRUE
  } else {
    changed <- FALSE
  }
  return(list(entities=entities, resources=resources, changed=changed))
}

NPend <- function(t, entities, resources, e) {
  # end nurse-practitioner classification
  resources[1] <- resources[1] + 1
  p <- runif(1)
  if (p < 0.3) {
    # send to nurse classification queue
    resources[5] <- resources[5] + 1
    entities$available[e] <- TRUE
    entities$eventtime[e] <- NA
    entities$nextevent[e] <- NA
    entities$activity[e] <- "NCQ"
    entities$queue_position[e] <- resources[5]
  } else if (p < 0.7) {
    # send to doctor queue                    
    resources[6] <- resources[6] + 1
    entities$available[e] <- TRUE
    entities$eventtime[e] <- NA
    entities$nextevent[e] <- NA
    entities$activity[e] <- "DQ"
    entities$queue_position[e] <- resources[6]
  } else {
    # send home, don't record time in system
    entities <- entities[-e,]
  }
  return(list(entities=entities, resources=resources))
}

NCstart <- function(t, entities, resources, e=NULL) {
  # start nurse classification
  if (resources[2] > 0 && resources[5] > 0) {
    # update nurse classification queue positions
    entities$queue_position[entities$activity == "NCQ"] <-
      entities$queue_position[entities$activity == "NCQ"] - 1
    # j is head of nurse classification queue
    j <- which(entities$activity == "NCQ" & entities$queue_position == 0)
    # update patient entity
    entities$available[j] <- FALSE
    entities$eventtime[j] <- t + NC.time()
    entities$nextevent[j] <- "NCend"
    entities$activity[j] <- "NC"
    entities$queue_position[j] <- NA
    # update resource numbers
    resources[2] <- resources[2] - 1   # one less N
    resources[5] <- resources[5] - 1   # one less NCQ
    changed <- TRUE
  } else {
    changed <- FALSE
  }
  return(list(entities=entities, resources=resources, changed=changed))
}

NCend <- function(t, entities, resources, e) {
  # end nurse classification
  resources[2] <- resources[2] + 1
  p <- runif(1)
  if (p < 0.1) {
    # send to doctor queue
    resources[6] <- resources[6] + 1
    entities$available[e] <- TRUE
    entities$eventtime[e] <- NA
    entities$nextevent[e] <- NA
    entities$activity[e] <- "DQ"
    entities$queue_position[e] <- resources[6]
  } else {
    # send home, record time in system
    cat(file=fname, append=TRUE, t - entities$arrival_time[e], "\n")
    entities <- entities[-e,]
  }
  return(list(entities=entities, resources=resources))
}

Dstart <- function(t, entities, resources, e=NULL) {
  # start doctor treatment
  if (resources[3] > 0 && resources[6] > 0) {
    # update doctor queue positions
    entities$queue_position[entities$activity == "DQ"] <-
      entities$queue_position[entities$activity == "DQ"] - 1
    # j is head of doctor queue
    j <- which(entities$activity == "DQ" & entities$queue_position == 0)
    # update patient entity
    entities$available[j] <- FALSE
    entities$eventtime[j] <- t + D.time()
    entities$nextevent[j] <- "Dend"
    entities$activity[j] <- "D"
    entities$queue_position[j] <- NA
    # update resource numbers
    resources[3] <- resources[3] - 1   # one less D
    resources[6] <- resources[6] - 1   # one less DQ
    changed <- TRUE
  } else {
    changed <- FALSE
  }
  return(list(entities=entities, resources=resources, changed=changed))
}

Dend <- function(t, entities, resources, e) {
  # end doctor treatment, send to nurse treatment queue
  resources[7] <- resources[7] + 1     # one more NTQ
  resources[3] <- resources[3] + 1     # one more D
  entities$available[e] <- TRUE
  entities$eventtime[e] <- NA
  entities$nextevent[e] <- NA
  entities$activity[e] <- "NTQ"
  entities$queue_position[e] <- resources[7]
  return(list(entities=entities, resources=resources))
}

NTstart <- function(t, entities, resources, e=NULL) {
  # start nurse treatment
  if (resources[2] > 0 && resources[7] > 0) {
    # update nurse treatment queue positions
    entities$queue_position[entities$activity == "NTQ"] <-
      entities$queue_position[entities$activity == "NTQ"] - 1
    # j is head of nurse treatment queue
    j <- which(entities$activity == "NTQ" & entities$queue_position == 0)
    # update patient entity
    entities$available[j] <- FALSE
    entities$eventtime[j] <- t + NT.time()
    entities$nextevent[j] <- "NTend"
    entities$activity[j] <- "NT"
    entities$queue_position[j] <- NA
    # update resource numbers
    resources[2] <- resources[2] - 1   # one less N
    resources[7] <- resources[7] - 1   # one less NTQ
    changed <- TRUE
  } else {
    changed <- FALSE
  }
  return(list(entities=entities, resources=resources, changed=changed))
}

NTend <- function(t, entities, resources, e) {
  # end nurse treatment, send home, record time in system
  resources[2] <- resources[2] + 1
  cat(file=fname, append=TRUE, t - entities$arrival_time[e], "\n")
  entities <- entities[-e,]
  return(list(entities=entities, resources=resources))
}

C.list <- c("NPstart", "NCstart", "Dstart", "NTstart")

# Initialise entities, resources and time
# entities are all patients, with these attributes:
#   activity in {Arr, NPQ, NP, NCQ, NC, DQ, D, NTQ, NT}
#   queue_position if in {NPQ, NCQ, NTQ}
#   and arrival_time
# resources are NP, N, D, NPQ, NCQ, DQ, NTQ
t <- 0
entities <- data.frame(available=FALSE,
                       eventtime=Arr.time(),
                       nextevent="Arr",
                       activity="Arr",
                       queue_position=NA,
                       arrival_time=NA,
                       stringsAsFactors=FALSE)
resources <- c(NP.num, N.num, D.num, 0, 0, 0, 0)


print.entities.resources <- function(t, entities, resources) {
  # print current system state
  cat(t, "\nNPQ  NP NCQ  NC  DQ   D NTQ  NT\n")
  x <- sum(entities$activity == "NC")
  cat(format(resources[4], width=3), 
      format(NP.num-resources[1], width=3),
      format(resources[5], width=3), 
      format(x, width=3), 
      format(resources[6], width=3), 
      format(D.num-resources[3], width=3),
      format(resources[7], width=3), 
      format(N.num-resources[2]-x, width=3), "\n")
  cat('\n')
}

# Control loop
cat(file=fname, "time in system\n")  # prepare output file
if (reporting) print.entities.resources(t, entities, resources)
while (t < t.end) {
  # A phase
  t <- min(entities$eventtime[!entities$available])
  B.index <- sort(which(entities$eventtime == t & !entities$available), decreasing=TRUE)
  # B phase
  for (i in B.index) {
    B <- entities$nextevent[i]
    x <- do.call(B, list(t, entities, resources, i))
    entities <- x$entities
    resources <- x$resources
    if (reporting) print.entities.resources(t, entities, resources)
  }
  # C phase
  C.activity <- TRUE
  while (C.activity) {
    C.activity <- FALSE
    for (C in C.list) {
      x <- do.call(C, list(t, entities, resources))
      if (x$changed) {
        entities <- x$entities
        resources <- x$resources
        C.activity <- TRUE
        if (reporting) print.entities.resources(t, entities, resources)
      }
    }
  }
}   

# Output
x <- scan(fname, skip=1)
hist(x)