RICM4: Évaluation de Performance

Table of Contents

Sitemap

---> misc
| ---> 2016
| ---> 2015
| ---> 2014
| ---> 2013
| ---> 2012
`--> Agenda

Informations Générales

Florence Perronnin est chargée des cours et Arnaud Legrand s'occupe des TDs.

Le planning avec les salles de cours est disponible ici.

La page de l'an dernier est ici.

Programme du cours

Semaine Cours (Jeudi, 8h00) TD (Vendredi, 11:30 ou exceptionnellement 8:00)
12-13 janvier Slides d'introduction Mesures de performance de code
19-20 janvier Introduction aux chaînes de Markov en temps discret Suite du TD précédent. Rendu DM pour le 3 février.
26-27 janvier Voir slides précédents Pas de TD.
2-3 février Fin des chaînes de Markov en temps discret. Illustration par simulation Modélisation et Simulation d'un Cache Web. Voir le corrigé détaillé de l'an dernier.
9-10 février Échange de cours. Rattrapage plus tard Page rank et marches aléatoires
16-17 février Loi Exponentielle + Poisson Fin du cours + Aloha
2-3 mars Chaine de Markov à temps continu, exemple de la M/M/1 Pas de cours (Projet)
9-10 mars Modélisation Markovienne du protocole Aloha Simulation à évènement discret d'une M/M/1 (version de base)
16-17 mars Files d'attentes Simulation à évènement discret d'une M/M/1 (exploitation)
22-23 mars Quick Simulation à évènement discret d'une M/M/1 (régression linéaire)
28-31 mars Cours Régression Linéaire Simulation à évènement discret d'une M/M/1 (réduction de variance)
7 avril   Révisions

TDs

Mesures de performance

On souhaite évaluer l'opportunité de remplacer l'implémentation du quicksort de la libc par une version parallèle afin que toutes les applications puissent bénéficier automatiquement des architectures multi-coeurs. L'objectif de ces séances est de répondre aux questions suivantes:

  • Quels sont les paramètres sur lesquels on peut jouer simplement pour améliorer les performances de ces différentes versions ?
  • À partir de quelle taille de tableau devient-il opportun d'utiliser la version parallèle ?

Votre mission pour la prochaine fois:

  1. Forkez le repos suivant https://github.com/alegrand/M2R-ParallelQuicksort puis faites un checkout sur votre machine. Ce repository comporte un code C qui réalise un quicksort séquentiel ou multi-threadé.
  2. Reproduire mes résultats sur votre machine. Un conseil, prenez un environnement unix (linux ou MacOsX) car sous windows, vous aurez des difficultés (ceci dit, l'an dernier, certains l'avaient fait et s'en étaient sortis) du genre: problème de noms de fichiers (le : est interdit sous windows!), de compilateurs (et oui, pas de compilateur C d'installé par défaut sous windows), de scripts (pas de bash ni de perl). Vous pouvez utiliser une VM mais dans ce cas là, vous aurez peut-être du mal à voir l'intérêt du parallélisme.
  3. Améliorer ce que vous voulez en prenant des notes sur ce que vous faites! Comme nous l'avons vu, il y a plein de points d'améliorations possibles:

    • le plan d'expérience
    • la façon dont les mesures sont faites (contrôle des paramètres, choix d'une configuration "intéressante")
    • l'analyse

    Donc n'hésitez pas à vous répartir le travail.

  4. Se préparer à présenter aux autres.

Quelques liens utiles:

Modélisation et Simulation d'un Cache Web

On dispose de \(N=100\) objets et on va étudier les performances d'un cache LRU de taille \(K\). Pour simplifier, on va supposer qu'un des objets est plus populaire (probabilité \(a\)) que les autres (probabilité \(b\)).

N = 100
a = 0.1
b = (1-a)/(N-1)
K = 20
T = 1000
  1. Générer une séquence d'accès et vérifier que la séquence respecte bien les propriétés attendues.

    Facile, avec la fonction sample!

    access = sample(c(1:N), T, prob=c(a,rep(b,N-1)), replace=TRUE)
    head(access)
    
    [1] 39 60 97 64  1 33
    

    On vérifie quand même que ça ressemble à ce à quoi on s'attendait!

    plot.ts(access)
    

    figure9036McI.png

    hist(access,breaks=N)
    

    figure9036ZmO.png

  2. Simuler la politique LRU en partant d'un cache vide. Quelques astuces de mise en œuvre:

    • 5 %in% X vous indique si la valeur 5 apparaît dans le vecteur X
    • (1:N)[X==5] vous indique la (les) position de la valeur 5 dans le vecteur X de longueur N.
    • X[-4] est le vecteur X privé de son 4ème élément.

    Allons-y, il suffit de prendre les choses par le bon bout.

    cache = sample(c(1:N), K, prob=c(a,rep(b,N-1)), replace=FALSE) # Un cache initial 
    miss=c()
    for(obj in access) {
        if(obj %in% cache) {
            pos = (1:K)[cache==obj];
            cache = c(obj,cache[-pos])
            miss=c(miss,0);
        } else {
            cache = c(obj,cache[-K])
            miss=c(miss,1);
        }
    }
    
    head(posA)
    head(miss)
    head(missA)
    
    [1] 16 18 20 36 38 40
    [1] 1 1 1 1 0 0
    [1] 0 0 0 0 0 0
    
    plot.ts(head(miss,n=40))
    

    figure9036m3I.png

  3. Améliorer votre code en stockant quand il y a un défaut de cache, quand ça concerne l'objet A et la position de l'objet A (pour débuger votre code…).

    Le code va vite devenir moche mais n'est pas compliqué pour autant.

    cache = sample(c(1:N), K, prob=c(a,rep(b,N-1)), replace=FALSE) # Un cache initial 
    miss=c()
    missA=c()
    posA=c();
    for(obj in access) {
        if(obj %in% cache) {
            pos = (1:K)[cache==obj];
            cache = c(obj,cache[-pos])
            miss=c(miss,0);
            missA=c(missA,0);
        } else {
            cache = c(obj,cache[-K])
            miss=c(miss,1);
            if(obj == 1) {
               missA=c(missA,1);
            } else {
               missA=c(missA,0);
            }
        }
        if(1 %in% cache) {
            posA = c(posA, (1:K)[cache==1]);
        } else {
            posA = c(posA, NA);
        }
    }
    head(posA)
    head(miss)
    head(missA)
    
    [1] NA NA NA NA  1  2
    [1] 1 0 1 1 1 1
    [1] 0 0 0 0 1 0
    

    Est-ce que la position de A est cohérente ? Oui, il monte, stagne de temps en temps, et repasse brutalement en position 1.

    plot.ts(head(posA,40))
    

    figure9036NWb.png

  4. Estimer la probabilité de défaut de cache et l'intervalle de confiance. Votre calcul est-il légitime ?

    #   mean(miss)
    err = sd(miss)/sqrt(length(miss))
    c(mean(miss)-2*err, mean(miss) +2*err)
    
    [1] 0.731937 0.786063
    

    Bof, pas terrible comme estimateur.

  5. Modéliser ce système par une chaîne de Markov.

    On a fait ça ensemble au tableau.

  6. Construire la matrice correspondante et l'utiliser pour calculer le régime stationnaire.

    Voir le corrigé de l'an dernier.

  7. En déduire la probabilité que A soit dans le cache. Comparer la valeur exacte avec l'estimation obtenue en simulation à évènements discrets.
  8. En déduire la probabilité qu'il y ait un défaut de cache. Comparer la valeur exacte avec l'estimation obtenue en simulation à évènements discrets.
  9. Tracer la probabilité de défaut de cache en fonction de la taille du cache.

Page rank et marches aléatoires

  • Google: notion de popularité, pas basée uniquement sur le contenu qui est facilement "forgeable" mais plutôt sur les interconnexions en tre les pages.
  • Intuitivement:

    1. plus on pointe vers la page i, plus elle est populaire. Donc la popularité est proportionnelle au degré entrant.
    2. plus le degré sortant d'une page est élevé, plus sa contribution à la popularité des pages devrait être faible
    3. Plus une page est populaire, plus sa contribution à la popularité des autres pages est élevée (notion "d'expert").

    Au final il semble y avoir un point compliqué, la définition de la popularité est récursive à cause du 3ème point. Une équation possible serait:

    \(\pi_i = \sum_{j} 1/|\sum_k L_{j,k}| \pi_j\) si \(L\) est la matrice des connexions avec des 0 et des 1

  • La popularité d'une page, c'est défini par la fréquence à la quelle les gens cliquent vers cette page, ce que Google va avoir du mal à tracer complètement… Et intuitivement, les surfeurs, se baladent sur le graphe du web. Quels sont les problèmes potentiels

    • Revenir plusieurs fois sur la même page. Dans le cadre d'une exploration aléatoire, ça n'est pas forcément un problème.
    • D'où part on ? Est-ce grave ? Il pourrait y avoir plusieurs composantes connexes et notre surfer serait bloqué juste dans une. S'il n'y a que quelques composantes, pourquoi pas, notre marche aléatoire pourrait nous donner des stats intéressantes mais si c'est plus fragmenté, ça va être compliqué.
    • Si on a une page sans lien vers d'autres pages, c'est un puits et tout surfeur va fatalement y tomber et y rester coincé.

    Il faudrait considérer que le graphe est fortement connexe mais ce n'est pas le cas. Une solution envisagée par les fondateurs de Google, était de se "téléporter". En moyenne au bout de 5 clicks, on repart au hasard de n'importe où. Ça revient à rajouter des probas très faibles vers tout le monde.

Au final, tout ça revient au même. Calculer des statistiques sur une marche aléatoire dans le graphe ou bien calculer la distributation stationnaire de la chaine de Markov associée (c'est à dire les valeurs/vecteurs propres de la distribution associée), c'est pareil.

Après discussions, on identifie 3 méthodes pour calculer le régime stationnaire:

  1. Faire une marche aléatoire sur le graphe associé et on calcule des statistiques à la fin sur la fréquence de passage dans chaque page. Intuitivement, la complexité est très très faible mais on a une estimation d'assez mauvaise qualité.
  2. Calculer les valeurs propres de la matrice associée (comme on a fait plus haut pour le cache). C'est très coûteux (\(n^3\)) pour un graphe de la taille du web mais on a un calcul exact
  3. Itérer la fonction! C'est à dire partir d'un vecteur \(\pi=(1/N,1/N,\dots,1/N)\).

    Que se passe-t-il sur la matrice ? Si \(A= MDM^{-1}\) alors \(A^k = MD^kM^{-1}\). Vu que notre matrice a une seule valeur propre à 1 et que les autres sont plus petites que 1, elles vont disparaître lors de l'itération et il ne restera plus que celle associée au vecteur propre 1.

Enfin, deux documents que vous pouvez vouloir regarder et dont je me suis inspiré pour cette séance: cette série de slides (lien d'origine) et cet article publié dans CACM en 2011 mais dont voici la version mise à jour et libre d'accès sur arxiv.

Simulation à évènement discrets de file d'attentes

L'objectif du TD est de comprendre la bonne façon d'organiser une simulation à évènement discret (sans reposer sur les hypothèse Markoviennes) avec la notion d'état du système, d'évènements possible et de mise à jour de l'état du système. Si vous lisiez des bouquins qui expliquent comment faire de la simulation à évènement discret, il n'est pas clair que ça vous aiderait vraiment. Le mieux est d'en écrire un ensemble pour que vous compreniez comment procéder. Pour celà, on va s'intéresser au cas simple d'une file d'attente.

Modélisation et principe de simulation

  • Paramètres du système
    • Taux d'arrivée dans le système lambda (en tâches par seconde)
    • Taux de service du système mu (en tâches par seconde)
    • Politique de service (on va coder FIFO là mais on essaye de coder tout ça de la façon la plus générique possible)
  • Description des variables d'état importantes:
    • Dates d'arrivées Arrival
    • Temps de service Service
    • Temps résiduel Remaining
    • Date de terminaison Completion
    • Date courante t
    • Client actuellement servi CurrClient (utile pour carractériser la politique de service utilisée)
  • Évolution possible
    • Soit une arrivée
    • Soit la terminaison d'une tâche
  • Structure du code:

    while(T) {
        dtA = ...  # temps jusqu'à la prochaine arrivée
        dtC = ...  # temps jusqu'à la prochaine terminaison
        if(is.na(dtA) & is.na(dtC)) {break;}
        dt = min(dtA,dtC)
        # Mettre à jour comme il faut.
    }
    

Pour la prochaine fois:

  • Repartir de votre propre version (si vous voulez bien comprendre comment ça marche) ou à défaut de celle donnée ci dessous (si vous renoncez à écrire ce code vous-même proprement…) et l'encapsuler de façon à pouvoir facilement controler les paramètres du système (lambda, mu)
  • Fixer le taux de service mu à 1 et étudier (je veux une courbe!) le temps de réponse en fonction du taux d'arrivée lambda. Idéalement, vous devriez obtenir quelque chose du genre:

    MM1.png

Principe: Une première version faite en TD

set.seed(37);
N = 10; # N <- 100; # N <<- 100;
lambda = .1;
mu = 1;
Arrival = cumsum(rexp(n=N, rate = lambda));
Service = rep(N,x=1/mu); # rexp(n = N, mu)
Remaining = rep(N, x=NA);
Completion = rep(N, x=NA);
t = 0;
CurrentTask = NA;
NextArrival = 1;
while (TRUE) {
    print(t);
    print(Arrival);
    print(Service);
    print(Remaining);
    dtA = NA;
    dtC = NA;
    if(length(Arrival[Arrival>t])>0) {
        dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
    }
    if(!is.na(CurrentTask)) {
        dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
    }
    if(is.na(dtA) & is.na(dtC)) {
        break;
    } 
    dt = min(dtA,dtC,na.rm=T)

    # Mettre à jour comme il faut.
    t = t + dt;
    if((NextArrival <=N) & (Arrival[NextArrival] == t)) {
        Remaining[NextArrival] = Service[NextArrival];
        NextArrival = NextArrival + 1;
    }
    if(!is.na(CurrentTask)) {
        Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
        if(Remaining[CurrentTask] <= 0) {
            Completion[CurrentTask] = t;
            Remaining[CurrentTask] = NA;
        }
        CurrentTask = NA;
    }
    WaitingList=(1:N)[!is.na(Remaining)];
    if(length(WaitingList)>0) {
        CurrentTask = head(WaitingList,n=1);
    }
}
print(Completion)
[1] 0
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 0.9928169
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1]  1 NA NA NA NA NA NA NA NA NA
[1] 1.992817
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 24.40118
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA  1 NA NA NA NA NA NA NA NA
[1] 25.40118
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 27.37712
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA  1 NA NA NA NA NA NA NA
[1] 28.37712
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 36.97357
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA  1 NA NA NA NA NA NA
[1] 37.97357
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 47.98916
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA  1 NA NA NA NA NA
[1] 48.98916
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 68.15102
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA  1 NA NA NA NA
[1] 69.15102
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 83.79521
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA  1 NA NA NA
[1] 84.79521
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 95.51369
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA  1 NA NA
[1] 96.51369
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 100.0268
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA  1 NA
[1] 101.0268
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
[1] 106.9737
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA  1
[1] 107.9737
 [1]   0.9928169  24.4011768  27.3771164  36.9735672  47.9891610  68.1510163
 [7]  83.7952085  95.5136902 100.0267969 106.9736605
 [1] 1 1 1 1 1 1 1 1 1 1
 [1] NA NA NA NA NA NA NA NA NA NA
 [1]   1.992817  25.401177  28.377116  37.973567  48.989161  69.151016
 [7]  84.795209  96.513690 101.026797 107.973660

Suite faite en TD: utilisation du code précédent pour étudier le temps de réponse d'une file M/G/1

set.seed(37);

MM1 = function(N = 100, lambda = .1, mu = 1, service="det", debug=F) {
    Inter = rexp(n=N, rate = lambda)

    Arrival = cumsum(Inter);
    if(service == "det") {
       Service = rep(N,x=1/mu); # rexp(n = N, mu)
    } else if (service == "exp") {
       Service = rexp(n = N, mu)
    } else {
       error("WTF!")
    } 
    Remaining = rep(N, x=NA);
    Completion = rep(N, x=NA);
    t = 0;
    CurrentTask = NA;
    NextArrival = 1;
    while (TRUE) {
        if(debug) {
            print(t);
            print(Arrival);
            print(Service);
            print(Remaining);
        }
        dtA = NA;
        dtC = NA;
        if(length(Arrival[Arrival>t])>0) {
            dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
        }
        if(!is.na(CurrentTask)) {
            dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
        }
        if(is.na(dtA) & is.na(dtC)) {
            break;
        } 
        dt = min(dtA,dtC,na.rm=T)
        
    # Mettre à jour comme il faut.
    t = t + dt;
        if((NextArrival <=N) & (Arrival[NextArrival] == t)) {
            Remaining[NextArrival] = Service[NextArrival];
            NextArrival = NextArrival + 1;
        }
        if(!is.na(CurrentTask)) {
            Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
            if(Remaining[CurrentTask] <= 0) {
                Completion[CurrentTask] = t;
                Remaining[CurrentTask] = NA;
            }
            CurrentTask = NA;
        }
        WaitingList=(1:N)[!is.na(Remaining)];
        if(length(WaitingList)>0) {
            CurrentTask = head(WaitingList,n=1);
        }
    }
    na.fail(Completion);
    if(debug) {print(Completion) ;}
    return(data.frame(N = N, lambda=lambda, mu = mu, service = service,
                      resp = mean(Completion - Arrival), 
                      resp_err = sd(Completion - Arrival)/sqrt(N)))
 }
MM1(N=100,lambda=.4,service="exp")
    N lambda mu service     resp  resp_err
1 100    0.4  1     exp 1.626987 0.1434172
df = data.frame();
for(service in c("det","exp")) {
    for(lambda in seq(from=.1,to=.9,by=.1)) {
        df = rbind(df,MM1(N=ceiling(40/(1-lambda)**2),lambda=lambda,service=service));
         # ceiling(40/(1-lambda)**2)   N = 1000
    }
}
df;
   lambda service     resp   resp_err    N
1     0.1     det 1.024441 0.01548081   50
2     0.2     det 1.039429 0.01596869   63
3     0.3     det 1.242695 0.05018003   82
4     0.4     det 1.316764 0.04892153  112
5     0.5     det 1.494866 0.06623087  160
6     0.6     det 1.812693 0.07703218  250
7     0.7     det 1.773051 0.04136567  445
8     0.8     det 2.396349 0.04671912 1001
9     0.9     det 5.554034 0.07156642 4001
10    0.1     exp 1.196517 0.18172348   50
11    0.2     exp 1.347830 0.17005655   63
12    0.3     exp 1.470158 0.14457873   82
13    0.4     exp 1.889351 0.17304625  112
14    0.5     exp 1.969257 0.15328482  160
15    0.6     exp 1.993911 0.11877594  250
16    0.7     exp 3.236034 0.12677742  445
17    0.8     exp 4.476521 0.15028374 1001
18    0.9     exp 8.750325 0.11068925 4001
library(ggplot2);
myrate = function(x) {1/(1-x)} # la formule théorique 1/(mu-lambda) pour une vraie MM1
ggplot(data=df, aes(x=lambda, y=resp, color=service)) + geom_point() + 
   geom_errorbar(aes(ymin=resp-2*resp_err, ymax=resp+2*resp_err),width=.05) + 
#    geom_smooth(method="lm",formula=y~I(1/(1-x))) + # this is not the model you should expect
    xlim(0,1)  + ylim(0,15)  +   
    stat_function(fun = myrate, geom = "line",color="blue")

MM1.png

Plusieurs remarques:

  • Le fait d'avoir une loi de service exponentielle plutôt que déterministe joue à la fois sur la variabilité et sur la moyenne du temps de réponse.
    • Avoir un temps de service plus variable va forcément induire un temps de réponse plus variable puisque le temps de réponse, c'est le temps de service plus le temps d'attente.
    • Mais avoir un temps de service plus variable a crée une variabilité plus importante du temps d'attente et quand on se retrouve à attendre longtemps derrière quelqu'un, on est rarement seul dans la queue. Du coup, le temps d'attente augmente pour tout le monde, d'où un temps de réponse moyen bien plus élevé.
    • Note: la variabilité des interarrivée aussi impacte le temps de réponse. Si on se mettait dans le cas complètement déterministe, tant que λ<μ, on n'a pas d'attente du tout et le temps de réponse est constant égal au temps de traitement…
  • Plus λ est élevé plus la variabilité est grande et donc plus l'estimation qu'on fait est à prendre avec précaution (l'errorbar est plus large). C'est d'autant plus vrai que pour appliquer rigoureusement le théorème central limite il faut:
    • Connaître σ alors que là je ne fais qu'une estimation avec sd et plus elle est élevée plus cette estimation est incertaine.
    • Il faut que les \(N\) observations soient indépendantes, ce qui n'est pas le cas du tout ici. Si un client a un temps d'attente élevé, il est très probable que le client d'après aussi ait un temps d'attente élevée… Le seul cas où ils seraient indépendants serait si le système s'est vidé. On devrait donc compter le nombre de fois où le système se vide, calculer le temps de réponse moyen des clients entre deux telles situations et appliquer le théorème central limite à ces nouveaux temps de réponse moyen… Ces instants où le système se vide s'appellent instant de regénération. Les errorbar vont alors énormément s'élargir et ce d'autant plus quand le système est chargé car la regénération va arriver de moins en moins souvent. Donc plus le système est chargé, plus estimé le temps de réponse devient difficile…
  • L'hypothèse d'un taux d'arrivée exponentiel est naturelle mais celle d'un taux de service exponentiel est très très discutable. Ceci dit, le comportement est assez proche et pour peu que la variance du système soit plutôt plus faible que celle d'une loi exponentielle, on obtient alors un genre de pire cas. En fait les estimations qu'on obtient avec l'hypothèse Markov/loi exponentiels sont assez robustes au changement de loi même s'il n'y a pas de méta-théorème qui le garantisse. Pour peu qu'on sache l'utiliser, c'est donc un outil très puissant pour évaluer les performances d'un système.
  • Ne pas oublier que pour une trajectoire/simulation de file d'attente donnée, on a une phase de remplissage du système, un état stationnaire, puis une phase où le système se vide. Clairement identifier ces différentes phases est difficile mais si on étudie des configurations bizarres (N=10 ou bien λ>μ), on obtiendra des valeurs de temps de réponse moyen et on calculera des intervalles de confiance mais qui n'ont tout simplement aucun sens… On peut parfaitement ne plus avoir que la partie remplissage et vidage du système. Une bonne façon de vérifier que les choses ont du sens c'est au moins de vérifier si les statistiques qu'on obtient se stabilisent bien quand on augmente N.

Application du modèle linéaire à l'étude de files d'attentes

Quand on voit un asymptote comme ça, ça sent le "you shall not pass". ☺ Il serait donc certainement plus logique de regarder 1/resp.

library(ggplot2);
ggplot(data=df, aes(x=lambda, y=1/resp, color=service)) + geom_point() + 
    geom_smooth(method="lm") +
    xlim(0,1)

MM1_inv.png

Ah oui, effectivement, ça a l'air déjà bien plus simple (tout bêtement linéaire) comme ça. Bon, mais du coup, on devrait arriver à en déduire les paramètres correspondants.

Commençons par extraire la partie pour le temps de service déterministe.

dd = df[df$service=="det",]
dd
  lambda service     resp   resp_err    N
1    0.1     det 1.024441 0.01548081   50
2    0.2     det 1.039429 0.01596869   63
3    0.3     det 1.242695 0.05018003   82
4    0.4     det 1.316764 0.04892153  112
5    0.5     det 1.494866 0.06623087  160
6    0.6     det 1.812693 0.07703218  250
7    0.7     det 1.773051 0.04136567  445
8    0.8     det 2.396349 0.04671912 1001
9    0.9     det 5.554034 0.07156642 4001
plot(data=dd, 1/resp ~ lambda)

MM1_inv_scatter.png

C'est comme avant mais sans la droite de régression. Faisons le fit nous-même.

summary(lm(data=dd, 1/resp ~ lambda))
Call:
lm(formula = 1/resp ~ lambda, data = dd)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.10657 -0.03271  0.01383  0.03286  0.09378 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.11280    0.04458   24.96 4.22e-08 ***
lambda      -0.91797    0.07921  -11.59 8.03e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06136 on 7 degrees of freedom
Multiple R-squared:  0.9505,	Adjusted R-squared:  0.9434 
F-statistic: 134.3 on 1 and 7 DF,  p-value: 8.033e-06

Le R square est élevé, les estimations bonnes. Super. On aurait donc envie de dire que \(resp= 1/(1.11-0.91*\lambda)\) (à prendre avec des pincettes! je vous laisse réfléchir à la confiance qu'on peut avoir dans de telles estimations et comment on pourrait améliorer les choses).

Vérifions que tout va bien du côté des résidus

par(mfrow=c(2,2)); plot(lm(data=dd, 1/resp ~ lambda)) ; par(mfrow=c(1,1));

MM1_inv_4plots.png

Pour mémoire, la première chose que vous avez faite, c'était sans inverser resp:

plot(data=dd, resp ~ lambda)

MM1_scatter.png

C'était évidemment pas linéaire mais ça n'empèche pas de faire la régression (qui va raconter n'importe quoi mais c'est pas forcément évident au début de comprendre pourquoi).

summary(lm(data=dd, resp ~ lambda))
Call:
lm(formula = resp ~ lambda, data = dd)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.9801 -0.5447 -0.2491  0.2651  2.0094 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.01722    0.70722  -0.024   0.9813  
lambda       3.95763    1.25677   3.149   0.0162 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9735 on 7 degrees of freedom
Multiple R-squared:  0.5862,	Adjusted R-squared:  0.5271 
F-statistic: 9.917 on 1 and 7 DF,  p-value: 0.01617

Le R squared est mauvais, évidemment mais c'est sans importance en fait, cette valeur ne veut rien dire car le modèle est complètement faux. On voit bien que les estimations sont très mauvaises et douteuses (pas de *** comme avant…).

Et quand on regarde ce qu'il se passe du côté des résidus, on peut un peu comprendre:

par(mfrow=c(2,2)); plot(lm(data=dd, resp ~ lambda)) ; par(mfrow=c(1,1));

MM1_4plots.png

Seul le premier plot est important. On voit que ça n'est absolument pas aléatoire et qu'il y a une sacré structure complètement ignorée. Le modèle est à changer!

Modèle linéaire et réduction de variance

  • Modèle linéaire

    Quels sont les problèmes dans la régression précédente?

    1. On ne peut pas calculer des intervalles de confiance n'importe comment sur des variables corrélées. Or les temps de réponses sont corrélés. Le plus simple est de repérer quand le système se vide complètement et de calculer le temps de réponse moyen sur ce batch. Ensuite, on pourra faire des statistiques sur "le temps de réponse d'un batch".
    2. On fait des régressions sur des valeurs agrégées et pas sur les données brutes. En général c'est une mauvaise idée. Attention, voulez vous regarder l'inverse du temps de réponse moyen ou la moyenne de l'inverse de temps de réponse ? C'est pas la même chose. Il faut donc bien réfléchir à votre modèle et à ce que vous voulez calculer.
    3. Le système est hétéroscédastique (quand λ augmente, la variabilité du temps de réponse augmente beaucoup). Il faut peut-être faire la régression un poil différemment.
    4. Il faut faire plus de mesures là où il y a beaucoup de variabilité et moins là où c'est stable. Par exemple: N=ceiling(40/(1-lambda)**2)

    Essayer de corriger le code précédant pour prendre en compte ces remarques et améliorer les choses.

  • Réduction de variance à l'aide de variables de contrôle.

    Présentation du principe de réduction de variance par variable de contrôle. Voir vieux TD. L'idée c'est que la variable qui nous intéresse (typiquement le temps de réponse) peut être assez variable mais que sa variabilité peut souvent s'expliquer par d'autres variables qu'on connaît bien. Dans le cas de la M/G/1, on pourrait par exemple utiliser le temps de service et le temps d'inter-arrivée qui imposent une variabilité.

    Essayer d'utiliser cette technique pour améliorer et mieux contrôler la qualité de l'estimation.

  • Réduction de variance à l'aide de variables antithétiques

    Présentation du principe de réduction de variance à l'aide de variables antithétiques. Voir vieux TD. L'idée est de faire des simulations de durée N/2, en contrôlant bien la graine, et d'utiliser \(U\) et \(1-U\) dans le générateur. Ces deux types de trajectoires sont d'une certaine façon correllées et ont plus de chances de se "compenser".

Une version avec tout ce qu'on a vu

  • Le code
    1. Tout d'abord, corrigeons cette histoire d'évaluation de l'intervalle de confiance. Pour être tout à fait honnête, quand j'ai commencé à écrire le code et à réfléchir sur si ce que je faisais était correct, je me suis rendu compte que c'était bien plus complexe que ce que j'imaginais. En me référant aux sources sur ce sujet (merci Jean-Marc!), i.e., le papier de 1975 sur la simulation régénérative par Lavenberg et Slutz, on trouve les explications sur comment calculer correctement l'intervalle de confiance.
    2. Ensuite, il faut déterminiser la simulation (pour pouvoir utiliser les variables antithétiques). Ça m'a permis de trouver la source des bugs qui survenaient de temps en temps. Une égalité sur des flottants (en plus, je me souviens très bien quand on l'a écrit avoir dit "ah, ça c'est un peu dangereux mais bon, on va faire avec pour l'instant"). Au passage, je vous rappelle que, comme dans la majorité des langages:

      3-2.9 == .1
      
      [1] FALSE
      
    3. Je rajoute la possibilité de renvoyer l'ensembles des variables d'état de la simulation (pour tenter de réduire la variance via les variables antithétiques ou via les variables de contrôle)
    set.seed(37);
    
    MM1 = function(N = 100, lambda = .1, mu = 1, service="det", 
                   seed=NA, antithetic_lambda=F, antithetic_mu=F, 
                   debug=F, output_state = F) {
        start.time <- Sys.time()
    
        if(!is.na(seed)) {
            set.seed(seed);
        }
    
        # Génération des arrivées
        U = runif(N);
        if(antithetic_lambda) { U=1-U; }
        Inter = -log(U)/lambda; # Inter = rexp(n=N, rate = lambda)
        Arrival = cumsum(Inter);
    
        # Génération des temps de service
        U = runif(N);
        if(antithetic_mu) { U=1-U; }
        if(service == "det") {
           Service = rep(N,x=1/mu); # rexp(n = N, mu)
        } else if (service == "exp") {
           Service = -log(U)/mu; # Service = rexp(n = N, mu)
        } else {
           error("WTF!")
        }
    
        # Initialisation des variables d'état de la simulation
        Remaining = rep(N, x=NA);
        Completion = rep(N, x=NA);
        t = 0;
        CurrentTask = NA;
        NextArrival = 1;
        CurrentBatch = c();
        BatchResponse = c();
        BatchSize = c();
    
        # Boucle principale
        while (TRUE) {
            if(debug) {
                print(paste("###########",t,"###########"));
                print(c("Arrival:", Arrival));
                print(c("Service:", Service));
                print(c("Remaining:", Remaining));
                print(c("Completion:", Completion));
                print(c("NextArrival", NextArrival));
            }
            dtA = NA;
            dtC = NA;
            if(length(Arrival[Arrival>t])>0) {
                dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
            }
            if(!is.na(CurrentTask)) {
                dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
            }
            if(is.na(dtA) & is.na(dtC)) {
                break;
            } 
            dt = min(dtA,dtC,na.rm=T)
            
        # Mise à jour des variables d'état
            t = t + dt;
            if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { # Arrivée
                #### WARNING : Ce qui buggait avant, c'est qu'on testait 
                ####     Arrival[NextArrival] == t    :(
                Remaining[NextArrival] = Service[NextArrival];
                CurrentBatch = c(CurrentBatch, NextArrival);
                NextArrival = NextArrival + 1;
            }
            if(!is.na(CurrentTask)) { # Terminaison
                Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
                if(Remaining[CurrentTask] <= 0) {
                    Completion[CurrentTask] = t;
                    Remaining[CurrentTask] = NA;
                }
                CurrentTask = NA;
            }
            # Sélection de la prochaine tâche à exécuter (FIFO)
            WaitingList=(1:N)[!is.na(Remaining)];
            if(length(WaitingList)>0) {
                CurrentTask = head(WaitingList,n=1);
            }
            if(is.na(CurrentTask)) { # Plus de tâche en cours d'exécution
                { ### Tracking batches
                    BatchResponse = c(BatchResponse,mean(Completion[CurrentBatch]-Arrival[CurrentBatch]))
                    BatchSize = c(BatchSize,length(CurrentBatch))
                    if(debug) {
                        print(c("CurrentBatch:", CurrentBatch));
                        print(c("BatchResponse:", BatchResponse));
                    }
                    CurrentBatch = c();
                }
            }
        }
        # Vérification minimale de la cohérence du résultat
        na.fail(Completion);
    
        if(debug) {print(Completion) ;}
    
        Response = Completion - Arrival
        { ## Computing statistics per batch
            Num_Indep_Measurements = length(BatchResponse);
            Q = mean(Response); # or sum(sigma)/sum(nu)
            sigma = (BatchResponse*BatchSize);
            nu = BatchSize
            V1 = var(sigma);
            V2 = var(nu);
            V12 = cov(nu, sigma)
            V = V1 - 2*Q*V12 + Q*Q*V2;
            ## print(c("Q:", Q));
            ## print(c("V1:", V1));
            ## print(c("V2:", V2));
            ## print(c("V12:", V12));
            ## print(c("V:", V));
            err = sqrt(V/Num_Indep_Measurements)/sqrt(mean(nu));
        }
    
        end.time <- Sys.time()
        time.taken <- end.time - start.time
        time.taken
    
        if(!output_state) {
           return(data.frame(N = N, lambda=lambda, mu=mu, service = service, 
                             resp = mean(Response), 
                             resp_err = sd(Response) /sqrt(N),
                             resp_size = N,
                             resp_err_batch = err,
                             resp_size_batch = Num_Indep_Measurements,
                             time_taken = time.taken))
        } else {
            return(data.frame(arrival = Arrival,
                              service = Service,
                              response = Response,
                              earlier_arrival = Inter,
                              time_taken = rep(time.taken,N))) # YUCK!
        }
     }
    

    Vérifions qu'on arrive toujours à obtenir des choses du même genre qu'avant:

    MM1(N=100,lambda=.4,service="exp", seed=42)
    
        N lambda mu service     resp resp_err resp_size resp_err_batch resp_size_batch      time_taken
    1 100    0.4  1     exp 1.636837 0.181725       100      0.4923027              57 0.01629901 secs
    
  • Comparaison des intervalles de confiance
    df = data.frame();
    for(service in c("det","exp")) {
        for(lambda in seq(from=.1,to=.9,by=.1)) {
            df = rbind(df,MM1(N=ceiling(160/(1-lambda)**2), lambda=lambda, 
                              service=service, seed=42));
        }
    }
    df;
    
           N lambda mu service     resp   resp_err resp_size resp_err_batch resp_size_batch       time_taken
    1    198    0.1  1     det 1.062065 0.01468198       198     0.01810787             174  0.02598739 secs
    2    250    0.2  1     det 1.133452 0.02019338       250     0.02801578             197  0.02644062 secs
    3    327    0.3  1     det 1.210554 0.02293535       327     0.03752081             230  0.03600883 secs
    4    445    0.4  1     det 1.300432 0.02404866       445     0.04575370             268  0.05901003 secs
    5    640    0.5  1     det 1.408054 0.02394448       640     0.05716385             347  0.10334420 secs
    6   1000    0.6  1     det 1.619144 0.02679712      1000     0.08121681             444  0.17916989 secs
    7   1778    0.7  1     det 1.898775 0.02610724      1778     0.12396933             604  0.44924521 secs
    8   4001    0.8  1     det 2.893992 0.03360215      4001     0.40735722             822  1.93757987 secs
    9  16001    0.9  1     det 5.122652 0.03724353     16001     1.24906572            1814 20.73342061 secs
    10   198    0.1  1     exp 1.222103 0.08242261       198     0.11062882             170  0.02579975 secs
    11   250    0.2  1     exp 1.380830 0.07831180       250     0.11687164             190  0.03336740 secs
    12   327    0.3  1     exp 1.430104 0.07252845       327     0.13069325             219  0.04917002 secs
    13   445    0.4  1     exp 1.918310 0.10522285       445     0.37295807             270  0.06772017 secs
    14   640    0.5  1     exp 2.033666 0.08003419       640     0.31915293             327  0.11357641 secs
    15  1000    0.6  1     exp 2.534484 0.07773870      1000     0.40575138             399  0.19646025 secs
    16  1778    0.7  1     exp 2.647708 0.06156654      1778     0.38229027             649  0.43771863 secs
    17  4001    0.8  1     exp 5.130352 0.08399659      4001     1.60501195             776  1.82564521 secs
    18 16001    0.9  1     exp 9.622424 0.08155676     16001     4.92104355            1786 21.01252937 secs
    

    Notons que pour bien faire, il faudrait au moins jeter le dernier batch qui est un peu dégénéré car on arrête d'injecter des tâches mais c'est négligeable car comme vous pouvez le voir, au final, on a quand même pas mal de batchs.

    library(ggplot2);
    myrate = function(x) {1/(1-x)} # la formule théorique 1/(mu-lambda) pour une vraie MM1
    ggplot(data=df, aes(x=lambda, y=resp, color=service)) + geom_point() + 
       geom_errorbar(aes(ymin=resp-2*resp_err, ymax=resp+2*resp_err),width=.02,alpha=.6) + 
       geom_errorbar(aes(ymin=resp-2*resp_err_batch, ymax=resp+2*resp_err_batch),width=.05) + 
        xlim(0,1)  + coord_cartesian(ylim = c(0,15)) +
        stat_function(fun = myrate, geom = "line",color="blue")
    

    MM1_new.png

    Les intervalles sont énormes et ce n'est pas dû uniquement au fait qu'on a bien moins de batchs que de tâches. C'est beaucoup dû au fait que le bon estimateur de la variance (qui pondère en prenant les batchs en compte) est très élevé. La variabilité d'une tâche à l'autre est vraiment très très importante quand on s'approche de la saturation du système.

  • Illustration de la réduction de variance grace aux variables antithétiques

    Tout d'abord, petite illustration du principe:

    f = function(x) { 1/(1+x) }
    U = runif(1000)
    V = runif(1000)
    var(f(U))
    var((f(U)+f(V))/2) # U et V sont indépendante donc la variance va être 2 fois plus faible
    var((f(U)+f(1-U))/2) # f(U) et f(1-U) sont correlés donc la variance est encore plus faible
    
    [1] 0.01970162
    [1] 0.0006373013
    [1] 0.009700396
    

    Le gain de variance est bien plus important que le fait d'avoir doublé le nombre de mesures. Mais comme certains l'ont remarqués, sur une fonction parfaitement symétrique, on ne gagnerait rien du tout.

    f = function(x) { (x-1/2)**2 }
    U = runif(1000)
    V = runif(1000)
    var(f(U))
    var((f(U)+f(V))/2) # U et V sont indépendante donc la variance va être 2 fois plus faible
    var((f(U)+f(1-U))/2) # f(U) et f(1-U) sont identiques donc la variance est exactement la même...
    
    [1] 0.00559895
    [1] 0.002784934
    [1] 0.00559895
    

    Voyons si ça aide pour notre simulation de file d'attente… Pour ça, on va avoir besoin de récupérer toute la trajectoire.

    head(MM1(N=20, lambda=.4, service="exp", seed=seed, output_state = T))
    
        arrival    service   response earlier_arrival  time_taken
    1 0.2226080 0.10089120 0.10089120       0.2226080 0.001753569
    2 0.3850868 1.97536865 1.97536865       0.1624788 0.001753569
    3 3.5132761 0.01117043 0.01117043       3.1281893 0.001753569
    4 3.9777521 0.05480658 0.05480658       0.4644760 0.001753569
    5 5.0866607 2.49571414 2.49571414       1.1089086 0.001753569
    6 6.7258271 0.66512007 1.52166786       1.6391663 0.001753569
    

    On a 2 variables antithétiques, essayons de toute les utiliser.

    lambda_control = .4;
    mu_control = 1;
    Ntot=10000
    seed = 42;
    
    df=data.frame()
    for(Ntot in c(1000,5000,10000)) {
        Base_Resp = MM1(N=Ntot, lambda=lambda_control, mu=mu_control, service="exp", seed=seed, output_state = T)
    
        Anti_Resp1 = MM1(N=Ntot/4,lambda=.4,service="exp", seed=seed, output_state = T)
        Anti_Resp2 = MM1(N=Ntot/4,lambda=.4,service="exp", seed=seed, output_state = T, antithetic_lambda = T)
        Anti_Resp3 = MM1(N=Ntot/4,lambda=.4,service="exp", seed=seed, output_state = T, antithetic_mu = T)
        Anti_Resp4 = MM1(N=Ntot/4,lambda=.4,service="exp", seed=seed, output_state = T, antithetic_lambda = T, antithetic_mu = T)
        Anti_Resp = (Anti_Resp1+Anti_Resp2+Anti_Resp3+Anti_Resp4)/4;
    
        df=rbind(df,data.frame(
                        n = Ntot,
                        mean_base = mean(Base_Resp$response),
                        se_base = sd(Base_Resp$response)/sqrt(length(Base_Resp$response)),
                        duration_base = head(Base_Resp$time_taken,n=1),
                        mean_anti = mean(Anti_Resp$response),
                        se_anti = sd(Anti_Resp$response)/sqrt(length(Anti_Resp$response)),
                        duration_anti = 4*head(Anti_Resp$time_taken,n=1)
                    ));
    
    }
    df
    
          n mean_base    se_base duration_base mean_anti    se_anti duration_anti
    1  1000  1.712526 0.05637287     0.1590347  1.668875 0.04010352    0.09941316
    2  5000  1.734233 0.02452186     2.2151918  1.727632 0.02113520    0.86967754
    3 10000  1.693142 0.01702518     7.5637808  1.646537 0.01487895    2.49059558
    

    On gagne d'autant plus que notre simulation est codée avec les pieds et que son temps d'exécution n'est pas linéaire… :( J'avoue ne pas trop savoir comment intégrer cette technique avec le calcul correct de l'intervalle de confiance lié à la regénération puisque les regénérations n'arrivent pas au même moment.

  • Illustration de la réduction de variance grace aux variables de contrôle
    df_control = MM1(N=1000, lambda=lambda_control, mu=mu_control, service="exp", seed=42, output_state = T)
    head(df_control)
    
        arrival   service  response earlier_arrival time_taken
    1 0.2226080 0.1645289 0.1645289       0.2226080  0.1640182
    2 0.3850868 2.7686552 2.7707053       0.1624788  0.1640182
    3 3.5132761 0.1986399 0.1986399       3.1281893  0.1640182
    4 3.9777521 0.6173715 0.6173715       0.4644760  0.1640182
    5 5.0866607 0.6951089 0.6951089       1.1089086  0.1640182
    6 6.7258271 3.8064331 3.8064331       1.6391663  0.1640182
    
    • Estimation naive

      Estimons (naïvement) l'espérance du temps de réponse:

      mean(df_control$response)
      sd(df_control$response)/sqrt(length(df_control$response))
      
      [1] 1.712526
      [1] 0.05637287
      
    • Exploitation du temps de service
      plot(df_control$service, df_control$response)
      

      MM1_cor_service.png

      cor(df_control$service, df_control$response)
      cov(df_control$service, df_control$response)
      
      [1] 0.5900534
      [1] 1.113425
      
      regularized_resp = df_control$response - (df_control$service - 1/mu_control)
      mean(regularized_resp)
      sd(regularized_resp)/sqrt(length(regularized_resp))
      
      [1] 1.67525
      [1] 0.04551393
      

      On peut donc très facilement gagner ≈ 20% sur la taille de l'intervalle de confiance en enlevant la variabilité du temps de service. Évidemment, ça ne marche que s'il y a de la variabilité sur le temps de service et que c'est ce qui explique le mieux la variabilité du temps de réponse. Quand λ est proche de 1, ce n'est clairement pas la cause la plus importante.

    • Exploitation du temps d'inter-arrivée

      Essayons de prendre en compte l'arrivée précédente.

      plot(df_control$earlier_arrival, df_control$response)
      

      MM1_cor_arrival.png

      cor(df_control$earlier_arrival, df_control$response)
      
      [1] -0.2243635
      

      La corrélation est bien plus faible. Pas clair que ça aide beaucoup mais essayons quand même. Je vous renvoie aux slides de Gersende FORT qui détaillent comment faire une bonne correction en exploitant intelligemment le lien entre les deux variables.

      b = cov(df_control$earlier_arrival, df_control$response) / (1/lambda_control**2)
      b
      
      [1] -0.1711974
      

      On peut alors construire un nouvel estimateur (vous pourrez essayez d'autres valeurs de b pour voir…):

      regularized_resp = df_control$response - (df_control$earlier_arrival - 1/lambda_control)
      mean(regularized_resp)
      sd(regularized_resp)/sqrt(length(regularized_resp))
      sd(df_control$response)/sqrt(length(df_control$response))
      
      [1] 1.586692
      [1] 0.1116894
      [1] 0.05637287
      

      Il y a l'air d'y avoir une réduction mais elle semble très faible.

      Évidemment, il y a d'autres variables avec qui le temps de réponse doit être bien plus correllée (par exemple, le nombre de tâches en attente quand une tâche arrive dans le système) mais on ne connaît pas grand chose sur leur statistique donc c'est difficile à exploiter.

DM

DM1: Mesures de performance de code

Rendu des étudiants

  • Anthony Geourjon (B)

    Bon, c'était pas bien fini… Il y a des données produites mais le lien vers les figures étaient cassées (et les figures n'étaient pas commitées de toutes façon) et du coup, il a fallu que je creuse pour comprendre.

    C'est très bien bien d'avoir essayé de randomiser les expériences. Par contre, c'est dommage que vous ne vous soyez pas débarassé de cette bouse de gnuplot pour passer en R, ça vous aurait simplifié la vie et vous aurait permi de repérer des problèmes. Histoire que vous en sachiez un peu plus sur vos mesures, voilà ce que j'ai fait à partir de vos données.

    library(ggplot2)
    library(dplyr)
    library(reshape2)
    df = read.csv("data/geourjoa-work-laptop_2017-02-09/measurements_09:05_wide.csv", header=T)
    head(df)
    df$XpId=1:length(df$Size)
    # Inspired by http://stackoverflow.com/questions/13250872/reshaping-data-to-plot-in-r-using-ggplot2 
    # or much better by http://seananderson.ca/2013/10/19/reshape.html
    df = df[!is.na(df$Size),] # il y a une valeur manquante au début, je ne sais pas pourquoi
    df2 <- melt(df, id.vars = c('Size','XpId'),
                variable.name = "Implementation", 
                value.name = "Time")
    head(df2)
    
        Size      Seq      Par     Libc
    1     NA 0.000000 0.000856 0.000004
    2 150000 0.062475 0.107621 0.053528
    3 150000 0.000000 0.000377 0.000002
    4 250000 0.099138 0.126824 0.089847
    5 250000 0.000000 0.000354 0.000002
    6 350000 0.150896 0.163016 0.130557
        Size XpId Implementation     Time
    1 150000    2            Seq 0.062475
    2 150000    3            Seq 0.000000
    3 250000    4            Seq 0.099138
    4 250000    5            Seq 0.000000
    5 350000    6            Seq 0.150896
    6 350000    7            Seq 0.000000
    
    ggplot(df2, aes(x=Size, y=Time, color=Implementation)) + geom_point()
    

    geourgon1.png

    Ouh là, il y a des points vraiment bizarres. Plus lentes, je comprendrais mais à ce point plus rapide, c'est vraiment étrange.

    head(df2[df2$Time<=.05,])
    
         Size XpId Implementation  Time
    2  150000    3            Seq 0e+00
    4  250000    5            Seq 0e+00
    6  350000    7            Seq 0e+00
    8  450000    9            Seq 1e-06
    10 550000   11            Seq 1e-06
    12 650000   13            Seq 1e-06
    

    Et ce n'a pas l'air d'être une série d'expériences bizarre, c'est un peu partout… Euh… encore qu'en y regardant de plus près, ça a l'air d'être exactement une expérience sur deux ?!? Honnêtement, je ne sais pas ce qu'il se passe sur votre machine! :) Et effectivement, quand je regarde dans le fichier de mesure initiale, ce phénomène est bien présent. Bon, y'a un problème mais admettons, j'enlève ces mesures très très bizarres.

    ggplot(df2[df2$XpId%%2==0,], aes(x=Size, y=Time, color=Implementation)) + 
        geom_point(alpha=.5) + geom_smooth(method="lm",formula=y~x*I(log(x)))
    

    geourgon2.png

    La version parallèle fini donc par surpasser les deux autres pour une taille de tableau supérieure à 65,000 sur cette machine et l'écart va aller en se creusant même si on semble loin du facteur 2 qu'on pourrait espérer. Je me demande quand même ce qu'il a bien pu se passer et pourquoi une mesure sur deux semble complètement aberrante…

    Ah!!! Je sais! :) J'ai ouvert votre "random.csv" et voilà à quoi il ressemble:

    1e+05
    150000
    2e+05
    250000
    3e+05
    350000
    4e+05
    450000
    

    La notation en E+05 est une notation pour de double, alors que le code de quicksort attend un int en entrée…

  • Maxime Chevalier (A+)
    • Très bien de bien noter les informations sur la machine.
    • Bon réflexe d'utiliser un moniteur ou un outil comme perf pour vérifier si vous utilisez bien la machine. Vous avez du vous rendre compte que ces informations n'étaient pas facile à interpréter.
    • Pensez à mettre la date plutôt que Test1, Test2, Test3…
    • Attention, dans Test3, vous utilisez à un moment geom_smooth (bien) et à un autre geom_line qui connecte les points par des lignes et n'est pas ce que vous voulez…
    • Pensez à bien commenter vos observations. Visiblement, la version parallèle est très très lente par rapport aux autres. J'ai l'impression que passer en -O3 n'a pas vraiment changé quoi que ce soit dans votre cas, non ?
    • Attention, je n'ai pas trouvé le code ayant généré les graphes présentés dans la section "Premiers résultats". Du coup, j'avais beaucoup de mal à savoir ce qu'ils représentaient exactement et comment les interpréter. Notamment, la partie avec la zone d'incertitude dont vous parlez, je ne suis pas sûr de savoir comment elle est calculée et ce qu'elle signifie. Je crains que ça soit une conséquence d'une erreur de syntaxe dans la commande R.
    • Ah, melt, super! Et la présentation après coup avec ggplot est bonne.
    • Effectivement, sur votre machine, l'algorithme parallèle n'est pas plus performant sur votre machine que le séquentiel en O3. Je peux vous assurer par contre que compiler en O3 ne générera pas de parallélisme permettant d'exploiter plusieurs cores.
  • Antoine Delise et Hugo Amodru A
    • la machine 1 a un nom qui manque d'originalité (Machine) mais 8 cores i7, c'est sympa pour ce type d'expériences.
    • Pénible les script shells pour lancer les expériences, non ? Un conseil: évitez de mélanger la génération avec le lanceur. Le rand du shell est souvent super pourri. Et utilisez le bon langage pour ce lancer les expériences (perl, python, …), ça vous simplifiera la vie.
    • "Un meilleur indice de visualisation: la moyenne". Argh! Mais pourquoi vous faites ça en awk ? Non seulement, c'est lent et peu précis mais c'est tellement facile à faire en R
    • Vous auriez vraiment du vous débarasser de cette cochonerie de gnuplot. C'est aussi tellement plus facile, fiable, souple en R. Par exemple, on a plusieurs couleurs/styles non coordonnés, on n'a pas d'intervalles de confiance, etc. Et c'est logique parcequ'avec gnuplot, c'est juste hyper pénible à faire… Bon, mais effectivement, sur les 2 machines, la version parallèle fini par accélérer les choses. On est loin d'un facteur 8 ou d'un facteur 4…
  • Timothée Lemaire (B+)
    • Vos commentaires sont OK même s'il n'y en a pas beaucoup.
    • Votre façon de calculer des moyennes est compliquée (et donc error-prone). Les data-frames et des outils comme dplyr permettent de calculer tout ça plus simplement. Je vous invite à regarder par exemple http://stackoverflow.com/questions/35953394/calculating-length-of-95-ci-using-dplyr.

      wget https://raw.githubusercontent.com/TimotheeLemaire/M2R-ParallelQuicksort/master/data/tim-ThinkPad-X220-Tablet_2017-01-20/measurements_02%3A11.csv -O /tmp/tim.csv
      
      library(dplyr)
      library(ggplot2)
      library(reshape2)
      df=read.csv("/tmp/tim.csv");
      df$XpId=1:length(df$Size)
      df2 <- melt(df, id.vars = c('Size','XpId'),
                  variable.name = "Implementation", 
                  value.name = "Time")
      head(df2)
      
      Attachement du package : ‘dplyr’
      
      The following objects are masked from ‘package:stats’:
      
          filter, lag
      
      The following objects are masked from ‘package:base’:
      
          intersect, setdiff, setequal, union
          Size XpId Implementation     Time
      1 200000    1            Seq 0.069050
      2 200000    2            Seq 0.044572
      3 200000    3            Seq 0.044920
      4 200000    4            Seq 0.044939
      5 200000    5            Seq 0.045102
      6 400000    6            Seq 0.096460
      
      ggplot(df2, aes(x=Size, y=Time, color=Implementation)) + 
          geom_point(alpha=.5) + geom_smooth(method="lm",formula=y~x*I(log(x)))
      

      tim1.png

      Wow, il y a trop de varabilité sur la version parallèle pour estimer quoi que ce soit correctement. Si vous aviez voulu calculer les intervalles de confiances correctement, voilà comment vous auriez pu faire (en vous inspirant de http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/ qui définit la fonction summarySE par exemple mais il y a moyen de faire plus simple et plus spécifique…):

      df2_summary <- summarySE(df2, measurevar="Time", groupvars=c("Size","Implementation"))
      ggplot(data=df2_summary, aes(x=Size, y=Time, color=Implementation)) + 
          geom_point(data=df2, alpha=.2) + 
          geom_errorbar(aes(ymin=Time-ci, ymax=Time+ci), width=.1,colour="black", width=100000) +
          geom_point(shape=21, fill="white") +
          geom_line(linetype=2)
      

      tim2.png

  • Thibaud Vegreville
  • Raphaël Fustes (B)

    Bon, sur votre machine, la version parallèle est incroyablement mauvaise. C'est bizarre mais pas absurde. Vous avez cherché mais vous ne m'indiquez pas ce que vous avez essayé, ni comment vous vous y êtes pris. Du coup, difficile de vous aider (et d'évaluer le temps que vous y avez passé…). Dans ces cas là, Il faut que vous essayiez vous même (des fois que le problème vienne de vous ;) sur une autre machine, que vous indiquiez bien quel est le type de votre machine, du système d'exploitation, etc.

  • Antoine Blanc (B)
    • Quand vous dites "Nous pouvons remarquer que les méilleurs performence sont obtenu en O2.", c'est étrange parce que quand je regarde vos courbes, c'est le O3 qui semble le plus efficace.
    • Pour générer au hasard, bien plus simple que ce que vous avez fait:

      sample(rep(c(100,1000,10000,100000,350000,500000),times=10))
      
       [1]   1000 100000    100 350000 500000 350000  10000   1000 350000 350000
      [11] 500000  10000   1000    100    100  10000    100  10000 350000 500000
      [21] 100000 350000   1000 100000 500000    100  10000 500000 350000    100
      [31] 500000  10000 350000 100000   1000 100000 100000 500000  10000   1000
      [41]   1000    100   1000 500000    100 350000 500000 100000  10000 350000
      [51] 500000  10000 100000  10000   1000 100000    100 100000    100   1000
      
    • Chez vous non plus la version parallèle est beaucoup moins efficace que la version séquentielle. Vous auriez pu essayer des tailles plus grandes, non?..
  • Mesh33 ??? (A)
    • Vous dites observer des "moyennes questionnables". C'est l'interprétation que vous en faites qui est questionnable et c'est exactement à cela que sert la notion d'intervalle de confiance. C'est d'autant plus important quand vous changez le nombre de tests au fil des expériences…
    • Même si j'ai beaucoup de mal à lire vos graphiques et à vous suivre, c'est intéressant d'avoir remarqué que votre tout premier test est plus long que les suivants. Ça peut parfaitement s'expliquer par la politique de DVFS de votre OS et par le turboboost, d'où l'importance de bien randomiser votre plan d'expérience…
    • "J’aimerais que l’échelle de temps reste fixe pour pouvoir comparer différents tests mais je ne sais pas faire." Plusieurs solutions. La plus simple consiste à rajouter un +ylim(0,10) de façon à imposer l'échelle sur l'axe des y. Mieux, vous fusionnez vos différentes data frame en indiquant préalablement le type d'expérience (dans une colonne supplémentaire) puis vous utilisez facet_wrap.
    • La libc ne profite effectivement pas des optimisation du compilateur puisque le code est dans une lib qui a été compilée une bonne fois pour toute.
    • Au final, il est difficile de savoir si Parellele est meilleure ou pas que Sequentiel, justement parce que vous n'avez pas d'intervalle de confiance ni d'information sur la variabilité.
    • "On peut déduire que O2 et O3 permettent de paralléliser le code séquentiel." Ou bien, vous pouvez en déduire qu'une fois O3 activé, la version parallèle n'apporte rien… C'est à mon avis le cas sur votre machine. À moins que vous n'ayez annoté votre code avec des pragmas OpenMP, les options de compilation ne vous permettront pas d'exploiter les différents coeurs de votre machine.
  • Ahmed Nassik (B-)

    Vous avez fait de nouvelles mesures et produit une figure mais il n'y a aucun commentaire… Dans vos tests, la version parallèle est toujours plus lente que les versions séquentielles. Vous auriez pu essayer des tailles plus grandes, non?..

DM2: Suite d'un des deux TDs au choix

Comme d'habitude, si vous êtes bloqués, vous avez le droit de chercher de l'aide auprès de vos collègues, du prof, du web, des TDs des années précédentes mais il faut le signaler, tout particulièrement si vous reprenez du code ou que vous vous en inspirez.

À rendre pour le 10 mars (mais si ça pose un problème, signalez le moi au plus tôt). Comme d'habitude, publication sur rpubs, envoi de l'URL par mail (uniquement à arnaud.legrand@imag.fr cette fois-ci) avec "[RICM4-EP]" dans le sujet du mail.

Sujet 1: Politiques de cache

On se met dans la même configuration que dans le TD sur la Modélisation et Simulation d'un Cache Web.

  1. Étudiez de la même façon (en simulation) la politique "Move-Ahead". Comme pour LRU, on suppose que notre cache est trié par date d'accès. Si un objet est déjà présent dans le cache, on le décale d'un cran vers la gauche. Si un objet n'est pas présent dans le cache, on le ramène en position K.
  2. Construisez une chaîne de Markov (spécifique à la taille du cache cette fois-ci, à la différence de ce qu'on a fait pour LRU) et retrouvez les probabilités précédentes (probabilité d'éviction de l'élément A et taux de défauts de cache).
  3. Comparez "Move-Ahead" à "Move-to-front/LRU" en fonction de la taille du cache.
  4. Faites une simulation à évènement discret pour LRU et pour Move-Ahead en supposant cette fois-ci que la popularité des objets suit une loi de puissance (i.e., \(P(i)\) proportionnel à \(\alpha^i\) avec \(\alpha \in ]0,1[\)).

La conclusion que vous devriez tirer de tout cela est que ces modèles, même très grossiers, vous donnent les bonnes tendances

Sujet 2: Comparaisons de méthodes pour calculer l'état stationnaire

Nous avons identifié trois méthodes pour calculer l'état stationnaire (marche aléatoire et statistique sur la fréquence de passage dans chaque état, itérer la fonction de transition à partir d'un vecteur de proba arbitraire, calculer directement les valeurs propres de la matrice).

Évaluez ces différentes méthodes sur:

  1. une ligne,
  2. un anneau,
  3. et une sucette.

Pour chaque cas, vous comparerez le temps de calcul et la qualité du résultat des trois méthodes. Vous interpréterez en quoi la forme du graphe influe sur la qualité de la solution.

Rendus des étudiants

  • Timothée Lemaire, DM1 (A+)

    Excellent travail. Le code est clair (même s'il manque un peu de fonctions) et les explications aussi. Les choix des conditions de simulation et de leur durée sont judicieux.

    • "les données de la simmulation ne semble vraiment pas significatifs. les résultats peuvent être faussé par les 2 défaut de cache résultant du choix pour l’initialisation du cache." Effectivement, ta simulation est biaisée par le début. Tu obtiens 0.002 avec un intervalle de confiance qui contient 0. La simulation n'est pas assez précise pour faire une estimation raisonnable. Et quand on voit la proba calculée avec la chaîne de Markov (1.2E-22), on comprend pourquoi. Il va falloir une simulation super longue pour faire une telle estimation et se débarasser de l'impact de l'effet initial.
  • Vincent Turrin, DM2 (http://rpubs.com/Mesh33/DM2_EP_VT) (B)

    Le code du début (simulation, calcul des probas ou des valeurs propres) est assez clair. Ce que tu as fait pour la partie mesure l'est beaucoup moins. Tu as une grosse fonction qui fait toutes les mesures plus le graphe, c'est une mauvaise idée. Il vaut mieux couper en petits bouts:

    • faire les mesures pour une méthode, en n'hésitant pas à faire plusieurs mesure (de suite dans un premier temps) de la même méthode pour avoir des mesures plus stables
    • les stoquer dans une data frame
    • visualiser jusqu'à ce qu'on trouve la bonne représentation, en mettant les bons labels aux graphes pour que le lecteur puisse comprendre de quoi il s'agit.
    • éventuellement refaire des mesures si on n'est pas satisfait (paramètres d'entrées, répétition, variabilité, …)
    • passer à la deuxième méthode
    • etc.

    Si tu as passé ton temps à utiliser cette énorme fonction, c'est certainement pour ça que ça prenait des plombes et que tu n'as pas pu aboutir. Passer au second use case sans avoir compris ce qui n'allait pas sur le premier n'a pas du aider.

  • Clément Rouquier, DM2 (http://rpubs.com/clement_rouquier/EP_ROUQUIER_DM2_SUJET2) A-
    • Burp. Indigeste ton code. Tu aimes ne pas me faciliter la tâche, hein ?
    • Argh, le copier coller immonde avec le même code de simulation et juste le code d'initialisation de la matrice qui change. Et d'ailleurs pourquoi inclure cette initialisation dans tes mesures ?
    • if(runif(1)>0.85){ #Permet de se sortir des buckets. 1/6 du temp on se TP à un random quoi qu'il arrive current = sample(c(1:N),1) }. WTF!?! Pourquoi tu te téléportes pendant ta marche aléatoire ? Ça fout toutes tes estimations en l'air.
    • Les mesures de temps et les graphes sont pas mal du tout (à ceci près que je ne vois pas du tout pourquoi tu utilise formula=y~x*I(log(x)). Je sens que j'ai ouvert la boite de Pandore en vous montrant ce truc là.
    • Tu évoques les problèmes de précision, mais tu ne les as pas du tout regardé. Le nombre d'itérations pour obtenir une bonne convergence ne sont évidemment pas du tout les même pour RW et pour FT.
    • Tu évoques l'influence de la topologie du graphe sur cette vitesse de convergence de RW. Tes intuitions sont bonnes. Dommage que tu n'ais pas fait les expériences correspondantes. Pour FT, la topologie a aussi un impact sur vitesse de convergence mais c'est bien plus rapide que pour FT. Quant à VP, effectivement, le temps de calcul ne dépend normalement pas de la structure du graphe si c'est une méthode directe qui est utilisée. Il peut y a voir des problèmes de stabilité numérique par contre.
    • Tiens, oui, c'est bizarre que tu as eu pour 4b. Tu as débranché ta machine ? :)
    • Effectivement, la structure n'a pas d'impact sur le temps de calcul.

    Regarde ce qu'a fait Hugo pour l'étude de la convergence.

  • Hugo Amodru-Favin, DM2 (http://rpubs.com/Stiglitz/257378) (A-)
    • Illustration de la vitesse de convergence vers l'état stationnaire. Très bien.
    • Par contre, tu n'as pas mesuré le temps réel de chaque méthode alors que c'était ça qui m'intéressait.
    • Je ne sais pas si tu l'as trouvé tout seul ou si c'est une réflexion collective mais bravo pour avoir remarqué qu'il y avait une formule générale.

    Regarde ce qu'a fait Jean-Baptiste pour l'étude de la performance des méthodes.

  • Nicolas Homberg, DM1 (http://rpubs.com/nicomy/DM2_EP) (B)
    • La comparaison (en Q3) de la dynamique de MA et LRU est une bonne idée. Mais je parlais plutôt des performances, i.e., des probas de défaut de cache sur A et des défauts de cache totaux.
    • Le reste est léger. Regarde le DM1 de Timothée Lemaire.
  • Nicolas Homberg, DM2 (http://rpubs.com/nicomy/DM2_EP) (C)
    • Tu fais des restarts fréquents et des trajectoires courtes. Pas sûr que ça soit une si bonne idée… Tu aurais du regarder l'impact de ces paramètres pour utiliser la bonne configuration.
    • Au final, tu ne compares pas les vitesses de convergence ni la cohérence des résultats entre les différentes méthodes, ni le temps de calcul.
  • Maxime Chevalier DM2 (http://rpubs.com/Chevamax/257196) A-
    • "On voit que pour un anneau, c’est très aléatoire. Augmentons le nombre de répétitions pour voir si on a atteind l’état stationnaire." Tu as regardé l'échelle ? ☺ Erm, même problème ailleurs.
    • "Cependant vu que la sucette est générée aléatoirement ont sait pas non plus comment elle est et donc comment interpreter." Il ne tenait qu'à toi de controler sa génération et de mettre par exemple la moitié des noeuds dans la ligne et l'autre moitié dans l'anneau…
    • Bon, c'est pas parfait mais tu as regardé la vitesse de convergence et évalué le temps de calcul.
  • Raphael Fustes DM2 (http://rpubs.com/FustesRpubs/257501) (A+)
    • Ton code manque de clarté et de structure mais fait ce qu'il faut.
    • Bon choix du nombre d'itérations quand tu compares RW et itération de la fonction de transition pour montrer que l'itération de la fonction de transition est bien plus efficace.

    Tu n'as pas regardé comment ça évoluait en fonction de N mais sinon, il y a tout ce qu'il faut.

  • Antoine Delise DM2 (http://rpubs.com/FustesRpubs/257501) (A+)
    • Le code est propre et clair.
    • L'étude de la vitesse de convergence est bonne. Effectivement, la convergence est plus lente mais ce n'est pas flagrant sur les courbes. Tu aurais du fixer l'échelle avec un ylim(0,NA).

    Beau travail.

  • Anthony Geourjon DM1 (B-)
    • Pour Q2 (calcul du régime stationnaire en passant par les valeurs propres de la matrice), si t'as pas compris comment faire en R, il faut regarder ce qu'on a fait dans le TD sur page rank. Et si tu n'as pas compris comment ça marchait et ce que ça signifiait, il faut venir me voir pour que je te fasses une explication adaptée.
    • En général, faire des set.seed dans les fonctions, c'est un peu dangereux (à moins d'avoir une vraie bonne raison de le faire, genre variables antithétiques et tout). Ça biaise et pas sûr que ça fasse vraiment ce que tu veux.
    • Dans tes graphes, l'utilisation de lines n'aide pas à savoir où sont les mesures et comme il n'y a pas d'intervalle de confiance, on a du mal à savoir s'il y a une différence entre la courbe rouge et la bleue.
    • Dans tes graphes, tu comptes le nombre de défauts de cache, qui dépend de la longueur de la trajectoire. Ça permet de mieux vérifier que c'est cohérent. Par exemple, si je comprends assez bien les deux premiers de Q3, je ne comprends pas grand chose aux deux suivants. D'ailleurs, c'est nvec et pas kvec que tu aurais du utiliser.
    • Dans l'ensemble, c'est assez fouilli, on a du mal à te suivre et tu n'as pas l'air bien sûr de ce que tu fais. C'est bien d'avoir essayé des configurations un peu différentes (par exemple plus d'items populaires) mais on a du mal à voir ce que tu en conclues.
  • Nassik DM1 (A-)
    • Dans l'ensemble, c'est assez clair.
    • La comparaison (en Q3) de la dynamique de MA et LRU est une bonne idée.
    • C'est une bonne idée d'avoir un seul graphique qui résume ta comparaison de MA et de LRU. Dommage qu'il n'y ait pas les intervalles de confiance associés. On voit bien qu'il y a quelque irrégularités qui n'ont pas lieu d'être.
    • Il manque la simulation d'une loi de puissance, non ?
  • Antoine Blanc DM2 (A-)
    • Vous avez considéré un N grand, c'est intéressant.
    • "Les autres formes une ligne non régulière. Ceci est surement du au fait qu’il n’y ai qu’un seul départ. Il faudrait faire plusieurs fois l’algorithme et faire la moyenne de toutes les courbes pour bien voire le résultat." Et bien pourquoi ne l'avez-vous pas fait ?
    • "Pour l’anneau, les points sont bien répandu sur le graphique. On ne peut pas déterminer la tendance de la courbe". Avez-vous vu l'échelle ? L'anneau est symétrique donc, toutes les valeurs doivent être égales. Le fait que
    • "Pour la sucette on constate qu’il y a tout les points (sauf un)". Sauf deux en fait…
    • Le calcul par itération est effectivement bien plus long que celui par valeurs propres mais vous avez fait bien trop d'itérations

Bibliographie