# queue.R
#
# Usage:
# > MM1(EX=1, ES=1)
# > MG1(EX=1)

service = function(task, xnew=Inf){
  stamp = cumsum(task)
  done = sum(stamp < xnew)  # the number of served tasks before the next customer arrives
  if(done > 0){
    task = diff(c(xnew, stamp[-(1:done)]))  # remaining tasks
    stamp = stamp[1:done]
  }else{
  	task[1] = task[1] - xnew
  	stamp = c()
  }
  list(done = done, stamp = stamp, task = task)
}

# test:
# > service(c(1,1,3,8,6,5,6), xnew=9)

MM1 = function(EX=1, ES=1, tmax=100, plt=TRUE){
  # M/M/1, based on the definition

  lam = 1.0 / EX
  mu = 1.0 / ES
  history = matrix(0, 1, 2)
  t = rexp(1, lam)  # arriving time of 1st customer
  S = rexp(1, mu)  # service time
  task = c(S)  # "queue"
  Q = 1  # the number of waiting customers
  while(t < tmax){
  	history = rbind(history, c(t, Q))
    xnew = rexp(1, lam)  # next customer
    
    ser = service(task, xnew)  # service

    task = ser$task
    m = length(ser$stamp)
    if(m > 0){
      h1 = t + ser$stamp
      h2 = Q - (1:m)
      history = rbind(history, cbind(h1, h2))
      Q = Q - m
    }
    t = t + xnew
    S = rexp(1, mu)
    task = c(task, S)
    Q = Q + 1
  }
  if(plt){
    plot(history, type="s", xlab="time", ylab="queue")
  }
  history
}

MM1.check = function(EX=1, ES=1, tmax=100, plt=TRUE){
  # M/M/1, based on the birth-death process
  
  lam = 1.0 / EX
  mu = 1.0 / ES
  history = matrix(0, 1, 2)
  t = 0
  Q = 0
  while(t < tmax){
    if(Q == 0){
      xnew = rexp(1, lam)
      t = t + xnew
      Q = Q + 1
    }else{
      event = rexp(1, lam + mu)
      t = t + event
      U = runif(1)
      if(U < lam / (lam + mu)){
        Q = Q + 1
      }else{
        Q = Q - 1
      }
    }
    history = rbind(history, c(t, Q))
  }
  if(plt){
    plot(history, type="s", xlab="time", ylab="queue")
  }
  history
}

MG1 = function(EX=1, rg=function() 1, tmax=100, plt=TRUE){
  # M/G/1, based on the definition

  lam = 1.0 / EX
  history = matrix(0, 1, 2)
  t = rexp(1, lam)  # arriving time of 1st customer
  S = rg()  # service time
  task = c(S)  # "queue"
  Q = 1  # the number of waiting customers
  while(t < tmax){
  	history = rbind(history, c(t, Q))
    xnew = rexp(1, lam)  # next customer
    
    ser = service(task, xnew)  # service

    task = ser$task
    m = length(ser$stamp)
    if(m > 0){
      h1 = t + ser$stamp
      h2 = Q - (1:m)
      history = rbind(history, cbind(h1, h2))
      Q = Q - m
    }
    t = t + xnew
    S = rg()
    task = c(task, S)
    Q = Q + 1
  }
  if(plt){
    plot(history, type="s", xlab="time", ylab="queue")
  }
  history
}

test1 = function(){
  set.seed(20170629)
  par(cex=1.5, mar=c(5,4,2,1), cex=1.5, lwd=2)
  hist1 = MM1(1, 0.5, 50, TRUE)
  dev.copy2eps(file="queue_MM1_1.eps")
  hist2 = MM1(1, 1.0, 50, TRUE)
  dev.copy2eps(file="queue_MM1_2.eps")
  hist3 = MM1(1, 1.5, 50, TRUE)
  dev.copy2eps(file="queue_MM1_3.eps")

  # MM1.check(1, 0.5, 50, TRUE)
  # MM1.check(1, 1.0, 50, TRUE)
  # MM1.check(1, 1.5, 50, TRUE)

  MG1(EX=1.0/0.5, tmax=50)
  dev.copy2eps(file="queue_MG1_1.eps")
  MG1(EX=1.0, tmax=50)
  dev.copy2eps(file="queue_MG1_2.eps")
  MG1(EX=1.0/1.5, tmax=50)
  dev.copy2eps(file="queue_MG1_3.eps")
}