# DES.r:  R routines for discrete-event simulation (DES), with an example

# each event will be represented by a vector; the first component will
# be the time the event is to occur; the second component will be the
# numerical code for the programmer-defined event type; the programmer
# may add application-specific components

# a list named "sim" holds the events list and other information; for
# convenience, sim has been stored as a global variable; some functions
# have side effects

# create "sim"
newsim <- function(numfields) {
   sim <<- list()
   sim$currtime <<- 0.0  # current simulated time
   sim$evnts <<- NULL  # event list
}

# insert event evnt into event list
insevnt <- function(evnt) {
   if (is.null(sim$evnts)) {
      sim$evnts <<- matrix(evnt,nrow=1)
      return()
   }
   # find insertion point
   inspt <- binsearch(sim$evnts[,1],evnt[1])
   # now "insert"
   if (inspt > 1) e <- rbind(sim$evnts[1:(inspt-1),],evnt)
   nrse <- nrow(sim$evnts) 
   if (inspt <= nrse)
      e <- rbind(evnt, sim$evnts[inspt:nrse,])
   sim$evnts <<- e
}

# schedule new event; evnttime is the time at which the event is to
# occur; evnttype is the event type; and appfields are the values of the
# programmer-defined fields, if any
schedevnt <- function(evnttime,evnttype,appfields=NULL) {
   evnt <- c(evnttime,evnttype,appfields)
   insevnt(evnt)  
}

# start to process next event (second half done by application
# programmer via call to reactevnt())
getnextevnt <- function() {
   head <- sim$evnts[1,]
   # delete head
   if (nrow(sim$evnts) == 1) sim$evnts <<- NULL else 
      sim$evnts <<- sim$evnts[-1,,drop=F]
   return(head)
}

# main loop of the simulation
mainloop <- function(maxsimtime) {
   while(sim$currtime < maxsimtime) {
      head <- getnextevnt()
      sim$currtime <<- head[1]  # update current simulated time
      reactevnt(head)  # process this event (programmer-supplied ftn)
   }
}

# binary search of insertion point of y in the sorted vector x; returns
# the position in x before which y should be inserted, with the value
# length(x)+1 if y is larger than x[length(x)]
binsearch <- function(x,y) {
   n <- length(x)
   lo <- 1
   hi <- n
   while(lo+1 < hi) {
      mid <- floor((lo+hi)/2)
      if (y == x[mid]) return(mid)
      if (y < x[mid]) hi <- mid else lo <- mid
   }
   if (y <= x[lo]) return(lo)
   if (y < x[hi]) return(hi)
   return(hi+1)
}

# application:  M/M/1 queue, arrival rate 0.5, service rate 1.0

# globals
# rates
arvrate <- 0.5
srvrate <- 1.0
# event types
arvtype <- 1
srvdonetype <- 2
# initialize server queue
srvq <- NULL  # will just consist of arrival times of queued jobs
# statistics
njobsdone <- NULL  # jobs done so far
totwait <- NULL  # total wait time so far

# event processing function requried by general DES code above
reactevnt <- function(head) {
   if (head[2] == arvtype) {  # arrival
      # if server free, start service, else add to queue
      if (length(srvq) == 0) {
         srvq <<- head[3]
         srvdonetime <- sim$currtime + rexp(1,srvrate)
         schedevnt(srvdonetime,srvdonetype,head[3])
      } else srvq <<- c(srvq,head[3])
      # generate next arrival
      arvtime <- sim$currtime + rexp(1,arvrate)
      schedevnt(arvtime,arvtype,arvtime)
   } else {  # service done
      # process job that just finished
      # do accounting
      njobsdone <<- njobsdone + 1
      totwait <<- totwait + sim$currtime - head[3]
      # remove from queue
      srvq <<- srvq[-1]
      # more still in the queue?
      if (length(srvq) > 0) {
         # schedule new service
         srvdonetime <- sim$currtime + rexp(1,srvrate)
         schedevnt(srvdonetime,srvdonetype,srvq[1])
      }
   }
}

mm1sim <- function() {
   srvq <<- vector(length=0)  
   njobsdone <<- 0
   totwait <<- 0.0
   # create simulation, 1 extra field (arrival time)
   newsim(1)
   # get things going, with first arrival event
   arvtime <- rexp(1,rate=arvrate)
   schedevnt(arvtime,arvtype,arvtime)
   mainloop(100000.0)
   return(totwait/njobsdone)
}