R Notebook

correlation
similarity measures
plotting similarity

This session is based on chapter 6 of Levshina’s How to do linguistics with R. One of the aims of this course is that you should be able to go into any chapter of this textbook and follow the methods

Install the Rling package from https://benjamins.com/sites/z.195/content/package.html (or search Rling download).

Load up everything we need

library(tidyverse)

package ‘tidyverse’ was built under R version 3.3.2Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
package ‘ggplot2’ was built under R version 3.3.2package ‘readr’ was built under R version 3.3.2package ‘purrr’ was built under R version 3.3.2package ‘dplyr’ was built under R version 3.3.2Conflicts with tidy packages ------------------------------------------------------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats

library(Rling)
data(ldt) # lexical decision time data

Check out the data:

head(ldt)

Plot word length against reaction time.

ldt %>% ggplot(aes(x=Length, y=Mean_RT)) + geom_point() + geom_smooth(method="lm")

Intro to linear regression

Linear regression assumes that the observations (red points) show an underlying relationship between x and y (the blue line), which is obscured by some noise/random/other factors (green lines, called residuals). The residuals are the deviations from perfect correlation.

To calculate the linear model without plotting it, use the lm function:

rr m <- lm(Mean_RT ~ Length, data=ldt)

This contains all the data. You can look at the residuals:

#residuals(m) %>% head()
m$residuals %>% head()

    marveled   persuaders      midmost       crutch resuspension efflorescent 
      19.595      102.747      146.269       41.993      175.249       -1.841

Or plot them (use a histogram)

tibble(residuals=residuals(m)) %>% ggplot(aes(x=residuals)) + geom_histogram(bins=50)

Measure the correlation (using least squares, minimising the squares of the residuals)

cor(ldt$Mean_RT, ldt$Length)

[1] 0.6147

#help(cor)
#cor(Mean_RT, Length, data=ldt)

Test the strength of this correlation. Basically, is zero within the 95% confidence interval?

cor.test(ldt$Mean_RT, ldt$Length, alternative="greater")


    Pearson's product-moment correlation

data:  ldt$Mean_RT and ldt$Length
t = 7.7, df = 98, p-value = 5e-12
alternative hypothesis: true correlation is greater than 0
95 percent confidence interval:
 0.5001 1.0000
sample estimates:
   cor 
0.6147

Very high probability that the true correlation is greater than zero.

Compare with some random data:

ran.data <- tibble(a=rnorm(20, 10, 1), b=rnorm(20, 10, 1))
ran.data %>% ggplot(aes(x=a, y=b)) + geom_point() + geom_smooth(method="lm")

There might seem to be a correlation:


    Pearson's product-moment correlation

data:  ran.data$a and ran.data$b
t = 0.46, df = 8, p-value = 0.7
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.5214  0.7179
sample estimates:
   cor 
0.1611

But do we believe this? (use two-sided alternative: testing that the value is not zero)

cor.test(ran.data$a, ran.data$b, alternative="two.sided")


    Pearson's product-moment correlation

data:  ran.data$a and ran.data$b
t = 0.21, df = 18, p-value = 0.8
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4028  0.4805
sample estimates:
    cor 
0.04828

(at the p-value cut-off of 0.05, one in twenty times this demo won’t work)

ran.data <- tibble(a=runif(100, 0, 20))
#head(ran.data)
ran.data$b <- ran.data$a + rnorm(100, 0, 10)
head(ran.data)

ran.data %>% ggplot(aes(x=a, y=b)) + geom_point() + geom_smooth(method="lm")

Discussion point:

don’t confuse effect size with statistical significance!

Outliers

False correlations with outliers

Add an outlier

ran.data <- tibble(a=runif(20, 0, 10), b=runif(20, 0, 10))
ran.data %>% 
  arrange(desc(a)) %>% 
  mutate(rank=1:n()) %>% 
  mutate(b=if_else(rank < 4, b+100, b)) -> ran.data.2
ran.data.2 %>% ggplot(aes(x=a, y=b)) + geom_point() + geom_smooth(method="lm")

cor(ran.data.2$a, ran.data.2$b)

[1] 0.6366

cor.test(ran.data.2$a, ran.data.2$b, alternative="two.sided")


    Pearson's product-moment correlation

data:  ran.data.2$a and ran.data.2$b
t = 3.5, df = 18, p-value = 0.003
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2702 0.8420
sample estimates:
   cor 
