Regression-based imputation for the Social Indicators Survey. See Chapter 17 in Regression and Other Stories.
dplyr + ggplot version.


Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")
library("dplyr")
library("ggplot2")
theme_set(bayesplot::theme_default(base_family = "sans"))

Load data

SIS <- read.csv(root("Imputation/data","SIS.csv"))
head(SIS)
  earnings retirement interest assistance other male over65 white immig educ_r
1     84.0        0.7     0.20          0   0.0    1      0     1     1      4
2       NA        0.0     0.00          0   0.0    0      0     0     1      4
3     27.5        0.0     0.16          0   0.0    1      0     0     1      2
4     85.0       12.0     5.00          0    NA    1      1     1     0      4
5    135.0        0.0     0.10          0   0.1    1      0     0     0      4
6      0.0         NA       NA         NA    NA    0      0     0     0      4
  workmos workhrs_top any_ssi any_welfare any_charity
1      12          40       0           0           0
2      12          40       0           0           0
3      10          40       0           0           0
4      11           8       0           0           0
5      12          40       0           0           0
6       0          40       0           0           0
summary(SIS)
    earnings         retirement        interest         assistance     
 Min.   :   0.00   Min.   :  0.00   Min.   :  0.000   Min.   : 0.0000  
 1st Qu.:   4.80   1st Qu.:  0.00   1st Qu.:  0.000   1st Qu.: 0.0000  
 Median :  28.00   Median :  0.00   Median :  0.000   Median : 0.0000  
 Mean   :  52.19   Mean   :  2.19   Mean   :  1.722   Mean   : 0.6821  
 3rd Qu.:  65.00   3rd Qu.:  0.00   3rd Qu.:  0.000   3rd Qu.: 0.0000  
 Max.   :3250.00   Max.   :300.00   Max.   :650.000   Max.   :36.4000  
 NA's   :241       NA's   :123      NA's   :277       NA's   :128      
     other             male            over65            white       
 Min.   : 0.000   Min.   :0.0000   Min.   :0.00000   Min.   :0.0000  
 1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000  
 Median : 0.000   Median :0.0000   Median :0.00000   Median :0.0000  
 Mean   : 0.577   Mean   :0.3618   Mean   :0.07728   Mean   :0.3185  
 3rd Qu.: 0.025   3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:1.0000  
 Max.   :96.000   Max.   :1.0000   Max.   :1.00000   Max.   :1.0000  
 NA's   :133                                                         
     immig            educ_r         workmos        workhrs_top   
 Min.   :0.0000   Min.   :1.000   Min.   : 0.000   Min.   : 0.00  
 1st Qu.:0.0000   1st Qu.:2.000   1st Qu.: 7.000   1st Qu.:26.00  
 Median :0.0000   Median :3.000   Median :12.000   Median :40.00  
 Mean   :0.3924   Mean   :2.747   Mean   : 8.985   Mean   :30.19  
 3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:12.000   3rd Qu.:40.00  
 Max.   :1.0000   Max.   :4.000   Max.   :12.000   Max.   :40.00  
                                                                  
    any_ssi         any_welfare       any_charity      
 Min.   :0.00000   Min.   :0.00000   Min.   :0.000000  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.000000  
 Median :0.00000   Median :0.00000   Median :0.000000  
 Mean   :0.03664   Mean   :0.03931   Mean   :0.009993  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.000000  
 Max.   :1.00000   Max.   :1.00000   Max.   :1.000000  
                                                       

Imputation helper functions
Create a completed data vector using imputations

impute <- function(a, a_impute) {
  ifelse(is.na(a), a_impute, a)
}

Top code function

topcode <- function(a, top) {
  ifelse(a>top, top, a)
}

se missing codes to NA

na_fix <- function(a) {
  ifelse(a<0 | a==999999, NA, a)
}
is_any <- function(a) {
  any_a <- ifelse(a>0, 1, 0)
  any_a[is.na(a)] <- 0
  any_a
}

1. Deterministic imputation

Impute 0 earnings using the logical rule (if worked 0 months and 0 hrs/wk)

SIS <- mutate(SIS,
  zero_earnings = (workhrs_top==0 & workmos==0),
  earnings_top = if_else(zero_earnings, 0, topcode(earnings, 100)))

Create a little dataset with all our redefined variables

SIS <- select(SIS,
              earnings, earnings_top, interest, male, over65,
              white, immig, educ_r, workmos, workhrs_top,
              any_ssi, any_welfare, any_charity)

Impute subset of earnings that are nonzero: linear scale

