Code and figures for ChileSchools example. See Chapter 21 in Regression and Other Stories.

Data are from

  • Chay, K. Y., McEwan, P. J., and Urquiola, M. (2005). The central role of noise in evaluating interventions that use test scores to rank schools. American Economic Review 95, 1237–1258.

Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("foreign")
library("arm")
library("rstanarm")
library("brms")
library("ivpack")

Load data

The outcomes in these analyses are the gain scores between 88 and 92.

chile <- read.csv(root("ChileSchools/data","chile.csv"))
print(head(chile), digits=3)
  p90       cmb_regn urban88  rule2 cutoff cutoff_cmb eligible read92 read88
1   0 Regions 2,5,10       1  3.415   49.4       49.4        0   57.0   54.4
2   0    Regions 6,8       1 31.815   43.4       43.4        0   85.5   79.1
3   0    Regions 6,8       1  4.445   43.4       43.4        0   52.0   47.8
4   0      Region 13       1 16.800   46.4       46.4        0   66.4   65.1
5   0 Regions 2,5,10       1  0.095   49.4       49.4        0   52.5   49.3
6   1      Region 13       1 -2.030   46.4       46.4        1   55.3   50.5
  math92 math88
1   67.1   51.3
2   84.7   71.4
3   53.3   47.9
4   56.4   61.3
5   50.9   49.7
6   59.0   38.2

Fit models

fit_1a <- stan_glm(read92 ~ eligible + rule2, data=chile, refresh=0)
fit_1b <- stan_glm(read92 ~ eligible + rule2, data=chile, subset = abs(rule2)<5, refresh=0)

fit_2b <- stan_glm(read92 ~ eligible + rule2 + eligible:rule2, data=chile, subset = abs(rule2)<5, refresh=0)

fit_3a <- stan_glm(read92 ~ eligible + rule2 + read88 + math88, data=chile, refresh=0)
fit_3b <- stan_glm(read92 ~ eligible + rule2 + read88 + math88, data=chile, subset = abs(rule2)<5, refresh=0)

chile$z_read88 <- (chile$read88 - mean(chile$read88))/sd(chile$read88)
chile$z_math88 <- (chile$math88 - mean(chile$math88))/sd(chile$math88)
fit_5a <- stan_glm(read92 ~ eligible + rule2 + z_read88 + z_math88 + eligible:z_read88, data=chile, refresh=0)
fit_5b <- stan_glm(read92 ~ eligible + rule2 + z_read88 + z_math88 + eligible:z_read88, data=chile, subset = abs(rule2)<5, refresh=0)

Plots

whiteline <- function(a, b, from=-1000, to=1000){
  curve(a + b*x, from=from, to=to, col="white", lwd=4, add=TRUE)
  curve(a + b*x, from=from, to=to, lwd=2, add=TRUE)
}
par(mar=c(3,3,2,1), mgp=c(1.7, .5, 0), tck=-.01)
yy <- chile$read92
xx <- chile$rule2
plot(y=yy, x=xx, type="n", 
     xlab="Assignment variable",
     ylab="Post-test", bty="l", yaxt="n")
axis(2, seq(0,100,20))
points(xx[chile$eligible==0],yy[chile$eligible==0],pch=20,
       cex=.7,col="gray60")
points(xx[chile$eligible==1],yy[chile$eligible==1],pch=20,
       cex=.7)
whiteline(coef(fit_1a)[1], coef(fit_1a)[3], from=0)
whiteline(coef(fit_1a)[1] + coef(fit_1a)[2], coef(fit_1a)[3], to=0)
abline(v=0, lwd=4, col="gray")
mtext("All the data", 3, 1)

par(mar=c(3,3,2,1), mgp=c(1.7, .5, 0), tck=-.01)
yy <- chile$read92
xx <- chile$rule2
plot(y=yy, x=xx, type="n", 
     xlab="Assignment variable",
     ylab="Post-test", bty="l", yaxt="n", xlim=c(-5,5), xaxs="i")
axis(2, seq(0,100,20))
points(xx[chile$eligible==0],yy[chile$eligible==0],pch=20,
       cex=.7,col="gray60")
points(xx[chile$eligible==1],yy[chile$eligible==1],pch=20,
       cex=.7)