0.6366

Parsed with column specification:
cols(
  wals.code = col_character(),
  name = col_character(),
  family = col_character(),
  F82A = col_integer(),
  F83A = col_integer(),
  F85A = col_integer(),
  F86A = col_integer(),
  F87A = col_integer(),
  F88A = col_integer(),
  F89A = col_integer(),
  F90A = col_integer()
)

rr head(wals)

Solution: remove outliers

Check ldt data for outliers in mean_RT
What might these mean?
Remove outliers and recalculate correlation

Assumptions of Pearson correlation

The sample is randomly selected from the population it represents
Both variables are interval-scaled
The sample size is large (> 30); or for smaller samples, both variables come from a bivariate normal distribution. This means: both variables are normally distributed, and their linear combination is also normally distributed.
The residual variance is homoscedastic). This means: the relationship between the variables is of equal strength across the entire range of both variables
The residuals are independent (no autocorrelation)

Distance

From here on is for Session #8

Data from The World Atlas of Linguistic Structures

wals <- read_tsv("./typ/wals_data.tsv")

Parsed with column specification:
cols(
  wals.code = col_character(),
  name = col_character(),
  family = col_character(),
  F82A = col_integer(),
  F83A = col_integer(),
  F85A = col_integer(),
  F86A = col_integer(),
  F87A = col_integer(),
  F88A = col_integer(),
  F89A = col_integer(),
  F90A = col_integer()
)

head(wals)

Check the data quality. How much is missing?

wals %>% group_by(wals.code) %>% 
  mutate(N=sum(!is.na(F82A), !is.na(F83A), !is.na(F85A), !is.na(F86A), !is.na(F87A), !is.na(F88A), !is.na(F89A), !is.na(F90A))) %>%
  ggplot(aes(x=N)) + geom_histogram(bins=10)

Let’s just take the complete samples

wals %>% group_by(wals.code) %>% 
  mutate(N=sum(!is.na(F82A), !is.na(F83A), !is.na(F85A), !is.na(F86A), !is.na(F87A), !is.na(F88A), !is.na(F89A), !is.na(F90A))) %>%
  filter(N==8) -> wals.554

554 rows.

Now let’s make a distance matrix. We start from an empty matrix, then fill it pair by pair

rr biplot(prcomp(as.dist(m)))

Now for the tricky bit:

rr prcomp(as.dist(m)) -> p p$x[1:8, 1:2]

        PC1     PC2
aar -12.592  34.078
abk -50.578   8.102
abu  26.597  20.667
abv -21.273  16.571
ace  25.308 -18.981
acn -39.788  20.811
adi   7.033  27.869
adz  26.597  20.667

biplot(prcomp(as.dist(m)))

prcomp(as.dist(m)) -> p
p$x[1:8, 1:2]

        PC1     PC2
aar -12.592  34.078
abk -50.578   8.102
abu  26.597  20.667
abv -21.273  16.571
ace  25.308 -18.981
acn -39.788  20.811
adi   7.033  27.869
adz  26.597  20.667

