Preamble

Standard packages

### Load packages
library(tidyverse) # Collection of all the good stuff like dplyr, ggplot2 ect.
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.4     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   2.0.1     ✓ forcats 0.5.1
── Conflicts ─────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
library(magrittr) # For extra-piping operators (eg. %<>%)

Attaching package: ‘magrittr’

The following object is masked from ‘package:purrr’:

    set_names

The following object is masked from ‘package:tidyr’:

    extract

Load data

data <- read_csv('https://sds-aau.github.io/SDS-master/M1/data/cities.csv')
Rows: 780 Columns: 25
── Column specification ─────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): place, alpha-2, region, sub-region
dbl (21): cost_nomad, cost_coworking, cost_expat, coffee_in_cafe, cost_beer, places_to_work, free_wifi_av...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
data %>% glimpse()
Rows: 780
Columns: 25
$ place                  <chr> "Budapest", "Chiang Mai", "Phuket", "Bangkok", "Ko Samui", "Ko Lanta", "Chia…
$ `alpha-2`              <chr> "HU", "TH", "TH", "TH", "TH", "TH", "TH", "TH", "TH", "TH", "TH", "TH", "TH"…
$ region                 <chr> "Europe", "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "A…
$ `sub-region`           <chr> "Eastern Europe", "South-eastern Asia", "South-eastern Asia", "South-eastern…
$ cost_nomad             <dbl> 1364, 777, 1012, 1197, 1352, 812, 1134, 1134, 866, 1232, 777, 1247, 1285, 12…
$ cost_coworking         <dbl> 152.41, 98.88, 155.43, 131.41, 169.56, 135.65, 195.38, 113.04, 172.39, 172.3…
$ cost_expat             <dbl> 1273, 780, 1714, 1158, 1347, 1016, 1119, 1100, 1483, 2173, 1152, 1092, 730, …
$ coffee_in_cafe         <dbl> 1.73, 0.85, 1.41, 2.12, 1.41, 1.41, 1.41, 2.12, 1.84, 1.55, 1.41, 1.41, 1.41…
$ cost_beer              <dbl> 1.73, 0.85, 1.41, 2.12, 1.41, 1.41, 1.41, 2.12, 1.84, 1.55, 1.41, 1.41, 1.41…
$ places_to_work         <dbl> 1.0, 0.8, 0.8, 1.0, 0.8, 0.4, 0.6, 0.4, 0.4, 0.4, 0.4, 0.2, 0.2, 0.6, 1.0, 0…
$ free_wifi_available    <dbl> 0.40, 0.60, 0.40, 1.00, 0.40, 0.20, 0.60, 0.60, 0.40, 0.40, 0.20, 0.40, 0.20…
$ internet_speed         <dbl> 31, 14, 14, 24, 15, 15, 12, 9, 0, 13, 2, 2, 3, 4, 15, 26, 16, 25, 7, 9, 5, 1…
$ freedom_score          <dbl> 0.6, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.8, 0…
$ peace_score            <dbl> 8.000000e-01, 4.000000e-01, 4.000000e-01, 4.000000e-01, 4.000000e-01, 4.0000…
$ safety                 <dbl> 0.60, 0.80, 0.80, 0.77, 0.80, 0.80, 0.60, 0.56, 0.60, 0.40, 0.60, 0.60, 0.80…
$ fragile_states_index   <dbl> 5.270000e+01, 7.880000e+01, 7.880000e+01, 7.880000e+01, 7.880000e+01, 7.8800…
$ press_freedom_index    <dbl> 28.17, 44.53, 44.53, 44.53, 44.53, 44.53, 44.53, 44.53, 44.53, 44.53, 44.53,…
$ female_friendly        <dbl> 1.0, 0.8, 0.6, 0.8, 0.8, 1.0, 0.6, 0.2, 0.6, 0.6, 0.6, 1.0, 0.6, 0.6, 1.0, 0…
$ lgbt_friendly          <dbl> 0.27, 0.60, 0.80, 0.80, 0.80, 0.80, 0.40, 1.00, 0.40, 0.40, 0.80, 0.80, 0.60…
$ friendly_to_foreigners <dbl> 0.60, 0.60, 0.60, 1.00, 1.00, 0.80, 1.00, 1.00, 0.80, 0.60, 0.80, 0.80, 0.80…
$ racism                 <dbl> 0.40, 0.40, 0.42, 0.42, 0.40, 0.40, 0.60, 0.40, 0.40, 0.60, 0.40, 0.40, 0.40…
$ leisure                <dbl> 0.80, 0.62, 0.60, 0.82, 0.80, 0.62, 0.80, 0.80, 0.60, 0.60, 0.60, 0.60, 0.60…
$ life_score             <dbl> 0.86, 0.75, 0.75, 0.72, 0.80, 0.73, 0.76, 0.66, 0.70, 0.66, 0.65, 0.67, 0.66…
$ nightlife              <dbl> 1.00, 0.40, 0.82, 1.00, 0.80, 0.43, 0.80, 0.80, 0.40, 0.60, 1.00, 0.40, 1.00…
$ weed                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0…

