Prepare analysis workflow

Set parameters

knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file(),
                     fig.width=15,
                     digit=5,
                     scipen=8)
options(readr.show_progress = FALSE,
        digits=5, 
        scipen=8,
        future.globals.maxSize = +Inf)

Set filepaths and parameters

project_dir <- rprojroot::find_rstudio_root_file()
if(is.null(project_dir)){
  project_dir <- getwd()
  warning(sprintf("No rstudio project root file  found. 
                  Setting project directory to current workflow.Rmd file location: %s. 
                  Override if needed.",
                  project_dir))
 
}
message(sprintf("Project directory: %s",
                project_dir))
Project directory: /home/rfarouni/Documents/index_hopping

Load libraries

library(rhdf5)
#library(DropletUtils) # install but not load
library(tidyverse)
library(matrixStats)
library(broom)
library(furrr)
library(tictoc)
library(data.table)
library(cowplot)
plan(multiprocess)

Load functions

code_dir <- file.path(project_dir, "code")
source(file.path(code_dir, "1_create_joined_counts_table.R"))
source(file.path(code_dir, "2_create_counts_by_outcome_table.R"))
source(file.path(code_dir, "3_estimate_sample_index_hopping_rate.R"))
source(file.path(code_dir, "4_compute_summary_statistics.R"))
source(file.path(code_dir, "5_reassign_hopped_reads.R"))
source(file.path(code_dir, "6_purge_phantom_molecules.R"))
source(file.path(code_dir, "7_call_cells.R"))
source(file.path(code_dir, "8_summarize_purge.R"))
source(file.path(code_dir, "9_plotting_functions.R"))

Load data

validation_output_dir <- file.path(project_dir, "data", "hiseq4000_validation")
data <- read_tsv(file.path(validation_output_dir,
                           "hiseq4000_inner_joined_with_labels.txt"))
Parsed with column specification:
cols(
  cell = col_character(),
  gene = col_double(),
  umi = col_double(),
  s1_nonplexed = col_double(),
  s2_nonplexed = col_double(),
  s1_plexed = col_double(),
  s2_plexed = col_double(),
  outcome = col_character(),
  label = col_character()
)
data
ABCDEFGHIJ0123456789
cell
<chr>
gene
<dbl>
umi
<dbl>
s1_nonplexed
<dbl>
s2_nonplexed
<dbl>
s1_plexed
<dbl>
s2_plexed
<dbl>
outcome
<chr>
label
<chr>
AAACCTGAGAAACCTA1412115900010101,0r,0
AAACCTGAGAAACCTA810884239230202,0r,0
AAACCTGAGAAACCTA260513854610101,0r,0
AAACCTGAGAAACCTA272585581250202,0r,0
AAACCTGAGAAACCTA221157712420202,0r,0
AAACCTGAGAAACCTA533925023310303,0r,0
AAACCTGAGAAACCTA1961930254430303,0r,0
AAACCTGAGAAACCTA2088464417930202,0r,0
AAACCTGAGAAACCTA2614104337820202,0r,0
AAACCTGAGAAACCTA417075174932401851185,1r,f
Next
123456
...
100
Previous
1-10 of 9,252,095 rows

Compute Index hopping rate

Estimates conditional on duplication level and label

summary_counts <-
    data %>%
  mutate(r=as.integer(rowSums(.[c(6,7)]))) %>%
  arrange(r) %>%
  group_by(r, label) %>%
  summarize_at(vars(matches("^s1_pl|s2_pl")), 
               list(~ sum(.)))  %>%
  mutate(s1_hopped=if_else(label %in% c("0,f", "r,f"), s2_plexed,0),
         s2_hopped=if_else(label %in% c("f,0", "f,r"), s1_plexed,0),
         s1_nonhopped=if_else(label %in% c("r,0", "r,f"), s1_plexed,0),
         s2_nonhopped=if_else(label %in% c("0,r", "f,r"), s2_plexed,0)) %>%
  select(-s1_plexed,- s2_plexed)
summary_counts
ABCDEFGHIJ0123456789
r
<int>
label
<chr>
s1_hopped
<dbl>
s2_hopped
<dbl>
s1_nonhopped
<dbl>
s2_nonhopped
<dbl>
10,f3266000
10,r0001863385
1f,00644600
1r,0006996670
20,f330000
20,r0001542816
2f,0012600
2f,r0503805038
2r,0006409760
2r,f2436024360
Next
123456
...
97
Previous
1-10 of 965 rows

Estimates conditional on duplication level

summary_counts_conditional <-
summary_counts  %>%
  group_by(r) %>%
  summarize_at(vars(matches("^s")), 
               list(~ sum(.)))%>%
  mutate(SIHR_12 = s1_hopped/(s1_hopped+s1_nonhopped),
         SIHR_21 = s2_hopped/(s2_hopped+s2_nonhopped),
         frac_s1= (s1_hopped+s1_nonhopped)/ (s1_hopped+s1_nonhopped+s2_hopped+s2_nonhopped))
summary_counts_conditional
ABCDEFGHIJ0123456789
r
<int>
s1_hopped
<dbl>
s2_hopped
<dbl>
s1_nonhopped
<dbl>
s2_nonhopped
<dbl>
SIHR_12
<dbl>
SIHR_21
<dbl>
frac_s1
<dbl>
13266644669966718633850.004646250.00344740.273221
22766516464341215478540.004280550.00332510.293825
31678332739970110231380.004180590.00324120.281108
4105126282382057966720.004392780.00328790.230374
574424631693117698320.004375060.00318920.180458
655126331457058256490.003767370.00317890.150077
755028141408019084530.003891020.00308800.134285
8557308314194710058450.003908660.00305570.123762
9559331914644711161280.003802570.00296490.116077
10501378914757912335010.003383310.00306230.106888
Next
123456
...
42
Previous
1-10 of 416 rows

Marginal estimates

summary_counts_marginal <- 
  summary_counts %>%
  ungroup() %>%
  summarize_at(vars(matches("^s")), 
               list(~ sum(.)))%>%
  mutate(SIHR_12 = s1_hopped/(s1_hopped+s1_nonhopped),
         SIHR_21 = s2_hopped/(s2_hopped+s2_nonhopped),
         SIHR=1-(s1_nonhopped+s2_nonhopped)/(s1_hopped+s1_nonhopped +s2_hopped+s2_nonhopped),
         frac_s1= (s1_hopped+s1_nonhopped)/ (s1_hopped+s1_nonhopped+s2_hopped+s2_nonhopped))
summary_counts_marginal
ABCDEFGHIJ0123456789
s1_hopped
<dbl>
s2_hopped
<dbl>
s1_nonhopped
<dbl>
s2_nonhopped
<dbl>
SIHR_12
<dbl>
SIHR_21
<dbl>
SIHR
<dbl>
frac_s1
<dbl>
31790423021091653829759760170.00345650.00302090.00325910.54687
1 row
p1 <- 
  ggplot(summary_counts_conditional) +
    geom_line(aes(x = r,
                    y = SIHR_12*100,
               colour="SIHR_12"))+
    geom_line(aes(x = r,
                    y = SIHR_21*100,
               colour="SIHR_21")) +
      geom_line(aes(x = r,
                    y = frac_s1,
               colour="frac_s1")) +
    geom_hline(yintercept = unlist(summary_counts_marginal[5:7])* c(100, 100, 100),
               linetype="dashed") +
    xlim(0,90) + 
  ylim(0,1) 
p1

NA

Molecules

summary_mol_counts <-
  data %>%
  mutate(r=as.integer(rowSums(.[c(6,7)]))) %>%
  arrange(r) %>%  
  mutate_at(vars(matches("^s")), 
            list(~ as.integer(.!=0)))  %>%
  group_by(r, label) %>%
  summarize_at(vars(matches("^s1_pl|^s2_pl")), 
               list(~ sum(.)))  %>%
    mutate(s1_phantom=if_else(label %in% c("0,f", "r,f"), s2_plexed,0L),
         s2_phantom=if_else(label %in% c("f,0", "f,r"), s1_plexed,0L),
         s1_real=if_else(label %in% c("r,0", "r,f"), s1_plexed,0L),
         s2_real=if_else(label %in% c("0,r", "f,r"), s2_plexed,0L)) %>%
  select(-s1_plexed,- s2_plexed) 
summary_mol_counts
ABCDEFGHIJ0123456789
r
<int>
label
<chr>
s1_phantom
<int>
s2_phantom
<int>
s1_real
<int>
s2_real
<int>
10,f3266000
10,r0001863385
1f,00644600
1r,0006996670
20,f165000
20,r000771408
2f,006300
2f,r0503805038
2r,0003204880
2r,f2436024360
Next
123456
...
97
Previous
1-10 of 965 rows
summary_mol_counts_marginal <-
  summary_mol_counts%>%
  ungroup() %>%
  summarize_at(vars(matches("^s")), 
               list(~ sum(.))) %>%
  mutate(ppm_12 = s1_phantom/(s1_phantom+s1_real),# prop hopped phantom molec
         ppm_21 = s2_phantom/(s2_phantom+s2_real),
         ppm_1 = s2_phantom/(s2_phantom+s1_real), # prop phantom molec
         ppm_2 = s1_phantom/(s1_phantom+s2_real),
         ppm=(s1_phantom+s2_phantom)/(s1_phantom+s1_real+s2_phantom+s2_real),
         frac_mol_s1= (s1_phantom+s1_real)/ (s1_phantom+s1_real+s2_phantom+s2_real))
summary_mol_counts_marginal
ABCDEFGHIJ0123456789
s1_phantom
<int>
s2_phantom
<int>
s1_real
<int>
s2_real
<int>
ppm_12
<dbl>
ppm_21
<dbl>
ppm_1
<dbl>
ppm_2
<dbl>
ppm
<dbl>
263741221658238933768527630.0994090.0313320.0848940.037060.0499
1 row | 1-9 of 10 columns
summary_mol_counts_conditional <-
  summary_mol_counts%>%
  group_by(r) %>%
  summarize_at(vars(matches("^s")), 
               list(~ sum(.))) %>%
  mutate(ppm_12 = s1_phantom/(s1_phantom+s1_real),
         ppm_21 = s2_phantom/(s2_phantom+s2_real),
         ppm=(s1_phantom+s2_phantom)/(s1_phantom+s1_real+s2_phantom+s2_real),
         frac_mol_s1= (s1_phantom+s1_real)/ (s1_phantom+s1_real+s2_phantom+s2_real))
summary_mol_counts_conditional
ABCDEFGHIJ0123456789
r
<int>
s1_phantom
<int>
s2_phantom
<int>
s1_real
<int>
s2_real
<int>
ppm_12
<dbl>
ppm_21
<dbl>
ppm
<dbl>
frac_mol_s1
<dbl>
13266644669966718633850.00464620.00344740.00377490.273221
2260151013229247764460.00799020.00652680.00695710.294041
3154732841337533421530.01143390.00950680.01004920.281443
49522593598021998250.01566980.01281010.01347030.230853
56932439340101544590.01996950.01554510.01634650.181121
65292598243761380470.02124070.01847200.01888860.150438
75322775201931301810.02566950.02087160.02151860.134857
85413028178131261160.02947590.02344670.02419690.124436
95483267163341243830.03246060.02559340.02639550.116805
104933731148081237290.03222010.02927190.02958790.107179
Next
123456
...
42
Previous
1-10 of 416 rows
p4 <- ggplot(summary_mol_counts_conditional) +
  geom_line(aes(x = r,
                  y = ppm_12,
             colour="ppm_12"))+
  geom_line(aes(x = r,
                  y =  ppm_21,
             colour=" ppm_21")) +
  #geom_hline(yintercept = unlist(summary_mol_counts_marginal[5:7]), linetype="dashed") +
  xlim(0,300) 
