#
# file: Assignment3_Rsolution.R
#
# Author: B. M. Bolstad
# Date: Oct 16, 2004
#
# Aim: More sophisticated queuing simulation code. Showing
#      a more complete solution than original code.
#
# Instructions: Copy and paste.
#
#
#

queue.simulation <- function(rate,ncustomers=1000,do.plots=TRUE,smooth.frac=0.1,do.summary.stats=TRUE,sim.data=FALSE){
# simulate the service times
  servicetimes <- rnorm(ncustomers,1.8,0.25)  # ie S
  arrivaltimes <- rexp(ncustomers,rate/60)
  
                                        # figure out when customers actually arrive
  cum.arrival <- cumsum(arrivaltimes)

# each customer could possibly be served as soon as the server
# becomes free. First customer is always served immediately.
# Otherwise the first possible time that they could start being served is
# the maximum of their arrival time and the time at which the the
# previous customer finished being served.


  start.servicetimes <- matrix(0,1000,1)
  
  start.servicetimes[1] <- cum.arrival[1]
  for (i in 2:1000){
    start.servicetimes[i] <- max(cum.arrival[i], start.servicetimes[i-1]+ servicetimes[i-1])
  }
  



# The time that the customer leaves the store is when they
# complete being served. which is the time that the start being served plus
# their service time

  leave.times <- start.servicetimes + servicetimes
  
  

#
# Queue length changes whenever anyone arrives or leaves
# Count the person being served as being a member of the
# queue. Queue starts out empty.
#

  how.many.changes <- 2*1000 + 1 #every customer joins and leaves plus the starting time
  
  queue.length <- matrix(0,how.many.changes,2) # store both times and lengths
  
  queue.length[1,] <- c(0,0)  # empty queue at beginning
  change.times <- sort(c(cum.arrival,leave.times))
  queue.length[2:how.many.changes,1] <- change.times
  

  for (i in 2:how.many.changes){
    if (is.element(change.times[i-1],cum.arrival)){
      queue.length[i,2] <- queue.length[i-1,2] +1
    } else {
      queue.length[i,2] <- queue.length[i-1,2] -1
    }   
  }
  
  colnames(queue.length) <- c("Time","Length")
  



#
# Helpful hints:
#
# 1. W is the difference between the starting service time and arrival time
# 2. T is the difference between the leaving time and arrival time
# 3. The period of time between when the previous customer left
#    and service for the next customer begins gives the idle time I.

  W <- start.servicetimes - cum.arrival
  Total.time <- leave.times - cum.arrival
  idle.time <- start.servicetimes - c(0,leave.times[1:(ncustomers-1)])

  if (do.plots){
    hist(W,prob=TRUE,main="Histogram (Smoothed) of time spent in queue")
    lines(density(W))
    plot(W,xlab="Customer Number",ylab="Total Time",main="Time spent in queue versus Customer",type="l")
    lines(lowess(1:ncustomers,W,f=smooth.frac),lwd=2,col="red")
    
    hist(Total.time,prob=TRUE,main="Histogram (Smoothed) of Total Time")
    lines(density(Total.time))
    plot(Total.time,xlab="Customer Number",ylab="Total Time",main="Total Time versus Customer",type="l")
    lines(lowess(1:ncustomers,Total.time,f=smooth.frac),lwd=2,col="red")
    
    plot(queue.length,ylab="Queue Length",xlab="Time",main="Queue Length versus Time",type="l")
    lines(lowess(queue.length,f=smooth.frac),lwd=2,col="red")

    plot(idle.time,type="l")
    lines(lowess(1:ncustomers,idle.time,f=smooth.frac),lwd=2,col="red")
    plot(cumsum(idle.time),type="l",xlab="Customer Number",ylab="Idle time",main="Cumulative Idle time")
  }

  summary.stats <- matrix(0,4,7)
  rownames(summary.stats) <- c("W","T","Idle time","Queue Length")
  colnames(summary.stats) <- c("Mean","Median","Stdev","Max","Min","LQ","UQ")
  summary.stats[1,1] <- mean(W)
  summary.stats[1,2] <- median(W)
  summary.stats[1,3] <- sd(W)
  summary.stats[1,4] <- max(W)
  summary.stats[1.5] <- min(W)
  summary.stats[1,6:7] <- quantile(W,c(0.25,0.75))
  summary.stats[2,1] <- mean(Total.time)
  summary.stats[2,2] <- median(Total.time)
  summary.stats[2,3] <- sd(Total.time)
  summary.stats[2,4] <- max(Total.time)
  summary.stats[2.5] <- min(Total.time)
  summary.stats[2,6:7] <- quantile(Total.time,c(0.25,0.75))
  summary.stats[3,1] <- mean(idle.time)
  summary.stats[3,2] <- median(idle.time)
  summary.stats[3,3] <- sd(idle.time)
  summary.stats[3,4] <- max(idle.time)
  summary.stats[3.5] <- min(idle.time)
  summary.stats[3,6:7] <- quantile(idle.time,c(0.25,0.75))
   summary.stats[4,1] <- mean( queue.length[,2])
  summary.stats[4,2] <- median( queue.length[,2])
  summary.stats[4,3] <- sd( queue.length[,2])
  summary.stats[4,4] <- max( queue.length[,2])
  summary.stats[4.5] <- min( queue.length[,2])
  summary.stats[4,6:7] <- quantile( queue.length[,2],c(0.25,0.75))
    
  if (do.summary.stats){
    print(summary.stats)
  }

  if (sim.data){
    X <- cbind(servicetimes, cum.arrival, start.servicetimes,leave.times,W,Total.time)
    colnames(X) <- c("Service Times","Arrival Time","Start Service Time","Leave Time","Waiting Time","Total Time Spent")
    rownames(X) <- 1:ncustomers


    list(times=X,queuelength=queue.length)
    
  }
  
  
}

#par(ask=FALSE)
par(ask=TRUE)
queue.simulation(rate=25)
par(ask=FALSE)
results <- queue.simulation(rate=25,sim.data=TRUE)


# you could later take the output stored in results and do further analysis