as_tibble(p$x[,1:2], rownames="wals_code") %>% 
  mutate(family=wals.554$family) %>%
  mutate(family=if_else(family %in% c("Austronesian", "Indo-European", "Niger-Congo"), family, "other")) %>%
  ggplot(aes(x=PC1, y=PC2, colour=family)) + geom_point() + scale_colour_manual(values=c("darkgreen", "blue", "brown", "gray"))

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKLSBjb3JyZWxhdGlvbgotIHNpbWlsYXJpdHkgbWVhc3VyZXMKLSBwbG90dGluZyBzaW1pbGFyaXR5CgpUaGlzIHNlc3Npb24gaXMgYmFzZWQgb24gY2hhcHRlciA2IG9mIExldnNoaW5hJ3MgKkhvdyB0byBkbyBsaW5ndWlzdGljcyB3aXRoIFIqLiBPbmUgb2YgdGhlIGFpbXMgb2YgdGhpcyBjb3Vyc2UgaXMgdGhhdCB5b3Ugc2hvdWxkIGJlIGFibGUgdG8gZ28gaW50byBhbnkgY2hhcHRlciBvZiB0aGlzIHRleHRib29rIGFuZCBmb2xsb3cgdGhlIG1ldGhvZHMgCgpJbnN0YWxsIHRoZSBgUmxpbmdgIHBhY2thZ2UgZnJvbSBodHRwczovL2JlbmphbWlucy5jb20vc2l0ZXMvei4xOTUvY29udGVudC9wYWNrYWdlLmh0bWwgKG9yIHNlYXJjaCAqUmxpbmcgZG93bmxvYWQqKS4KCkxvYWQgdXAgZXZlcnl0aGluZyB3ZSBuZWVkCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShSbGluZykKZGF0YShsZHQpICMgbGV4aWNhbCBkZWNpc2lvbiB0aW1lIGRhdGEKYGBgCgpDaGVjayBvdXQgdGhlIGRhdGE6CmBgYHtyfQpoZWFkKGxkdCkKYGBgClBsb3Qgd29yZCBsZW5ndGggYWdhaW5zdCByZWFjdGlvbiB0aW1lLgoKYGBge3J9CmxkdCAlPiUgZ2dwbG90KGFlcyh4PUxlbmd0aCwgeT1NZWFuX1JUKSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX3Ntb290aChtZXRob2Q9ImxtIikKYGBgCgpJbnRybyB0byBsaW5lYXIgcmVncmVzc2lvbgoKIVtdKC4vMjIwcHgtTGluZWFyX2xlYXN0X3NxdWFyZXNfZXhhbXBsZTIucG5nKQpMaW5lYXIgcmVncmVzc2lvbiBhc3N1bWVzIHRoYXQgdGhlIG9ic2VydmF0aW9ucyAocmVkIHBvaW50cykgc2hvdyBhbiB1bmRlcmx5aW5nIHJlbGF0aW9uc2hpcCBiZXR3ZWVuIHggYW5kIHkgKHRoZSBibHVlIGxpbmUpLCB3aGljaCBpcyBvYnNjdXJlZCBieSBzb21lIG5vaXNlL3JhbmRvbS9vdGhlciBmYWN0b3JzIChncmVlbiBsaW5lcywgY2FsbGVkICoqcmVzaWR1YWxzKiopLiBUaGUgcmVzaWR1YWxzIGFyZSB0aGUgZGV2aWF0aW9ucyBmcm9tIHBlcmZlY3QgY29ycmVsYXRpb24uCgpUbyBjYWxjdWxhdGUgdGhlIGxpbmVhciBtb2RlbCB3aXRob3V0IHBsb3R0aW5nIGl0LCB1c2UgdGhlIGxtIGZ1bmN0aW9uOgoKYGBge3J9Cm0gPC0gbG0oTWVhbl9SVCB+IExlbmd0aCwgZGF0YT1sZHQpCiNuYW1lcyhtKQojbSRmaXR0ZWQudmFsdWVzCmBgYAoKVGhpcyBjb250YWlucyBhbGwgdGhlIGRhdGEuIFlvdSBjYW4gbG9vayBhdCB0aGUgcmVzaWR1YWxzOgpgYGB7cn0KI3Jlc2lkdWFscyhtKSAlPiUgaGVhZCgpCm0kcmVzaWR1YWxzICU+JSBoZWFkKCkKYGBgCk9yIHBsb3QgdGhlbSAodXNlIGEgaGlzdG9ncmFtKQpgYGB7cn0KdGliYmxlKHJlc2lkdWFscz1yZXNpZHVhbHMobSkpICU+JSBnZ3Bsb3QoYWVzKHg9cmVzaWR1YWxzKSkgKyBnZW9tX2hpc3RvZ3JhbShiaW5zPTUwKQpgYGAKCk1lYXN1cmUgdGhlIGNvcnJlbGF0aW9uICh1c2luZyAqbGVhc3Qgc3F1YXJlcyosIG1pbmltaXNpbmcgdGhlIHNxdWFyZXMgb2YgdGhlIHJlc2lkdWFscykKCmBgYHtyfQpjb3IobGR0JE1lYW5fUlQsIGxkdCRMZW5ndGgpCiNoZWxwKGNvcikKI2NvcihNZWFuX1JULCBMZW5ndGgsIGRhdGE9bGR0KQpgYGAKClRlc3QgdGhlIHN0cmVuZ3RoIG9mIHRoaXMgY29ycmVsYXRpb24uIEJhc2ljYWxseSwgaXMgemVybyB3aXRoaW4gdGhlIDk1JSBjb25maWRlbmNlIGludGVydmFsPwpgYGB7cn0KY29yLnRlc3QobGR0JE1lYW5fUlQsIGxkdCRMZW5ndGgsIGFsdGVybmF0aXZlPSJncmVhdGVyIikKYGBgClZlcnkgaGlnaCBwcm9iYWJpbGl0eSB0aGF0IHRoZSB0cnVlIGNvcnJlbGF0aW9uIGlzIGdyZWF0ZXIgdGhhbiB6ZXJvLgoKQ29tcGFyZSB3aXRoIHNvbWUgcmFuZG9tIGRhdGE6CmBgYHtyfQpyYW4uZGF0YSA8LSB0aWJibGUoYT1ybm9ybSgyMCwgMTAsIDEpLCBiPXJub3JtKDIwLCAxMCwgMSkpCnJhbi5kYXRhICU+JSBnZ3Bsb3QoYWVzKHg9YSwgeT1iKSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX3Ntb290aChtZXRob2Q9ImxtIikKYGBgCgpUaGVyZSBtaWdodCBzZWVtIHRvIGJlIGEgY29ycmVsYXRpb246CmBgYHtyfQpjb3IocmFuLmRhdGEkYSwgcmFuLmRhdGEkYikKYGBgCkJ1dCBkbyB3ZSBiZWxpZXZlIHRoaXM/ICh1c2UgdHdvLXNpZGVkIGFsdGVybmF0aXZlOiB0ZXN0aW5nIHRoYXQgdGhlIHZhbHVlIGlzIG5vdCB6ZXJvKQpgYGB7cn0KY29yLnRlc3QocmFuLmRhdGEkYSwgcmFuLmRhdGEkYiwgYWx0ZXJuYXRpdmU9InR3by5zaWRlZCIpCmBgYAooYXQgdGhlIHAtdmFsdWUgY3V0LW9mZiBvZiAwLjA1LCBvbmUgaW4gdHdlbnR5IHRpbWVzIHRoaXMgZGVtbyB3b24ndCB3b3JrKQoKYGBge3J9CnJhbi5kYXRhIDwtIHRpYmJsZShhPXJ1bmlmKDEwMCwgMCwgMjApKQojaGVhZChyYW4uZGF0YSkKcmFuLmRhdGEkYiA8LSByYW4uZGF0YSRhICsgcm5vcm0oMTAwLCAwLCAxMCkKaGVhZChyYW4uZGF0YSkKYGBgCmBgYHtyfQpyYW4uZGF0YSAlPiUgZ2dwbG90KGFlcyh4PWEsIHk9YikpICsgZ2VvbV9wb2ludCgpICsgZ2VvbV9zbW9vdGgobWV0aG9kPSJsbSIpCmBgYApgYGB7cn0KY29yLnRlc3QocmFuLmRhdGEkYSwgcmFuLmRhdGEkYikKYGBgCgoKRGlzY3Vzc2lvbiBwb2ludDoKCi0gZG9uJ3QgY29uZnVzZSBlZmZlY3Qgc2l6ZSB3aXRoIHN0YXRpc3RpY2FsIHNpZ25pZmljYW5jZSEKCiMjIE91dGxpZXJzCgotIEZhbHNlIGNvcnJlbGF0aW9ucyB3aXRoIG91dGxpZXJzCgpBZGQgYW4gb3V0bGllcgpgYGB7cn0KcmFuLmRhdGEgPC0gdGliYmxlKGE9cnVuaWYoMjAsIDAsIDEwKSwgYj1ydW5pZigyMCwgMCwgMTApKQpyYW4uZGF0YSAlPiUgCiAgYXJyYW5nZShkZXNjKGEpKSAlPiUgCiAgbXV0YXRlKHJhbms9MTpuKCkpICU+JSAKICBtdXRhdGUoYj1pZl9lbHNlKHJhbmsgPCA0LCBiKzEwMCwgYikpIC0+IHJhbi5kYXRhLjIKcmFuLmRhdGEuMiAlPiUgZ2dwbG90KGFlcyh4PWEsIHk9YikpICsgZ2VvbV9wb2ludCgpICsgZ2VvbV9zbW9vdGgobWV0aG9kPSJsbSIpCmBgYAoKYGBge3J9CmNvcihyYW4uZGF0YS4yJGEsIHJhbi5kYXRhLjIkYikKYGBgCgpgYGB7cn0KY29yLnRlc3QocmFuLmRhdGEuMiRhLCByYW4uZGF0YS4yJGIsIGFsdGVybmF0aXZlPSJ0d28uc2lkZWQiKQpgYGAKYGBge3J9CnJhbi5kYXRhLjIgJT4lIGZpbHRlcihiIDwgNTApIC0+IHJhbi5kYXRhLjMKY29yLnRlc3QocmFuLmRhdGEuMyRhLCByYW4uZGF0YS4zJGIpCmBgYAoKKipTb2x1dGlvbjoqKiByZW1vdmUgb3V0bGllcnMKCi0gQ2hlY2sgbGR0IGRhdGEgZm9yIG91dGxpZXJzIGluIG1lYW5fUlQKLSBXaGF0IG1pZ2h0IHRoZXNlIG1lYW4/Ci0gUmVtb3ZlIG91dGxpZXJzIGFuZCByZWNhbGN1bGF0ZSBjb3JyZWxhdGlvbgoKIyMgQXNzdW1wdGlvbnMgb2YgUGVhcnNvbiBjb3JyZWxhdGlvbgoKLSBUaGUgc2FtcGxlIGlzIHJhbmRvbWx5IHNlbGVjdGVkIGZyb20gdGhlIHBvcHVsYXRpb24gaXQgcmVwcmVzZW50cwotIEJvdGggdmFyaWFibGVzIGFyZSBpbnRlcnZhbC1zY2FsZWQKLSBUaGUgc2FtcGxlIHNpemUgaXMgbGFyZ2UgKD4gMzApOyBvciBmb3Igc21hbGxlciBzYW1wbGVzLCBib3RoIHZhcmlhYmxlcyBjb21lIGZyb20gYSBiaXZhcmlhdGUgbm9ybWFsIGRpc3RyaWJ1dGlvbi4gVGhpcyBtZWFuczogYm90aCB2YXJpYWJsZXMgYXJlIG5vcm1hbGx5IGRpc3RyaWJ1dGVkLCBhbmQgdGhlaXIgbGluZWFyIGNvbWJpbmF0aW9uIGlzIGFsc28gbm9ybWFsbHkgZGlzdHJpYnV0ZWQuCi0gVGhlIHJlc2lkdWFsIHZhcmlhbmNlIGlzICpob21vc2NlZGFzdGljKikuIFRoaXMgbWVhbnM6IHRoZSByZWxhdGlvbnNoaXAgYmV0d2VlbiB0aGUgdmFyaWFibGVzIGlzIG9mIGVxdWFsIHN0cmVuZ3RoIGFjcm9zcyB0aGUgZW50aXJlIHJhbmdlIG9mIGJvdGggdmFyaWFibGVzCi0gVGhlIHJlc2lkdWFscyBhcmUgaW5kZXBlbmRlbnQgKG5vIGF1dG9jb3JyZWxhdGlvbikKCiMjIERpc3RhbmNlCgoqKkZyb20gaGVyZSBvbiBpcyBmb3IgU2Vzc2lvbiAjOCoqCgpEYXRhIGZyb20gW1RoZSBXb3JsZCBBdGxhcyBvZiBMaW5ndWlzdGljIFN0cnVjdHVyZXNdKGh0dHBzOi8vd2Fscy5pbmZvKQoKYGBge3J9CndhbHMgPC0gcmVhZF90c3YoIi4vdHlwL3dhbHNfZGF0YS50c3YiKQpoZWFkKHdhbHMpCmBgYAoKQ2hlY2sgdGhlIGRhdGEgcXVhbGl0eS4gSG93IG11Y2ggaXMgbWlzc2luZz8KCmBgYHtyfQp3YWxzICU+JSBncm91cF9ieSh3YWxzLmNvZGUpICU+JSAKICBtdXRhdGUoTj1zdW0oIWlzLm5hKEY4MkEpLCAhaXMubmEoRjgzQSksICFpcy5uYShGODVBKSwgIWlzLm5hKEY4NkEpLCAhaXMubmEoRjg3QSksICFpcy5uYShGODhBKSwgIWlzLm5hKEY4OUEpLCAhaXMubmEoRjkwQSkpKSAlPiUKICBnZ3Bsb3QoYWVzKHg9TikpICsgZ2VvbV9oaXN0b2dyYW0oYmlucz0xMCkKYGBgCkxldCdzIGp1c3QgdGFrZSB0aGUgY29tcGxldGUgc2FtcGxlcwoKYGBge3J9CndhbHMgJT4lIGdyb3VwX2J5KHdhbHMuY29kZSkgJT4lIAogIG11dGF0ZShOPXN1bSghaXMubmEoRjgyQSksICFpcy5uYShGODNBKSwgIWlzLm5hKEY4NUEpLCAhaXMubmEoRjg2QSksICFpcy5uYShGODdBKSwgIWlzLm5hKEY4OEEpLCAhaXMubmEoRjg5QSksICFpcy5uYShGOTBBKSkpICU+JQogIGZpbHRlcihOPT04KSAtPiB3YWxzLjU1NApgYGAKCjU1NCByb3dzLgoKTm93IGxldCdzIG1ha2UgYSBkaXN0YW5jZSBtYXRyaXguIFdlIHN0YXJ0IGZyb20gYW4gZW1wdHkgbWF0cml4LCB0aGVuIGZpbGwgaXQgcGFpciBieSBwYWlyCgpgYGB7cn0KbSA8LSBtYXRyaXgoTkEsIG5yb3c9NTU0LCBuY29sPTU1NCkKd2Fscy5jb2RlczwtIHdhbHMuNTU0ICU+JSBwdWxsKHdhbHMuY29kZSkKcm93bmFtZXMobSkgPC0gd2Fscy5jb2Rlcwpjb2xuYW1lcyhtKSA8LSB3YWxzLmNvZGVzCm1bMTo4LDE6OF0KYGBgCk5vdyBmb3IgdGhlIHRyaWNreSBiaXQ6CmBgYHtyfQoKd2Fscy41NTQgJT4lIHNlbGVjdCgiRjgyQSIsICJGODNBIiwgIkY4NUEiLCAiRjg2QSIsICJGODdBIiwgIkY4OEEiLCAiRjg5QSIsICJGOTBBIikgJT4lIGFzLmRhdGEuZnJhbWUoKSAtPiB3YWxzLmRmCnJvdy5uYW1lcyh3YWxzLmRmKSA8LSB3YWxzLmNvZGVzCgppID0gMApmb3IgKGxnMSBpbiB3YWxzLmNvZGVzKXsKICBpID0gaSArIDEKICBpZiAoaSAlJSAxMCA9PSAwKSBjYXQocGFzdGUoaSwgbGcxLCAiXG4iKSkKICBmb3IgKGxnMiBpbiB3YWxzLmNvZGVzKXsKICAgIGlmIChsZzEgPCBsZzIpewogICAgICBsZzEudmFsdWVzIDwtIHdhbHMuZGZbbGcxLCBdCiAgICAgIGxnMi52YWx1ZXMgPC0gd2Fscy5kZltsZzIsIF0KICAgICAgZCA8LSBzdW0obGcxLnZhbHVlcyAhPSBsZzIudmFsdWVzKQogICAgICBtW2xnMSwgbGcyXSA8LSBkCiAgICAgIG1bbGcyLCBsZzFdIDwtIGQKICAgIH0KICB9Cn0KbVsxOjgsMTo4XQpgYGAKCmBgYHtyfQpiaXBsb3QocHJjb21wKGFzLmRpc3QobSkpKQpgYGAKCmBgYHtyfQpwcmNvbXAoYXMuZGlzdChtKSkgLT4gcApwJHhbMTo4LCAxOjJdCmBgYAoKYGBge3J9CmFzX3RpYmJsZShwJHhbLDE6Ml0sIHJvd25hbWVzPSJ3YWxzX2NvZGUiKSAlPiUgCiAgbXV0YXRlKGZhbWlseT13YWxzLjU1NCRmYW1pbHkpICU+JQogIG11dGF0ZShmYW1pbHk9aWZfZWxzZShmYW1pbHkgJWluJSBjKCJBdXN0cm9uZXNpYW4iLCAiSW5kby1FdXJvcGVhbiIsICJOaWdlci1Db25nbyIpLCBmYW1pbHksICJvdGhlciIpKSAlPiUKICBnZ3Bsb3QoYWVzKHg9UEMxLCB5PVBDMiwgY29sb3VyPWZhbWlseSkpICsgZ2VvbV9wb2ludCgpICsgc2NhbGVfY29sb3VyX21hbnVhbCh2YWx1ZXM9YygiZGFya2dyZWVuIiwgImJsdWUiLCAiYnJvd24iLCAiZ3JheSIpKQpgYGAKCgo=