The Children’s television workshop experiment.

## Data input

W <- rep(c(0, 1), 8)
G <- as.vector(rbind(1:8, 1:8))
X <- c(12.9, 12, 15.1, 12.3, 16.8, 17.2, 15.8, 18.9, 13.9, 15.3, 14.5, 16.6, 17, 16, 15.8, 20.1)
Y <- c(54.6, 60.6, 56.5, 55.5, 75.2, 84.8, 75.6, 101.9, 55.3, 70.6, 59.3, 78.4, 87, 84.2, 73.7, 108.6)
data <- data.frame(G, W, X, Y)
data

library(ggpubr)
Loading required package: ggplot2
temp <- data.frame("Control" = data$Y[data$W == 0], "Treatment" = data$Y[data$W == 1])
temp
ggpaired(data = temp, cond1 = "Control", cond2 = "Treatment", fill = "condition", palette = "jco", line.color = "gray", 
         ylab = "Post-test Score")


temp <- data.frame("Control" = data$X[data$W == 0], "Treatment" = data$X[data$W == 1])
ggpaired(data = temp, cond1 = "Control", cond2 = "Treatment", fill = "condition", palette = "jco", line.color = "gray", 
         ylab = "Pre-test Score")


testStat <- function(Y0, Y1, w, pair_idx, method = "mean diff") {
 n <- length(Y0)
 if (method == "mean diff")
   val <- abs(mean(Y0[w == 1]) - mean(Y0[w == 0])) ## replace abs(mean(Y1[w == 1]) - mean(Y0[w == 0])-tau)
 else if (method == "rank diff"){
  r <- rank(Y0)
  val <- abs(mean(r[w == 1]) - mean(r[w == 0]))
 } else if (method == "within-pair rank diff") {
   pairs <- unique(pair_idx)
   v <- sapply(pairs, function(idx) {
     Yt <- Y0[pair_idx == idx & w == 1]
     Yc <- Y0[pair_idx == idx & w == 0]
     return((Yt > Yc) - (Yt < Yc))
   })
   val <- abs(2 * sum(v)/length(Y0))
 }
  return(val)
}

find.p.value <- function(tau, yobs, w, pair_idx, seed = 0, method = "mean diff", 
                         B = 10000, error = 1e-10){
  ## error is to fix our problems with floating point errors
  n <- length(yobs)
  y0 <- rep(NA, n)
  y1 <- rep(NA, n)
  
  y1[w == 1] <- yobs[w == 1]
  y0[w == 0] <- yobs[w == 0]
  
  y1[w == 0] <- yobs[w == 0] + tau
  y0[w == 1] <- yobs[w == 1] - tau
  
  t.obs <-  testStat(y0, y1, w, pair_idx, method = method)
  
  
  set.seed(seed)
  pairs <- unique(pair_idx)
  t.stat <- sapply(1:B, function(i){
    w.iter <- rep(NA, n)
    for (i in 1:length(pairs)) {
      temp.idx <- which(pair_idx == pairs[i])
      w.iter[sample(temp.idx)] <- c(0, 1)
    }
    return(testStat(y0, y1, w.iter, pair_idx, method = method))
  })

  # subtracting error fixes floating point errors
  p.val <- sum(t.stat >= (t.obs - error))/B
  res <- c(t.obs, p.val)
  names(res) <- paste(method, c("test stats", "p-value"))
  return(res)
}


find.p.value(0, data$Y, data$W, data$G)
mean diff test stats    mean diff p-value 
             13.4250               0.0332 
find.p.value(0, data$Y, data$W, data$G, method = "rank diff")
rank diff test stats    rank diff p-value 
              3.7500               0.0332 
find.p.value(0, data$Y, data$W, data$G, method = "within-pair rank diff")
within-pair rank diff test stats    within-pair rank diff p-value 
                             0.5                              0.3 