whiteline(coef(fit_1b)[1], coef(fit_1b)[3], from=0)
whiteline(coef(fit_1b)[1] + coef(fit_1b)[2], coef(fit_1b)[3], to=0)
abline(v=0, lwd=4, col="gray")
mtext("Restricting to data near the cutoff", 3, 1)

par(mar=c(3,3,2,1), mgp=c(1.7, .5, 0), tck=-.01)
yy <- chile$read92
xx <- chile$rule2
plot(y=yy, x=xx, type="n", 
     xlab="Assignment variable",
     ylab="Ppst-test", bty="l", yaxt="n", xlim=c(-5,5), xaxs="i")
axis(2, seq(0,100,20))
points(xx[chile$eligible==0],yy[chile$eligible==0],pch=20,
       cex=.7,col="gray60")
points(xx[chile$eligible==1],yy[chile$eligible==1],pch=20,
       cex=.7)
whiteline(coef(fit_2b)[1], coef(fit_2b)[3], from=0)
whiteline(coef(fit_2b)[1] + coef(fit_2b)[2], coef(fit_2b)[3] + coef(fit_2b)[4], to=0)
abline(v=0, lwd=4, col="gray")
mtext("Model with interaction", 3, 1)

par(mar=c(3,3,2,1), mgp=c(1.7, .5, 0), tck=-.01)
yy <- chile$read92 - coef(fit_3b)[4] * (chile$read88 - mean(chile$read88)) - coef(fit_3b)[5] * (chile$math88 - mean(chile$math88))
xx <- chile$rule2
plot(y=yy, x=xx, type="n", 
     xlab="Assignment variable",
     ylab="Adjusted outcome", bty="l", yaxt="n", xlim=c(-5,5), xaxs="i")
axis(2, seq(-100,100,20))
points(xx[chile$eligible==0],yy[chile$eligible==0],pch=20,
       cex=.7,col="gray60")
points(xx[chile$eligible==1],yy[chile$eligible==1],pch=20,
       cex=.7)
whiteline(coef(fit_3b)[1] + coef(fit_3b)[4] *mean(chile$read88) + coef(fit_3b)[5] *mean(chile$math88), coef(fit_3b)[3], from=0)
whiteline(coef(fit_3b)[1] + coef(fit_3b)[4] *mean(chile$read88) + coef(fit_3b)[5] *mean(chile$math88) + coef(fit_3b)[2], coef(fit_3b)[3], to=0)
abline(v=0, lwd=4, col="gray")
mtext("Adjusting for pre-test scores", 3, 1)

par(mar=c(3,3,2,1), mgp=c(1.7, .5, 0), tck=-.01)
yy <- chile$read92 - coef(fit_3b)[4] * (chile$read88 - mean(chile$read88)) - coef(fit_3b)[5] * (chile$math88 - mean(chile$math88))
xx <- chile$rule2
n_bins <- 20
n <- length(xx)
halfwidth <- 5
cutpoints <- c(seq(-halfwidth, halfwidth, length=n_bins+1))
xx_bin <- rep(NA, n_bins)
yy_bin <- rep(NA, n_bins)
for (i in 1:n_bins){
  keep <- xx > cutpoints[i] & xx <= cutpoints[i+1]
  xx_bin[i] <- mean(xx[keep])
  yy_bin[i] <- mean(yy[keep])
}
plot(y=yy_bin, x=xx_bin, type="n", 
     xlab="Assignment variable",
     ylab="Adjusted outcome", bty="l", yaxt="n", xlim=c(-halfwidth, halfwidth), xaxs="i")
axis(2, seq(-100,100,1))
points(xx_bin[xx_bin>0], yy_bin[xx_bin>0],pch=20,
       cex=2, col="gray60")
points(xx_bin[xx_bin<0], yy_bin[xx_bin<0],pch=20,
       cex=2)
whiteline(coef(fit_3b)[1] + coef(fit_3b)[4] *mean(chile$read88) + coef(fit_3b)[5] *mean(chile$math88), coef(fit_3b)[3], from=0)
whiteline(coef(fit_3b)[1] + coef(fit_3b)[4] *mean(chile$read88) + coef(fit_3b)[5] *mean(chile$math88) + coef(fit_3b)[2], coef(fit_3b)[3], to=0)
abline(v=0, lwd=4, col="gray")
mtext("Binned averages with same regression lines", 3, 1)