EDA

# Variables for descriptives
vars.desc <- c("cost_nomad", "places_to_work", "freedom_score", "friendly_to_foreigners", "life_score")
Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

First, lets look at a classical correlation matrix.

ggcorr(data[,vars.desc], label = TRUE, label_size = 3, label_round = 2, label_alpha = TRUE)

Dimionality Reduction

library(FactoMineR)
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
library(factoextra)
Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

Do the PCA

res_pca <- data %>%
  column_to_rownames('place') %>%
  select_if(is_numeric) %>%
  PCA(scale.unit = TRUE, graph =FALSE)

Do a screeplot

res_pca %>% 
  fviz_screeplot()
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Plot the variable loadings

res_pca %>%
  fviz_pca_var(alpha.var = "cos2",
               col.var = "contrib",
               gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
               repel = TRUE) 

Bonus: Plot variables plus observations

res_pca %>% glimpse()
List of 5
 $ eig : num [1:21, 1:3] 8.29 1.98 1.28 1.12 1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:21] "comp 1" "comp 2" "comp 3" "comp 4" ...
  .. ..$ : chr [1:3] "eigenvalue" "percentage of variance" "cumulative percentage of variance"
 $ var :List of 4
  ..$ coord  : num [1:21, 1:5] 0.659 0.431 0.6 0.751 0.751 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ cor    : num [1:21, 1:5] 0.659 0.431 0.6 0.751 0.751 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ cos2   : num [1:21, 1:5] 0.434 0.186 0.36 0.564 0.564 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ contrib: num [1:21, 1:5] 5.24 2.24 4.34 6.81 6.81 ...
  .. ..- attr(*, "dimnames")=List of 2
 $ ind :List of 4
  ..$ coord  : num [1:780, 1:5] 0.549 -1.894 -1.332 0.288 -0.504 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ cos2   : num [1:780, 1:5] 0.01766 0.27572 0.15585 0.00399 0.01865 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ contrib: num [1:780, 1:5] 0.00467 0.05545 0.02742 0.00128 0.00392 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ dist   : Named num [1:780] 4.13 3.61 3.37 4.55 3.69 ...
  .. ..- attr(*, "names")= chr [1:780] "Budapest" "Chiang Mai" "Phuket" "Bangkok" ...
 $ svd :List of 3
  ..$ vs: num [1:21] 2.88 1.41 1.13 1.06 1 ...
  ..$ U : num [1:780, 1:5] 0.1908 -0.6577 -0.4625 0.0999 -0.1749 ...
  ..$ V : num [1:21, 1:5] 0.229 0.15 0.208 0.261 0.261 ...
 $ call:List of 9
  ..$ row.w     : num [1:780] 0.00128 0.00128 0.00128 0.00128 0.00128 ...
  ..$ col.w     : num [1:21] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ scale.unit: logi TRUE
  ..$ ncp       : num 5
  ..$ centre    : num [1:21] 2331.7 210.3 1880.9 3.3 3.3 ...
  ..$ ecart.type: num [1:21] 1117.41 173.96 1265.34 1.98 1.98 ...
  ..$ X         :'data.frame':  780 obs. of  21 variables:
  .. ..$ cost_nomad            : num [1:780] 1364 777 1012 1197 1352 ...
  .. ..$ cost_coworking        : num [1:780] 152.4 98.9 155.4 131.4 169.6 ...
  .. ..$ cost_expat            : num [1:780] 1273 780 1714 1158 1347 ...
  .. ..$ coffee_in_cafe        : num [1:780] 1.73 0.85 1.41 2.12 1.41 1.41 1.41 2.12 1.84 1.55 ...
  .. ..$ cost_beer             : num [1:780] 1.73 0.85 1.41 2.12 1.41 1.41 1.41 2.12 1.84 1.55 ...
  .. ..$ places_to_work        : num [1:780] 1 0.8 0.8 1 0.8 0.4 0.6 0.4 0.4 0.4 ...
  .. ..$ free_wifi_available   : num [1:780] 0.4 0.6 0.4 1 0.4 0.2 0.6 0.6 0.4 0.4 ...
  .. ..$ internet_speed        : num [1:780] 31 14 14 24 15 15 12 9 0 13 ...
  .. ..$ freedom_score         : num [1:780] 0.6 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...
  .. ..$ peace_score           : num [1:780] 0.8 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 ...
  .. ..$ safety                : num [1:780] 0.6 0.8 0.8 0.77 0.8 0.8 0.6 0.56 0.6 0.4 ...
  .. ..$ fragile_states_index  : num [1:780] 52.7 78.8 78.8 78.8 78.8 78.8 78.8 78.8 78.8 78.8 ...
  .. ..$ press_freedom_index   : num [1:780] 28.2 44.5 44.5 44.5 44.5 ...
  .. ..$ female_friendly       : num [1:780] 1 0.8 0.6 0.8 0.8 1 0.6 0.2 0.6 0.6 ...
  .. ..$ lgbt_friendly         : num [1:780] 0.27 0.6 0.8 0.8 0.8 0.8 0.4 1 0.4 0.4 ...
  .. ..$ friendly_to_foreigners: num [1:780] 0.6 0.6 0.6 1 1 0.8 1 1 0.8 0.6 ...
  .. ..$ racism                : num [1:780] 0.4 0.4 0.42 0.42 0.4 0.4 0.6 0.4 0.4 0.6 ...
  .. ..$ leisure               : num [1:780] 0.8 0.62 0.6 0.82 0.8 0.62 0.8 0.8 0.6 0.6 ...
  .. ..$ life_score            : num [1:780] 0.86 0.75 0.75 0.72 0.8 0.73 0.76 0.66 0.7 0.66 ...
  .. ..$ nightlife             : num [1:780] 1 0.4 0.82 1 0.8 0.43 0.8 0.8 0.4 0.6 ...
  .. ..$ weed                  : num [1:780] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ row.w.init: num [1:780] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ call      : language PCA(X = ., scale.unit = TRUE, graph = FALSE)
 - attr(*, "class")= chr [1:2] "PCA" "list"