find.p.value.wrong <- function(tau, yobs, w, pair_idx, seed = 0, method = "mean diff", 
                         B = 10000, error = 1e-10){
  ## error is to fix our problems with floating point errors
  n <- length(yobs)
  y0 <- rep(NA, n)
  y1 <- rep(NA, n)
  
  y1[w == 1] <- yobs[w == 1]
  y0[w == 0] <- yobs[w == 0]
  
  y1[w == 0] <- yobs[w == 0] + tau
  y0[w == 1] <- yobs[w == 1] - tau
  
  t.obs <-  testStat(y0, y1, w, pair_idx, method)
  
  
  set.seed(seed)
  pairs <- unique(pair_idx)
  t.stat <- sapply(1:B, function(i){
    w.iter <- sample(w)
    return(testStat(y0, y1, w.iter, pair_idx, method = method))
  })

  # subtracting error fixes floating point errors
  p.val <- sum(t.stat >= (t.obs - error))/B
  res <- c(t.obs, p.val)
  names(res) <- paste(method, c("test stats", "p-value"))
  return(res)
}


find.p.value.wrong(0, data$Y, data$W, data$G)
mean diff test stats    mean diff p-value 
             13.4250               0.1058 
find.p.value.wrong(0, data$Y, data$W, data$G, method = "rank diff")
rank diff test stats    rank diff p-value 
              3.7500               0.1282 

We get larger p-values from complete randomization mechanism.

## CI from Fisher
pvals <- sapply(0:30, function(tau) find.p.value(tau, data$Y, data$W, data$G)[2])

The 95% CI from Fisher’s exact p-values is [1.9, 25.6]

## CI from Neyman
pairs <- unique(data$G)
tau_hat_vec <- sapply(pairs, function(idx) {
  Yt <- data$Y[data$G == idx & data$W == 1]
  Yc <- data$Y[data$G == idx & data$W == 0]
  return(Yt - Yc)
})

tau_hat <- mean(tau_hat_vec)
sd_tau_hat <- sd(tau_hat_vec)/sqrt(length(tau_hat_vec))
  
## 95% CI
print(c(tau_hat- 1.96 * sd_tau_hat, tau_hat + 1.96 * sd_tau_hat))
[1]  4.337779 22.512221
## CI from Neyman

tau_hat <- mean(data$Y[data$W == 1]) - mean(data$Y[data$W == 0])

st <- sd(data$Y[data$W == 1])
sc <- sd(data$Y[data$W == 0])
varhat_tauhat <- st^2/sum(data$W == 1) + sc^2/sum(data$W == 0)

print(c(tau_hat - 1.96 * sqrt(varhat_tauhat), tau_hat + 1.96 * sqrt(varhat_tauhat)))
[1] -1.957377 28.807377
## need to create pair-level data
pairs <- unique(data$G)
X_bar <- mean(data$X)
group_data <- sapply(pairs, function(idx) {
  Yt <- data$Y[data$G == idx & data$W == 1]
  Yc <- data$Y[data$G == idx & data$W == 0]
  Xt <- data$X[data$G == idx & data$W == 1]
  Xc <- data$X[data$G == idx & data$W == 0]
  Y_diff <- Yt - Yc
  X_diff <- Xt - Xc
  X_mean <- (Xt + Xc)/2 - X_bar
  return(c(X_mean, X_diff, Y_diff))
})
group_data
        [,1]    [,2]   [,3]    [,4]    [,5]    [,6]    [,7]    [,8]
[1,] -3.1875 -1.9375 1.3625  1.7125 -1.0375 -0.0875  0.8625  2.3125
[2,] -0.9000 -2.8000 0.4000  3.1000  1.4000  2.1000 -1.0000  4.3000
[3,]  6.0000 -1.0000 9.6000 26.3000 15.3000 19.1000 -2.8000 34.9000
group_data <- data.frame(t(group_data))
colnames(group_data) <- c("X_mean", "X_diff", "Y_diff")
group_data

lm_result1 <- lm(Y_diff ~ X_diff, data =group_data)
summary(lm_result1) ## substantially smaller sd of the estimate