p4

 #ggsave("phantom_molecules_validation.pdf", p4,  width =8, height = 5)

Examine extent of contamination in cells

summary_mol_counts_cell<-
  data %>%
  filter(label!="NA") %>%
  mutate_at(vars(matches("^s")), 
            list(~ as.integer(.!=0)))  %>%
  group_by(cell, label) %>%
  summarize_at(vars(matches("^s1_pl|^s2_pl")), 
               list(~ sum(.)))  %>%
    mutate(s1_phantom=if_else(label %in% c("f,0", "f,r"), s1_plexed,0L),
         s2_phantom=if_else(label %in% c("0,f", "r,f"), s2_plexed,0L),
         s1_real=if_else(label %in% c("r,0", "r,f"), s1_plexed,0L),
         s2_real=if_else(label %in% c("0,r", "f,r"), s2_plexed,0L)) %>%
  select(-s1_plexed,- s2_plexed) %>%
  group_by(cell) %>%
  summarize_at(vars(matches("^s")), 
               list(~ sum(.))) %>% 
    mutate_at(vars(matches("^s")), 
            list(nonempty= ~ as.integer(.!=0)))  %>%
  unite(label,matches("nonempty"), sep=",") %>%
  mutate(cell_status= 
           case_when(label %in% c("0,0,0,1","0,0,1,0") ~ "real",
                     label %in% c("0,0,1,1") ~ "real",
                     label %in% c("1,0,0,0","0,1,0,0", "1,1,0,0") ~  "phantom",
                     label %in% c("1,0,0,1","0,1,1,0") ~  "phantom",
                     TRUE ~  "contaminated"))  %>%
  mutate(s1_total =(s1_phantom+s1_real),
         s2_total =(s2_phantom+s2_real),
         s1_ppm = s1_phantom/(s1_total),
         s2_ppm = s2_phantom/(s2_total))
  
summary_mol_counts_cell
ABCDEFGHIJ0123456789
cell
<chr>
s1_phantom
<int>
s2_phantom
<int>
s1_real
<int>
s2_real
<int>
label
<chr>
cell_status
<chr>
s1_total
<int>
AAACCTGAGAAACCTA011000,1,1,0phantom10
AAACCTGAGAAACGCC00010,0,0,1real0
AAACCTGAGAAAGTGG00900,0,1,0real9
AAACCTGAGAACAACT000240,0,0,1real0
AAACCTGAGAACTCGG000260,0,0,1real0
AAACCTGAGAACTGTA00010,0,0,1real0
AAACCTGAGAAGCCCA000560,0,0,1real0
AAACCTGAGAAGGCCT100351,0,0,1phantom1
AAACCTGAGAAGGTGA0012220,0,1,1real12
AAACCTGAGAATGTGT000120,0,0,1real0
Next
123456
...
100
Previous
1-10 of 244,201 rows | 1-8 of 11 columns
cell_status_tally <-
  summary_mol_counts_cell %>%
  group_by(cell_status,label) %>%
  tally(sort=TRUE)
cell_status_tally
ABCDEFGHIJ0123456789
cell_status
<chr>
label
<chr>
n
<int>
real0,0,0,198256
real0,0,1,067769
phantom1,0,0,126246
phantom0,1,1,025291
real0,0,1,115105
contaminated0,1,1,15561
contaminated1,0,1,14407
contaminated1,1,1,11502
phantom1,0,0,030
phantom0,1,0,015
Next
12
Previous
1-10 of 14 rows
n_affected_cells <- 
  cell_status_tally %>%
  ungroup() %>%
  filter(cell_status!="real") %>% 
  summarise(n=sum(n)) %>% 
  pull(n)
n_affected_cells <- n_affected_cells + 1502 +6
n_affected_cells
[1] 64579
n_total_cells <-
  summary_mol_counts_cell %>%
  mutate_at(vars(matches("_total")), 
            list(~ as.integer(.!=0))) %>%
  ungroup() %>%
  summarise_at(vars(matches("_total")), list(~sum(.))) %>%
  mutate(cells_total= s1_total+s2_total)%>%
  pull(cells_total)
n_total_cells 
[1] 322321

Proportion of affected cells

p_affected_cells <-  n_affected_cells/n_total_cells
p_affected_cells
[1] 0.20036

For each cell-barcode, plot the number of phantom molecules against the number of total molecules associated with it.

Sample 1 plot

p2 <- ggplot(summary_mol_counts_cell %>%
               filter(cell_status!="real" & s1_phantom >0 )) +
  geom_point(aes(x = s1_total,
                  y = s1_phantom,
                 colour=cell_status)) +
  scale_x_log10() +
  scale_y_log10()
p2

 #ggsave("phantom_molecules_validation.pdf", p2,  width =8, height = 5)

Sample 2 plot

p3 <- ggplot(summary_mol_counts_cell %>%
               filter(cell_status!="real" & s2_phantom >0 )) +
  geom_point(aes(x = s2_total,
                  y = s2_phantom,
                 colour=cell_status)) +
  scale_x_log10() +
  scale_y_log10()
p3

#ggsave("phantom_cells_validation_s2.pdf", p3,  width =8, height = 5)

Run workflow on multiplexed data

read_counts <- 
  data %>% 
  filter(label!="NA")%>%
  select(-ends_with("nonplexed")) %>%
  set_names(c("cell", "gene", "umi", "s1", "s2", "outcome", "label"))
read_counts
ABCDEFGHIJ0123456789
cell
<chr>
gene
<dbl>
umi
<dbl>
s1
<dbl>
s2
<dbl>
outcome
<chr>
label
<chr>
AAACCTGAGAAACCTA14121159000101,0r,0
AAACCTGAGAAACCTA8108842392202,0r,0
AAACCTGAGAAACCTA2605138546101,0r,0
AAACCTGAGAAACCTA2725855812202,0r,0
AAACCTGAGAAACCTA2211577124202,0r,0
AAACCTGAGAAACCTA5339250233303,0r,0
AAACCTGAGAAACCTA19619302544303,0r,0
AAACCTGAGAAACCTA20884644179202,0r,0
AAACCTGAGAAACCTA26141043378202,0r,0
AAACCTGAGAAACCTA41707517491851185,1r,f
Next
123456
...
100
Previous
1-10 of 9,252,095 rows
S <- 2
sample_names <- colnames(read_counts)[4:(S+3)]
sample_names
[1] "s1" "s2"
tic("Step 2: creating outcome counts datatable with grouping vars")
outcome_counts <- create_outcome_counts(read_counts%>%
                                          select(-label), 
                                        sample_names,  
                                        min_frac=0.8)
toc()
Step 2: creating outcome counts datatable with grouping vars: 1.591 sec elapsed
outcome_counts
ABCDEFGHIJ0123456789
outcome
<chr>
n
<int>
r
<int>
k_chimera
<int>
s1
<dbl>
s2
<dbl>
s_maxprop
<dbl>
0,1186665111012
1,070611311101
0,277157321022
2,032055121201
1,1747422110
0,333891131032
3,013224831301
1,2329432120
2,1149532210
0,419724441042
Next
123456
...
100
Previous
1-10 of 1,974 rows
tic("Step 3: creating a chimera counts datatable and estimating hopping rate")
  fit_out <-
    estimate_hopping_rate(
      outcome_counts,
      S
    )
  toc()
Step 3: creating a chimera counts datatable and estimating hopping rate: 0.083 sec elapsed
fit_out 
$glm_estimates
ABCDEFGHIJ0123456789
max_r
<int>
phat
<dbl>
phat_low
<dbl>
phat_high
<dbl>
SIHR
<dbl>
SBIHR
<dbl>
5080.99680.996780.996810.00320420.0056074
1 row

$chimera_counts
ABCDEFGHIJ0123456789
r
<dbl>
1
<dbl>
2
<dbl>
n_cugs
<dbl>
n_chimeras
<dbl>
p_chimeras
<dbl>
pcum_cugs
<dbl>
phat_chimeras
<dbl>
phat_chimeras_low
<dbl>
125727640257276400.0000000.278070.00320420.0032162
210921247474109959874740.0067970.396920.00639820.0064221
3471159478947594847890.0100620.448360.00958190.0096177
4256106353325963935330.0136070.476430.01275540.0128030
5185339313118847031310.0166130.496800.01591870.0159780
6159296312716242331270.0192520.514350.01907190.0191428
7147067330715037433070.0219920.530610.02221510.0222975
8140360356914392935690.0247970.546160.02534810.0254420
9136902381514071738150.0271110.561370.02847110.0285764
10134313422413853742240.0304900.576350.03158410.0317007
Next
123456
...
42
Previous
1-10 of 416 rows | 1-9 of 10 columns
NA
  # compute_molecular_complexity_profile
  tic("Step 4: compute molecular complexity profile and other summary statistics")
  summary_stats <-
    compute_summary_stats(
      outcome_counts,
      fit_out$glm_estimates$phat,
      sample_names
    )
  toc()
Step 4: compute molecular complexity profile and other summary statistics: 0.16 sec elapsed
summary_stats
$summary_estimates
ABCDEFGHIJ0123456789
n_reads
<int>
n_cugs
<dbl>
n_pm
<dbl>
u
<dbl>
g
<dbl>
frm
<dbl>
p_chimeras
<dbl>
16817796092520954836590.0513830.0008922318.1770.051383
1 row

$marginal
ABCDEFGHIJ0123456789
sample
<chr>
n_reads
<dbl>
prop_reads
<dbl>
n_molecs
<int>
p_molecs
<dbl>
FRM
<dbl>
s1918840390.5463526109950.2684135.191
s2762939210.4536571165040.7315910.721
2 rows

$conditional
ABCDEFGHIJ0123456789
r
<int>
n_cugs
<dbl>
p_cugs
<dbl>
s1
<dbl>
s2
<dbl>
max_hop_rate
<dbl>
125727640.278073668720.2744570.725543030.27445697
210995980.118848541870.2949150.705084950.29491505
34759480.051442186880.2822630.717736670.28226333
42596390.028062725250.2318920.768107830.23189217
51884700.020370521490.1822830.817717410.18228259
61624230.017555267210.1522140.847786340.15221366
71503740.016252967570.1364360.863563990.13643601
81439290.015556368580.1259560.874043800.12595620
91407170.015209203970.1182560.881743740.11825626
101385370.014973581660.1092620.890738210.10926179
Next
123456
...
42
Previous
1-10 of 416 rows

$pi_r_hat
ABCDEFGHIJ0123456789
r
<int>
s1
<dbl>
s2
<dbl>
10.2730020.72699773429
20.2935920.70640769902
30.2808590.71914102194
40.2301630.76983706695
50.1802330.81976661133
60.1499710.85002947592
70.1340910.86590889238
80.1235440.87645629409
90.1157940.88420589375
100.1067420.89325838412
Next
123456
...
42
Previous
1-10 of 416 rows
NA

