Some preliminary notes

Three important notes about this week’s lab exercise:

  1. Since this is week 5, please name the file you submit LabWeek5_firstname_lastname.extension.

  2. The next several weeks1 will involve further exploration of ideas that we cover in class rather than analysis of data sets.

  3. This week’s exercise won’t require you to answer any questions. It will simply ask you to explore (through simulations) some properties of genetic drift. If you show me the explorations I ask for, you’ll get full credit for this exercise, i.e., 10 out of 10 possible points. The exercises for weeks 6 and 7 will be similar, but I’ll probably ask you to stretch yourself a bit and explain the patterns you see.

The idea we want to explore2

In lecture we’ll see that the variance effective size of a population is defined as the size of an ideal population that has the same variance in allele frequency among offspring populations as the actual population we’re interested in, i.e.,

\[ N_e^{(v)} = \frac{p(1-p)}{2\widehat{Var(p)}} \]

We will also see that for a population with separate sexes

\[ N_e \approx \frac{4N_fN_m}{N_f + N_m} \]

We’re going to explore how well the \(N_e\) we expect based on the number of reproducing females and males predicts the \(N_e\) we infer from the observed variance in allele frequency among offspring populations.

The basic simulation

To explore this relationship we will simulate the process of mating and reproduction in a finite population with separate sexes. Here’s the process:

  1. We calculate genotype frequencies in the base population given the allele frequency \(p_0\).

  2. We construct \(N_f\) female genotypes and \(N_m\) male genotypes at random given those genotype frequencies.

  3. We construct allele frequencies in an offspring population by picking one female and one male at random, averaging the allele frequency within them to produce the allele frequency in their offspring, and repeating until we’ve constructed \(N_{pop}\) offspring. Then we take the average allele frequency across all of these offspring and we have the allele frequency in the offspring population.

  4. We then repeat that process to produce \(N_{samp}\) different offspring populations.

Here’s an example using \(p_0 = 0.5\), $N_f = 10, \(N_m = 20\), \(N_{pop} = 100\), and \(N_{samp} = 100\).

library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.4     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   2.0.1     ✓ forcats 0.5.1
── Conflicts ──────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
library(ggplot2)

rm(list = ls())

## NOTE: Don't include this line in your code. I include it so that if 
## you run the R code yourself, you'll get the same results that I do.
## The random number generator in R (and in other programming languages
## for that matter) is a pseudo-random number generator. The numbers it
## produces appear random, but it you give it the same "seed", it will
## generate the same sequence of numbers.
##
set.seed(1234)

simulate <- function(n_f, n_m, p_0, n_samp, n_pop) {
  ## genotype frequenies in Hardy-Weinberg
  ##
  x <- c(p_0^2, 2*p_0*(1-p_0), (1-p_0)^2)
  p <- numeric(n_samp)
  for (i in 1:n_samp) {
    ## construct female genotypes at random given H-W
    ##
    f <- numeric(n_f)
    for (j in 1:n_f) {
      tmp <- rmultinom(1, 1, x)
      f[j] <- which(tmp[, 1] == 1)
    }
    ## construct male genotypes at random given H-W
    ##
    m <- numeric(n_m)
    for (j in 1:n_m) {
      tmp <- rmultinom(1, 1, x)
      m[j] <- which(tmp[, 1] == 1)
    }
    ## construct offspring as average of mom and dad
    ##
    p[i] <- 0
    for (j in 1:n_pop) {
      p_f <- (sample(f, 1) - 1)/2
      p_m <- (sample(m, 1) - 1)/2
      p[i] <- p[i] + (p_f + p_m)/2
    }
    p[i] <- p[i]/n_pop
  }
  return(p)
}

result <- simulate(n_f = 10, n_m = 20, p_0 = 0.5, n_samp = 100, 
                   n_pop = 100)

p <- ggplot(tibble(p = result), aes(x = p)) + 
  geom_histogram(binwidth = 0.01) + 
  theme_bw() +
  ggtitle("One replication of sampling")
p

As you can see, there’s a lot of variation in allele frequency among the offspring populations, but both the mean and the variance of the allele frequency are pretty close to what we’d expect.

p_0 <- 0.5
n_e <- 4*10*20/(10 + 20)

dat <- tibble(Category = c("Expected", "Observed"),
              Mean = c(p_0, round(mean(result), 3)),
              Variance = c(round(p_0*(1-p_0)/(2*n_e), 6),
                           round(var(result), 6)))
dat

Extending the simulation

What we’ve done so far is nice, but we’ve only looked at the results from one simulated base population. Maybe we were lucky in our choice of base population. Let’s see what this looks like if we look at the mean and variance from 1000 simulated base populations. Here we don’t want to look at 1000 graphs like the one above, so I’ll keep track of the mean and variance of allele frequencies from each base population and display histograms of the means and variances and see where the values fall relative to the expectation.3

