Correction des exercices du chapitre 4

Exercice 4.1 : Factorielle

Construction de la fonction factorielle avec la fonction prod:

factorielle <- function(n) {
  if (n<0) stop("l'entier doit être strictement positif")
  if (n==0) return(1)
  if (floor(n)!=n) warning(paste("arrondi de",n,"en",floor(n)))
  resultat <- prod(1:n)
  return(resultat)
}

Construction de la fonction factorielle avec une boucle for:

factorielle <- function(n) {
  if (n<0) stop("l'entier doit être strictement positif")
  if (n==0) return(1)
  if (floor(n)!=n){
    warning(paste("arrondi de",n,"en",floor(n)))
    n <- floor(n)
  }
  resultat <- 1
  for (i in 1:n) resultat <- resultat*i
  return(resultat)
}
factorielle(5)

[1] 120

factorielle(4.23)

[1] 24

Exercice 4.2 : Ventilation

Construisons la fonction de ventilation:

ventilation <- function(Xqual,p=0.05) {
  if (!is.factor(Xqual)) stop("Xqual doit etre un facteur \n")
  modalites <- levels(Xqual)
  if (length(modalites)<=1) stop("pas assez de niveaux \n")
  tabl <- table(Xqual)
  selecti <- (tabl/sum(tabl))<p
  if (!any(selecti)) return(Xqual) else {
   lesquels <- modalites[!selecti]
   prov <- factor(Xqual[(Xqual%in%lesquels)],levels=lesquels)
   prov <- table(prov)
   proba <- prov/sum(prov)
   for (j in modalites[selecti]) {
    ## tirages dans les modalites au hasard et remplacement
    if (length(lesquels)==1){
     warning("1 seule modalite\n")
     Xqual[Xqual==j]<-lesquels
    } else Xqual[Xqual==j]<-sample(lesquels,sum(Xqual==j),
             replace=TRUE,prob=proba)
   }
   Xqualvent <- factor(as.character(Xqual))
  }
  return(Xqualvent)
 }

Nous appliquons la fonction de la question précédente à chaque colonne du tableau qui est un facteur:

ventil.tab <- function (tab, seuil=0.05) {
 for (i in 1:ncol(tab)) {
  if (is.factor(tab[,i])) tab[,i]<-ventilation(tab[,i],p=seuil)
 }
 return(tab)
}

Exercice 4.3 : Ventilation sur facteur ordonné

La fonction est identique ou presque à ce qui a été vu dans la correction de l’exercice 2. Nous ajoutons la sortie (return) et quelques contrôles :

ventilation.ordonnee <- function(Xqual,p=0.05) {
  if (!is.ordered(Xqual)) stop("Xqual doit etre un ordonne \n")
  modalites <- levels(Xqual)
  if (length(modalites)<=1) stop("pas assez de niveaux \n")
  tabl <- table(Xqual)
  selecti <- (tabl/sum(tabl))<p
  if (!any(selecti)) return(Xqual) else {
    numero <- which(selecti)
    while(any((tabl/sum(tabl))<p)) {
      Xqual[,i] <- ventilation(Xqual[,i])
   }
 }
 return(Xqual)
}

Nous appliquons la fonction de la question précédente à chaque colonne du tableau qui est un facteur ordonné:

ventil.ordonnee.tab <- function (tab, seuil=0.05) {
 for (i in 1:ncol(tab)) {
  if (is.ordered(tab[,i])) tab[,i]<-ventilation.ordonnee(tab[,i],seuil)
 }
 return(tab)
}

Exercice 4.4 : Parallélisation

max_inc <- function(seq) {
  dseq <- seq[-1] > seq[-length(seq)]; cur_l <- 1+dseq[1]; max_l <- cur_l
  for (i in 1:(length(dseq)-1)) {
    if (dseq[i] && (dseq[i] == dseq[i+1])) {
    cur_l <- cur_l + 1;
    } else { cur_l <- 1+dseq[i+1] }
  max_l <- max(cur_l, max_l)
  }
  return(max_l)
}

Utilisation de sapply et table

n <- 1000
k <- 100000

system.time({
  set.seed(234)
  res <- sapply(1:k, function(aux){max_inc(runif(n))})
})

   user  system elapsed 
 80.474   0.152  81.267

table(res)

res
    4     5     6     7     8     9    10    11 
   82 30280 53828 13624  1905   248    32     1