res_pca %>%
  fviz_pca_biplot() 

alpha.ind = “cos2”, col.ind = “contrib”, gradient.cols = c(“#00AFBB”, “#E7B800”, “#FC4E07”), geom = “point”

Clustering

Determine number of clusters

data %>% 
  drop_na() %>%
  column_to_rownames('place') %>%
  select_if(is_numeric) %>%
  scale() %>%
  fviz_nbclust(hcut, method = "wss")  

Do the clustering

hc <- data %>%
  column_to_rownames('place') %>%
  select_if(is_numeric) %>%
  hcut(hc_func = "hclust", 
       k = 3, 
       stand = TRUE)

Visualize clusters

hc %>% 
  fviz_cluster(data = data %>% select_if(is_numeric))  

Where do we find the clusters?

hc %>%
  glimpse()
List of 12
 $ merge      : int [1:779, 1:2] -460 -493 -547 -463 -459 -487 -475 -578 -453 -539 ...
 $ height     : num [1:779] 0.157 0.161 0.169 0.201 0.22 ...
 $ order      : int [1:780] 111 122 28 95 128 140 141 145 146 73 ...
 $ labels     : chr [1:780] "Budapest" "Chiang Mai" "Phuket" "Bangkok" ...
 $ method     : chr "ward.D2"
 $ call       : language stats::hclust(d = x, method = hc_method)
 $ dist.method: chr "euclidean"
 $ cluster    : Named int [1:780] 1 1 1 1 1 1 1 1 2 2 ...
  ..- attr(*, "names")= chr [1:780] "Budapest" "Chiang Mai" "Phuket" "Bangkok" ...
 $ nbclust    : num 3
 $ silinfo    :List of 3
  ..$ widths         :'data.frame': 780 obs. of  3 variables:
  .. ..$ cluster  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ neighbor : num [1:780] 2 3 2 3 2 2 2 2 2 2 ...
  .. ..$ sil_width: num [1:780] 0.228 0.216 0.215 0.213 0.212 ...
  ..$ clus.avg.widths: num [1:3] 0.0542 0.2006 0.1811
  ..$ avg.width      : num 0.144
 $ size       : int [1:3] 257 200 323
 $ data       : num [1:780, 1:21] -0.866 -1.39 -1.18 -1.015 -0.876 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:780] "Budapest" "Chiang Mai" "Phuket" "Bangkok" ...
  .. ..$ : chr [1:21] "cost_nomad" "cost_coworking" "cost_expat" "coffee_in_cafe" ...
  ..- attr(*, "scaled:center")= Named num [1:21] 2331.7 210.3 1880.9 3.3 3.3 ...
  .. ..- attr(*, "names")= chr [1:21] "cost_nomad" "cost_coworking" "cost_expat" "coffee_in_cafe" ...
  ..- attr(*, "scaled:scale")= Named num [1:21] 1118.12 174.07 1266.15 1.98 1.98 ...
  .. ..- attr(*, "names")= chr [1:21] "cost_nomad" "cost_coworking" "cost_expat" "coffee_in_cafe" ...
 - attr(*, "class")= chr [1:2] "hclust" "hcut"