p_mean <- numeric(1000)
p_var <- numeric(1000)
for (i in 1:1000) {
  p <- simulate(n_f = 10, n_m = 20, p_0 = 0.5, n_samp = 100, 
                n_pop = 100)
  p_mean[i] <- mean(p)
  p_var[i] <- var(p)
}

p_dat <- tibble(p_mean = p_mean, p_var = p_var)
p <- ggplot(p_dat, aes(x = p_mean)) +
  geom_histogram(bins = 25) +
  geom_vline(xintercept = 0.5, linetype = "dashed", 
             color = "salmon") +
  ggtitle("Distribution of mean allele frequencies") +
  theme_bw()
p

p <- ggplot(p_dat, aes(x = p_var)) +
  geom_histogram(bins = 25) +
  geom_vline(xintercept = p_0*(1-p_0)/(2*n_e), linetype = "dashed",
             color = "salmon") +
  ggtitle("Distribution of variance in allele frequencies") +
  theme_bw()
p

There are two striking features of these distributions:

  1. There’s a reasonable amount of variability in both the mean and the variance of allele frequency, and there’s a lot more variability in the variance than in the mean.

  2. The average mean allele frequency across the 1000 simulations is pretty close to the expectation (the vertical dashed line), but the mean variance in allele frequency is noticeably larger in the simulations than we expect.

Exercise 5

The lab exercise this week asks you to explore this discrepancy a bit further. To do so you can use the following function.4

run_simulation <- function(p_0, n_samp, n_pop) {
  dat <- data.frame(n_f = NA, n_m = NA, n_e = NA, n_e_obs = NA, 
                    var_p = NA, var_p_obs = NA)
  for (n_f in c(5, 10, 25, 50, 100)) {
    for (n_m in c(5, 10, 25, 50, 100)) {
      p <- simulate(n_f = n_f, n_m = n_m, p_0 = 0.5, n_samp = n_samp, 
                    n_pop = n_pop)
      ## calculate Ne and Var(p) from formulas
      ##
      n_e <- 4.0*n_f*n_m/(n_f + n_m)
      var_p <- p_0*(1 - p_0)/(2*n_e)
      ## calculate observed Ne from observed variance
      ##
      n_e_obs <- 0.5*0.5/(2*var(p))
      dat <- add_row(dat, n_f = n_f, n_m = n_m, n_e = n_e, 
                     n_e_obs = n_e_obs, var_p = var_p, 
                     var_p_obs = var(p))
    }
  }
  dat <- subset(dat, !is.na(n_f))

  p <- ggplot(dat, aes(x = n_e, y = n_e_obs)) +
    geom_point() +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed",
                color = "salmon") +
    theme_bw()
  print(p)
  p <- ggplot(dat, aes(x = var_p, y = var_p_obs)) +
    geom_point() +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed",
                color = "salmon") +
    theme_bw()
  print(p)
}
  1. Execute run_simulation() with \(p_0 = 0.5\), \(n_{samp} = 100\), and \(n_obs = 100\). The dashed line is the 1:1 line. If the simulations matched expectations, all of the points would lie on the line (with a bit of statistical sampling error). With these parameters, you’ll see that there’s a large mismatch between the observed and expected \(N_e\). There’s a smaller mismatch between the observed and expected variance, but it’s consistently in one direction.

  2. Try some different combinations of \(p_0\), \(n_{samp}\), and \(n_{obs}\) and report your results. I am looking for a minimum of 10 combinations, but if you’d like to try more, feel free. Either way include the plots that your simulations produce. Treat your exploration of different combinations of parameters like an experiment, i.e., only change one at a time (at least at first) so that you can isolate the effect of changes in that parameter. If you’re feeling ambitious you can try to look at combinations, but that could become challenging pretty quickly. See if you can find a combination of parameters where the fit of observed to expected is reasonably good.

  3. Summarize any overall pattern you see. In particular, how do different choices for the parameters affect the extent of the discrepancy you see, if at all?

We’ll discuss any general patterns that emerge and explore explanations for them in lecture on 5 October.


  1. Until we get to week 8↩︎

  2. More precisely, the idea I’m assigning you to explore.↩︎

  3. Note: If you run the following block of code, don’t be worried if you don’t see anything happening. It takes between 3-5 minutes to run on my MacBook↩︎

  4. Note that this function calls simulate(), so you’ll need to execute simulate() in your R notebook before you run this function.↩︎