Additional models

Noncompliance rates

mean(chile$p90[chile$eligible==0 & abs(chile$rule2) < 5])  # .053
[1] 0.05333333
mean(chile$p90[chile$eligible==1 & abs(chile$rule2) < 5]) # .606
[1] 0.6061453

IV model on restricted dataset

# with brms as in IV section earlier in the chapter
set.seed(1234)
chile$diff_read92 <- chile$read92 - chile$read88
rd_f1 <- bf(p90 ~ eligible)
rd_f2 <- bf(diff_read92 ~ p90)
IV_brm_rd <- brm(formula=rd_f1 + rd_f2, data = chile[abs(chile$rule2) < 5,])
print(IV_brm_rd, digits=3)
 Family: MV(gaussian, gaussian) 
  Links: mu = identity; sigma = identity
         mu = identity; sigma = identity 
Formula: p90 ~ eligible 
         diff_read92 ~ p90 
   Data: chile[abs(chile$rule2) < 5, ] (Number of observations: 883) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS
p90_Intercept           0.054     0.016    0.023    0.084 1.000     4812
diffread92_Intercept   10.944     0.341   10.271   11.594 1.001     2611
p90_eligible            0.552     0.025    0.502    0.600 1.000     4484
diffread92_p90          5.065     0.876    3.347    6.749 1.002     2329
                     Tail_ESS
p90_Intercept            2873
diffread92_Intercept     2584
p90_eligible             2840
diffread92_p90           2602

Family Specific Parameters: 
                 Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
sigma_p90           0.357     0.009    0.341    0.375 1.001     4718     2808
sigma_diffread92    7.123     0.176    6.786    7.484 1.001     3478     3244

Residual Correlations: 
                       Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS
rescor(p90,diffread92)   -0.182     0.053   -0.282   -0.073 1.001     2524
                       Tail_ESS
rescor(p90,diffread92)     2733

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
# t.e. is 5.08 with s.e. of .90

Compare to IV regression by two-stage least squares

summary(ivreg(formula = diff_read92 ~ p90 | eligible, data=chile, subset = abs(rule2) < 5))

Call:
ivreg(formula = diff_read92 ~ p90 | eligible, data = chile, subset = abs(rule2) < 
    5)

Residuals:
     Min       1Q   Median       3Q      Max 
