Éval de perf, TD2: Visualisation de données, application à l’analyse d’une file M/GI/1 (RICM4, 2014)

L’objectif est d’apprendre à visualiser des données de simulation en R, et en particulier la cohérence d’une simulation. Pour celà, nous utiliserons une simulation équationnelle d’une fille d’attente M/GI/1 avec un ordonnancement FIFO. Le processus d’arrivée est un processus de Poisson de taux (débit), les clients ont un temps de service de moyenne 1 pris comme unité de temps de référence.

Pour faire nos visualisations, nous utiliserons ggplot2. La doc regorge d’exemple et celle du R graphics cookbook également. En cas de soucis, google et stackoverflow sont vos amis. Enfin, les données n’étant pas toujours structurées comme on le souhaiterait, la bibliothèque reshape pourra s’avérer utile…

Dans la première section, je vous fourni le code effectuant la simulation (le code est un peu long car j’ai rajouté plein d’informations que vous pourrez visualiser) et dans le seconde, je visualise une simulation. Votre objectif est de refaire ces mêmes visualisations avec ggplot2.

Simulation de la file

Pour simuler une file FIFO, il suffit d’utiliser la méthode équationnelle à partir de la relation de récurrence sur les temps de réponse (équation de /Lindley/)

library(ggplot2)
library(reshape)
library(gridExtra)

## Loading required package: grid

set.seed(42);

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)
         )
}

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;
      }
    }

    list(jobs = df, events = data.frame(date=date,count=count),
         workevolution = data.frame(date=datework,work=work),
         sys_vide = sys_vide)
}    

Analyse d’une simulation.

Commençons par générer une trajectoire.

df <- Response(n=50, lambda=.7,"exp",1,0)
df$jobs$id=1:length(df$jobs$arrival)

Maintenant, regardons comment évoluent les départs et les arrivées.

ggplot(data=df$jobs) + 
  geom_point(aes(x=id,y=arrival)) + 
  geom_point(aes(x=id,y=completion),color="red") +
  xlab("Job Id") +
  ylab("Date") + 
  ggtitle("Arrival and departure of jobs in a M/M/1 queue")

Par contre, pas facile de mettre une légende dans ces conditions. Vous aurez donc intérêt à reformater vos données comme ceci:

df$jobs$response_time = df$jobs$completion-df$jobs$arrival
d <- melt(df$jobs,id=c("id"))

Cela vous permettra de facilement faire ceci:

ggplot(data=d[d$variable %in% c("arrival","completion"),],
       aes(x=id,y=value, color=variable)) + 
  geom_point() + scale_color_brewer(palette="Set1") +
  xlab("Job Id") +
  ylab("Date") + 
  ggtitle("Arrival and departure of jobs in a M/M/1 queue")

ou encore ceci qui est peut-être plus lisible ? Attention, ce qui est non-trivial ici, c’est que je ne stacke pas service et response_time mais que je les supperpose (avec même un peu de transparence pour détecter le cas aberrant où service serait plus grand que response_time ) de façon à percevoir le temps d’attente…

ggplot(data=d[d$variable %in% c("service","response_time"),],
       aes(x=factor(id),y=value, fill=variable)) + 
  geom_bar(data=d[d$variable %in% c("response_time"),],stat="identity", position="identity") +
  geom_bar(data=d[d$variable %in% c("service"),],stat="identity", position="identity", alpha=.8)+
  xlab("Job Id") +
  ylab("Duration") + 
  ggtitle("Response and service time of jobs in a M/M/1 queue")

Mais mettre le temps en ordonné, c’est un peu étrange en fait. Si on souhaite plutôt voir l’évolution au cours du temps, il faut le mettre en abcisse. Essayons plutôt ceci:

p1 <- ggplot(data=df$jobs) + geom_errorbarh(aes(y=id,xmin=arrival,xmax=completion,x=arrival)) + 
  xlab("Date") +
  ylab("Job Id") + 
  ggtitle("Job lifespan in a M/M/1 queue")
p1

Pas très pratique ce job id. Et si on regardait la taille de la file.

p2 <- ggplot(data=df$events) + geom_step(aes(x=date,y=count)) + 
  xlab("Date") +
  ylab("Queue length") + 
  ggtitle("Queue length in a FIFO M/M/1 queue")
p2

Ou encore, la quantité de travail au cours du temps

p3 <- ggplot(data=df$workevolution) + geom_line(aes(x=date,y=work)) + 
  xlab("Date") +
  ylab("Work") + 
  ggtitle("Work evolution in a FIFO M/M/1 queue")
p3

On peut alors même combiner les trois représentations.

grid.arrange(p1,p2,p3) # heights=c(1/2, 1/4, 1/4)