Hierarchical model for Rats experiment (BDA3, p. 102)

ggplot2, grid, and gridExtra are used for plotting, tidyr for manipulating data frames

library(ggplot2)
theme_set(theme_minimal())
library(gridExtra)
library(tidyr)
library(latex2exp)

Data

y <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,
        2,1,5,2,5,3,2,7,7,3,3,2,9,10,4,4,4,4,4,4,4,10,4,4,4,5,11,12,
        5,5,6,5,6,6,6,6,16,15,15,9,4)
n <- c(20,20,20,20,20,20,20,19,19,19,19,18,18,17,20,20,20,20,19,19,18,18,25,24,
       23,20,20,20,20,20,20,10,49,19,46,27,17,49,47,20,20,13,48,50,20,20,20,20,
       20,20,20,48,19,19,19,22,46,49,20,20,23,19,22,20,20,20,52,46,47,24,14)

Evaluate densities in grid

x <- seq(0.0001, 0.9999, length.out = 1000)

Helper function to evaluate density over observations

bdens <- function(n, y, x)
    dbeta(x, y+1, n-y+1)

Separate model

df_sep <- mapply(bdens, n, y, MoreArgs = list(x = x)) %>%
  as.data.frame() %>% cbind(x) %>% 
  pivot_longer(cols = !x, names_to = "ind", values_to = "p")

Plot the separate model

