Q1

set.seed(10)
library(ggplot2)
library(reshape)
  # I will only deal with SRPT_pmtn and SPT_pmtn
FileSRPT <- function(n=1000,lambda=0.5,policy="SRPT_pmtn") {
  # simulates a M/GI/1 srpt queue with different preemption policy
  # parameters:
  #    n :  total number of jobs
  #    lambda : arrival rate
  #    policy: SRPT_pmtn, SPT_pmtn, SPT
  # return value:
  #    vector with response time of each task assuming the queue is initially empty
  
  A <- rexp(n,lambda)         # inter arrival
  t1 <- cumsum(A)             # arrival dates
  t2 <- rep(NA,n)  # completion dates
  S<-rexp(n)
  
  #### Variables that define the state of the queue
  t = 0               # time node of each event ocurring
  remaining = rep(NA,n)  # how much work remains to do for each task
  running = NA        # index of the currently running task
  waiting = c()       # stack with tasks which have arrived and have not been completed yet
  next_arrival = 1    # index of the next task to arrive
  
  get_task=function(temp=rep(0,n)){#choose a the right task
    running<<-switch(policy,
                     SRPT_pmtn=waiting[which.min(remaining)], #which.min() to get the index of mininum element
                     SPT_pmtn=waiting[which.min(temp)]
    )
  }
  
  #### A few useful local functions 
  run_task = function() { # runs the last task of the waiting list
    if(length(waiting)>0) {
      temp <- remaining
      if(policy == "SRPT_pmtn" & !is.na(running)){ 
        temp[running]<-S[running] 
        temp[waiting(length(waiting))] <- S[waiting(length(waiting))]
      }
      if(is.na(running))
        running<<-waiting[length(waiting)]
      else
        running <<- get_task(temp)
          
      remaining[running]<<-min(S[running],remaining[running],na.rm=T) 
      waiting <<- waiting[-length(waiting)]        
      
    }
  }
  
  #no change
  push_task = function() { 
        if(!is.na(running)) {
        waiting <<- c(waiting,running)
        }
        running <<- NA      
        waiting <<- c(waiting,next_arrival)
        next_arrival <<- next_arrival+1 
        if(is.na(running)) { run_task() }
    
  }
  
  #### Main simulation loop
  #no change
  while(TRUE) { 
    # Look for next event
    dt = NA
    if(next_arrival <=n) { dt = min(dt,(t1[next_arrival]-t), na.rm=T) }
    if(!is.na(running))  { dt = min(dt,remaining[running], na.rm=T)   }
    if(is.na(dt)) { break }
    
    # Update state
    t=t+dt 
    if(!is.na(running)) {
        remaining[running] = remaining[running] - dt
        if(remaining[running]<=0) { 
          t2[running]<- t 
            running = NA
            run_task()  
        }      
    }
    if((next_arrival<=n) & (t==t1[next_arrival])) { 
          push_task()
    }
  }
  
  t2-t1
} 

Q2

Service <- function(n=1,typeservice,x,y) {
# genere un temps de service
  switch(typeservice,
         det = seq(from=1,to=1,length.out=n),
         uni = runif(n,x,y),
         gamma = rgamma(n,shape=x,scale=y),
         exp = rexp(n,x)
         )
}