hc$cluster %>%
  head()
  Budapest Chiang Mai     Phuket    Bangkok   Ko Samui   Ko Lanta 
         1          1          1          1          1          1 

Add them to dataset

data[,"cluster"] <- hc$cluster

Inspect clusters per region

table(data$cluster, data$region)
   
    Africa Americas Asia Europe Oceania
  1     17       67   90     83       0
  2     23       15  145     16       1
  3      5      139   27    136      16

Also add PCA to orignal data

data[,"pca1"] <- res_pca$ind$coord[,1]
data[,"pca2"] <- res_pca$ind$coord[,2]

Component mean per cluster

data %>%
  group_by(cluster) %>%
  summarise(pca1 = pca1 %>% mean(),
            pca2 = pca2 %>% mean())

Bonus: add trips data

trips <- read_csv('https://sds-aau.github.io/SDS-master/M1/data/trips.csv')
New names:
* `` -> ...1
Rows: 46510 Columns: 11
── Column specification ─────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (6): username, country, country_code, country_slug, place, place_slug
dbl  (3): ...1, latitude, longitude
date (2): date_end, date_start

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
trips %>% glimpse()
Rows: 46,510
Columns: 11
$ ...1         <dbl> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, …
$ username     <chr> "@lewellenmichael", "@lewellenmichael", "@lewellenmichael", "@lewellenmichael", "@wayl…
$ country      <chr> "Mexico", "Mexico", "Mexico", "Jordan", "China", "Vietnam", "Hong Kong", "China", "Chi…
$ country_code <chr> "MX", "MX", "MX", "JO", "CN", "VN", "HK", "CN", "CN", "CN", "TH", "MY", "KH", "VN", "I…
$ country_slug <chr> "mexico", "mexico", "mexico", "jordan", "china", "vietnam", "hong-kong", "china", "chi…
$ date_end     <date> 2018-06-15, 2018-06-03, 2017-11-05, 2017-08-07, 2017-03-18, 2017-02-16, 2016-09-01, 2…
$ date_start   <date> 2018-06-04, 2018-05-31, 2017-11-01, 2017-07-24, 2017-02-17, 2016-09-02, 2016-08-02, 2…
$ latitude     <dbl> 21, 19, 21, 31, 40, 10, 22, 22, 22, 18, 7, 3, 11, 10, 13, 26, 27, 27, 28, 28, 19, 11, …
$ longitude    <dbl> -101, -99, -86, 35, 122, 106, 114, 114, 113, 109, 98, 101, 104, 106, 80, 75, 78, 78, 7…
$ place        <chr> "Guanajuato", "Mexico City", "Cancun", "Amman", "Yingkou", "Ho Chi Minh City", "Shenzh…
$ place_slug   <chr> "mexico", "mexico-city-mexico", "cancun-mexico", "amman-jordan", "china", "ho-chi-minh…

Add number of trips per city

data %<>%
  left_join(trips %>% count(place, sort = TRUE, name = 'n_city'), by = 'place')

Check most popular cities per cluster

data %>%
  select(place, cluster, n_city) %>%
  group_by(cluster) %>%
  arrange(desc(n_city)) %>%
  slice(1:5) %>%
  ungroup() 

Count cluster popularity

data %>%
  count(cluster, wt = n_city)

To finish up, lets plot it in a map, simplest way possible.

geo_merge <- trips %>%
  select(place, longitude, latitude) %>%
  distinct(place, .keep_all = TRUE)