labs1 <- paste('posterior of', c('theta_j', 'theta_71'))
plot_sep <- ggplot(data = df_sep) +
  geom_line(aes(x = x, y = p, color = (ind=='V71'), group = ind)) +
  labs(x = expression(theta), y = '', title = 'Separate model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = c('blue','red'), labels = labs1) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))
# The last one is for emphasize colored red
plot_sep

Pooled model

df_pool <- data.frame(x = x, p = dbeta(x, sum(y)+1, sum(n)-sum(y)+1))

Create a plot for the separate model

plot_pool <- ggplot(data = df_pool) +
  geom_line(aes(x = x, y = p, color = '1')) +
  labs(x = expression(theta), y = '', title = 'Pooled model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = 'red', labels = 'Posterior of common theta') +
  theme(legend.background = element_blank(), legend.position = c(0.7,0.9))

Plot both separate and pooled model

grid.arrange(plot_sep, plot_pool)

Compute the marginal posterior of alpha and beta in hierarchical model Use grid

A <- seq(0.5, 6, length.out = 100)
B <- seq(3, 33, length.out = 100)

Make vectors that contain all pairwise combinations of A and B

cA <- rep(A, each = length(B))
cB <- rep(B, length(A))

Use logarithms for numerical accuracy!

lpfun <- function(a, b, y, n) log(a+b)*(-5/2) +
  sum(lgamma(a+b)-lgamma(a)-lgamma(b)+lgamma(a+y)+lgamma(b+n-y)-lgamma(a+b+n))
lp <- mapply(lpfun, cA, cB, MoreArgs = list(y, n))

Subtract maximum value to avoid over/underflow in exponentation

df_marg <- data.frame(x = cA, y = cB, p = exp(lp - max(lp)))

Create a plot of the marginal posterior density

title1 <- TeX('The marginal of $\\alpha$ and $\\beta$')
ggplot(data = df_marg, aes(x = x, y = y)) +
  geom_raster(aes(fill = p, alpha = p), interpolate = T) +
  geom_contour(aes(z = p), colour = 'black', size = 0.2) +
  coord_cartesian(xlim = c(1,5), ylim = c(4, 26)) +
  labs(x = TeX('$\\alpha$'), y = TeX('$\\beta$'), title = title1) +
  scale_fill_gradient(low = 'yellow', high = 'red', guide = F) +
  scale_alpha(range = c(0, 1), guide = F)
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.
## It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

Sample from the grid (with replacement)

nsamp <- 100
samp_indices <- sample(length(df_marg$p), size = nsamp,
                       replace = T, prob = df_marg$p/sum(df_marg$p))
samp_A <- cA[samp_indices[1:nsamp]]
samp_B <- cB[samp_indices[1:nsamp]]
df_psamp <- mapply(function(a, b, x) dbeta(x, a, b),
                  samp_A, samp_B, MoreArgs = list(x = x)) %>%
  as.data.frame() %>% cbind(x) %>% 
  pivot_longer(cols = !x, names_to = "ind", values_to = "p")

Create plot for samples from the distribution of distributions Beta(alpha,beta), that is, plot Beta(alpha,beta) using posterior samples of alpha and beta

# helper function to convert ind to numeric for subsetting
indtonum <- function(x) strtoi(substring(x,2))
title2 <- TeX('Beta($\\alpha,\\beta$) given posterior draws of $\\alpha$ and $\\beta$')
plot_psamp <- ggplot(data = subset(df_psamp, indtonum(ind) <= 20)) +
  geom_line(aes(x = x, y = p, group = ind), color='forestgreen') +
  labs(x = expression(theta), y = '', title = title2) +
  scale_y_continuous(breaks = NULL)

The average of above distributions, is the predictive distribution for a new theta, and also the prior distribution for theta_j

df_psampmean <- spread(df_psamp, ind, p) %>% subset(select = -x) %>%
    rowMeans() %>% data.frame(x = x, p = .)

Create plot for samples from the predictive distribution for new theta

title3 <- TeX('Population distribution (prior) for $\\theta_j$')
plot_psampmean <- ggplot(data = df_psampmean) +
  geom_line(aes(x = x, y = p), color='forestgreen') +
  labs(x = expression(theta), y = '', title = title3) +
  scale_y_continuous(breaks = NULL)

Combine the plots

grid.arrange(plot_psamp, plot_psampmean)

And finally compare the separate model and hierarchical model (using every seventh sample for clarity) Create plot for the separate model

plot_sep7 <- ggplot(data = subset(df_sep, indtonum(ind)%%7==0)) +
  geom_line(aes(x = x, y = p, color = (ind=='V49'), group = ind)) +
  labs(x = expression(theta), y = '', title = 'Separate model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = c('blue', 'red'), guide = F) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))

The hierarchical model

Note that these marginal posteriors for theta_j are more narrow than in the separate model case, due to the borrowed information from the other theta_j’s

Average density over samples (of a and b) for each (n,y)-pair at each point x

bdens2 <- function(n, y, a, b, x)
  rowMeans(mapply(dbeta, a + y, n - y + b, MoreArgs = list(x = x)))
df_hier <- mapply(bdens2, n, y, MoreArgs = list(samp_A, samp_B, x)) %>%
  as.data.frame() %>% cbind(x) %>% 
  pivot_longer(cols = !x, names_to = "ind", values_to = "p")

Create plot for the hierarchical model

plot_hier7 <- ggplot(data = subset(df_hier, indtonum(ind)%%7==0)) +
  geom_line(aes(x = x, y = p, color = (ind=='V49'), group = ind)) +
  labs(x = expression(theta), y = '', title = 'Hierarchical model', color = '') +
  scale_color_manual(values = c('blue', 'red'), guide = F) +
  scale_y_continuous(breaks = NULL) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))

Combine the plots

grid.arrange(plot_sep7, plot_hier7)
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.
## It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

90% intervals for separate models

qq_separate<-data.frame(id=1:length(n), n=n, y=y,
               q10=qbeta(0.1,y+1,n-y+1), q50=qbeta(0.5,y+1,n-y+1), q90=qbeta(0.9,y+1,n-y+1))

90% intervals for hierarchical model

qh <- function(q, n, y) colMeans(mapply(function(q, n, y, a, b)
    mapply(qbeta, q, a + y, n - y + b), q, n, y, MoreArgs = list(samp_A, samp_B)))
qq_hier <- data.frame(id=1:length(n), n=n, y=y,
                      q10 = qh(0.1, n, y), q50 = qh(0.5, n, y), q90 = qh(0.9, n, y))

pooled mean

m_pooled <- (sum(y)+1)/(sum(n)+2)

plot

plot_sep_int <- qq_separate %>%
  ggplot(aes(x=n,y=q50,ymin=q10,ymax=q90)) +
  geom_pointrange(color="blue", alpha=0.2) +
  geom_hline(yintercept=m_pooled, linetype="dashed")+
  lims(x=c(0,60), y=c(0,0.51))+
  labs(x="N", y="Probability of polyps", title="Separate")+
  annotate("text", x=0, y=m_pooled, label = "Pooled mean", vjust=-0.2, hjust=0.3)
plot_hier_int <- qq_hier %>%
  ggplot(aes(x=n,y=q50,ymin=q10,ymax=q90)) +
  geom_pointrange(color="blue", alpha=0.2) +
  geom_hline(yintercept=m_pooled, linetype="dashed")+
  lims(x=c(0,60), y=c(0,0.51))+
  labs(x="N", y="Probability of polyps", title="Hierarchical")+
  annotate("text", x=0, y=m_pooled, label = "Pooled mean", vjust=-0.2, hjust=0.3)
grid.arrange(plot_sep_int, plot_hier_int)

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDUuMSIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIEhpZXJhcmNoaWNhbCBtb2RlbCBmb3IgUmF0cyBleHBlcmltZW50IChCREEzLCBwLiAxMDIpCgpnZ3Bsb3QyLCBncmlkLCBhbmQgZ3JpZEV4dHJhIGFyZSB1c2VkIGZvciBwbG90dGluZywgdGlkeXIgZm9yCm1hbmlwdWxhdGluZyBkYXRhIGZyYW1lcwoKYGBge3Igc2V0dXAsIG1lc3NhZ2U9RkFMU0UsIGVycm9yPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KGdncGxvdDIpCnRoZW1lX3NldCh0aGVtZV9taW5pbWFsKCkpCmxpYnJhcnkoZ3JpZEV4dHJhKQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KGxhdGV4MmV4cCkKYGBgCgpEYXRhCgpgYGB7ciB9CnkgPC0gYygwLDAsMCwwLDAsMCwwLDAsMCwwLDAsMCwwLDAsMSwxLDEsMSwxLDEsMSwxLDIsMiwyLDIsMiwyLDIsMiwKICAgICAgICAyLDEsNSwyLDUsMywyLDcsNywzLDMsMiw5LDEwLDQsNCw0LDQsNCw0LDQsMTAsNCw0LDQsNSwxMSwxMiwKICAgICAgICA1LDUsNiw1LDYsNiw2LDYsMTYsMTUsMTUsOSw0KQpuIDwtIGMoMjAsMjAsMjAsMjAsMjAsMjAsMjAsMTksMTksMTksMTksMTgsMTgsMTcsMjAsMjAsMjAsMjAsMTksMTksMTgsMTgsMjUsMjQsCiAgICAgICAyMywyMCwyMCwyMCwyMCwyMCwyMCwxMCw0OSwxOSw0NiwyNywxNyw0OSw0NywyMCwyMCwxMyw0OCw1MCwyMCwyMCwyMCwyMCwKICAgICAgIDIwLDIwLDIwLDQ4LDE5LDE5LDE5LDIyLDQ2LDQ5LDIwLDIwLDIzLDE5LDIyLDIwLDIwLDIwLDUyLDQ2LDQ3LDI0LDE0KQpgYGAKCkV2YWx1YXRlIGRlbnNpdGllcyBpbiBncmlkCgpgYGB7ciB9CnggPC0gc2VxKDAuMDAwMSwgMC45OTk5LCBsZW5ndGgub3V0ID0gMTAwMCkKYGBgCgpIZWxwZXIgZnVuY3Rpb24gdG8gZXZhbHVhdGUgZGVuc2l0eSBvdmVyIG9ic2VydmF0aW9ucwoKYGBge3IgfQpiZGVucyA8LSBmdW5jdGlvbihuLCB5LCB4KQogICAgZGJldGEoeCwgeSsxLCBuLXkrMSkKYGBgCgpTZXBhcmF0ZSBtb2RlbAoKYGBge3IgfQpkZl9zZXAgPC0gbWFwcGx5KGJkZW5zLCBuLCB5LCBNb3JlQXJncyA9IGxpc3QoeCA9IHgpKSAlPiUKICBhcy5kYXRhLmZyYW1lKCkgJT4lIGNiaW5kKHgpICU+JSAKICBwaXZvdF9sb25nZXIoY29scyA9ICF4LCBuYW1lc190byA9ICJpbmQiLCB2YWx1ZXNfdG8gPSAicCIpCmBgYAoKUGxvdCB0aGUgc2VwYXJhdGUgbW9kZWwKCmBgYHtyIH0KbGFiczEgPC0gcGFzdGUoJ3Bvc3RlcmlvciBvZicsIGMoJ3RoZXRhX2onLCAndGhldGFfNzEnKSkKcGxvdF9zZXAgPC0gZ2dwbG90KGRhdGEgPSBkZl9zZXApICsKICBnZW9tX2xpbmUoYWVzKHggPSB4LCB5ID0gcCwgY29sb3IgPSAoaW5kPT0nVjcxJyksIGdyb3VwID0gaW5kKSkgKwogIGxhYnMoeCA9IGV4cHJlc3Npb24odGhldGEpLCB5ID0gJycsIHRpdGxlID0gJ1NlcGFyYXRlIG1vZGVsJywgY29sb3IgPSAnJykgKwogIHNjYWxlX3lfY29udGludW91cyhicmVha3MgPSBOVUxMKSArCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGMoJ2JsdWUnLCdyZWQnKSwgbGFiZWxzID0gbGFiczEpICsKICB0aGVtZShsZWdlbmQuYmFja2dyb3VuZCA9IGVsZW1lbnRfYmxhbmsoKSwgbGVnZW5kLnBvc2l0aW9uID0gYygwLjgsMC45KSkKIyBUaGUgbGFzdCBvbmUgaXMgZm9yIGVtcGhhc2l6ZSBjb2xvcmVkIHJlZApwbG90X3NlcApgYGAKClBvb2xlZCBtb2RlbAoKYGBge3IgfQpkZl9wb29sIDwtIGRhdGEuZnJhbWUoeCA9IHgsIHAgPSBkYmV0YSh4LCBzdW0oeSkrMSwgc3VtKG4pLXN1bSh5KSsxKSkKYGBgCgpDcmVhdGUgYSBwbG90IGZvciB0aGUgc2VwYXJhdGUgbW9kZWwKCmBgYHtyIH0KcGxvdF9wb29sIDwtIGdncGxvdChkYXRhID0gZGZfcG9vbCkgKwogIGdlb21fbGluZShhZXMoeCA9IHgsIHkgPSBwLCBjb2xvciA9ICcxJykpICsKICBsYWJzKHggPSBleHByZXNzaW9uKHRoZXRhKSwgeSA9ICcnLCB0aXRsZSA9ICdQb29sZWQgbW9kZWwnLCBjb2xvciA9ICcnKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcyA9IE5VTEwpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gJ3JlZCcsIGxhYmVscyA9ICdQb3N0ZXJpb3Igb2YgY29tbW9uIHRoZXRhJykgKwogIHRoZW1lKGxlZ2VuZC5iYWNrZ3JvdW5kID0gZWxlbWVudF9ibGFuaygpLCBsZWdlbmQucG9zaXRpb24gPSBjKDAuNywwLjkpKQpgYGAKClBsb3QgYm90aCBzZXBhcmF0ZSBhbmQgcG9vbGVkIG1vZGVsCgpgYGB7ciB9CmdyaWQuYXJyYW5nZShwbG90X3NlcCwgcGxvdF9wb29sKQpgYGAKCkNvbXB1dGUgdGhlIG1hcmdpbmFsIHBvc3RlcmlvciBvZiBhbHBoYSBhbmQgYmV0YSBpbiBoaWVyYXJjaGljYWwgbW9kZWwKVXNlIGdyaWQKCmBgYHtyIH0KQSA8LSBzZXEoMC41LCA2LCBsZW5ndGgub3V0ID0gMTAwKQpCIDwtIHNlcSgzLCAzMywgbGVuZ3RoLm91dCA9IDEwMCkKYGBgCgpNYWtlIHZlY3RvcnMgdGhhdCBjb250YWluIGFsbCBwYWlyd2lzZSBjb21iaW5hdGlvbnMgb2YgQSBhbmQgQgoKYGBge3IgfQpjQSA8LSByZXAoQSwgZWFjaCA9IGxlbmd0aChCKSkKY0IgPC0gcmVwKEIsIGxlbmd0aChBKSkKYGBgCgpVc2UgbG9nYXJpdGhtcyBmb3IgbnVtZXJpY2FsIGFjY3VyYWN5IQoKYGBge3IgfQpscGZ1biA8LSBmdW5jdGlvbihhLCBiLCB5LCBuKSBsb2coYStiKSooLTUvMikgKwogIHN1bShsZ2FtbWEoYStiKS1sZ2FtbWEoYSktbGdhbW1hKGIpK2xnYW1tYShhK3kpK2xnYW1tYShiK24teSktbGdhbW1hKGErYituKSkKbHAgPC0gbWFwcGx5KGxwZnVuLCBjQSwgY0IsIE1vcmVBcmdzID0gbGlzdCh5LCBuKSkKYGBgCgpTdWJ0cmFjdCBtYXhpbXVtIHZhbHVlIHRvIGF2b2lkIG92ZXIvdW5kZXJmbG93IGluIGV4cG9uZW50YXRpb24KCmBgYHtyIH0KZGZfbWFyZyA8LSBkYXRhLmZyYW1lKHggPSBjQSwgeSA9IGNCLCBwID0gZXhwKGxwIC0gbWF4KGxwKSkpCmBgYAoKQ3JlYXRlIGEgcGxvdCBvZiB0aGUgbWFyZ2luYWwgcG9zdGVyaW9yIGRlbnNpdHkKCmBgYHtyIH0KdGl0bGUxIDwtIFRlWCgnVGhlIG1hcmdpbmFsIG9mICRcXGFscGhhJCBhbmQgJFxcYmV0YSQnKQpnZ3Bsb3QoZGF0YSA9IGRmX21hcmcsIGFlcyh4ID0geCwgeSA9IHkpKSArCiAgZ2VvbV9yYXN0ZXIoYWVzKGZpbGwgPSBwLCBhbHBoYSA9IHApLCBpbnRlcnBvbGF0ZSA9IFQpICsKICBnZW9tX2NvbnRvdXIoYWVzKHogPSBwKSwgY29sb3VyID0gJ2JsYWNrJywgc2l6ZSA9IDAuMikgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygxLDUpLCB5bGltID0gYyg0LCAyNikpICsKICBsYWJzKHggPSBUZVgoJyRcXGFscGhhJCcpLCB5ID0gVGVYKCckXFxiZXRhJCcpLCB0aXRsZSA9IHRpdGxlMSkgKwogIHNjYWxlX2ZpbGxfZ3JhZGllbnQobG93ID0gJ3llbGxvdycsIGhpZ2ggPSAncmVkJywgZ3VpZGUgPSBGKSArCiAgc2NhbGVfYWxwaGEocmFuZ2UgPSBjKDAsIDEpLCBndWlkZSA9IEYpCmBgYAoKU2FtcGxlIGZyb20gdGhlIGdyaWQgKHdpdGggcmVwbGFjZW1lbnQpCgpgYGB7ciB9Cm5zYW1wIDwtIDEwMApzYW1wX2luZGljZXMgPC0gc2FtcGxlKGxlbmd0aChkZl9tYXJnJHApLCBzaXplID0gbnNhbXAsCiAgICAgICAgICAgICAgICAgICAgICAgcmVwbGFjZSA9IFQsIHByb2IgPSBkZl9tYXJnJHAvc3VtKGRmX21hcmckcCkpCnNhbXBfQSA8LSBjQVtzYW1wX2luZGljZXNbMTpuc2FtcF1dCnNhbXBfQiA8LSBjQltzYW1wX2luZGljZXNbMTpuc2FtcF1dCmRmX3BzYW1wIDwtIG1hcHBseShmdW5jdGlvbihhLCBiLCB4KSBkYmV0YSh4LCBhLCBiKSwKICAgICAgICAgICAgICAgICAgc2FtcF9BLCBzYW1wX0IsIE1vcmVBcmdzID0gbGlzdCh4ID0geCkpICU+JQogIGFzLmRhdGEuZnJhbWUoKSAlPiUgY2JpbmQoeCkgJT4lIAogIHBpdm90X2xvbmdlcihjb2xzID0gIXgsIG5hbWVzX3RvID0gImluZCIsIHZhbHVlc190byA9ICJwIikKYGBgCgpDcmVhdGUgcGxvdCBmb3Igc2FtcGxlcyBmcm9tIHRoZSBkaXN0cmlidXRpb24gb2YgZGlzdHJpYnV0aW9ucwpCZXRhKGFscGhhLGJldGEpLCB0aGF0IGlzLCBwbG90IEJldGEoYWxwaGEsYmV0YSkgdXNpbmcgcG9zdGVyaW9yCnNhbXBsZXMgb2YgYWxwaGEgYW5kIGJldGEKCmBgYHtyIH0KIyBoZWxwZXIgZnVuY3Rpb24gdG8gY29udmVydCBpbmQgdG8gbnVtZXJpYyBmb3Igc3Vic2V0dGluZwppbmR0b251bSA8LSBmdW5jdGlvbih4KSBzdHJ0b2koc3Vic3RyaW5nKHgsMikpCnRpdGxlMiA8LSBUZVgoJ0JldGEoJFxcYWxwaGEsXFxiZXRhJCkgZ2l2ZW4gcG9zdGVyaW9yIGRyYXdzIG9mICRcXGFscGhhJCBhbmQgJFxcYmV0YSQnKQpwbG90X3BzYW1wIDwtIGdncGxvdChkYXRhID0gc3Vic2V0KGRmX3BzYW1wLCBpbmR0b251bShpbmQpIDw9IDIwKSkgKwogIGdlb21fbGluZShhZXMoeCA9IHgsIHkgPSBwLCBncm91cCA9IGluZCksIGNvbG9yPSdmb3Jlc3RncmVlbicpICsKICBsYWJzKHggPSBleHByZXNzaW9uKHRoZXRhKSwgeSA9ICcnLCB0aXRsZSA9IHRpdGxlMikgKwogIHNjYWxlX3lfY29udGludW91cyhicmVha3MgPSBOVUxMKQpgYGAKClRoZSBhdmVyYWdlIG9mIGFib3ZlIGRpc3RyaWJ1dGlvbnMsIGlzIHRoZSBwcmVkaWN0aXZlIGRpc3RyaWJ1dGlvbgpmb3IgYSBuZXcgdGhldGEsIGFuZCBhbHNvIHRoZSBwcmlvciBkaXN0cmlidXRpb24gZm9yIHRoZXRhX2oKCmBgYHtyIH0KZGZfcHNhbXBtZWFuIDwtIHNwcmVhZChkZl9wc2FtcCwgaW5kLCBwKSAlPiUgc3Vic2V0KHNlbGVjdCA9IC14KSAlPiUKICAgIHJvd01lYW5zKCkgJT4lIGRhdGEuZnJhbWUoeCA9IHgsIHAgPSAuKQpgYGAKCkNyZWF0ZSBwbG90IGZvciBzYW1wbGVzIGZyb20gdGhlIHByZWRpY3RpdmUgZGlzdHJpYnV0aW9uIGZvciBuZXcgdGhldGEKCmBgYHtyIH0KdGl0bGUzIDwtIFRlWCgnUG9wdWxhdGlvbiBkaXN0cmlidXRpb24gKHByaW9yKSBmb3IgJFxcdGhldGFfaiQnKQpwbG90X3BzYW1wbWVhbiA8LSBnZ3Bsb3QoZGF0YSA9IGRmX3BzYW1wbWVhbikgKwogIGdlb21fbGluZShhZXMoeCA9IHgsIHkgPSBwKSwgY29sb3I9J2ZvcmVzdGdyZWVuJykgKwogIGxhYnMoeCA9IGV4cHJlc3Npb24odGhldGEpLCB5ID0gJycsIHRpdGxlID0gdGl0bGUzKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcyA9IE5VTEwpCmBgYAoKQ29tYmluZSB0aGUgcGxvdHMKCmBgYHtyIH0KZ3JpZC5hcnJhbmdlKHBsb3RfcHNhbXAsIHBsb3RfcHNhbXBtZWFuKQpgYGAKCkFuZCBmaW5hbGx5IGNvbXBhcmUgdGhlIHNlcGFyYXRlIG1vZGVsIGFuZCBoaWVyYXJjaGljYWwgbW9kZWwKKHVzaW5nIGV2ZXJ5IHNldmVudGggc2FtcGxlIGZvciBjbGFyaXR5KQpDcmVhdGUgcGxvdCBmb3IgdGhlIHNlcGFyYXRlIG1vZGVsCgpgYGB7ciB9CnBsb3Rfc2VwNyA8LSBnZ3Bsb3QoZGF0YSA9IHN1YnNldChkZl9zZXAsIGluZHRvbnVtKGluZCklJTc9PTApKSArCiAgZ2VvbV9saW5lKGFlcyh4ID0geCwgeSA9IHAsIGNvbG9yID0gKGluZD09J1Y0OScpLCBncm91cCA9IGluZCkpICsKICBsYWJzKHggPSBleHByZXNzaW9uKHRoZXRhKSwgeSA9ICcnLCB0aXRsZSA9ICdTZXBhcmF0ZSBtb2RlbCcsIGNvbG9yID0gJycpICsKICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzID0gTlVMTCkgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCdibHVlJywgJ3JlZCcpLCBndWlkZSA9IEYpICsKICB0aGVtZShsZWdlbmQuYmFja2dyb3VuZCA9IGVsZW1lbnRfYmxhbmsoKSwgbGVnZW5kLnBvc2l0aW9uID0gYygwLjgsMC45KSkKYGBgCgpUaGUgaGllcmFyY2hpY2FsIG1vZGVsCgpOb3RlIHRoYXQgdGhlc2UgbWFyZ2luYWwgcG9zdGVyaW9ycyBmb3IgdGhldGFfaiBhcmUgbW9yZSBuYXJyb3cgdGhhbgppbiB0aGUgc2VwYXJhdGUgbW9kZWwgY2FzZSwgZHVlIHRvIHRoZSBib3Jyb3dlZCBpbmZvcm1hdGlvbiBmcm9tIHRoZQpvdGhlciB0aGV0YV9qJ3MKCkF2ZXJhZ2UgZGVuc2l0eSBvdmVyIHNhbXBsZXMgKG9mIGEgYW5kIGIpIGZvciBlYWNoIChuLHkpLXBhaXIKYXQgZWFjaCBwb2ludCB4CgpgYGB7ciB9CmJkZW5zMiA8LSBmdW5jdGlvbihuLCB5LCBhLCBiLCB4KQogIHJvd01lYW5zKG1hcHBseShkYmV0YSwgYSArIHksIG4gLSB5ICsgYiwgTW9yZUFyZ3MgPSBsaXN0KHggPSB4KSkpCmRmX2hpZXIgPC0gbWFwcGx5KGJkZW5zMiwgbiwgeSwgTW9yZUFyZ3MgPSBsaXN0KHNhbXBfQSwgc2FtcF9CLCB4KSkgJT4lCiAgYXMuZGF0YS5mcmFtZSgpICU+JSBjYmluZCh4KSAlPiUgCiAgcGl2b3RfbG9uZ2VyKGNvbHMgPSAheCwgbmFtZXNfdG8gPSAiaW5kIiwgdmFsdWVzX3RvID0gInAiKQpgYGAKCkNyZWF0ZSBwbG90IGZvciB0aGUgaGllcmFyY2hpY2FsIG1vZGVsCgpgYGB7ciB9CnBsb3RfaGllcjcgPC0gZ2dwbG90KGRhdGEgPSBzdWJzZXQoZGZfaGllciwgaW5kdG9udW0oaW5kKSUlNz09MCkpICsKICBnZW9tX2xpbmUoYWVzKHggPSB4LCB5ID0gcCwgY29sb3IgPSAoaW5kPT0nVjQ5JyksIGdyb3VwID0gaW5kKSkgKwogIGxhYnMoeCA9IGV4cHJlc3Npb24odGhldGEpLCB5ID0gJycsIHRpdGxlID0gJ0hpZXJhcmNoaWNhbCBtb2RlbCcsIGNvbG9yID0gJycpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gYygnYmx1ZScsICdyZWQnKSwgZ3VpZGUgPSBGKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcyA9IE5VTEwpICsKICB0aGVtZShsZWdlbmQuYmFja2dyb3VuZCA9IGVsZW1lbnRfYmxhbmsoKSwgbGVnZW5kLnBvc2l0aW9uID0gYygwLjgsMC45KSkKYGBgCgpDb21iaW5lIHRoZSBwbG90cwoKYGBge3IgfQpncmlkLmFycmFuZ2UocGxvdF9zZXA3LCBwbG90X2hpZXI3KQpgYGAKCjkwJSBpbnRlcnZhbHMgZm9yIHNlcGFyYXRlIG1vZGVscwoKYGBge3IgfQpxcV9zZXBhcmF0ZTwtZGF0YS5mcmFtZShpZD0xOmxlbmd0aChuKSwgbj1uLCB5PXksCiAgICAgICAgICAgICAgIHExMD1xYmV0YSgwLjEseSsxLG4teSsxKSwgcTUwPXFiZXRhKDAuNSx5KzEsbi15KzEpLCBxOTA9cWJldGEoMC45LHkrMSxuLXkrMSkpCmBgYAoKOTAlIGludGVydmFscyBmb3IgaGllcmFyY2hpY2FsIG1vZGVsCgpgYGB7ciB9CnFoIDwtIGZ1bmN0aW9uKHEsIG4sIHkpIGNvbE1lYW5zKG1hcHBseShmdW5jdGlvbihxLCBuLCB5LCBhLCBiKQogICAgbWFwcGx5KHFiZXRhLCBxLCBhICsgeSwgbiAtIHkgKyBiKSwgcSwgbiwgeSwgTW9yZUFyZ3MgPSBsaXN0KHNhbXBfQSwgc2FtcF9CKSkpCnFxX2hpZXIgPC0gZGF0YS5mcmFtZShpZD0xOmxlbmd0aChuKSwgbj1uLCB5PXksCiAgICAgICAgICAgICAgICAgICAgICBxMTAgPSBxaCgwLjEsIG4sIHkpLCBxNTAgPSBxaCgwLjUsIG4sIHkpLCBxOTAgPSBxaCgwLjksIG4sIHkpKQpgYGAKCnBvb2xlZCBtZWFuCgpgYGB7ciB9Cm1fcG9vbGVkIDwtIChzdW0oeSkrMSkvKHN1bShuKSsyKQpgYGAKCnBsb3QKCmBgYHtyIH0KcGxvdF9zZXBfaW50IDwtIHFxX3NlcGFyYXRlICU+JQogIGdncGxvdChhZXMoeD1uLHk9cTUwLHltaW49cTEwLHltYXg9cTkwKSkgKwogIGdlb21fcG9pbnRyYW5nZShjb2xvcj0iYmx1ZSIsIGFscGhhPTAuMikgKwogIGdlb21faGxpbmUoeWludGVyY2VwdD1tX3Bvb2xlZCwgbGluZXR5cGU9ImRhc2hlZCIpKwogIGxpbXMoeD1jKDAsNjApLCB5PWMoMCwwLjUxKSkrCiAgbGFicyh4PSJOIiwgeT0iUHJvYmFiaWxpdHkgb2YgcG9seXBzIiwgdGl0bGU9IlNlcGFyYXRlIikrCiAgYW5ub3RhdGUoInRleHQiLCB4PTAsIHk9bV9wb29sZWQsIGxhYmVsID0gIlBvb2xlZCBtZWFuIiwgdmp1c3Q9LTAuMiwgaGp1c3Q9MC4zKQpwbG90X2hpZXJfaW50IDwtIHFxX2hpZXIgJT4lCiAgZ2dwbG90KGFlcyh4PW4seT1xNTAseW1pbj1xMTAseW1heD1xOTApKSArCiAgZ2VvbV9wb2ludHJhbmdlKGNvbG9yPSJibHVlIiwgYWxwaGE9MC4yKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0PW1fcG9vbGVkLCBsaW5ldHlwZT0iZGFzaGVkIikrCiAgbGltcyh4PWMoMCw2MCksIHk9YygwLDAuNTEpKSsKICBsYWJzKHg9Ik4iLCB5PSJQcm9iYWJpbGl0eSBvZiBwb2x5cHMiLCB0aXRsZT0iSGllcmFyY2hpY2FsIikrCiAgYW5ub3RhdGUoInRleHQiLCB4PTAsIHk9bV9wb29sZWQsIGxhYmVsID0gIlBvb2xlZCBtZWFuIiwgdmp1c3Q9LTAuMiwgaGp1c3Q9MC4zKQpncmlkLmFycmFuZ2UocGxvdF9zZXBfaW50LCBwbG90X2hpZXJfaW50KQpgYGAKCg==