This John Snow knew something!

1 Visualizzare i dati

La Visualizzazione dei dati permettere di individuare trend, connessioni e carpire informazioni dai dati che in forma tabulare non sono ovvie.

I comandi che seguono aiutano in tali indagini. Tutti questi comandi possono essere personalizzati usando diversi colori, personalizzando le label degli assi o affiancando le immagini.

Un pacchetto che permette non solo di produrre immagini (un po’ più belle), ma anche di utilizzare funzioni avanzate è ggplot2 che però non tratteremo in questo corso.

Oltre ad introdurre alcune funzioni che permettono di ottenere dei grafici, esploreremo anche alcuni dataset.

1.1 Box plot

A lezione abbiamo già incontrato i boxplot. Usando R è molto semplice produrre dei boxplot usando gli appositi comandi. Vediamoli usando il dataset birthwt, contenuto in MASS.

library(MASS)
boxplot(birthwt$bwt)

La seguente notazione permette di produrre più boxplot relativi a bwt, divisi rispetto alla variabile race.

boxplot(bwt ~ race, data = birthwt )

Esercizio Si investighi la variabile lwt rispetto a race e smoke.

Esercizio Si noti che il pacchetto car contiene una funzione chiamata Boxplot(), con la maiuscola, che vicino ai potenziali outliers indica l’indice dell’elemento. Si replichino i precedenti boxplot usando tale funzione.

1.2 Istogrammi e scatter plost

Usando il comando hist() è possibile ottenere degli istogrammi per le variabili, passando gli argomenti in modo analogo ai boxplot.

Analogamente, con i comandi plot(x,y) or scatterplot(x,y) è possibile rappresentare i dati come punti.

Vedremo che questo può essere molto utile per ottenere degli indizi sulle relazioni tra variabili.

2 BodyTemperature

BodyTemperature è un dataset che contiene la misurazione di quattro variabili (sesso, età, frequenza cardiaca e temperatura corporea ascellare) di 100 pazienti. Possiamo leggere importare il file specificando invece che il percorso, l’indirizzo url.

BodyTemperature <- read.csv(url('http://extras.springer.com/2012/978-1-4614-1301-1/BodyTemperature.txt'))

Se però visualizziamo questo dataset qualosa non torna. Il comando giusto specifica quale separatore viene considerato: in questo caso lo spazio, che è diverso da quello di default di read.csv. In modo analogo avremmo potuto leggere i dati usando read.table(). Tuttavia, controllaimo bene l’output per capire se effettivamente il dataset è stato letto nel modo giusto!

(BodyTemperature <- read.table(url('http://extras.springer.com/2012/978-1-4614-1301-1/BodyTemperature.txt')) )
BodyTemperature <- read.csv(url('http://extras.springer.com/2012/978-1-4614-1301-1/BodyTemperature.txt'), sep= " ")

Esploriamo un po’ il dataset per capirne il contenuto e le variabili

str(BodyTemperature)
'data.frame':   100 obs. of  5 variables:
 $ Gender     : Factor w/ 2 levels "F","M": 2 2 2 1 1 2 1 1 1 2 ...
 $ Age        : int  33 32 42 33 26 37 32 45 31 49 ...
 $ HeartRate  : int  69 72 68 75 68 79 71 73 77 81 ...
 $ Temperature: num  97 98.8 96.2 97.8 98.8 ...
 $ TempC      : num  36.1 37.1 35.7 36.6 37.1 ...
head(BodyTemperature)
names(BodyTemperature)
[1] "Gender"      "Age"         "HeartRate"   "Temperature"
[5] "TempC"      

Solo a scopo didattico, vediamo l’applicazione di uno dei comandi visti: levels():

BodyTemperature$GenderLong <- BodyTemperature$Gender
levels(BodyTemperature$GenderLong) <- c("Female", "Male")
summary(BodyTemperature)
 Gender      Age          HeartRate      Temperature    
 F:51   Min.   :21.00   Min.   :61.00   Min.   : 96.20  
 M:49   1st Qu.:33.75   1st Qu.:69.00   1st Qu.: 97.70  
        Median :37.00   Median :73.00   Median : 98.30  
        Mean   :37.62   Mean   :73.66   Mean   : 98.33  
        3rd Qu.:42.00   3rd Qu.:78.00   3rd Qu.: 98.90  
        Max.   :50.00   Max.   :87.00   Max.   :101.30  
     TempC      
 Min.   :35.67  
 1st Qu.:36.50  
 Median :36.83  
 Mean   :36.85  
 3rd Qu.:37.17  
 Max.   :38.50  
range(BodyTemperature[,-1])
[1]  21.0 101.3

In modo analogo possiamo cercare di ottenere la distanza interquartile per i gli elementi del dataset con IQR(BodyTemperature[,2:4]), tuttavia questo ci da errore.

Possiamo aggirare questa limitazione usando il comando apply(). La stessa funzione la possiamo usare per individuare, ad esempio, il valore minimo rispetto alle variabili (numeriche) e quale indice lo assume.