data %<>%
  left_join(geo_merge , by = 'place')

Load a worldmap geom

library(ggmap)
Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
Please cite ggmap if you use it! See citation("ggmap") for details.

Attaching package: ‘ggmap’

The following object is masked from ‘package:magrittr’:

    inset
mapWorld <- borders("world", colour = "gray50", fill = "gray50")

GThat’s how it looks

mapWorld
mapping: group = ~group, x = ~long, y = ~lat 
geom_polygon: na.rm = FALSE, rule = evenodd
stat_identity: na.rm = FALSE
position_identity 

Add it to an empty ggplot surface

mp <- ggplot() +   
  mapWorld 

That’s how it looks so far

mp

Add a geom with the cities as points

nomad_map <- mp + 
  geom_point(data = data, aes(x = longitude, y = latitude, col = factor(cluster))) 
nomad_map

Or do a density plot of popular nomad cities

mp + 
  stat_density2d(data = trips, 
                 aes(x = longitude, y = latitude, fill = stat(nlevel), col = stat(nlevel) ), 
                 alpha = 0.2, size = 0.2, bins = 10, geom = "polygon") +
  scale_fill_gradient(low = "skyblue", high = "red") +
  scale_color_gradient(low = "skyblue", high = "red")

LS0tCnRpdGxlOiAiSW4tQ2xhc3MgRXhlcmNpc2U6IFVNTCAmIERpZ2l0YWwgTm9tYWRzKSIKYXV0aG9yOiAiRGFuaWVsIFMuIEhhaW4gKGRzaEBidXNpbmVzcy5hYXUuZGspIgpkYXRlOiAiVXBkYXRlZCBgciBmb3JtYXQoU3lzLnRpbWUoKSwgJyVCICVkLCAlWScpYCIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICBjb2RlX2ZvbGRpbmc6IHNob3cKICAgIGRmX3ByaW50OiBwYWdlZAogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDIKICAgIHRvY19mbG9hdDoKICAgICAgY29sbGFwc2VkOiBmYWxzZQogICAgdGhlbWU6IGZsYXRseQotLS0KCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQojIEtuaXRyIG9wdGlvbnMKIyMjIEdlbmVyaWMgcHJlYW1ibGUKcm0obGlzdD1scygpKTsgZ3JhcGhpY3Mub2ZmKCkgClN5cy5zZXRlbnYoTEFORyA9ICJlbiIpICMgRm9yIGVuZ2xpc2ggbGFuZ3VhZ2UKb3B0aW9ucyhzY2lwZW4gPSA1KSAjIFRvIGRlYWN0aXZhdGUgYW5ub3lpbmcgc2NpZW50aWZpYyBudW1iZXIgbm90YXRpb24KCiMgcm0obGlzdD1scygpKTsgZ3JhcGhpY3Mub2ZmKCkgIyBnZXQgcmlkIG9mIGV2ZXJ5dGhpbmcgaW4gdGhlIHdvcmtzcGFjZQppZiAoIXJlcXVpcmUoImtuaXRyIikpIGluc3RhbGwucGFja2FnZXMoImtuaXRyIik7IGxpYnJhcnkoa25pdHIpICMgRm9yIGRpc3BsYXkgb2YgdGhlIG1hcmtkb3duCgojIyMgS25pdHIgb3B0aW9ucwprbml0cjo6b3B0c19jaHVuayRzZXQod2FybmluZz1GQUxTRSwKICAgICAgICAgICAgICAgICAgICAgbWVzc2FnZT1GQUxTRSwKICAgICAgICAgICAgICAgICAgICAgZmlnLmFsaWduPSJjZW50ZXIiCiAgICAgICAgICAgICAgICAgICAgICkKYGBgCgojIFByZWFtYmxlCgojIyBTdGFuZGFyZCBwYWNrYWdlcwoKYGBge3J9CiMjIyBMb2FkIHBhY2thZ2VzCmxpYnJhcnkodGlkeXZlcnNlKSAjIENvbGxlY3Rpb24gb2YgYWxsIHRoZSBnb29kIHN0dWZmIGxpa2UgZHBseXIsIGdncGxvdDIgZWN0LgpsaWJyYXJ5KG1hZ3JpdHRyKSAjIEZvciBleHRyYS1waXBpbmcgb3BlcmF0b3JzIChlZy4gJTw+JSkKYGBgCgojIyBMb2FkIGRhdGEKCmBgYHtyfQpkYXRhIDwtIHJlYWRfY3N2KCdodHRwczovL3Nkcy1hYXUuZ2l0aHViLmlvL1NEUy1tYXN0ZXIvTTEvZGF0YS9jaXRpZXMuY3N2JykKZGF0YSAlPiUgZ2xpbXBzZSgpCmBgYAoKIyBFREEKCmBgYHtyfQojIFZhcmlhYmxlcyBmb3IgZGVzY3JpcHRpdmVzCnZhcnMuZGVzYyA8LSBjKCJjb3N0X25vbWFkIiwgInBsYWNlc190b193b3JrIiwgImZyZWVkb21fc2NvcmUiLCAiZnJpZW5kbHlfdG9fZm9yZWlnbmVycyIsICJsaWZlX3Njb3JlIikKYGBgCgpgYGB7cix3YXJuaW5nPUZBTFNFLGVjaG89RkFMU0V9CmxpYnJhcnkoR0dhbGx5KQpgYGAKCkZpcnN0LCBsZXRzIGxvb2sgYXQgYSBjbGFzc2ljYWwgY29ycmVsYXRpb24gbWF0cml4LgoKYGBge3J9CmdnY29ycihkYXRhWyx2YXJzLmRlc2NdLCBsYWJlbCA9IFRSVUUsIGxhYmVsX3NpemUgPSAzLCBsYWJlbF9yb3VuZCA9IDIsIGxhYmVsX2FscGhhID0gVFJVRSkKYGBgCgojIERpbWlvbmFsaXR5IFJlZHVjdGlvbgoKYGBge3J9CmxpYnJhcnkoRmFjdG9NaW5lUikKbGlicmFyeShmYWN0b2V4dHJhKQpgYGAKCkRvIHRoZSBQQ0EKCmBgYHtyfQpyZXNfcGNhIDwtIGRhdGEgJT4lCiAgY29sdW1uX3RvX3Jvd25hbWVzKCdwbGFjZScpICU+JQogIHNlbGVjdF9pZihpc19udW1lcmljKSAlPiUKICBQQ0Eoc2NhbGUudW5pdCA9IFRSVUUsIGdyYXBoID1GQUxTRSkKYGBgCgpEbyBhIHNjcmVlcGxvdAoKYGBge3IsZmlnLmFsaWduPSdjZW50ZXInfQpyZXNfcGNhICU+JSAKICBmdml6X3NjcmVlcGxvdCgpCmBgYAoKUGxvdCB0aGUgdmFyaWFibGUgbG9hZGluZ3MKCmBgYHtyLGZpZy53aWR0aD01LGZpZy5oZWlnaHQ9NSxmaWcuYWxpZ249J2NlbnRlcid9CnJlc19wY2EgJT4lCiAgZnZpel9wY2FfdmFyKGFscGhhLnZhciA9ICJjb3MyIiwKICAgICAgICAgICAgICAgY29sLnZhciA9ICJjb250cmliIiwKICAgICAgICAgICAgICAgZ3JhZGllbnQuY29scyA9IGMoIiMwMEFGQkIiLCAiI0U3QjgwMCIsICIjRkM0RTA3IiksCiAgICAgICAgICAgICAgIHJlcGVsID0gVFJVRSkgCmBgYAoKCgpCb251czogUGxvdCB2YXJpYWJsZXMgcGx1cyBvYnNlcnZhdGlvbnMKCmBgYHtyfQpyZXNfcGNhICU+JSBnbGltcHNlKCkKYGBgCgoKYGBge3IsLGZpZy53aWR0aD0xNSxmaWcuaGVpZ2h0PTEwLGZpZy5hbGlnbj0nY2VudGVyJ30KcmVzX3BjYSAlPiUKICBmdml6X3BjYV9iaXBsb3QoKSAKYGBgCgphbHBoYS5pbmQgPSAiY29zMiIsCiAgICAgICAgICAgICAgICAgIGNvbC5pbmQgPSAiY29udHJpYiIsCiAgICAgICAgICAgICAgICAgIGdyYWRpZW50LmNvbHMgPSBjKCIjMDBBRkJCIiwgIiNFN0I4MDAiLCAiI0ZDNEUwNyIpLAogICAgICAgICAgICAgICAgICBnZW9tID0gInBvaW50IgogICAgICAgICAgICAgICAgICAKCiMgQ2x1c3RlcmluZwoKCkRldGVybWluZSBudW1iZXIgb2YgY2x1c3RlcnMKCmBgYHtyLGZpZy5hbGlnbj0nY2VudGVyJ30KZGF0YSAlPiUgCiAgZHJvcF9uYSgpICU+JQogIGNvbHVtbl90b19yb3duYW1lcygncGxhY2UnKSAlPiUKICBzZWxlY3RfaWYoaXNfbnVtZXJpYykgJT4lCiAgc2NhbGUoKSAlPiUKICBmdml6X25iY2x1c3QoaGN1dCwgbWV0aG9kID0gIndzcyIpICAKYGBgCgpEbyB0aGUgY2x1c3RlcmluZwoKYGBge3J9CmhjIDwtIGRhdGEgJT4lCiAgY29sdW1uX3RvX3Jvd25hbWVzKCdwbGFjZScpICU+JQogIHNlbGVjdF9pZihpc19udW1lcmljKSAlPiUKICBoY3V0KGhjX2Z1bmMgPSAiaGNsdXN0IiwgCiAgICAgICBrID0gMywgCiAgICAgICBzdGFuZCA9IFRSVUUpCmBgYAoKVmlzdWFsaXplIGNsdXN0ZXJzCgpgYGB7ciwsZmlnLndpZHRoPTE1LGZpZy5oZWlnaHQ9MTAsZmlnLmFsaWduPSdjZW50ZXInfQpoYyAlPiUgCiAgZnZpel9jbHVzdGVyKGRhdGEgPSBkYXRhICU+JSBzZWxlY3RfaWYoaXNfbnVtZXJpYykpICAKYGBgCgpXaGVyZSBkbyB3ZSBmaW5kIHRoZSBjbHVzdGVycz8KCmBgYHtyfQpoYyAlPiUKICBnbGltcHNlKCkKYGBgCgpgYGB7cn0KaGMkY2x1c3RlciAlPiUKICBoZWFkKCkKYGBgCgpBZGQgdGhlbSB0byBkYXRhc2V0CgpgYGB7cn0KZGF0YVssImNsdXN0ZXIiXSA8LSBoYyRjbHVzdGVyCmBgYAoKSW5zcGVjdCBjbHVzdGVycyBwZXIgcmVnaW9uCgpgYGB7cn0KdGFibGUoZGF0YSRjbHVzdGVyLCBkYXRhJHJlZ2lvbikKYGBgCgpBbHNvIGFkZCBQQ0EgdG8gb3JpZ25hbCBkYXRhCgpgYGB7cn0KZGF0YVssInBjYTEiXSA8LSByZXNfcGNhJGluZCRjb29yZFssMV0KZGF0YVssInBjYTIiXSA8LSByZXNfcGNhJGluZCRjb29yZFssMl0KYGBgCgpDb21wb25lbnQgbWVhbiBwZXIgY2x1c3RlcgoKYGBge3J9CmRhdGEgJT4lCiAgZ3JvdXBfYnkoY2x1c3RlcikgJT4lCiAgc3VtbWFyaXNlKHBjYTEgPSBwY2ExICU+JSBtZWFuKCksCiAgICAgICAgICAgIHBjYTIgPSBwY2EyICU+JSBtZWFuKCkpCmBgYApCb251czogYWRkIHRyaXBzIGRhdGEKCmBgYHtyfQp0cmlwcyA8LSByZWFkX2NzdignaHR0cHM6Ly9zZHMtYWF1LmdpdGh1Yi5pby9TRFMtbWFzdGVyL00xL2RhdGEvdHJpcHMuY3N2JykKdHJpcHMgJT4lIGdsaW1wc2UoKQpgYGAKCkFkZCBudW1iZXIgb2YgdHJpcHMgcGVyIGNpdHkKCmBgYHtyfQpkYXRhICU8PiUKICBsZWZ0X2pvaW4odHJpcHMgJT4lIGNvdW50KHBsYWNlLCBzb3J0ID0gVFJVRSwgbmFtZSA9ICduX2NpdHknKSwgYnkgPSAncGxhY2UnKQpgYGAKCkNoZWNrIG1vc3QgcG9wdWxhciBjaXRpZXMgcGVyIGNsdXN0ZXIKCmBgYHtyfQpkYXRhICU+JQogIHNlbGVjdChwbGFjZSwgY2x1c3Rlciwgbl9jaXR5KSAlPiUKICBncm91cF9ieShjbHVzdGVyKSAlPiUKICBhcnJhbmdlKGRlc2Mobl9jaXR5KSkgJT4lCiAgc2xpY2UoMTo1KSAlPiUKICB1bmdyb3VwKCkgCmBgYApDb3VudCBjbHVzdGVyIHBvcHVsYXJpdHkKCmBgYHtyfQpkYXRhICU+JQogIGNvdW50KGNsdXN0ZXIsIHd0ID0gbl9jaXR5KQpgYGAKClRvIGZpbmlzaCB1cCwgbGV0cyBwbG90IGl0IGluIGEgbWFwLCBzaW1wbGVzdCB3YXkgcG9zc2libGUuCgpgYGB7cn0KZ2VvX21lcmdlIDwtIHRyaXBzICU+JQogIHNlbGVjdChwbGFjZSwgbG9uZ2l0dWRlLCBsYXRpdHVkZSkgJT4lCiAgZGlzdGluY3QocGxhY2UsIC5rZWVwX2FsbCA9IFRSVUUpCmBgYAoKYGBge3J9CmRhdGEgJTw+JQogIGxlZnRfam9pbihnZW9fbWVyZ2UgLCBieSA9ICdwbGFjZScpCmBgYAoKTG9hZCBhIHdvcmxkbWFwIGdlb20KCmBgYHtyLCBmaWcuaGVpZ2h0PTEwLCBmaWcud2lkdGg9MTV9CmxpYnJhcnkoZ2dtYXApCm1hcFdvcmxkIDwtIGJvcmRlcnMoIndvcmxkIiwgY29sb3VyID0gImdyYXk1MCIsIGZpbGwgPSAiZ3JheTUwIikKYGBgCgpHVGhhdCdzIGhvdyBpdCBsb29rcwoKYGBge3J9Cm1hcFdvcmxkCmBgYAoKQWRkIGl0IHRvIGFuIGVtcHR5IGdncGxvdCBzdXJmYWNlCgpgYGB7cn0KbXAgPC0gZ2dwbG90KCkgKyAgIAogIG1hcFdvcmxkIApgYGAKCgpUaGF0J3MgaG93IGl0IGxvb2tzIHNvIGZhcgoKYGBge3IsLGZpZy53aWR0aD0xNSxmaWcuaGVpZ2h0PTEwLGZpZy5hbGlnbj0nY2VudGVyJ30KbXAKYGBgCgpBZGQgYSBnZW9tIHdpdGggdGhlIGNpdGllcyBhcyBwb2ludHMKCmBgYHtyfQpub21hZF9tYXAgPC0gbXAgKyAKICBnZW9tX3BvaW50KGRhdGEgPSBkYXRhLCBhZXMoeCA9IGxvbmdpdHVkZSwgeSA9IGxhdGl0dWRlLCBjb2wgPSBmYWN0b3IoY2x1c3RlcikpKSAKYGBgCgpgYGB7ciwsZmlnLndpZHRoPTE1LGZpZy5oZWlnaHQ9MTAsZmlnLmFsaWduPSdjZW50ZXInfQpub21hZF9tYXAKYGBgCgpPciBkbyBhIGRlbnNpdHkgcGxvdCBvZiBwb3B1bGFyIG5vbWFkIGNpdGllcwoKYGBge3IsLGZpZy53aWR0aD0xNSxmaWcuaGVpZ2h0PTEwLGZpZy5hbGlnbj0nY2VudGVyJ30KbXAgKyAKICBzdGF0X2RlbnNpdHkyZChkYXRhID0gdHJpcHMsIAogICAgICAgICAgICAgICAgIGFlcyh4ID0gbG9uZ2l0dWRlLCB5ID0gbGF0aXR1ZGUsIGZpbGwgPSBzdGF0KG5sZXZlbCksIGNvbCA9IHN0YXQobmxldmVsKSApLCAKICAgICAgICAgICAgICAgICBhbHBoYSA9IDAuMiwgc2l6ZSA9IDAuMiwgYmlucyA9IDEwLCBnZW9tID0gInBvbHlnb24iKSArCiAgc2NhbGVfZmlsbF9ncmFkaWVudChsb3cgPSAic2t5Ymx1ZSIsIGhpZ2ggPSAicmVkIikgKwogIHNjYWxlX2NvbG9yX2dyYWRpZW50KGxvdyA9ICJza3libHVlIiwgaGlnaCA9ICJyZWQiKQpgYGAKCgogICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAg