Some preliminary notes
A couple of important notes about this week’s lab exercise:
The next several weeks will involve further exploration of ideas
that we cover in class rather than analysis of data sets.
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 Nick the explorations I ask for, you’ll
get full credit for this exercise, i.e., 10 out of 10 possible points.
The exercises for the next couple of weeks will be similar, but I’ll ask
you to stretch yourself a bit and explain the patterns you see.
The idea we want to explore
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:
We calculate genotype frequencies in the base population given
the allele frequency \(p_0\).
We construct \(N_f\) female
genotypes and \(N_m\) male genotypes at
random given those genotype frequencies.
We construct allele frequencies in an offspring population by
picking one female and one male at random, using Mendel’s rules to
construct one 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.
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\).
## The options line is optional. (Do you like the pun?) It merely turns
## off the verbose messages we'd otherwise get when loading the
## tidyverse
##
options(tidyverse.quiet = TRUE)
library(tidyverse)
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 if you give it the same "seed", it will
## generate the same sequence of numbers.
##
set.seed(1234)
## Homozygous A1A1 if x == 1
## Heterozygous A1A2 if x == 2
## Homozygous A2A2 if x == 3
##
## Return value is the number of A1 alleles in the gamete
##
mendel <- function(x) {
if (x == 1) {
retval <- 1
} else if (x == 3) {
retval <- 0
} else {
if (runif(1) < 0.5) {
retval <- 0
} else {
retval <- 1
}
}
return(retval)
}
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 <- mendel(sample(f, 1))
p_m <- mendel(sample(m, 1))
p[i] <- p[i] + p_f + p_m
}
p[i] <- p[i]/(2*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 the table below shows that 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.
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:
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.
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.
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)
}
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.
Try some different combinations of \(p_0\), \(n_{samp}\), and \(n_{pop}\) 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.
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 3 October.
LS0tCnRpdGxlOiAiRXhwbG9yaW5nIHRoZSBwcm9wZXJ0aWVzIG9mIGVmZmVjdGl2ZSBwb3B1bGF0aW9uIHNpemUiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMjIFNvbWUgcHJlbGltaW5hcnkgbm90ZXMKCkEgY291cGxlIG9mIGltcG9ydGFudCBub3RlcyBhYm91dCB0aGlzIHdlZWsncyBsYWIgZXhlcmNpc2U6CgoxLiBUaGUgbmV4dCBzZXZlcmFsIHdlZWtzXltVbnRpbCB3ZSBnZXQgdG8gdGhlIDE3dGggb2YgT2N0b2Jlcl0gd2lsbCBpbnZvbHZlIGZ1cnRoZXIgZXhwbG9yYXRpb24gb2YgaWRlYXMgdGhhdCB3ZSBjb3ZlciBpbiBjbGFzcyByYXRoZXIgdGhhbiBhbmFseXNpcyBvZiBkYXRhIHNldHMuCgoyLiBUaGlzIHdlZWsncyBleGVyY2lzZSB3b24ndCByZXF1aXJlIHlvdSB0byBhbnN3ZXIgYW55IHF1ZXN0aW9ucy4gSXQgd2lsbCBzaW1wbHkgYXNrIHlvdSB0byBleHBsb3JlICh0aHJvdWdoIHNpbXVsYXRpb25zKSBzb21lIHByb3BlcnRpZXMgb2YgZ2VuZXRpYyBkcmlmdC4gSWYgeW91IHNob3cgTmljayB0aGUgZXhwbG9yYXRpb25zIEkgYXNrIGZvciwgeW91J2xsIGdldCBmdWxsIGNyZWRpdCBmb3IgdGhpcyBleGVyY2lzZSwgaS5lLiwgMTAgb3V0IG9mIDEwIHBvc3NpYmxlIHBvaW50cy4gVGhlIGV4ZXJjaXNlcyBmb3IgdGhlIG5leHQgY291cGxlIG9mIHdlZWtzIHdpbGwgYmUgc2ltaWxhciwgYnV0IEknbGwgYXNrIHlvdSB0byBzdHJldGNoIHlvdXJzZWxmIGEgYml0IGFuZCBleHBsYWluIHRoZSBwYXR0ZXJucyB5b3Ugc2VlLgoKIyMgVGhlIGlkZWEgd2Ugd2FudCB0byBleHBsb3JlXltNb3JlIHByZWNpc2VseSwgdGhlIGlkZWEgSSdtICphc3NpZ25pbmcqIHlvdSB0byBleHBsb3JlLl0KCkluIGxlY3R1cmUgd2UnbGwgc2VlIHRoYXQgdGhlIHZhcmlhbmNlIGVmZmVjdGl2ZSBzaXplIG9mIGEgcG9wdWxhdGlvbiBpcyBkZWZpbmVkIGFzIHRoZSBzaXplIG9mIGFuIGlkZWFsIHBvcHVsYXRpb24gdGhhdCBoYXMgdGhlIHNhbWUgdmFyaWFuY2UgaW4gYWxsZWxlIGZyZXF1ZW5jeSBhbW9uZyBvZmZzcHJpbmcgcG9wdWxhdGlvbnMgYXMgdGhlIGFjdHVhbCBwb3B1bGF0aW9uIHdlJ3JlIGludGVyZXN0ZWQgaW4sIGkuZS4sCgokJApOX2Veeyh2KX0gPSBcZnJhY3twKDEtcCl9ezJcd2lkZWhhdHtWYXIocCl9fQokJAoKV2Ugd2lsbCBhbHNvIHNlZSB0aGF0IGZvciBhIHBvcHVsYXRpb24gd2l0aCBzZXBhcmF0ZSBzZXhlcwoKJCQKTl9lIFxhcHByb3ggXGZyYWN7NE5fZk5fbX17Tl9mICsgTl9tfQokJAoKV2UncmUgZ29pbmcgdG8gZXhwbG9yZSBob3cgd2VsbCB0aGUgJE5fZSQgd2UgZXhwZWN0IGJhc2VkIG9uIHRoZSBudW1iZXIgb2YgcmVwcm9kdWNpbmcgZmVtYWxlcyBhbmQgbWFsZXMgcHJlZGljdHMgdGhlICROX2UkIHdlIGluZmVyIGZyb20gdGhlIG9ic2VydmVkIHZhcmlhbmNlIGluIGFsbGVsZSBmcmVxdWVuY3kgYW1vbmcgb2Zmc3ByaW5nIHBvcHVsYXRpb25zLgoKIyMjIFRoZSBiYXNpYyBzaW11bGF0aW9uCgpUbyBleHBsb3JlIHRoaXMgcmVsYXRpb25zaGlwIHdlIHdpbGwgc2ltdWxhdGUgdGhlIHByb2Nlc3Mgb2YgbWF0aW5nIGFuZCByZXByb2R1Y3Rpb24gaW4gYSBmaW5pdGUgcG9wdWxhdGlvbiB3aXRoIHNlcGFyYXRlIHNleGVzLiBIZXJlJ3MgdGhlIHByb2Nlc3M6CgoxLiBXZSBjYWxjdWxhdGUgZ2Vub3R5cGUgZnJlcXVlbmNpZXMgaW4gdGhlIGJhc2UgcG9wdWxhdGlvbiBnaXZlbiB0aGUgYWxsZWxlIGZyZXF1ZW5jeSAkcF8wJC5eW0FuZCBhc3N1bWluZyB0aGF0IHRoZSBnZW5vdHlwZXMgYXJlIGluIEhhcmR5LVdlaW5iZXJnIHByb3BvcnRpb25zLl0KCjIuIFdlIGNvbnN0cnVjdCAkTl9mJCBmZW1hbGUgZ2Vub3R5cGVzIGFuZCAkTl9tJCBtYWxlIGdlbm90eXBlcyBhdCByYW5kb20gZ2l2ZW4gdGhvc2UgZ2Vub3R5cGUgZnJlcXVlbmNpZXMuXltBc3N1bWluZyB0aGF0IHRoZSBnZW5vdHlwZSBmcmVxdWVuY2llcyBhcmUgdGhlIHNhbWUgaW4gYm90aCBzZXhlcy5dCgozLiBXZSBjb25zdHJ1Y3QgYWxsZWxlIGZyZXF1ZW5jaWVzIGluIGFuIG9mZnNwcmluZyBwb3B1bGF0aW9uIGJ5IHBpY2tpbmcgb25lIGZlbWFsZSBhbmQgb25lIG1hbGUgYXQgcmFuZG9tLCB1c2luZyBNZW5kZWwncyBydWxlcyB0byBjb25zdHJ1Y3Qgb25lIG9mZnNwcmluZywgYW5kIHJlcGVhdGluZyB1bnRpbCB3ZSd2ZSBjb25zdHJ1Y3RlZCAkTl97cG9wfSQgb2Zmc3ByaW5nLiBUaGVuIHdlIHRha2UgdGhlIGF2ZXJhZ2UgYWxsZWxlIGZyZXF1ZW5jeSBhY3Jvc3MgYWxsIG9mIHRoZXNlIG9mZnNwcmluZyBhbmQgd2UgaGF2ZSB0aGUgYWxsZWxlIGZyZXF1ZW5jeSBpbiB0aGUgb2Zmc3ByaW5nIHBvcHVsYXRpb24uCgo0LiBXZSB0aGVuIHJlcGVhdCB0aGF0IHByb2Nlc3MgdG8gcHJvZHVjZSAkTl97c2FtcH0kIGRpZmZlcmVudCBvZmZzcHJpbmcgcG9wdWxhdGlvbnMuIAoKSGVyZSdzIGFuIGV4YW1wbGUgdXNpbmcgJHBfMCA9IDAuNSQsICROX2YgPSAxMCQsICROX20gPSAyMCQsICROX3twb3B9ID0gMTAwJCwgYW5kICROX3tzYW1wfSA9IDEwMCQuCgpgYGB7cn0KIyMgVGhlIG9wdGlvbnMgbGluZSBpcyBvcHRpb25hbC4gKERvIHlvdSBsaWtlIHRoZSBwdW4/KSBJdCBtZXJlbHkgdHVybnMKIyMgb2ZmIHRoZSB2ZXJib3NlIG1lc3NhZ2VzIHdlJ2Qgb3RoZXJ3aXNlIGdldCB3aGVuIGxvYWRpbmcgdGhlIAojIyB0aWR5dmVyc2UKIyMKb3B0aW9ucyh0aWR5dmVyc2UucXVpZXQgPSBUUlVFKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShnZ3Bsb3QyKQoKcm0obGlzdCA9IGxzKCkpCgojIyBOT1RFOiBEb24ndCBpbmNsdWRlIHRoaXMgbGluZSBpbiB5b3VyIGNvZGUuIEkgaW5jbHVkZSBpdCBzbyB0aGF0IGlmIAojIyB5b3UgcnVuIHRoZSBSIGNvZGUgeW91cnNlbGYsIHlvdSdsbCBnZXQgdGhlIHNhbWUgcmVzdWx0cyB0aGF0IEkgZG8uCiMjIFRoZSByYW5kb20gbnVtYmVyIGdlbmVyYXRvciBpbiBSIChhbmQgaW4gb3RoZXIgcHJvZ3JhbW1pbmcgbGFuZ3VhZ2VzCiMjIGZvciB0aGF0IG1hdHRlcikgaXMgYSBwc2V1ZG8tcmFuZG9tIG51bWJlciBnZW5lcmF0b3IuIFRoZSBudW1iZXJzIGl0CiMjIHByb2R1Y2VzIGFwcGVhciByYW5kb20sIGJ1dCBpZiB5b3UgZ2l2ZSBpdCB0aGUgc2FtZSAic2VlZCIsIGl0IHdpbGwKIyMgZ2VuZXJhdGUgdGhlIHNhbWUgc2VxdWVuY2Ugb2YgbnVtYmVycy4KIyMKc2V0LnNlZWQoMTIzNCkKCiMjIEhvbW96eWdvdXMgQTFBMSBpZiB4ID09IDEKIyMgSGV0ZXJvenlnb3VzIEExQTIgaWYgeCA9PSAyCiMjIEhvbW96eWdvdXMgQTJBMiBpZiB4ID09IDMKIyMKIyMgUmV0dXJuIHZhbHVlIGlzIHRoZSBudW1iZXIgb2YgQTEgYWxsZWxlcyBpbiB0aGUgZ2FtZXRlCiMjCm1lbmRlbCA8LSBmdW5jdGlvbih4KSB7CiAgaWYgKHggPT0gMSkgewogICAgcmV0dmFsIDwtIDEKICB9IGVsc2UgaWYgKHggPT0gMykgewogICAgcmV0dmFsIDwtIDAKICB9IGVsc2UgewogICAgaWYgKHJ1bmlmKDEpIDwgMC41KSB7CiAgICAgIHJldHZhbCA8LSAwCiAgICB9IGVsc2UgewogICAgICByZXR2YWwgPC0gMQogICAgfQogIH0KICByZXR1cm4ocmV0dmFsKQp9CgpzaW11bGF0ZSA8LSBmdW5jdGlvbihuX2YsIG5fbSwgcF8wLCBuX3NhbXAsIG5fcG9wKSB7CiAgIyMgZ2Vub3R5cGUgZnJlcXVlbmllcyBpbiBIYXJkeS1XZWluYmVyZwogICMjCiAgeCA8LSBjKHBfMF4yLCAyKnBfMCooMS1wXzApLCAoMS1wXzApXjIpCiAgcCA8LSBudW1lcmljKG5fc2FtcCkKICBmb3IgKGkgaW4gMTpuX3NhbXApIHsKICAgICMjIGNvbnN0cnVjdCBmZW1hbGUgZ2Vub3R5cGVzIGF0IHJhbmRvbSBnaXZlbiBILVcKICAgICMjCiAgICBmIDwtIG51bWVyaWMobl9mKQogICAgZm9yIChqIGluIDE6bl9mKSB7CiAgICAgIHRtcCA8LSBybXVsdGlub20oMSwgMSwgeCkKICAgICAgZltqXSA8LSB3aGljaCh0bXBbLCAxXSA9PSAxKQogICAgfQogICAgIyMgY29uc3RydWN0IG1hbGUgZ2Vub3R5cGVzIGF0IHJhbmRvbSBnaXZlbiBILVcKICAgICMjCiAgICBtIDwtIG51bWVyaWMobl9tKQogICAgZm9yIChqIGluIDE6bl9tKSB7CiAgICAgIHRtcCA8LSBybXVsdGlub20oMSwgMSwgeCkKICAgICAgbVtqXSA8LSB3aGljaCh0bXBbLCAxXSA9PSAxKQogICAgfQogICAgIyMgY29uc3RydWN0IG9mZnNwcmluZyBhcyBhdmVyYWdlIG9mIG1vbSBhbmQgZGFkCiAgICAjIwogICAgcFtpXSA8LSAwCiAgICBmb3IgKGogaW4gMTpuX3BvcCkgewogICAgICBwX2YgPC0gbWVuZGVsKHNhbXBsZShmLCAxKSkKICAgICAgcF9tIDwtIG1lbmRlbChzYW1wbGUobSwgMSkpCiAgICAgIHBbaV0gPC0gcFtpXSArIHBfZiArIHBfbQogICAgfQogICAgcFtpXSA8LSBwW2ldLygyKm5fcG9wKQogIH0KICByZXR1cm4ocCkKfQoKcmVzdWx0IDwtIHNpbXVsYXRlKG5fZiA9IDEwLCBuX20gPSAyMCwgcF8wID0gMC41LCBuX3NhbXAgPSAxMDAsIAogICAgICAgICAgICAgICAgICAgbl9wb3AgPSAxMDApCgpwIDwtIGdncGxvdCh0aWJibGUocCA9IHJlc3VsdCksIGFlcyh4ID0gcCkpICsgCiAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSAwLjAxKSArIAogIHRoZW1lX2J3KCkgKwogIGdndGl0bGUoIk9uZSByZXBsaWNhdGlvbiBvZiBzYW1wbGluZyIpCnAKYGBgCgpBcyB5b3UgY2FuIHNlZSwgdGhlcmUncyBhIGxvdCBvZiB2YXJpYXRpb24gaW4gYWxsZWxlIGZyZXF1ZW5jeSBhbW9uZyB0aGUgb2Zmc3ByaW5nIHBvcHVsYXRpb25zLCBidXQgdGhlIHRhYmxlIGJlbG93IHNob3dzIHRoYXQgYm90aCB0aGUgbWVhbiBhbmQgdGhlIHZhcmlhbmNlIG9mIHRoZSBhbGxlbGUgZnJlcXVlbmN5IGFyZSBwcmV0dHkgY2xvc2UgdG8gd2hhdCB3ZSdkIGV4cGVjdC4KCmBgYHtyfQpwXzAgPC0gMC41Cm5fZSA8LSA0KjEwKjIwLygxMCArIDIwKQoKZGF0IDwtIHRpYmJsZShDYXRlZ29yeSA9IGMoIkV4cGVjdGVkIiwgIk9ic2VydmVkIiksCiAgICAgICAgICAgICAgTWVhbiA9IGMocF8wLCByb3VuZChtZWFuKHJlc3VsdCksIDMpKSwKICAgICAgICAgICAgICBWYXJpYW5jZSA9IGMocm91bmQocF8wKigxLXBfMCkvKDIqbl9lKSwgNiksCiAgICAgICAgICAgICAgICAgICAgICAgICAgIHJvdW5kKHZhcihyZXN1bHQpLCA2KSkpCmRhdApgYGAKCiMjIyBFeHRlbmRpbmcgdGhlIHNpbXVsYXRpb24KCldoYXQgd2UndmUgZG9uZSBzbyBmYXIgaXMgbmljZSwgYnV0IHdlJ3ZlIG9ubHkgbG9va2VkIGF0IHRoZSByZXN1bHRzIGZyb20gb25lIHNpbXVsYXRlZCBiYXNlIHBvcHVsYXRpb24uIE1heWJlIHdlIHdlcmUgbHVja3kgaW4gb3VyIGNob2ljZSBvZiBiYXNlIHBvcHVsYXRpb24uIExldCdzIHNlZSB3aGF0IHRoaXMgbG9va3MgbGlrZSBpZiB3ZSBsb29rIGF0IHRoZSBtZWFuIGFuZCB2YXJpYW5jZSBmcm9tIDEwMDAgc2ltdWxhdGVkIGJhc2UgcG9wdWxhdGlvbnMuIEhlcmUgd2UgZG9uJ3Qgd2FudCB0byBsb29rIGF0IDEwMDAgZ3JhcGhzIGxpa2UgdGhlIG9uZSBhYm92ZSwgc28gSSdsbCBrZWVwIHRyYWNrIG9mIHRoZSBtZWFuIGFuZCB2YXJpYW5jZSBvZiBhbGxlbGUgZnJlcXVlbmNpZXMgZnJvbSBlYWNoIGJhc2UgcG9wdWxhdGlvbiBhbmQgZGlzcGxheSBoaXN0b2dyYW1zIG9mIHRoZSBtZWFucyBhbmQgdmFyaWFuY2VzIGFuZCBzZWUgd2hlcmUgdGhlIHZhbHVlcyBmYWxsIHJlbGF0aXZlIHRvIHRoZSBleHBlY3RhdGlvbi5eW05vdGU6IElmIHlvdSBydW4gdGhlIGZvbGxvd2luZyBibG9jayBvZiBjb2RlLCBkb24ndCBiZSB3b3JyaWVkIGlmIHlvdSBkb24ndCBzZWUgYW55dGhpbmcgaGFwcGVuaW5nLiBJdCB0YWtlcyAzLTUgbWludXRlcyB0byBydW4gb24gbXkgTWFjQm9va10KCmBgYHtyfQpwX21lYW4gPC0gbnVtZXJpYygxMDAwKQpwX3ZhciA8LSBudW1lcmljKDEwMDApCmZvciAoaSBpbiAxOjEwMDApIHsKICBwIDwtIHNpbXVsYXRlKG5fZiA9IDEwLCBuX20gPSAyMCwgcF8wID0gMC41LCBuX3NhbXAgPSAxMDAsIAogICAgICAgICAgICAgICAgbl9wb3AgPSAxMDApCiAgcF9tZWFuW2ldIDwtIG1lYW4ocCkKICBwX3ZhcltpXSA8LSB2YXIocCkKfQoKcF9kYXQgPC0gdGliYmxlKHBfbWVhbiA9IHBfbWVhbiwgcF92YXIgPSBwX3ZhcikKcCA8LSBnZ3Bsb3QocF9kYXQsIGFlcyh4ID0gcF9tZWFuKSkgKwogIGdlb21faGlzdG9ncmFtKGJpbnMgPSAyNSkgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAuNSwgbGluZXR5cGUgPSAiZGFzaGVkIiwgCiAgICAgICAgICAgICBjb2xvciA9ICJzYWxtb24iKSArCiAgZ2d0aXRsZSgiRGlzdHJpYnV0aW9uIG9mIG1lYW4gYWxsZWxlIGZyZXF1ZW5jaWVzIikgKwogIHRoZW1lX2J3KCkKcApwIDwtIGdncGxvdChwX2RhdCwgYWVzKHggPSBwX3ZhcikpICsKICBnZW9tX2hpc3RvZ3JhbShiaW5zID0gMjUpICsKICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSBwXzAqKDEtcF8wKS8oMipuX2UpLCBsaW5ldHlwZSA9ICJkYXNoZWQiLAogICAgICAgICAgICAgY29sb3IgPSAic2FsbW9uIikgKwogIGdndGl0bGUoIkRpc3RyaWJ1dGlvbiBvZiB2YXJpYW5jZSBpbiBhbGxlbGUgZnJlcXVlbmNpZXMiKSArCiAgdGhlbWVfYncoKQpwCmBgYAoKVGhlcmUgYXJlIHR3byBzdHJpa2luZyBmZWF0dXJlcyBvZiB0aGVzZSBkaXN0cmlidXRpb25zOgoKMS4gVGhlcmUncyBhIHJlYXNvbmFibGUgYW1vdW50IG9mIHZhcmlhYmlsaXR5IGluIGJvdGggdGhlIG1lYW4gYW5kIHRoZSB2YXJpYW5jZSBvZiBhbGxlbGUgZnJlcXVlbmN5LCBhbmQgdGhlcmUncyBhIGxvdCBtb3JlIHZhcmlhYmlsaXR5IGluIHRoZSB2YXJpYW5jZSB0aGFuIGluIHRoZSBtZWFuLgoKMi4gVGhlIGF2ZXJhZ2UgbWVhbiBhbGxlbGUgZnJlcXVlbmN5IGFjcm9zcyB0aGUgMTAwMCBzaW11bGF0aW9ucyBpcyBwcmV0dHkgY2xvc2UgdG8gdGhlIGV4cGVjdGF0aW9uICh0aGUgdmVydGljYWwgZGFzaGVkIGxpbmUpLCBidXQgdGhlIG1lYW4gdmFyaWFuY2UgaW4gYWxsZWxlIGZyZXF1ZW5jeSBpcyBub3RpY2VhYmx5IGxhcmdlciBpbiB0aGUgc2ltdWxhdGlvbnMgdGhhbiB3ZSBleHBlY3QuCgojIyBFeGVyY2lzZSA1CgpUaGUgbGFiIGV4ZXJjaXNlIHRoaXMgd2VlayBhc2tzIHlvdSB0byBleHBsb3JlIHRoaXMgZGlzY3JlcGFuY3kgYSBiaXQgZnVydGhlci4gVG8gZG8gc28geW91IGNhbiB1c2UgdGhlIGZvbGxvd2luZyBmdW5jdGlvbi5eW05vdGUgdGhhdCB0aGlzIGZ1bmN0aW9uIGNhbGxzIGJvdGggYG1lbmRlbCgpYCBhbmQgYHNpbXVsYXRlKClgLCBzbyB5b3UnbGwgbmVlZCB0byBleGVjdXRlIGJvdGggaW4geW91ciBSIG5vdGVib29rIGJlZm9yZSB5b3UgcnVuIHRoaXMgZnVuY3Rpb24uXQoKYGBge3J9CnJ1bl9zaW11bGF0aW9uIDwtIGZ1bmN0aW9uKHBfMCwgbl9zYW1wLCBuX3BvcCkgewogIGRhdCA8LSBkYXRhLmZyYW1lKG5fZiA9IE5BLCBuX20gPSBOQSwgbl9lID0gTkEsIG5fZV9vYnMgPSBOQSwgCiAgICAgICAgICAgICAgICAgICAgdmFyX3AgPSBOQSwgdmFyX3Bfb2JzID0gTkEpCiAgZm9yIChuX2YgaW4gYyg1LCAxMCwgMjUsIDUwLCAxMDApKSB7CiAgICBmb3IgKG5fbSBpbiBjKDUsIDEwLCAyNSwgNTAsIDEwMCkpIHsKICAgICAgcCA8LSBzaW11bGF0ZShuX2YgPSBuX2YsIG5fbSA9IG5fbSwgcF8wID0gMC41LCBuX3NhbXAgPSBuX3NhbXAsIAogICAgICAgICAgICAgICAgICAgIG5fcG9wID0gbl9wb3ApCiAgICAgICMjIGNhbGN1bGF0ZSBOZSBhbmQgVmFyKHApIGZyb20gZm9ybXVsYXMKICAgICAgIyMKICAgICAgbl9lIDwtIDQuMCpuX2Yqbl9tLyhuX2YgKyBuX20pCiAgICAgIHZhcl9wIDwtIHBfMCooMSAtIHBfMCkvKDIqbl9lKQogICAgICAjIyBjYWxjdWxhdGUgb2JzZXJ2ZWQgTmUgZnJvbSBvYnNlcnZlZCB2YXJpYW5jZQogICAgICAjIwogICAgICBuX2Vfb2JzIDwtIDAuNSowLjUvKDIqdmFyKHApKQogICAgICBkYXQgPC0gYWRkX3JvdyhkYXQsIG5fZiA9IG5fZiwgbl9tID0gbl9tLCBuX2UgPSBuX2UsIAogICAgICAgICAgICAgICAgICAgICBuX2Vfb2JzID0gbl9lX29icywgdmFyX3AgPSB2YXJfcCwgCiAgICAgICAgICAgICAgICAgICAgIHZhcl9wX29icyA9IHZhcihwKSkKICAgIH0KICB9CiAgZGF0IDwtIHN1YnNldChkYXQsICFpcy5uYShuX2YpKQoKICBwIDwtIGdncGxvdChkYXQsIGFlcyh4ID0gbl9lLCB5ID0gbl9lX29icykpICsKICAgIGdlb21fcG9pbnQoKSArCiAgICBnZW9tX2FibGluZShzbG9wZSA9IDEsIGludGVyY2VwdCA9IDAsIGxpbmV0eXBlID0gImRhc2hlZCIsCiAgICAgICAgICAgICAgICBjb2xvciA9ICJzYWxtb24iKSArCiAgICB0aGVtZV9idygpCiAgcHJpbnQocCkKICBwIDwtIGdncGxvdChkYXQsIGFlcyh4ID0gdmFyX3AsIHkgPSB2YXJfcF9vYnMpKSArCiAgICBnZW9tX3BvaW50KCkgKwogICAgZ2VvbV9hYmxpbmUoc2xvcGUgPSAxLCBpbnRlcmNlcHQgPSAwLCBsaW5ldHlwZSA9ICJkYXNoZWQiLAogICAgICAgICAgICAgICAgY29sb3IgPSAic2FsbW9uIikgKwogICAgdGhlbWVfYncoKQogIHByaW50KHApCn0KYGBgCgoxLiBFeGVjdXRlIGBydW5fc2ltdWxhdGlvbigpYCB3aXRoICRwXzAgPSAwLjUkLCAkbl97c2FtcH0gPSAxMDAkLCBhbmQgCiRuX3tvYnN9ID0gMTAwJC4gVGhlIGRhc2hlZCBsaW5lIGlzIHRoZSAxOjEgbGluZS4gSWYgdGhlIHNpbXVsYXRpb25zIG1hdGNoZWQgZXhwZWN0YXRpb25zLCBhbGwgb2YgdGhlIHBvaW50cyB3b3VsZCBsaWUgb24gdGhlIGxpbmUgKHdpdGggYSBiaXQgb2Ygc3RhdGlzdGljYWwgc2FtcGxpbmcgZXJyb3IpLiBXaXRoIHRoZXNlIHBhcmFtZXRlcnMsIHlvdSdsbCBzZWUgdGhhdCB0aGVyZSdzIGEgbGFyZ2UgbWlzbWF0Y2ggYmV0d2VlbiB0aGUgb2JzZXJ2ZWQgYW5kIGV4cGVjdGVkICROX2UkLiBUaGVyZSdzIGEgc21hbGxlciBtaXNtYXRjaCBiZXR3ZWVuIHRoZSBvYnNlcnZlZCBhbmQgZXhwZWN0ZWQgdmFyaWFuY2UsIGJ1dCBpdCdzIGNvbnNpc3RlbnRseSBpbiBvbmUgZGlyZWN0aW9uLgoKMi4gVHJ5IHNvbWUgZGlmZmVyZW50IGNvbWJpbmF0aW9ucyBvZiAkcF8wJCwgJG5fe3NhbXB9JCwgYW5kICRuX3twb3B9JCBhbmQgcmVwb3J0IHlvdXIgcmVzdWx0cy4gSSBhbSBsb29raW5nIGZvciBhIG1pbmltdW0gb2YgMTAgY29tYmluYXRpb25zLCBidXQgaWYgeW91J2QgbGlrZSB0byB0cnkgbW9yZSwgZmVlbCBmcmVlLiBFaXRoZXIgd2F5IGluY2x1ZGUgdGhlIHBsb3RzIHRoYXQgeW91ciBzaW11bGF0aW9ucyBwcm9kdWNlLiBUcmVhdCB5b3VyIGV4cGxvcmF0aW9uIG9mIGRpZmZlcmVudCBjb21iaW5hdGlvbnMgb2YgcGFyYW1ldGVycyBsaWtlIGFuIGV4cGVyaW1lbnQsIGkuZS4sIG9ubHkgY2hhbmdlIG9uZSBhdCBhIHRpbWUgKGF0IGxlYXN0IGF0IGZpcnN0KSBzbyB0aGF0IHlvdSBjYW4gaXNvbGF0ZSB0aGUgZWZmZWN0IG9mIGNoYW5nZXMgaW4gdGhhdCBwYXJhbWV0ZXIuIElmIHlvdSdyZSBmZWVsaW5nIGFtYml0aW91cyB5b3UgY2FuIHRyeSB0byBsb29rIGF0IGNvbWJpbmF0aW9ucywgYnV0IHRoYXQgY291bGQgYmVjb21lIGNoYWxsZW5naW5nIHByZXR0eSBxdWlja2x5LiBTZWUgaWYgeW91IGNhbiBmaW5kIGEgY29tYmluYXRpb24gb2YgcGFyYW1ldGVycyB3aGVyZSB0aGUgZml0IG9mIG9ic2VydmVkIHRvIGV4cGVjdGVkIGlzIHJlYXNvbmFibHkgZ29vZC4KCjMuIFN1bW1hcml6ZSBhbnkgb3ZlcmFsbCBwYXR0ZXJuIHlvdSBzZWUuIEluIHBhcnRpY3VsYXIsIGhvdyBkbyBkaWZmZXJlbnQgY2hvaWNlcyBmb3IgdGhlIHBhcmFtZXRlcnMgYWZmZWN0IHRoZSBleHRlbnQgb2YgdGhlIGRpc2NyZXBhbmN5IHlvdSBzZWUsIGlmIGF0IGFsbD8KCldlJ2xsIGRpc2N1c3MgYW55IGdlbmVyYWwgcGF0dGVybnMgdGhhdCBlbWVyZ2UgYW5kIGV4cGxvcmUgZXhwbGFuYXRpb25zIGZvciB0aGVtIGluIGxlY3R1cmUgb24gMyBPY3RvYmVyLgo=