#FIFO
Response <- function(n,lambda,typeservice,x,y) {
    # genere une trajectoire de temps de reponse pour une file M/GI/1 FIFO
    # parametres :
    #    n :  taille de la trajectoire
    #    lambda : debit d'arrivee
    #    typeservice : type de service (det uni gamma exp)
    #    x ,y : parametres de la loi de service
    # valeur de retour :
    #    vecteur de n valeurs successives du temps de
    #    réponse dans une M/GI/1 FIFO initialement vide
    
    # generation des arrivees et des temps de service
    A <- rexp(n,lambda)
    t1 <- cumsum(A)
    S <- Service(n,typeservice,x,y)
    
    # initialisation des structures (pour éviter une allocation dynamique et gagner du temps)
    t2 <- seq(0,0,length.out=n)

    # Par curiosité comptons également le nombre de fois où le système se vide
    sys_vide = 1
    
    # initialisation du premier terme
    t2[1] <- t1[1]+ S[1]
    
    # calcul de la trajectoire equation de recurrence de Lindley
    for (i in 2:n) {
        t2[i] <- S[i] + max(t1[i], t2[i-1])
        if(t1[i]>t2[i-1]) {
          sys_vide = sys_vide + 1
        }
    }
        
    df=data.frame(arrival=t1,completion=t2,service=S)

    # Une itération sur ce que l'on vient de calculer pour calculer l'évolution de la
    # taille de la file d'attente au cours du temps.
    t1[n+1] = t2[n]+1;       # Un hack pour sortir de ma boucle
    date = rep.int(0,2*n);     # À initialiser absolument (performance + nom déjà utilisé) !
    count = rep.int(0,2*n);   
    i1 = i2 = 1;
    while(i1<=n | i2 <=n) {
      if(t1[i1]<t2[i2]) {
        date[i1+i2] <- t1[i1];
        count[i1+i2]<- count[i1+i2-1]+1;
        i1 <- i1+1;
      } else {
        date[i1+i2] <- t2[i2];
        count[i1+i2]<- count[i1+i2-1]-1;        
        i2 <- i2+1;
      }
    }
            
    # Une deuxième itération pour avoir l'évolution du travail au cours du temps.
    datework = rep.int(0,2*n);
    work =rep.int(0,2*n);
    i1 = i2 = 1;
    while(i1<=n | i2 <=n) {
      if(t1[i1]<t2[i2]) {
        i = 2*(i1+i2-1);
        datework[i] <- t1[i1];
        if(work[i-1]>0) {
           work[i] <- work[i-1] - (datework[i]-datework[i-1])
        } else {
           work[i] <- 0
        }
        if(work[i]<1e-15) { work[i]<-0 } # Fix floating point precision issue :( 
        datework[i+1] <- datework[i];
        work[i+1] <- work[i] + S[i1];
        
        i1 <- i1+1;
      } else {
        i = 2*(i1+i2-1);
        datework[i] <- t2[i2];
        work[i] <- work[i-1] - (datework[i]-datework[i-1]);
        datework[i+1] <- datework[i];
        work[i+1] <- work[i];

        i2 <- i2+1;
      }
    }

    df$completion - df$arrival
}    

df_fifo <- Response(n=1000, lambda=0.5,"exp",1,0)

df_srpt<-FileSRPT(n=1000,policy="SRPT_pmtn")

df_spt<-FileSRPT(n=1000,policy="SPT_pmtn")



dfr<-data.frame(SRPT_pmtn=df_srpt,SPT_pmtn=df_spt, FIFO=df_fifo) 
dfr$id=1:length(dfr$SRPT_pmtn)
d <- melt(dfr,id=c("id"))

ggplot(data=d[d$variable %in% c("SRPT_pmtn","SPT_pmtn", "FIFO"),],
       aes(x=id,y=value, color=variable)) + 
  geom_line(size=1) + scale_color_brewer(palette="Set1") +  xlab("Job Id") +
  ylab("Response Time") + 
  ggtitle("Response time")

sprintf("Average response time for SRPT_pmtn, SPT_pmtn and FIFO are %s, %s, %s",mean(df_srpt), mean(df_spt), mean(df_fifo) )

## [1] "Average response time for SRPT_pmtn, SPT_pmtn and FIFO are 1.75063251254203, 2.0920669595711, 1.84957006630721"

Q3

sprintf("Max response time for SRPT_pmtn, SPT_pmtn and FIFO are %s, %s, %s",max(df_srpt), max(df_spt), max(df_fifo) )

## [1] "Max response time for SRPT_pmtn, SPT_pmtn and FIFO are 28.2120779462152, 37.2354891224178, 14.3260307890545"