LS0tCnRpdGxlOiAiRXhwbG9yaW5nIHRoZSBwcm9wZXJ0aWVzIG9mIGVmZmVjdGl2ZSBwb3B1bGF0aW9uIHNpemUiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMjIFNvbWUgcHJlbGltaW5hcnkgbm90ZXMKClRocmVlIGltcG9ydGFudCBub3RlcyBhYm91dCB0aGlzIHdlZWsncyBsYWIgZXhlcmNpc2U6CgoxLiBTaW5jZSB0aGlzIGlzIHdlZWsgNSwgcGxlYXNlIG5hbWUgdGhlIGZpbGUgeW91IHN1Ym1pdCBgTGFiV2VlazVfZmlyc3RuYW1lX2xhc3RuYW1lLmV4dGVuc2lvbmAuCgoyLiBUaGUgbmV4dCBzZXZlcmFsIHdlZWtzXltVbnRpbCB3ZSBnZXQgdG8gd2VlayA4XSB3aWxsIGludm9sdmUgZnVydGhlciBleHBsb3JhdGlvbiBvZiBpZGVhcyB0aGF0IHdlIGNvdmVyIGluIGNsYXNzIHJhdGhlciB0aGFuIGFuYWx5c2lzIG9mIGRhdGEgc2V0cy4KCjMuIFRoaXMgd2VlaydzIGV4ZXJjaXNlIHdvbid0IHJlcXVpcmUgeW91IHRvIGFuc3dlciBhbnkgcXVlc3Rpb25zLiBJdCB3aWxsIHNpbXBseSBhc2sgeW91IHRvIGV4cGxvcmUgKHRocm91Z2ggc2ltdWxhdGlvbnMpIHNvbWUgcHJvcGVydGllcyBvZiBnZW5ldGljIGRyaWZ0LiBJZiB5b3Ugc2hvdyBtZSB0aGUgZXhwbG9yYXRpb25zIEkgYXNrIGZvciwgeW91J2xsIGdldCBmdWxsIGNyZWRpdCBmb3IgdGhpcyBleGVyY2lzZSwgaS5lLiwgMTAgb3V0IG9mIDEwIHBvc3NpYmxlIHBvaW50cy4gVGhlIGV4ZXJjaXNlcyBmb3Igd2Vla3MgNiBhbmQgNyB3aWxsIGJlIHNpbWlsYXIsIGJ1dCBJJ2xsIHByb2JhYmx5IGFzayB5b3UgdG8gc3RyZXRjaCB5b3Vyc2VsZiBhIGJpdCBhbmQgZXhwbGFpbiB0aGUgcGF0dGVybnMgeW91IHNlZS4KCiMjIFRoZSBpZGVhIHdlIHdhbnQgdG8gZXhwbG9yZV5bTW9yZSBwcmVjaXNlbHksIHRoZSBpZGVhIEknbSAqYXNzaWduaW5nKiB5b3UgdG8gZXhwbG9yZS5dCgpJbiBsZWN0dXJlIHdlJ2xsIHNlZSB0aGF0IHRoZSB2YXJpYW5jZSBlZmZlY3RpdmUgc2l6ZSBvZiBhIHBvcHVsYXRpb24gaXMgZGVmaW5lZCBhcyB0aGUgc2l6ZSBvZiBhbiBpZGVhbCBwb3B1bGF0aW9uIHRoYXQgaGFzIHRoZSBzYW1lIHZhcmlhbmNlIGluIGFsbGVsZSBmcmVxdWVuY3kgYW1vbmcgb2Zmc3ByaW5nIHBvcHVsYXRpb25zIGFzIHRoZSBhY3R1YWwgcG9wdWxhdGlvbiB3ZSdyZSBpbnRlcmVzdGVkIGluLCBpLmUuLAoKJCQKTl9lXnsodil9ID0gXGZyYWN7cCgxLXApfXsyXHdpZGVoYXR7VmFyKHApfX0KJCQKCldlIHdpbGwgYWxzbyBzZWUgdGhhdCBmb3IgYSBwb3B1bGF0aW9uIHdpdGggc2VwYXJhdGUgc2V4ZXMKCiQkCk5fZSBcYXBwcm94IFxmcmFjezROX2ZOX219e05fZiArIE5fbX0KJCQKCldlJ3JlIGdvaW5nIHRvIGV4cGxvcmUgaG93IHdlbGwgdGhlICROX2UkIHdlIGV4cGVjdCBiYXNlZCBvbiB0aGUgbnVtYmVyIG9mIHJlcHJvZHVjaW5nIGZlbWFsZXMgYW5kIG1hbGVzIHByZWRpY3RzIHRoZSAkTl9lJCB3ZSBpbmZlciBmcm9tIHRoZSBvYnNlcnZlZCB2YXJpYW5jZSBpbiBhbGxlbGUgZnJlcXVlbmN5IGFtb25nIG9mZnNwcmluZyBwb3B1bGF0aW9ucy4KCiMjIyBUaGUgYmFzaWMgc2ltdWxhdGlvbgoKVG8gZXhwbG9yZSB0aGlzIHJlbGF0aW9uc2hpcCB3ZSB3aWxsIHNpbXVsYXRlIHRoZSBwcm9jZXNzIG9mIG1hdGluZyBhbmQgcmVwcm9kdWN0aW9uIGluIGEgZmluaXRlIHBvcHVsYXRpb24gd2l0aCBzZXBhcmF0ZSBzZXhlcy4gSGVyZSdzIHRoZSBwcm9jZXNzOgoKMS4gV2UgY2FsY3VsYXRlIGdlbm90eXBlIGZyZXF1ZW5jaWVzIGluIHRoZSBiYXNlIHBvcHVsYXRpb24gZ2l2ZW4gdGhlIGFsbGVsZSBmcmVxdWVuY3kgJHBfMCQuCgoyLiBXZSBjb25zdHJ1Y3QgJE5fZiQgZmVtYWxlIGdlbm90eXBlcyBhbmQgJE5fbSQgbWFsZSBnZW5vdHlwZXMgYXQgcmFuZG9tIGdpdmVuIHRob3NlIGdlbm90eXBlIGZyZXF1ZW5jaWVzLgoKMy4gV2UgY29uc3RydWN0IGFsbGVsZSBmcmVxdWVuY2llcyBpbiBhbiBvZmZzcHJpbmcgcG9wdWxhdGlvbiBieSBwaWNraW5nIG9uZSBmZW1hbGUgYW5kIG9uZSBtYWxlIGF0IHJhbmRvbSwgYXZlcmFnaW5nIHRoZSBhbGxlbGUgZnJlcXVlbmN5IHdpdGhpbiB0aGVtIHRvIHByb2R1Y2UgdGhlIGFsbGVsZSBmcmVxdWVuY3kgaW4gdGhlaXIgb2Zmc3ByaW5nLCBhbmQgcmVwZWF0aW5nIHVudGlsIHdlJ3ZlIGNvbnN0cnVjdGVkICROX3twb3B9JCBvZmZzcHJpbmcuIFRoZW4gd2UgdGFrZSB0aGUgYXZlcmFnZSBhbGxlbGUgZnJlcXVlbmN5IGFjcm9zcyBhbGwgb2YgdGhlc2Ugb2Zmc3ByaW5nIGFuZCB3ZSBoYXZlIHRoZSBhbGxlbGUgZnJlcXVlbmN5IGluIHRoZSBvZmZzcHJpbmcgcG9wdWxhdGlvbi4KCjQuIFdlIHRoZW4gcmVwZWF0IHRoYXQgcHJvY2VzcyB0byBwcm9kdWNlICROX3tzYW1wfSQgZGlmZmVyZW50IG9mZnNwcmluZyBwb3B1bGF0aW9ucy4gCgpIZXJlJ3MgYW4gZXhhbXBsZSB1c2luZyAkcF8wID0gMC41JCwgJE5fZiA9IDEwLCAkTl9tID0gMjAkLCAkTl97cG9wfSA9IDEwMCQsIGFuZCAkTl97c2FtcH0gPSAxMDAkLgoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGdncGxvdDIpCgpybShsaXN0ID0gbHMoKSkKCiMjIE5PVEU6IERvbid0IGluY2x1ZGUgdGhpcyBsaW5lIGluIHlvdXIgY29kZS4gSSBpbmNsdWRlIGl0IHNvIHRoYXQgaWYgCiMjIHlvdSBydW4gdGhlIFIgY29kZSB5b3Vyc2VsZiwgeW91J2xsIGdldCB0aGUgc2FtZSByZXN1bHRzIHRoYXQgSSBkby4KIyMgVGhlIHJhbmRvbSBudW1iZXIgZ2VuZXJhdG9yIGluIFIgKGFuZCBpbiBvdGhlciBwcm9ncmFtbWluZyBsYW5ndWFnZXMKIyMgZm9yIHRoYXQgbWF0dGVyKSBpcyBhIHBzZXVkby1yYW5kb20gbnVtYmVyIGdlbmVyYXRvci4gVGhlIG51bWJlcnMgaXQKIyMgcHJvZHVjZXMgYXBwZWFyIHJhbmRvbSwgYnV0IGl0IHlvdSBnaXZlIGl0IHRoZSBzYW1lICJzZWVkIiwgaXQgd2lsbAojIyBnZW5lcmF0ZSB0aGUgc2FtZSBzZXF1ZW5jZSBvZiBudW1iZXJzLgojIwpzZXQuc2VlZCgxMjM0KQoKc2ltdWxhdGUgPC0gZnVuY3Rpb24obl9mLCBuX20sIHBfMCwgbl9zYW1wLCBuX3BvcCkgewogICMjIGdlbm90eXBlIGZyZXF1ZW5pZXMgaW4gSGFyZHktV2VpbmJlcmcKICAjIwogIHggPC0gYyhwXzBeMiwgMipwXzAqKDEtcF8wKSwgKDEtcF8wKV4yKQogIHAgPC0gbnVtZXJpYyhuX3NhbXApCiAgZm9yIChpIGluIDE6bl9zYW1wKSB7CiAgICAjIyBjb25zdHJ1Y3QgZmVtYWxlIGdlbm90eXBlcyBhdCByYW5kb20gZ2l2ZW4gSC1XCiAgICAjIwogICAgZiA8LSBudW1lcmljKG5fZikKICAgIGZvciAoaiBpbiAxOm5fZikgewogICAgICB0bXAgPC0gcm11bHRpbm9tKDEsIDEsIHgpCiAgICAgIGZbal0gPC0gd2hpY2godG1wWywgMV0gPT0gMSkKICAgIH0KICAgICMjIGNvbnN0cnVjdCBtYWxlIGdlbm90eXBlcyBhdCByYW5kb20gZ2l2ZW4gSC1XCiAgICAjIwogICAgbSA8LSBudW1lcmljKG5fbSkKICAgIGZvciAoaiBpbiAxOm5fbSkgewogICAgICB0bXAgPC0gcm11bHRpbm9tKDEsIDEsIHgpCiAgICAgIG1bal0gPC0gd2hpY2godG1wWywgMV0gPT0gMSkKICAgIH0KICAgICMjIGNvbnN0cnVjdCBvZmZzcHJpbmcgYXMgYXZlcmFnZSBvZiBtb20gYW5kIGRhZAogICAgIyMKICAgIHBbaV0gPC0gMAogICAgZm9yIChqIGluIDE6bl9wb3ApIHsKICAgICAgcF9mIDwtIChzYW1wbGUoZiwgMSkgLSAxKS8yCiAgICAgIHBfbSA8LSAoc2FtcGxlKG0sIDEpIC0gMSkvMgogICAgICBwW2ldIDwtIHBbaV0gKyAocF9mICsgcF9tKS8yCiAgICB9CiAgICBwW2ldIDwtIHBbaV0vbl9wb3AKICB9CiAgcmV0dXJuKHApCn0KCnJlc3VsdCA8LSBzaW11bGF0ZShuX2YgPSAxMCwgbl9tID0gMjAsIHBfMCA9IDAuNSwgbl9zYW1wID0gMTAwLCAKICAgICAgICAgICAgICAgICAgIG5fcG9wID0gMTAwKQoKcCA8LSBnZ3Bsb3QodGliYmxlKHAgPSByZXN1bHQpLCBhZXMoeCA9IHApKSArIAogIGdlb21faGlzdG9ncmFtKGJpbndpZHRoID0gMC4wMSkgKyAKICB0aGVtZV9idygpICsKICBnZ3RpdGxlKCJPbmUgcmVwbGljYXRpb24gb2Ygc2FtcGxpbmciKQpwCmBgYAoKQXMgeW91IGNhbiBzZWUsIHRoZXJlJ3MgYSBsb3Qgb2YgdmFyaWF0aW9uIGluIGFsbGVsZSBmcmVxdWVuY3kgYW1vbmcgdGhlIG9mZnNwcmluZyBwb3B1bGF0aW9ucywgYnV0IGJvdGggdGhlIG1lYW4gYW5kIHRoZSB2YXJpYW5jZSBvZiB0aGUgYWxsZWxlIGZyZXF1ZW5jeSBhcmUgcHJldHR5IGNsb3NlIHRvIHdoYXQgd2UnZCBleHBlY3QuCgpgYGB7cn0KcF8wIDwtIDAuNQpuX2UgPC0gNCoxMCoyMC8oMTAgKyAyMCkKCmRhdCA8LSB0aWJibGUoQ2F0ZWdvcnkgPSBjKCJFeHBlY3RlZCIsICJPYnNlcnZlZCIpLAogICAgICAgICAgICAgIE1lYW4gPSBjKHBfMCwgcm91bmQobWVhbihyZXN1bHQpLCAzKSksCiAgICAgICAgICAgICAgVmFyaWFuY2UgPSBjKHJvdW5kKHBfMCooMS1wXzApLygyKm5fZSksIDYpLAogICAgICAgICAgICAgICAgICAgICAgICAgICByb3VuZCh2YXIocmVzdWx0KSwgNikpKQpkYXQKYGBgCgojIyMgRXh0ZW5kaW5nIHRoZSBzaW11bGF0aW9uCgpXaGF0IHdlJ3ZlIGRvbmUgc28gZmFyIGlzIG5pY2UsIGJ1dCB3ZSd2ZSBvbmx5IGxvb2tlZCBhdCB0aGUgcmVzdWx0cyBmcm9tIG9uZSBzaW11bGF0ZWQgYmFzZSBwb3B1bGF0aW9uLiBNYXliZSB3ZSB3ZXJlIGx1Y2t5IGluIG91ciBjaG9pY2Ugb2YgYmFzZSBwb3B1bGF0aW9uLiBMZXQncyBzZWUgd2hhdCB0aGlzIGxvb2tzIGxpa2UgaWYgd2UgbG9vayBhdCB0aGUgbWVhbiBhbmQgdmFyaWFuY2UgZnJvbSAxMDAwIHNpbXVsYXRlZCBiYXNlIHBvcHVsYXRpb25zLiBIZXJlIHdlIGRvbid0IHdhbnQgdG8gbG9vayBhdCAxMDAwIGdyYXBocyBsaWtlIHRoZSBvbmUgYWJvdmUsIHNvIEknbGwga2VlcCB0cmFjayBvZiB0aGUgbWVhbiBhbmQgdmFyaWFuY2Ugb2YgYWxsZWxlIGZyZXF1ZW5jaWVzIGZyb20gZWFjaCBiYXNlIHBvcHVsYXRpb24gYW5kIGRpc3BsYXkgaGlzdG9ncmFtcyBvZiB0aGUgbWVhbnMgYW5kIHZhcmlhbmNlcyBhbmQgc2VlIHdoZXJlIHRoZSB2YWx1ZXMgZmFsbCByZWxhdGl2ZSB0byB0aGUgZXhwZWN0YXRpb24uXltOb3RlOiBJZiB5b3UgcnVuIHRoZSBmb2xsb3dpbmcgYmxvY2sgb2YgY29kZSwgZG9uJ3QgYmUgd29ycmllZCBpZiB5b3UgZG9uJ3Qgc2VlIGFueXRoaW5nIGhhcHBlbmluZy4gSXQgdGFrZXMgYmV0d2VlbiAzLTUgbWludXRlcyB0byBydW4gb24gbXkgTWFjQm9va10KCmBgYHtyfQpwX21lYW4gPC0gbnVtZXJpYygxMDAwKQpwX3ZhciA8LSBudW1lcmljKDEwMDApCmZvciAoaSBpbiAxOjEwMDApIHsKICBwIDwtIHNpbXVsYXRlKG5fZiA9IDEwLCBuX20gPSAyMCwgcF8wID0gMC41LCBuX3NhbXAgPSAxMDAsIAogICAgICAgICAgICAgICAgbl9wb3AgPSAxMDApCiAgcF9tZWFuW2ldIDwtIG1lYW4ocCkKICBwX3ZhcltpXSA8LSB2YXIocCkKfQoKcF9kYXQgPC0gdGliYmxlKHBfbWVhbiA9IHBfbWVhbiwgcF92YXIgPSBwX3ZhcikKcCA8LSBnZ3Bsb3QocF9kYXQsIGFlcyh4ID0gcF9tZWFuKSkgKwogIGdlb21faGlzdG9ncmFtKGJpbnMgPSAyNSkgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAuNSwgbGluZXR5cGUgPSAiZGFzaGVkIiwgCiAgICAgICAgICAgICBjb2xvciA9ICJzYWxtb24iKSArCiAgZ2d0aXRsZSgiRGlzdHJpYnV0aW9uIG9mIG1lYW4gYWxsZWxlIGZyZXF1ZW5jaWVzIikgKwogIHRoZW1lX2J3KCkKcApwIDwtIGdncGxvdChwX2RhdCwgYWVzKHggPSBwX3ZhcikpICsKICBnZW9tX2hpc3RvZ3JhbShiaW5zID0gMjUpICsKICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSBwXzAqKDEtcF8wKS8oMipuX2UpLCBsaW5ldHlwZSA9ICJkYXNoZWQiLAogICAgICAgICAgICAgY29sb3IgPSAic2FsbW9uIikgKwogIGdndGl0bGUoIkRpc3RyaWJ1dGlvbiBvZiB2YXJpYW5jZSBpbiBhbGxlbGUgZnJlcXVlbmNpZXMiKSArCiAgdGhlbWVfYncoKQpwCmBgYAoKVGhlcmUgYXJlIHR3byBzdHJpa2luZyBmZWF0dXJlcyBvZiB0aGVzZSBkaXN0cmlidXRpb25zOgoKMS4gVGhlcmUncyBhIHJlYXNvbmFibGUgYW1vdW50IG9mIHZhcmlhYmlsaXR5IGluIGJvdGggdGhlIG1lYW4gYW5kIHRoZSB2YXJpYW5jZSBvZiBhbGxlbGUgZnJlcXVlbmN5LCBhbmQgdGhlcmUncyBhIGxvdCBtb3JlIHZhcmlhYmlsaXR5IGluIHRoZSB2YXJpYW5jZSB0aGFuIGluIHRoZSBtZWFuLgoKMi4gVGhlIGF2ZXJhZ2UgbWVhbiBhbGxlbGUgZnJlcXVlbmN5IGFjcm9zcyB0aGUgMTAwMCBzaW11bGF0aW9ucyBpcyBwcmV0dHkgY2xvc2UgdG8gdGhlIGV4cGVjdGF0aW9uICh0aGUgdmVydGljYWwgZGFzaGVkIGxpbmUpLCBidXQgdGhlIG1lYW4gdmFyaWFuY2UgaW4gYWxsZWxlIGZyZXF1ZW5jeSBpcyBub3RpY2VhYmx5IGxhcmdlciBpbiB0aGUgc2ltdWxhdGlvbnMgdGhhbiB3ZSBleHBlY3QuCgojIyBFeGVyY2lzZSA1CgpUaGUgbGFiIGV4ZXJjaXNlIHRoaXMgd2VlayBhc2tzIHlvdSB0byBleHBsb3JlIHRoaXMgZGlzY3JlcGFuY3kgYSBiaXQgZnVydGhlci4gVG8gZG8gc28geW91IGNhbiB1c2UgdGhlIGZvbGxvd2luZyBmdW5jdGlvbi5eW05vdGUgdGhhdCB0aGlzIGZ1bmN0aW9uIGNhbGxzIGBzaW11bGF0ZSgpYCwgc28geW91J2xsIG5lZWQgdG8gZXhlY3V0ZSBgc2ltdWxhdGUoKWAgaW4geW91ciBSIG5vdGVib29rIGJlZm9yZSB5b3UgcnVuIHRoaXMgZnVuY3Rpb24uXQoKYGBge3J9CnJ1bl9zaW11bGF0aW9uIDwtIGZ1bmN0aW9uKHBfMCwgbl9zYW1wLCBuX3BvcCkgewogIGRhdCA8LSBkYXRhLmZyYW1lKG5fZiA9IE5BLCBuX20gPSBOQSwgbl9lID0gTkEsIG5fZV9vYnMgPSBOQSwgCiAgICAgICAgICAgICAgICAgICAgdmFyX3AgPSBOQSwgdmFyX3Bfb2JzID0gTkEpCiAgZm9yIChuX2YgaW4gYyg1LCAxMCwgMjUsIDUwLCAxMDApKSB7CiAgICBmb3IgKG5fbSBpbiBjKDUsIDEwLCAyNSwgNTAsIDEwMCkpIHsKICAgICAgcCA8LSBzaW11bGF0ZShuX2YgPSBuX2YsIG5fbSA9IG5fbSwgcF8wID0gMC41LCBuX3NhbXAgPSBuX3NhbXAsIAogICAgICAgICAgICAgICAgICAgIG5fcG9wID0gbl9wb3ApCiAgICAgICMjIGNhbGN1bGF0ZSBOZSBhbmQgVmFyKHApIGZyb20gZm9ybXVsYXMKICAgICAgIyMKICAgICAgbl9lIDwtIDQuMCpuX2Yqbl9tLyhuX2YgKyBuX20pCiAgICAgIHZhcl9wIDwtIHBfMCooMSAtIHBfMCkvKDIqbl9lKQogICAgICAjIyBjYWxjdWxhdGUgb2JzZXJ2ZWQgTmUgZnJvbSBvYnNlcnZlZCB2YXJpYW5jZQogICAgICAjIwogICAgICBuX2Vfb2JzIDwtIDAuNSowLjUvKDIqdmFyKHApKQogICAgICBkYXQgPC0gYWRkX3JvdyhkYXQsIG5fZiA9IG5fZiwgbl9tID0gbl9tLCBuX2UgPSBuX2UsIAogICAgICAgICAgICAgICAgICAgICBuX2Vfb2JzID0gbl9lX29icywgdmFyX3AgPSB2YXJfcCwgCiAgICAgICAgICAgICAgICAgICAgIHZhcl9wX29icyA9IHZhcihwKSkKICAgIH0KICB9CiAgZGF0IDwtIHN1YnNldChkYXQsICFpcy5uYShuX2YpKQoKICBwIDwtIGdncGxvdChkYXQsIGFlcyh4ID0gbl9lLCB5ID0gbl9lX29icykpICsKICAgIGdlb21fcG9pbnQoKSArCiAgICBnZW9tX2FibGluZShzbG9wZSA9IDEsIGludGVyY2VwdCA9IDAsIGxpbmV0eXBlID0gImRhc2hlZCIsCiAgICAgICAgICAgICAgICBjb2xvciA9ICJzYWxtb24iKSArCiAgICB0aGVtZV9idygpCiAgcHJpbnQocCkKICBwIDwtIGdncGxvdChkYXQsIGFlcyh4ID0gdmFyX3AsIHkgPSB2YXJfcF9vYnMpKSArCiAgICBnZW9tX3BvaW50KCkgKwogICAgZ2VvbV9hYmxpbmUoc2xvcGUgPSAxLCBpbnRlcmNlcHQgPSAwLCBsaW5ldHlwZSA9ICJkYXNoZWQiLAogICAgICAgICAgICAgICAgY29sb3IgPSAic2FsbW9uIikgKwogICAgdGhlbWVfYncoKQogIHByaW50KHApCn0KYGBgCgoxLiBFeGVjdXRlIGBydW5fc2ltdWxhdGlvbigpYCB3aXRoICRwXzAgPSAwLjUkLCAkbl97c2FtcH0gPSAxMDAkLCBhbmQgCiRuX29icyA9IDEwMCQuIFRoZSBkYXNoZWQgbGluZSBpcyB0aGUgMToxIGxpbmUuIElmIHRoZSBzaW11bGF0aW9ucyBtYXRjaGVkIGV4cGVjdGF0aW9ucywgYWxsIG9mIHRoZSBwb2ludHMgd291bGQgbGllIG9uIHRoZSBsaW5lICh3aXRoIGEgYml0IG9mIHN0YXRpc3RpY2FsIHNhbXBsaW5nIGVycm9yKS4gV2l0aCB0aGVzZSBwYXJhbWV0ZXJzLCB5b3UnbGwgc2VlIHRoYXQgdGhlcmUncyBhIGxhcmdlIG1pc21hdGNoIGJldHdlZW4gdGhlIG9ic2VydmVkIGFuZCBleHBlY3RlZCAkTl9lJC4gVGhlcmUncyBhIHNtYWxsZXIgbWlzbWF0Y2ggYmV0d2VlbiB0aGUgb2JzZXJ2ZWQgYW5kIGV4cGVjdGVkIHZhcmlhbmNlLCBidXQgaXQncyBjb25zaXN0ZW50bHkgaW4gb25lIGRpcmVjdGlvbi4KCjIuIFRyeSBzb21lIGRpZmZlcmVudCBjb21iaW5hdGlvbnMgb2YgJHBfMCQsICRuX3tzYW1wfSQsIGFuZCAkbl97b2JzfSQgYW5kIHJlcG9ydCB5b3VyIHJlc3VsdHMuIEkgYW0gbG9va2luZyBmb3IgYSBtaW5pbXVtIG9mIDEwIGNvbWJpbmF0aW9ucywgYnV0IGlmIHlvdSdkIGxpa2UgdG8gdHJ5IG1vcmUsIGZlZWwgZnJlZS4gRWl0aGVyIHdheSBpbmNsdWRlIHRoZSBwbG90cyB0aGF0IHlvdXIgc2ltdWxhdGlvbnMgcHJvZHVjZS4gVHJlYXQgeW91ciBleHBsb3JhdGlvbiBvZiBkaWZmZXJlbnQgY29tYmluYXRpb25zIG9mIHBhcmFtZXRlcnMgbGlrZSBhbiBleHBlcmltZW50LCBpLmUuLCBvbmx5IGNoYW5nZSBvbmUgYXQgYSB0aW1lIChhdCBsZWFzdCBhdCBmaXJzdCkgc28gdGhhdCB5b3UgY2FuIGlzb2xhdGUgdGhlIGVmZmVjdCBvZiBjaGFuZ2VzIGluIHRoYXQgcGFyYW1ldGVyLiBJZiB5b3UncmUgZmVlbGluZyBhbWJpdGlvdXMgeW91IGNhbiB0cnkgdG8gbG9vayBhdCBjb21iaW5hdGlvbnMsIGJ1dCB0aGF0IGNvdWxkIGJlY29tZSBjaGFsbGVuZ2luZyBwcmV0dHkgcXVpY2tseS4gU2VlIGlmIHlvdSBjYW4gZmluZCBhIGNvbWJpbmF0aW9uIG9mIHBhcmFtZXRlcnMgd2hlcmUgdGhlIGZpdCBvZiBvYnNlcnZlZCB0byBleHBlY3RlZCBpcyByZWFzb25hYmx5IGdvb2QuCgozLiBTdW1tYXJpemUgYW55IG92ZXJhbGwgcGF0dGVybiB5b3Ugc2VlLiBJbiBwYXJ0aWN1bGFyLCBob3cgZG8gZGlmZmVyZW50IGNob2ljZXMgZm9yIHRoZSBwYXJhbWV0ZXJzIGFmZmVjdCB0aGUgZXh0ZW50IG9mIHRoZSBkaXNjcmVwYW5jeSB5b3Ugc2VlLCBpZiBhdCBhbGw/CgpXZSdsbCBkaXNjdXNzIGFueSBnZW5lcmFsIHBhdHRlcm5zIHRoYXQgZW1lcmdlIGFuZCBleHBsb3JlIGV4cGxhbmF0aW9ucyBmb3IgdGhlbSBpbiBsZWN0dXJlIG9uIDUgT2N0b2Jlci4=