apply(BodyTemperature[,-1], 2, IQR )
        Age   HeartRate Temperature 
       8.25        9.00        1.20 
apply(BodyTemperature[,-1], 2, min )
        Age   HeartRate Temperature 
       21.0        61.0        96.2 
apply(BodyTemperature[,-1], 2, which.min )
        Age   HeartRate Temperature 
         37          28           3 

Possiamo notare che la temperatura non è espressa in gradi Celsius (°C), bensì in Fahrenheit (F). Usando delle semplici operazioni di base possiamo trasformare le temperature in gradi Celsius e salvare questi nuovi dati nella colonna che chiamiamo TempC.

fromFtoC <- function(dataF) {
return((dataF - 32)*(5/9)) 
}
BodyTemperature$TempC <- (BodyTemperature$Temperature -32)*5/9
head(BodyTemperature)

2.1 Visualizzare il dataset BodyTemperature

Passiamo ora a visualizzare il contenuto del dateset.

hist(BodyTemperature$TempC )

boxplot(TempC ~ Gender, data = BodyTemperature)

Vediamo ora come affiancare diversi plot e usiamo alcuni dei parametri delle funzioni di plot.

par(mfrow=c(1, 2))
hist(BodyTemperature$TempC, main = "", xlab = "°C")
boxplot(TempC ~ Gender, data = BodyTemperature, xlab = "Gender", ylab = "Temperature [°C]")

par(mfrow=c(1, 2))
boxplot(BodyTemperature$TempC , ylab = "Temperature [°C]")
boxplot(TempC ~ Gender, data = BodyTemperature, xlab = "Gender")

par(mfrow=c(1, 2))
boxplot(BodyTemperature$TempC , ylab = "Temperature [°C]")
boxplot(TempC ~ Gender, data = BodyTemperature, xlab = "Gender")
title("Boxplots for Temperature", outer = T, line = -2)

par(mfrow=c(2, 2))
hist((BodyTemperature$TempC) , xlab = "Temperature [°C]", main = "")
boxplot(TempC ~ Gender, data = BodyTemperature, xlab = "Gender")
hist(BodyTemperature$TempC[BodyTemperature$Gender == "F"],freq = FALSE ,main = "Female", col = "pink", xlab = "Temperature [°C]")
hist(BodyTemperature$TempC[BodyTemperature$Gender == "M"],freq = FALSE ,main = "Male", col = "skyblue" , xlab = "Temperature [°C]")

Esistono altri comandi che si possono abbinare a quelli di visualizzazione visti fin qui, ad esempio abline(), text() e lines().

Esercizio Usare l’help per capire le funzionalità dei comandi precedenti e si usino per aggiungere informazioni ai grafici precedenti.

3 Pima.tr2

library(MASS) # necessario per accedere al dataset
str(Pima.tr2)
'data.frame':   300 obs. of  8 variables:
 $ npreg: int  5 7 5 0 0 5 3 1 3 2 ...
 $ glu  : int  86 195 77 165 107 97 83 193 142 128 ...
 $ bp   : int  68 70 82 76 60 76 58 50 80 78 ...
 $ skin : int  28 33 41 43 25 27 31 16 15 37 ...
 $ bmi  : num  30.2 25.1 35.8 47.9 26.4 35.6 34.3 25.9 32.4 43.3 ...
 $ ped  : num  0.364 0.163 0.156 0.259 0.133 ...
 $ age  : int  24 55 35 26 23 52 25 24 63 31 ...
 $ type : Factor w/ 2 levels "No","Yes": 1 2 1 1 1 2 1 1 1 2 ...
dim(Pima.tr2)
[1] 300   8
names(Pima.tr2)
[1] "npreg" "glu"   "bp"    "skin"  "bmi"   "ped"   "age"   "type" 
?Pima.tr2
(Pima.tr2)
summary(Pima.tr2)
     npreg             glu              bp              skin      
 Min.   : 0.000   Min.   : 56.0   Min.   : 38.00   Min.   : 7.00  
 1st Qu.: 1.000   1st Qu.:101.0   1st Qu.: 64.00   1st Qu.:21.00  
 Median : 3.000   Median :121.0   Median : 72.00   Median :29.00  
 Mean   : 3.787   Mean   :123.7   Mean   : 72.32   Mean   :29.15  
 3rd Qu.: 6.000   3rd Qu.:142.0   3rd Qu.: 80.00   3rd Qu.:36.00  
 Max.   :14.000   Max.   :199.0   Max.   :114.00   Max.   :99.00  
                                  NA's   :13       NA's   :98     
      bmi             ped              age        type    
 Min.   :18.20   Min.   :0.0780   Min.   :21.0   No :194  
 1st Qu.:27.10   1st Qu.:0.2367   1st Qu.:24.0   Yes:106  
 Median :32.00   Median :0.3360   Median :29.0            
 Mean   :32.05   Mean   :0.4357   Mean   :33.1            
 3rd Qu.:36.50   3rd Qu.:0.5867   3rd Qu.:40.0            
 Max.   :52.90   Max.   :2.2880   Max.   :72.0            
 NA's   :3                                                
