# Stochasticke simulacne metody
# Radoslav Harman, KAMS FMFI UK
# Cvicenie 5 (Simulovanie systemu hromadnej obsluhy) - strucny zaznam prikazov pre R


sho <- function(delta, Delta, mu, H, graph=TRUE) {
# Simulovanie jedneho dna specialneho systemu hromadnej obsluhy
# s jednou linkou, poissonovskym tokom netrpezlivych zakaznikov
# a s konstantnou dobou obsluhy. (Podrobnejsie na prednaske.)

  no <- 0     # celkovy pocet obsluzenych zakaznikov
  t <- 0      # cas prichodu dosial posledneho zakaznika do systemu
  l <- 0      # cas dokoncenia obsluhy dosial posledneho zakaznika, ktory sa zaradil do fronty
  tvec <- c() # doby prichodov jednotlivych zakaznikov (len kvoli grafu; takato forma tvorenia vektora je neefektivna)
  nvec <- c() # pocty zakaznikov vo fronte v casoch prichodu (len kvoli grafu)
  rvec <- c() # vektor urcujuci rozhodnutie zakaznika zoradit sa do fronty (len kvoli grafu)

  while(TRUE) {

    t <- t + rexp(1, rate=1/mu)
    if(t > Delta) break  
    if(graph) tvec <- c(tvec, t)

    # Vypocitajme aktualny pocet zakaznikov vo fronte
    n <- 0; if(l > t) n <- ceiling((l - t) / delta)
    if(graph) nvec <- c(nvec, n)

    if(runif(1) < exp(-H * n)){

      # Zakaznik sa rozhodol zaradit do fronty
      if(graph) rvec <- c(rvec, TRUE)
      no <- no + 1    
      if(l > t) {
        l <- l + delta
      } else {
        l <- t + delta
      }
    } else {
      if(graph) rvec <- c(rvec, FALSE)
    }
  }

  if(graph){
    plot(tvec[rvec], nvec[rvec], pch=19, cex=0.6, xlab="time", ylab="size of the queue")
    points(tvec[!rvec], nvec[!rvec], pch=19, col="red", cex=0.6)
  }

  return(no)  
}

sho(1, 600, 1, log(2)/10)  # "vyvazena" doba obsluhy a tok zakaznikov
sho(2, 600, 1, log(2)/10)  # prilis pomala obsluha
sho(1, 600, 2, log(2)/10)  # malo zakaznikov
sho(1, 600, 1, log(2)/3)   # velmi netrpezlivi zakaznici
sho(2, 600, 1, log(2)/50)  # pomala obsluha, ale trpezlivi zakaznici


sim.sho <- function(n, delta, Delta, mu, H, graph=TRUE) {
# Simulovanie specialneho systemu hromadnej obsluhy
# so sekvencnym pocitanim intervalu spolahlivosti
# n ... pocet simulacnych behov (Podrobnejsie na prednaske.)

  L <- 0; Q <- 0
  Z <- rep(0, n)           # Sluzi len na to, aby sme mohli zobrazit histogram
  ISd <- ISu <- rep(0,n)   # Sluzi na vykreslenie priebehu dolnych a hornych hranic IS
  kvantil <- qnorm(0.975)

  for(it in 1:n) {
    print(it)
    z <- sho(delta, Delta, mu, H, graph); Z[it]<-z
    L <- L + z; Q <- Q + z^2; priemer <- L / it
    poldlzka <- kvantil * sqrt((Q / it - priemer^2) / it)
    ISd[it] <- priemer - poldlzka
    ISu[it] <- priemer + poldlzka
  }

  hist(Z); windows()
  plot(ISd, ylim=c(min(ISd), max(ISu)), type="l", xlab="simulation run", ylab="IS")
  lines(ISu); lines((ISd + ISu) / 2, col="red")

  return(list(Z=Z, estimate=priemer, lower=ISd[n], upper=ISu[n]))
}

sim.sho(200, 1, 600, 1, log(2)/10)
sim.sho(200, 2, 600, 1, log(2)/10)
sim.sho(200, 1, 600, 2, log(2)/10)
sim.sho(200, 1, 600, 1, log(2)/3) 
sim.sho(200, 2, 600, 1, log(2)/50)
sim.sho(2000, 2, 600, 1, log(2)/50, FALSE)