-33.8471  -4.3416   0.1009   4.6948  25.7337 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.9411     0.3421  31.981  < 2e-16 ***
p90           5.0871     0.8814   5.771 1.09e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.109 on 881 degrees of freedom
Multiple R-Squared: 0.009344,   Adjusted R-squared: 0.00822 
Wald test: 33.31 on 1 and 881 DF,  p-value: 1.088e-08 
# est is 5.09 with s.e. of .88
LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogQ2hpbGVTY2hvb2xzIgphdXRob3I6ICJBbmRyZXcgR2VsbWFuLCBKZW5uaWZlciBIaWxsLCBBa2kgVmVodGFyaSIKZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdGhlbWU6IHJlYWRhYmxlCiAgICB0b2M6IHRydWUKICAgIHRvY19kZXB0aDogMgogICAgdG9jX2Zsb2F0OiB0cnVlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCi0tLQpDb2RlIGFuZCBmaWd1cmVzIGZvciBDaGlsZVNjaG9vbHMgZXhhbXBsZS4gU2VlIENoYXB0ZXIgMjEgaW4KUmVncmVzc2lvbiBhbmQgT3RoZXIgU3Rvcmllcy4KCkRhdGEgYXJlIGZyb20KCi0gQ2hheSwgSy4gWS4sIE1jRXdhbiwgUC4gSi4sIGFuZCBVcnF1aW9sYSwgTS4gKDIwMDUpLiBUaGUgY2VudHJhbAogIHJvbGUgb2Ygbm9pc2UgaW4gZXZhbHVhdGluZyBpbnRlcnZlbnRpb25zIHRoYXQgdXNlIHRlc3Qgc2NvcmVzIHRvCiAgcmFuayBzY2hvb2xzLiBBbWVyaWNhbiBFY29ub21pYyBSZXZpZXcgOTUsIDEyMzfigJMxMjU4LgoKLS0tLS0tLS0tLS0tLQoKCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQobWVzc2FnZT1GQUxTRSwgZXJyb3I9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGNvbW1lbnQ9TkEpCiMgc3dpdGNoIHRoaXMgdG8gVFJVRSB0byBzYXZlIGZpZ3VyZXMgaW4gc2VwYXJhdGUgZmlsZXMKc2F2ZWZpZ3MgPC0gRkFMU0UKYGBgCgojIyMjIExvYWQgcGFja2FnZXMKCmBgYHtyIH0KbGlicmFyeSgicnByb2pyb290IikKcm9vdDwtaGFzX2ZpbGUoIi5ST1MtRXhhbXBsZXMtcm9vdCIpJG1ha2VfZml4X2ZpbGUoKQpsaWJyYXJ5KCJmb3JlaWduIikKbGlicmFyeSgiYXJtIikKbGlicmFyeSgicnN0YW5hcm0iKQpsaWJyYXJ5KCJicm1zIikKbGlicmFyeSgiaXZwYWNrIikKYGBgCgojIyMjIExvYWQgZGF0YQoKVGhlIG91dGNvbWVzIGluIHRoZXNlIGFuYWx5c2VzIGFyZSB0aGUgKmdhaW4gc2NvcmVzKiBiZXR3ZWVuIDg4IGFuZCA5Mi4gCgpgYGB7ciB9CmNoaWxlIDwtIHJlYWQuY3N2KHJvb3QoIkNoaWxlU2Nob29scy9kYXRhIiwiY2hpbGUuY3N2IikpCnByaW50KGhlYWQoY2hpbGUpLCBkaWdpdHM9MykKYGBgCgojIyBGaXQgbW9kZWxzCgpgYGB7ciB9CmZpdF8xYSA8LSBzdGFuX2dsbShyZWFkOTIgfiBlbGlnaWJsZSArIHJ1bGUyLCBkYXRhPWNoaWxlLCByZWZyZXNoPTApCmZpdF8xYiA8LSBzdGFuX2dsbShyZWFkOTIgfiBlbGlnaWJsZSArIHJ1bGUyLCBkYXRhPWNoaWxlLCBzdWJzZXQgPSBhYnMocnVsZTIpPDUsIHJlZnJlc2g9MCkKCmZpdF8yYiA8LSBzdGFuX2dsbShyZWFkOTIgfiBlbGlnaWJsZSArIHJ1bGUyICsgZWxpZ2libGU6cnVsZTIsIGRhdGE9Y2hpbGUsIHN1YnNldCA9IGFicyhydWxlMik8NSwgcmVmcmVzaD0wKQoKZml0XzNhIDwtIHN0YW5fZ2xtKHJlYWQ5MiB+IGVsaWdpYmxlICsgcnVsZTIgKyByZWFkODggKyBtYXRoODgsIGRhdGE9Y2hpbGUsIHJlZnJlc2g9MCkKZml0XzNiIDwtIHN0YW5fZ2xtKHJlYWQ5MiB+IGVsaWdpYmxlICsgcnVsZTIgKyByZWFkODggKyBtYXRoODgsIGRhdGE9Y2hpbGUsIHN1YnNldCA9IGFicyhydWxlMik8NSwgcmVmcmVzaD0wKQoKY2hpbGUkel9yZWFkODggPC0gKGNoaWxlJHJlYWQ4OCAtIG1lYW4oY2hpbGUkcmVhZDg4KSkvc2QoY2hpbGUkcmVhZDg4KQpjaGlsZSR6X21hdGg4OCA8LSAoY2hpbGUkbWF0aDg4IC0gbWVhbihjaGlsZSRtYXRoODgpKS9zZChjaGlsZSRtYXRoODgpCmZpdF81YSA8LSBzdGFuX2dsbShyZWFkOTIgfiBlbGlnaWJsZSArIHJ1bGUyICsgel9yZWFkODggKyB6X21hdGg4OCArIGVsaWdpYmxlOnpfcmVhZDg4LCBkYXRhPWNoaWxlLCByZWZyZXNoPTApCmZpdF81YiA8LSBzdGFuX2dsbShyZWFkOTIgfiBlbGlnaWJsZSArIHJ1bGUyICsgel9yZWFkODggKyB6X21hdGg4OCArIGVsaWdpYmxlOnpfcmVhZDg4LCBkYXRhPWNoaWxlLCBzdWJzZXQgPSBhYnMocnVsZTIpPDUsIHJlZnJlc2g9MCkKYGBgCgojIyBQbG90cwoKYGBge3IgfQp3aGl0ZWxpbmUgPC0gZnVuY3Rpb24oYSwgYiwgZnJvbT0tMTAwMCwgdG89MTAwMCl7CiAgY3VydmUoYSArIGIqeCwgZnJvbT1mcm9tLCB0bz10bywgY29sPSJ3aGl0ZSIsIGx3ZD00LCBhZGQ9VFJVRSkKICBjdXJ2ZShhICsgYip4LCBmcm9tPWZyb20sIHRvPXRvLCBsd2Q9MiwgYWRkPVRSVUUpCn0KCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJDaGlsZVNjaG9vbHMvZmlncyIsInJkX2NoaWxlXzFhLnBkZiIpLCBoZWlnaHQ9My42LCB3aWR0aD01LCBjb2xvcm1vZGVsPSJncmF5IikKYGBgCmBgYHtyIH0KcGFyKG1hcj1jKDMsMywyLDEpLCBtZ3A9YygxLjcsIC41LCAwKSwgdGNrPS0uMDEpCnl5IDwtIGNoaWxlJHJlYWQ5Mgp4eCA8LSBjaGlsZSRydWxlMgpwbG90KHk9eXksIHg9eHgsIHR5cGU9Im4iLCAKICAgICB4bGFiPSJBc3NpZ25tZW50IHZhcmlhYmxlIiwKICAgICB5bGFiPSJQb3N0LXRlc3QiLCBidHk9ImwiLCB5YXh0PSJuIikKYXhpcygyLCBzZXEoMCwxMDAsMjApKQpwb2ludHMoeHhbY2hpbGUkZWxpZ2libGU9PTBdLHl5W2NoaWxlJGVsaWdpYmxlPT0wXSxwY2g9MjAsCiAgICAgICBjZXg9LjcsY29sPSJncmF5NjAiKQpwb2ludHMoeHhbY2hpbGUkZWxpZ2libGU9PTFdLHl5W2NoaWxlJGVsaWdpYmxlPT0xXSxwY2g9MjAsCiAgICAgICBjZXg9LjcpCndoaXRlbGluZShjb2VmKGZpdF8xYSlbMV0sIGNvZWYoZml0XzFhKVszXSwgZnJvbT0wKQp3aGl0ZWxpbmUoY29lZihmaXRfMWEpWzFdICsgY29lZihmaXRfMWEpWzJdLCBjb2VmKGZpdF8xYSlbM10sIHRvPTApCmFibGluZSh2PTAsIGx3ZD00LCBjb2w9ImdyYXkiKQptdGV4dCgiQWxsIHRoZSBkYXRhIiwgMywgMSkKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmRldi5vZmYoKQoKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgcGRmKHJvb3QoIkNoaWxlU2Nob29scy9maWdzIiwicmRfY2hpbGVfMWIucGRmIiksIGhlaWdodD0zLjYsIHdpZHRoPTUsIGNvbG9ybW9kZWw9ImdyYXkiKQpgYGAKYGBge3IgfQpwYXIobWFyPWMoMywzLDIsMSksIG1ncD1jKDEuNywgLjUsIDApLCB0Y2s9LS4wMSkKeXkgPC0gY2hpbGUkcmVhZDkyCnh4IDwtIGNoaWxlJHJ1bGUyCnBsb3QoeT15eSwgeD14eCwgdHlwZT0ibiIsIAogICAgIHhsYWI9IkFzc2lnbm1lbnQgdmFyaWFibGUiLAogICAgIHlsYWI9IlBvc3QtdGVzdCIsIGJ0eT0ibCIsIHlheHQ9Im4iLCB4bGltPWMoLTUsNSksIHhheHM9ImkiKQpheGlzKDIsIHNlcSgwLDEwMCwyMCkpCnBvaW50cyh4eFtjaGlsZSRlbGlnaWJsZT09MF0seXlbY2hpbGUkZWxpZ2libGU9PTBdLHBjaD0yMCwKICAgICAgIGNleD0uNyxjb2w9ImdyYXk2MCIpCnBvaW50cyh4eFtjaGlsZSRlbGlnaWJsZT09MV0seXlbY2hpbGUkZWxpZ2libGU9PTFdLHBjaD0yMCwKICAgICAgIGNleD0uNykKd2hpdGVsaW5lKGNvZWYoZml0XzFiKVsxXSwgY29lZihmaXRfMWIpWzNdLCBmcm9tPTApCndoaXRlbGluZShjb2VmKGZpdF8xYilbMV0gKyBjb2VmKGZpdF8xYilbMl0sIGNvZWYoZml0XzFiKVszXSwgdG89MCkKYWJsaW5lKHY9MCwgbHdkPTQsIGNvbD0iZ3JheSIpCm10ZXh0KCJSZXN0cmljdGluZyB0byBkYXRhIG5lYXIgdGhlIGN1dG9mZiIsIDMsIDEpCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQpkZXYub2ZmKCkKCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJDaGlsZVNjaG9vbHMvZmlncyIsInJkX2NoaWxlXzJiLnBkZiIpLCBoZWlnaHQ9My42LCB3aWR0aD01LCBjb2xvcm1vZGVsPSJncmF5IikKYGBgCmBgYHtyIH0KcGFyKG1hcj1jKDMsMywyLDEpLCBtZ3A9YygxLjcsIC41LCAwKSwgdGNrPS0uMDEpCnl5IDwtIGNoaWxlJHJlYWQ5Mgp4eCA8LSBjaGlsZSRydWxlMgpwbG90KHk9eXksIHg9eHgsIHR5cGU9Im4iLCAKICAgICB4bGFiPSJBc3NpZ25tZW50IHZhcmlhYmxlIiwKICAgICB5bGFiPSJQcHN0LXRlc3QiLCBidHk9ImwiLCB5YXh0PSJuIiwgeGxpbT1jKC01LDUpLCB4YXhzPSJpIikKYXhpcygyLCBzZXEoMCwxMDAsMjApKQpwb2ludHMoeHhbY2hpbGUkZWxpZ2libGU9PTBdLHl5W2NoaWxlJGVsaWdpYmxlPT0wXSxwY2g9MjAsCiAgICAgICBjZXg9LjcsY29sPSJncmF5NjAiKQpwb2ludHMoeHhbY2hpbGUkZWxpZ2libGU9PTFdLHl5W2NoaWxlJGVsaWdpYmxlPT0xXSxwY2g9MjAsCiAgICAgICBjZXg9LjcpCndoaXRlbGluZShjb2VmKGZpdF8yYilbMV0sIGNvZWYoZml0XzJiKVszXSwgZnJvbT0wKQp3aGl0ZWxpbmUoY29lZihmaXRfMmIpWzFdICsgY29lZihmaXRfMmIpWzJdLCBjb2VmKGZpdF8yYilbM10gKyBjb2VmKGZpdF8yYilbNF0sIHRvPTApCmFibGluZSh2PTAsIGx3ZD00LCBjb2w9ImdyYXkiKQptdGV4dCgiTW9kZWwgd2l0aCBpbnRlcmFjdGlvbiIsIDMsIDEpCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQpkZXYub2ZmKCkKCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJDaGlsZVNjaG9vbHMvZmlncyIsInJkX2NoaWxlXzNiLnBkZiIpLCBoZWlnaHQ9My42LCB3aWR0aD01LCBjb2xvcm1vZGVsPSJncmF5IikKYGBgCmBgYHtyIH0KcGFyKG1hcj1jKDMsMywyLDEpLCBtZ3A9YygxLjcsIC41LCAwKSwgdGNrPS0uMDEpCnl5IDwtIGNoaWxlJHJlYWQ5MiAtIGNvZWYoZml0XzNiKVs0XSAqIChjaGlsZSRyZWFkODggLSBtZWFuKGNoaWxlJHJlYWQ4OCkpIC0gY29lZihmaXRfM2IpWzVdICogKGNoaWxlJG1hdGg4OCAtIG1lYW4oY2hpbGUkbWF0aDg4KSkKeHggPC0gY2hpbGUkcnVsZTIKcGxvdCh5PXl5LCB4PXh4LCB0eXBlPSJuIiwgCiAgICAgeGxhYj0iQXNzaWdubWVudCB2YXJpYWJsZSIsCiAgICAgeWxhYj0iQWRqdXN0ZWQgb3V0Y29tZSIsIGJ0eT0ibCIsIHlheHQ9Im4iLCB4bGltPWMoLTUsNSksIHhheHM9ImkiKQpheGlzKDIsIHNlcSgtMTAwLDEwMCwyMCkpCnBvaW50cyh4eFtjaGlsZSRlbGlnaWJsZT09MF0seXlbY2hpbGUkZWxpZ2libGU9PTBdLHBjaD0yMCwKICAgICAgIGNleD0uNyxjb2w9ImdyYXk2MCIpCnBvaW50cyh4eFtjaGlsZSRlbGlnaWJsZT09MV0seXlbY2hpbGUkZWxpZ2libGU9PTFdLHBjaD0yMCwKICAgICAgIGNleD0uNykKd2hpdGVsaW5lKGNvZWYoZml0XzNiKVsxXSArIGNvZWYoZml0XzNiKVs0XSAqbWVhbihjaGlsZSRyZWFkODgpICsgY29lZihmaXRfM2IpWzVdICptZWFuKGNoaWxlJG1hdGg4OCksIGNvZWYoZml0XzNiKVszXSwgZnJvbT0wKQp3aGl0ZWxpbmUoY29lZihmaXRfM2IpWzFdICsgY29lZihmaXRfM2IpWzRdICptZWFuKGNoaWxlJHJlYWQ4OCkgKyBjb2VmKGZpdF8zYilbNV0gKm1lYW4oY2hpbGUkbWF0aDg4KSArIGNvZWYoZml0XzNiKVsyXSwgY29lZihmaXRfM2IpWzNdLCB0bz0wKQphYmxpbmUodj0wLCBsd2Q9NCwgY29sPSJncmF5IikKbXRleHQoIkFkanVzdGluZyBmb3IgcHJlLXRlc3Qgc2NvcmVzIiwgMywgMSkKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmRldi5vZmYoKQoKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgcGRmKHJvb3QoIkNoaWxlU2Nob29scy9maWdzIiwicmRfY2hpbGVfNGIucGRmIiksIGhlaWdodD0zLjYsIHdpZHRoPTUsIGNvbG9ybW9kZWw9ImdyYXkiKQpgYGAKYGBge3IgfQpwYXIobWFyPWMoMywzLDIsMSksIG1ncD1jKDEuNywgLjUsIDApLCB0Y2s9LS4wMSkKeXkgPC0gY2hpbGUkcmVhZDkyIC0gY29lZihmaXRfM2IpWzRdICogKGNoaWxlJHJlYWQ4OCAtIG1lYW4oY2hpbGUkcmVhZDg4KSkgLSBjb2VmKGZpdF8zYilbNV0gKiAoY2hpbGUkbWF0aDg4IC0gbWVhbihjaGlsZSRtYXRoODgpKQp4eCA8LSBjaGlsZSRydWxlMgpuX2JpbnMgPC0gMjAKbiA8LSBsZW5ndGgoeHgpCmhhbGZ3aWR0aCA8LSA1CmN1dHBvaW50cyA8LSBjKHNlcSgtaGFsZndpZHRoLCBoYWxmd2lkdGgsIGxlbmd0aD1uX2JpbnMrMSkpCnh4X2JpbiA8LSByZXAoTkEsIG5fYmlucykKeXlfYmluIDwtIHJlcChOQSwgbl9iaW5zKQpmb3IgKGkgaW4gMTpuX2JpbnMpewogIGtlZXAgPC0geHggPiBjdXRwb2ludHNbaV0gJiB4eCA8PSBjdXRwb2ludHNbaSsxXQogIHh4X2JpbltpXSA8LSBtZWFuKHh4W2tlZXBdKQogIHl5X2JpbltpXSA8LSBtZWFuKHl5W2tlZXBdKQp9CnBsb3QoeT15eV9iaW4sIHg9eHhfYmluLCB0eXBlPSJuIiwgCiAgICAgeGxhYj0iQXNzaWdubWVudCB2YXJpYWJsZSIsCiAgICAgeWxhYj0iQWRqdXN0ZWQgb3V0Y29tZSIsIGJ0eT0ibCIsIHlheHQ9Im4iLCB4bGltPWMoLWhhbGZ3aWR0aCwgaGFsZndpZHRoKSwgeGF4cz0iaSIpCmF4aXMoMiwgc2VxKC0xMDAsMTAwLDEpKQpwb2ludHMoeHhfYmluW3h4X2Jpbj4wXSwgeXlfYmluW3h4X2Jpbj4wXSxwY2g9MjAsCiAgICAgICBjZXg9MiwgY29sPSJncmF5NjAiKQpwb2ludHMoeHhfYmluW3h4X2JpbjwwXSwgeXlfYmluW3h4X2JpbjwwXSxwY2g9MjAsCiAgICAgICBjZXg9MikKd2hpdGVsaW5lKGNvZWYoZml0XzNiKVsxXSArIGNvZWYoZml0XzNiKVs0XSAqbWVhbihjaGlsZSRyZWFkODgpICsgY29lZihmaXRfM2IpWzVdICptZWFuKGNoaWxlJG1hdGg4OCksIGNvZWYoZml0XzNiKVszXSwgZnJvbT0wKQp3aGl0ZWxpbmUoY29lZihmaXRfM2IpWzFdICsgY29lZihmaXRfM2IpWzRdICptZWFuKGNoaWxlJHJlYWQ4OCkgKyBjb2VmKGZpdF8zYilbNV0gKm1lYW4oY2hpbGUkbWF0aDg4KSArIGNvZWYoZml0XzNiKVsyXSwgY29lZihmaXRfM2IpWzNdLCB0bz0wKQphYmxpbmUodj0wLCBsd2Q9NCwgY29sPSJncmF5IikKbXRleHQoIkJpbm5lZCBhdmVyYWdlcyB3aXRoIHNhbWUgcmVncmVzc2lvbiBsaW5lcyIsIDMsIDEpCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQpkZXYub2ZmKCkKYGBgCgojIyBBZGRpdGlvbmFsIG1vZGVscwojIyMjIE5vbmNvbXBsaWFuY2UgcmF0ZXMKCmBgYHtyIH0KbWVhbihjaGlsZSRwOTBbY2hpbGUkZWxpZ2libGU9PTAgJiBhYnMoY2hpbGUkcnVsZTIpIDwgNV0pICAjIC4wNTMKbWVhbihjaGlsZSRwOTBbY2hpbGUkZWxpZ2libGU9PTEgJiBhYnMoY2hpbGUkcnVsZTIpIDwgNV0pICMgLjYwNgpgYGAKCiMjIyMgSVYgbW9kZWwgb24gcmVzdHJpY3RlZCBkYXRhc2V0CgpgYGB7ciB9CiMgd2l0aCBicm1zIGFzIGluIElWIHNlY3Rpb24gZWFybGllciBpbiB0aGUgY2hhcHRlcgpzZXQuc2VlZCgxMjM0KQpjaGlsZSRkaWZmX3JlYWQ5MiA8LSBjaGlsZSRyZWFkOTIgLSBjaGlsZSRyZWFkODgKcmRfZjEgPC0gYmYocDkwIH4gZWxpZ2libGUpCnJkX2YyIDwtIGJmKGRpZmZfcmVhZDkyIH4gcDkwKQpgYGAKYGBge3IgcmVzdWx0cz0naGlkZSd9CklWX2JybV9yZCA8LSBicm0oZm9ybXVsYT1yZF9mMSArIHJkX2YyLCBkYXRhID0gY2hpbGVbYWJzKGNoaWxlJHJ1bGUyKSA8IDUsXSkKYGBgCmBgYHtyIH0KcHJpbnQoSVZfYnJtX3JkLCBkaWdpdHM9MykKIyB0LmUuIGlzIDUuMDggd2l0aCBzLmUuIG9mIC45MApgYGAKCiMjIyMgQ29tcGFyZSB0byBJViByZWdyZXNzaW9uIGJ5IHR3by1zdGFnZSBsZWFzdCBzcXVhcmVzCgpgYGB7ciB9CnN1bW1hcnkoaXZyZWcoZm9ybXVsYSA9IGRpZmZfcmVhZDkyIH4gcDkwIHwgZWxpZ2libGUsIGRhdGE9Y2hpbGUsIHN1YnNldCA9IGFicyhydWxlMikgPCA1KSkKIyBlc3QgaXMgNS4wOSB3aXRoIHMuZS4gb2YgLjg4CmBgYAoK