which(is.na(Pima.tr2))
  [1]  804  834  836  853  854  863  871  885  887  888  889  898  899
 [14] 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113
 [27] 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126
 [40] 1127 1128 1129 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140
 [53] 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153
 [66] 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166
 [79] 1167 1168 1169 1170 1172 1173 1174 1175 1176 1177 1178 1179 1180
 [92] 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193
[105] 1194 1195 1196 1197 1198 1199 1200 1413 1430 1468
# amico di apply
lapply(lapply(Pima.tr2 , is.na),which)
$npreg
integer(0)

$glu
integer(0)

$bp
 [1] 204 234 236 253 254 263 271 285 287 288 289 298 299

$skin
 [1] 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
[17] 217 218 219 220 221 222 223 224 225 226 227 228 229 231 232 233
[33] 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249
[49] 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
[65] 266 267 268 269 270 272 273 274 275 276 277 278 279 280 281 282
[81] 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
[97] 299 300

$bmi
[1] 213 230 268

$ped
integer(0)

$age
integer(0)

$type
integer(0)

Cosa possiamo fare con i diversi NA?

Un’idea ovviamente è quella di rimuovere le righe dove compare un NA

data_no_Na <- na.omit(Pima.tr2)
dim(data_no_Na)
[1] 200   8

Abbiamo rimosso un terzo del dataset…..

Rimuovere così tanti dati impoverisce significativamente il dataset, quindi magari si possono trovare soluzioni meno radicali.

Ad esempio, possiamo notare come gran parte dei NA siano nella variabile skin… se questa non ci interessa, possiamo evitare di rimuovere quelle righe!

data_no_Na <- na.omit(Pima.tr2[,-4])
dim(data_no_Na)
[1] 284   7

In questo modo abbiamo rimosso molte meno righe. Si noti che alcune funzioni (vedi mean) hanno specifici argomenti per gestire i NA e si può quindi evitare di rimuoverli dal dataset ma gestirli i vari casi dalle funzioni

In alcuni casi, si preferisce non rimuovere alcun dato, ma piuttosto si sostituiscono i valori mancanti con informazioni prese dai dati correnti.

Ad esempio, si possono sostituire con la media, la mediana o si possono anche definire modelli più complessi per definire i sostituti.

3.1 Visualizzare il dataset Pima.tr2

Usando gli strumenti di visualizzazione visti, si esplorino le relazioni tra le variabili del dataset.

Si faccia attenzione agli NA e si confrontino i risultati nei casi in cui:

  • vengono rimosse tutte le righe contenenti NA;
  • vengono solo le righe contenenti NA nelle variabili di interesse;
  • vengono sostituiti gli NA con la media o la mediana della variabile;