Call:
lm(formula = Y_diff ~ X_diff, data = group_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4238 -1.2954 -0.2577  2.0825  5.0432 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.9943     1.4096   6.381 0.000697 ***
X_diff        5.3705     0.5992   8.964 0.000108 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.734 on 6 degrees of freedom
Multiple R-squared:  0.9305,    Adjusted R-squared:  0.9189 
F-statistic: 80.34 on 1 and 6 DF,  p-value: 0.0001077
lm_result2 <- lm(Y_diff ~ X_diff + X_mean, data =group_data)
summary(lm_result2) 

Call:
lm(formula = Y_diff ~ X_diff + X_mean, data = group_data)

Residuals:
      1       2       3       4       5       6       7       8 
-0.4555  5.1236  0.1038  1.1179 -2.6106 -1.9925 -4.4907  3.2041 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.5220     1.4745   5.780 0.002181 ** 
X_diff        5.9431     0.8138   7.303 0.000754 ***
X_mean       -1.0298     0.9969  -1.033 0.348978    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.713 on 5 degrees of freedom
Multiple R-squared:  0.9427,    Adjusted R-squared:  0.9198 
F-statistic: 41.15 on 2 and 5 DF,  p-value: 0.0007848
LS0tCnRpdGxlOiAiUiBleGFtcGxlMzogdGhlIENoaWxkcmVu4oCZcyB0ZWxldmlzaW9uIHdvcmtzaG9wIGV4cGVyaW1lbnQuIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGUgQ2hpbGRyZW7igJlzIHRlbGV2aXNpb24gd29ya3Nob3AgZXhwZXJpbWVudC4gIAoKYGBge3J9CiMjIERhdGEgaW5wdXQKClcgPC0gcmVwKGMoMCwgMSksIDgpCkcgPC0gYXMudmVjdG9yKHJiaW5kKDE6OCwgMTo4KSkKWCA8LSBjKDEyLjksIDEyLCAxNS4xLCAxMi4zLCAxNi44LCAxNy4yLCAxNS44LCAxOC45LCAxMy45LCAxNS4zLCAxNC41LCAxNi42LCAxNywgMTYsIDE1LjgsIDIwLjEpClkgPC0gYyg1NC42LCA2MC42LCA1Ni41LCA1NS41LCA3NS4yLCA4NC44LCA3NS42LCAxMDEuOSwgNTUuMywgNzAuNiwgNTkuMywgNzguNCwgODcsIDg0LjIsIDczLjcsIDEwOC42KQpkYXRhIDwtIGRhdGEuZnJhbWUoRywgVywgWCwgWSkKZGF0YQoKbGlicmFyeShnZ3B1YnIpCgp0ZW1wIDwtIGRhdGEuZnJhbWUoIkNvbnRyb2wiID0gZGF0YSRZW2RhdGEkVyA9PSAwXSwgIlRyZWF0bWVudCIgPSBkYXRhJFlbZGF0YSRXID09IDFdKQp0ZW1wCmdncGFpcmVkKGRhdGEgPSB0ZW1wLCBjb25kMSA9ICJDb250cm9sIiwgY29uZDIgPSAiVHJlYXRtZW50IiwgZmlsbCA9ICJjb25kaXRpb24iLCBwYWxldHRlID0gImpjbyIsIGxpbmUuY29sb3IgPSAiZ3JheSIsIAogICAgICAgICB5bGFiID0gIlBvc3QtdGVzdCBTY29yZSIpCgp0ZW1wIDwtIGRhdGEuZnJhbWUoIkNvbnRyb2wiID0gZGF0YSRYW2RhdGEkVyA9PSAwXSwgIlRyZWF0bWVudCIgPSBkYXRhJFhbZGF0YSRXID09IDFdKQpnZ3BhaXJlZChkYXRhID0gdGVtcCwgY29uZDEgPSAiQ29udHJvbCIsIGNvbmQyID0gIlRyZWF0bWVudCIsIGZpbGwgPSAiY29uZGl0aW9uIiwgcGFsZXR0ZSA9ICJqY28iLCBsaW5lLmNvbG9yID0gImdyYXkiLCAKICAgICAgICAgeWxhYiA9ICJQcmUtdGVzdCBTY29yZSIpCmBgYAoKCi0gQ29tcHV0ZSBGaXNoZXIncyBleGFjdCBwLXZhbHVlCgpgYGB7cn0KCnRlc3RTdGF0IDwtIGZ1bmN0aW9uKFkwLCBZMSwgdywgcGFpcl9pZHgsIG1ldGhvZCA9ICJtZWFuIGRpZmYiKSB7CiBuIDwtIGxlbmd0aChZMCkKIGlmIChtZXRob2QgPT0gIm1lYW4gZGlmZiIpCiAgIHZhbCA8LSBhYnMobWVhbihZMFt3ID09IDFdKSAtIG1lYW4oWTBbdyA9PSAwXSkpICMjIHJlcGxhY2UgYWJzKG1lYW4oWTFbdyA9PSAxXSkgLSBtZWFuKFkwW3cgPT0gMF0pLXRhdSkKIGVsc2UgaWYgKG1ldGhvZCA9PSAicmFuayBkaWZmIil7CiAgciA8LSByYW5rKFkwKQogIHZhbCA8LSBhYnMobWVhbihyW3cgPT0gMV0pIC0gbWVhbihyW3cgPT0gMF0pKQogfSBlbHNlIGlmIChtZXRob2QgPT0gIndpdGhpbi1wYWlyIHJhbmsgZGlmZiIpIHsKICAgcGFpcnMgPC0gdW5pcXVlKHBhaXJfaWR4KQogICB2IDwtIHNhcHBseShwYWlycywgZnVuY3Rpb24oaWR4KSB7CiAgICAgWXQgPC0gWTBbcGFpcl9pZHggPT0gaWR4ICYgdyA9PSAxXQogICAgIFljIDwtIFkwW3BhaXJfaWR4ID09IGlkeCAmIHcgPT0gMF0KICAgICByZXR1cm4oKFl0ID4gWWMpIC0gKFl0IDwgWWMpKQogICB9KQogICB2YWwgPC0gYWJzKDIgKiBzdW0odikvbGVuZ3RoKFkwKSkKIH0KICByZXR1cm4odmFsKQp9CgpmaW5kLnAudmFsdWUgPC0gZnVuY3Rpb24odGF1LCB5b2JzLCB3LCBwYWlyX2lkeCwgc2VlZCA9IDAsIG1ldGhvZCA9ICJtZWFuIGRpZmYiLCAKICAgICAgICAgICAgICAgICAgICAgICAgIEIgPSAxMDAwMCwgZXJyb3IgPSAxZS0xMCl7CiAgIyMgZXJyb3IgaXMgdG8gZml4IG91ciBwcm9ibGVtcyB3aXRoIGZsb2F0aW5nIHBvaW50IGVycm9ycwogIG4gPC0gbGVuZ3RoKHlvYnMpCiAgeTAgPC0gcmVwKE5BLCBuKQogIHkxIDwtIHJlcChOQSwgbikKICAKICB5MVt3ID09IDFdIDwtIHlvYnNbdyA9PSAxXQogIHkwW3cgPT0gMF0gPC0geW9ic1t3ID09IDBdCiAgCiAgeTFbdyA9PSAwXSA8LSB5b2JzW3cgPT0gMF0gKyB0YXUKICB5MFt3ID09IDFdIDwtIHlvYnNbdyA9PSAxXSAtIHRhdQogIAogIHQub2JzIDwtICB0ZXN0U3RhdCh5MCwgeTEsIHcsIHBhaXJfaWR4LCBtZXRob2QgPSBtZXRob2QpCiAgCiAgCiAgc2V0LnNlZWQoc2VlZCkKICBwYWlycyA8LSB1bmlxdWUocGFpcl9pZHgpCiAgdC5zdGF0IDwtIHNhcHBseSgxOkIsIGZ1bmN0aW9uKGkpewogICAgdy5pdGVyIDwtIHJlcChOQSwgbikKICAgIGZvciAoaSBpbiAxOmxlbmd0aChwYWlycykpIHsKICAgICAgdGVtcC5pZHggPC0gd2hpY2gocGFpcl9pZHggPT0gcGFpcnNbaV0pCiAgICAgIHcuaXRlcltzYW1wbGUodGVtcC5pZHgpXSA8LSBjKDAsIDEpCiAgICB9CiAgICByZXR1cm4odGVzdFN0YXQoeTAsIHkxLCB3Lml0ZXIsIHBhaXJfaWR4LCBtZXRob2QgPSBtZXRob2QpKQogIH0pCgogICMgc3VidHJhY3RpbmcgZXJyb3IgZml4ZXMgZmxvYXRpbmcgcG9pbnQgZXJyb3JzCiAgcC52YWwgPC0gc3VtKHQuc3RhdCA+PSAodC5vYnMgLSBlcnJvcikpL0IKICByZXMgPC0gYyh0Lm9icywgcC52YWwpCiAgbmFtZXMocmVzKSA8LSBwYXN0ZShtZXRob2QsIGMoInRlc3Qgc3RhdHMiLCAicC12YWx1ZSIpKQogIHJldHVybihyZXMpCn0KCgpmaW5kLnAudmFsdWUoMCwgZGF0YSRZLCBkYXRhJFcsIGRhdGEkRykKZmluZC5wLnZhbHVlKDAsIGRhdGEkWSwgZGF0YSRXLCBkYXRhJEcsIG1ldGhvZCA9ICJyYW5rIGRpZmYiKQpmaW5kLnAudmFsdWUoMCwgZGF0YSRZLCBkYXRhJFcsIGRhdGEkRywgbWV0aG9kID0gIndpdGhpbi1wYWlyIHJhbmsgZGlmZiIpCgpgYGAKCgotIFdoYXQgaWYgd2UgdXNlIGlnbm9yZSB0aGUgcGFpcndpc2UgcmFuZG9taXphdGlvbiBtZWNoYW5pc20gYW5kIHRyZWF0IGl0IGFzIGNvbXBsZXRlIHJhbmRvbWl6YXRpb24/CgpgYGB7cn0KZmluZC5wLnZhbHVlLndyb25nIDwtIGZ1bmN0aW9uKHRhdSwgeW9icywgdywgcGFpcl9pZHgsIHNlZWQgPSAwLCBtZXRob2QgPSAibWVhbiBkaWZmIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICBCID0gMTAwMDAsIGVycm9yID0gMWUtMTApewogICMjIGVycm9yIGlzIHRvIGZpeCBvdXIgcHJvYmxlbXMgd2l0aCBmbG9hdGluZyBwb2ludCBlcnJvcnMKICBuIDwtIGxlbmd0aCh5b2JzKQogIHkwIDwtIHJlcChOQSwgbikKICB5MSA8LSByZXAoTkEsIG4pCiAgCiAgeTFbdyA9PSAxXSA8LSB5b2JzW3cgPT0gMV0KICB5MFt3ID09IDBdIDwtIHlvYnNbdyA9PSAwXQogIAogIHkxW3cgPT0gMF0gPC0geW9ic1t3ID09IDBdICsgdGF1CiAgeTBbdyA9PSAxXSA8LSB5b2JzW3cgPT0gMV0gLSB0YXUKICAKICB0Lm9icyA8LSAgdGVzdFN0YXQoeTAsIHkxLCB3LCBwYWlyX2lkeCwgbWV0aG9kKQogIAogIAogIHNldC5zZWVkKHNlZWQpCiAgcGFpcnMgPC0gdW5pcXVlKHBhaXJfaWR4KQogIHQuc3RhdCA8LSBzYXBwbHkoMTpCLCBmdW5jdGlvbihpKXsKICAgIHcuaXRlciA8LSBzYW1wbGUodykKICAgIHJldHVybih0ZXN0U3RhdCh5MCwgeTEsIHcuaXRlciwgcGFpcl9pZHgsIG1ldGhvZCA9IG1ldGhvZCkpCiAgfSkKCiAgIyBzdWJ0cmFjdGluZyBlcnJvciBmaXhlcyBmbG9hdGluZyBwb2ludCBlcnJvcnMKICBwLnZhbCA8LSBzdW0odC5zdGF0ID49ICh0Lm9icyAtIGVycm9yKSkvQgogIHJlcyA8LSBjKHQub2JzLCBwLnZhbCkKICBuYW1lcyhyZXMpIDwtIHBhc3RlKG1ldGhvZCwgYygidGVzdCBzdGF0cyIsICJwLXZhbHVlIikpCiAgcmV0dXJuKHJlcykKfQoKCmZpbmQucC52YWx1ZS53cm9uZygwLCBkYXRhJFksIGRhdGEkVywgZGF0YSRHKQpmaW5kLnAudmFsdWUud3JvbmcoMCwgZGF0YSRZLCBkYXRhJFcsIGRhdGEkRywgbWV0aG9kID0gInJhbmsgZGlmZiIpCmBgYAogV2UgZ2V0IGxhcmdlciBwLXZhbHVlcyBmcm9tIGNvbXBsZXRlIHJhbmRvbWl6YXRpb24gbWVjaGFuaXNtLgogCiAtIENvbXBhcmUgOTVcJSBDSSBvZiBGaXNoZXIgVi5TLiBOZXltYW4KIApgYGB7cn0KIyMgQ0kgZnJvbSBGaXNoZXIKcHZhbHMgPC0gc2FwcGx5KDA6MzAsIGZ1bmN0aW9uKHRhdSkgZmluZC5wLnZhbHVlKHRhdSwgZGF0YSRZLCBkYXRhJFcsIGRhdGEkRylbMl0pCm5hbWVzKHB2YWxzKSA8LSAwOjMwCnB2YWxzCgpMQiA8LSBzYXBwbHkoc2VxKDEsIDIsIDAuMSksIGZ1bmN0aW9uKHRhdSkgZmluZC5wLnZhbHVlKHRhdSwgZGF0YSRZLCBkYXRhJFcsIGRhdGEkRylbMl0pCm5hbWVzKExCKSA8LSBzZXEoMSwgMiwgMC4xKQpMQgoKVUIgPC0gc2FwcGx5KHNlcSgyNSwgMjYsIDAuMSksIGZ1bmN0aW9uKHRhdSkgZmluZC5wLnZhbHVlKHRhdSwgZGF0YSRZLCBkYXRhJFcsIGRhdGEkRylbMl0pCm5hbWVzKFVCKSA8LSBzZXEoMjUsIDI2LCAwLjEpClVCCmBgYAogVGhlIDk1XCUgQ0kgZnJvbSBGaXNoZXIncyBleGFjdCBwLXZhbHVlcyBpcyBbMS45LCAyNS42XQogCgpgYGB7cn0KIyMgQ0kgZnJvbSBOZXltYW4KcGFpcnMgPC0gdW5pcXVlKGRhdGEkRykKdGF1X2hhdF92ZWMgPC0gc2FwcGx5KHBhaXJzLCBmdW5jdGlvbihpZHgpIHsKICBZdCA8LSBkYXRhJFlbZGF0YSRHID09IGlkeCAmIGRhdGEkVyA9PSAxXQogIFljIDwtIGRhdGEkWVtkYXRhJEcgPT0gaWR4ICYgZGF0YSRXID09IDBdCiAgcmV0dXJuKFl0IC0gWWMpCn0pCgp0YXVfaGF0IDwtIG1lYW4odGF1X2hhdF92ZWMpCnNkX3RhdV9oYXQgPC0gc2QodGF1X2hhdF92ZWMpL3NxcnQobGVuZ3RoKHRhdV9oYXRfdmVjKSkKICAKIyMgOTUlIENJCnByaW50KGModGF1X2hhdC0gMS45NiAqIHNkX3RhdV9oYXQsIHRhdV9oYXQgKyAxLjk2ICogc2RfdGF1X2hhdCkpCmBgYAoKLSBXaGF0IGRvZXMgTmV5bWFuJ3MgQ0kgbG9vayBsaWtlIGlmIHdlIHRyZWF0IHRoZSBkYXRhIGFzIGZyb20gY29tcGxldGVseSByYW5kb21pemVkIGV4cGVyaW1lbnQKCmBgYHtyfQojIyBDSSBmcm9tIE5leW1hbgoKdGF1X2hhdCA8LSBtZWFuKGRhdGEkWVtkYXRhJFcgPT0gMV0pIC0gbWVhbihkYXRhJFlbZGF0YSRXID09IDBdKQoKc3QgPC0gc2QoZGF0YSRZW2RhdGEkVyA9PSAxXSkKc2MgPC0gc2QoZGF0YSRZW2RhdGEkVyA9PSAwXSkKdmFyaGF0X3RhdWhhdCA8LSBzdF4yL3N1bShkYXRhJFcgPT0gMSkgKyBzY14yL3N1bShkYXRhJFcgPT0gMCkKCnByaW50KGModGF1X2hhdCAtIDEuOTYgKiBzcXJ0KHZhcmhhdF90YXVoYXQpLCB0YXVfaGF0ICsgMS45NiAqIHNxcnQodmFyaGF0X3RhdWhhdCkpKQpgYGAKCi0gTGluZWFyIHJlZ3Jlc3Npb24gcmVzdWx0cwoKYGBge3J9CiMjIG5lZWQgdG8gY3JlYXRlIHBhaXItbGV2ZWwgZGF0YQpwYWlycyA8LSB1bmlxdWUoZGF0YSRHKQpYX2JhciA8LSBtZWFuKGRhdGEkWCkKZ3JvdXBfZGF0YSA8LSBzYXBwbHkocGFpcnMsIGZ1bmN0aW9uKGlkeCkgewogIFl0IDwtIGRhdGEkWVtkYXRhJEcgPT0gaWR4ICYgZGF0YSRXID09IDFdCiAgWWMgPC0gZGF0YSRZW2RhdGEkRyA9PSBpZHggJiBkYXRhJFcgPT0gMF0KICBYdCA8LSBkYXRhJFhbZGF0YSRHID09IGlkeCAmIGRhdGEkVyA9PSAxXQogIFhjIDwtIGRhdGEkWFtkYXRhJEcgPT0gaWR4ICYgZGF0YSRXID09IDBdCiAgWV9kaWZmIDwtIFl0IC0gWWMKICBYX2RpZmYgPC0gWHQgLSBYYwogIFhfbWVhbiA8LSAoWHQgKyBYYykvMiAtIFhfYmFyCiAgcmV0dXJuKGMoWF9tZWFuLCBYX2RpZmYsIFlfZGlmZikpCn0pCmdyb3VwX2RhdGEKCmdyb3VwX2RhdGEgPC0gZGF0YS5mcmFtZSh0KGdyb3VwX2RhdGEpKQpjb2xuYW1lcyhncm91cF9kYXRhKSA8LSBjKCJYX21lYW4iLCAiWF9kaWZmIiwgIllfZGlmZiIpCmdyb3VwX2RhdGEKCmxtX3Jlc3VsdDEgPC0gbG0oWV9kaWZmIH4gWF9kaWZmLCBkYXRhID1ncm91cF9kYXRhKQpzdW1tYXJ5KGxtX3Jlc3VsdDEpICMjIHN1YnN0YW50aWFsbHkgc21hbGxlciBzZCBvZiB0aGUgZXN0aW1hdGUKCmxtX3Jlc3VsdDIgPC0gbG0oWV9kaWZmIH4gWF9kaWZmICsgWF9tZWFuLCBkYXRhID1ncm91cF9kYXRhKQpzdW1tYXJ5KGxtX3Jlc3VsdDIpIApgYGA=