fit_imp_1 <- stan_glm(
  earnings ~ male + over65 + white + immig + educ_r +
              workmos + workhrs_top + any_ssi +
              any_welfare + any_charity,
  data = SIS,
  subset = earnings > 0,
  refresh = 0
)
print(fit_imp_1)
stan_glm
 family:       gaussian [identity]
 formula:      earnings ~ male + over65 + white + immig + educ_r + workmos + 
       workhrs_top + any_ssi + any_welfare + any_charity
 observations: 988
 predictors:   11
------
            Median MAD_SD
(Intercept) -87.0   30.0 
male         -0.5    9.1 
over65      -40.3   39.9 
white        26.2    9.9 
immig       -12.3    9.3 
educ_r       22.6    4.3 
workmos       5.4    2.1 
workhrs_top   0.8    0.6 
any_ssi     -34.3   38.6 
any_welfare -11.6   24.5 
any_charity -29.8   40.7 

Auxiliary parameter(s):
      Median MAD_SD
sigma 132.2    3.1 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
SIS_predictors <- select(SIS, -starts_with("earnings"))
pred_1 <- colMeans(posterior_linpred(fit_imp_1, newdata = SIS_predictors))
SIS <- mutate(SIS,
  earnings_imp_1 = impute(earnings, pred_1))

Impute subset of earnings that are nonzero: square root scale and topcoding

fit_imp_2 <- stan_glm(
  sqrt(earnings_top) ~ male + over65 + white + immig +
                       educ_r + workmos + workhrs_top + any_ssi +
                       any_welfare + any_charity,
  data = SIS,
  subset = earnings > 0,
  refresh = 0
)
print(fit_imp_2)
stan_glm
 family:       gaussian [identity]
 formula:      sqrt(earnings_top) ~ male + over65 + white + immig + educ_r + 
       workmos + workhrs_top + any_ssi + any_welfare + any_charity
 observations: 988
 predictors:   11
------
            Median MAD_SD
(Intercept) -1.7    0.4  
male         0.3    0.1  
over65      -1.4    0.6  
white        1.0    0.1  
immig       -0.6    0.1  
educ_r       0.8    0.1  
workmos      0.3    0.0  
workhrs_top  0.1    0.0  
any_ssi     -0.9    0.5  
any_welfare -1.3    0.4  
any_charity -1.2    0.6  

Auxiliary parameter(s):
      Median MAD_SD
sigma 2.0    0.0   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

point predictions

pred_2_sqrt <- colMeans(posterior_linpred(fit_imp_2, newdata = SIS_predictors))
pred_2 <- topcode(pred_2_sqrt^2, 100)
SIS <- mutate(SIS,
  pred_2 = pred_2,
  earnings_imp_2 = impute(earnings_top, pred_2))

2. Random imputation

Linear scale (use fitted model fit_imp_1)

pred_3 <- posterior_predict(fit_imp_1, newdata = SIS_predictors, draws = 1)
SIS <- mutate(SIS,
  earnings_imp_3 = impute(earnings, pred_3))

Square root scale and topcoding (use fitted model fit_imp_2)

pred_4_sqrt <- posterior_predict(fit_imp_2, newdata = SIS_predictors, draws = 1)
pred_4 <- topcode(pred_4_sqrt^2, 100)
SIS <- mutate(SIS,
  earnings_imp_4 = impute(earnings_top, pred_4))

3. Histograms and scatterplots of data and imputations