Set the trade-off ratio cost cutoff (torc). The parameter torc represents the number of real molecules one is willing to incorrectly discard in order to correctly purge one phantom molecule. Since discarding a large proportion of the data is undesirable, reasonable values of torc are expected to be within the range of 1-5.

torc <- 3 
tic("Step 5: reassign read counts, determine cutoff, and mark retained observations")
  outcome_counts <-
    reassign_reads_and_mark_retained_observations(
      outcome_counts,
      summary_stats,
      sample_names,
      fit_out,
      torc
    )
  # get the tradoff ratio cutoff
  summary_stats <- get_threshold(outcome_counts, summary_stats)
  toc()
Step 5: reassign read counts, determine cutoff, and mark retained observations: 1.176 sec elapsed
tic("Step 6: Purge and save read counts datatable to disk")
read_counts <-
  left_join(read_counts %>%
    select(outcome, cell, umi, gene, sample_names, label),
  outcome_counts,
  by = c("outcome")
  ) %>%
  select(-outcome)
toc()
Step 6: Purge and save read counts datatable to disk: 2.166 sec elapsed

Compare the SIHR estimates with ground truth estimates

p5 <- 
  ggplot(summary_counts_conditional) +
  geom_line(aes(x = r,
                  y = SIHR_12,
             colour="12"))+
  geom_line(aes(x = r,
                  y = SIHR_21,
             colour="21"))  +
geom_hline(aes(yintercept =  summary_counts_marginal$SIHR, 
               colour="true mean"), 
           linetype="solid",
           size=.5)   +
geom_hline(aes(yintercept =  fit_out$glm_estimates$SIHR, 
               colour="estimate"),
           linetype="solid",
           size=.1)   +
    geom_linerange(data=summary_counts_conditional, 
                  aes(x=r,
                      ymax=1-fit_out$glm_estimates$phat_low,
                      ymin=1-fit_out$glm_estimates$phat_high,
                      colour="estimate"), 
                  size=.5)+ 
  xlim(1,210) +
  ylim(0.002,0.005)
#ggsave("index_hopping_rate_200.pdf", p5, width=9, height=6)
p5

Determine the number of false positives and false negatives

read_counts <-
    read_counts %>%
  ungroup() %>%
  arrange(-qr) %>%
  mutate(t= case_when(
      label %in% c("f,r","0,r") ~ 2,
      label %in% c("r,f","r,0") ~ 1    ),
    f= case_when(
      label %in% c("f,r","f,0") ~ 1,
      label %in% c("r,f","0,f") ~ 2    )) %>%
  mutate(tp= if_else( t == s, 1L, 0L, missing =0L),
         fp= if_else( f == s, 1L, 0L, missing =0L),
         tn= if_else( f != s, 1L, 0L, missing =0L),
         fn= if_else( t != s, 1L, 0L, missing =0L),
         tp0= if_else( t == 0, 1L, 0L, missing =0L), #0 if  predict all molecules to be phantom
         fn0= if_else( t != 0, 1L, 0L, missing =0L),
         tn0= if_else( f != 0, 1L, 0L, missing =0L),
         fp0= if_else( f == 0, 1L, 0L, missing =0L),
         tp_max= if_else( t == s_maxprop, 1L, 0L, missing =0L),
         fp_max= if_else( f == s_maxprop, 1L, 0L, missing =0L),
         tn_max= if_else( f != s_maxprop, 1L, 0L, missing =0L),
         fn_max= if_else( t != s_maxprop, 1L, 0L, missing =0L))

The maximum read fraction method

false_counts_maxprop <-
  read_counts %>%
  summarize_at(vars(c("tp_max", "fp_max", "tn_max", "fn_max")),
            list( ~ sum(.))) %>%
  set_names(c("tp", "fp", "tn", "fn"))
false_counts_maxprop
ABCDEFGHIJ0123456789
tp
<int>
fp
<int>
tn
<int>
fn
<int>
92260061000047539916094
1 row

The tor method

false_counts_min_cutoff <-
  read_counts %>%
  summarize_at(vars(c("tp", "fp", "tn", "fn")),
            list( ~ sum(.))) 
false_counts_min_cutoff
ABCDEFGHIJ0123456789
tp
<int>
fp
<int>
tn
<int>
fn
<int>
9239510125854728142590
1 row

No purging

false_counts_nopurging <- 
  read_counts%>%
  summarize(n_cugs = n(),
            n_real= sum(t>0, na.rm = TRUE),
            n_fantom = sum(f>0, na.rm = TRUE),
            n_mol=n_real+n_fantom,
            g =n_cugs- n_real,
            u = n_mol-n_cugs,
            tp=n_mol-g,
            fp=u+g,
            tn=0,
            fn=0)
false_counts_nopurging
ABCDEFGHIJ0123456789
n_cugs
<int>
n_real
<int>
n_fantom
<int>
n_mol
<int>
g
<int>
u
<int>
tp
<int>
fp
<int>
tn
<dbl>
fn
<dbl>
9252095924210048539997274999995475404971750448539900
1 row

TOR cutoff

read_counts <-
  read_counts %>%
  mutate_at(vars(c("tp", "fp", "tn", "fn","tp0", "fp0", "tn0", "fn0")),
            list(cum= ~ cumsum(.)))  %>%
  mutate_at(vars(c("tp_cum", "fp_cum", "tn_cum", "fn_cum")), 
            list( ~ (last(.)-lag(., default =0)))) %>%
  mutate(tp_t=tp_cum + tp0_cum,
         fp_t=fp_cum + fp0_cum,
         tn_t=tn_cum + tn0_cum,
         fn_t=fn_cum + fn0_cum,
         fpm= first(fp_t)- fp_t, 
         fnm= fn_t-first(fn_t),
         tor_true= fnm/fpm)
false_counts_tor_cutoff <-
  read_counts%>%
  filter(retain) %>%
  slice(1)%>%
  select(c("s1", "s2", "qr", "tor", "tp_t", "fp_t", "tn_t", "fn_t", "fpm", "fnm", "tor_true"))  
false_counts_tor_cutoff
ABCDEFGHIJ0123456789
s1
<dbl>
s2
<dbl>
qr
<dbl>
tor
<dbl>
tp_t
<dbl>
fp_t
<dbl>
tn_t
<dbl>
fn_t
<dbl>
fpm
<dbl>
fnm
<dbl>
100.00848752.4222923441610116475283768524695094
1 row | 1-10 of 11 columns

Create comparison datatable

false_counts_dt <-
  bind_rows(
    list(no_purging=false_counts_nopurging %>%
           select(c("tp", "fp", "tn", "fn")),
         no_cutoff=false_counts_min_cutoff,
         tor_cutoff=
           false_counts_tor_cutoff %>%
           select(c("tp_t", "fp_t", "tn_t", "fn_t")) %>%
           set_names(c("tp", "fp", "tn", "fn")),
         max_frac=false_counts_maxprop),
    .id="approach") %>%
  select(approach, fp,fn, tp, tn) %>%
  mutate(fpr=fp/false_counts_nopurging$n_fantom,
         fnr=fn/false_counts_nopurging$n_real)
false_counts_dt
ABCDEFGHIJ0123456789
approach
<chr>
fp
<dbl>
fn
<dbl>
tp
<dbl>
tn
<dbl>
fpr
<dbl>
fnr
<dbl>
no_purging4853990971750401.0000000.00000000
no_cutoff12585259092395104728140.0259270.00028024
tor_cutoff10116768592344164752830.0208410.00083152
max_frac100001609492260064753990.0206020.00174138
4 rows

Plots

Datatable for plotting

classification_curves <-
  read_counts %>% 
  group_by(qr) %>%
  slice(1L) %>%
  ungroup() %>%
  select( qr, qs, tor, retain, fp_t, FP, fn_t, FN, tp_t, tn_t, TP, TN, FPm, FNm, fpm, fnm, tor_true,o,r) %>%
  mutate(fpr=fp_t/false_counts_nopurging$n_fantom,
         fnr=fn_t/false_counts_nopurging$n_real)
classification_curves
ABCDEFGHIJ0123456789
qr
<dbl>
qs
<dbl>
tor
<dbl>
retain
<lgl>
fp_t
<dbl>
FP
<dbl>
fn_t
<dbl>
FN
<dbl>
tp_t
<dbl>
tn_t
<dbl>
0.0000e+000.0000455.6849TRUE00475585247782594486249485400
2.2204e-165.0792440.7896TRUE00461737146221424624730485399
6.2172e-1518.0053438.2251TRUE00459352445952634648577485399
7.7716e-1518.9606438.1748TRUE00459299745947364649104485400
9.3259e-1519.7432438.1740TRUE00459298945947284649112485400
3.7970e-1425.8058438.1738TRUE00459298745947264649114485400
4.4853e-1426.5276438.1691TRUE00459293745946764649164485400
5.3069e-1427.2566437.9079TRUE10459020145919394651900485399
7.5495e-1428.7849423.4035TRUE20443818044399184803920485398
1.5612e-1241.9349420.2246TRUE20440486344066004837238485397
Next
12345
Previous
1-10 of 45 rows | 1-10 of 21 columns

Preformance Plots

p_tradeoff <-  
    ggplot(classification_curves) + 
    geom_point(
      aes(x = FPm,
          y = FNm),
      size=.5)+
    geom_line(
      aes(x = FPm,
          y = FPm,
          colour="1")
    ) +
    geom_line(
      aes(x = FPm,
          y = 2*FPm,
          colour="2"))+
    geom_line(
      aes(x = FPm,
          y = 3*FPm,
          colour="3"))+
    geom_line(
      aes(x = FPm,
          y = 4*FPm,
          colour="4"))+
    geom_line(
      aes(x = FPm,
          y = 5*FPm,
          colour="5"))+
    geom_line(
      aes(x = FPm,
          y = 9*FPm,
          colour="9"))+
    scale_y_log10() +
    theme_bw()  +
      theme(
        legend.title = element_text(face = "bold")) + 
      scale_color_discrete(name = "TORC") +
    labs(x="Marginal Decrease in False Positives (reduce phantom molecs) ",
         y="Marginal Increase in False Negatives (discard real molecs)") 
    
#ggsave(file.path(figures_dir, "validation_tradeoff.pdf"), p_tradeoff, width=9, height=6)
p_tradeoff

p6 <-
  ggplot(classification_curves)  + 
    geom_point(
             aes(x = fp_t,
                  y = fn_t,
                 colour="true")) +
  geom_line(
             aes(x = fp_t,
                  y = fn_t,
                 colour="true")) +
    geom_line(
             aes(x = FP,
                  y = FN,
                 colour="predicted"))+
      geom_point(
             aes(x = FP,
                  y = FN,
                 colour="predicted"))+ 
  geom_point(data=false_counts_dt ,
             aes(x = fp,
                  y = fn,
                  shape=approach),
             size=2) +
    labs(x="False Positive Count",
       y="False Negative Count")  + 
    scale_y_sqrt() + 
  scale_x_sqrt() 
#ggsave("peformance_groundtruth.pdf", p6, width=9, height=6)
 p6

p7 <- ggplot(false_counts_dt
              %>% filter(approach %in% c("no_cutoff", "tor_cutoff", "max_frac")))  + 
  geom_point(
             aes(x = fp ,
                  y = fn,
                  color=approach),
             size=2) 
