RICM4: Évaluation de Performance 2016

Modélisation et principe de simulation

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

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

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

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)

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

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

plot(data=dd, resp ~ lambda)

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

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

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

  

  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)
    

    

    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)
    

    

    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.

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

  

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

  

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

    

    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)
    

    

• 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?..

Sujet 1: Politiques de cache

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

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 (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