Observed earnings (excluding 0's)

qplot(earnings_top, data=filter(SIS, !is.na(earnings)), breaks=seq(0,100,10),
      fill=I("white"), col=I("black")) +
    xlab("earnings") + ggtitle("Observed earnings (excluding 0's)")

Random imputation of earnings

qplot(earnings_imp_4, data = filter(SIS, !is.na(earnings)),
     breaks=seq(0,100,10), fill=I("white"), col=I("black")) +
    xlab("earnings") + ggtitle("Random imputation of earnings")

Deterministic imputation scatter plot

ggplot() + 
    geom_point(aes(x=earnings_imp_2, y=earnings_imp_2), shape=19,
               data = filter(SIS, !is.na(earnings))) +
    geom_point(aes(x=pred_2, y=earnings_top), shape=20, color="darkgrey",
               data = filter(SIS, !is.na(earnings))) +
    xlab("Regression prediction") + ylab("Earnings") +
    ggtitle("Deterministic imputation")

Random imputation scatter plot

ggplot() + 
    geom_point(aes(x=earnings_imp_2, y=earnings_imp_4), shape=19,
               data = filter(SIS, !is.na(earnings))) +
    geom_point(aes(x=pred_2, y=earnings_top), shape=20, color="darkgrey",
               data = filter(SIS, !is.na(earnings))) +
    xlab("Regression prediction") + ylab("Earnings") +
    ggtitle("Random imputation")

4. Two-stage imputation model

Fit the 2 models

fit_positive <- stan_glm((earnings>0) ~ male + over65 + white + immig +
  educ_r + any_ssi + any_welfare + any_charity,
  data=SIS, family=binomial(link=logit), refresh=0)
print(fit_positive, digits=2)
stan_glm
 family:       binomial [logit]
 formula:      (earnings > 0) ~ male + over65 + white + immig + educ_r + any_ssi + 
       any_welfare + any_charity
 observations: 1260
 predictors:   9
------
            Median MAD_SD
(Intercept)  0.35   0.26 
male         0.31   0.18 
over65      -4.06   0.34 
white       -0.10   0.21 
immig        0.18   0.18 
educ_r       0.56   0.09 
any_ssi     -2.48   0.38 
any_welfare -0.57   0.35 
any_charity  0.52   0.88 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
# 2nd model was alraedy fitted before
fit_positive_sqrt <- fit_imp_2
print(fit_positive_sqrt)
stan_glm
 family:       gaussian [identity]
 formula:      sqrt(earnings_top) ~ male + over65 + white + immig + educ_r + 
       workmos + workhrs_top + any_ssi + any_welfare + any_charity
 observations: 988
 predictors:   11
------
            Median MAD_SD
(Intercept) -1.7    0.4  
male         0.3    0.1  
over65      -1.4    0.6  
white        1.0    0.1  
immig       -0.6    0.1  
educ_r       0.8    0.1  
workmos      0.3    0.0  
workhrs_top  0.1    0.0  
any_ssi     -0.9    0.5  
any_welfare -1.3    0.4  
any_charity -1.2    0.6  

Auxiliary parameter(s):
      Median MAD_SD
sigma 2.0    0.0   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Predict the sign and then the earnings (if positive)

pred_pos <- posterior_predict(fit_positive, newdata = SIS_predictors, draws = 1)
pred_ifpos_sqrt <- posterior_predict(fit_positive_sqrt, newdata = SIS_predictors, draws = 1)
pred_ifpos <- topcode(pred_ifpos_sqrt^2, 100)
SIS <- mutate(SIS,
 earnings_imp = impute(earnings, pred_pos*pred_ifpos))

5. Iterative regression imputation

Starting values

random_imp <- function (a){
  missing <- is.na(a)
  n_missing <- sum(missing)
  a_obs <- a[!missing]
  imputed <- a
  imputed[missing] <- sample(a_obs, n_missing)
  imputed
}
SIS <- mutate(SIS,
 interest_imp = random_imp(interest),
 earnings_imp = random_imp(earnings))

Simplest regression imputation

n_sims <- 10
for (s in 1:n_sims) {
  # Predict earnings
  output <- capture.output(
      fit <- stan_glm(earnings ~ interest_imp + male + over65 + white +
                            immig + educ_r + workmos + workhrs_top + any_ssi +
                            any_welfare + any_charity,
                        data = SIS, cores = 1, open_progress = FALSE, refresh=0))
  pred1 <- posterior_predict(fit, draws = 1)
  # Impute earnings
  SIS <- mutate(SIS,
    earnings_imp = impute(earnings, pred1))
  # Predict interest
  output <- capture.output(
    fit <- stan_glm(interest ~ earnings_imp + male + over65 + white +
                          immig + educ_r + workmos + workhrs_top + any_ssi +
                          any_welfare + any_charity,
                      data = SIS, cores = 1, open_progress = FALSE, refresh=0))
  pred2 <- posterior_predict(fit, draws = 1)
  # Impute interest
  SIS <- mutate(SIS,
    interest_imp = impute(interest, pred2))
}
LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogSW1wdXRhdGlvbiIKYXV0aG9yOiAiQW5kcmV3IEdlbG1hbiwgSmVubmlmZXIgSGlsbCwgQWtpIFZlaHRhcmkiCmRhdGU6ICJgciBmb3JtYXQoU3lzLkRhdGUoKSlgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIHRoZW1lOiByZWFkYWJsZQogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDIKICAgIHRvY19mbG9hdDogdHJ1ZQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQotLS0KUmVncmVzc2lvbi1iYXNlZCBpbXB1dGF0aW9uIGZvciB0aGUgU29jaWFsIEluZGljYXRvcnMgU3VydmV5LiBTZWUKQ2hhcHRlciAxNyBpbiBSZWdyZXNzaW9uIGFuZCBPdGhlciBTdG9yaWVzLjxicj4KZHBseXIgKyBnZ3Bsb3QgdmVyc2lvbi4gCgotLS0tLS0tLS0tLS0tCgoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRSwgY29tbWVudD1OQSkKIyBzd2l0Y2ggdGhpcyB0byBUUlVFIHRvIHNhdmUgZmlndXJlcyBpbiBzZXBhcmF0ZSBmaWxlcwpzYXZlZmlncyA8LSBGQUxTRQpgYGAKCioqTG9hZCBwYWNrYWdlcyoqCgpgYGB7ciB9CmxpYnJhcnkoInJwcm9qcm9vdCIpCnJvb3Q8LWhhc19maWxlKCIuUk9TLUV4YW1wbGVzLXJvb3QiKSRtYWtlX2ZpeF9maWxlKCkKbGlicmFyeSgicnN0YW5hcm0iKQpsaWJyYXJ5KCJkcGx5ciIpCmxpYnJhcnkoImdncGxvdDIiKQp0aGVtZV9zZXQoYmF5ZXNwbG90Ojp0aGVtZV9kZWZhdWx0KGJhc2VfZmFtaWx5ID0gInNhbnMiKSkKYGBgCgoqKkxvYWQgZGF0YSoqCgpgYGB7ciB9ClNJUyA8LSByZWFkLmNzdihyb290KCJJbXB1dGF0aW9uL2RhdGEiLCJTSVMuY3N2IikpCmhlYWQoU0lTKQpzdW1tYXJ5KFNJUykKYGBgCgoqKkltcHV0YXRpb24gaGVscGVyIGZ1bmN0aW9ucyoqPC9icj4KQ3JlYXRlIGEgY29tcGxldGVkIGRhdGEgdmVjdG9yIHVzaW5nIGltcHV0YXRpb25zCgpgYGB7ciB9CmltcHV0ZSA8LSBmdW5jdGlvbihhLCBhX2ltcHV0ZSkgewogIGlmZWxzZShpcy5uYShhKSwgYV9pbXB1dGUsIGEpCn0KYGBgCgpUb3AgY29kZSBmdW5jdGlvbgoKYGBge3IgfQp0b3Bjb2RlIDwtIGZ1bmN0aW9uKGEsIHRvcCkgewogIGlmZWxzZShhPnRvcCwgdG9wLCBhKQp9CmBgYAoKc2UgbWlzc2luZyBjb2RlcyB0byBOQQoKYGBge3IgfQpuYV9maXggPC0gZnVuY3Rpb24oYSkgewogIGlmZWxzZShhPDAgfCBhPT05OTk5OTksIE5BLCBhKQp9CmlzX2FueSA8LSBmdW5jdGlvbihhKSB7CiAgYW55X2EgPC0gaWZlbHNlKGE+MCwgMSwgMCkKICBhbnlfYVtpcy5uYShhKV0gPC0gMAogIGFueV9hCn0KYGBgCgojIyMgMS4gIERldGVybWluaXN0aWMgaW1wdXRhdGlvbgoqKkltcHV0ZSAwIGVhcm5pbmdzIHVzaW5nIHRoZSBsb2dpY2FsIHJ1bGUgKGlmIHdvcmtlZCAwIG1vbnRocyBhbmQgMCBocnMvd2spKioKCmBgYHtyIH0KU0lTIDwtIG11dGF0ZShTSVMsCiAgemVyb19lYXJuaW5ncyA9ICh3b3JraHJzX3RvcD09MCAmIHdvcmttb3M9PTApLAogIGVhcm5pbmdzX3RvcCA9IGlmX2Vsc2UoemVyb19lYXJuaW5ncywgMCwgdG9wY29kZShlYXJuaW5ncywgMTAwKSkpCmBgYAoKKipDcmVhdGUgYSBsaXR0bGUgZGF0YXNldCB3aXRoIGFsbCBvdXIgcmVkZWZpbmVkIHZhcmlhYmxlcyoqCgpgYGB7ciB9ClNJUyA8LSBzZWxlY3QoU0lTLAogICAgICAgICAgICAgIGVhcm5pbmdzLCBlYXJuaW5nc190b3AsIGludGVyZXN0LCBtYWxlLCBvdmVyNjUsCiAgICAgICAgICAgICAgd2hpdGUsIGltbWlnLCBlZHVjX3IsIHdvcmttb3MsIHdvcmtocnNfdG9wLAogICAgICAgICAgICAgIGFueV9zc2ksIGFueV93ZWxmYXJlLCBhbnlfY2hhcml0eSkKYGBgCgoqKkltcHV0ZSBzdWJzZXQgb2YgZWFybmluZ3MgdGhhdCBhcmUgbm9uemVybzogIGxpbmVhciBzY2FsZSoqCgpgYGB7ciB9CmZpdF9pbXBfMSA8LSBzdGFuX2dsbSgKICBlYXJuaW5ncyB+IG1hbGUgKyBvdmVyNjUgKyB3aGl0ZSArIGltbWlnICsgZWR1Y19yICsKICAgICAgICAgICAgICB3b3JrbW9zICsgd29ya2hyc190b3AgKyBhbnlfc3NpICsKICAgICAgICAgICAgICBhbnlfd2VsZmFyZSArIGFueV9jaGFyaXR5LAogIGRhdGEgPSBTSVMsCiAgc3Vic2V0ID0gZWFybmluZ3MgPiAwLAogIHJlZnJlc2ggPSAwCikKcHJpbnQoZml0X2ltcF8xKQpTSVNfcHJlZGljdG9ycyA8LSBzZWxlY3QoU0lTLCAtc3RhcnRzX3dpdGgoImVhcm5pbmdzIikpCnByZWRfMSA8LSBjb2xNZWFucyhwb3N0ZXJpb3JfbGlucHJlZChmaXRfaW1wXzEsIG5ld2RhdGEgPSBTSVNfcHJlZGljdG9ycykpClNJUyA8LSBtdXRhdGUoU0lTLAogIGVhcm5pbmdzX2ltcF8xID0gaW1wdXRlKGVhcm5pbmdzLCBwcmVkXzEpKQpgYGAKCioqSW1wdXRlIHN1YnNldCBvZiBlYXJuaW5ncyB0aGF0IGFyZSBub256ZXJvOiAgc3F1YXJlIHJvb3Qgc2NhbGUgYW5kIHRvcGNvZGluZyoqCgpgYGB7ciB9CmZpdF9pbXBfMiA8LSBzdGFuX2dsbSgKICBzcXJ0KGVhcm5pbmdzX3RvcCkgfiBtYWxlICsgb3ZlcjY1ICsgd2hpdGUgKyBpbW1pZyArCiAgICAgICAgICAgICAgICAgICAgICAgZWR1Y19yICsgd29ya21vcyArIHdvcmtocnNfdG9wICsgYW55X3NzaSArCiAgICAgICAgICAgICAgICAgICAgICAgYW55X3dlbGZhcmUgKyBhbnlfY2hhcml0eSwKICBkYXRhID0gU0lTLAogIHN1YnNldCA9IGVhcm5pbmdzID4gMCwKICByZWZyZXNoID0gMAopCnByaW50KGZpdF9pbXBfMikKYGBgCgpwb2ludCBwcmVkaWN0aW9ucwoKYGBge3IgfQpwcmVkXzJfc3FydCA8LSBjb2xNZWFucyhwb3N0ZXJpb3JfbGlucHJlZChmaXRfaW1wXzIsIG5ld2RhdGEgPSBTSVNfcHJlZGljdG9ycykpCnByZWRfMiA8LSB0b3Bjb2RlKHByZWRfMl9zcXJ0XjIsIDEwMCkKU0lTIDwtIG11dGF0ZShTSVMsCiAgcHJlZF8yID0gcHJlZF8yLAogIGVhcm5pbmdzX2ltcF8yID0gaW1wdXRlKGVhcm5pbmdzX3RvcCwgcHJlZF8yKSkKYGBgCgojIyMgMi4gIFJhbmRvbSBpbXB1dGF0aW9uCioqTGluZWFyIHNjYWxlICh1c2UgZml0dGVkIG1vZGVsIGZpdF9pbXBfMSkqKgoKYGBge3IgfQpwcmVkXzMgPC0gcG9zdGVyaW9yX3ByZWRpY3QoZml0X2ltcF8xLCBuZXdkYXRhID0gU0lTX3ByZWRpY3RvcnMsIGRyYXdzID0gMSkKU0lTIDwtIG11dGF0ZShTSVMsCiAgZWFybmluZ3NfaW1wXzMgPSBpbXB1dGUoZWFybmluZ3MsIHByZWRfMykpCmBgYAoKKipTcXVhcmUgcm9vdCBzY2FsZSBhbmQgdG9wY29kaW5nICh1c2UgZml0dGVkIG1vZGVsIGZpdF9pbXBfMikqKgoKYGBge3IgfQpwcmVkXzRfc3FydCA8LSBwb3N0ZXJpb3JfcHJlZGljdChmaXRfaW1wXzIsIG5ld2RhdGEgPSBTSVNfcHJlZGljdG9ycywgZHJhd3MgPSAxKQpwcmVkXzQgPC0gdG9wY29kZShwcmVkXzRfc3FydF4yLCAxMDApClNJUyA8LSBtdXRhdGUoU0lTLAogIGVhcm5pbmdzX2ltcF80ID0gaW1wdXRlKGVhcm5pbmdzX3RvcCwgcHJlZF80KSkKYGBgCgojIyMgMy4gIEhpc3RvZ3JhbXMgYW5kIHNjYXR0ZXJwbG90cyBvZiBkYXRhIGFuZCBpbXB1dGF0aW9ucwoKT2JzZXJ2ZWQgZWFybmluZ3MgKGV4Y2x1ZGluZyAwJ3MpCgpgYGB7ciB9CnFwbG90KGVhcm5pbmdzX3RvcCwgZGF0YT1maWx0ZXIoU0lTLCAhaXMubmEoZWFybmluZ3MpKSwgYnJlYWtzPXNlcSgwLDEwMCwxMCksCiAgICAgIGZpbGw9SSgid2hpdGUiKSwgY29sPUkoImJsYWNrIikpICsKICAgIHhsYWIoImVhcm5pbmdzIikgKyBnZ3RpdGxlKCJPYnNlcnZlZCBlYXJuaW5ncyAoZXhjbHVkaW5nIDAncykiKQpgYGAKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBnZ3NhdmUocm9vdCgiSW1wdXRhdGlvbi9maWdzIiwiaW1wdXRlX2hpc3QyX2dnLnBkZiIpLCB3aWR0aCA9IDUsIGhlaWdodCA9IDQpCnFwbG90KGVhcm5pbmdzX2ltcF8yLCBkYXRhID0gZmlsdGVyKFNJUywgIWlzLm5hKGVhcm5pbmdzKSksCiAgICAgYnJlYWtzPXNlcSgwLDEwMCwxMCksIGZpbGw9SSgid2hpdGUiKSwgY29sPUkoImJsYWNrIikpICsKICAgIHhsYWIoImVhcm5pbmdzIikgKyBnZ3RpdGxlKCJEZXRlcm1pbmlzdGljIGltcHV0YXRpb24gb2YgZWFybmluZ3MiKQpgYGAKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBnZ3NhdmUocm9vdCgiSW1wdXRhdGlvbi9maWdzIiwiaW1wdXRlX2hpc3QzX2dnLnBkZiIpLCB3aWR0aCA9IDUsIGhlaWdodCA9IDQpCmBgYAoKUmFuZG9tIGltcHV0YXRpb24gb2YgZWFybmluZ3MKCmBgYHtyIH0KcXBsb3QoZWFybmluZ3NfaW1wXzQsIGRhdGEgPSBmaWx0ZXIoU0lTLCAhaXMubmEoZWFybmluZ3MpKSwKICAgICBicmVha3M9c2VxKDAsMTAwLDEwKSwgZmlsbD1JKCJ3aGl0ZSIpLCBjb2w9SSgiYmxhY2siKSkgKwogICAgeGxhYigiZWFybmluZ3MiKSArIGdndGl0bGUoIlJhbmRvbSBpbXB1dGF0aW9uIG9mIGVhcm5pbmdzIikKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgZ2dzYXZlKHJvb3QoIkltcHV0YXRpb24vZmlncyIsImltcHV0ZV9oaXN0NF9nZy5wZGYiKSwgd2lkdGggPSA1LCBoZWlnaHQgPSA0KQpgYGAKCkRldGVybWluaXN0aWMgaW1wdXRhdGlvbiBzY2F0dGVyIHBsb3QKCmBgYHtyIH0KZ2dwbG90KCkgKyAKICAgIGdlb21fcG9pbnQoYWVzKHg9ZWFybmluZ3NfaW1wXzIsIHk9ZWFybmluZ3NfaW1wXzIpLCBzaGFwZT0xOSwKICAgICAgICAgICAgICAgZGF0YSA9IGZpbHRlcihTSVMsICFpcy5uYShlYXJuaW5ncykpKSArCiAgICBnZW9tX3BvaW50KGFlcyh4PXByZWRfMiwgeT1lYXJuaW5nc190b3ApLCBzaGFwZT0yMCwgY29sb3I9ImRhcmtncmV5IiwKICAgICAgICAgICAgICAgZGF0YSA9IGZpbHRlcihTSVMsICFpcy5uYShlYXJuaW5ncykpKSArCiAgICB4bGFiKCJSZWdyZXNzaW9uIHByZWRpY3Rpb24iKSArIHlsYWIoIkVhcm5pbmdzIikgKwogICAgZ2d0aXRsZSgiRGV0ZXJtaW5pc3RpYyBpbXB1dGF0aW9uIikKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgZ2dzYXZlKHJvb3QoIkltcHV0YXRpb24vZmlncyIsImltcHV0ZV9zY2F0MV9nZy5wZGYiKSwgd2lkdGggPSA1LCBoZWlnaHQgPSA0KQpgYGAKClJhbmRvbSBpbXB1dGF0aW9uIHNjYXR0ZXIgcGxvdAoKYGBge3IgfQpnZ3Bsb3QoKSArIAogICAgZ2VvbV9wb2ludChhZXMoeD1lYXJuaW5nc19pbXBfMiwgeT1lYXJuaW5nc19pbXBfNCksIHNoYXBlPTE5LAogICAgICAgICAgICAgICBkYXRhID0gZmlsdGVyKFNJUywgIWlzLm5hKGVhcm5pbmdzKSkpICsKICAgIGdlb21fcG9pbnQoYWVzKHg9cHJlZF8yLCB5PWVhcm5pbmdzX3RvcCksIHNoYXBlPTIwLCBjb2xvcj0iZGFya2dyZXkiLAogICAgICAgICAgICAgICBkYXRhID0gZmlsdGVyKFNJUywgIWlzLm5hKGVhcm5pbmdzKSkpICsKICAgIHhsYWIoIlJlZ3Jlc3Npb24gcHJlZGljdGlvbiIpICsgeWxhYigiRWFybmluZ3MiKSArCiAgICBnZ3RpdGxlKCJSYW5kb20gaW1wdXRhdGlvbiIpCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIGdnc2F2ZShyb290KCJJbXB1dGF0aW9uL2ZpZ3MiLCJpbXB1dGVfc2NhdDJfZ2cucGRmIiksIHdpZHRoID0gNSwgaGVpZ2h0ID0gNCkKYGBgCgojIyMgNC4gVHdvLXN0YWdlIGltcHV0YXRpb24gbW9kZWwKKipGaXQgdGhlIDIgbW9kZWxzKioKCmBgYHtyIH0KZml0X3Bvc2l0aXZlIDwtIHN0YW5fZ2xtKChlYXJuaW5ncz4wKSB+IG1hbGUgKyBvdmVyNjUgKyB3aGl0ZSArIGltbWlnICsKICBlZHVjX3IgKyBhbnlfc3NpICsgYW55X3dlbGZhcmUgKyBhbnlfY2hhcml0eSwKICBkYXRhPVNJUywgZmFtaWx5PWJpbm9taWFsKGxpbms9bG9naXQpLCByZWZyZXNoPTApCnByaW50KGZpdF9wb3NpdGl2ZSwgZGlnaXRzPTIpCiMgMm5kIG1vZGVsIHdhcyBhbHJhZWR5IGZpdHRlZCBiZWZvcmUKZml0X3Bvc2l0aXZlX3NxcnQgPC0gZml0X2ltcF8yCnByaW50KGZpdF9wb3NpdGl2ZV9zcXJ0KQpgYGAKCioqUHJlZGljdCB0aGUgc2lnbiBhbmQgdGhlbiB0aGUgZWFybmluZ3MgKGlmIHBvc2l0aXZlKSoqCgpgYGB7ciB9CnByZWRfcG9zIDwtIHBvc3Rlcmlvcl9wcmVkaWN0KGZpdF9wb3NpdGl2ZSwgbmV3ZGF0YSA9IFNJU19wcmVkaWN0b3JzLCBkcmF3cyA9IDEpCnByZWRfaWZwb3Nfc3FydCA8LSBwb3N0ZXJpb3JfcHJlZGljdChmaXRfcG9zaXRpdmVfc3FydCwgbmV3ZGF0YSA9IFNJU19wcmVkaWN0b3JzLCBkcmF3cyA9IDEpCnByZWRfaWZwb3MgPC0gdG9wY29kZShwcmVkX2lmcG9zX3NxcnReMiwgMTAwKQpTSVMgPC0gbXV0YXRlKFNJUywKIGVhcm5pbmdzX2ltcCA9IGltcHV0ZShlYXJuaW5ncywgcHJlZF9wb3MqcHJlZF9pZnBvcykpCmBgYAoKIyMjIDUuICBJdGVyYXRpdmUgcmVncmVzc2lvbiBpbXB1dGF0aW9uCioqU3RhcnRpbmcgdmFsdWVzKioKCmBgYHtyIH0KcmFuZG9tX2ltcCA8LSBmdW5jdGlvbiAoYSl7CiAgbWlzc2luZyA8LSBpcy5uYShhKQogIG5fbWlzc2luZyA8LSBzdW0obWlzc2luZykKICBhX29icyA8LSBhWyFtaXNzaW5nXQogIGltcHV0ZWQgPC0gYQogIGltcHV0ZWRbbWlzc2luZ10gPC0gc2FtcGxlKGFfb2JzLCBuX21pc3NpbmcpCiAgaW1wdXRlZAp9ClNJUyA8LSBtdXRhdGUoU0lTLAogaW50ZXJlc3RfaW1wID0gcmFuZG9tX2ltcChpbnRlcmVzdCksCiBlYXJuaW5nc19pbXAgPSByYW5kb21faW1wKGVhcm5pbmdzKSkKYGBgCgoqKlNpbXBsZXN0IHJlZ3Jlc3Npb24gaW1wdXRhdGlvbioqCgpgYGB7ciB9Cm5fc2ltcyA8LSAxMApmb3IgKHMgaW4gMTpuX3NpbXMpIHsKICAjIFByZWRpY3QgZWFybmluZ3MKICBvdXRwdXQgPC0gY2FwdHVyZS5vdXRwdXQoCiAgICAgIGZpdCA8LSBzdGFuX2dsbShlYXJuaW5ncyB+IGludGVyZXN0X2ltcCArIG1hbGUgKyBvdmVyNjUgKyB3aGl0ZSArCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBpbW1pZyArIGVkdWNfciArIHdvcmttb3MgKyB3b3JraHJzX3RvcCArIGFueV9zc2kgKwogICAgICAgICAgICAgICAgICAgICAgICAgICAgYW55X3dlbGZhcmUgKyBhbnlfY2hhcml0eSwKICAgICAgICAgICAgICAgICAgICAgICAgZGF0YSA9IFNJUywgY29yZXMgPSAxLCBvcGVuX3Byb2dyZXNzID0gRkFMU0UsIHJlZnJlc2g9MCkpCiAgcHJlZDEgPC0gcG9zdGVyaW9yX3ByZWRpY3QoZml0LCBkcmF3cyA9IDEpCiAgIyBJbXB1dGUgZWFybmluZ3MKICBTSVMgPC0gbXV0YXRlKFNJUywKICAgIGVhcm5pbmdzX2ltcCA9IGltcHV0ZShlYXJuaW5ncywgcHJlZDEpKQogICMgUHJlZGljdCBpbnRlcmVzdAogIG91dHB1dCA8LSBjYXB0dXJlLm91dHB1dCgKICAgIGZpdCA8LSBzdGFuX2dsbShpbnRlcmVzdCB+IGVhcm5pbmdzX2ltcCArIG1hbGUgKyBvdmVyNjUgKyB3aGl0ZSArCiAgICAgICAgICAgICAgICAgICAgICAgICAgaW1taWcgKyBlZHVjX3IgKyB3b3JrbW9zICsgd29ya2hyc190b3AgKyBhbnlfc3NpICsKICAgICAgICAgICAgICAgICAgICAgICAgICBhbnlfd2VsZmFyZSArIGFueV9jaGFyaXR5LAogICAgICAgICAgICAgICAgICAgICAgZGF0YSA9IFNJUywgY29yZXMgPSAxLCBvcGVuX3Byb2dyZXNzID0gRkFMU0UsIHJlZnJlc2g9MCkpCiAgcHJlZDIgPC0gcG9zdGVyaW9yX3ByZWRpY3QoZml0LCBkcmF3cyA9IDEpCiAgIyBJbXB1dGUgaW50ZXJlc3QKICBTSVMgPC0gbXV0YXRlKFNJUywKICAgIGludGVyZXN0X2ltcCA9IGltcHV0ZShpbnRlcmVzdCwgcHJlZDIpKQp9CmBgYAoK