#### Code by Matt Asher. Published at StatisticsBlog.com ####

#### CONFIG ####
# Where do you want the animation stored?
setwd("path/to/directory/")

# Number of slots to fill
numbSlots = 5

# Total time intervals in the simulation
intervals = 250

# Probability that a new person will show up during an interval
# Note, a maximum of one new person can show up during an interval
p = 0.1

# Average time each person takes at the teller, using discretized exponential distribution 
# Times will be augmented by one, so that everyone takes at least 1 interval to serve
meanServiceTime = 60

#### INITIALIZATION ####
queueLengths = rep(0, intervals)
slots = list()
waitTimes = c()
leavingTimes = c()
queue = list()
arrivalTimes = c()
frontOfLineWaits = c()
id = 1

#### Libraries ####
# Use the proto library to treat people like objects in traditional oop
library("proto")

#### Functions ####
# R is missing a nice way to do ++, so we use this
inc <- function(x) {
  eval.parent(substitute(x <- x + 1))
}

# Main object, really a "proto" function
# Changed based on suggestion by Johann (http://www.statisticsblog.com/2011/10/waiting-in-line-waiting-on-r/comment-page-1/#comment-12272)
person <- proto(
	new = function (.,id) {
		# Set id and how muct teller time this person will need
		.$proto(intervalsNeeded = floor(rexp(1, 1/meanServiceTime)) + 1,personId = id)
	},
	
	intervalArrived = 0,
	intervalAttended = NULL,
	
	intervalsWaited = 0,
	intervalsWaitedAtHeadOfQueue = 0,
)

# Fill our slots with Nullos (don't Google that term, NSFL)
nullPerson = person$new(0)

for(k in 1:numbSlots) {
	slots[[k]] = nullPerson
}

#### Main loop ####
for(i in 1:intervals) {
	# Check if anyone is leaving the slots
	for(j in 1:numbSlots) {
	
		if(slots[[j]]$personId) {
			if( (slots[[j]]$intervalAttended + slots[[j]]$intervalsNeeded) == i) {
				# They are leaving the queue, slot to 0
				slots[[j]] = nullPerson
				leavingTimes = c(leavingTimes, i)
			}
		}
	}

	# See if a new person is to be added to the queue
	if(runif(1) < p) {
		newPerson = person$new(id)
		inc(id)
		newPerson$intervalArrived = i
		queue = c(queue, newPerson)
		arrivalTimes  = c(arrivalTimes, i)
	}

	# Can we place someone into a slot?
	for(j in 1:numbSlots) {
		# Is this slot free
		if(!slots[[j]]$personId) {
			if(length(queue) > 0) {

				slots[[j]] = queue[[1]]
				slots[[j]]$intervalAttended = i
				waitTimes = c(waitTimes, queue[[1]]$intervalsWaited)
				
				# Only interested in these if person waited 1 or more intevals at front of line
				if(queue[[1]]$intervalsWaitedAtHeadOfQueue) {
					frontOfLineWaits = c(frontOfLineWaits, queue[[1]]$intervalsWaitedAtHeadOfQueue)
				}

				# Remove placed person from queue
				queue[[1]] = NULL
			}
		}
	}

	# Everyone left in the queue has now waited one more interval to be served
	if(length(queue)) {
		for(j in 1:length(queue)) {
			inc(queue[[j]]$intervalsWaited)
		}

		# The (possibley new) person at the front of the queue has had to wait there one more interval.
		inc(queue[[1]]$intervalsWaitedAtHeadOfQueue) 
	}

	# End of the interval, what is the state of things
	queueLengths[i] = length(queue);
	
	# Create a snapshot of this interval.
	png(file=paste("snapshot-",paste(sprintf( "%04d", i),sep="",collapse=""),".png",sep=""), width=550, height=550)
	plot.new()
	plot(xlim=c(0, numbSlots), ylim=c(0, 20), 1, 1, pch=18, col=(slots[[1]]$personId %% 153), cex=sqrt((slots[[1]]$intervalsNeeded) ), main=paste("Interval ",i,sep=""), xlab="Queueing area", ylab="Slots")

	# Add the rest of the slots
	if(numbSlots > 1) {
		for(sl in 1:numbSlots) {
			points(sl,1, pch=18, col=(slots[[sl]]$personId %% 153), cex=log(slots[[sl]]$intervalsNeeded+1))
		}
	}
	
	# Show the line, if there is one
	if(length(queue)) {
		for(j in 1:length(queue)) {
			points((numbSlots/2),(j+1), pch=18, col=(queue[[j]]$personId %% 153), cex=log(queue[[j]]$intervalsNeeded+1))
		}
	}
	
	dev.off()
}

# Create an animated gif of all of this. You will need ImageMagick installed for this
system('"path/to/convert" -delay 15 snapshot-*.png animation.gif')

#### Output ####
# par(ask=T) # So we can see plots one by one
plot(queueLengths, type="o", col="blue", pch=20, main="Queue lengths over time", xlab="Interval", ylab="Queue length")
#plot(waitTimes, type="o", col="blue", pch=20, main="Wait times", xlab="Person", ylab="Wait time")
#if(length(frontOfLineWaits)) { hist(frontOfLineWaits, col="blue", breaks=10) }

# Optional removal of all snapshot files. Careful with this.
# file.remove(list.files(pattern=".png"))