Utilisation du package parallel

require(parallel)
system.time({
  cl <- makeCluster(detectCores()-1)    
  clusterExport(cl, varlist=c("n","k","max_inc"))
  res <- parSapply(cl, 1:k, function(aux){max_inc(runif(n))})
  stopCluster(cl)
})

   user  system elapsed 
  0.137   0.150  50.554

table(res)

res
    4     5     6     7     8     9    10    11    12 
   78 30530 53413 13790  1909   253    23     2     2

Pour accelérer encore les temps de calcul, il faut utiliser le package Rcpp.

#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]
int max_inc_rcpp(NumericVector x) {
  int max_l = 0;
  int cur_l = 0;
  bool cur_g, prev_g;
  cur_g = x[1] > x[0];
  cur_l = 1 + cur_g;
  max_l = cur_l;
  for (int i = 1; i < (x.length()-1); i++) {
    prev_g = cur_g;
    cur_g = x[i+1] > x[i];
    if (cur_g && (prev_g == cur_g)) {
      cur_l = cur_l + 1;
    } else {
      cur_l = 1 + cur_g;
    }
    if (cur_l > max_l) { max_l = cur_l;};
  }
  return max_l;
}

system.time(res <- sapply(1:k, function(i) { max_inc_rcpp(runif(n)) }))

   user  system elapsed 
  4.650   0.020   4.698

Exercice 4.5 :

Simulons deux variables x et y:

set.seed(1234)
x <- rnorm(50)
epsilon <- rnorm(50)
y <- 2+ 0.5*x + epsilon

On calule le coefficient entre x et y et le test associé, puis on compte le nombre de fois où le coefficient de corrélation entre x et y permuté est plus grand que le coefficient de corrélation observé:

cor.xy=cor(x,y)
cor.test(x,y)


    Pearson's product-moment correlation

data:  x and y
t = 2.5274, df = 48, p-value = 0.01484
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.07114562 0.56697017
sample estimates:
      cor 
0.3427066

nbsim=100000

system.time({
set.seed(234)
som=0
for (i in 1:nbsim) som=som+(cor(x,y[sample(length(y))])>cor.xy)
})

   user  system elapsed 
  6.197   0.003   6.243

print(som/nbsim*2)

[1] 0.01476

La probabilité critique obtenue par le test classique, 0.01484, est proche de celle calculée par permutation 0.01476.

On utilise le package parallel et la fonction parLapply

require(parallel)
cl <- makeCluster(3)
clusterExport(cl, varlist=c("x","y","cor.xy"))
system.time({
  res <- parSapply(cl, 1:nbsim, function(aux){1*(cor(x,y[sample(length(y))])>cor.xy)})
})

   user  system elapsed 
  0.057   0.011   4.178

stopCluster(cl)
mean(res)*2

[1] 0.0142

On multiplie par 2 le pourcentage obtenu car le test est bilatéral. On note que le temps de calcul est en gros divisé par 2.6.

On relance les lignes de code avec \(rchisq\) pour générer les aléas:

set.seed(1234)
x <- rchisq(50, df=1)
epsilon <- rchisq(50, df=1)
y <- 2+ 0.5*x + epsilon

cor.xy=cor(x,y)
cor.test(x,y)


    Pearson's product-moment correlation

data:  x and y
t = 2.7855, df = 48, p-value = 0.007628
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1056607 0.5901113
sample estimates:
      cor 
0.3730372

nbsim=100000

system.time({
set.seed(234)
som=0
for (i in 1:nbsim) som=som+(cor(x,y[sample(length(y))])>cor.xy)
})

   user  system elapsed 
  6.560   0.008   6.609

print(som/nbsim*2)

[1] 0.0252

require(parallel)
cl <- makeCluster(3)
clusterExport(cl, varlist=c("x","y","cor.xy"))
system.time({
  set.seed(234)
  res <- parSapply(cl, 1:nbsim, function(aux){1*(cor(x,y[sample(length(y))])>cor.xy)})
})

   user  system elapsed 
  0.062   0.016   4.228

stopCluster(cl)
2*mean(res)

[1] 0.02406

La probabilité critique obtenue avec le test classique, 0.007628, est cette fois assez différente de celle obtenue par le test de permutation. En effet, l’hypothèse de normalité n’est pas vérifiée ici et le test classique n’est pas adaptée. On aura donc plus confiance dans l’estimation de la probabilité critique obtenue par le test de permutation : 0.025

