Simulating Research Designs for Social Science pt 2
DesignLibrary - Beginner Guide & Case Study
John-Henry Pezzuto
This is part two of a DeclareDesign tutorial series. You can find part 1 here.
DesignLibrary is a online repository of ready-made designs for experimental simulating using DeclareDesign. It is basically DeclareDesign with training wheels, which makes it a great jumping off point if you already know what type of experimental design you want to use. It also has really great documentation! I’ve found it especially useful if you want to plan or replicate a simple design. You can find a full list of the designs here. Here are a few designs I selected arbitrarily:
- Binary Instrumental Variable -
binary_iv_designer - Block-Cluster-Two-Arm Design -
block_cluster_two_arm_designer - Cluster Sampling -
cluster_sampling_designer - Mediation Analysis Design -
mediation_analysis_designer - Pretest Post-test Design -
pretest_posttest_designer - Randomized Response -
randomized_response_designer - Two-by-Two Factorial -
two_by_two_designer
Trying it out
Here they we are interested in replicating the results from the famous Brescoll & Uhlmann, 2008 on gender, emotion, and status conferrel (the amount of status you think someone deserves to have).
The high level summary of this study is that men who get angry in the workplace are conferred more status than men who get sad in the workplace. This is kind of intuitive, right? People who get angry in the workplace are like bosses of course they’re high status! However, this study found that when women get angry in the workplace they are conferred less status compared to sad women.
This sets us up for a 2 (gender: male, female) x 2 (emotion: anger sadness) study design where:
| 0 | 1 | ||
|---|---|---|---|
| B | 0 | male anger | male sadness |
| 1 | female anger | female sadness |
To fill the above tables with numbers, we can take the results directly from Brescoll & Uhlmann, 2008 where they report their findings in Table 1 on page 270. We can plug their findings for status conferral directly into our table.
| 0 | 1 | ||
|---|---|---|---|
| B | 0 | 6.47 | 4.05 |
| 1 | 3.75 | 5.02 |
Enter DesignLibrary
Now we can begin to start thinking about simulation! Here we plug-in the above table into the two_by_two_designer function from DesignLibrary, and uhh that’s basically it! A lot easier than DeclareDesign right? We can learn more about the arguments for two_by_two_designer by using ?two_by_two_designer
library(tidyverse)
library(DesignLibrary)
library(DeclareDesign)
brescoll_uhlmann <-
two_by_two_designer(mean_A0B0 = 6.47, # equivalent to 'outcome_means = c(6.47, 3.75, 4.05, 5.02)'
mean_A0B1 = 3.75,
mean_A1B0 = 4.05,
mean_A1B1 = 5.02,
outcome_sds = c(2.25, 1.77, 1.61, 1.8))How many participants do we need?
The next goal is explore the right N for our design. We can use the above brescoll_uhlmann together with redesign function to evaluate the design with different N’s. I could’ve done more granular increments of 50, but I use 100 to help lower compiling time.
brescoll_uhlmann_sims <-
brescoll_uhlmann %>%
redesign(N = seq(100, 2000, by = 100)) %>% # from 100 to 2000 in increments of 100
diagnose_design() %>%
get_diagnosands()Alternatively we can use the more elegant and generalized approach with expand_design and the two_by_two_designer base function. This is recommended by the folks at DeclareDesign because it is considered safer code to run.
brescoll_uhlmann_sims <-
two_by_two_designer %>%
expand_design(mean_A0B0 = 6.47,
mean_A0B1 = 3.75,
mean_A1B0 = 4.05,
mean_A1B1 = 5.02,
outcome_sds = list(c(2.25, 1.77, 1.61, 1.8)),
N = seq(100, 2000, by = 100)) %>%
diagnose_design() %>%
get_diagnosands()Here we can take an ugly peek at our data. Each different N we evaluated takes up 3 different rows. A row evaluating the “A”, “B”, “A:B” aka “emotion”, “gender” and the “interaction”. Because we are replicating a previous study, they did not change the means in the assumption.
glimpse(brescoll_uhlmann_sims)Observations: 60
Variables: 24
$ design_label <chr> "design_1", "design_1", "design_1", "design_…
$ N <dbl> 100, 100, 100, 200, 200, 200, 300, 300, 300,…
$ estimand_label <chr> "ate_A", "ate_B", "interaction", "ate_A", "a…
$ estimator_label <chr> "No_Interaction", "No_Interaction", "Interac…
$ term <chr> "A", "B", "A:B", "A", "B", "A:B", "A", "B", …
$ bias <dbl> -0.0190363867, -0.0056538398, 0.0451473134, …
$ `se(bias)` <dbl> 0.014596484, 0.017010385, 0.032738626, 0.010…
$ rmse <dbl> 0.3779281, 0.3805540, 0.7541199, 0.2603786, …
$ `se(rmse)` <dbl> 0.012729110, 0.012772565, 0.022977813, 0.008…
$ power <dbl> 0.228, 0.492, 0.990, 0.430, 0.756, 1.000, 0.…
$ `se(power)` <dbl> 0.018172284, 0.022759835, 0.004664762, 0.021…
$ coverage <dbl> 0.990, 0.986, 0.974, 0.986, 0.982, 0.964, 0.…
$ `se(coverage)` <dbl> 0.004939390, 0.005572624, 0.006692646, 0.005…
$ mean_estimate <dbl> -0.6058400, -0.8858187, 3.7354927, -0.604670…
$ `se(mean_estimate)` <dbl> 0.016631210, 0.020047839, 0.038163507, 0.012…
$ sd_estimate <dbl> 0.4312540, 0.4298720, 0.8602923, 0.2856536, …
$ `se(sd_estimate)` <dbl> 0.014984936, 0.014688257, 0.024194064, 0.009…
$ mean_se <dbl> 0.4631760, 0.4631760, 0.8454376, 0.3273779, …
$ `se(mean_se)` <dbl> 0.0014485996, 0.0014485996, 0.0026268929, 0.…
$ type_s_rate <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ `se(type_s_rate)` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ mean_estimand <dbl> -0.5868036, -0.8801648, 3.6903453, -0.581334…
$ `se(mean_estimand)` <dbl> 0.008287417, 0.008452199, 0.016137426, 0.005…
$ n_sims <dbl> 500, 500, 500, 500, 500, 500, 500, 500, 500,…Visualize our results
And next we visualize to see when it will be well powered!
brescoll_uhlmann_sims %>%
mutate(estimand_label = case_when(estimand_label == "ate_A" ~ "Emotion",
estimand_label == "ate_B" ~ "Gender",
estimand_label == "interaction" ~ "Emotion x Gender")) %>%
ggplot(aes(N, power,
col = estimand_label,
group = estimand_label,
ymin = power - 1.96*`se(power)`,
ymax = power + 1.96*`se(power)`)) +
geom_smooth(position = "dodge", stat = "identity") +
geom_hline(yintercept = 0.8) +
scale_y_continuous(labels = scales::percent) +
theme_bw()
We see that this design will need about 500 participants to be well-powered or about 125 in each of the four cells. Surprisingly, we see that the interaction is actually the largest effect. Next we’re interested in what N exactly do emotion and gender reach the power the 80% power mark.
What is the minimum number of participants needed?
brescoll_uhlmann_sims %>%
mutate(estimand_label = case_when(estimand_label == "ate_A" ~ "Emotion",
estimand_label == "ate_B" ~ "Gender",
estimand_label == "interaction" ~ "Emotion x Gender")) %>%
group_by(estimand_label) %>%
filter(power > .8) %>%
filter(N == min(N)) %>%
select(estimand_label, term, N, power, `se(power)`) %>%
kable() %>%
kable_styling(position = "center")| estimand_label | term | N | power | se(power) |
|---|---|---|---|---|
| Emotion x Gender | A:B | 100 | 0.990 | 0.0046648 |
| Gender | B | 300 | 0.926 | 0.0116498 |
| Emotion | A | 500 | 0.830 | 0.0168127 |
Stealing the DesignLibrary code
Lastly I wanted to note how you could customize DesignLibrary code for your own designs, let’s say you found a code in DesignLibrary that is almost right for your needs, but not quite. To see the internals in a DesignLibrary function use get_design_code.
In order to take the code for a copy pasting you, can cat(paste(code, collapse = "\n")). This is equivalent to print_table in the blpl package. Note that the ##’s won’t show up in your R console.
twoxtwo_code <-
get_design_code(two_by_two_designer()) # generalized 2x2, note the parenthesis on the inside function
bu_code <- get_design_code(brescoll_uhlmann) # the example we just made
cat(paste(bu_code, collapse = "\n"))N <- 100
prob_A <- 0.5
prob_B <- 0.5
weight_A <- 0.5
weight_B <- 0.5
mean_A0B0 <- 6.47
mean_A0B1 <- 3.75
mean_A1B0 <- 4.05
mean_A1B1 <- 5.02
sd_i <- 1
outcome_sds <- c(2.25, 1.77, 1.61, 1.8)
population <- declare_population(N, u = rnorm(N, sd = sd_i))
potential_outcomes <- declare_potential_outcomes(Y_A_0_B_0 = mean_A0B0 +
u + rnorm(N, sd = outcome_sds[1]), Y_A_0_B_1 = mean_A0B1 +
u + rnorm(N, sd = outcome_sds[2]), Y_A_1_B_0 = mean_A1B0 +
u + rnorm(N, sd = outcome_sds[3]), Y_A_1_B_1 = mean_A1B1 +
u + rnorm(N, sd = outcome_sds[4]))
estimand_1 <- declare_estimand(ate_A = weight_B * mean(Y_A_1_B_1 -
Y_A_0_B_1) + (1 - weight_B) * mean(Y_A_1_B_0 - Y_A_0_B_0))
estimand_2 <- declare_estimand(ate_B = weight_A * mean(Y_A_1_B_1 -
Y_A_1_B_0) + (1 - weight_A) * mean(Y_A_0_B_1 - Y_A_0_B_0))
estimand_3 <- declare_estimand(interaction = mean((Y_A_1_B_1 -
Y_A_1_B_0) - (Y_A_0_B_1 - Y_A_0_B_0)))
assign_A <- declare_assignment(prob = prob_A, assignment_variable = A)
assign_B <- declare_assignment(prob = prob_B, assignment_variable = B,
blocks = A)
reveal_Y <- declare_reveal(Y_variables = Y, assignment_variables = c(A,
B))
estimator_1 <- declare_estimator(Y ~ A + B, model = lm_robust,
term = c("A", "B"), estimand = c("ate_A", "ate_B"), label = "No_Interaction")
estimator_2 <- declare_estimator(Y ~ A + B + A:B, model = lm_robust,
term = "A:B", estimand = "interaction", label = "Interaction")
two_by_two_design <- population + potential_outcomes + estimand_1 +
estimand_2 + estimand_3 + assign_A + assign_B + reveal_Y +
estimator_1 + estimator_2LS0tCnRpdGxlOiAiU2ltdWxhdGluZyBSZXNlYXJjaCBEZXNpZ25zIGZvciBTb2NpYWwgU2NpZW5jZSBwdCAyIgpzdWJ0aXRsZTogIkRlc2lnbkxpYnJhcnkgLSBCZWdpbm5lciBHdWlkZSAmIENhc2UgU3R1ZHkiCmF1dGhvcjogIkpvaG4tSGVucnkgUGV6enV0byIKb3V0cHV0OiAKICBybWFya2Rvd246Omh0bWxfbm90ZWJvb2s6CiAgICBzbWFydDogZmFsc2UKICAgIHRoZW1lOiBzcGFjZWxhYgogICAgbnVtYmVyX3NlY3Rpb25zOiBubwogICAgdG9jOiB0cnVlCiAgICB0b2NfZmxvYXQ6CiAgICAgIGNvbGxhcHNlZDogZmFsc2UKICAgIHRvY19kZXB0aDogNAotLS0KCgpUaGlzIGlzIHBhcnQgdHdvIG9mIGEgRGVjbGFyZURlc2lnbiBbdHV0b3JpYWwgc2VyaWVzXShodHRwczovL3JwdWJzLmNvbS9qb2huLWhlbnJ5KS4gWW91IGNhbiBbZmluZCBwYXJ0IDEgaGVyZV0oaHR0cHM6Ly9ycHVicy5jb20vam9obi1oZW5yeS9kZWNsYXJlLWRlc2lnbi0xKS4KCltEZXNpZ25MaWJyYXJ5XShodHRwczovL2RlY2xhcmVkZXNpZ24ub3JnL2xpYnJhcnkvKSBpcyBhIG9ubGluZSByZXBvc2l0b3J5IG9mIHJlYWR5LW1hZGUgZGVzaWducyBmb3IgZXhwZXJpbWVudGFsIHNpbXVsYXRpbmcgdXNpbmcgRGVjbGFyZURlc2lnbi4gSXQgaXMgYmFzaWNhbGx5IERlY2xhcmVEZXNpZ24gd2l0aCB0cmFpbmluZyB3aGVlbHMsIHdoaWNoIG1ha2VzIGl0IGEgZ3JlYXQganVtcGluZyBvZmYgcG9pbnQgaWYgeW91IGFscmVhZHkga25vdyB3aGF0IHR5cGUgb2YgZXhwZXJpbWVudGFsIGRlc2lnbiB5b3Ugd2FudCB0byB1c2UuIEl0IGFsc28gaGFzIHJlYWxseSBncmVhdCBkb2N1bWVudGF0aW9uISBJJ3ZlIGZvdW5kIGl0IGVzcGVjaWFsbHkgdXNlZnVsIGlmIHlvdSB3YW50IHRvIHBsYW4gb3IgcmVwbGljYXRlIGEgc2ltcGxlIGRlc2lnbi4gWW91IGNhbiBmaW5kIGEgZnVsbCBsaXN0IG9mIHRoZSBkZXNpZ25zIFtoZXJlXShodHRwczovL2RlY2xhcmVkZXNpZ24ub3JnL2xpYnJhcnkvKS4gSGVyZSBhcmUgYSBmZXcgZGVzaWducyBJIHNlbGVjdGVkIGFyYml0cmFyaWx5OgoKKiBbQmluYXJ5IEluc3RydW1lbnRhbCBWYXJpYWJsZV0oaHR0cHM6Ly9kZWNsYXJlZGVzaWduLm9yZy9saWJyYXJ5L2FydGljbGVzL2JpbmFyeV9pdi5odG1sKSAtIGBiaW5hcnlfaXZfZGVzaWduZXJgCiogW0Jsb2NrLUNsdXN0ZXItVHdvLUFybSBEZXNpZ25dKGh0dHBzOi8vZGVjbGFyZWRlc2lnbi5vcmcvbGlicmFyeS9hcnRpY2xlcy9ibG9ja19jbHVzdGVyX3R3b19hcm0uaHRtbCkgLSBgYmxvY2tfY2x1c3Rlcl90d29fYXJtX2Rlc2lnbmVyYAoqIFtDbHVzdGVyIFNhbXBsaW5nXShodHRwczovL2RlY2xhcmVkZXNpZ24ub3JnL2xpYnJhcnkvYXJ0aWNsZXMvY2x1c3Rlcl9zYW1wbGluZy5odG1sKSAtIGBjbHVzdGVyX3NhbXBsaW5nX2Rlc2lnbmVyYAoqIFtNZWRpYXRpb24gQW5hbHlzaXMgRGVzaWduXShodHRwczovL2RlY2xhcmVkZXNpZ24ub3JnL2xpYnJhcnkvYXJ0aWNsZXMvbWVkaWF0aW9uX2FuYWx5c2lzLmh0bWwpIC0gYG1lZGlhdGlvbl9hbmFseXNpc19kZXNpZ25lcmAKKiBbUHJldGVzdCBQb3N0LXRlc3QgRGVzaWduXShodHRwczovL2RlY2xhcmVkZXNpZ24ub3JnL2xpYnJhcnkvYXJ0aWNsZXMvcHJldGVzdF9wb3N0dGVzdC5odG1sKSAtIGBwcmV0ZXN0X3Bvc3R0ZXN0X2Rlc2lnbmVyYAoqIFtSYW5kb21pemVkIFJlc3BvbnNlXShodHRwczovL2RlY2xhcmVkZXNpZ24ub3JnL2xpYnJhcnkvYXJ0aWNsZXMvcmFuZG9taXplZF9yZXNwb25zZS5odG1sKSAtIGByYW5kb21pemVkX3Jlc3BvbnNlX2Rlc2lnbmVyYAoqIFtUd28tYnktVHdvIEZhY3RvcmlhbF0oaHR0cHM6Ly9kZWNsYXJlZGVzaWduLm9yZy9saWJyYXJ5L2FydGljbGVzL3R3b19ieV90d28uaHRtbCkgLSBgdHdvX2J5X3R3b19kZXNpZ25lcmAKCgojIyMgVHJ5aW5nIGl0IG91dAoKSGVyZSB0aGV5IHdlIGFyZSBpbnRlcmVzdGVkIGluIHJlcGxpY2F0aW5nIHRoZSByZXN1bHRzIGZyb20gdGhlIGZhbW91cyBbQnJlc2NvbGwgJiBVaGxtYW5uLCAyMDA4XShodHRwczovL3d3dy5uY2JpLm5sbS5uaWguZ292L3B1Ym1lZC8xODMxNTgwMCkgb24gZ2VuZGVyLCBlbW90aW9uLCBhbmQgc3RhdHVzIGNvbmZlcnJlbCAodGhlIGFtb3VudCBvZiBzdGF0dXMgeW91IHRoaW5rIHNvbWVvbmUgZGVzZXJ2ZXMgdG8gaGF2ZSkuIAoKVGhlIGhpZ2ggbGV2ZWwgc3VtbWFyeSBvZiB0aGlzIHN0dWR5IGlzIHRoYXQgbWVuIHdobyBnZXQgYW5ncnkgaW4gdGhlIHdvcmtwbGFjZSBhcmUgY29uZmVycmVkIG1vcmUgc3RhdHVzIHRoYW4gbWVuIHdobyBnZXQgc2FkIGluIHRoZSB3b3JrcGxhY2UuIFRoaXMgaXMga2luZCBvZiBpbnR1aXRpdmUsIHJpZ2h0PyBQZW9wbGUgd2hvIGdldCBhbmdyeSBpbiB0aGUgd29ya3BsYWNlIGFyZSBsaWtlIGJvc3NlcyBvZiBjb3Vyc2UgdGhleSdyZSBoaWdoIHN0YXR1cyEgSG93ZXZlciwgdGhpcyBzdHVkeSBmb3VuZCB0aGF0IHdoZW4gd29tZW4gZ2V0IGFuZ3J5IGluIHRoZSB3b3JrcGxhY2UgdGhleSBhcmUgY29uZmVycmVkIF9sZXNzXyBzdGF0dXMgY29tcGFyZWQgdG8gc2FkIHdvbWVuLgoKVGhpcyBzZXRzIHVzIHVwIGZvciBhIDIgKGdlbmRlcjogbWFsZSwgZmVtYWxlKSB4IDIgKGVtb3Rpb246IGFuZ2VyIHNhZG5lc3MpIHN0dWR5IGRlc2lnbiB3aGVyZToKYGBge3IgZnJhbWVfbnVtYmVycywgZWNobz1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShrbml0cikKbGlicmFyeShrYWJsZUV4dHJhKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKCnRpYmJsZShsYWJlbCA9IGMoIkIiLCAiQiIpLAogICAgICAgbnVtID0gYygwLCAxKSwKICAgICAgIGAwYCA9IGMoIm1hbGUgYW5nZXIiLCAiZmVtYWxlIGFuZ2VyIiksCiAgICAgICBgMWAgPSBjKCJtYWxlIHNhZG5lc3MiLCAiZmVtYWxlIHNhZG5lc3MiKSkgJT4lIAogIGthYmxlKGFsaWduID0gJ2MnLCBjb2wubmFtZXMgPSBjKCIgIiwiICIsICIwIiwgIjEiKSwgbGluZXNlcCA9ICdcXGFkZGxpbmVzcGFjZSBcXGFkZGxpbmVzcGFjZScpICU+JSAKICBrYWJsZV9zdHlsaW5nKHBvc2l0aW9uID0gImNlbnRlciIpICU+JSAKICBhZGRfaGVhZGVyX2Fib3ZlKGMoIiAiID0gMiwgIkEiID0gMiksIGxpbmUgPSBUUlVFLCBleHRyYV9jc3MgPSApICU+JSAKICBjb2x1bW5fc3BlYyhjb2x1bW4gPSAxLCBib2xkID0gVCwgd2lkdGggPSAiMWVtIiwgYm9yZGVyX3JpZ2h0ID0gRiwgCiAgICAgICAgICAgICAgZXh0cmFfY3NzID0gImJvcmRlci1yaWdodDoxLjNweCBzb2xpZCBsaWdodGdyZXk7IikgJT4lIAogIGNvbHVtbl9zcGVjKGNvbHVtbiA9IDIsIGJvbGQgPSBULCB3aWR0aCA9ICIxZW0iLCBib3JkZXJfcmlnaHQgPSBGKSAlPiUgCiAgY29sdW1uX3NwZWMoY29sdW1uID0gMywgd2lkdGggPSAiN2VtIiwgYm9yZGVyX3JpZ2h0ID0gRikgJT4lIAogIGNvbHVtbl9zcGVjKGNvbHVtbiA9IDQsIHdpZHRoID0gIjdlbSIsIGJvcmRlcl9yaWdodCA9IEYpICU+JQogIHJvd19zcGVjKHJvdyA9IDEpICU+JSAKICBjb2xsYXBzZV9yb3dzKGNvbHVtbnMgPSAxKQpgYGAKClRvIGZpbGwgdGhlIGFib3ZlIHRhYmxlcyB3aXRoIG51bWJlcnMsIHdlIGNhbiB0YWtlIHRoZSByZXN1bHRzIGRpcmVjdGx5IGZyb20gQnJlc2NvbGwgJiBVaGxtYW5uLCAyMDA4IHdoZXJlIHRoZXkgcmVwb3J0IHRoZWlyIGZpbmRpbmdzIGluIFRhYmxlIDEgb24gcGFnZSAyNzAuIFdlIGNhbiBwbHVnIHRoZWlyIGZpbmRpbmdzIGZvciBzdGF0dXMgY29uZmVycmFsIGRpcmVjdGx5IGludG8gb3VyIHRhYmxlLgoKCiFbXSh0YWJsZV8xX2Rlc2lnbl9saWJyYXJ5LnBuZykKCgpgYGB7ciwgcGx1Z2luX251bWJlcnMsIGVjaG89RkFMU0V9CnRpYmJsZShsYWJlbCA9IGMoIkIiLCAiQiIpLAogICAgICAgbnVtID0gYygwLCAxKSwKICAgICAgIGAwYCA9IGMoNi40NywgMy43NSksCiAgICAgICBgMWAgPSBjKCA0LjA1LCA1LjAyKSkgJT4lIAogIGthYmxlKGFsaWduID0gJ2MnLCAsIGNvbC5uYW1lcyA9IGMoIiAiLCIgIiwgIjAiLCAiMSIpLCBsaW5lc2VwID0gJ1xcYWRkbGluZXNwYWNlIFxcYWRkbGluZXNwYWNlJykgJT4lCiAga2FibGVfc3R5bGluZyhwb3NpdGlvbiA9ICJjZW50ZXIiKSAlPiUgCiAgYWRkX2hlYWRlcl9hYm92ZShjKCIgIiA9IDIsICJBIiA9IDIpLCBsaW5lID0gVFJVRSwgZXh0cmFfY3NzID0gKSAlPiUgCiAgY29sdW1uX3NwZWMoY29sdW1uID0gMSwgYm9sZCA9IFQsIHdpZHRoID0gIjFlbSIsIGJvcmRlcl9yaWdodCA9IEYsIAogICAgICAgICAgICAgIGV4dHJhX2NzcyA9ICJib3JkZXItcmlnaHQ6MS4zcHggc29saWQgbGlnaHRncmV5OyIpICU+JSAKICBjb2x1bW5fc3BlYyhjb2x1bW4gPSAyLCBib2xkID0gVCwgd2lkdGggPSAiMWVtIiwgYm9yZGVyX3JpZ2h0ID0gRikgJT4lIAogIGNvbHVtbl9zcGVjKGNvbHVtbiA9IDMsIHdpZHRoID0gIjdlbSIsIGJvcmRlcl9yaWdodCA9IEYpICU+JSAKICBjb2x1bW5fc3BlYyhjb2x1bW4gPSA0LCB3aWR0aCA9ICI3ZW0iLCBib3JkZXJfcmlnaHQgPSBGKSAlPiUKICByb3dfc3BlYyhyb3cgPSAxKSAlPiUgCiAgY29sbGFwc2Vfcm93cyhjb2x1bW5zID0gMSkKYGBgCgoKIyBFbnRlciBEZXNpZ25MaWJyYXJ5CgpOb3cgd2UgY2FuIGJlZ2luIHRvIHN0YXJ0IHRoaW5raW5nIGFib3V0IHNpbXVsYXRpb24hIEhlcmUgd2UgcGx1Zy1pbiB0aGUgYWJvdmUgdGFibGUgaW50byB0aGUgYHR3b19ieV90d29fZGVzaWduZXJgIGZ1bmN0aW9uIGZyb20gRGVzaWduTGlicmFyeSwgYW5kIHVoaCB0aGF0J3MgYmFzaWNhbGx5IGl0ISBBIGxvdCBlYXNpZXIgdGhhbiBEZWNsYXJlRGVzaWduIHJpZ2h0PyBXZSBjYW4gbGVhcm4gbW9yZSBhYm91dCB0aGUgYXJndW1lbnRzIGZvciBgdHdvX2J5X3R3b19kZXNpZ25lcmAgYnkgdXNpbmcgYD90d29fYnlfdHdvX2Rlc2lnbmVyYAoKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KERlc2lnbkxpYnJhcnkpCmxpYnJhcnkoRGVjbGFyZURlc2lnbikKCmJyZXNjb2xsX3VobG1hbm4gPC0gCiAgdHdvX2J5X3R3b19kZXNpZ25lcihtZWFuX0EwQjAgPSA2LjQ3LCAgIyBlcXVpdmFsZW50IHRvICdvdXRjb21lX21lYW5zID0gYyg2LjQ3LCAzLjc1LCA0LjA1LCA1LjAyKScKICAgICAgICAgICAgICAgICAgICAgIG1lYW5fQTBCMSA9IDMuNzUsIAogICAgICAgICAgICAgICAgICAgICAgbWVhbl9BMUIwID0gNC4wNSwgCiAgICAgICAgICAgICAgICAgICAgICBtZWFuX0ExQjEgPSA1LjAyLAogICAgICAgICAgICAgICAgICAgICAgb3V0Y29tZV9zZHMgPSBjKDIuMjUsIDEuNzcsIDEuNjEsIDEuOCkpCmBgYAoKCiMjIEhvdyBtYW55IHBhcnRpY2lwYW50cyBkbyB3ZSBuZWVkPwoKVGhlIG5leHQgZ29hbCBpcyBleHBsb3JlIHRoZSByaWdodCBOIGZvciBvdXIgZGVzaWduLiBXZSBjYW4gdXNlIHRoZSBhYm92ZSBgYnJlc2NvbGxfdWhsbWFubmAgdG9nZXRoZXIgd2l0aCBgcmVkZXNpZ25gIGZ1bmN0aW9uIHRvIGV2YWx1YXRlIHRoZSBkZXNpZ24gd2l0aCBkaWZmZXJlbnQgTidzLiBJIGNvdWxkJ3ZlIGRvbmUgbW9yZSBncmFudWxhciBpbmNyZW1lbnRzIG9mIDUwLCBidXQgSSB1c2UgMTAwIHRvIGhlbHAgbG93ZXIgY29tcGlsaW5nIHRpbWUuCmBgYHtyIHJlZGVzaWduLCBldmFsPUZBTFNFfQpicmVzY29sbF91aGxtYW5uX3NpbXMgPC0gCiAgYnJlc2NvbGxfdWhsbWFubiAlPiUgCiAgcmVkZXNpZ24oTiA9IHNlcSgxMDAsIDIwMDAsIGJ5ID0gMTAwKSkgJT4lICMgZnJvbSAxMDAgdG8gMjAwMCBpbiBpbmNyZW1lbnRzIG9mIDEwMAogIGRpYWdub3NlX2Rlc2lnbigpICU+JSAKICBnZXRfZGlhZ25vc2FuZHMoKQpgYGAKCgpBbHRlcm5hdGl2ZWx5IHdlIGNhbiB1c2UgdGhlIG1vcmUgZWxlZ2FudCBhbmQgZ2VuZXJhbGl6ZWQgYXBwcm9hY2ggd2l0aCBgZXhwYW5kX2Rlc2lnbmAgYW5kIHRoZSBgdHdvX2J5X3R3b19kZXNpZ25lcmAgYmFzZSBmdW5jdGlvbi4gVGhpcyBpcyByZWNvbW1lbmRlZCBieSB0aGUgZm9sa3MgYXQgRGVjbGFyZURlc2lnbiBiZWNhdXNlIGl0IGlzIGNvbnNpZGVyZWQgc2FmZXIgY29kZSB0byBydW4uCmBgYHtyIGV4cGFuZF9kZXNpZ24sIGV2YWw9RkFMU0V9CmJyZXNjb2xsX3VobG1hbm5fc2ltcyA8LSAKICB0d29fYnlfdHdvX2Rlc2lnbmVyICU+JSAKICBleHBhbmRfZGVzaWduKG1lYW5fQTBCMCA9IDYuNDcsIAogICAgICAgICAgICAgICAgbWVhbl9BMEIxID0gMy43NSwgCiAgICAgICAgICAgICAgICBtZWFuX0ExQjAgPSA0LjA1LCAKICAgICAgICAgICAgICAgIG1lYW5fQTFCMSA9IDUuMDIsCiAgICAgICAgICAgICAgICBvdXRjb21lX3NkcyA9IGxpc3QoYygyLjI1LCAxLjc3LCAxLjYxLCAxLjgpKSwKICAgICAgICAgICAgICAgIE4gPSBzZXEoMTAwLCAyMDAwLCBieSA9IDEwMCkpICAlPiUgIAogIGRpYWdub3NlX2Rlc2lnbigpICU+JSAKICBnZXRfZGlhZ25vc2FuZHMoKQpgYGAKCmBgYHtyLCByZWFkX2RhdGEsIGVjaG89RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CiMgd3JpdGVfY3N2KGJyZXNjb2xsX3VobG1hbm5fc2ltcywgImJyZXNjb2xsX3VobG1hbm5fc2ltcy5jc3YiKQpicmVzY29sbF91aGxtYW5uX3NpbXMgPC0gcmVhZF9jc3YoImJyZXNjb2xsX3VobG1hbm5fc2ltcy5jc3YiKQpgYGAKCkhlcmUgd2UgY2FuIHRha2UgYW4gdWdseSBwZWVrIGF0IG91ciBkYXRhLiBFYWNoIGRpZmZlcmVudCBOIHdlIGV2YWx1YXRlZCB0YWtlcyB1cCAzIGRpZmZlcmVudCByb3dzLiBBIHJvdyBldmFsdWF0aW5nIHRoZSAiQSIsICJCIiwgIkE6QiIgYWthICJlbW90aW9uIiwgImdlbmRlciIgYW5kIHRoZSAiaW50ZXJhY3Rpb24iLiBCZWNhdXNlIHdlIGFyZSByZXBsaWNhdGluZyBhIHByZXZpb3VzIHN0dWR5LCB0aGV5IGRpZCBub3QgY2hhbmdlIHRoZSBtZWFucyBpbiB0aGUgYXNzdW1wdGlvbi4KYGBge3J9CmdsaW1wc2UoYnJlc2NvbGxfdWhsbWFubl9zaW1zKQpgYGAKIyMgVmlzdWFsaXplIG91ciByZXN1bHRzCgpBbmQgbmV4dCB3ZSB2aXN1YWxpemUgdG8gc2VlIHdoZW4gaXQgd2lsbCBiZSB3ZWxsIHBvd2VyZWQhCmBgYHtyLCB3YXJuaW5nPUZBTFNFfQpicmVzY29sbF91aGxtYW5uX3NpbXMgJT4lCiAgbXV0YXRlKGVzdGltYW5kX2xhYmVsID0gY2FzZV93aGVuKGVzdGltYW5kX2xhYmVsID09ICJhdGVfQSIgfiAiRW1vdGlvbiIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGVzdGltYW5kX2xhYmVsID09ICJhdGVfQiIgfiAiR2VuZGVyIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZXN0aW1hbmRfbGFiZWwgPT0gImludGVyYWN0aW9uIiB+ICJFbW90aW9uIHggR2VuZGVyIikpICU+JSAKICBnZ3Bsb3QoYWVzKE4sIHBvd2VyLCAKICAgICAgICAgICAgIGNvbCA9IGVzdGltYW5kX2xhYmVsLCAKICAgICAgICAgICAgIGdyb3VwID0gZXN0aW1hbmRfbGFiZWwsIAogICAgICAgICAgICAgeW1pbiA9IHBvd2VyIC0gMS45Nipgc2UocG93ZXIpYCwgCiAgICAgICAgICAgICB5bWF4ID0gcG93ZXIgKyAxLjk2KmBzZShwb3dlcilgKSkgKyAKICBnZW9tX3Ntb290aChwb3NpdGlvbiA9ICJkb2RnZSIsIHN0YXQgPSAiaWRlbnRpdHkiKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMC44KSArCiAgc2NhbGVfeV9jb250aW51b3VzKGxhYmVscyA9IHNjYWxlczo6cGVyY2VudCkgKwogIHRoZW1lX2J3KCkKYGBgCldlIHNlZSB0aGF0IHRoaXMgZGVzaWduIHdpbGwgbmVlZCBhYm91dCA1MDAgcGFydGljaXBhbnRzIHRvIGJlIHdlbGwtcG93ZXJlZCBvciBhYm91dCAxMjUgaW4gZWFjaCBvZiB0aGUgZm91ciBjZWxscy4gU3VycHJpc2luZ2x5LCB3ZSBzZWUgdGhhdCB0aGUgaW50ZXJhY3Rpb24gaXMgYWN0dWFsbHkgdGhlIGxhcmdlc3QgZWZmZWN0LiBOZXh0IHdlJ3JlIGludGVyZXN0ZWQgaW4gd2hhdCBOIGV4YWN0bHkgZG8gZW1vdGlvbiBhbmQgZ2VuZGVyIHJlYWNoIHRoZSBwb3dlciB0aGUgODAlIHBvd2VyIG1hcmsuCgoKIyMgV2hhdCBpcyB0aGUgbWluaW11bSBudW1iZXIgb2YgcGFydGljaXBhbnRzIG5lZWRlZD8KYGBge3IgbWluX25fdGFibGV9CmJyZXNjb2xsX3VobG1hbm5fc2ltcyAlPiUKICAgIG11dGF0ZShlc3RpbWFuZF9sYWJlbCA9IGNhc2Vfd2hlbihlc3RpbWFuZF9sYWJlbCA9PSAiYXRlX0EiIH4gIkVtb3Rpb24iLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBlc3RpbWFuZF9sYWJlbCA9PSAiYXRlX0IiIH4gIkdlbmRlciIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGVzdGltYW5kX2xhYmVsID09ICJpbnRlcmFjdGlvbiIgfiAiRW1vdGlvbiB4IEdlbmRlciIpKSAlPiUgCiAgZ3JvdXBfYnkoZXN0aW1hbmRfbGFiZWwpICU+JSAKICBmaWx0ZXIocG93ZXIgPiAuOCkgJT4lIAogIGZpbHRlcihOID09IG1pbihOKSkgJT4lIAogIHNlbGVjdChlc3RpbWFuZF9sYWJlbCwgdGVybSwgTiwgcG93ZXIsIGBzZShwb3dlcilgKSAlPiUgCiAga2FibGUoKSAlPiUgCiAga2FibGVfc3R5bGluZyhwb3NpdGlvbiA9ICJjZW50ZXIiKQpgYGAKCgoKIyBTdGVhbGluZyB0aGUgRGVzaWduTGlicmFyeSBjb2RlCgpMYXN0bHkgSSB3YW50ZWQgdG8gbm90ZSBob3cgeW91IGNvdWxkIGN1c3RvbWl6ZSBEZXNpZ25MaWJyYXJ5IGNvZGUgZm9yIHlvdXIgb3duIGRlc2lnbnMsIGxldCdzIHNheSB5b3UgZm91bmQgYSBjb2RlIGluIERlc2lnbkxpYnJhcnkgdGhhdCBpcyBhbG1vc3QgcmlnaHQgZm9yIHlvdXIgbmVlZHMsIGJ1dCBub3QgcXVpdGUuIFRvIHNlZSB0aGUgaW50ZXJuYWxzIGluIGEgRGVzaWduTGlicmFyeSBmdW5jdGlvbiB1c2UgYGdldF9kZXNpZ25fY29kZWAuIAoKSW4gb3JkZXIgdG8gdGFrZSB0aGUgY29kZSBmb3IgYSBjb3B5IHBhc3RpbmcgeW91LCBjYW4gYGNhdChwYXN0ZShjb2RlLCBjb2xsYXBzZSA9ICJcbiIpKWAuIFRoaXMgaXMgZXF1aXZhbGVudCB0byBgcHJpbnRfdGFibGVgIGluIHRoZSBgYmxwbGAgcGFja2FnZS4gTm90ZSB0aGF0IHRoZSAjIydzIHdvbid0IHNob3cgdXAgaW4geW91ciBSIGNvbnNvbGUuCmBgYHtyfQp0d294dHdvX2NvZGUgPC0gCiAgZ2V0X2Rlc2lnbl9jb2RlKHR3b19ieV90d29fZGVzaWduZXIoKSkgIyBnZW5lcmFsaXplZCAyeDIsIG5vdGUgdGhlIHBhcmVudGhlc2lzIG9uIHRoZSBpbnNpZGUgZnVuY3Rpb24KCgpidV9jb2RlIDwtIGdldF9kZXNpZ25fY29kZShicmVzY29sbF91aGxtYW5uKSAjIHRoZSBleGFtcGxlIHdlIGp1c3QgbWFkZQoKY2F0KHBhc3RlKGJ1X2NvZGUsIGNvbGxhcHNlID0gIlxuIikpCmBgYAoKCg==