LS0tCnRpdGxlOiAiVmlzdWFsaXp6YXJlIGkgZGF0aSIKYXV0aG9yOiAiRmVkZXJpY28gUmVhbGkiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogeWVzCiAgICB0b2NfZGVwdGg6IDIKICAgIHRvY19mbG9hdDoKICAgICAgY29sbGFwc2VkOiBubwogICAgICBzbW9vdGhfc2Nyb2xsOiB5ZXMKICAgIG51bWJlcl9zZWN0aW9uczogdHJ1ZQogICAgdGhlbWU6IHlldGkKICAgIGhpZ2hsaWdodDogdGFuZ28KLS0tCgohW10oLi4vLi4vbWVkaWEvQ2hvbGVyYS5qcGcpClRoaXMgW0pvaG4gU25vd10oaHR0cHM6Ly9pdC53aWtpcGVkaWEub3JnL3dpa2kvSm9obl9Tbm93XyhtZWRpY28pKSBrbmV3IHNvbWV0aGluZyEKCiMgVmlzdWFsaXp6YXJlIGkgZGF0aQoKCkxhIFZpc3VhbGl6emF6aW9uZSBkZWkgZGF0aSBwZXJtZXR0ZXJlIGRpIGluZGl2aWR1YXJlIHRyZW5kLCBjb25uZXNzaW9uaSBlIGNhcnBpcmUgaW5mb3JtYXppb25pIGRhaSBkYXRpIGNoZSBpbiBmb3JtYSB0YWJ1bGFyZSBub24gc29ubyBvdnZpZS4KCkkgY29tYW5kaSBjaGUgc2VndW9ubyBhaXV0YW5vIGluIHRhbGkgaW5kYWdpbmkuIFR1dHRpIHF1ZXN0aSBjb21hbmRpIHBvc3Nvbm8gZXNzZXJlIHBlcnNvbmFsaXp6YXRpIHVzYW5kbyBkaXZlcnNpIGNvbG9yaSwgcGVyc29uYWxpenphbmRvIGxlIGxhYmVsIGRlZ2xpIGFzc2kgbyBhZmZpYW5jYW5kbyBsZSBpbW1hZ2luaS4KClVuIHBhY2NoZXR0byBjaGUgcGVybWV0dGUgbm9uIHNvbG8gZGkgcHJvZHVycmUgaW1tYWdpbmkgKHVuIHBvJyBwacO5IGJlbGxlKSwgbWEgYW5jaGUgZGkgdXRpbGl6emFyZSBmdW56aW9uaSBhdmFuemF0ZSDDqCBbZ2dwbG90Ml0oaHR0cDovL2dncGxvdDIub3JnLykgY2hlIHBlcsOyIG5vbiB0cmF0dGVyZW1vIGluIHF1ZXN0byBjb3Jzby4KCk9sdHJlIGFkIGludHJvZHVycmUgYWxjdW5lIGZ1bnppb25pIGNoZSBwZXJtZXR0b25vIGRpIG90dGVuZXJlIGRlaSBncmFmaWNpLCBlc3Bsb3JlcmVtbyBhbmNoZSBhbGN1bmkgZGF0YXNldC4KCgojIyBCb3ggcGxvdAoKQSBsZXppb25lIGFiYmlhbW8gZ2nDoCBpbmNvbnRyYXRvIGkgYm94cGxvdC4gVXNhbmRvIFIgw6ggbW9sdG8gc2VtcGxpY2UgcHJvZHVycmUgZGVpIGJveHBsb3QgdXNhbmRvIGdsaSBhcHBvc2l0aSBjb21hbmRpLiBWZWRpYW1vbGkgdXNhbmRvIGlsIGRhdGFzZXQgKmJpcnRod3QqLCBjb250ZW51dG8gaW4gTUFTUy4KCmBgYHtyfQpsaWJyYXJ5KE1BU1MpCmJveHBsb3QoYmlydGh3dCRid3QpCmBgYAoKTGEgc2VndWVudGUgbm90YXppb25lIHBlcm1ldHRlIGRpIHByb2R1cnJlIHBpw7kgYm94cGxvdCByZWxhdGl2aSBhIGJ3dCwgZGl2aXNpIHJpc3BldHRvIGFsbGEgdmFyaWFiaWxlIHJhY2UuCgpgYGB7cn0KYm94cGxvdChid3QgfiByYWNlLCBkYXRhID0gYmlydGh3dCApCmBgYAoKX0VzZXJjaXppb18KU2kgaW52ZXN0aWdoaSBsYSB2YXJpYWJpbGUgbHd0IHJpc3BldHRvIGEgcmFjZSBlIHNtb2tlLgoKX0VzZXJjaXppb18KU2kgbm90aSBjaGUgaWwgcGFjY2hldHRvICpjYXIqIGNvbnRpZW5lIHVuYSBmdW56aW9uZSBjaGlhbWF0YSBgQm94cGxvdCgpYCwgY29uIGxhIG1haXVzY29sYSwgY2hlIHZpY2lubyBhaSBwb3RlbnppYWxpIG91dGxpZXJzIGluZGljYSBsJ2luZGljZSBkZWxsJ2VsZW1lbnRvLiBTaSByZXBsaWNoaW5vIGkgcHJlY2VkZW50aSBib3hwbG90IHVzYW5kbyB0YWxlIGZ1bnppb25lLgoKIyMgSXN0b2dyYW1taSBlIHNjYXR0ZXIgcGxvc3QKClVzYW5kbyBpbCBjb21hbmRvIGBoaXN0KClgIMOoIHBvc3NpYmlsZSBvdHRlbmVyZSBkZWdsaSBpc3RvZ3JhbW1pIHBlciBsZSB2YXJpYWJpbGksIHBhc3NhbmRvIGdsaSBhcmdvbWVudGkgaW4gbW9kbyBhbmFsb2dvIGFpIGJveHBsb3QuCgpBbmFsb2dhbWVudGUsIGNvbiBpIGNvbWFuZGkgYHBsb3QoeCx5KWAgb3IgYHNjYXR0ZXJwbG90KHgseSlgIMOoIHBvc3NpYmlsZSByYXBwcmVzZW50YXJlIGkgZGF0aSBjb21lIHB1bnRpLiAKClZlZHJlbW8gY2hlIHF1ZXN0byBwdcOyIGVzc2VyZSBtb2x0byB1dGlsZSBwZXIgb3R0ZW5lcmUgZGVnbGkgaW5kaXppIHN1bGxlIHJlbGF6aW9uaSB0cmEgdmFyaWFiaWxpLgoKIyBCb2R5VGVtcGVyYXR1cmUKCkJvZHlUZW1wZXJhdHVyZSDDqCB1biBkYXRhc2V0IGNoZSBjb250aWVuZSBsYSBtaXN1cmF6aW9uZSBkaSBxdWF0dHJvIHZhcmlhYmlsaSAoc2Vzc28sIGV0w6AsIGZyZXF1ZW56YSBjYXJkaWFjYSBlIHRlbXBlcmF0dXJhIGNvcnBvcmVhIGFzY2VsbGFyZSkgZGkgMTAwIHBhemllbnRpLiBQb3NzaWFtbyBsZWdnZXJlIGltcG9ydGFyZSBpbCBmaWxlIHNwZWNpZmljYW5kbyBpbnZlY2UgY2hlIGlsIHBlcmNvcnNvLCBsJ2luZGlyaXp6byB1cmwuCgpgYGB7ciwgZXZhbCA9IEZBTFNFfQpCb2R5VGVtcGVyYXR1cmUgPC0gcmVhZC5jc3YodXJsKCdodHRwOi8vZXh0cmFzLnNwcmluZ2VyLmNvbS8yMDEyLzk3OC0xLTQ2MTQtMTMwMS0xL0JvZHlUZW1wZXJhdHVyZS50eHQnKSkKYGBgClNlIHBlcsOyIHZpc3VhbGl6emlhbW8gcXVlc3RvIGRhdGFzZXQgcXVhbG9zYSBub24gdG9ybmEuCklsIGNvbWFuZG8gZ2l1c3RvIHNwZWNpZmljYSBxdWFsZSBzZXBhcmF0b3JlIHZpZW5lIGNvbnNpZGVyYXRvOiBpbiBxdWVzdG8gY2FzbyBsbyBzcGF6aW8sIGNoZSDDqCBkaXZlcnNvIGRhIHF1ZWxsbyBkaSBkZWZhdWx0IGRpIGByZWFkLmNzdmAuIEluIG1vZG8gYW5hbG9nbyBhdnJlbW1vIHBvdHV0byBsZWdnZXJlIGkgZGF0aSB1c2FuZG8gYHJlYWQudGFibGUoKWAuIFR1dHRhdmlhLCBjb250cm9sbGFpbW8gYmVuZSBsJ291dHB1dCBwZXIgY2FwaXJlIHNlIGVmZmV0dGl2YW1lbnRlIGlsIGRhdGFzZXQgw6ggc3RhdG8gbGV0dG8gbmVsIG1vZG8gZ2l1c3RvIQoKYGBge3J9CihCb2R5VGVtcGVyYXR1cmUgPC0gcmVhZC50YWJsZSh1cmwoJ2h0dHA6Ly9leHRyYXMuc3ByaW5nZXIuY29tLzIwMTIvOTc4LTEtNDYxNC0xMzAxLTEvQm9keVRlbXBlcmF0dXJlLnR4dCcpKSApCmBgYAoKYGBge3J9CkJvZHlUZW1wZXJhdHVyZSA8LSByZWFkLmNzdih1cmwoJ2h0dHA6Ly9leHRyYXMuc3ByaW5nZXIuY29tLzIwMTIvOTc4LTEtNDYxNC0xMzAxLTEvQm9keVRlbXBlcmF0dXJlLnR4dCcpLCBzZXA9ICIgIikKYGBgCgoKRXNwbG9yaWFtbyB1biBwbycgaWwgZGF0YXNldCBwZXIgY2FwaXJuZSBpbCBjb250ZW51dG8gZSBsZSB2YXJpYWJpbGkKYGBge3J9CnN0cihCb2R5VGVtcGVyYXR1cmUpCmhlYWQoQm9keVRlbXBlcmF0dXJlKQpuYW1lcyhCb2R5VGVtcGVyYXR1cmUpCmBgYApTb2xvIGEgc2NvcG8gZGlkYXR0aWNvLCB2ZWRpYW1vIGwnYXBwbGljYXppb25lIGRpIHVubyBkZWkgY29tYW5kaSB2aXN0aTogYGxldmVscygpYDoKYGBge3IsIGV2YWw9RkFMU0V9CkJvZHlUZW1wZXJhdHVyZSRHZW5kZXJMb25nIDwtIEJvZHlUZW1wZXJhdHVyZSRHZW5kZXIKbGV2ZWxzKEJvZHlUZW1wZXJhdHVyZSRHZW5kZXJMb25nKSA8LSBjKCJGZW1hbGUiLCAiTWFsZSIpCmBgYAoKCgpgYGB7cn0Kc3VtbWFyeShCb2R5VGVtcGVyYXR1cmUpCgpyYW5nZShCb2R5VGVtcGVyYXR1cmVbLC0xXSkKYGBgCgpJbiBtb2RvIGFuYWxvZ28gcG9zc2lhbW8gY2VyY2FyZSBkaSBvdHRlbmVyZSBsYSBkaXN0YW56YSBpbnRlcnF1YXJ0aWxlIHBlciBpIGdsaSBlbGVtZW50aSBkZWwgZGF0YXNldCBjb24gYElRUihCb2R5VGVtcGVyYXR1cmVbLDI6NF0pYCwgdHV0dGF2aWEgcXVlc3RvIGNpIGRhIGVycm9yZS4KClBvc3NpYW1vIGFnZ2lyYXJlIHF1ZXN0YSBsaW1pdGF6aW9uZSB1c2FuZG8gaWwgY29tYW5kbyBgYXBwbHkoKWAuIExhIHN0ZXNzYSBmdW56aW9uZSBsYSBwb3NzaWFtbyB1c2FyZSBwZXIgaW5kaXZpZHVhcmUsIGFkIGVzZW1waW8sIGlsICB2YWxvcmUgbWluaW1vIHJpc3BldHRvIGFsbGUgdmFyaWFiaWxpIChudW1lcmljaGUpIGUgcXVhbGUgaW5kaWNlIGxvIGFzc3VtZS4KCmBgYHtyfQphcHBseShCb2R5VGVtcGVyYXR1cmVbLC0xXSwgMiwgSVFSICkKYXBwbHkoQm9keVRlbXBlcmF0dXJlWywtMV0sIDIsIG1pbiApCmFwcGx5KEJvZHlUZW1wZXJhdHVyZVssLTFdLCAyLCB3aGljaC5taW4gKQpgYGAKClBvc3NpYW1vIG5vdGFyZSBjaGUgbGEgdGVtcGVyYXR1cmEgbm9uIMOoIGVzcHJlc3NhIGluIGdyYWRpIENlbHNpdXMgKMKwQyksIGJlbnPDrCBpbiBGYWhyZW5oZWl0IChGKS4gVXNhbmRvIGRlbGxlIHNlbXBsaWNpIG9wZXJhemlvbmkgZGkgYmFzZSBwb3NzaWFtbyB0cmFzZm9ybWFyZSBsZSB0ZW1wZXJhdHVyZSBpbiBncmFkaSBDZWxzaXVzIGUgc2FsdmFyZSBxdWVzdGkgbnVvdmkgZGF0aSBuZWxsYSBjb2xvbm5hIGNoZSBjaGlhbWlhbW8gKlRlbXBDKi4KCmBgYHtyfQpmcm9tRnRvQyA8LSBmdW5jdGlvbihkYXRhRikgewpyZXR1cm4oKGRhdGFGIC0gMzIpKig1LzkpKSAKfQpCb2R5VGVtcGVyYXR1cmUkVGVtcEMgPC0gKEJvZHlUZW1wZXJhdHVyZSRUZW1wZXJhdHVyZSAtMzIpKjUvOQoKaGVhZChCb2R5VGVtcGVyYXR1cmUpCmBgYAoKYGBge3IsIGVjaG89IEZBTFNFfQojIHJhbmdlKHggPC0gc29ydChyb3VuZChzdGF0czo6cm5vcm0oMTApIC0gMS4yLCAxKSkpCiMgaWYoYW55KHggPCAwKSkgY2F0KCJ4IGNvbnRhaW5zIG5lZ2F0aXZlIHZhbHVlc1xuIikKCiNPY2NoaW8gYWQgdXNhcmUgcm91bmQKYGBgCgojIyBWaXN1YWxpenphcmUgaWwgZGF0YXNldCBCb2R5VGVtcGVyYXR1cmUKClBhc3NpYW1vIG9yYSBhIHZpc3VhbGl6emFyZSBpbCBjb250ZW51dG8gZGVsIGRhdGVzZXQuCgpgYGB7cn0KaGlzdChCb2R5VGVtcGVyYXR1cmUkVGVtcEMgKQpgYGAKCmBgYHtyfQpib3hwbG90KFRlbXBDIH4gR2VuZGVyLCBkYXRhID0gQm9keVRlbXBlcmF0dXJlKQpgYGAKClZlZGlhbW8gb3JhIGNvbWUgYWZmaWFuY2FyZSBkaXZlcnNpIHBsb3QgZSB1c2lhbW8gYWxjdW5pIGRlaSBwYXJhbWV0cmkgZGVsbGUgZnVuemlvbmkgZGkgcGxvdC4KCmBgYHtyfQpwYXIobWZyb3c9YygxLCAyKSkKaGlzdChCb2R5VGVtcGVyYXR1cmUkVGVtcEMsIG1haW4gPSAiIiwgeGxhYiA9ICLCsEMiKQpib3hwbG90KFRlbXBDIH4gR2VuZGVyLCBkYXRhID0gQm9keVRlbXBlcmF0dXJlLCB4bGFiID0gIkdlbmRlciIsIHlsYWIgPSAiVGVtcGVyYXR1cmUgW8KwQ10iKQpgYGAKCgoKYGBge3J9CnBhcihtZnJvdz1jKDEsIDIpKQpib3hwbG90KEJvZHlUZW1wZXJhdHVyZSRUZW1wQyAsIHlsYWIgPSAiVGVtcGVyYXR1cmUgW8KwQ10iKQpib3hwbG90KFRlbXBDIH4gR2VuZGVyLCBkYXRhID0gQm9keVRlbXBlcmF0dXJlLCB4bGFiID0gIkdlbmRlciIpCmBgYAoKYGBge3J9CnBhcihtZnJvdz1jKDEsIDIpKQpib3hwbG90KEJvZHlUZW1wZXJhdHVyZSRUZW1wQyAsIHlsYWIgPSAiVGVtcGVyYXR1cmUgW8KwQ10iKQpib3hwbG90KFRlbXBDIH4gR2VuZGVyLCBkYXRhID0gQm9keVRlbXBlcmF0dXJlLCB4bGFiID0gIkdlbmRlciIpCnRpdGxlKCJCb3hwbG90cyBmb3IgVGVtcGVyYXR1cmUiLCBvdXRlciA9IFQsIGxpbmUgPSAtMikKYGBgCgoKYGBge3IsIGZpZy53aWR0aCA9IDEwLCBmaWcuaGVpZ2h0ID0gMTB9CnBhcihtZnJvdz1jKDIsIDIpKQpoaXN0KChCb2R5VGVtcGVyYXR1cmUkVGVtcEMpICwgeGxhYiA9ICJUZW1wZXJhdHVyZSBbwrBDXSIsIG1haW4gPSAiIikKYm94cGxvdChUZW1wQyB+IEdlbmRlciwgZGF0YSA9IEJvZHlUZW1wZXJhdHVyZSwgeGxhYiA9ICJHZW5kZXIiKQpoaXN0KEJvZHlUZW1wZXJhdHVyZSRUZW1wQ1tCb2R5VGVtcGVyYXR1cmUkR2VuZGVyID09ICJGIl0sZnJlcSA9IEZBTFNFICxtYWluID0gIkZlbWFsZSIsIGNvbCA9ICJwaW5rIiwgeGxhYiA9ICJUZW1wZXJhdHVyZSBbwrBDXSIpCmhpc3QoQm9keVRlbXBlcmF0dXJlJFRlbXBDW0JvZHlUZW1wZXJhdHVyZSRHZW5kZXIgPT0gIk0iXSxmcmVxID0gRkFMU0UgLG1haW4gPSAiTWFsZSIsIGNvbCA9ICJza3libHVlIiAsIHhsYWIgPSAiVGVtcGVyYXR1cmUgW8KwQ10iKQpgYGAKCkVzaXN0b25vIGFsdHJpIGNvbWFuZGkgY2hlIHNpIHBvc3Nvbm8gYWJiaW5hcmUgYSBxdWVsbGkgZGkgdmlzdWFsaXp6YXppb25lIHZpc3RpIGZpbiBxdWksIGFkIGVzZW1waW8gYGFibGluZSgpYCwgYHRleHQoKWAgZSBgbGluZXMoKWAuIAoKX0VzZXJjaXppb18gVXNhcmUgbCdoZWxwIHBlciBjYXBpcmUgbGUgZnVuemlvbmFsaXTDoCBkZWkgY29tYW5kaSBwcmVjZWRlbnRpIGUgc2kgdXNpbm8gcGVyIGFnZ2l1bmdlcmUgaW5mb3JtYXppb25pIGFpIGdyYWZpY2kgcHJlY2VkZW50aS4KCiMgUGltYS50cjIgCgoKCmBgYHtyfQpsaWJyYXJ5KE1BU1MpICMgbmVjZXNzYXJpbyBwZXIgYWNjZWRlcmUgYWwgZGF0YXNldApzdHIoUGltYS50cjIpCmRpbShQaW1hLnRyMikKbmFtZXMoUGltYS50cjIpCj9QaW1hLnRyMgooUGltYS50cjIpCnN1bW1hcnkoUGltYS50cjIpCmBgYAoKYGBge3J9CndoaWNoKGlzLm5hKFBpbWEudHIyKSkKIyBhbWljbyBkaSBhcHBseQpsYXBwbHkobGFwcGx5KFBpbWEudHIyICwgaXMubmEpLHdoaWNoKQpgYGAKCkNvc2EgcG9zc2lhbW8gZmFyZSBjb24gaSBkaXZlcnNpIE5BPwoKVW4naWRlYSBvdnZpYW1lbnRlIMOoIHF1ZWxsYSBkaSByaW11b3ZlcmUgbGUgcmlnaGUgZG92ZSBjb21wYXJlIHVuIE5BCgpgYGB7cn0KZGF0YV9ub19OYSA8LSBuYS5vbWl0KFBpbWEudHIyKQoKZGltKGRhdGFfbm9fTmEpCmBgYAoKQWJiaWFtbyByaW1vc3NvIHVuIHRlcnpvIGRlbCBkYXRhc2V0Li4uLi4KClJpbXVvdmVyZSBjb3PDrCB0YW50aSBkYXRpIGltcG92ZXJpc2NlIHNpZ25pZmljYXRpdmFtZW50ZSBpbCBkYXRhc2V0LCBxdWluZGkgbWFnYXJpIHNpIHBvc3Nvbm8gdHJvdmFyZSBzb2x1emlvbmkgbWVubyByYWRpY2FsaS4KCkFkIGVzZW1waW8sIHBvc3NpYW1vIG5vdGFyZSBjb21lIGdyYW4gcGFydGUgZGVpIE5BIHNpYW5vIG5lbGxhIHZhcmlhYmlsZSAqc2tpbiouLi4gc2UgcXVlc3RhIG5vbiBjaSBpbnRlcmVzc2EsIHBvc3NpYW1vIGV2aXRhcmUgZGkgcmltdW92ZXJlIHF1ZWxsZSByaWdoZSEKCmBgYHtyfQpkYXRhX25vX05hIDwtIG5hLm9taXQoUGltYS50cjJbLC00XSkKCmRpbShkYXRhX25vX05hKQpgYGAKCkluIHF1ZXN0byBtb2RvIGFiYmlhbW8gcmltb3NzbyBtb2x0ZSBtZW5vIHJpZ2hlLgpTaSBub3RpIGNoZSBhbGN1bmUgZnVuemlvbmkgKHZlZGkgbWVhbikgaGFubm8gc3BlY2lmaWNpIGFyZ29tZW50aSBwZXIgZ2VzdGlyZSBpIE5BIGUgc2kgcHXDsiBxdWluZGkgZXZpdGFyZSBkaSByaW11b3ZlcmxpIGRhbCBkYXRhc2V0IG1hIGdlc3RpcmxpIGkgdmFyaSBjYXNpIGRhbGxlIGZ1bnppb25pCgoKSW4gYWxjdW5pIGNhc2ksIHNpIHByZWZlcmlzY2Ugbm9uIHJpbXVvdmVyZSBhbGN1biBkYXRvLCBtYSBwaXV0dG9zdG8gc2kgc29zdGl0dWlzY29ubyBpIHZhbG9yaSBtYW5jYW50aSBjb24gaW5mb3JtYXppb25pIHByZXNlIGRhaSBkYXRpIGNvcnJlbnRpLgoKQWQgZXNlbXBpbywgc2kgcG9zc29ubyBzb3N0aXR1aXJlIGNvbiBsYSBtZWRpYSwgbGEgbWVkaWFuYSBvIHNpIHBvc3Nvbm8gYW5jaGUgZGVmaW5pcmUgbW9kZWxsaSBwacO5IGNvbXBsZXNzaSBwZXIgZGVmaW5pcmUgaSBzb3N0aXR1dGkuCgojIyBWaXN1YWxpenphcmUgaWwgZGF0YXNldCBQaW1hLnRyMgoKVXNhbmRvIGdsaSBzdHJ1bWVudGkgZGkgdmlzdWFsaXp6YXppb25lIHZpc3RpLCBzaSBlc3Bsb3Jpbm8gbGUgcmVsYXppb25pIHRyYSBsZSB2YXJpYWJpbGkgZGVsIGRhdGFzZXQuCgpTaSBmYWNjaWEgYXR0ZW56aW9uZSBhZ2xpIE5BIGUgc2kgY29uZnJvbnRpbm8gaSByaXN1bHRhdGkgbmVpIGNhc2kgaW4gY3VpOgoKLSB2ZW5nb25vIHJpbW9zc2UgdHV0dGUgbGUgcmlnaGUgY29udGVuZW50aSBOQTsKLSB2ZW5nb25vIHNvbG8gbGUgcmlnaGUgY29udGVuZW50aSBOQSBuZWxsZSB2YXJpYWJpbGkgZGkgaW50ZXJlc3NlOwotIHZlbmdvbm8gc29zdGl0dWl0aSBnbGkgTkEgY29uIGxhIG1lZGlhIG8gbGEgbWVkaWFuYSBkZWxsYSB2YXJpYWJpbGU7CgpgYGB7ciwgaW5jbHVkZT1GQUxTRX0KcGxvdChnbHUgfiBibWksIGRhdGEgPSBQaW1hLnRyMikKcGxvdCggUGltYS50cjIkYm1pLCBQaW1hLnRyMiRnbHUpCnBsb3QoYnAgfiBibWksIGRhdGEgPSBQaW1hLnRyMikKcGxvdChucHJlZyB+IGJtaSwgZGF0YSA9IFBpbWEudHIyKQpwbG90KG5wcmVnIH4gYWdlLCBkYXRhID0gUGltYS50cjIpCmJveHBsb3QobnByZWcgfiB0eXBlLCBkYXRhID0gUGltYS50cjIpCmhpc3QoUGltYS50cjIkZ2x1KQojT2NjaGlvIGcgbyBHIHN1IGdsdWNvc2UgKGNvbmNlbnRyYXRpb24pIC0+IGltcG9ydGFudGUgc2FwZXJlIGxlIHVuaXTDoCBkaSBtaXN1cmEhCmBgYAoKCgpgYGB7ciwgaW5jbHVkZT1GQUxTRX0KcGFyKG1mcm93PWMoMSwgMikpCmhpc3QoUGltYS50cjJbIFBpbWEudHIyJHR5cGUgPT0gIlllcyIgICwiZ2x1Il0sIG1haW4gPSAiRGlhYmV0ZSIsICAgIHhsYWIgPSAiR2x1Y29zZSIpCmhpc3QoUGltYS50cjJbIFBpbWEudHIyJHR5cGUgPT0gIk5vIiAgICwiZ2x1Il0sIG1haW4gPSAiTm8gRGlhYmV0ZSIsIHhsYWIgPSAiR2x1Y29zZSIpCgpgYGAKCgoKCgoK