LS0tCnRpdGxlOiAiQ29ycmVjdGlvbiBkZXMgZXhlcmNpY2VzIGR1IGNoYXBpdHJlIDQiCmF1dGhvcjogIkh1c3NvbiBldCBhbC4iCmRhdGU6ICIwOS8wOS8yMDE4IgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogeWVzCiAgICB0b2NfZGVwdGg6IDMKICAgIHRvY19mbG9hdDogeWVzCiAgaHRtbF9kb2N1bWVudDoKICAgIHRvYzogeWVzCiAgICB0b2NfZGVwdGg6ICczJwogICAgdG9jX2Zsb2F0OiB5ZXMKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCBjYWNoZSA9IFRSVUUpCmBgYAoKIyMgRXhlcmNpY2UgNC4xIDogRmFjdG9yaWVsbGUKCjEuIENvbnN0cnVjdGlvbiBkZSBsYSBmb25jdGlvbiBmYWN0b3JpZWxsZSBhdmVjIGxhIGZvbmN0aW9uICpwcm9kKjoKCmBgYHtyLG1lc3NhZ2U9RkFMU0Usd2FybmluZz1GQUxTRX0KZmFjdG9yaWVsbGUgPC0gZnVuY3Rpb24obikgewogIGlmIChuPDApIHN0b3AoImwnZW50aWVyIGRvaXQgw6p0cmUgc3RyaWN0ZW1lbnQgcG9zaXRpZiIpCiAgaWYgKG49PTApIHJldHVybigxKQogIGlmIChmbG9vcihuKSE9bikgd2FybmluZyhwYXN0ZSgiYXJyb25kaSBkZSIsbiwiZW4iLGZsb29yKG4pKSkKICByZXN1bHRhdCA8LSBwcm9kKDE6bikKICByZXR1cm4ocmVzdWx0YXQpCn0KYGBgCgoyLiBDb25zdHJ1Y3Rpb24gZGUgbGEgZm9uY3Rpb24gZmFjdG9yaWVsbGUgYXZlYyB1bmUgYm91Y2xlICpmb3IqOgoKYGBge3IsbWVzc2FnZT1GQUxTRSx3YXJuaW5nPUZBTFNFfQpmYWN0b3JpZWxsZSA8LSBmdW5jdGlvbihuKSB7CiAgaWYgKG48MCkgc3RvcCgibCdlbnRpZXIgZG9pdCDDqnRyZSBzdHJpY3RlbWVudCBwb3NpdGlmIikKICBpZiAobj09MCkgcmV0dXJuKDEpCiAgaWYgKGZsb29yKG4pIT1uKXsKICAgIHdhcm5pbmcocGFzdGUoImFycm9uZGkgZGUiLG4sImVuIixmbG9vcihuKSkpCiAgICBuIDwtIGZsb29yKG4pCiAgfQogIHJlc3VsdGF0IDwtIDEKICBmb3IgKGkgaW4gMTpuKSByZXN1bHRhdCA8LSByZXN1bHRhdCppCiAgcmV0dXJuKHJlc3VsdGF0KQp9CmZhY3RvcmllbGxlKDUpCmZhY3RvcmllbGxlKDQuMjMpCmBgYAojIyBFeGVyY2ljZSA0LjIgOiBWZW50aWxhdGlvbgoKMS4gQ29uc3RydWlzb25zIGxhIGZvbmN0aW9uIGRlIHZlbnRpbGF0aW9uOgpgYGB7cixtZXNzYWdlPUZBTFNFLHdhcm5pbmc9RkFMU0V9CnZlbnRpbGF0aW9uIDwtIGZ1bmN0aW9uKFhxdWFsLHA9MC4wNSkgewogIGlmICghaXMuZmFjdG9yKFhxdWFsKSkgc3RvcCgiWHF1YWwgZG9pdCBldHJlIHVuIGZhY3RldXIgXG4iKQogIG1vZGFsaXRlcyA8LSBsZXZlbHMoWHF1YWwpCiAgaWYgKGxlbmd0aChtb2RhbGl0ZXMpPD0xKSBzdG9wKCJwYXMgYXNzZXogZGUgbml2ZWF1eCBcbiIpCiAgdGFibCA8LSB0YWJsZShYcXVhbCkKICBzZWxlY3RpIDwtICh0YWJsL3N1bSh0YWJsKSk8cAogIGlmICghYW55KHNlbGVjdGkpKSByZXR1cm4oWHF1YWwpIGVsc2UgewogICBsZXNxdWVscyA8LSBtb2RhbGl0ZXNbIXNlbGVjdGldCiAgIHByb3YgPC0gZmFjdG9yKFhxdWFsWyhYcXVhbCVpbiVsZXNxdWVscyldLGxldmVscz1sZXNxdWVscykKICAgcHJvdiA8LSB0YWJsZShwcm92KQogICBwcm9iYSA8LSBwcm92L3N1bShwcm92KQogICBmb3IgKGogaW4gbW9kYWxpdGVzW3NlbGVjdGldKSB7CiAgICAjIyB0aXJhZ2VzIGRhbnMgbGVzIG1vZGFsaXRlcyBhdSBoYXNhcmQgZXQgcmVtcGxhY2VtZW50CiAgICBpZiAobGVuZ3RoKGxlc3F1ZWxzKT09MSl7CiAgICAgd2FybmluZygiMSBzZXVsZSBtb2RhbGl0ZVxuIikKICAgICBYcXVhbFtYcXVhbD09al08LWxlc3F1ZWxzCiAgICB9IGVsc2UgWHF1YWxbWHF1YWw9PWpdPC1zYW1wbGUobGVzcXVlbHMsc3VtKFhxdWFsPT1qKSwKICAgICAgICAgICAgIHJlcGxhY2U9VFJVRSxwcm9iPXByb2JhKQogICB9CiAgIFhxdWFsdmVudCA8LSBmYWN0b3IoYXMuY2hhcmFjdGVyKFhxdWFsKSkKICB9CiAgcmV0dXJuKFhxdWFsdmVudCkKIH0KYGBgCjIuIE5vdXMgYXBwbGlxdW9ucyBsYSBmb25jdGlvbiBkZSBsYSBxdWVzdGlvbiBwcsOpY8OpZGVudGUgw6AgY2hhcXVlIGNvbG9ubmUgZHUgdGFibGVhdSBxdWkgZXN0IHVuIGZhY3RldXI6CmBgYHtyLG1lc3NhZ2U9RkFMU0Usd2FybmluZz1GQUxTRX0KdmVudGlsLnRhYiA8LSBmdW5jdGlvbiAodGFiLCBzZXVpbD0wLjA1KSB7CiBmb3IgKGkgaW4gMTpuY29sKHRhYikpIHsKICBpZiAoaXMuZmFjdG9yKHRhYlssaV0pKSB0YWJbLGldPC12ZW50aWxhdGlvbih0YWJbLGldLHA9c2V1aWwpCiB9CiByZXR1cm4odGFiKQp9CmBgYAoKIyMgRXhlcmNpY2UgNC4zIDogVmVudGlsYXRpb24gc3VyIGZhY3RldXIgb3Jkb25uw6kKCjEuIExhIGZvbmN0aW9uIGVzdCBpZGVudGlxdWUgb3UgcHJlc3F1ZSDDoCBjZSBxdWkgYSDDqXTDqSB2dSBkYW5zIGxhIGNvcnJlY3Rpb24gZGUgbCdleGVyY2ljZSAyLgpOb3VzIGFqb3V0b25zIGxhIHNvcnRpZSAoKnJldHVybiopIGV0IHF1ZWxxdWVzIGNvbnRyw7RsZXMgOgpgYGB7cixtZXNzYWdlPUZBTFNFLHdhcm5pbmc9RkFMU0V9CnZlbnRpbGF0aW9uLm9yZG9ubmVlIDwtIGZ1bmN0aW9uKFhxdWFsLHA9MC4wNSkgewogIGlmICghaXMub3JkZXJlZChYcXVhbCkpIHN0b3AoIlhxdWFsIGRvaXQgZXRyZSB1biBvcmRvbm5lIFxuIikKICBtb2RhbGl0ZXMgPC0gbGV2ZWxzKFhxdWFsKQogIGlmIChsZW5ndGgobW9kYWxpdGVzKTw9MSkgc3RvcCgicGFzIGFzc2V6IGRlIG5pdmVhdXggXG4iKQogIHRhYmwgPC0gdGFibGUoWHF1YWwpCiAgc2VsZWN0aSA8LSAodGFibC9zdW0odGFibCkpPHAKICBpZiAoIWFueShzZWxlY3RpKSkgcmV0dXJuKFhxdWFsKSBlbHNlIHsKICAgIG51bWVybyA8LSB3aGljaChzZWxlY3RpKQogICAgd2hpbGUoYW55KCh0YWJsL3N1bSh0YWJsKSk8cCkpIHsKICAgICAgWHF1YWxbLGldIDwtIHZlbnRpbGF0aW9uKFhxdWFsWyxpXSkKICAgfQogfQogcmV0dXJuKFhxdWFsKQp9CmBgYAoyLiBOb3VzIGFwcGxpcXVvbnMgbGEgZm9uY3Rpb24gZGUgbGEgcXVlc3Rpb24gcHLDqWPDqWRlbnRlIMOgIGNoYXF1ZSBjb2xvbm5lIGR1IHRhYmxlYXUgcXVpIGVzdCB1biBmYWN0ZXVyIG9yZG9ubsOpOgpgYGB7cixtZXNzYWdlPUZBTFNFLHdhcm5pbmc9RkFMU0V9CnZlbnRpbC5vcmRvbm5lZS50YWIgPC0gZnVuY3Rpb24gKHRhYiwgc2V1aWw9MC4wNSkgewogZm9yIChpIGluIDE6bmNvbCh0YWIpKSB7CiAgaWYgKGlzLm9yZGVyZWQodGFiWyxpXSkpIHRhYlssaV08LXZlbnRpbGF0aW9uLm9yZG9ubmVlKHRhYlssaV0sc2V1aWwpCiB9CiByZXR1cm4odGFiKQp9CmBgYAoKIyMgRXhlcmNpY2UgNC40IDogUGFyYWxsw6lsaXNhdGlvbgoKYGBge3J9Cm1heF9pbmMgPC0gZnVuY3Rpb24oc2VxKSB7CiAgZHNlcSA8LSBzZXFbLTFdID4gc2VxWy1sZW5ndGgoc2VxKV07IGN1cl9sIDwtIDErZHNlcVsxXTsgbWF4X2wgPC0gY3VyX2wKICBmb3IgKGkgaW4gMToobGVuZ3RoKGRzZXEpLTEpKSB7CiAgICBpZiAoZHNlcVtpXSAmJiAoZHNlcVtpXSA9PSBkc2VxW2krMV0pKSB7CiAgICBjdXJfbCA8LSBjdXJfbCArIDE7CiAgICB9IGVsc2UgeyBjdXJfbCA8LSAxK2RzZXFbaSsxXSB9CiAgbWF4X2wgPC0gbWF4KGN1cl9sLCBtYXhfbCkKICB9CiAgcmV0dXJuKG1heF9sKQp9CmBgYAoKMS4gVXRpbGlzYXRpb24gZGUgc2FwcGx5IGV0IHRhYmxlCgpgYGB7cn0KbiA8LSAxMDAwCmsgPC0gMTAwMDAwCmBgYAoKYGBge3J9CnN5c3RlbS50aW1lKHsKICBzZXQuc2VlZCgyMzQpCiAgcmVzIDwtIHNhcHBseSgxOmssIGZ1bmN0aW9uKGF1eCl7bWF4X2luYyhydW5pZihuKSl9KQp9KQp0YWJsZShyZXMpCmBgYAoyLiBVdGlsaXNhdGlvbiBkdSBwYWNrYWdlIHBhcmFsbGVsCgpgYGB7cixtZXNzYWdlPUZBTFNFLHdhcm5pbmc9RkFMU0V9CnJlcXVpcmUocGFyYWxsZWwpCnN5c3RlbS50aW1lKHsKICBjbCA8LSBtYWtlQ2x1c3RlcihkZXRlY3RDb3JlcygpLTEpICAgIAogIGNsdXN0ZXJFeHBvcnQoY2wsIHZhcmxpc3Q9YygibiIsImsiLCJtYXhfaW5jIikpCiAgcmVzIDwtIHBhclNhcHBseShjbCwgMTprLCBmdW5jdGlvbihhdXgpe21heF9pbmMocnVuaWYobikpfSkKICBzdG9wQ2x1c3RlcihjbCkKfSkKdGFibGUocmVzKQpgYGAKCjMuIFBvdXIgYWNjZWzDqXJlciBlbmNvcmUgbGVzIHRlbXBzIGRlIGNhbGN1bCwgaWwgZmF1dCB1dGlsaXNlciBsZSBwYWNrYWdlIFJjcHAuCgpgYGB7UmNwcH0KI2luY2x1ZGUgPFJjcHAuaD4KdXNpbmcgbmFtZXNwYWNlIFJjcHA7CgoKLy8gW1tSY3BwOjpleHBvcnRdXQppbnQgbWF4X2luY19yY3BwKE51bWVyaWNWZWN0b3IgeCkgewogIGludCBtYXhfbCA9IDA7CiAgaW50IGN1cl9sID0gMDsKICBib29sIGN1cl9nLCBwcmV2X2c7CiAgY3VyX2cgPSB4WzFdID4geFswXTsKICBjdXJfbCA9IDEgKyBjdXJfZzsKICBtYXhfbCA9IGN1cl9sOwogIGZvciAoaW50IGkgPSAxOyBpIDwgKHgubGVuZ3RoKCktMSk7IGkrKykgewogICAgcHJldl9nID0gY3VyX2c7CiAgICBjdXJfZyA9IHhbaSsxXSA+IHhbaV07CiAgICBpZiAoY3VyX2cgJiYgKHByZXZfZyA9PSBjdXJfZykpIHsKICAgICAgY3VyX2wgPSBjdXJfbCArIDE7CiAgICB9IGVsc2UgewogICAgICBjdXJfbCA9IDEgKyBjdXJfZzsKICAgIH0KICAgIGlmIChjdXJfbCA+IG1heF9sKSB7IG1heF9sID0gY3VyX2w7fTsKICB9CiAgcmV0dXJuIG1heF9sOwp9CmBgYAoKYGBge3J9CnN5c3RlbS50aW1lKHJlcyA8LSBzYXBwbHkoMTprLCBmdW5jdGlvbihpKSB7IG1heF9pbmNfcmNwcChydW5pZihuKSkgfSkpCmBgYAoKIyMgRXhlcmNpY2UgNC41IDogCgoxLiBTaW11bG9ucyBkZXV4IHZhcmlhYmxlcyB4IGV0IHk6CgpgYGB7cn0Kc2V0LnNlZWQoMTIzNCkKeCA8LSBybm9ybSg1MCkKZXBzaWxvbiA8LSBybm9ybSg1MCkKeSA8LSAyKyAwLjUqeCArIGVwc2lsb24KYGBgCgoyLiBPbiBjYWx1bGUgbGUgY29lZmZpY2llbnQgZW50cmUgeCBldCB5IGV0IGxlIHRlc3QgYXNzb2Npw6ksIHB1aXMgb24gY29tcHRlIGxlIG5vbWJyZSBkZSBmb2lzIG/DuSBsZSBjb2VmZmljaWVudCBkZSAgY29ycsOpbGF0aW9uIGVudHJlIHggZXQgeSBwZXJtdXTDqSBlc3QgcGx1cyBncmFuZCBxdWUgbGUgY29lZmZpY2llbnQgZGUgY29ycsOpbGF0aW9uIG9ic2VydsOpOgoKYGBge3J9CmNvci54eT1jb3IoeCx5KQpjb3IudGVzdCh4LHkpCm5ic2ltPTEwMDAwMAoKc3lzdGVtLnRpbWUoewpzZXQuc2VlZCgyMzQpCnNvbT0wCmZvciAoaSBpbiAxOm5ic2ltKSBzb209c29tKyhjb3IoeCx5W3NhbXBsZShsZW5ndGgoeSkpXSk+Y29yLnh5KQp9KQpwcmludChzb20vbmJzaW0qMikKYGBgCkxhIHByb2JhYmlsaXTDqSBjcml0aXF1ZSBvYnRlbnVlIHBhciBsZSB0ZXN0IGNsYXNzaXF1ZSwgMC4wMTQ4NCwgZXN0IHByb2NoZSBkZSBjZWxsZSBjYWxjdWzDqWUgcGFyIHBlcm11dGF0aW9uIDAuMDE0NzYuCgozLiBPbiB1dGlsaXNlIGxlIHBhY2thZ2UgcGFyYWxsZWwgZXQgbGEgZm9uY3Rpb24gcGFyTGFwcGx5CgpgYGB7cixtZXNzYWdlPUZBTFNFLHdhcm5pbmc9RkFMU0V9CnJlcXVpcmUocGFyYWxsZWwpCmNsIDwtIG1ha2VDbHVzdGVyKDMpCmNsdXN0ZXJFeHBvcnQoY2wsIHZhcmxpc3Q9YygieCIsInkiLCJjb3IueHkiKSkKc3lzdGVtLnRpbWUoewogIHJlcyA8LSBwYXJTYXBwbHkoY2wsIDE6bmJzaW0sIGZ1bmN0aW9uKGF1eCl7MSooY29yKHgseVtzYW1wbGUobGVuZ3RoKHkpKV0pPmNvci54eSl9KQp9KQpzdG9wQ2x1c3RlcihjbCkKbWVhbihyZXMpKjIKYGBgCk9uIG11bHRpcGxpZSBwYXIgMiBsZSBwb3VyY2VudGFnZSBvYnRlbnUgY2FyIGxlIHRlc3QgZXN0IGJpbGF0w6lyYWwuCk9uIG5vdGUgcXVlIGxlIHRlbXBzIGRlIGNhbGN1bCBlc3QgZW4gZ3JvcyBkaXZpc8OpIHBhciAyLjYuCgo0LiBPbiByZWxhbmNlIGxlcyBsaWduZXMgZGUgY29kZSBhdmVjICRyY2hpc3EkIHBvdXIgZ8OpbsOpcmVyIGxlcyBhbMOpYXM6CgpgYGB7cn0Kc2V0LnNlZWQoMTIzNCkKeCA8LSByY2hpc3EoNTAsIGRmPTEpCmVwc2lsb24gPC0gcmNoaXNxKDUwLCBkZj0xKQp5IDwtIDIrIDAuNSp4ICsgZXBzaWxvbgoKY29yLnh5PWNvcih4LHkpCmNvci50ZXN0KHgseSkKbmJzaW09MTAwMDAwCgpzeXN0ZW0udGltZSh7CnNldC5zZWVkKDIzNCkKc29tPTAKZm9yIChpIGluIDE6bmJzaW0pIHNvbT1zb20rKGNvcih4LHlbc2FtcGxlKGxlbmd0aCh5KSldKT5jb3IueHkpCn0pCnByaW50KHNvbS9uYnNpbSoyKQoKcmVxdWlyZShwYXJhbGxlbCkKY2wgPC0gbWFrZUNsdXN0ZXIoMykKY2x1c3RlckV4cG9ydChjbCwgdmFybGlzdD1jKCJ4IiwieSIsImNvci54eSIpKQpzeXN0ZW0udGltZSh7CiAgc2V0LnNlZWQoMjM0KQogIHJlcyA8LSBwYXJTYXBwbHkoY2wsIDE6bmJzaW0sIGZ1bmN0aW9uKGF1eCl7MSooY29yKHgseVtzYW1wbGUobGVuZ3RoKHkpKV0pPmNvci54eSl9KQp9KQpzdG9wQ2x1c3RlcihjbCkKMiptZWFuKHJlcykKYGBgCkxhIHByb2JhYmlsaXTDqSBjcml0aXF1ZSBvYnRlbnVlIGF2ZWMgbGUgdGVzdCBjbGFzc2lxdWUsIDAuMDA3NjI4LCBlc3QgY2V0dGUgZm9pcyBhc3NleiBkaWZmw6lyZW50ZSBkZSBjZWxsZSBvYnRlbnVlIHBhciBsZSB0ZXN0IGRlIHBlcm11dGF0aW9uLiBFbiBlZmZldCwgbCdoeXBvdGjDqHNlIGRlIG5vcm1hbGl0w6kgbidlc3QgcGFzIHbDqXJpZmnDqWUgaWNpIGV0IGxlIHRlc3QgY2xhc3NpcXVlIG4nZXN0IHBhcyBhZGFwdMOpZS4gT24gYXVyYSBkb25jIHBsdXMgY29uZmlhbmNlIGRhbnMgbCdlc3RpbWF0aW9uIGRlIGxhIHByb2JhYmlsaXTDqSBjcml0aXF1ZSBvYnRlbnVlIHBhciBsZSB0ZXN0IGRlIHBlcm11dGF0aW9uIDogMC4wMjU=