#ggsave("peformance_zoom.pdf", p7, width=9, height=6)
p7

sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] cowplot_0.9.4      data.table_1.12.2  tictoc_1.0         furrr_0.1.0        future_1.13.0      broom_0.5.2        matrixStats_0.54.0 forcats_0.4.0     
 [9] stringr_1.4.0      dplyr_0.8.1        purrr_0.3.2        readr_1.3.1        tidyr_0.8.3        tibble_2.1.2       ggplot2_3.1.1      tidyverse_1.2.1   
[17] rhdf5_2.28.0      

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5 xfun_0.7         listenv_0.7.0    haven_2.1.0      lattice_0.20-38  colorspace_1.4-1 generics_0.0.2   yaml_2.2.0      
 [9] rlang_0.3.4      pillar_1.4.1     glue_1.3.1       withr_2.1.2      modelr_0.1.4     readxl_1.3.1     plyr_1.8.4       munsell_0.5.0   
[17] gtable_0.3.0     cellranger_1.1.0 rvest_0.3.4      codetools_0.2-16 labeling_0.3     knitr_1.23       parallel_3.6.0   Rcpp_1.0.1      
[25] scales_1.0.0     backports_1.1.4  jsonlite_1.6     hms_0.4.2        digest_0.6.19    stringi_1.4.3    grid_3.6.0       rprojroot_1.3-2 
[33] cli_1.1.0        tools_3.6.0      magrittr_1.5     lazyeval_0.2.2   crayon_1.3.4     pkgconfig_2.0.2  MASS_7.3-51.1    xml2_1.2.0      
[41] lubridate_1.7.4  assertthat_0.2.1 httr_1.4.0       rstudioapi_0.10  Rhdf5lib_1.6.0   R6_2.4.0         globals_0.12.4   nlme_3.1-140    
[49] compiler_3.6.0
LS0tCnRpdGxlOiAiUGhhbnRvbSBQdXJnZSIKc3VidGl0bGU6ICJWYWxpZGF0aW9uIEFuYWx5c2lzOiBQYXJ0IElJIgphdXRob3I6IAotIG5hbWU6IFJpY2sgRmFyb3VuaQogIGFmZmlsaWF0aW9uOgogIC0gJmNydWsgR8Opbm9tZSBRdcOpYmVjIElubm92YXRpb24gQ2VudHJlLCBNY0dpbGwgVW5pdmVyc2l0eSwgTW9udHJlYWwsIENhbmFkYQpkYXRlOiAnYHIgZm9ybWF0KFN5cy5EYXRlKCksICIlWS0lQi0lZCIpYCcKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICBkZl9wcmludDogcGFnZWQKICAgIGNvZGVfZm9sZGluZzogc2hvdwogICAgdG9jOiBubwogICAgdG9jX2Zsb2F0OiAKICAgICAgY29sbGFwc2VkOiBmYWxzZQogICAgICBzbW9vdGhfc2Nyb2xsOiBmYWxzZQotLS0KCiMgUHJlcGFyZSBhbmFseXNpcyB3b3JrZmxvdwoKIyMjIFNldCBwYXJhbWV0ZXJzCgpgYGB7ciBzZXR1cH0Ka25pdHI6Om9wdHNfa25pdCRzZXQocm9vdC5kaXIgPSBycHJvanJvb3Q6OmZpbmRfcnN0dWRpb19yb290X2ZpbGUoKSwKICAgICAgICAgICAgICAgICAgICAgZmlnLndpZHRoPTE1LAogICAgICAgICAgICAgICAgICAgICBkaWdpdD01LAogICAgICAgICAgICAgICAgICAgICBzY2lwZW49OCkKb3B0aW9ucyhyZWFkci5zaG93X3Byb2dyZXNzID0gRkFMU0UsCiAgICAgICAgZGlnaXRzPTUsIAogICAgICAgIHNjaXBlbj04LAogICAgICAgIGZ1dHVyZS5nbG9iYWxzLm1heFNpemUgPSArSW5mKQpgYGAKCgojIyMgU2V0IGZpbGVwYXRocyBhbmQgcGFyYW1ldGVycwoKYGBge3J9CnByb2plY3RfZGlyIDwtIHJwcm9qcm9vdDo6ZmluZF9yc3R1ZGlvX3Jvb3RfZmlsZSgpCgppZihpcy5udWxsKHByb2plY3RfZGlyKSl7CiAgcHJvamVjdF9kaXIgPC0gZ2V0d2QoKQogIHdhcm5pbmcoc3ByaW50ZigiTm8gcnN0dWRpbyBwcm9qZWN0IHJvb3QgZmlsZSAgZm91bmQuIAogICAgICAgICAgICAgICAgICBTZXR0aW5nIHByb2plY3QgZGlyZWN0b3J5IHRvIGN1cnJlbnQgd29ya2Zsb3cuUm1kIGZpbGUgbG9jYXRpb246ICVzLiAKICAgICAgICAgICAgICAgICAgT3ZlcnJpZGUgaWYgbmVlZGVkLiIsCiAgICAgICAgICAgICAgICAgIHByb2plY3RfZGlyKSkKIAp9Cm1lc3NhZ2Uoc3ByaW50ZigiUHJvamVjdCBkaXJlY3Rvcnk6ICVzIiwKICAgICAgICAgICAgICAgIHByb2plY3RfZGlyKSkKYGBgCgojIyMgTG9hZCBsaWJyYXJpZXMKCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmxpYnJhcnkocmhkZjUpCiNsaWJyYXJ5KERyb3BsZXRVdGlscykgIyBpbnN0YWxsIGJ1dCBub3QgbG9hZApsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShtYXRyaXhTdGF0cykKbGlicmFyeShicm9vbSkKbGlicmFyeShmdXJycikKbGlicmFyeSh0aWN0b2MpCmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeShjb3dwbG90KQpwbGFuKG11bHRpcHJvY2VzcykKYGBgCgoKIyMjIExvYWQgZnVuY3Rpb25zCgoKYGBge3IgbWVzc2FnZT1GQUxTRX0KY29kZV9kaXIgPC0gZmlsZS5wYXRoKHByb2plY3RfZGlyLCAiY29kZSIpCnNvdXJjZShmaWxlLnBhdGgoY29kZV9kaXIsICIxX2NyZWF0ZV9qb2luZWRfY291bnRzX3RhYmxlLlIiKSkKc291cmNlKGZpbGUucGF0aChjb2RlX2RpciwgIjJfY3JlYXRlX2NvdW50c19ieV9vdXRjb21lX3RhYmxlLlIiKSkKc291cmNlKGZpbGUucGF0aChjb2RlX2RpciwgIjNfZXN0aW1hdGVfc2FtcGxlX2luZGV4X2hvcHBpbmdfcmF0ZS5SIikpCnNvdXJjZShmaWxlLnBhdGgoY29kZV9kaXIsICI0X2NvbXB1dGVfc3VtbWFyeV9zdGF0aXN0aWNzLlIiKSkKc291cmNlKGZpbGUucGF0aChjb2RlX2RpciwgIjVfcmVhc3NpZ25faG9wcGVkX3JlYWRzLlIiKSkKc291cmNlKGZpbGUucGF0aChjb2RlX2RpciwgIjZfcHVyZ2VfcGhhbnRvbV9tb2xlY3VsZXMuUiIpKQpzb3VyY2UoZmlsZS5wYXRoKGNvZGVfZGlyLCAiN19jYWxsX2NlbGxzLlIiKSkKc291cmNlKGZpbGUucGF0aChjb2RlX2RpciwgIjhfc3VtbWFyaXplX3B1cmdlLlIiKSkKc291cmNlKGZpbGUucGF0aChjb2RlX2RpciwgIjlfcGxvdHRpbmdfZnVuY3Rpb25zLlIiKSkKYGBgCgojIyMgTG9hZCBkYXRhCmBgYHtyfQp2YWxpZGF0aW9uX291dHB1dF9kaXIgPC0gZmlsZS5wYXRoKHByb2plY3RfZGlyLCAiZGF0YSIsICJoaXNlcTQwMDBfdmFsaWRhdGlvbiIpCmBgYAoKYGBge3J9CmRhdGEgPC0gcmVhZF90c3YoZmlsZS5wYXRoKHZhbGlkYXRpb25fb3V0cHV0X2RpciwKICAgICAgICAgICAgICAgICAgICAgICAgICAgImhpc2VxNDAwMF9pbm5lcl9qb2luZWRfd2l0aF9sYWJlbHMudHh0IikpCmRhdGEKYGBgCgoKIyBDb21wdXRlIEluZGV4IGhvcHBpbmcgcmF0ZQoKIyMgIEVzdGltYXRlcyBjb25kaXRpb25hbCBvbiBkdXBsaWNhdGlvbiBsZXZlbCBhbmQgbGFiZWwKCmBgYHtyfQpzdW1tYXJ5X2NvdW50cyA8LQogICAgZGF0YSAlPiUKICBtdXRhdGUocj1hcy5pbnRlZ2VyKHJvd1N1bXMoLltjKDYsNyldKSkpICU+JQogIGFycmFuZ2UocikgJT4lCiAgZ3JvdXBfYnkociwgbGFiZWwpICU+JQogIHN1bW1hcml6ZV9hdCh2YXJzKG1hdGNoZXMoIl5zMV9wbHxzMl9wbCIpKSwgCiAgICAgICAgICAgICAgIGxpc3QofiBzdW0oLikpKSAgJT4lCiAgbXV0YXRlKHMxX2hvcHBlZD1pZl9lbHNlKGxhYmVsICVpbiUgYygiMCxmIiwgInIsZiIpLCBzMl9wbGV4ZWQsMCksCiAgICAgICAgIHMyX2hvcHBlZD1pZl9lbHNlKGxhYmVsICVpbiUgYygiZiwwIiwgImYsciIpLCBzMV9wbGV4ZWQsMCksCiAgICAgICAgIHMxX25vbmhvcHBlZD1pZl9lbHNlKGxhYmVsICVpbiUgYygiciwwIiwgInIsZiIpLCBzMV9wbGV4ZWQsMCksCiAgICAgICAgIHMyX25vbmhvcHBlZD1pZl9lbHNlKGxhYmVsICVpbiUgYygiMCxyIiwgImYsciIpLCBzMl9wbGV4ZWQsMCkpICU+JQogIHNlbGVjdCgtczFfcGxleGVkLC0gczJfcGxleGVkKQoKc3VtbWFyeV9jb3VudHMKYGBgCgojIyBFc3RpbWF0ZXMgY29uZGl0aW9uYWwgb24gZHVwbGljYXRpb24gbGV2ZWwKCmBgYHtyfQpzdW1tYXJ5X2NvdW50c19jb25kaXRpb25hbCA8LQpzdW1tYXJ5X2NvdW50cyAgJT4lCiAgZ3JvdXBfYnkocikgJT4lCiAgc3VtbWFyaXplX2F0KHZhcnMobWF0Y2hlcygiXnMiKSksIAogICAgICAgICAgICAgICBsaXN0KH4gc3VtKC4pKSklPiUKICBtdXRhdGUoU0lIUl8xMiA9IHMxX2hvcHBlZC8oczFfaG9wcGVkK3MxX25vbmhvcHBlZCksCiAgICAgICAgIFNJSFJfMjEgPSBzMl9ob3BwZWQvKHMyX2hvcHBlZCtzMl9ub25ob3BwZWQpLAogICAgICAgICBmcmFjX3MxPSAoczFfaG9wcGVkK3MxX25vbmhvcHBlZCkvIChzMV9ob3BwZWQrczFfbm9uaG9wcGVkK3MyX2hvcHBlZCtzMl9ub25ob3BwZWQpKQpzdW1tYXJ5X2NvdW50c19jb25kaXRpb25hbApgYGAKCgoKIyMgTWFyZ2luYWwgZXN0aW1hdGVzCgpgYGB7cn0Kc3VtbWFyeV9jb3VudHNfbWFyZ2luYWwgPC0gCiAgc3VtbWFyeV9jb3VudHMgJT4lCiAgdW5ncm91cCgpICU+JQogIHN1bW1hcml6ZV9hdCh2YXJzKG1hdGNoZXMoIl5zIikpLCAKICAgICAgICAgICAgICAgbGlzdCh+IHN1bSguKSkpJT4lCiAgbXV0YXRlKFNJSFJfMTIgPSBzMV9ob3BwZWQvKHMxX2hvcHBlZCtzMV9ub25ob3BwZWQpLAogICAgICAgICBTSUhSXzIxID0gczJfaG9wcGVkLyhzMl9ob3BwZWQrczJfbm9uaG9wcGVkKSwKICAgICAgICAgU0lIUj0xLShzMV9ub25ob3BwZWQrczJfbm9uaG9wcGVkKS8oczFfaG9wcGVkK3MxX25vbmhvcHBlZCArczJfaG9wcGVkK3MyX25vbmhvcHBlZCksCiAgICAgICAgIGZyYWNfczE9IChzMV9ob3BwZWQrczFfbm9uaG9wcGVkKS8gKHMxX2hvcHBlZCtzMV9ub25ob3BwZWQrczJfaG9wcGVkK3MyX25vbmhvcHBlZCkpCnN1bW1hcnlfY291bnRzX21hcmdpbmFsCmBgYAoKCmBgYHtyLCBmaWcuaGVpZ2h0PTEwfQpwMSA8LSAKICBnZ3Bsb3Qoc3VtbWFyeV9jb3VudHNfY29uZGl0aW9uYWwpICsKICAgIGdlb21fbGluZShhZXMoeCA9IHIsCiAgICAgICAgICAgICAgICAgICAgeSA9IFNJSFJfMTIqMTAwLAogICAgICAgICAgICAgICBjb2xvdXI9IlNJSFJfMTIiKSkrCiAgICBnZW9tX2xpbmUoYWVzKHggPSByLAogICAgICAgICAgICAgICAgICAgIHkgPSBTSUhSXzIxKjEwMCwKICAgICAgICAgICAgICAgY29sb3VyPSJTSUhSXzIxIikpICsKICAgICAgZ2VvbV9saW5lKGFlcyh4ID0gciwKICAgICAgICAgICAgICAgICAgICB5ID0gZnJhY19zMSwKICAgICAgICAgICAgICAgY29sb3VyPSJmcmFjX3MxIikpICsKICAgIGdlb21faGxpbmUoeWludGVyY2VwdCA9IHVubGlzdChzdW1tYXJ5X2NvdW50c19tYXJnaW5hbFs1OjddKSogYygxMDAsIDEwMCwgMTAwKSwKICAgICAgICAgICAgICAgbGluZXR5cGU9ImRhc2hlZCIpICsKICAgIHhsaW0oMCw5MCkgKyAKICB5bGltKDAsMSkgCnAxCiAKYGBgCgoKIyMgTW9sZWN1bGVzCgpgYGB7cn0Kc3VtbWFyeV9tb2xfY291bnRzIDwtCiAgZGF0YSAlPiUKICBtdXRhdGUocj1hcy5pbnRlZ2VyKHJvd1N1bXMoLltjKDYsNyldKSkpICU+JQogIGFycmFuZ2UocikgJT4lICAKICBtdXRhdGVfYXQodmFycyhtYXRjaGVzKCJecyIpKSwgCiAgICAgICAgICAgIGxpc3QofiBhcy5pbnRlZ2VyKC4hPTApKSkgICU+JQogIGdyb3VwX2J5KHIsIGxhYmVsKSAlPiUKICBzdW1tYXJpemVfYXQodmFycyhtYXRjaGVzKCJeczFfcGx8XnMyX3BsIikpLCAKICAgICAgICAgICAgICAgbGlzdCh+IHN1bSguKSkpICAlPiUKICAgIG11dGF0ZShzMV9waGFudG9tPWlmX2Vsc2UobGFiZWwgJWluJSBjKCIwLGYiLCAicixmIiksIHMyX3BsZXhlZCwwTCksCiAgICAgICAgIHMyX3BoYW50b209aWZfZWxzZShsYWJlbCAlaW4lIGMoImYsMCIsICJmLHIiKSwgczFfcGxleGVkLDBMKSwKICAgICAgICAgczFfcmVhbD1pZl9lbHNlKGxhYmVsICVpbiUgYygiciwwIiwgInIsZiIpLCBzMV9wbGV4ZWQsMEwpLAogICAgICAgICBzMl9yZWFsPWlmX2Vsc2UobGFiZWwgJWluJSBjKCIwLHIiLCAiZixyIiksIHMyX3BsZXhlZCwwTCkpICU+JQogIHNlbGVjdCgtczFfcGxleGVkLC0gczJfcGxleGVkKSAKc3VtbWFyeV9tb2xfY291bnRzIApgYGAKYGBge3J9CnN1bW1hcnlfbW9sX2NvdW50c19tYXJnaW5hbCA8LQogIHN1bW1hcnlfbW9sX2NvdW50cyU+JQogIHVuZ3JvdXAoKSAlPiUKICBzdW1tYXJpemVfYXQodmFycyhtYXRjaGVzKCJecyIpKSwgCiAgICAgICAgICAgICAgIGxpc3QofiBzdW0oLikpKSAlPiUKICBtdXRhdGUocHBtXzEyID0gczFfcGhhbnRvbS8oczFfcGhhbnRvbStzMV9yZWFsKSwjIHByb3AgaG9wcGVkIHBoYW50b20gbW9sZWMKICAgICAgICAgcHBtXzIxID0gczJfcGhhbnRvbS8oczJfcGhhbnRvbStzMl9yZWFsKSwKICAgICAgICAgcHBtXzEgPSBzMl9waGFudG9tLyhzMl9waGFudG9tK3MxX3JlYWwpLCAjIHByb3AgcGhhbnRvbSBtb2xlYwogICAgICAgICBwcG1fMiA9IHMxX3BoYW50b20vKHMxX3BoYW50b20rczJfcmVhbCksCiAgICAgICAgIHBwbT0oczFfcGhhbnRvbStzMl9waGFudG9tKS8oczFfcGhhbnRvbStzMV9yZWFsK3MyX3BoYW50b20rczJfcmVhbCksCiAgICAgICAgIGZyYWNfbW9sX3MxPSAoczFfcGhhbnRvbStzMV9yZWFsKS8gKHMxX3BoYW50b20rczFfcmVhbCtzMl9waGFudG9tK3MyX3JlYWwpKQpzdW1tYXJ5X21vbF9jb3VudHNfbWFyZ2luYWwKYGBgCgoKYGBge3J9CnN1bW1hcnlfbW9sX2NvdW50c19jb25kaXRpb25hbCA8LQogIHN1bW1hcnlfbW9sX2NvdW50cyU+JQogIGdyb3VwX2J5KHIpICU+JQogIHN1bW1hcml6ZV9hdCh2YXJzKG1hdGNoZXMoIl5zIikpLCAKICAgICAgICAgICAgICAgbGlzdCh+IHN1bSguKSkpICU+JQogIG11dGF0ZShwcG1fMTIgPSBzMV9waGFudG9tLyhzMV9waGFudG9tK3MxX3JlYWwpLAogICAgICAgICBwcG1fMjEgPSBzMl9waGFudG9tLyhzMl9waGFudG9tK3MyX3JlYWwpLAogICAgICAgICBwcG09KHMxX3BoYW50b20rczJfcGhhbnRvbSkvKHMxX3BoYW50b20rczFfcmVhbCtzMl9waGFudG9tK3MyX3JlYWwpLAogICAgICAgICBmcmFjX21vbF9zMT0gKHMxX3BoYW50b20rczFfcmVhbCkvIChzMV9waGFudG9tK3MxX3JlYWwrczJfcGhhbnRvbStzMl9yZWFsKSkKc3VtbWFyeV9tb2xfY291bnRzX2NvbmRpdGlvbmFsCmBgYAoKYGBge3IsIGZpZy5oZWlnaHQ9MTB9CnA0IDwtIGdncGxvdChzdW1tYXJ5X21vbF9jb3VudHNfY29uZGl0aW9uYWwpICsKICBnZW9tX2xpbmUoYWVzKHggPSByLAogICAgICAgICAgICAgICAgICB5ID0gcHBtXzEyLAogICAgICAgICAgICAgY29sb3VyPSJwcG1fMTIiKSkrCiAgZ2VvbV9saW5lKGFlcyh4ID0gciwKICAgICAgICAgICAgICAgICAgeSA9ICBwcG1fMjEsCiAgICAgICAgICAgICBjb2xvdXI9IiBwcG1fMjEiKSkgKwogICNnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSB1bmxpc3Qoc3VtbWFyeV9tb2xfY291bnRzX21hcmdpbmFsWzU6N10pLCBsaW5ldHlwZT0iZGFzaGVkIikgKwogIHhsaW0oMCwzMDApIApwNAogI2dnc2F2ZSgicGhhbnRvbV9tb2xlY3VsZXNfdmFsaWRhdGlvbi5wZGYiLCBwNCwgIHdpZHRoID04LCBoZWlnaHQgPSA1KQpgYGAKCgoKIyMgRXhhbWluZSBleHRlbnQgb2YgY29udGFtaW5hdGlvbiBpbiBjZWxscwoKYGBge3J9CnN1bW1hcnlfbW9sX2NvdW50c19jZWxsPC0KICBkYXRhICU+JQogIGZpbHRlcihsYWJlbCE9Ik5BIikgJT4lCiAgbXV0YXRlX2F0KHZhcnMobWF0Y2hlcygiXnMiKSksIAogICAgICAgICAgICBsaXN0KH4gYXMuaW50ZWdlciguIT0wKSkpICAlPiUKICBncm91cF9ieShjZWxsLCBsYWJlbCkgJT4lCiAgc3VtbWFyaXplX2F0KHZhcnMobWF0Y2hlcygiXnMxX3BsfF5zMl9wbCIpKSwgCiAgICAgICAgICAgICAgIGxpc3QofiBzdW0oLikpKSAgJT4lCiAgICBtdXRhdGUoczFfcGhhbnRvbT1pZl9lbHNlKGxhYmVsICVpbiUgYygiZiwwIiwgImYsciIpLCBzMV9wbGV4ZWQsMEwpLAogICAgICAgICBzMl9waGFudG9tPWlmX2Vsc2UobGFiZWwgJWluJSBjKCIwLGYiLCAicixmIiksIHMyX3BsZXhlZCwwTCksCiAgICAgICAgIHMxX3JlYWw9aWZfZWxzZShsYWJlbCAlaW4lIGMoInIsMCIsICJyLGYiKSwgczFfcGxleGVkLDBMKSwKICAgICAgICAgczJfcmVhbD1pZl9lbHNlKGxhYmVsICVpbiUgYygiMCxyIiwgImYsciIpLCBzMl9wbGV4ZWQsMEwpKSAlPiUKICBzZWxlY3QoLXMxX3BsZXhlZCwtIHMyX3BsZXhlZCkgJT4lCiAgZ3JvdXBfYnkoY2VsbCkgJT4lCiAgc3VtbWFyaXplX2F0KHZhcnMobWF0Y2hlcygiXnMiKSksIAogICAgICAgICAgICAgICBsaXN0KH4gc3VtKC4pKSkgJT4lIAogICAgbXV0YXRlX2F0KHZhcnMobWF0Y2hlcygiXnMiKSksIAogICAgICAgICAgICBsaXN0KG5vbmVtcHR5PSB+IGFzLmludGVnZXIoLiE9MCkpKSAgJT4lCiAgdW5pdGUobGFiZWwsbWF0Y2hlcygibm9uZW1wdHkiKSwgc2VwPSIsIikgJT4lCiAgbXV0YXRlKGNlbGxfc3RhdHVzPSAKICAgICAgICAgICBjYXNlX3doZW4obGFiZWwgJWluJSBjKCIwLDAsMCwxIiwiMCwwLDEsMCIpIH4gInJlYWwiLAogICAgICAgICAgICAgICAgICAgICBsYWJlbCAlaW4lIGMoIjAsMCwxLDEiKSB+ICJyZWFsIiwKICAgICAgICAgICAgICAgICAgICAgbGFiZWwgJWluJSBjKCIxLDAsMCwwIiwiMCwxLDAsMCIsICIxLDEsMCwwIikgfiAgInBoYW50b20iLAogICAgICAgICAgICAgICAgICAgICBsYWJlbCAlaW4lIGMoIjEsMCwwLDEiLCIwLDEsMSwwIikgfiAgInBoYW50b20iLAogICAgICAgICAgICAgICAgICAgICBUUlVFIH4gICJjb250YW1pbmF0ZWQiKSkgICU+JQogIG11dGF0ZShzMV90b3RhbCA9KHMxX3BoYW50b20rczFfcmVhbCksCiAgICAgICAgIHMyX3RvdGFsID0oczJfcGhhbnRvbStzMl9yZWFsKSwKICAgICAgICAgczFfcHBtID0gczFfcGhhbnRvbS8oczFfdG90YWwpLAogICAgICAgICBzMl9wcG0gPSBzMl9waGFudG9tLyhzMl90b3RhbCkpCiAgCnN1bW1hcnlfbW9sX2NvdW50c19jZWxsCmBgYApgYGB7cn0KY2VsbF9zdGF0dXNfdGFsbHkgPC0KICBzdW1tYXJ5X21vbF9jb3VudHNfY2VsbCAlPiUKICBncm91cF9ieShjZWxsX3N0YXR1cyxsYWJlbCkgJT4lCiAgdGFsbHkoc29ydD1UUlVFKQpjZWxsX3N0YXR1c190YWxseQpgYGAKCmBgYHtyfQpuX2FmZmVjdGVkX2NlbGxzIDwtIAogIGNlbGxfc3RhdHVzX3RhbGx5ICU+JQogIHVuZ3JvdXAoKSAlPiUKICBmaWx0ZXIoY2VsbF9zdGF0dXMhPSJyZWFsIikgJT4lIAogIHN1bW1hcmlzZShuPXN1bShuKSkgJT4lIAogIHB1bGwobikKbl9hZmZlY3RlZF9jZWxscyA8LSBuX2FmZmVjdGVkX2NlbGxzICsgMTUwMiArNgpuX2FmZmVjdGVkX2NlbGxzCmBgYAoKCmBgYHtyfQpuX3RvdGFsX2NlbGxzIDwtCiAgc3VtbWFyeV9tb2xfY291bnRzX2NlbGwgJT4lCiAgbXV0YXRlX2F0KHZhcnMobWF0Y2hlcygiX3RvdGFsIikpLCAKICAgICAgICAgICAgbGlzdCh+IGFzLmludGVnZXIoLiE9MCkpKSAlPiUKICB1bmdyb3VwKCkgJT4lCiAgc3VtbWFyaXNlX2F0KHZhcnMobWF0Y2hlcygiX3RvdGFsIikpLCBsaXN0KH5zdW0oLikpKSAlPiUKICBtdXRhdGUoY2VsbHNfdG90YWw9IHMxX3RvdGFsK3MyX3RvdGFsKSU+JQogIHB1bGwoY2VsbHNfdG90YWwpCm5fdG90YWxfY2VsbHMgCmBgYAoKUHJvcG9ydGlvbiBvZiBhZmZlY3RlZCBjZWxscwoKYGBge3J9CnBfYWZmZWN0ZWRfY2VsbHMgPC0gIG5fYWZmZWN0ZWRfY2VsbHMvbl90b3RhbF9jZWxscwpwX2FmZmVjdGVkX2NlbGxzCmBgYAoKRm9yIGVhY2ggY2VsbC1iYXJjb2RlLCBwbG90IHRoZSBudW1iZXIgb2YgcGhhbnRvbSBtb2xlY3VsZXMgYWdhaW5zdCB0aGUgbnVtYmVyIG9mIHRvdGFsIG1vbGVjdWxlcyBhc3NvY2lhdGVkIHdpdGggaXQuCgpTYW1wbGUgMSBwbG90CgpgYGB7ciwgZmlnLmhlaWdodD0xMH0KcDIgPC0gZ2dwbG90KHN1bW1hcnlfbW9sX2NvdW50c19jZWxsICU+JQogICAgICAgICAgICAgICBmaWx0ZXIoY2VsbF9zdGF0dXMhPSJyZWFsIiAmIHMxX3BoYW50b20gPjAgKSkgKwogIGdlb21fcG9pbnQoYWVzKHggPSBzMV90b3RhbCwKICAgICAgICAgICAgICAgICAgeSA9IHMxX3BoYW50b20sCiAgICAgICAgICAgICAgICAgY29sb3VyPWNlbGxfc3RhdHVzKSkgKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgc2NhbGVfeV9sb2cxMCgpCgpwMgogI2dnc2F2ZSgicGhhbnRvbV9tb2xlY3VsZXNfdmFsaWRhdGlvbi5wZGYiLCBwMiwgIHdpZHRoID04LCBoZWlnaHQgPSA1KQpgYGAKClNhbXBsZSAyIHBsb3QKCgpgYGB7ciwgZmlnLmhlaWdodD0xMH0KcDMgPC0gZ2dwbG90KHN1bW1hcnlfbW9sX2NvdW50c19jZWxsICU+JQogICAgICAgICAgICAgICBmaWx0ZXIoY2VsbF9zdGF0dXMhPSJyZWFsIiAmIHMyX3BoYW50b20gPjAgKSkgKwogIGdlb21fcG9pbnQoYWVzKHggPSBzMl90b3RhbCwKICAgICAgICAgICAgICAgICAgeSA9IHMyX3BoYW50b20sCiAgICAgICAgICAgICAgICAgY29sb3VyPWNlbGxfc3RhdHVzKSkgKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgc2NhbGVfeV9sb2cxMCgpCgpwMwojZ2dzYXZlKCJwaGFudG9tX2NlbGxzX3ZhbGlkYXRpb25fczIucGRmIiwgcDMsICB3aWR0aCA9OCwgaGVpZ2h0ID0gNSkKYGBgCgojIFJ1biB3b3JrZmxvdyBvbiBtdWx0aXBsZXhlZCBkYXRhCgpgYGB7cn0KcmVhZF9jb3VudHMgPC0gCiAgZGF0YSAlPiUgCiAgZmlsdGVyKGxhYmVsIT0iTkEiKSU+JQogIHNlbGVjdCgtZW5kc193aXRoKCJub25wbGV4ZWQiKSkgJT4lCiAgc2V0X25hbWVzKGMoImNlbGwiLCAiZ2VuZSIsICJ1bWkiLCAiczEiLCAiczIiLCAib3V0Y29tZSIsICJsYWJlbCIpKQpyZWFkX2NvdW50cwpgYGAKCgpgYGB7cn0KUyA8LSAyCnNhbXBsZV9uYW1lcyA8LSBjb2xuYW1lcyhyZWFkX2NvdW50cylbNDooUyszKV0Kc2FtcGxlX25hbWVzCmBgYAoKCgpgYGB7cn0KdGljKCJTdGVwIDI6IGNyZWF0aW5nIG91dGNvbWUgY291bnRzIGRhdGF0YWJsZSB3aXRoIGdyb3VwaW5nIHZhcnMiKQoKb3V0Y29tZV9jb3VudHMgPC0gY3JlYXRlX291dGNvbWVfY291bnRzKHJlYWRfY291bnRzJT4lCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNlbGVjdCgtbGFiZWwpLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNhbXBsZV9uYW1lcywgIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWluX2ZyYWM9MC44KQp0b2MoKQoKCm91dGNvbWVfY291bnRzCmBgYAoKYGBge3J9CnRpYygiU3RlcCAzOiBjcmVhdGluZyBhIGNoaW1lcmEgY291bnRzIGRhdGF0YWJsZSBhbmQgZXN0aW1hdGluZyBob3BwaW5nIHJhdGUiKQogIGZpdF9vdXQgPC0KICAgIGVzdGltYXRlX2hvcHBpbmdfcmF0ZSgKICAgICAgb3V0Y29tZV9jb3VudHMsCiAgICAgIFMKICAgICkKICB0b2MoKQpmaXRfb3V0IApgYGAKCgoKYGBge3J9CiAgIyBjb21wdXRlX21vbGVjdWxhcl9jb21wbGV4aXR5X3Byb2ZpbGUKICB0aWMoIlN0ZXAgNDogY29tcHV0ZSBtb2xlY3VsYXIgY29tcGxleGl0eSBwcm9maWxlIGFuZCBvdGhlciBzdW1tYXJ5IHN0YXRpc3RpY3MiKQogIHN1bW1hcnlfc3RhdHMgPC0KICAgIGNvbXB1dGVfc3VtbWFyeV9zdGF0cygKICAgICAgb3V0Y29tZV9jb3VudHMsCiAgICAgIGZpdF9vdXQkZ2xtX2VzdGltYXRlcyRwaGF0LAogICAgICBzYW1wbGVfbmFtZXMKICAgICkKICB0b2MoKQpzdW1tYXJ5X3N0YXRzCmBgYAoKU2V0IHRoZSB0cmFkZS1vZmYgcmF0aW8gY29zdCBjdXRvZmYgKCp0b3JjKikuIFRoZSBwYXJhbWV0ZXIgKnRvcmMqIHJlcHJlc2VudHMgdGhlIG51bWJlciBvZiByZWFsIG1vbGVjdWxlcyBvbmUgaXMgd2lsbGluZyB0byBpbmNvcnJlY3RseSBkaXNjYXJkIGluIG9yZGVyIHRvIGNvcnJlY3RseSBwdXJnZSBvbmUgcGhhbnRvbSBtb2xlY3VsZS4gU2luY2UgZGlzY2FyZGluZyBhIGxhcmdlIHByb3BvcnRpb24gb2YgdGhlIGRhdGEgaXMgdW5kZXNpcmFibGUsIHJlYXNvbmFibGUgdmFsdWVzIG9mICp0b3JjKiBhcmUgZXhwZWN0ZWQgdG8gYmUgd2l0aGluIHRoZSByYW5nZSBvZiAxLTUuCgpgYGB7cn0KdG9yYyA8LSAzIApgYGAKCgpgYGB7cn0KdGljKCJTdGVwIDU6IHJlYXNzaWduIHJlYWQgY291bnRzLCBkZXRlcm1pbmUgY3V0b2ZmLCBhbmQgbWFyayByZXRhaW5lZCBvYnNlcnZhdGlvbnMiKQoKICBvdXRjb21lX2NvdW50cyA8LQogICAgcmVhc3NpZ25fcmVhZHNfYW5kX21hcmtfcmV0YWluZWRfb2JzZXJ2YXRpb25zKAogICAgICBvdXRjb21lX2NvdW50cywKICAgICAgc3VtbWFyeV9zdGF0cywKICAgICAgc2FtcGxlX25hbWVzLAogICAgICBmaXRfb3V0LAogICAgICB0b3JjCiAgICApCiAgIyBnZXQgdGhlIHRyYWRvZmYgcmF0aW8gY3V0b2ZmCiAgc3VtbWFyeV9zdGF0cyA8LSBnZXRfdGhyZXNob2xkKG91dGNvbWVfY291bnRzLCBzdW1tYXJ5X3N0YXRzKQoKICB0b2MoKQpgYGAKCgoKYGBge3J9CnRpYygiU3RlcCA2OiBQdXJnZSBhbmQgc2F2ZSByZWFkIGNvdW50cyBkYXRhdGFibGUgdG8gZGlzayIpCgpyZWFkX2NvdW50cyA8LQogIGxlZnRfam9pbihyZWFkX2NvdW50cyAlPiUKICAgIHNlbGVjdChvdXRjb21lLCBjZWxsLCB1bWksIGdlbmUsIHNhbXBsZV9uYW1lcywgbGFiZWwpLAogIG91dGNvbWVfY291bnRzLAogIGJ5ID0gYygib3V0Y29tZSIpCiAgKSAlPiUKICBzZWxlY3QoLW91dGNvbWUpCgp0b2MoKQoKYGBgCgoKIyMjIENvbXBhcmUgdGhlIFNJSFIgZXN0aW1hdGVzIHdpdGggZ3JvdW5kIHRydXRoIGVzdGltYXRlcwoKYGBge3IsIGZpZy53aWR0aD0xMH0KcDUgPC0gCiAgZ2dwbG90KHN1bW1hcnlfY291bnRzX2NvbmRpdGlvbmFsKSArCiAgZ2VvbV9saW5lKGFlcyh4ID0gciwKICAgICAgICAgICAgICAgICAgeSA9IFNJSFJfMTIsCiAgICAgICAgICAgICBjb2xvdXI9IjEyIikpKwogIGdlb21fbGluZShhZXMoeCA9IHIsCiAgICAgICAgICAgICAgICAgIHkgPSBTSUhSXzIxLAogICAgICAgICAgICAgY29sb3VyPSIyMSIpKSAgKwpnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gIHN1bW1hcnlfY291bnRzX21hcmdpbmFsJFNJSFIsIAogICAgICAgICAgICAgICBjb2xvdXI9InRydWUgbWVhbiIpLCAKICAgICAgICAgICBsaW5ldHlwZT0ic29saWQiLAogICAgICAgICAgIHNpemU9LjUpICAgKwpnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gIGZpdF9vdXQkZ2xtX2VzdGltYXRlcyRTSUhSLCAKICAgICAgICAgICAgICAgY29sb3VyPSJlc3RpbWF0ZSIpLAogICAgICAgICAgIGxpbmV0eXBlPSJzb2xpZCIsCiAgICAgICAgICAgc2l6ZT0uMSkgICArCiAgICBnZW9tX2xpbmVyYW5nZShkYXRhPXN1bW1hcnlfY291bnRzX2NvbmRpdGlvbmFsLCAKICAgICAgICAgICAgICAgICAgYWVzKHg9ciwKICAgICAgICAgICAgICAgICAgICAgIHltYXg9MS1maXRfb3V0JGdsbV9lc3RpbWF0ZXMkcGhhdF9sb3csCiAgICAgICAgICAgICAgICAgICAgICB5bWluPTEtZml0X291dCRnbG1fZXN0aW1hdGVzJHBoYXRfaGlnaCwKICAgICAgICAgICAgICAgICAgICAgIGNvbG91cj0iZXN0aW1hdGUiKSwgCiAgICAgICAgICAgICAgICAgIHNpemU9LjUpKyAKICB4bGltKDEsMjEwKSArCiAgeWxpbSgwLjAwMiwwLjAwNSkKI2dnc2F2ZSgiaW5kZXhfaG9wcGluZ19yYXRlXzIwMC5wZGYiLCBwNSwgd2lkdGg9OSwgaGVpZ2h0PTYpCnA1CmBgYAoKCgojIERldGVybWluZSB0aGUgbnVtYmVyIG9mIGZhbHNlIHBvc2l0aXZlcyBhbmQgZmFsc2UgbmVnYXRpdmVzCgpgYGB7cn0KcmVhZF9jb3VudHMgPC0KICAgIHJlYWRfY291bnRzICU+JQogIHVuZ3JvdXAoKSAlPiUKICBhcnJhbmdlKC1xcikgJT4lCiAgbXV0YXRlKHQ9IGNhc2Vfd2hlbigKICAgICAgbGFiZWwgJWluJSBjKCJmLHIiLCIwLHIiKSB+IDIsCiAgICAgIGxhYmVsICVpbiUgYygicixmIiwiciwwIikgfiAxICAgICksCiAgICBmPSBjYXNlX3doZW4oCiAgICAgIGxhYmVsICVpbiUgYygiZixyIiwiZiwwIikgfiAxLAogICAgICBsYWJlbCAlaW4lIGMoInIsZiIsIjAsZiIpIH4gMiAgICApKSAlPiUKICBtdXRhdGUodHA9IGlmX2Vsc2UoIHQgPT0gcywgMUwsIDBMLCBtaXNzaW5nID0wTCksCiAgICAgICAgIGZwPSBpZl9lbHNlKCBmID09IHMsIDFMLCAwTCwgbWlzc2luZyA9MEwpLAogICAgICAgICB0bj0gaWZfZWxzZSggZiAhPSBzLCAxTCwgMEwsIG1pc3NpbmcgPTBMKSwKICAgICAgICAgZm49IGlmX2Vsc2UoIHQgIT0gcywgMUwsIDBMLCBtaXNzaW5nID0wTCksCiAgICAgICAgIHRwMD0gaWZfZWxzZSggdCA9PSAwLCAxTCwgMEwsIG1pc3NpbmcgPTBMKSwgIzAgaWYgIHByZWRpY3QgYWxsIG1vbGVjdWxlcyB0byBiZSBwaGFudG9tCiAgICAgICAgIGZuMD0gaWZfZWxzZSggdCAhPSAwLCAxTCwgMEwsIG1pc3NpbmcgPTBMKSwKICAgICAgICAgdG4wPSBpZl9lbHNlKCBmICE9IDAsIDFMLCAwTCwgbWlzc2luZyA9MEwpLAogICAgICAgICBmcDA9IGlmX2Vsc2UoIGYgPT0gMCwgMUwsIDBMLCBtaXNzaW5nID0wTCksCiAgICAgICAgIHRwX21heD0gaWZfZWxzZSggdCA9PSBzX21heHByb3AsIDFMLCAwTCwgbWlzc2luZyA9MEwpLAogICAgICAgICBmcF9tYXg9IGlmX2Vsc2UoIGYgPT0gc19tYXhwcm9wLCAxTCwgMEwsIG1pc3NpbmcgPTBMKSwKICAgICAgICAgdG5fbWF4PSBpZl9lbHNlKCBmICE9IHNfbWF4cHJvcCwgMUwsIDBMLCBtaXNzaW5nID0wTCksCiAgICAgICAgIGZuX21heD0gaWZfZWxzZSggdCAhPSBzX21heHByb3AsIDFMLCAwTCwgbWlzc2luZyA9MEwpKSAKYGBgCgoKIyMjIFRoZSBtYXhpbXVtIHJlYWQgZnJhY3Rpb24gbWV0aG9kCgoKYGBge3J9CmZhbHNlX2NvdW50c19tYXhwcm9wIDwtCiAgcmVhZF9jb3VudHMgJT4lCiAgc3VtbWFyaXplX2F0KHZhcnMoYygidHBfbWF4IiwgImZwX21heCIsICJ0bl9tYXgiLCAiZm5fbWF4IikpLAogICAgICAgICAgICBsaXN0KCB+IHN1bSguKSkpICU+JQogIHNldF9uYW1lcyhjKCJ0cCIsICJmcCIsICJ0biIsICJmbiIpKQpmYWxzZV9jb3VudHNfbWF4cHJvcApgYGAKCiMjIyBUaGUgdG9yIG1ldGhvZAoKCmBgYHtyfQpmYWxzZV9jb3VudHNfbWluX2N1dG9mZiA8LQogIHJlYWRfY291bnRzICU+JQogIHN1bW1hcml6ZV9hdCh2YXJzKGMoInRwIiwgImZwIiwgInRuIiwgImZuIikpLAogICAgICAgICAgICBsaXN0KCB+IHN1bSguKSkpIApmYWxzZV9jb3VudHNfbWluX2N1dG9mZgpgYGAKCiMjIyBObyBwdXJnaW5nCgpgYGB7cn0KZmFsc2VfY291bnRzX25vcHVyZ2luZyA8LSAKICByZWFkX2NvdW50cyU+JQogIHN1bW1hcml6ZShuX2N1Z3MgPSBuKCksCiAgICAgICAgICAgIG5fcmVhbD0gc3VtKHQ+MCwgbmEucm0gPSBUUlVFKSwKICAgICAgICAgICAgbl9mYW50b20gPSBzdW0oZj4wLCBuYS5ybSA9IFRSVUUpLAogICAgICAgICAgICBuX21vbD1uX3JlYWwrbl9mYW50b20sCiAgICAgICAgICAgIGcgPW5fY3Vncy0gbl9yZWFsLAogICAgICAgICAgICB1ID0gbl9tb2wtbl9jdWdzLAogICAgICAgICAgICB0cD1uX21vbC1nLAogICAgICAgICAgICBmcD11K2csCiAgICAgICAgICAgIHRuPTAsCiAgICAgICAgICAgIGZuPTApCgpmYWxzZV9jb3VudHNfbm9wdXJnaW5nCmBgYAoKCiMjIyAgVE9SIGN1dG9mZgoKYGBge3J9CnJlYWRfY291bnRzIDwtCiAgcmVhZF9jb3VudHMgJT4lCiAgbXV0YXRlX2F0KHZhcnMoYygidHAiLCAiZnAiLCAidG4iLCAiZm4iLCJ0cDAiLCAiZnAwIiwgInRuMCIsICJmbjAiKSksCiAgICAgICAgICAgIGxpc3QoY3VtPSB+IGN1bXN1bSguKSkpICAlPiUKICBtdXRhdGVfYXQodmFycyhjKCJ0cF9jdW0iLCAiZnBfY3VtIiwgInRuX2N1bSIsICJmbl9jdW0iKSksIAogICAgICAgICAgICBsaXN0KCB+IChsYXN0KC4pLWxhZyguLCBkZWZhdWx0ID0wKSkpKSAlPiUKICBtdXRhdGUodHBfdD10cF9jdW0gKyB0cDBfY3VtLAogICAgICAgICBmcF90PWZwX2N1bSArIGZwMF9jdW0sCiAgICAgICAgIHRuX3Q9dG5fY3VtICsgdG4wX2N1bSwKICAgICAgICAgZm5fdD1mbl9jdW0gKyBmbjBfY3VtLAogICAgICAgICBmcG09IGZpcnN0KGZwX3QpLSBmcF90LCAKICAgICAgICAgZm5tPSBmbl90LWZpcnN0KGZuX3QpLAogICAgICAgICB0b3JfdHJ1ZT0gZm5tL2ZwbSkKYGBgCgoKCmBgYHtyfQpmYWxzZV9jb3VudHNfdG9yX2N1dG9mZiA8LQogIHJlYWRfY291bnRzJT4lCiAgZmlsdGVyKHJldGFpbikgJT4lCiAgc2xpY2UoMSklPiUKICBzZWxlY3QoYygiczEiLCAiczIiLCAicXIiLCAidG9yIiwgInRwX3QiLCAiZnBfdCIsICJ0bl90IiwgImZuX3QiLCAiZnBtIiwgImZubSIsICJ0b3JfdHJ1ZSIpKSAgCmZhbHNlX2NvdW50c190b3JfY3V0b2ZmCmBgYAojIyBDcmVhdGUgY29tcGFyaXNvbiBkYXRhdGFibGUKCmBgYHtyfQpmYWxzZV9jb3VudHNfZHQgPC0KICBiaW5kX3Jvd3MoCiAgICBsaXN0KG5vX3B1cmdpbmc9ZmFsc2VfY291bnRzX25vcHVyZ2luZyAlPiUKICAgICAgICAgICBzZWxlY3QoYygidHAiLCAiZnAiLCAidG4iLCAiZm4iKSksCiAgICAgICAgIG5vX2N1dG9mZj1mYWxzZV9jb3VudHNfbWluX2N1dG9mZiwKICAgICAgICAgdG9yX2N1dG9mZj0KICAgICAgICAgICBmYWxzZV9jb3VudHNfdG9yX2N1dG9mZiAlPiUKICAgICAgICAgICBzZWxlY3QoYygidHBfdCIsICJmcF90IiwgInRuX3QiLCAiZm5fdCIpKSAlPiUKICAgICAgICAgICBzZXRfbmFtZXMoYygidHAiLCAiZnAiLCAidG4iLCAiZm4iKSksCiAgICAgICAgIG1heF9mcmFjPWZhbHNlX2NvdW50c19tYXhwcm9wKSwKICAgIC5pZD0iYXBwcm9hY2giKSAlPiUKICBzZWxlY3QoYXBwcm9hY2gsIGZwLGZuLCB0cCwgdG4pICU+JQogIG11dGF0ZShmcHI9ZnAvZmFsc2VfY291bnRzX25vcHVyZ2luZyRuX2ZhbnRvbSwKICAgICAgICAgZm5yPWZuL2ZhbHNlX2NvdW50c19ub3B1cmdpbmckbl9yZWFsKQpmYWxzZV9jb3VudHNfZHQKYGBgCgoKCiMjIFBsb3RzCgojIyMgRGF0YXRhYmxlIGZvciBwbG90dGluZwoKCmBgYHtyfQpjbGFzc2lmaWNhdGlvbl9jdXJ2ZXMgPC0KICByZWFkX2NvdW50cyAlPiUgCiAgZ3JvdXBfYnkocXIpICU+JQogIHNsaWNlKDFMKSAlPiUKICB1bmdyb3VwKCkgJT4lCiAgc2VsZWN0KCBxciwgcXMsIHRvciwgcmV0YWluLCBmcF90LCBGUCwgZm5fdCwgRk4sIHRwX3QsIHRuX3QsIFRQLCBUTiwgRlBtLCBGTm0sIGZwbSwgZm5tLCB0b3JfdHJ1ZSxvLHIpICU+JQogIG11dGF0ZShmcHI9ZnBfdC9mYWxzZV9jb3VudHNfbm9wdXJnaW5nJG5fZmFudG9tLAogICAgICAgICBmbnI9Zm5fdC9mYWxzZV9jb3VudHNfbm9wdXJnaW5nJG5fcmVhbCkKY2xhc3NpZmljYXRpb25fY3VydmVzCmBgYAoKIyMjIFByZWZvcm1hbmNlIFBsb3RzCgoKCmBgYHtyIGZpZy5oZWlnaHQ9NywgZmlnLndpZHRoPTEwLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpwX3RyYWRlb2ZmIDwtICAKICAgIGdncGxvdChjbGFzc2lmaWNhdGlvbl9jdXJ2ZXMpICsgCiAgICBnZW9tX3BvaW50KAogICAgICBhZXMoeCA9IEZQbSwKICAgICAgICAgIHkgPSBGTm0pLAogICAgICBzaXplPS41KSsKICAgIGdlb21fbGluZSgKICAgICAgYWVzKHggPSBGUG0sCiAgICAgICAgICB5ID0gRlBtLAogICAgICAgICAgY29sb3VyPSIxIikKICAgICkgKwogICAgZ2VvbV9saW5lKAogICAgICBhZXMoeCA9IEZQbSwKICAgICAgICAgIHkgPSAyKkZQbSwKICAgICAgICAgIGNvbG91cj0iMiIpKSsKICAgIGdlb21fbGluZSgKICAgICAgYWVzKHggPSBGUG0sCiAgICAgICAgICB5ID0gMypGUG0sCiAgICAgICAgICBjb2xvdXI9IjMiKSkrCiAgICBnZW9tX2xpbmUoCiAgICAgIGFlcyh4ID0gRlBtLAogICAgICAgICAgeSA9IDQqRlBtLAogICAgICAgICAgY29sb3VyPSI0IikpKwogICAgZ2VvbV9saW5lKAogICAgICBhZXMoeCA9IEZQbSwKICAgICAgICAgIHkgPSA1KkZQbSwKICAgICAgICAgIGNvbG91cj0iNSIpKSsKICAgIGdlb21fbGluZSgKICAgICAgYWVzKHggPSBGUG0sCiAgICAgICAgICB5ID0gOSpGUG0sCiAgICAgICAgICBjb2xvdXI9IjkiKSkrCiAgICBzY2FsZV95X2xvZzEwKCkgKwogICAgdGhlbWVfYncoKSAgKwogICAgICB0aGVtZSgKICAgICAgICBsZWdlbmQudGl0bGUgPSBlbGVtZW50X3RleHQoZmFjZSA9ICJib2xkIikpICsgCiAgICAgIHNjYWxlX2NvbG9yX2Rpc2NyZXRlKG5hbWUgPSAiVE9SQyIpICsKICAgIGxhYnMoeD0iTWFyZ2luYWwgRGVjcmVhc2UgaW4gRmFsc2UgUG9zaXRpdmVzIChyZWR1Y2UgcGhhbnRvbSBtb2xlY3MpICIsCiAgICAgICAgIHk9Ik1hcmdpbmFsIEluY3JlYXNlIGluIEZhbHNlIE5lZ2F0aXZlcyAoZGlzY2FyZCByZWFsIG1vbGVjcykiKSAKICAgIAoKI2dnc2F2ZShmaWxlLnBhdGgoZmlndXJlc19kaXIsICJ2YWxpZGF0aW9uX3RyYWRlb2ZmLnBkZiIpLCBwX3RyYWRlb2ZmLCB3aWR0aD05LCBoZWlnaHQ9NikKCnBfdHJhZGVvZmYKYGBgCgpgYGB7ciBmaWcuaGVpZ2h0PTZ9CnA2IDwtCiAgZ2dwbG90KGNsYXNzaWZpY2F0aW9uX2N1cnZlcykgICsgCiAgICBnZW9tX3BvaW50KAogICAgICAgICAgICAgYWVzKHggPSBmcF90LAogICAgICAgICAgICAgICAgICB5ID0gZm5fdCwKICAgICAgICAgICAgICAgICBjb2xvdXI9InRydWUiKSkgKwoKICBnZW9tX2xpbmUoCiAgICAgICAgICAgICBhZXMoeCA9IGZwX3QsCiAgICAgICAgICAgICAgICAgIHkgPSBmbl90LAogICAgICAgICAgICAgICAgIGNvbG91cj0idHJ1ZSIpKSArCiAgICBnZW9tX2xpbmUoCiAgICAgICAgICAgICBhZXMoeCA9IEZQLAogICAgICAgICAgICAgICAgICB5ID0gRk4sCiAgICAgICAgICAgICAgICAgY29sb3VyPSJwcmVkaWN0ZWQiKSkrCiAgICAgIGdlb21fcG9pbnQoCiAgICAgICAgICAgICBhZXMoeCA9IEZQLAogICAgICAgICAgICAgICAgICB5ID0gRk4sCiAgICAgICAgICAgICAgICAgY29sb3VyPSJwcmVkaWN0ZWQiKSkrIAoKICBnZW9tX3BvaW50KGRhdGE9ZmFsc2VfY291bnRzX2R0ICwKICAgICAgICAgICAgIGFlcyh4ID0gZnAsCiAgICAgICAgICAgICAgICAgIHkgPSBmbiwKICAgICAgICAgICAgICAgICAgc2hhcGU9YXBwcm9hY2gpLAogICAgICAgICAgICAgc2l6ZT0yKSArCiAgICBsYWJzKHg9IkZhbHNlIFBvc2l0aXZlIENvdW50IiwKICAgICAgIHk9IkZhbHNlIE5lZ2F0aXZlIENvdW50IikgICsgCiAgICBzY2FsZV95X3NxcnQoKSArIAogIHNjYWxlX3hfc3FydCgpIAojZ2dzYXZlKCJwZWZvcm1hbmNlX2dyb3VuZHRydXRoLnBkZiIsIHA2LCB3aWR0aD05LCBoZWlnaHQ9NikKIHA2CmBgYAoKCgpgYGB7ciBmaWcuaGVpZ2h0PTh9CnA3IDwtIGdncGxvdChmYWxzZV9jb3VudHNfZHQKICAgICAgICAgICAgICAlPiUgZmlsdGVyKGFwcHJvYWNoICVpbiUgYygibm9fY3V0b2ZmIiwgInRvcl9jdXRvZmYiLCAibWF4X2ZyYWMiKSkpICArIAogIGdlb21fcG9pbnQoCiAgICAgICAgICAgICBhZXMoeCA9IGZwICwKICAgICAgICAgICAgICAgICAgeSA9IGZuLAogICAgICAgICAgICAgICAgICBjb2xvcj1hcHByb2FjaCksCiAgICAgICAgICAgICBzaXplPTIpIAoKCiNnZ3NhdmUoInBlZm9ybWFuY2Vfem9vbS5wZGYiLCBwNywgd2lkdGg9OSwgaGVpZ2h0PTYpCnA3CmBgYApgYGB7cn0Kc2Vzc2lvbkluZm8oKQpgYGAKCg==