From 833a3408c9985fe227ef114f601dd952653c8cb6 Mon Sep 17 00:00:00 2001
From: Jennifer Blanc <jgblanc@uchicago.edu>
Date: Wed, 20 Mar 2019 15:19:31 -0500
Subject: [PATCH 1/3] Add Coalescent Vignette

---
 analysis/The Ancestral Process.Rmd  | 195 ++++++++++++++++++++++++++++
 analysis/The Coalescent Process.Rmd | 193 +++++++++++++++++++++++++++
 2 files changed, 388 insertions(+)
 create mode 100644 analysis/The Ancestral Process.Rmd
 create mode 100644 analysis/The Coalescent Process.Rmd

diff --git a/analysis/The Ancestral Process.Rmd b/analysis/The Ancestral Process.Rmd
new file mode 100644
index 0000000..e71d295
--- /dev/null
+++ b/analysis/The Ancestral Process.Rmd	
@@ -0,0 +1,195 @@
+---
+title: "The Ancestral Process"
+author: "Jennifer Blanc"
+date: "3/14/2019"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+rm(list = ls())
+library(reshape2)
+library(ggplot2)
+library(expm)
+library(ggpubr)
+set.seed(12)
+```
+
+### Pre-requisites 
+
+* The Coalescent Process
+* Continuous Time Markov Chains 
+
+### Ancestral Process  
+
+In coalescent theory, the ancestral process is the number of active ancestral lineages at time t. Recall that we can think of the coalescent process as modeling the relationship back in time between a sample of haploid individuals at a single neutral locus, in a population of constant size 2N. In the previous vignette, we showed the time until two individuals coalesce, $T_2$, is an exponential random variable with parameter 1. We then extended this idea to show that for a sample of size n, each coalescent coalescent time, $T_n, T_{n-1}, ... T_2$, is an independent random variable where $T_i \sim exp{{i}\choose{2}}$. So each time a coalescent event happens, we are transitioning from $i$ to $i-1$ lineages after waiting $T_i \sim exp{{i}\choose{2}}$ amount of time.  
+
+Keeping this in mind, we can define the ancestral process as a continuous time Markov chain where the states are the number of active lineages (1,2, .., n) in generation t:  
+
+$$A_n(t) := \text{number of distinct lineages in generation t of a sample of size n at time 0} $$
+
+### Deriving Kingman's Coalescent
+
+We will first re-derive the results from the previous vignette while explicitly considering the ancestral process as a continuous time Markov chain. Mainly, we will show that the coalescent process describes the ancestral process for a sample of fixed size n in the limit as N approaches infinity in the Wright-Fisher model. Remember that in the wright-fisher model we assume a fixed population size (in this case N) and that the j parents of i individuals are sampled with replacement from the N possible parents in the previous generation. This process is described by the expression below where i is the number of active lineages in the current generation and j is the number of ancestral lineages in the previous generation. 
+
+$$G_{i,j} := P(\text{i individuals have j separate parents)}= \frac{N(N-1)...(N-j+1)\delta_k^{(j)}}{N^i} \\ 1 \le j \le i$$
+
+This probability can be broken into three parts:  
+$\bullet \ N(N-1)...(N-j+1)$ is a descending factorial that describes the number of ways that you can choose j parents  
+$\bullet \ \delta_k^{(j)}$ is a Stirling number of the second kind, it describes the number of ways to assign k individuals to their j parents  
+$\bullet \ N^i$ is the total number of ways of assigning the number of i individuals to their parents (there are N possible parents) 
+
+Taken together, $G_{i,j}$ describes the probability of i individuals having j separate ancestors (remember that $j \le i$ because individuals cannot have more than one parent). Since our CTMC is defined as the number of lineages at a given time, $G_{ij}$ describes the transition probabilities for our CTMC:  
+
+$$P(A_n(t+1) = i|A_n(t) = j) = G_{ij}$$
+
+When we look closer at this transition probability, we can see that Kingman's coalescent does not apply when we have a fixed population size N. Right now, i individuals can have anywhere from i to 1 distinct parents. However, we know that Kingman's coalescent only allows only one coalescent event per generation, in other words, i individuals can have only i or i - 1 parents. We will now show that for large N, we recover Kingman's coalescent.    
+
+First let's consider the case where i individuals have i - 1 parents.  
+
+$$G_{i,i-1} = \delta_i^{(i-1)}\frac{N(N-1)...N(N-k+1)}{N^i} =  {{i}\choose{2}}\frac{1}{N}  + O(N^{-2})$$
+
+
+We can arrive at the above statement by first recalling that $\delta_k^{(k-1)} = {{k}\choose{2}}$ and recognizing that the other terms are proportional to $1/N^2$ and will decay quickly when N is is very large. 
+
+Next, we will consider the other two possibilities. First that i individuals have fewer than i - 1 parents and the case where they have exactly i parents:  
+
+$$G_{i,j} = \delta_i^{(j)}\frac{N(N-1)...N(N-j+1)}{N^i}= O(N^{-2}) \ \text{for} \ j <i-1$$
+
+$$G_{i,i} =  \delta_i^{(i)}N^{-i} N(N-1)..(N-i+1) = 1-{{i}\choose{2}}\frac{1}{N} + O(N^{-2})$$  
+
+Now we can show formally, that in the limit as N tends to infinity, the ancestral process converges to the continuous-time coalescent process. In order to do this we are again going to re-scale time so that one unit of t corresponds to N generations.   
+
+$$P(T_i^{(N)} > t) = G_{i,i}^{[Nt]} \rightarrow e^{-{{i}\choose{2}}t} \ \text{as} \ N \rightarrow \infty$$
+
+Here the notation $[Nt]$ means only the integer part of Nt. So while t could be any value greater than zero, the probability of $G_{i,i}^{[Nt]}$ only makes sense for whole generations. We have shown that the probability of a coalescent event, $T_i$ happens after time t happens with probability $e^{-{{i}\choose{2}}}$. This is that same coalescent time that we arrived at in the last vignette.   
+
+### Properties of the Ancestral Process  
+
+In the last section, we showed that as N tends to infinity, the ancestral process converges to the coalescent process. We can use the ancestral process transition probabilities, which we defined above to write the general transition matrix for this process:  
+
+$$G_{i,j}^{Nt} = (I + N^{-1}Q + O(N^{-2}))^{Nt} \rightarrow e^{Qt}, \ \text{as} \ N \rightarrow \infty$$
+
+Thus, the ancestral process is approximated by the Markov chain $A_n(t)$ whose behave is determined by the rate matrix Q and $A_n(0) = n$. Here Q is is an n x n matrix with non-zero entries given by:  
+
+$$q_{i,i} = -{{i}\choose{2}}, q_{i, i-1} = {{i}\choose{2}}, i = n, n-1, ..., 2 $$
+
+$$Q =
+  \left[ {\begin{array}{ccccc}
+   0 &  &  & \\
+   {{i}\choose{2}} & -{{i}\choose{2}} &  & \\
+   0 & {{i}\choose{2}} & -{{i}\choose{2}} &  \\
+   0 & 0 & {{i}\choose{2}} & -{{i}\choose{2}} \\
+   \vdots & & & & \ddots \\
+  \end{array} } \right]$$
+
+Intuitively, we can think about $q_{i,i-1}$ as the rate individuals coalesce, going from i to i-1 lineages. Therefore, $q_{i,i}$ is the rate that process is leaving the state corresponding to i=1. Since these are the only two options, we know the $A_n(t)$ is pure death process that starts at $A_n(0)$ and decreases by only one jump at a time. We can use Q to calculate the distribution of $A_n(t)$.  
+
+One way we can do this is write $Q = RDL$, where D is a diagonal matrix of eigenvalues $\lambda_i = - {{i}\choose{2}}$ and L and R are matrices of right and left eigenvalues of Q, normalized so that RL = LR = I. Using this method and some calculation, its possible to solve for the the distribution of $A_n(t)$:    
+
+$$G_{n,j} = P(A_n(t) = j) = \sum\limits_{i=j}^n e^{-i(i-1)t/2}\frac{(2i-1)(-1)^{i-j}j_{(i-1)}n_{[i]}}{j!(i-j)!n_{(i)}}  \\ s_{(n)} = s(s+1)...(s + n-1) \\ s_{[n]} = s(s-1) ... (s -n +1) \\ s_{(0)} = s_{[0]} = 1$$  
+
+The expectation of this distribution: 
+
+$$E[A_n(t)] = \sum\limits_{k=1}^n e^{-i(i-1)t/2}\frac{(2i-1)n_{[i]}}{n_{(i)}}$$
+
+Now we have a both a distribution and its expectation for the number of ancestral lineages for a sample of size n and different times, t. Lets plots the mean number of ancestors over time for different values of n:  
+
+
+```{r}
+up <- function(a,n) { # Function for s(n) calculations
+  x <- seq(a, (a+n-1))
+  return(prod(x))
+}
+
+down <- function(a,n) { #Function for s[n] calculations
+  x <- seq(a, (a-n+1))
+  return(prod(x))
+}
+
+
+mean_ancestral_lineages <- function(t, n) { 
+  tmp <- rep(NA, n)
+  for (i in 1:n) { # Evaluate espression for every value of i
+    tmp[i] <- exp(-i * (i-1) * (t /2)) * ((((2*i)-1) * down(n,i)) / up(n,i))
+  }
+  return(sum(tmp)) # Return sum over values of i 
+}
+```
+
+```{r}
+# Calculate the mean number of ancestral lineages for n = 5, 10, 50 
+x <- seq(0, 2, 0.01)
+df <- as.data.frame(cbind(x,(sapply(x,FUN = mean_ancestral_lineages, n = 5)), (sapply(x,FUN = mean_ancestral_lineages, n = 10)), (sapply(x,FUN = mean_ancestral_lineages, n = 50))))
+dat <- melt(df, id = "x")
+
+# Plot the data 
+ggplot(data = dat, aes(x = x, y = value, color = variable)) + geom_line() + theme_bw() + xlab("time") + ylab("Expecent Number Ancestral Lineages") + scale_color_manual(name = "n", labels = c("10", "50", "100"), values = c("red4", "navy", "darkgreen")) + geom_hline(yintercept = 1, linetype = 2)
+```
+
+
+Again, we can see that most coalescent events happen quickly as the number of active lineages quickly decreases and that most of the time is take up by the last coalescent event.  
+
+
+### Example  
+
+Let's do an example for a sample size of n = 5. Recalling that the ancestral process is a pure death CTMC with death rate of ${{i}\choose{2}}$, we can write down the Q matrix for our example:  
+
+$$Q =
+  \left[ {\begin{array}{ccccc}
+   0 &   \\
+   1 & -1 &  \\
+   0 & 2 & -2 \\
+   0 & 0 & 6 & -6 \\
+   0 & 0 & 0 & 10 & -10
+  \end{array} } \right]$$  
+
+
+Now, we can use the fact the the transition matrix, $G_{i,j}(t)$, is equal to $e^{Qt}$ when $N \rightarrow \infty$ to calculate the transition matrix for different time points.  
+
+```{r}
+# Enter Q matrix 
+Q <- matrix(c(0,0,0,0,0,1,-1,0,0,0,0,2,-2,0,0,0,0,6,-6,0,0,0,0,10,-10), nrow = 5, byrow = T)
+
+# Calculate transition matrix at different time points 
+P0.1 <- expm(Q * 0.1)
+P0.5 <- expm(Q * 0.5)
+P1 <- expm(Q * 1)
+P10 <- expm(Q * 10)
+```
+
+
+Let's visualize the transition matrix with heat maps for each time point:  
+```{r, fig.width=12}
+# Prepare data for plotting 
+dat_plot1 <- cbind(melt(P0.1, varnames = names(dimnames(P0.1)), na.rm = FALSE, as.is = FALSE, value.name = "value"), 0.1)
+dat_plot2 <- cbind(melt(P0.5, varnames = names(dimnames(P0.5)), na.rm = FALSE, as.is = FALSE, value.name = "value"), 0.5)
+dat_plot3 <- cbind(melt(P1, varnames = names(dimnames(P1)), na.rm = FALSE, as.is = FALSE, value.name = "value"), 1)
+dat_plot4 <- cbind(melt(P10, varnames = names(dimnames(P10)), na.rm = FALSE, as.is = FALSE, value.name = "value"), 10)
+
+# Make a heat map for each tranisition matrix
+p1 <- ggplot(data = dat_plot1, aes(x=Var2, y=Var1, fill=value)) + 
+  geom_tile() + theme_classic() + scale_fill_gradient2(low = "white", high = "navy") + geom_text(aes(Var2, Var1, label = round(value,2)), color = "black", size = 4) + scale_x_reverse() + scale_y_reverse() + xlab("j") + ylab("i")
+p2 <- ggplot(data = dat_plot2, aes(x=Var2, y=Var1, fill=value)) + 
+  geom_tile() + theme_classic() + scale_fill_gradient2(low = "white", high = "navy") + geom_text(aes(Var2, Var1, label = round(value,2)), color = "black", size = 4) + scale_x_reverse() + scale_y_reverse() + xlab("j") + ylab("i")
+p3 <- ggplot(data = dat_plot3, aes(x=Var2, y=Var1, fill=value)) + 
+  geom_tile() + theme_classic() + scale_fill_gradient2(low = "white", high = "navy") + geom_text(aes(Var2, Var1, label = round(value,2)), color = "black", size = 4) + scale_x_reverse() + scale_y_reverse() + xlab("j") + ylab("i")
+p4 <- ggplot(data = dat_plot4, aes(x=Var2, y=Var1, fill=value)) + 
+  geom_tile() + theme_classic() + scale_fill_gradient2(low = "white", high = "navy") + geom_text(aes(Var2, Var1, label = round(value,2)), color = "black", size = 4) + scale_x_reverse() + scale_y_reverse() + xlab("j") + ylab("i")
+
+# View plots 
+ggarrange(p1,p2,p3,p4, common.legend = T, ncol = 4, labels = c("0.1", "0.5", "1", "10"), hjust = -.3)
+```
+
+In each transition matrix, each cell in the lower diagonal matrix represents the probability of transitioning from i lineages to j lineages for the given time period. If we focus on the fifth row in each transition matrix, we can see that for t = 0.1 the highest probability is that one coalescent event will happen and the process will go from 5 lineages to 4. If we look at t = 0.5 and t = 1, the highest probability is to go from i = 5 to j = 3 and j = 2, respectively. Unsurprisingly, given that the coalescent is a pure death process, by t = 10 all of 5 of our samples will have coalesced and there will be only 1 lineage left. 
+
+```{r}
+sessionInfo()
+```
+
+
+
+
+
+
+
diff --git a/analysis/The Coalescent Process.Rmd b/analysis/The Coalescent Process.Rmd
new file mode 100644
index 0000000..95b1f07
--- /dev/null
+++ b/analysis/The Coalescent Process.Rmd	
@@ -0,0 +1,193 @@
+---
+title: "The Coalescent Process"
+author: "Jennifer Blanc"
+date: "3/13/2019"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+rm(list = ls())
+set.seed(12)
+library(ggplot2)
+library(reshape2)
+```
+
+
+### Pre-requisites 
+
+* Wright-Fisher Model 
+* Intro probability  
+
+### The Coalescent Process
+
+In population genetics, we are often interested in modeling the relationship between a sample of genetic variants. The coalescent works backwards in time to model the genelogical relationship between the alleles in a sample. Coalescent theory was pioneered by John Kingman in the 1980s and remains one of the most widely used population genetics models. To understand how the coalescent works, it is helpful to first think about the wright-fisher model where we have idealized population with 2N haploid individuals, a neutral genetic locus, a constant population size of 2N, and offspring being randomly assigned to parents. While these might sound like idealized and unrealistic conditions, we can still derive interesting and applicable results from this model. Furthermore, it is the foundation on which a variety of more complicated models are built.     
+
+### Time to Coalescence  
+
+To understand how this model works, the first question we can ask is about the time it takes for two individuals in this population of 2N to share the same parent. In other words, how long does it take them to coalesce? Since we are thinking backwards in time and randomly assigning offspring to parents, the probability that the two individuals shared the same parents is $1/2N$. So the probability that they coalesce in the previous generation is:  
+
+$$P_C = \frac{1}{2N}$$  
+
+By the law of total probability we know that the probability that the two individuals don't coalesce is:  
+
+$$P_{NC} = 1 - \frac{1}{2N}$$  
+
+So if we want to know the probability that these two individuals coalesce in the previous $k+1$ generations, it is just the probability that they do not coalesce for $k$ generations and then do coalesce in the $k+1$ generation:  
+
+$$P_{C,k+1} = \big(1 - \frac{1}{2N}\big)^k\big(\frac{1}{2N}\big)$$  
+
+If we consider a real population, 2N is likely to be large and we can use the Taylor series expansion for e $e^{(\frac{-1}{2N})} \approx(1-\frac{1}{2N})$to approximate the above expression:  
+
+$$P_{C,k+1} = \frac{1}{2N}e^{\frac{-k}{2N}}$$
+
+Looking at this expression, we can recognize it as the probability distribution of an exponential random variable with parameter $\frac{1}{2N}$. Immediately we can derive some interesting results from the expression, namely the expectation and variance of the time to coalescent for two individuals. Calling this time, T:  
+
+$$E[T] = \frac{1}{1/2N} =  2N$$
+$$Var(T) = \big(\frac{1}{1/2N}\big)^2 = 4N^2$$ 
+
+
+So the time to coalescence for two individuals is 2N generations is exponentially distributed with a mean of 2N. We can easily plot this distribution and the expected coalescence time: 
+
+```{r}
+N <- 1000 
+
+f <- function(t){ # Time to coalescent for 2 individuals
+  return(((1/(2*N))* exp(-t *(1/(2*N)))))
+}
+
+# Plot
+curve(f,1,6000, ylab = "P(coalescencence)", xlab = "t")
+abline(v = (2*N), col = "red")
+```
+
+
+## Multiple Coalescence  
+
+We can extend these ideas to model the coalescent times for a sample of size n in a population. Instead of just asking when two individuals coalesce, we instead want to model the time until each pairwise coalescent event until we have reached a single lineage, the common ancestor of all the individuals in our sample. In other would we are interested in in $T_n, T_{n-1}, ... , T_2$ where $T_n$ is the coalescent time of 2 of the n samples, $T_{n-1}$ is the coalescence time of the 2 of the remaining $n-1$ lineages, etc. $T_2$ is the last coalescent event of the 2 remaining lineages into 1.  
+
+First we will consider the case where $n=3$. There are now ${3}\choose{2}$ pairs of lineages, any of which could coalesce. So we can re-write the probability of that 2 of the 3 samples coalesce as:  
+
+$$P_C = \frac{{3}\choose{2}}{2N} = \frac{3}{2N}$$
+
+Just as before, the probability that the first coalescent event happens in $t+1$ generations can be written as:  
+$$P_{C, k+1} = \frac{3}{2N}e^{-\frac{3k}{2N}}$$
+
+Note, that after this coalescent event of 2 of our 3 lineages, we now have only 2 lineages left and the time coalescence for those two is the same as in the first section. We can now extend this idea to an arbitrary number of of samples, n.   
+
+$$P_C = \frac{{n}\choose{2}}{2N}$$
+
+$$P_{C,k+1} = \frac{{n}\choose{2}}{2N}e^{-\frac{{{n}\choose{2}}}{2N}k}$$
+
+We have now written the probability of a coalescent event in a sample of size n in $k+1$ generations. Again, this is an exponentially distributed random variable with parameter $\frac{{n}\choose{2}}{2N}$. Using this distribution we can write can calculate the expected value of the time to first coalescent event in a sample of size n as:  
+
+$$E[T_n] = \frac{2N}{{n}\choose{2}}$$  
+
+## Simulating Coalescent Times  
+
+Finally, we will simulate the coalescent times for an entire sample. Before we do that, we are going to re-scale our expression slightly. Since we are typically working with large populations, the expected time for two alleles to coalesce ($2N$) is usually a very large number, which can make looking at coalescent times on the scale of generations difficult. Instead, we are going to set $k = 2N \mathrm{ generations} \Longleftrightarrow t = 1$. So each $t$ time-step backwards, actually represents 2N generations. We can now re-write the probability of a coalescent event before $t +1$ as:   
+
+$$P_{C,t+1} = {{n}\choose{2}}e^{-{{n}\choose{2}}}t$$
+
+This is still an exponential distribution that can be written as:  
+
+$$T_n \sim exp{{n}\choose{2}} = exp\big(\frac{n(n-1)}{2}\big)$$
+
+Keeping this in mind, lets simulate the re-scaled coalescent times for a given population:    
+```{r}
+simulate_cumulative_colescent_times <- function(n) {
+  times <- rep(NA, n) 
+  times[n] <- rexp(1, (n*(n-1))/2) # First coalescent time is distributed exp(n(n-1)/2)
+  for (k in (n-1):2) { # After each coal. event we have 1 fewer lineage - so lets loop backwards to last coal. event T2
+    times[k] <- rexp(1, (k*(k-1))/2) + times[k+1] # Coal event happend after a exp(n(n-1)/2) distributed waiting time, add to previous time
+  }
+  times <- c(times[2:n], 0) # Account for starting time of 0
+  return(times)
+}
+
+```
+
+```{r}
+time_sim <- simulate_cumulative_colescent_times(100) # Simulate coal. times for 100 individuals
+num_individuals <-seq(1, 100, 1)
+data <- as.data.frame(cbind(time_sim, num_individuals))
+ggplot(data) + geom_point(aes(x = time_sim, y=num_individuals)) + ylab("Number of Lineages") + xlab("t") + theme_bw()
+```
+
+
+In this plot we are looking backward in time on the x-axis in rescaled coalescent time. Each point is a coalescent event that happens after an exponentially distributed waiting period. So our coalescent times $T_n, T_{n-1}, ..., T_2$ are the distances between each point. From the plot we can see that most coalescence events happen very fast and the longest waiting time occurs before the last coalescent event. We can also see this from our expression for expected waiting time, $E[T_n] = \frac{2N}{{n}\choose{2}}$, the larger the $n$, the shorter the waiting period until the first coalescent event.   
+
+
+## Time to MRCA  
+
+Another question we may as is how long does it take for all of my samples to coalesce? Lets try to answer this question using our simulation.   
+
+```{r}
+mean(replicate(10000, simulate_cumulative_colescent_times(100)[1]))
+```
+
+According to our simulation, the re-scaled time to the most recent common ancestor of all 100 of our individuals is 1.9888. We can actually solve for the $T_{MRCA}$ and check how well our simulation matches with it.  
+
+$$T_{MRCA} = T_n + T_{n-1} + ...T_2$$
+
+$$E[T_{MRCA}] = \sum\limits_{i=2}^n E[T_i] = \sum\limits_{i=2}^n \frac{2}{i(i-1)} = 2\big(\frac{n-1}{n}\big)$$
+
+Using the identity $\sum\limits_{i=2}^x \frac{1}{i(i-1)} = \frac{x-1}{x}$:  
+
+$$= 2\big(\frac{n-1}{n}\big)$$
+
+Plugging in 100 for n:  
+
+$$T_{MRCA} = 2(99/100) = 1.98$$
+
+Our simulation is approximating the time to the most recent common ancestor well. Remembering that we have re-scaled time, we can recognize that no matter the size of the sample we take, all of the individuals will coalesce in before 2*2N generations.   
+
+
+## Tree length  
+
+$L^n$ is the total length of of all the branches in our coalescent process. In other words, it is the sum of all of the coalescent times multiplied the number of active lineages at that time:  
+
+$$L^n :=  2T_2 + 3T_3 + ..nT_n$$  
+
+Let's use our simulation to estimate the expected value of $L^n$:  
+```{r}
+mean(replicate(1000, sum(abs(diff(simulate_cumulative_colescent_times(100))) * seq(2,100))))
+```
+
+It is fairly straight forward to find the expect value for $L^n$:  
+
+$$E[L^n] = \sum\limits_{i=2}^nkE[T_i] = \sum\limits_{i=2}^{n}\frac{2}{1-i}$$
+
+Plugging in n=100:  
+
+$$E[L^{100}] = \sum\limits_{i=2}^{100} \frac{2}{i-1} = 10.35476 $$
+
+Again, our simulation matches the expected total tree length really well. Finally, we will plot $E[T_{MRCA}]$ and $E[L^n]$ for different values of n:  
+
+```{r}
+L <- function(n) { # Calculated total tree length
+  if (n == 1) {
+    return(0)
+  } else{
+   return(sum(2 / (seq(2,n) -1)))
+  }
+}
+
+TMRCA <- function(n) { # Calculate TMRCA 
+  return(2 * ((n-1) / n))
+}
+
+n <- seq(1, 20) # Get Ln and TMRCA for values of n 1 - 20 
+L_n <- sapply(n, L)
+TMRCA_n <- sapply(n, TMRCA)
+
+data <- as.data.frame(cbind(n, L_n, TMRCA_n)) # Prepare data and plot! 
+dat <- melt(data, id = "n")
+ggplot(dat, aes(x = n, y = value, color = variable)) + geom_point() + geom_line() + theme_classic() + scale_color_manual(name = " ", labels = c("E[Ln]", "E[TMRCA]"), values  = c("red4", "navy")) + ylab("Time") + theme(legend.position="bottom")
+```
+
+```{r}
+sessionInfo()
+```
+
+

From 218a36c64b1596345ae6ab5b2581aafe83659900 Mon Sep 17 00:00:00 2001
From: Jennifer Blanc <jgblanc@uchicago.edu>
Date: Fri, 5 Apr 2019 09:27:22 -0700
Subject: [PATCH 2/3] New Coalescent Rmd and Picture

Hi Matthew,
I'm not sure I'm doing this right but I updated the Coalescent page and added pictures plus another paragraph that walks through the example in the picture. I wasn't sure  if I was allowed to include pictures from the internet...I included where I found them from. Let me know what you think and anything you think would make the page better!
Thank you,
Jennifer
---
 analysis/Cont_Coal_Pic.jpg          | Bin 0 -> 13591 bytes
 analysis/Discrete_Coal_Pic.jpg      | Bin 0 -> 55338 bytes
 analysis/The Coalescent Process.Rmd |  29 ++++++++++++++++++----------
 3 files changed, 19 insertions(+), 10 deletions(-)
 create mode 100644 analysis/Cont_Coal_Pic.jpg
 create mode 100644 analysis/Discrete_Coal_Pic.jpg

diff --git a/analysis/Cont_Coal_Pic.jpg b/analysis/Cont_Coal_Pic.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..211a9774c8d5c6e085902f88d5b926f2badbd2dd
GIT binary patch
literal 13591
zcmb_?2Ut{1((WNhR1hRbB?w3kqvQdVEIEUKN(RYsNF!Nt5>zDTj7ZJ`0s_KF7Lc4X
zLmXt7xug5-7xv%X|K8nu=IJwip3~i@tGcSX>a9k9LoWc=mF1P>0W2&4Fva`;^db-~
z=WS~R0II6M9RL6bfGb$!01l>xh4}$kbikEgbpTMtqW`_Fg~j<>8*Bgww*_#2Yh#4@
z$6Vybn48e=|2SD#e~g%g{rA>b+F3Y%uLDOv`=NIMQ4Jd>cPBR+C+C~Od;);zeI-@g
zpQB^0Uv>6h)$3!7(Nl522i(mUgl~2b_XE%$0Wt#M2um9Wiv_?Y!@?oMLVp67F$7$}
z`gQ&4hPh#3<KSMw!zUmlBE~eRy$)bw;oxB7;#|4%^D|ifnD+oK*%k6zg0grwG@s(L
zxKIcMCS((^-mmJQ)EYiu6Mp6zL`XzMO+!n^ew%}n>y8NM?mba4@dpp(<P{W^l(ltq
z_4Ex4jm#}9t*mWq?cCfwJiWXjKEWZOVc{<$A`_EdC8xZ8^ENdnH!r`Su&B7?Lv>AU
zU427iQ)kzw?w;Pyef=Y&W8)K(-=?M)mzGyn*VZ>Sx8R3I$0w&}i1Qym>B8{yZ)9Qa
zf1~Ul=pw_=g^i1ggNy%@E-Y-XpM;a)Ub!WRM=q<0|J3CMi%=i|#r=frst!U{VXXto
zXRgCURBR%P?C_tY{i5vuj<BHrh_c@a`-83--~oV#S)41lczBp4!p9>Zx=utyNJvCM
zLVEQ&H3ba~H3bzF9V0t49X%@p71d3io2<7vxbAS#GV==Xatg3>a&e+DgJ9v~;}a1O
z-5?^m!AVC&$NB$!p=U8-Nrj#UuHs-}gb9ZX00Tb?B000L|0o=cHb73bJFl1O<Su*D
z@IFg(eI$f;PBB&{fu@x7OTf(_QI?uefI?I5zb;Vp*I!=;(mV?Gn8Mz<k;9zW)J!zT
ziaA$lum!ugaG9NDDVN06-r{(B^?FUVvB9rIZTu)7yIUVuVM1B2I8KYDAD?^Z5~iKD
zupSAVU6+N`@<Jv<{G2@f2p73qpLHIhflqYZMIJTPio>gS-leLiCtP>va_!m&@I}i9
z%;*hFHWJK?jpq8?c*jz>3ZH=*?g{Q6-8xUZy?Xe4Jz)NIkareIR+_%pvey`&2tuTl
zL6xy5pMPp@PY?|-#V@>NhgnS%RbD(Q%mNL0mo8CdKV3xdyQ2Z(t#%s+Ek^Mr;*>z9
zYa=Y{5k!p+3w-npL5hsNGqib)wPCirYqiN=g=Gb`;F)_>_Xje*lS^^p3<S__Ds=HO
zOta@+V>#0uE8|-TB)l%-RgX>E7N|d9zPY*{%H}hYfoczh)CWm2W^eiMj-A!980vTX
zgjJNUD>gJm9@-?T(KE1c&YY=<;jY#9cFkXnI}dZatVRQb_((sMxOynesTB>x`p??|
zTSxNt$!H+dcQIk~5gK>{L$Wc$;rvMGZ*4B;?ZAts2n+X7@S5;%1E2xs`J<VA=yJQ)
z^b7H_;nY>=ziSeQ&~Kw<6oH*bqk&FYOlNX5(0hhFy_iSr&z@$Hl!`=m$Vw*<Ks#Fx
zK5rI4zgVNlNfC(GC~*6F{T!xkz<w<V#f*3=53P0v*FhI=qJg*w=<)^1*MA@BwRfR$
z<sF#8+x6&^mW|$Y8M$3m^{r61pn<bjckiaqh-tVAxQc>?@$lMo-*x&F1TQO$mT5^`
zIAD8wT)%+5V{>V7<>ZdGAkT`2bcn9USszpcF2Kk)BYe#Z)6l?kG~lE}ug>Xxaur&}
z>V=XxBiJ<>%H%6r=rY%}RZiQ<dsm;TkmHyO>L(3H(ZJ}OkUhU<X!aK#nj!S(Xkm2L
zzrFm+^z{@V<5tm=+-QL4tc$haEgBHaL@0unm&_0-G>|dxiv}1k+Q$AJdjJnZ%mgs(
z{ns31KB)8)61UAfqD2EC!6<2V3={XychB14^Ub?5jCwYMuI>M>i2{;){%UOMKczqe
zUpH0wapk{hq4~cMhxOlcK%zIFzBju<A4Pc<rhx4FO<aGR|4C5vUm;a%N8L%uJM41}
zLltXD5sHr5U#zeB?BKsWx>1A74SaFWeXZiUo2LRFKlN2!u&wl|^oy&c<$U%r)d%2E
z$MPo);cC|?n%=|dGHX-@H)(JDO#*e67tg|c{{?e3Cwi3=6Pv_a{!2x%{{O(_7jA!B
zrqh3+&P~_&)#^re=ZHCxMH#k^v|dzS5<ibfxNYtAGU7Q42ykN0QfSIt!4zin@GtXC
zuxDlL3%MOH_E75+i*Oa-sf(MmgLS(Fx>yjJWly&s=?H2kqHcd}ka{|@t@|}Ex(45B
z!jL`s$sTMVUS2qG(3h4ruY*xghJoV~cj-3}jPwDUPvgN#Iaj;?9YTwwL5%_mC?_qI
zke5Z#N~qo~CH%0F<84JrQs3NdX^Kgf4;XQdgiz6Q)9M*hKui`hm^qnltcb-q<FICg
z&k|u2=y*27y%+HnKStHOQt9B(VDFvYcv%^!lsk&BV)wDxpf$xob4`K+$xuVQYxWyH
zXIH{gV&c%|y2wQ;3ndDnvf(%LOdH?s_`k&U!Y(8ZMvRoa@bT5a=RFuj$<4qnu5Qi$
zaA~M0X*3JVJmi~$ezVp?1EVqB9y)-g9L2AqZLv9>t{QIv4F&vQi|VvnM}C^&hAC4I
zQ(H$i4>#w`1a4WXwR-hX2%k!rh3*df8AP6SNVANdkP?1xOTR2$EAP6)EzfvAKI^VN
zZ^}(F<h985#douW?KZsbou{|q3j%p1!TfjfjK4~YJ)t}lkG_{(X*e!2aP8S(*@-c?
z{kHtg8frnFqHkifdbo6RS}XQdE1F$*##Sn=@^*HsKKs$b9p)pEAh+8E!+j7-_la^5
zbDK0=J#q&WbBoStf0}vLe&8u+HG}1#6&76C`%;qHF=2)S>2EvCsl$cpkDZ%sb}``(
zM%}bwy7PKx+F=3pRfp%aqO457c}9uGkBLyrWuB5e5|LXq9#*lxW4FnUIM1yZ<P~5z
zE+naaJI(u>Q2i9FE4T8qzTz0bQ{m<KkC=qo;B^IEUzECxrt8gPbx`!I;20lp97ug9
zj$q5nZpOkuS>B(qP4`l<3>HdfeQF)VgY{U+^;#G7W7|GTqOLA@H8@eSMHlq)lYG?F
z?fI3L0i;{jng^CGTI5#M=FiWM87enqklYQFYfM+sKm_ThL8t;Wg?RY1NCI$c#zFmf
ziUzPM)$s~S&T8%rY>8=V`u3_f7;)xQOWgm%Z~Em!k)lWb(Vxm9(O)b?M%8DY(vG(u
zFT0rNy0%{fub8x;f%ik@XrPcja;pf)iR8qTYXEJ5qBGWJoRN!29|3Wq3{K;&46`W0
z?4z(OTuy|RZ&^GbO2j9`?a(sW+-61YW;C`6rxi>>4S+E#6oQWK1tXnGOfVL1GDNds
z!cnVxM^`NDuEnFU?MCddX93z_3l`>^Z?=t^b72E%k2(5~U-TYFL&eW$0~irei>Pwz
z3L{YR-rGE77*ih&j7N@^^i96ZV{V4{L6|6GTV#PU4j<{w)05)pbG2xBCeWl+oqu%Y
zwe30|R5ImB@V)|#B^W~A^nOPk3ZpD3&qzL+{6y;k#`a#BL8zUo-kmrdGr4=qknEi=
zp4?Re;0g0%(QcphE`G9srAj;lyE7HK#{+9@Is8pu5``Q!TxA@xU#w#_Wz=eiz83m6
zKpf%bm~6tqh(r!miZrgJl4j>K+DkCz49iix;LM2xxajkFS>3Zpg%I7`^KGG7##Q9=
zCfMCRq0D?x_!4vhw0Zqu-GKsyGL`7qouVD!Bjs8O7!F3j+`x17BR3DnV&7Fo2Vi#l
zvq7Iwdu_!Yj`X0*Rpr{|X8wnexNgAR>;}T^Q&fC{$|(ePfe*gmf{wHb2@&9Wtu;a;
zQuR;kgHAYKS=r<s<%;%wZ-MX0D=cgIiH?J7k<10(g{^er`jRrdI0lkVdZ2V{=f#EJ
zxB|ARMcosdF4Zz3Im_Z6{-om11wHzo3ey7_JR1&*kWnv^sVt(@p>p$_&FU?>@-Los
z8rMctZ+&D~Aazv&u=?1r1|ea!OOx=+g|B|p5};gY@9JVq?zh%Sd%#(iG0eZmuCklE
z80iJIUCboraZ|T+<y*~=XSv>A_Rbhc^Sg9*cIl{yTVxJ&l4MCupD35Ao}WLr88)hx
z)XcLyMSdwUjC<qiEJ`F{vG=<AG0bhCSj?P#XA{Y2wh@B{t|GW)56%+!((|X-4;88^
z{X&mcV2|5`VFfNTBTaITnPw$4xMsJh_IxygLUi;#bJ5%wNr{xB^_d9G-*#3YsVy+r
z5@dDKX%4)t``rAdfGB8fy5gPh!HkZ-{sV<~UDfw~$mO7c0j4SD0@4M4`$<)6@8CyK
z247FE2`$dvYOPP4nw35b*UF0vuv_Fad+_c?m~Qe!qOmY_XMn`PdZe|Dh=~$MqXDuf
zA|as;-oJ<AOiRn1@OeGKL^&-}`TSyGKcwPy@1)S2I(A{Wg?G|t{L5S4l&BP0S0Cev
z_N_Y=N@Jn~;coE+&6J>yiF7+<H$7i+O*zW9V-%W-Ho&0n*xty1T4VF8B~Dh`#}}M^
zC|zsO?Z^B*CW=$uLN=_TsfTibGTU)2Pyqs@-6E4Mo7<Az#T}vIDGC(Pw1h2PShXv7
z1Vo?W>Mz5A2CR5Z<EB6OaTmmJbWQueq!H*WncI!7LiUy0c~}nK?hel1DeB>O`2c*@
z4bwSLD2Ehwrig=IO3E!Zl3Zz-w+Ug&SL!{SI-e}CYU%BMa6!CO5MEL^fS1*7`7RCH
zt)An^>|!mQH%-mY_tdq}?A#<X_v`V2n~YcipCZKpH|sFf|FDN+|CtEqnMZ4l``_Op
zrst}tWRIn{(dlDLygM$(kf>h3iLKCK>f_J7273u5L#n!!ZKpJu%uhXc;A8us!kMGz
z8xvFix@Td?$AL9WGssmv%YOcPF8pSc!U4mem#@cRpQu9g$98$3xYy5pb4*`h5ggt`
z70g(+TWI)w3+^}Xlrl@FF~3tMl$!|9EKk%$uIQ(%ztzj_hD~XCM$V$b3WJei=8NUL
z4m|pNkkn_=w#JR%ub$-#Vg)7fN*p3^!@D!hQbg{H%6G0C-LlDvlnbTM66^-9R84<C
zs+UwmEe@({3qsn>>k?h4T-GSH?gwwV#Pt*-3cIJ7Lze?<cV&4WeLA<VN}p#`)Wuz}
zxVp>`IM)D{K3Yo2qZ(eRJ4kM}_$+Flw2SiI-_)lu)354__~a{bU%YQ8_Q3)y;^Z#;
zLV;V(i}|GrOe|8W)Um(vrpF5MExYrSnB{s!K%{2uw=xIfF7ewhOpZ-hI(rN$coL0g
zawK0sAuk33ZrABcZK^#U+9>0t)hz7ea2eB-QyQTNO9_^_=9s+3B*HwCc>@W~jZ3AM
zJuMD*dd%Dbpn*XB!m5vhkhm^&L^=ZYsf`6JB5zhLa~Zh;-6_dTw^x~Mi32n5H@QRK
zqKLIebzC-;vG!#Oy~^Tb?vsdC+*oFM<K#!r+ljMo@VTC7U<|8^`6D@$X#Y-IwgcmC
z^l3HUMnO#O(=CFIl8|z?FM=YiwOrbJkxZ0hVO5J4#0T<Lp@ttfm}#COuFKh5ryj)x
zQ)@jawDV!?Lj#JVe1?J>d6}1Q?yZ}~l8S6`E9J=(Rn}(HN}&Nw%nAh;{l)H%|5a#%
zoP})i9MfH77%IE8<4FGQ;DvU37>7a%15Z>vuaU3vv`C$PHttRhEmb;Qd1M4y`Eq8!
z3T91F=ig^+4#%rU19dHrl+Zw=N=5Z4Ikd(c6Md1MUv8p-_<*xkz~f!0Tl;l)QiO5~
zD>OI%V?Cr(NKszK!ovgCi&lf+cu+@tKB5`%vorHKp&Z7nfR8Rl(ZJ2ieV>o;`52_a
zB8<>3HP_pUlPz4!p3~HL#p7;aNXspGXM@01{{=#}1v;1U4LSK&^eH_YJG*(3o~aY{
zUaP!=1HE>vkKR9wQL%ObR-;>@Nt?tcKU4_zj<irO#5R^s8_wH`2P9p-PuS8cSX)tx
zhmbw0_@L-2um)yV>+WEY-Gmj!IDPtr20rIc8e9-!Jg8St_#zqrw}KBz=gvw?;ll6u
z^}0G|W2>cKSszwtU)X>8l!?e|+6F)6xKIWo4$wg6T;@6czVGm5?b0)Kq?U<C=KG|}
zd}gFK6t=I93_fe!7Jk1i#L2^*gqns6a-g{Qnqc3QP$O|DN$~nvAxidQ?l{~Cj=^N!
zs-Vol7>6zXmD!&rkwsGNRaLad5lfS5n5IC7{pE`@z^DZ3mu_5=)7oDYC;`P)C3<i+
zyvIys@7Gi&Up`m{Q+T5R(^2>8yY92kHn^uAkaN|AU>ioIQBCLaCJi`eITtLFIQTGp
zi&oA_jQ!xt>eel8S_<dxaoaMbik#WeOBt}TY<-)f&?RiCoO@mEsE0&`AvSr32{<@I
zvSd(HjbQvu8*j-Ti06WEOZ0BwH@s)S4p0_DwfuPa81WnfBjFS5E{&r)@lRhpe8+dT
zvYhDKeRjEZFoFDClo-%}?O(-#L)w9i=DBd;MStqx-cT768t9Ln<ocJa-T0R*J;8rE
zVWFAuIP-KeA_5Fo@|XzdR6qU8>gN9I#r%qIn4qYh@0vMddYcUmbhA$yUJ!z(F*rbP
zI1+V?kr%M<TyVMfn+!h;5+RuRe*LE~z<(EpuUBA<c9~O5WaIo4$q7%nX^%b6^-ox;
z+iuwAK01&KlOcZ`>^yhz%g1lnf0>$6S^>%kEV6}8xjwdoU2CK!7Hg<CDM>=MyLldZ
zGRe=?P8r*^v4&AZf(fp~c(@dXnke%SlVm#lRKd2hR=dk-&ppT(8t{#J)g)`YHq<3z
zq0jN3uq|OD-Aay?Nfgu%b?|<Zga#;Hn^^B)cwno5(Oc%<_4YwI;+-VoPXfFAo7T*q
zWs!K{$Rhd)`e*E~2Vh}X#)UB#`US-&m2`;gVzcwRK4JO!G5nhE$5-Rj@m6nN=zak7
zps)Piy=1<%5&)XooMVSmp2W0{F;UE3>nW=1ht^;mp?C-_^(@BZjr;HWOc#_0Nq%Q|
z#jTq?HNP+MC{pO(Dc8Rz<kN9%wR4sL+j6a|DuyO}6c_(_w0}5ev7vuj>n6rK_HP#w
z_8p}A?hNwSQ>U_hv1jFRO87VzXvQeN!}JQv^EdZ_6g}}*3!D%p+iT`y#EFxHuoGUK
z`6$k4t&0h<F&3CdK+3J#jrkuVD&tzfb&dxP3zy^)k}cm-y3$qnrV7{Ep;45LQ$kh@
zNFyQ5u8?>%S}D!8Wna9l_!VJC1-7bjOuDEqIst1)q~^l;7h`M7yfe-Wzw+_0Q5S{T
z?6Z4M-+M@K*FpX}fo2DHa0X;K3t_k<nl<_y$z38p;ZyM%+o(h(VV%c@f($i*uXZUf
z@6w{B(c)p9WU7|$q1KkLIrRe#B+P7z%Z0_+o~p9+jP!gZe|UG8$2MtfR|C%uVcK}y
z&60LO=e=!EJ+3!PpkLWeLG}Z#mgo2?Oauj4`rxO!#K_&{1D)|`{#0AfL#Be=J(E<c
zE6cXb@j?wts#F=W6hxg0As3COvNjLZ5Kr{y?rL$ODugST)i6;7@mJa2*FD;2%#GB(
z{r!+4UcS4xUyX1x-fpw&?d*>0lFOaQ)t{EkUUSkhfNdUX>1}GQN(f!vLlkRuNfzN}
ze4AEabL+PD6AhLbMTyv_h?CWS(v!crfoEZP=oRzr{N)o@tc73`Sm!*!;36r>paVq2
z&2;uVeRCX#w2|0D)#7+UK)xkjNJ*5o8{;fM@dk5|W+u`qHHkC#V^v;SYN>q3k9!tQ
zL(}<1>1Nu8({i1PRAbg41BiUScrqx-XS+V}P&#m0dAW|0Ffh_J(r8^OA>Zc$E^FF`
z)Go@xOWxk^zJhsH%XvZflX;Ydx+Z#?Jnq}?8e>_4?k2~*$3+9m<ZVH(UTG(7Kw4jl
zpY<3Q|Crf0fY>Rhe(M3>g<h<+N?t{=j<jqql;_sP71>kcT(W94T62tHLr<TLZ!Fc8
zU2gM*lRSQxziT!<i45VMoR1_rIa^AXv+#=f`5>>2aWZ~HZdugVr`Q_sxD<Puf5S?&
zx#6jAo#fp=-0~DoRukg%-F8BV+-*dj(fiB>4eTeG;CHmj%pqhQ@;P8163S-{j$(7q
zmk&zJ@E3jeciXRj=Q{GBFCa`n;&8-0uVOx+i|Smzz*(xR&#&)uDk_}Lfq7nh4{(Lf
zgT{HycP>~zW_#<z@?O?W_U)cCrN@#s<T{zTD3aKAHmYY`uWTi1R$3R?`UyMyDmduL
zzv>)D7yl8uP|^Pf$SWJsYXsPn;2&4=BTgL59TpWPC|-izdBr$b%6FwTT>H3MS8TPM
zP+wFIrT0@CmgFj;Cw^1u#5W<x7GZn+DFN6E!uBLnHYu9HRXA>umukBwHAI>dE=@N)
zCKdBq=MZ))uBeB!9Ty(GtH$dKlRM7{;LEc^1JxaFN|ytJ5eLlOKD`k&m#P>~Y>ao8
z?iOB7QR!N><2>yh?8JB4377B;NTqH8XOG)=xZ18^y{D&Ea0rEQ`<uoh#|oDEddPYL
zYZ!=cZyigAL`kYtxI2So6v>~fQfxOO^f8Ymhi^R@Ca1fnz5=i<VB3Eoed!p$y=AR5
zJhIm%UpK?PnSQ_DRc-i%Vn;_WM>gz{c8FN*hIz!Nk~;&2Qm;?GE||kBUrOKfeyu)I
z*oFoanp<M5MjR4k>!ZBCfBE*^lh2@<fm{@arS0|bakB}BdSrhmZqS_IVNTrGk4){k
z7;e?Nk+NGK-bX9os)RCs$<*>$Z#BWei%>&teM^(3?PyAxaH`yB32RucxZ*dt(9ZO@
z3+`W8yB0g=Ak&(^n=$ts(9oAc*bFO4d65f_w7mu;ElCl<%trF7#%>6L<+YDkJ)ihv
zHl>P}O1!iNCBj9oiuW0{K&;LSlw!CR*WJ36QsKtT4C)7Smq@WCZ^sjYVce9&>6=z$
z);JO`&wT=i<Ox=;#u<IBwdZwHt{<!E3pw>7Zn_6V9|aLe1`*s0{(Li%1E2goBGk|`
zan@4=<7NkEapnZe5ExT-`dN*|ybfj$5{{^3GoT=8xad8mLE1NQ1e-oxR!3y`sU4ig
zO4~T#PZuyMn0x~VVfI1naxQ+K%UFgU2MrhF0l6cIN@CkC#e>uIXLwfy?$kG%uYoh{
z4YalE8J=Zo%+6R_+{~t={EGSDNW%Ra9TogbYo(Xz3#GO9AJ$tsBr%?OV7Ar=`c^*r
zZ@p(w9GPTDu^j%vYjuw6j7=*z2tM$FZa(|I$+2iPVZ*RX?BTHA&_G9%^(i<^WQO62
zceV*eZ!?TfZa3^lZcVGzFV2(8-@QzQ-hL+jN^<#d{(Le&t8=Ll_A=G%v^}_%+rwn!
z<(I0!u@A_P0iwOjYUF2opSv;6>J)_aR4<b?`#$XhE*e%=%9|D1i?OAb<qNIidJZEG
zGp*g`xju4i2|($#dX0vZhk8xuBFjOb$)BN%e6}2g1{S8^|CUYtQ!ZBos}Bi!?(Mk`
zKBoHxJ{p;F!LP-IyBD^OF%z}0khZe=kRElrU?1tdN2o5*XOhYv+PQ2sRVG*6`w(0_
z#5nHzq1XxspIxS7!*5_wYK7@GA{f_btz5d&ma)2ztf^6$Gfmt<M$IQ^-RXk0YF$zO
zM|Bi+MDpF@4&RZP58^7&aPi()kOYi;Zi{TAw5Ulic&9+$LYIIjeY<8?1M61cM+pV!
z^WmBqe`e(mqh&hf<#~ioL@7gKazwdcnL8uV&iE8?Wo_4Ek*>vDS+Atov#WIX#n#eQ
zJ&DAM8b@3*scHy7XyA$|WJE;!8(v#o{Iwa2htDm%aH+JUNX-YI`CU0c`ar`hC)i#O
zje#CSyx|)oywQ|XlIU-+l7^h?Oj_Q5$u^9<nUeCF7-#DSv2O}Pi)Xj7FR7YmF-O)j
zZ&8-;on<K0Eyqr8koB>R=U(R1#%J_yRc?3iQ0&_8*E1eR!vpKaU=pc9J*B)lW+%|E
zC%f(Gd|WnV>JJ~tQq<4^N9_^J*TSPb=pEz-v`oZ*h=#UApRT#52&iyQGYZIBtxM5F
z>nU`VE1Th{Fb5;=W3x23md_vbn+(*c<D^(gPxa>ebkHs}LSIf}+h^PuGmjm21QC?4
zB!0Pyvwnpgkow?!&6cn-OyQcJu-56A$ZIpOb}<RT(kq6^A+sx>{p;$>Dkfqe?o0X6
z0gJI`&7LlHl=QY8Xh3JCvtl@U<-l{rpC(y2j=QRjuBnp5y?jeNSPn1mySaO%5W)sr
zhqPWYBV$fw<)6w5=yt;-)Jp<9zIh+%pgteO7J8P~$%W)u56`TcjYSS2R12V1-nYA4
zmsH{%>F2#>G3w<xQFZoY-fgI&wc<`Xc(o;SCALvFa4l1gGGQYWKc6Bz9{|PTO9`6V
z9(-YR>TY9))4S0RoZ5Tt=wM=GOTzQit$<s0y$S~b`R;|!3;6k@oxSe$6a26Shdwyy
zQpU`}DwCLxJz4Tx%JuugB_95wv4Y~g??<HpXOsl)k*?50@puwQo@IZ$51l`&dzrdy
zX@>a$CSEo7uKbk25I`;#FzndAa|3mSdsMHbw+HgI?Y6hCk!Brn%k+>B>7&e!Ks4}T
zvXIwa-&mzdR0s{Ys6}1yAIy^EL~^r!CF0C^@w0$`%ZB}~e!+4SSbcu?wB{I^d%!E8
zQjDp2AJZ+-%S2WWXHTr#uig^yn>gd@kl(`UeLjK|!8nGVN2o6GdQ_StchPe=CJB(F
zwlP%BsT&nN5e>f-UTX68fryyVOq-@F4+p7A)jS=hhz<*o9TGd2Wpj&__qQ9`A|CK|
z<h3x&dW)OzIs8t)%+25pS-h5kCH&4pI7zch>8ALH%}Hv`oFy?{<Z43l{Zg*)%_QT%
zLDvhHj={>XjGN-mj1gets2#6qcO?gh=Dgu7q6VoUV7-0asHULuMr5+#;`)xglsNcE
ziW{`nCX8s!Mg*Lmc5*uL^OXk}vPTN8$Z?z*_Q4qYgYa34M2}skDgn6Di|2c7ybtoC
z!c{p8=Y2g$Q&b$C1oL$L4rq#pAoQ^W%QsYSpNRB~r`>FQH(t|r)Ow&EvRC1OvjnCL
z<SXS6__(}UXZqb!@<77h$A^yVTB3KJ`eTpN!n5MEwDiFujZlKaPne8M-uIt4&LbUt
zHw6a=wxiDp56EA8xrORtV8f@8sm0q=R6i<`y$a{DeT^GT33V78<+aL$k8i&LvwoH)
z_0CjngV#qqiz$eyl?ZmRu!=5fkQou=-jI3X8Bc2VAkHRkq=jSQEJTZ+$!>Jcq_y|-
zDvUCP^_!BDkkI}R@pD{f;5_2Ew>h>p-|6MTfYT1;Mf(>IR`>io$K-@L$m<?1)5~LC
z?AezoKQ^ei$2iY1$`p6wXL3oYF7!#5XpTlzPN**93mWW}p<WaQ&!ye7zL@y&kjW@v
z)Kf5SFcWoNsX<w;%1U!8xkJ5wN4I=-%I7D*`1$yLLX+0WJ!Rks9GSQEV?SK-9_{C*
zhNA+UBDrO7$`L|8?^8C#ofiNN;8`hhP>57-yJMePe7<|ULiIkQqY<AO6Q$1UZg$p|
z=4kFf5X^JR`2l+a!#2FES4Y0uvwjtYMok-z`BA%W2=l(<?Kdw17Te(?&^z-7o<*R{
ztDBSOSJpxxZ+4c$nBRbMA&#&0s*S%i=Oz;goZvp0%s|pL_$qqDV6Rq7r@Mc&Zc0>O
zj!Z&CPA-KuGr7-y=;3lttsjJbl;RfA!bvs4kn?*W_`m1VR6Hw;|C(WY0FSK!CWk~}
zi>6QEKZr$<)O<OQAC!EAN7?EuJ)+h?oODt-6(PBed**M{s1w=Yb*^#G4Ae#|u^2t1
zxlcMfr$oc?)nXZ+A^J%Cxw%0o(qs8;3h!;2Pp2v0c!GrFM_))kPZYrCQC!BwsTHl=
zm!e*SeCw|&zIa(-g<4(fMZr?{ic1abiT!<DiZ(rWsi&Mc;`>$7QO_XZ>_JPyeWJ6`
zKOk`tl7^-8ej1HqBF$BN&9<(mO(w6rxkT&?HE_4K77}4Y_XMwI%g-ycz-!O|nWYJS
z=X;Z>(h4WbCYf>uYgQWEn8IzC`o4Yv!d&*l0V;_xG~|iAm?!I6LYy9WeYymh+By|Y
zQAHw>Gtc95X7z!@ZBSARD|lCJ|MAYrM8|J#1C}A`X3FU`VvTpoq1I#O^$EK<j)}??
z$=WBPl-U}yFv`}Nr0<G96$!#7s;eQSBo!v{HRM#<8Npg&w+|<OSl^g&VEUHXh`jPg
z3eir*KyS*N^+1ie#Jjjx=9?hN$tXW2xcoc+W$29u^j~6%ztcytidf2Z;w40wJ*W+K
zxXEDsb+WqBXrkhS=!gA<a;H*c?^?;0m#t;t<Jo=}3bL9C|5p}h;7)a#M+NmVU*7W;
zS`o4|RS3(+ERWtf@EhC!kU*SB;@wZ<r#hC}AGnZ|HJ&4&YgQ}JJiB{5D?eVSqrMom
zGU%y?x)0|eeATOP@Mv_6XRA`cZ$V|^qosVdbv%sxmhytWQL;I~hUqL;ACTI`6$y*m
z7e~<@#NxqU>!dXiX+)z^r9%NikNc!E3v4`B&VQtj6iDfiVPX-=(n9(^*4MQ9Q<|?*
z&_F34l3~=Sw_{kH@6p)^sO!d9@JHlFe0ar%t}n68&dx$z$TLw*O@%BR>@u%>F*l_T
z2V<?toLwKN;zvg4l3WT!D}`FX-<WB6*DmU&zc3jy&A?xPL1q_QtrE1Al89l8<@?p+
z<qZv^j(Rd(9}@9VAyl<FC)HcdR%8JK)4kr!8DTYeeI==?D?YJr?IDC;yvTT>j&IU6
z;Db_zzU)S6McQCYWl#`C$4LJ%#hCDt(bks6nr2sqq3|cxjoHBeewxA7yY$mtX~vPc
zBMm=J(TB{CnfH_(+bJmVgZe;FH3nU>7<2_iH8(`tvn7^?jP!BcBnuXA30b~!w?d$|
z@(og~4o)3~eEQVQ(x0h5Z%Q%Fm#^G~XSk<%()FP%xaE%0bbw}!HL$QrP(IZgQ7PJ|
zRnWTWu@PH98gA4ZbD{zxH=Xaz<;&UV!LqaTl+y@&wQo{r_k!?R2FdXJSZyW#=8y0|
z@FF#y0S{|^B!P1JAaim}kG#)9PNtuZa@$=-@FE5)jfY*{!`#i#00xg+P7Yr#g)!UC
zkxcOI$VWOZk1>35Fi&1l6^X_t^1er$-dpm70oyxA%hhv3+`vXNq^Yjf-!SPy!uI~W
z;)I9Tsv*B)1^Jx3+;kwz<02dme4#d1$mdEEQ|^y#xavzvA_GN!cPDOa5*N)V@s3|8
zUMyOAG#Ho!97mEfSg&4rDIGixz~4(kJ26wbmHER71LwD+lzyUpx`;PVv<<uqJRq1;
zUCTYog-BKiY6=Y$pn*eTG_aD1bdC1@?v6S6;sl-nAJ3y?p-7AgkJZ5dW0-?j6_SqP
zm{Bu#|2!yXZ%nT`LjLIIF+f$!7}RK>4~$sFoE*z~k8`{6@jn|+O91|!`G|LATM}CF
z2ism@jlHI=6h|a4G1n`&EHcuOL5l3iK(!z7S0nFfl_z~$Krtsn&_ECAb#Iu`xLjRh
z(TuG)PvQW3%+PC*__z9=WVbz6lqoDiZ>cTGu*-{sLz<W)%UT4ji=zD5J-ByOEgD~I
z9QmD%`eul$_3~hyJIn96a_dr3et}8XDkm3gXphXaAMzok5LzgviE@yGvA*#|)I+d!
ztc2V3h$$jrVQPnUdjT3EF%Cw-t<tok$3<rJu_$^B3B7+_|5pug4CJ5v4?GkTlc=sp
zrqPUJ2U&QqH&hD1Gm+}j92<0^T@P%e#|3sTcgT)q(SW~hCH*W?W)#U<W@0<M5$DT@
z%utxhO<pk(Eyn7s$6HnUuE6!ene@SGJ3gPJN{RRPd4$j5q;72|UooUam93j~E-UX*
z;iavHGl`A^ZI8U;>!!5WM^5|fbX2qR{hxoBcvJ8x7_C+<2h1N}(&eA6QSUH!el(B>
zh4*9DiwV3)4?F#<L?>(*KPVZ*o#c+AkIg)YS-<#8C&6XtkCgu@0kq2#Ekn0{9`5Xs
z_uH!-&YH={g`gTS{<ohF>RvM<un?i3-!wWYXa9Nsd#6KgiW0S256%@y`uKonVvlfM
zVB%*aSqoeW7Qx{}3^Ruoz<6A@&Qg<FLS>+DpLwkG>wa5o(~Nb=@Cb{jhIOM3%f;1O
zPj~!Tgwrx!E!ZdjO)@buDf}1Fq=7eJz<Gi>XdQI39T|iM&M{6A114j*RgF2b`p?Dr
zpZCckDX4CdEvzWbG5x0Wf{(?fLn?xS`QwE2KT{O0GRfyyFeIxvibD$O|A$f*T0L4S
zq^2GlKnNGr_J$huyFUEkH;BWvh36`Ko`BE~DCD3#oD2fRU7ta!#fn;_N;&>O@)hPM
zy;*+Q_f}f%Zpi&^vo?G-LD(>g5j8WE|C8aTPRIbIQ|9;k7+=s-vi%+uS-airG%!l(
zq2C&?&Sb{1qr(Irb?oV;&@@c{xoP90f6J){xyapGago83yA2QAA8CP|hL0RyRE9Qn
z51ypwjZEjXkPi}7y?9GOrRC!?NwZw;Ce$AGdPzMbuCs15WX?YQi4A+ez^Q6XTn$aT
zSNyz}mzbny7OC4+F5xsEec~K%nC?WagH}JROY6bq+hX?gQ4`ag<cI`@{=0a3Pnad{
zx#RE<&KM*iEryrJYL>7501YMkx~F6>j4%6s9!M-@H3nRm?yr+J;Pq+90ytxzuPBP#
z*=4%ij;)7VWkX9oH%Y&>_RTjoOBdhM3>a{G-aliQLCP7Ebtvtu9RCgbeTyL+s{GZm
z)jNN<qWpOC3$=4wC4+vW9@dw;D~uPA#kYiaoC$8mld9t~3D3nJidN8p>Fw^JIL1X(
zIl4wvq{>^U1=Kw}t7IizX{9|a^f}v`I}Z$_TuZQ?#w|e^?C+V89y5A%QKck+b}wto
zzk7LIrDK^SChjRKl084nzu;5U@&=glqSSCGuLL&Q8)W?8%x<#8Ggy4}IdZnMN|%k)
zBv(49>=*d*6M(@|bem3z1fp;CZK~<^_-~GBP~D(#Q;`EeuM#3D6}Q|lnDnn`@fw8&
zj-l682945Kb<^WNAXuQFX^*B*Q|9V7|6-kI>3<xshoXPWNC&x0MBW$19F**dlF~Lf
zU`NGeZ&Dj2&u-{4J$m0$wd?V`Dp77b_U_6hclpSkRrBal##;YUrxm_OZAWw;Z6Cv%
zsaNZ?i*aV16$>Z(^7E+rS^#qP!okDFPjgueWP*35UOmoFwKWIHuYb~iE0!fsfmPs{
z=F87L62aY;-sJquHrGf)_G(lJ-BqqWa9G0f!P~aPu@QLEtP&l)s?XZMe6N&}4HM~Z
zYu3-bkc`th7%D7AOo5D{*S_0K@7CpHdu;G!8#wMf2<j8!-hVb66R6Cj8(>M>9V5=n
z<Zb-PCpIZ>LyxX9ywZuirEYAs?K=Gp;g&?6s(7PE66ALSd&E|DPae5MVlS8%j=YfE
zRyT1SSW;^*9$0KuwkM7epPuG^mHHJ@N%%Id-wN09h=NaMegBD@fsuZA%YBR~*%j$K
zmWXmNRv*u4yI}n;)8lSD?0jG7v%Z89v!JQ$ZjgVJ?J~{heHl{$$TFi1Ha%fp#srsj
zy{|3f0}&wst2l6Hw-PPmf#MJZ^ep}=ozY(oZ2b&A|6j^Krn!ELdVYrx{)=Z)^wj?W
Dff9e(

literal 0
HcmV?d00001

diff --git a/analysis/Discrete_Coal_Pic.jpg b/analysis/Discrete_Coal_Pic.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..a6371beb3ea98d65d3f458c8f8bbe4f9614ed857
GIT binary patch
literal 55338
zcmdpf1zZ(*zxU8c3IYOB0t$jicZ-xL-6<j>Al+R8N{65{3J8+Y-5t{9p}RYzdFP<+
zn)|qW@AKSyU;G?6pR)&;8D{?T`{tewo{j_OZ;ObF05C8Bpb7p1PA33AVO#x&0C4v%
zKmz~(GJpVs4ZwlVV8DL>3?YDU_8b6i!(95~xhxFz=l8$@fVVyX|M@+t;8*YiTMhig
z^v74Yc$hzbVm$1Zcf)kY!+m)Uiw6C^(-nY2`jLsHiNzxmQ!?fo3;@SX@w@QQpALSU
zJ->4H^n5R|--kdT2Y$iz+{@*?o6e_I02VSV8ca7F3^@Rc1p|i#b6N+Gf;)f!bM|re
zZQv&uSU7kDL?q;MD5&5Y^3MaXFmQ0N@NftS&|d@N1papb9t#2cGNT~k1sQE5ax)wz
zx8PUE6gRW0ab>&KDVgt^yPreByNFLfc;zY;^)(t6RyKAHPA;Ka!Xlz#;<x4G6%>_}
zRaA8z=swhYq;FthX=QC=YiIBA#M8_B=`){@(C1;{FJ4AOzD`I?dXt>;HZ><VFTbF$
zsJNu2wywURv8lPGyQjCWe_-$fWMXn^dS-TReqnKAb8CBN_v7CF0d&8><N0Z~z|Wud
z?7RKK0{06R9v%)J3A$e}ur|;g$AU+=%!r6BD1)SJc7dG94H@TV@T=_Va}><7>$vyL
zyHM~bSthP*KzHqI&wg*m-2c&@{jp<z?$;0?1R#Qc90YhoMDS09gouoC9tGvxITRcW
zOtkYCaq#gk;^5&C5?>)DyhK5ShevjujN&TQHJWP#q;w2))C^asuTh_ZKLiE|2?+%m
z<pK)I1!_V(Lh8T$aXJE?mL#VT01XZXJelCI0Dj=$jSqD^>{kys7VE4_$8UKQGHm%L
zPU2E*jZY*{n_0!i^<@-%FMG*~wqaV4hROI=jS;*y+BmwCdjn}>kOH&m)rmaJs=U&Z
zJ=94gKNhWL*|u<#%DE+FrB7RRQXtfV47H>$0Yr=jxEWdGc|mxGd;mvIru!D=CE9EA
z^3#N7WU%5@IDkx~@VzQrGIu8AERKOA3p?$Nv|2vYf<ljFn_cr!DeecAs5koH;%9k6
z%$;MZc*%X3nEQPI>emUew()+=Ja35%`(j+-HdtV;;G*Y`ly{rt$6gEN7`fy{nshhv
zfB*XuST%PKx@$hH(r4Z)D3mI{-VLvC5Lyo?zka(A`6e#4>S6a<hKV*3qJ*8Ia7Ukb
z+C7iX&Wm_HG@Nrj$Zer5!42gqg8=^p0c;uRA&sGttA{F+5zPKwMcN1QAHYBJ|J3`a
z<1tS=e>l}Jf#ZCa>2Q6QT|(be;FgqP;|iw1o^l!HWPi|l+OA;p(PHvS+9^PpD{%_E
zYtuR)xBs}T7eWm;7kExTmU=I2!&h>np)!7W&qc6nf>C8i!~IG`#bfMrn5C<egN-df
zb8|)(p)LB6_>{k>%w9~ETJ`u|u-u)V$x4RH9Q%>#p8F-C?#}B3HIwCerd&$R2jc6o
zg&PSCmFvaE!6tOZM;SEeB6>a}+Gu;#O<F}14Yu65?&2&k@U=A+@O$OsS$hhamoE%0
z+H?&rT1#sjzorYM)DNj8o3kJhua*{}HJn_T4RzXRJtUnBTRBelP;92`31PEzR|y!-
z^Aw~CnX)lTWkRqR!?KD)Zr#9~IO2otp`*()qE~)0yl8Fxnu9Vd>FE>fdnPog;?2BH
zX?>1HQ~n5imxuTGCyjlRh2!&R*<w;OUxjeKbc+sQcw){RH?h<dkKzoojvgL1Wr%Jv
zq2-g2xc(~eVZs~V9`j7QYIQxNN)9rQfOELbwsR9F9u|j^^Q!Wl!}+n<dS)(`2(4bs
z7fkEY^k|rn-~#%h9Nz3H)=p7*U{*^!yHXmmWSxxbx5iqcwyIJ8G8RT-ADD8h;^@%W
z&{NG>&sw=Q9%G~N9I`-hC-P!3b-hA1$4Yax<<TV;n6*f!XjSVtBfZBcLqm1hBIlY>
z;$fPXs?w~a?vxkyZXc>GUAxs^C_}m4R9u^Y-crjK?^7&*B&ZJa(qQwrWc(af^CKL;
zbhBLuH#UR#qQx?T%}({G;IYpddegN>>38JfO7phd?XMAYE^-oYm>L|s9GIQ?D8lk^
zRLJT9ywW44%!IXrG<%PQp)91xmB%*}cWWSv{KY(m)CUY5JDFF{Nv04I)0NGddC)j^
zHSm`yn>Fe#8W!gkN~np^y@nS>xo)@(z#OGSDRoA3g~c)Q|FbPy7ff4Ja5mPm)DbZn
z9#Mh8g`lOSOdKne3vRCrc)YhRxJc`uZDd6{9LxIL3IBl$*H(;{vJWwhLli&Lj69(0
zW<xW#=Zk@0Q$bvkQ`)scdMDKBbpwTz0s!OAjTP|xZrZ#gX=KOCo_ShVW^_%hLUkB&
zPCD4&umHXx(fcCjcqN@&Gmi$4?2vY8%h<`1;_r91i?-(NQoRli;#}misYc9ttFT|}
z)a6f581!ZP4o^}?;}5Kts%J~pP{YdK6_f<7N)fj3`7s|GK64j(JB!X(Hg|U5;!ifT
z8;M`M&75r?p)(>Yc#1_tot@PQAZ@-j(z<E7sOVUCj9`kWoz}X$<Vq=eC0>>(yK(5z
zvRF0M+m3g*V}EYg>ug3z=hdF}bMuLob1BpRWW(fA5xpf$IHT{i#IcwKuHg}xeXx&Q
zeZ11$-;Mgn4XehN@E8DGXs-5cQ(+|OwMlg<0S^hwDvgU^Z<(tiP|y-N6$I4zVcSO^
z!e|_u+61}WFRe94QKh|rGq<|JGal<BDrp}e>UiGZAWSRj#Rz^^lpT6s5SB)KOG#km
zV#Xq7vT+%SVg9>fwINi)l#}U=^o03Z{{EtuSGeLISr%+XG@!lR62<@ZA+xS{GO?e&
z|4a&wbyLsBQIOuB)s9tSmI-2_tqW%!R3D2HATkT!8JA6Q0GR0;?;Ufe;!xH1C}fI!
zY#*G|=@Vq!sEMPKj($2Efb|@kv+5oUan%rdYsVVNLkGg{+ALkCokE>Uk|o8RR=9mE
zRENAr{0y{J&oS;8BeGwc7wkx>1_IN{_}|v2&vRKP8>7@a7!kL?eb5lXoYSJf<d)4Y
z3KjP;d?|-YH%zTb4KsE9n4`I7j%tf_$HhapN2`|bvHkVZWimIq$sGk@6qs!?j4<0I
zydEFU&2zmn&uJ=eecWQVM;S+}ReG^90Eh)h2qfz21Qv7+2D?ki#oeQ=Xxhod-o#v1
zvNaZz5U76@eLRw!&5Psb)*;hV!5QCeAG`APr8KNVv&iA1U?T0rZ8g57k>S}gI;f^0
z`jN_CmsEs+w7qOjsc*{nc2Q9-w))e7*${-`Y>y*hZ`{#6>&2EWk-X45p>27}{VY4u
z06yOnoumr}2N%3&b7xr&{$d*|S^r#x@bj31H&=L{%_ly<un&m&;adOSyn)l~iUPvc
z9yE)5r$B093SUggx?N^}&5X%oy3B;b(1x7A1hi8?b@E=&qMgSb6V74iRrTmUvQGZ?
zrafPDd(}PEu7}T*{GLeXu-==<TJue56t`u5o;!r>mg(z;N&$eh<ypSdx!Nt2tqWbk
zxJZJhT|bN^D1@0Ww9^FM8HD2vyrVV28m`Ns=AhuXCt^pjk+v-(fS|qPjN8dN;J=fB
zDbY(Q_{5>~*7KBh0R%6B-HaA$xT{i_^PBqX%HoL&E4pxNJ|h`HjhD>fXr_QoJCZkV
z`dN&aRSWPS?NPcZL#rb*8Ui`-d`0(!AIS(15R!Vc7$Gu6UJcexrUmd-hZ5DYm+~z3
zdP^#N-xw-Rj;ZuK-a$1rCogf+6>u`G^67r<N|s;f(+{C_u$9Fpmwtw0fGdsR3R49K
zNq+4E&{EXLZPFT&?#!`quD$Co*V!xY<DpGkyH}hWCL+mnn+OGH-P4NO6U?hU1vL1#
zb#%5fbBgIbGn1<h?7a}4x8TDy@fq%97P0_}qQ3moiv9dc?T;<Sv@`4*Hyl!Jm6S%S
zE%a#oq8~ixz<~)~?`C;XrE1^g7FPGljfkje|NXs1MfM8V-Z;5#+~m%O{0nTw9bv1y
zNUv5pH#CrH2IR(_c!gkFbkSg@Nm&mD6lf^3b8e?L)uN>{A;+UcN+ubmDL2~akcC_%
z!s;cx$xUids;T!@d2|m|)=5fCF;4L1DUb-6<5McyCh;y5s<#OzeDAsYGV;kqMqed@
z(YA{w*&RXxgn(H<Z~jm4z<>7eyjZW~U-!`WJ}GgTZ)jZ5&^6aW)sS-A-a$jX_M8DJ
z$Z=pv!}t*lEowu;3H7lF|JD@yiL0S?@hL#6E@U)bZY&Xz;K|oWvl{H~_nd=(cQIju
zg60&cd697B=kjFZ`YC{6IXEr7@h(C|Ey#j)W0f|h&vIzZTQzMFLeAUhoqW&Oh)F7R
zmA4$z*xP5G4|QZn{fM`xqQeYrh){Y|%?ko<K|t8hT_THHT`j5eweSp0<ay7FAND}z
z&=^LC(C;SeC#E0jAFC(7UM#S0kHcM!V+cQo1aE~@#ODpkatYOxN?j-%!^l_c+rDLY
z#LEn4p&#NHXSOo$V9BHON@dUYc4?QAqS}(}GIDdI67t*Sx@~^dh?|f``S<Szli3#f
zx(zqBG(~bnGChffXiWuyX{Bfb>x;2)y>SX%UCEuZWFxA2d`gExLp?pOCQDbC<$HFN
z(I6T-8Ovnj%_cNW#RqR$wfnA!h)RQ7j};{R(0Vu{aZ#}(|A|dF)49t*rlQhLn6J{p
zW=mrH=UjKF`YbK=w>hJZX?i#}D<LFft>911;8)JB=(!;6KFOa`bhun?9VZm)<x0<E
z(A2a@pyrN@L|!%}_pWT6B=I0Kz@N{X5pnU;caO%#VB*e)4_JQ`ujsdZb7#uZ7~lE%
zJ}#w*I1Uq*J8~1bE6oB|SZ0Pc@r<9l*aPS(Ei`X1W0#q%V>iw#zN0}*=!Eg`6u6z&
z?wGafUAQ@Ko4F*_f88Xnu9vJRa5<CjkU=B4bM93<I>6Crpy~x#9bI6Dbl@~c8~fK{
zcE5QK;08~2?BnJK9&AMseH=8@%s$dwjK%Hy_)h&pO9B<a6$}JG_lUmZhzz2A3S`C|
zFlHYpsi?pB;Yw74f@L4zwEKtE@oTQ&%X7H+Oqe%>b#cbUunG_9U>RB(VxAH|BMKek
zt5R}NC!Qu9%v>LK4|nqw;v~J%GZ^b5A~2V-Fm@4~HvOTuaQ{iRO{dhl>q<{y6lMJP
zPVjb+-td!f4H#*g7Z%_R$Q(CMsfc!Yxn!6zh*|MezhRGFWzdz1WFZ~B&M}p<!#ez7
zC@CWmuM3vK6OZ)(08Zq;48n20xn<2>=6o8^=^`L{hnRcZL#a3QXt-kGQk~KCp~6$;
ztn;wx370mG=oT^WIP+L#LcYB-`atiCuGL9}u^pQca<l$OAAK=6mdo#2W3tBjz7e@W
z6@Y4JiIpZ+Wr_$KY}|r;8ti;?8bTVH$<q?haJ7RB+blI`8gRlVN6(YTiiA^|yd=Ni
zNk+njHkYB=&qvm5PUTz6X2jJ{8oOz#$U;(4yA+tDYb&L+#d6O?W-0(&@7$l<+MCM~
zFoQf*;k}*pzI8P!<F|+Se=sZepx>7&RClAv#dE{g)sd{js0$EZ1E;piOpLE*w!jzq
z%TmG=WEQ=aYwqY;vT0l9hlwUwg`q2om+zMepyc?OG_Dt>@G|ag^p_g`CcKzzhX?QJ
zqhyXnI?4m{0NlO|Th=QMXG6aG&V5y7wMwBY17u?9bya>$*~tO|%VN{9`;GT#i^__<
zx*-z25qKOr8g0v_Ko(P%m4m*~8{4T!TI4sj6Gg!&Z-%GQw&5|gZj<%x&{VS1ml*Jf
zj7sGzd-`kuG%bvkq&1`Q+l^C=sw?ULB@D^W8pn8zB66R(UX$!aLu6iL%oc$Wd*Ptb
z&~vZrM~N8VV2?dwZWaw7^?Ow)Z*+KOI%1^^3-e&R+~@J@f8iz-{#5dou_$uOIihg}
zAYk0cYi90aAUlu0zV<+N#lcXv?;DNbM@k!ukrz;nZj!#l!B-LBWGuoT&hR^LvB)4m
zv3Lq#E}9vtV#K7Xy;U|bN+=fQ^+>tLuez^OTO$0??}l)ochLLN>F^B^?LxWghtA$`
z90KO{_K^p<qCJz{mQi$Ey;vrCuxfRf9~Eu}cFNV8kC?CO1c{XCOF%{kx@RI1#9u!M
zQ<&2BQY3&^YP6x#Pc^uBG5eNt?UlL^Re@y<?LZE7xCT?04o&E84$ir0UEf=(8Q?>%
zEQaEjS`fdOz{b&?Kf^CCIOe<ntYJ7*>Md4^JG)^jva+x9)rK7KAy$3u9x0h`U$K8!
zo7i_#ccUiS{~i27AC;;+ES^wL*NY5<uOc5c#UuH*9dpc~M#r@c^qrtb-3|Q;!0`6T
z+;nw9hD%bU4a=MKGe4)`rz@GZE#<;&a7zU%<QCeKdr-=fu`90>x!*9#!IZoxIO!Sw
zDzdj{6^~^jI8&DkQS=u{H1TI^OcW!~!QpvXe6)_O+m_`Z0@s_VNOY>Vt4e>5vA{@v
z$<nYD14WS{9)CUUzE#=67V{-%5_1dkJ1-)LJIlrPLYwftg^4k|5;jXV1#$KH^&+GI
zaG2UrO8GqyZULoL%>4c}dl5iK<KTPwN~iFw*hTXO%ncBSYpi)2r)RMkinBEHv<my{
zK<NvQc@eHJSFw)NU#!~Bh8|@ehq$Acf@C`xAL&udKpZY~i8VnF<KJDj|1O|`UbMV8
zGSZf$r#-)QUxW=^PfD4KrKVt#q7fzB(Yps6-353s7$M=W3Q7t~?(XK+GwEGQFs<EA
zSzF6^Ob2oJFy(&vNR|O^K>#l`V`rrx&RDNJmgBN!76rh6Y*jmEJGmrZY)Fr+zo~im
z<J~sjCtA+-&ir9)W7Sd?)V+>|sZMVP(gk;K#nv4-+!$+cC3cvKDwmx|KO#+TwYzF+
zOsP(OI6g(Xe%*<N<B0I#Nb25{ykXpRlSj2-&ye0Y@xCDSw`O}mc-JeuTcueK_X%w{
zBko|rdM0k<$H}WHvdnZ@UMpj>#{S$cft$p3@yK*^NgL%5$in(5pl2)Fb+T-j+rhg>
zd<wjf3)N_-)JJ&X!OCr@98OK19VUY$087}Hljq*h&+Q-<PE7$LgJTM$lF(a8p)4K_
z-grDC6Y3Y@4eQEua&hkDyb;4GFd%h|KX^==o+&Y5_;#Hy_LZA1lNYdw8F_4&@gBg5
z^Q05y*$^ZMv6JN(NaF^OCLh66nXSlZ4A<%#U8!&LC1T8+wdR45E=azAeWvtK9o-0j
zc+(p#-bvVQiFK{+u5{u<H7)d^%C-?gYjDh-=Gqc~`}%ddu4MkHBX3RH%toy|_txsj
znEZ4J94zIQ6pwfqbb}A7f;xu8<%!@+t|>kofCA>ll4MUSr{#-m<|)k$;zDi$1TLss
zrvNhBOP8hO6TipKVf{*1*7JM!NTsB!aw3EPJEZr{`Teu5*zd{=_$LGV!8h|UyUZgN
z8pdtx)Ul1C_D>f_a>y;4XD;AVHw(Fij;8WsGi4oJv*&MW-Mh)WekgJ9{uD^iES21;
zIl5R{P`u}m>=ZjmP{8u4tPnkKKWl)Pwzw%-Y0CMHE%&h|`mRtKX}!+LH6tC~^x2h@
z>I|PBkQ^?Oc;I?J>x-&L5XuEVNo}!1Vq&x1&eT%s$*LeW2r)Py=Q$ybsO(PPbg=@_
zT&3}mlEQla^F0bF=@ykcue{E7s<)Z2#V>^Vfwnsc=02u9*>~*WG~>|KxEhmKpv!wv
zA-6woFUxX2T$yrSRNBE~-YIp*PBAN4%aOZ}RJjXlZvChYy)RaL)2T83<I0%N+VYW5
zlE!><YQ}MyH5!v{SPe~_2-T~Dp67JwR=2dXC5?L05(lPc%lu8j+1N_1hxIPMg6sRF
zG1bu;47Bjgbl<qtP{nGbSLy9TzA7XjS_)QEJ_WKkvpbx;@|Q^KSKm%(XgMIfk4|MW
z3y~g|MpUY0QQ<*KzV4ImM+t%n(U&#v^0`qvU>2oP=jVIoEzTQfFvrmzzgk&X=-g1N
z`~l0!`-D(X`Vb&j!#}_hz_pCIJ)rn5V6#Dxk_iKO5%bkHt9jXgT=t(pg4reMGrHLa
z;gVJI5*|C`zQK8Q(Oo4{r<XUt-`LrS=t`uMPbXw#V`3S#qUZCVon2(YdqM7@D+={%
zSN_x9ua+fd5z!0)M@IOk58{7D!`4@)E$Tqo&dZT_-cAj%w=Sk8Cd4;LC0J`OTO1qk
zUaZc_yi$=i=`W(T8-XY9DK3vvOQzTJ!ECR$^R|^h=Yr83vmb%nYyji&WNiR9y`0AT
zm=aZppwW494jCP-cmQ+e*2)<@Gi($VQ|_N5s1I*?GXaJ|I^J2KF2XB|ekc<NCJJ{G
zYiT4kO30>wk4^p-WBAn7a~_sTae$0=Pi#W)7(8J`0(2OTIOq>Ha?z!+6MWF15IRE{
zVCrO8+FAduYuWz^EqYo_AR)a#B1pK8vpU_UwtTz#_^35Z8YUcJI4y4SI>#rHa7{GN
ziRtdHFTOvoJb;M2e0Il*e|OaH*bpbE4tIz38Q+LGd}^C2ympvB_jzt?PO6IE<J-I0
zpAI$erk&5$BtuXUVV~Ml>m_-4WV@u1|1FDol3!PSxpxck<)}qbbr|2HP8ZvEQvt2L
z0~UJ5;RQizw~zyHzMqKQi{KdgW)7{Yuqi@1uMsFut>jr#-V5!W8$?-K`r}^OrTj9P
z-$s89ep{ByqSyS-dud$80^Y!`EY)+oJGv@-vSM>c{iWsTMlVEj*v_Z3WucTkDu)dY
zK3WQL|M!3h;rsWY#BVfslV9kNY5GSX;%{sj=6k)i+hVg7cQ%xsUfFn)aW&boLGA9!
z?j`G)%f9^A;DU=?;S1HH^#GWiS6Sga788bX@CNKyQt)wuffTGs<8m{vF9cr-u3>T)
z<v;@{iY{pv+Z|86@1Rgswpuo+C|JI~Q22rnciaArSIv;L_bdV;p7a^xPOB62-C+K9
z=*wUp*2_Se;TTNVk*0e0I1sIn@dln+QBaxU_(pMDGSP%Z=yr2fIuAL-<%q)mu_{m0
zGTpd@ZMi(TClFAl1J11BOV1n0{2z6Pra5RyL}?(HqBHbU9$GbFr!hXl@=|13{)CZs
zVM|Z?e#vKVbbxSu1v;dPK`7v&z>`%Y9-iyh@QsN-#;!dptr0KG%<Lg$F;>=lO4MQ<
zBv8JhEni>h3!C-VzZyXm@hCf8e-Qxg9ek3n?qV2nztFuE$&~pj*bC)d6L|%1iNjqb
z>Ed#q>&mS>AaZEaxO`{A^xR<`ApWrmR6eo;VL^m6;JM0!nWA2AhU$-%sFzm;>_2)k
zEi!HDSKRpN?=YcD@-WfPZiOc2b&|M7#!9)L6W8LLV6T#VXVmh?^`ojxf#(>Gg>aXi
zZ<EdS-g6c08zt?X?6)N5mMcWKUtSO}1GAj@%6ls-IJCX*x?ZTLk6Oz4`)WlLTjZ`J
zsLYmARDN6gkq4F4x%xS_I(GwwUpPHEf2B;cB$Zkk^A%v3T2EaNE#oCflggB=a6nZV
zQ3?u9zXH*}kUj#Ccts`mW<A6g9r8hy?%%j-h+{nv5kLQlwJK{A>0YK#$uwH2@|gsw
z37sNxl8eAb59vwj4&nG3*-5>P*ECL^;?-*%BM`glGY*9)VDAfpW=vq)yp%>pUMd)7
zMe5S0)=MfmV#_mN_#?OTyLX7s{K>QV*FGNcr+uI9dGE?(CSxwt6lAOKdCR&ta7K$>
zd}$XPs4LCnf-;hYwsBk#@z%NP1TS~vmi4B-%mQ|YKJfq}Z^z3I`24iQsmUa*D+0u0
zbyCk_3HBVr$JT>2a-_`9`5@~siM3$Txpg{(Fna`qWsPUyJvw3TqUdwvwJlV9Tz?Ko
z6PXH0PYx0h)r=qT+dF4eK?eRB+6^>c+|{sGYDu0^=ZV>){afs@Hrrjsg&_GAh^Uu0
zL?@rFe|&N5dLzRz8~J!3CjaRD%c(WQcJBMGybWKualSBYXJ2!8R89f%g6va3(?Tcq
z5Z$n!-*-f=rlevQHHuUI!4&1Is(7?Y&;3*2$)!8|J3Xg>2pj+HCEAcZ)kR0l<dl2h
zZ05@2F<Ys+WFdg?z&7_p*huH3Y_^eimvA86<;U5T6rDoc5g|GKDS`TJe5~rW!03LA
zBpruwI`JLe%c@R@Iwzrq(fgWHo%SRnn7EyebXI&wJ-rOq-$+-H2ySU*SjoPYxI~$n
z?#IwlOXdT-J)>g%Y^*^l)~uRXaGN4~H!XkGKFoj3HJ_??+f#p=An@RBrp1gypx_y2
zY^R+*sHOibMbH92-z)~?Vh@S7SzzLk@K+>$=5P6rAE19a+*_!h0)d%i1?!rHSK7L2
z;&j4n#HvMZ6QM<>9Er75uFuS$0u)pe*A3}wGgyfbMJ~`@n<HyEa-yY@1(|2O#Dmz7
z3_cIWQ$Rc)q>9bvd7+Gc_a42Ck;VE)$wS^cre%F0JL60GTUwc1N=8swMGPeL5tJWb
zNhUWO-&t2WD7&;=NnQ~8R&}wTab&QG(z6ld@jC+=A(;aai)_CXZDzc3wVmEveTC3=
zzpue27h!1GqO+Iy%ACIn!EaFZ<=Ed*)^X&CPn*fTZZZfGz)K@&P|4=^rBjq;2Z^J3
zxO0J+09~!#uHbLe<ESAbISG`hnSKmHXdXMbGm#ZBPU}nZp-N}4;{8dOUg#n8s2$o&
zmfwUdj;?$~qQ_Ah(lCjAy7^FyW?q=t*1JK4mX_F)7v6TJd-Owf52jZ57;uO;b;8Wo
z%%bVrURv|EPgo*6!$oc-xy^i1pdxi#x}k9#_tY}2oPbrKvM`|9>rPIQu-rUU^)UC$
z4@mK2yY{}s$};lEaUeJoJgFs9?H|5PUVKN8SGuG&C|Dgbst|0kBWo~7w#^zbWM=r^
zr{>8B<G`Hue_A=qPJzie?<4f|dB*0Aa~b#Y+@5(Y8U|kJib@edqT)C#^Kfe<vPZBb
ztv5pHsRLp1(Z!^_1SVU%?ft~U0zF!y0h%@yDPwo$Ms=9ZPIz!FsrF~DD)D}jS9xAp
zWkG49%!0*&F|5_^;`OSd5yogYw+2Moygj~A#DdN%-4(<ds>2^3T;<izS~+dzSL@6r
zJB}=2T!A#52_@I;#wRRNNAOzf8Ce*%N;lV=qHSHob9x+)2%O&|2RbtsaxD2p0>iI;
zsDYn7KpTv{S2op9ra+jAuK<(|zK?WYl=9<Wn@#>eITW`3WpM`CsJN4W!gC#Ge$R5j
zg}}yu{2^uaj1|xPUKexQewl6UxFz*F*fG*9Txxc0Q;72W34XsdmLtA2#C9r~KBw?9
z54o0>(3;Ge(p<{%Tek@!9on>GsnmV@YXoqy3CVxBP5#c75lb=X`_FMhan~QpgbkkQ
z6;LtRfSai${8kpnZYBiRY@f4h$DQ_DRf5IjrAXq0DJC(4MSnTX=IkX2{TnwPC{Os?
zFlSf6mbD7^tt1u1baTFb=krj^nDXY^Tuzx-*Qx%etWIZ&MC9))5<j&S{har=qW!<c
z`-9%GUt-0RBJQs^;?J7YC&vr&=Ytqr_S=dC!6%M5RuJlH`@@{|<tYB(h`^pE8fYH9
zwgF$k+?Uix{~;-F)VbF9D3M`qhRmFu0SI3T^&1EHQp-yVv-3h`7jO9>4_K5&1GLkI
zdvbRTKdh`|hSpBb3L4t4ulQZ#L}}wtBxt0HhUIpvz!i<$<cZ$3BEBLSI$#~Xd+}~*
zQD@dvAInud$P<1BYK9zFa~)I0Vm18F!}kev_|7BdmM@fjH##J0?&E+YUFDp=S?SFZ
zHHJS-7w9eDPC6a3fjDo}yRH)Vs)`C3LrO|cOL6|dkxMaDH(Ax4E`v<M*6n05Mdkp7
zQIcK#b<<@{mZTLE!v&YRLQ4JXk3(om7^>0Eczeq<d|5SLTLB+CLglUsZNg9}Z(rS^
zx}8?~(ugjxs5n^1i*>lY5kT3ZyQowU4zUrN^{KfQaT%%^watjw!v(6=pliu}fN~B0
zzSCUkRtOs?CNn@OB~6_KeMPMF0k++7VIF@<u|`vW@tScGE}c`LPczKU0zE13#YLv;
zrVXVLgX1uox~~^^yNB*<G_ojS-Px!$Bgeg4+(a;sOj7Y2G$nFD)y~weD4Q=npnu_{
z>L>TU9Xs36R;&OJ+|y_&J|VAl)0hKk{mg9#(H_bzIy^=~`xF2TYtYXU^!^D-&*m1}
z%Zn$_EE3?~OG9&~of=Bc>$5_Q8hl;vm)m?L5v{%(r~M)jS0(@-=LqKt9OAFj&IA(M
zYED%hc(y5>FX|!`ageCO{H%MmMq9X6K~F3^nOOtT7qQ)GG&J&IkkO&R#-7YfOJZTG
zyXQ_iTrk^E)~~>YUUCEH@6Wo|8dUb;!iaf0<esewo>;ECO47Ca15GUtdzXzVyD1_U
zaW3hWdZN|?Ja~Ew4N|S21)&$;2tt22+LSXDdwdw0uE>bUFMe^f5o)K4?=A&hQM}(t
zftId6=r$PK%)13zRzKj@Xnygmf^IR;v-*egvXO;QN^S&t_WX_O``q8qfc|tm4Sc@s
zhfbN)s$~&H8AkBjXQ8`M!_PV8pb(5Fe0|h<i~KPek={^-{fMnzeL3%u)rL#0CR}An
z<d^^+eVRt&g~#mIiNoahhAfuxV~fb(+NZk6*l+gn>vLx@#~&ZP;$O?kTRAoyamL<b
zS6YwBKU^YVt=Ldks-Y?jS+YCtiXW1c=;(#UrL!7P9$&O~Wqs8UGdm9TL!yHzF+smz
zPs~IiiAV9;qD8wSAd!V>8`wyKXxMCGxq>o#F*`&kJb+Uns+DT!(UHwP6kxXn<UXPN
z=MRsj4ZR~DNaL~H`v8Y!)-2_d{|iz>mPk%6D3Q<_iiP7-T#`9|4@NjF5sxi<2y!c1
z5Fh1ff-8c=|0kjcid4-D^0d)0Rb@5a`kgks{)LZ8kpe!IvbOQ(xbK}<q0bECDTZih
z+zu8bT(!d=hf%mZkQe38rFLj^yJfP4>`rJcQ^^fpF%I0>_JLF2f|YqN!Q&%ZL#nn~
zuV@fqop|h|`I%HSl2G9;XCNC=<s3V~K)`9&I9-vMF4@yaYgRjG0NzFG7->l}Dfxq7
z@`rjOUF=$TiF_<uKz3o5?Z7EfvA{H4af4c3n1~vea;qNx&Q{Y|%7U8oS@Jb4;dk2a
zMq9R0*4RsrF9uUD{9JmP{w?WACC2r=!^Q>qYXso~wbpB`miirgr@-ZbSlX`i^zrO`
zA4vBVFXV>aTlj#Cu0Dc)y7g}Gb!idU*VN-MvK>6qA`fEa16UW^t6=sb0aGFec_kYL
zh?=;W1$iTyVR1|#AW!#p(3jKrWB^W0gq`PAEj6i(rv2Bya*JIgrm}G#qpTD7`sJ^-
zBzs_mENA^;vXLQnJ-i}QD%m_Bq*|aDqzq{^(Fb3;M?)R};0paHrij*xb!)-sFOefE
zVSB-wT2Rcfdi;fv&yw^<IVA;j%h04zpfQ20pev0&)Jdg7JOVnY)S<R6<YlMSS=#rk
zw3o(U45Z#Jp!Cc$lWs?;&$WCc@#x$ItyG(}D(^lSyBH0oxIehB>~hXOf*sDL?8-Bb
zw;eRn)eaJUm36I;1?rj0U4+npp1C)VhanH*Xou}Fg@#^Dmdkzv7ST>{_I1P)9?gKB
zyJ|%O#6Q4d66MO3+>;WVZ4Jc^m7dA+o>}Z_vhjfUk4T+|3TvB&%{v-whZFs&kcy{)
z{N8WA@#SsP@?+^M+0!hxI++l=5isALz~y}z>%Al?U69?!#BgSlQ~U-orkx?i)+M#3
zKM~`sXvRb9(F0E$83EGI$s0J<S<;IYqy~IfO1qDQCW!CoJJWT;IfXwY5SKm$TG`u0
zIC+9ZDbEW|SxlOET1j1;XWJ#Df41hxZiU<3K47iwTk^z;F(%U+{vEG<mXxH_q}~%N
z`U|#kx8HSFT-CBQKECcUO={3lGL6U?WzZ5QdCb*`h}G-zh_fr*CEJnN`1Z<4C-yb6
zQ^0{lf9?b(&|vih%NF4XJ#maJ$r%z`o+`DQA+inHSO#I{nTBJ<b(WFfP5~75hErf$
z<w$&A>J(_k?RO;Tdecc?Ci-UdNT6iHX_+47J#j7H(uKozK7CSz^6ZZeok(xHguy7;
zPac-IRV=m2P<znr^!j)a0_C1;xNM3LhriZBCdd1B7|(9ac_K)k>|nZ|RT}G`c!8Qu
ztme2XvT?ZNKe9RNwxj^9u&W9kpz#4^&5_$m*CqkN7jgAa*iR%j6@Y`crXZ$<rb?jZ
zg5<x-AngBZW7Qky!bCN#4`Q*>DvV(d_QffpNmW<3Sqn5sng(yagOt&1yxv9;5#7qX
z`FAD<&^>Zt8@Gz!H`KT99e$wtmIPG{?i`zl6;$Q~^Tg-q&@^v`h;s1dkIWlBF5wM+
z=c%1<t1u-zkK-3Or~U-{Z=hA^hV5n|%o3JCorl5?uR?V<F4!}K-_8xAQ9skh7T&F{
z-jrQSxT>e)ieMM+n>Hk2#GPm%tc}waG%V&R_s~3n7p5`eONxV)gu;Xhl76$^))iM5
ze=KmqMA}|wV=BYC|C?lv-od~HG;g^!LLIvfvTuYqGd);P&ouc;&-}txTLk;;v)6MS
zf#52!kAARKzv$Qt=}c6TX&SDARm<XN2jCV<he{o2#f^=V#Nz#oZ7f_79Q|F;><!26
zniyF|!VWjlnP?JcKWV?2i-GIPF}(>@m?q@yvM5?Q^kll_7kKTH%6*jZ(g=C_pOP2J
z#D&gaBXFV&QkagL#WhWls@mRkpp8RaU{_T~HPNv!#96OH8c`B}0;?zFMw;C^3Bq)}
z6lZ}Be3&T<eo4+xfe!B{YEpcJJ+Hf_$8=0!)#?{heQ*$opElmMO}><L;ktX6a*b?i
zIeleZ)6mLaM7EC8{?Hr0_wG!^g6L`ZlTtehRcgt<RciY`SUON>qtqI#r#1M3GA3zz
zF)}4r?7f~#?KV)uKOyG7J>kQ=lR*MGX}PuO-s+b|)grzKvy-_dgx8wz;S#dqC8cL1
zl$vg0(2M_|Ju~>xM|M1v0X6AavFABdr3_OfeWczz7B)%xyU#LI2rNg^aM%lKDF7X$
z>;~Sz{IjfyZ0;<tT#_l>FSNOmV+zH0eJf=vT#MXK*77!!J3X3?3`74gLvb>+gNZ57
zTnpxscGukTr2f9wG8NGSOiabi=>N4s|K*In;(x_IL>|MBLBR{^DZ3h`gDqiZguIza
zsefN&-~AegEX>p!4EuqvEM@=cz?ls(i$2SpkkqorXA;ET>P!Z!#9`FneSYnbNpU(d
z>l^Rhm=+=Vjl;#B#i+77U@KJ4?3NcmF_DMK<vvULw|2|*Pq0%W;w_1y!|n99*Sb}q
zVxke->eS83yK+J<xr^sxs<L2Z23kQrLy|$TBMTmuFLV7QnIu1c0g_Zii@uX*L!JE2
z3JqyV>$<ItSjC2CQ}+61M4@BID0Lg#qV4sL&T+j@bcVtLp`+?&I%5S&XDp|_5}l_$
zGh2!X61GH_XRJ8g(0c7d`}_dzbCv~-o?t=gpu=k|-O^vuh*l<;*~IRmax_r@w@2sU
z1uafZN4mY#F|54(KwsU+na3JRs#wH9Hze=#Sn;|L9;B90O~D#gFzsc$vW1iA&&)HT
z+n(!POD>?~@huOchriGgz#JGyiH!g#oC2PD8xAKvFHVF^)-!v?Kh{PZp<gX@NL!*$
zH!P`p%Q{kfqswz74=qIM`Uf5`mEr@2J|^fI(xy%Ud|Qnpy)4M^{K>6?+=x5w%1L)W
z=niudz>&GVx5r#=O%R?~;L%CaYPO*Ap`)kpm9=@%CVmk^JA}veW=^hE^HhOAJf_vp
zIA7>%>4A7IO&RGf1JQ$X6}756fs-HQ*E6%q-Dles2soOlO{w81CKC1OUGsK=;b@-;
z!E$g0Z$(!3rkgk-yyZ(yq?L+qbLlp+_+#=b!Tcc;lPWK+0~%a|y$C>6MF5ut>hH-8
zSOxt(H#jI#_?}qv^;ws7$Yz(H=`}$=>NRWES)))_7+H~T=-_2Cm%qF*rUM2ISb!<w
zw+5FCOC1u&40+_Pnw7Al@!yVO0c4$V@`*qI@=$z6Y)UhI86bTYj|9`0FKqJ-)w7Op
z6B6?@H+*fylr>MsqJ@abZ3lNQeV`JH^=5E{kt4hDDZp&0KXl|0wDcZ8`|9<hLJ=rv
z0)2gx<%BKmsF##8(kGSg`5E03bVb0keg~cL4GBNb_&9}y{o!~WW=<9e`6W&Z4}-|F
zdSj_!B?tLqgzRh(!V+I7$tFsRl*WpHwm31Xm8(;d;U`0~IYl{7FHo+xW%>rnVj$S%
z=AQy}BkGAFYWygW<}HQ^&~A4M#2*_>U0Ta-d{*y3dy}wvxl&8Ox^gcn4cD(#e?m=Z
z|K$VLYElsI2_6hMWqeRuI6>`&ZfV9O=rOzq_2vDD+=S}e(<Y3xED#X86$qe0>9^SJ
z6riV|ok%nUPl&x~>G>K3T?#J}9a!!SWv*X~lY&3^`1)=?16jq?9EMt2km6!~_rF8%
z%<_7aV~Z)Nf~$frhD<>g$mnTt&>mqn9<`89%tV@R>pXr6Bz-8rrz+|&u>f;5KRefj
z#j!$!Ke?RaF}Xp3SwgZ;+gEwr=qS>$?VQ&2n26OxD=o7gNI2<5Lq6hQqcW{zZs9Hy
zqbjmvGglWhQUl*#TiuLIQuQHk0&esxQ)Y*mjvB)rrolBh`2!M(i~WRoB+fBbL<?j3
zGtTjw{DhrZd}Zqj{%zjU>*24j-*~<;EnY(I#IbN<oAP?YI(K40aj5>lf`=1QSTKIx
zscBxS2MOZJOBlMIN+_A3ou;-Aaq6{5c1jLA0Pmg#(>=oA^)xLuMGUP?AOviG;ja6^
zo%ve=F8u^@N9_vz68nk2{R;O0>bmACWNu<LX4m1V?qIjq<}h2&Y^Azo<{Wv9mj6M-
z_}Y#w`v;8IJ~SAiUxGJffs<=qWpboGc~XRXi1_rZMoa{B#@F@ky`DwntO*e%3+EMm
zY>v!Vti|hSJTS4Y)DVbd8+i&yKkOmPx4mIGe`Sb4aSM2URlus1p5O!{x)bq59UXoB
z3GSZ1(bO99yl?t(p%flhY#af1EgG&>BY#Z-8EJ`=2&*&F{3xZ9UVSX%hG*ce?kVu}
zimUK_srNWC*A#~No?zJRi?1nVobWiZ#ga1EhEpFaO6Wg7AxpC2!`Nu_c6@YWsFo|*
z-Zc+GMV(p7u4gg|m)Rpg4BDp!5PtFilm5S>iQqo-f%@i6!9+<iuP8;u4RB~A{2M`?
zGspMe%!y)4Kn1ZB;!nJFo8?`a2FmbB2-b|cSa+TO7QC!!B3$DgqMa<6KLOgWOlt{*
z%qEvha$srm2$W<#{wntmMHleN^5_`hqZS-cbTReUK_}wd9WTM46VS$L=R9Hx|8zPJ
zMMRO~(1`dMKIIN4=;}K7)1`x^O3<`_mONlOGZTGr>F6!l3+Ep{e>Q^<N}wq9JcP6Q
z>)}?>Q<LS>aN8vG*?-vdZ30WB3hF<s``&-Z9WqCOyJ_NMweP5xCZ4L59o}Phh=$oB
zgy$@qaOR~uOMmp)6k6{ph76B@OdtO6>KuTwPk4gAJx~!j3Rz%5C=OB`E#uZX7w!3w
zh!!wwo(4os*=tg^G)2Q=+Y=JscH8erIN)T`DPYay3gXuI9mP?HB3-@&Z%hE^cOr^^
z&{K2agZiG9SV{M|ngeQ-a<bRuE`?60r>=J5J^9B5cJkKR)Sk4X5KZb4^paTAi}R+;
zPmrKVcFn-#+9Z88RjU;8(RQ3dU;YsV6n;S|T<Q_(hzB^S4cK7{3r<HT!Y6&dCKe8r
z{^(oyzSf|Thp-~K7;32fJA&@tl$ynLpo0qBS<Fn>F5p>`a$U%<+D)FcXN1bZ75983
zaa?T|JYw?~zsj<I*QoJzmi^z7$NtPQ`~!04AB8^t?h^11!ZUvt>0Hvk+`HYR`>1cc
z_}PQOafZuYphMaJ=k*$26x?s~Z=epY<Jad1gB(4&XTErsX5Ae2!c1Ii5Wkc`!keg=
z+byF@%3;QtC^O>HSfWPpXlNff;7`Klm0o<!5`*mRtm=&{Bf7lR;P<#s<d}q@N?}Bc
z-!>9S-d955Z+2Zg_yp1I5Wh(&`Ep=`B^<WjlyKmYMz&8EBe8DUK&vkNcJ!?zcJ+<e
z<K7KbUPI&D7bEg-1Fn=zL%S|gvy4}kjVOytg2~-kzg{h3|HWz%`gyfj?H;XJD6}EO
zJ+goE<a-+~Jkf#QdE|r2+Jk7-evf2>@I}Q-LyuNvr0y-^WA;}lug?ZMOkMKKSznCk
zI}xXsHuf5kLh#+@AgF;0?dcupi{y*L3r+l>-yiSmz=1tV=2GH=<LvGy=WF;}1)jFy
zITW$c@4ZNQ?K8u>^(9DDRyoC&KVFFeQ5m&Os|{B@$uPfL5%=hELQ$;6z@y=0)FQvB
z0ax*Qdrw=FI*BEEOy@2V%>FWwW<+1M2)dZ^-ZD{y_Sq}GB!~)O7DIdy@o1QcizRFo
zJUYa4)$cJHWjl>1itaKmW0lvnFSez)SRzv*G+ZpVSE@ELsnMiZ-#&QZk;l!8Cq(9T
z(vK?dHWem{atMt0HzIY$dXri3o$Je6G6^KbwL<kh6iY{_h-{dcl_BA4=Ny*2RtYm#
z*k3_*bZ3Lamw*s`XQfQB8=>YdrNN1g-$6L~DF7FbY4Tq{YV;G%E)$47n8%i<K$m}D
zcS*~11zt%U+`=M0o*eZ(OtSr~{*vKNnBNL@Xw@8nk$Lntqz4ONYUUDYl2usqkG{p?
zcxZI%y;6Ba-n#+^-=&vy89i$fEN-~Yi5mjg*40Ogms*RiXHb@ugh5ZiO~xwNRG}<=
zc0$X-T}W71{d0i`Dlc72b^9Pw+V?S8Cs!9|$eBLuDfUZ6c%yTI7AjD2Qo<i|3gD5X
z{h)fjUrP7j0PYM6DW^YFlR|-1w$)zFL_VoJ!0azwF>V=?TGvc6(!;Jd!SzdhS`Hu^
zDoJG*_N<WS=X%~Re0K9S^@Knw$qg$lRD;8p)NCOs9nh#Z=j{VJ4|K0xr5I7W^09q8
zf=bi|NAx!&r$G7Kz{f0~pF2ZUV%+|FkL3SWS0<=FhJW812aI#<rev;a=M{uZ&I$+M
z4_I)M04Zv6*?8t7$C-!)e#Uk}UebK8%XWSqh=ck1<?#8QxQ!E`EFPr<6_sioL)`pm
z2{bJKl&1{9>-K&SU$3m_D98ID+>)t+4gYEHoA?mYR*K1u@tnMn#YK6j-KGh&puEK>
zPaM|o9vd=HRKk;ieGp3!)@8ORSe3<g29uvZZPlAlQ{j6^#4>l}vQ3(yuH&@JvB8n0
z`$$NI_r=mY#+g=nKRd#Rht>dO1JJyRC=O16d4AQk{8Cyc{%ZR$a7%H^Kz)yV8@QaI
z!PhV*mFLsnj_Wge1r(q?VN-02`yA=gR`10_ZYS_S%Rgpl=L^u-E`89nBA)%|GBf?~
zh@N1b?WM!X$2hIgieDu79$G`O0O^<5QZeZ>X<L!0+6zKc1Mf752QBQ#?(lO^kJ<7u
z`<LxxmVm`@1)lnoE`cYg-A6JJJUj=E-D*S3VF89=53eT|3KC3Z0`rgH0(`6Eau<i6
zsnw%_*Eo;rPQlb74*Qbk$Mrd-b*GaMa2bJBAC%BeMiB98)4o3}PSDgYi9fwV7Ak}3
zsR?rmyb`^23RDS$+SAbboTnu{+OtOQ?I=3-j_Mc7oliebaa#lwHOHg>%65R^Ui43_
zed>f2Z&?33T>W1p2tIYj^gmn6p!yQO!MLV&S~KWPp;Wwo4=R7fe=;3Ty~7>WRR20O
z$pI)uJSr*AE9l`dgU$gwrckxS*c{Ow7AY&~AIK%Y3uq(<K51ds#vg(q94TlBC+q54
zFodJAd_CS+7C*Q8OfGppaF)S|<;)6q{3C>e7B}b*hHwVuW>S9=!fET!frfChzz~l4
zN4(Uhn;+tRCPbX2o{M_~`Ma~Dxiz>@_&MoMoLkOxq)*h)H>E50e4j2?dwxE39(w23
zUSAD%m^^%?Gclv@F>$3#=H+MWS^Ie+XlUG?FAgK@%Z*=l6igUe$x|OAC|+w=xK&G*
zq;OL%wN&+R|H$_^?boi<-9JGjwDRif`h);mRwpTQXeoin#w9zxzWaEzbwC<*!;uTB
zy+;peAhBV=k7Z*&*GR69LwoHeky^c2{heI_{;)4)S5zzNs344l=Gc7N?m%s+mhflI
z4$;=IrM|5{_rmD~KpmfY0hgYgPc>>6=nxRt*?#@G>)}np?^S606J!omC4U1kd@nWq
zBgF92-m9WqH?Zce@G0pH{QmEHNAgYS8TzK787o>w2B-_H^HUzsY!GFY{>In$7b1YR
zEl)!eqNYDW0qvo9;5B)RCW<+XU><|_(U&5Ulz;1A{0^)B8F5he5*JK^%`<%JSPZCA
zGW&&%P;!B&OIetik{$c9{D<h24gs^JTXdg1zTvFp<eyZs<};OyX#8Nn=9RDR=b+`y
ztj|Hq>C2xF>GZeq-RD#D+sZ|L)|lr1q<;FE#qM#ZuxSE6RSPA=9P6FPMAa03a06f7
z`jcHHzKh+)c|~}m4YU2%CZ1xJWV*Q?8_wY4Io<EKcwNuhdmAj2z@f(hzJ&8SinG#=
z$RA^R<qhW9u@4Sa1qm(U;!wj|^H?VPn3c#dHC<EiyF!!J4XzFT6~9L$e2*^w8IkZ)
zt)0z!!B8I|PvGf<@ONtil$rGSlD8_0n$oJ}?jYe!Uzt1wmU)@#;WI9%p2O9ZZEVMQ
zaPtYAU!Y|E66?em9?R}VM|}tS-UaERNK8yrq(0-RC(G=@8o@f0*7zL-Q)BRe8AQvC
zKA$=(AmZ@m;zfpqu5-qTUiu^)7ozdBE%Gp*Z$<8aoU<loUm2aD%hil2!!f*CXL1KS
zW1Haf8Xq94qD@A%Z}{eyhTQke#*Yz)CJXDOP=@rENL>5F&qpZ)Z)A_Rq7N0CD{n<_
zsk#91QMH;h8*A=5rQu6rUB>a@O-&N`;SU~Oi;XrSx@V66F&@+CSAZq)=$gUzZMzt7
zfd)RS-`L^DQ9SdxvGnN1AmfIQ8s?sJ=)x`1cBo`IlN)$aFhTZ=M@}Bj_Y|HPydnXD
zuU_9b_J<+_u#2GHn(H_bPZO<oCF8zwyibF|En;x;S8VA%8y(5NjE*U`|G3mG%yUZS
z6u@yVjN9LcIt2t+OC#rwoXijk-FJ#_JPg1&BuXn0L1B^}f6ym7DFO<9-t?mT1TM%T
z&ZPP-DwW%Y1Nx&yTJ9gM2lz29*K&u2ZZte+51%NAV?CHkvf;#8JTmF+Kypc|u4uBN
znvjtB&ZXv&_`l7j_IKdn56vDN7`&;p-}=Q}Y4CcNg`*2y4BavkF2zbuJjYDQ;(b`Y
z_be61uhnp()Mh^_L2?Y*G4Tla`sJAt=<pu2)PXzSuV|_asW;LGjfpx2<Y)l(_Z1p%
z^}eeZ`S)&!-X9y1n}62jq1^wu%OmOkGra7AmH=&lNc17VzfB+|bl8Zrq--DK`Db(7
zQ+CCOj<A!?-R<pp+l;XvHR7*yv-I|w%(J*1)~Iq5rI#Qy0LkZ-sT8?r&HiV@oym6E
zow=4%;E{vN4R4>-yyN=~hj;GYM16YFGPRbA0CB>wCj>g-{H&8wj51eNL!Ht$?nE@q
z4>KIoJnn=M#*s$)9P5H3FYB4?s9jpSFrsh&(KDa5+ffCy5f|Pm7RPmKGZ{c_UKTlr
z&zWK38_VWBaSFH_?-Xkv-Wh*t6wgZ;mYYXve2%ZG5AzA_RC78#EB#YB$TRw~Vx#N#
z0=0A1i4x2P!ltf^V2CfnOh@W@Rd5XM<!bfQ8L?2a(&>{0Ol~*%y`3QwexM#rWYrQJ
z++FT5LFyXlmKQZvsKu6?xhNNvA+nRh^MM;?naPePM4Y1uuC9V$&!KS)iv?XyIe#xA
zU4IQ>yFVE6&JN2b^<C?}tE-L9YP3$decbZe6t8IA$;I*>$6LyK-h#z!?0(aD#$t*l
zyA_IZi*mg4aQ%<K@;d6DJh=X<%k|fzw=Xz{H9#ltpUWg_U#6RyKc*Y9zcDbL={Ns#
z1LN;=v{cY`1c@LEEh^jOI_$j;#TW9G55s%x4$(_DgKFcRyPX-?G9{T^frHfD)zOA-
ziKw^O)D#-+t6<w&aMr|K)-1ikf)XiA-w1?^rP>09#A#ZF31vN^SX1Ox#61NP!=_D(
z551OQTv8ekwLy(DQ%#_J=Sww%^G8+jH7zjB{#6Y#kMDWJ$+*LiD($(*`KoqWuE>tu
zu7$P}Ut1;?f;kvyIU(vcanvy_;3Jmfa0!v9yY%W)hZzt}#>vZj*GVm5VSPdMZg`9V
zazP$gP>gi=$yUQLV(v@kY1{oGzdS4H6u@3H7+2L!lYq8Wi%?djgZV*0Fw6ng6Q9AM
zZ<qQ%qCr1{L+O_d7ah*feM3IfqIvsUizee^me20CgwO7_BJ_V;{Qb+sW8}Y<0sqU$
z@~8dU(9`C3`?Y_Scx2xB8Lr64nco$u!S(8!Q1ZVB)`|1Sn4kZ7MdW9J67vO`;djN~
zkRrU72Hj^R^UoM)p-}DKr-Bq>*UZVk5pq=hR0MwA@PAlIHC0iSu&^Q*hBqu8X$NpO
zS&tYOvYc4U^Sd9RT9D<J+gdkV_!7<w0z1LSP!is+kw`jl#Wu?&?ktAyUQ<_1f3&4a
z_K6)4!<nYN=U@9RJ2C=gNBqy&5mGA+IOsLUc>N$heQa?3y#@A!Zs{{<fXY~Jp^-y<
z`s86n*1E6$T&=1lw4B5UT27(>7V6bz<1a|zOj*9x3zAKJb*)eH#aPF*(e7V#i9&uI
zkK_bR+Pp1)*9-8|Rr=?(|4;&r@@WnQ<<n1<PSE_F_2>MZ8HEHk<1iK`DyXCeo&pN?
z{{It&@z9Qau#Bi~TSi#0)8+J2&I=RcI8AwIPp_u1=W?iI-ay7|$(;v}AxCk>2^(2P
zUXSLyIJJXGuHi5-1qs12A@kNh*sH89Afa8K*?pE)qnlHDHvpzu7DEW`gBuWNN?V)G
zLh4tEHcn{bIV8DZ;yt-<lu&wE*NcC3yd>D~9?*t#k4qXUdk?NUQqRrO*PmO2{p#rC
z%<5Q25VwhG{J1N}zObFi@{vyEUd)D?^1D+Y=)r-XnJD@32X&UDmn7aP{OepD6;c5e
zJZ{^it#tG?st?8}KQb0=p^%m#5)bhU6wgaCYu2vMRomOR9$hnfwLqgv{mky670g_4
zcJa`XE3YeDK48W2i{h*QU}xMPK`xH(f?TR1&>&Yzf}4-p(2@kIUQtD{J7k1!J68-#
zX;BIXhUzhC`>V22_m*ynw70V0?X!b|2%t&i+6^^JmT)n8tD>%Q>s*NKLnYyZAT!O!
zD=r)CN*}PoJ8BBDKnq<aCI~7)F2AVZaM^6qOtYoaE60I;ppzIXL)GgrL2MQKP974V
zn~H*%y3;v-QPt@=p;NSCzsz-pTct7=T@h}?>6j!#gbnpj>P@gcZ^_iS$x87^-|FCb
z2L0dbxV!B|nul$xAzngX&f*H#a;WUG!wv|;9r?y=UxkZmU`^gGQJG)W-{s#yZQ6p=
zs}!57d337}MMO7$eWdj0eOIyxR32fC82p6|v6R@k8<if>z`8&Cy0owk&1(!-c-FM<
zvrqCoKXe>FaAfhnvW^#<L#As+eNPmLiO=G{;N+u*ee;W+j9-_V|EwqDFWTGwPfAgs
zP3L!|D1u+Ao~)n)bVcj25tzrzhzp-)ttPHGQ$1ZSs8lE0UQjB?@zi$bKBs_X51j8G
z@$@5vM)J%k|6hCO9gp??|Nl!Wg-F@EB$ZLg-YO(AlD!f_S2nrGR>%zD*eQEv@0A&u
z6|%P|*)#Ney{`*(ijL0bobT`Z{r-M`^iSv9x?S^qy<gAg<8i;sHG@hXG3Gw4rG)`3
z^?h(}Qx;>91<o5R9{5_3rLQJs8}IqH*K(+9Uo{OW?69yqTT{CvMp`%Z@tO#l&x6`p
z%Nphqs;jTwCY|L$CVmrtO#EgR@5Dz(hCHF_H~`z;c)Qda-|SfOuCG^nnILof46yWX
zW`LiL?HB^y1%-cin0W}7fxhKrqHidsx4!6CHlx>}<{2ZJj0lls$Ywj$-iM=QCq(6J
zWJ!F|o|n*TzCV2bbajQkywr57PqJSA^(8-%bZ2ET;{7zUQNrffl?}pdc0{SH@2!t5
zIedj39PzO_%b+gaW|p^JAE&@MLRXE={qV#&*jz_>sFuCep>2YI&_eJ()&6YXHoCRh
zAZhy6Aff#&70c(V2$EDk9+`oI)_RGr64cKJwB5%#MFxn?Q#T+!LWcP$Nh!=m$8A}s
znvXDVN;uM(eCpMjbIS##o_%W_$2jrw$Z4328_;3ku+wg!W`~$Ht!WFOxKB2jTtHq~
z8FytUxPGr`T;u0;>634^c2bVe*tu0+ji`kwGY%R%w}7#;9cKQ1s)k|-xj>n;usPR}
zd<OwVgZ_)t0=OigwJt%#$sYu!&SXrX>%?y~o3)Q92<}lXZ#)whDL|73GfS*gmXPtI
zmhon}GuiT`z57?Coxdur_x|gpot-FY-=W5L3S36%<2kG&H)SLnWO<6TpGo5=R9%t~
zjz4dGM36}WjJdSH-xSDU_7wEubStWUSn!S!EpFDti%qbHl|j8H$QlHlYIJuU4*f1%
z_X|<XuETY|apTJ6zk6TRP7kfIb<Q3B<@C{7tdq=15vNa>I(<~_ac110Kqh~qKss-L
zTOR*O5VN@`v;F+zqAZYeD;^5IQE35fP&{-IiCT;Z0;ok)JTx#+O0lC3G?~Gnrl`AN
zkGG6E&7H#{F<Eo1t!-lc)27EedJ2_q9gu)}an8oCMgbDt@S*fbp*KBmnB=XEfT+RW
z2#A(#1w@}Bc_YH#^F~JC<MlL>E>~ObW`2|sDhfSqNNAUz0McIJ(9XAxI$G}Blkt?+
z;&r>@V*vKfj)c8iWd=uo0kHS&l5L3wDaW&NeD6Pk5cLfJQ9sW*!<VSbVmTS**fFHL
z2e4HT)3>ODYUDItzG&ba2lc<B5c%E9V}>=_MXl<~M}u;fx&lwkbzYf^8v1dUX7@N@
zox(WDj!x#E=7?h>d|<DQdjFaPBV6!JmY6yHp^0s)NF>`;k>qp%Sw#XxR`~>ua2sS5
zM2d2X>T8n&t)8d3a<b;~3=~=v-qQ-8U}I{U&-CQZpU5F>@;(Eh1A*1R53(BE98I?J
zzVmGy@gqABLaB5>*A)-l`a<XNngc0K+NvVG^K$bTZ*tUdSv%VQBLEYU7__^T#!kAg
ze=)CXb?VvE5B!^AqC)n8>}V=SJ{E*#Rq)&E5jd}IuH0zDjIK;72zuvALU+ui4>#Hh
z#PeB;xXU|K@&hu=0`8sb<C#IvGBLW~75bQ;bjAigsP1T3D1+gWl9CDbJUFWF3j>M#
zI1$Zd)R(1Fa0D=tVWvHJOYZ%EMbBOlc5k+V5N}N0NN=P^<D~V(3s@Px^=m<Jc1v;T
z!=6ESSMu>p7s{n~>Tz|y#Z#Ch6Hna3ocYIM&1joZ5@KL2{Q8gte?KL&7186^wHaIW
zmICz$C3uY1`pHWXlFgrz2bwcQoXP0_xR7q9-`8)%=J4Vsu6i`g%6y0Wf!&}AcS+*I
zte#f*3UgCdXoPbX*u8+?+D4F~xytrPXNoX_JdT3=YEqOI-pDzc5`tIb!TO8k$%^y(
z1eAm2Z`dH0&X&r^;$^E$pXY|kXjn3Mm2_O&^Q7l<fAlH{x<{HT8CZWI@TjwF%M;b8
zY`eCL;qD3VKXO^i=Ucy2H~~a3lzl5ycZ$}pS+XU@^jdShQ%0id=ITR2#Qlc8cpD<=
z12~NpT(UiG`RTpN)i;I3TXKaYZ*LN?#gQyIS=pEpDi*jy6}L=6cu>OFO45FP=9Si&
z!$w*;{#q()H-htgJ}QRg6i4nV;`d(yZJ4kO-<c`dqX1_kLQ=4+7*%xcpHnTa<~;<{
zMXTbz{+8{}rX7DaT^Lx-nw%+`pNM{#(WS0Zrf(2or5JS}ah}c%7MRU_gWbYNN$Hib
zDTDDPz!M?wJ2lUG{gB1gyf;I|S#Oh47zXwhPK&3wd*xC0$%oV~mK<f3A*VD)f<laE
zvdZlu>ONL}3)7WDVY;y{THh}N1__!@^_f%28LJI^a=9Aob~vlkY7w!Do$ZR%fGASA
z7yc)rh%p0_C<2Zz*I1gBE7bI8+>qRk<>pPI$ov`Z`p;isLPMgREzk3E-CBD@Kc9d=
zYust2=&h0KE$LrQ1YZ2SKgZK#Cdocw?V9DBx8*eM^_j?7{X%5L&!#oZWqw^<$?3(3
z`;D$vfq}Ro`u&L(%MqV5Kk+yA&!@yzvOR8J0x+vgPYp-<gZ_9{XfImdeUF-n1E?7-
z#L57f;l7V1D)r>27dYk3HqTd7k4WTOo)nP;PCd}+UkIEBmlffndto1)w<wgX6&LR6
zU0U9xP#UAilNJkyOSHIlkqM3C68|Szrn|jBAUWzN?GdI`Oawe(Y>|!Bv(SRg;(%rV
z3DKEx5kRmAr!QZpUgcg=>?`Or#78){_B4un5oW70bJQ`ujXiOgD**iiud&*)8pYt-
zw7j*bRU%y13jG*QHn|s~0MC^*GTS~~fObCo%l!d+PD4y}>*LcW!v<6@_V%v^EaVeB
z>#k)wCVj%K%(Wx3rh$kjeEw4j=2;$^xTT9^5wJ42VN^j_bJ<~!sHca@s#TY&RfaTu
z&%u4C2>JUXZO`aDb-C8Kj_&q+id-SJ!rP61${Qb~C+<&V_grybQwOZs?mq^}GTvl&
zz!rZ(CT(Dpc9~)P@z?S1WEfAsL%$TA6WRdvi8!h^%t;nIr%-msbo?<OvuQT{4Q6ur
zFPl#M7T+%tptAs<LW}NuS|LQch7$xhs|(0{XZYFhxU2hVI!k@1X%3EMupupboJng4
zA3i_)mN1~DYXfcsU5ri^i@c@wzIOk8@<d7?0O<*VS^X7P0Mb)1O1j`B$+6c0{W%uA
zSl^q9PsX1kL7ixmqxa@UQ}3sh2Xsb2-)TP6g>4`MEI%JN^qtkox?K0-4KT}XzAeb=
zDF0$XCLg6Yt;dC2kmVp3Wa$YHsh&o=8S9%k3zCi_Yv-WEkdN7?CJhBjpMDFJQa*43
zKozh%85tj-SJtD;MWXXQSHp0XQ0%;(d!BTFofpml5h8Oq(RtC(+g9NR6nQ1@Kr2?g
z-SZ89$J&b=-={?Zzo0`i5-eU=eYqv+mZ$+^B&BSHC&QKAb3$~O0!ZrL_gL$l3Uj;1
zTGybIWi!~@0XT8He*$oh9xuK&El~nGFhKXZz})efU&CR2=)BG_wdpJbVQhfbtLZc2
zSlo&ROj=i-dRi3Zl}eAphx9KVm}@N#209jNJ|-<?<u>=bf(JK_3Os#yU5BO;c-nl1
zbu!{V$#O}m#T(x8tc>&Gq|J~fsN#G-^XOjagS&V>k2yH24UD~$d5*DqHa?tlPiETZ
zSTz7e-OOaoAZEe2B~#I&Y0ujC2y^j%(15{5HeiBQ$oCl8L%W{D=D8)a<dbz0S$mZ`
zE(T1Gro5ASK8HUo!G86m;=VHp$j;nXSigJsKBH@A7cxsAc*`wu%yRl8_ezm2U?jZd
zU4>kBi=;6UUe)adqs&6}l8VZO{VKoFTlQO_JN2%@!wENfbabtm&_PwE4(l`tcEo9p
zA@q-pOY<WIMcB$r-dcB$Dd;N^AeycA(a<cn9ovgSxHKT4t%93Mipr>idV2c`e-6N%
zJYlMSaw!qH`S!a}qRALDyFg|!e9Yy~UI2aQk<UPn9FuAE9tX@Sc$rEzQ<V}vb}T1P
z9ku$dbam7@rs#rMN^bfVlo{E|IpXhEhlF1&eDSfp*G9knc9^+o<eh^Fy5fwG)=Kc<
z6BNa;eQ2wl^hBuvYA5w2$!9lrJ>^9?0z8tL;x&1KFPn4=*fZo9MIxSJcx22C(rs$3
z0<YO&>!p22f5fom<)v`llS)wAVGOjqveN!ssm;9gJ--Mcrqe&zLXH6c5^zRr`Q5y9
zS;9BkjJh7UD4j*^hTv^t8*VPJMl_g1x{U&N_w%$hg3iFn8~+!&4W#Dne~uukFm@mx
z!Cu*OyrmFr0P(pY@&3{7_WHZIaBGac&ZB!D*dCx69?i1KLP%B{9zwKysVb&yv!(h#
z&lew4J~Mvn^QreH3dBmB#QRl~EU9{Vn)aL7W!z>;QBtRbS9Od@GfCUZ2i?$+?KIwt
zQyO1+pWV@#E($>}=yIirvA9u(gvq2bT!f=8!$Q-0Xg>^Q$(ZBD$0u@FaTPsMjBc^A
zODv*%g5U!cc9tZ?9lbbEj4GX43vg@xW4&%cuT`C$J4ffPw}*xeS1|1Jhi%d{P{QUP
zbf^NtMIM3L%0kP__n*(E+|16AEMjQcE4$)X2@~t2NFXc7Wn|RK%?sd&^t9mdGQzhC
zj^AGptui*Ft~3ssFvfe-5Ha#X)hbZwYWRNdqu_WiT&_3JR*Yj=slTLgczz(|eDN1@
zM5ugvrbbeD-Td6=Bd=_x*=5g{jexb-`*j^EdKI9Z;&-t{q)NU?n)t*pn4SLUmLqG2
zYa#5k*%vY!z7~3AjsSYZlQYLX_U(xY)=eHb%i>dsl6!r4lKQOuT)UyRl|R&ic{X!m
zVobWLO9*awV(|z>VhsZ%){lU;rl)A*D;8;iu3VL8DJLLn8&kL<Qn=(2Zz;8$-9wsV
z7qtFSq8wCm&L#`C<19uTo<GZUlrhISxp%}a{|(rg^1HO=ul>#4*lzuwr};5e-rfJW
z-GcJ$-FtBcpZhYiKuJnr0Dy64IKdClohZ`J{-&|7Y~av_vf`lSK#(dCg$d{-vHMTj
zp1c#o!SfMWZ5tehJA>yC*ygsiqija*DBt9Y(JoeW->Ak20iX|^<xuNnqw>{tisgYK
zeu<2n7a3k1wX~oGb?5jLEc(qS{kbn!lJ7#`52-k}bNHCSde3rP3tL3+!Ish<efHV)
zV)vWn8b$u&H!MQB*WGLUa-u~{jR}$q{kpC|NH7acdyM2uYth;ZqD&9F4jGaT4=bX{
zt_07HQO*ufMh{;a=|tS=Uf21O%#HhTwKxmoCU=kZrMHE-;@{q+)paVRIWc+;2K`&v
znU6Fa&u8K7hv6&!AO^dU(TI;1MCVWAX_fH;W_!|}eE)vjt6n4E2FddK8zh%6SL;6c
zs$%7m_ORB~rte9dXB}<f9j)LMt97Z*o_bfX+PJ86kE(6iMi*v*{J<Cm=@+s>NS8Hm
z-dkSf%WC_$@K|B#!VPe&yX8nvg4&xUe}All2%=Iy+nc?oI%p$Gm8?ykni-!-_W{1E
zQm<8cx&%GQ=&SQC^TvH|T=DmLnXk5N%jUgUijs~j8~$Q{KKK+oK7b{d%}W-S404*@
z>E}~tX;6_Mk-sW^&-XnrOk-LxCF&%NO{<N~9=r+Jc#+<(#bwr^m6Q;G=$F-3zbq0_
z{^lGzTM*r=MamvX>TE?2;hLFPL}0F6oa0sM5fp9Y7XDiD(tPYGWrgt|;=3{-X5)$1
zZ290d?-0IM-k(6y7-3HHscbv*e2Or{-?6O8p|b_r!0#X6U{<b8XR=MWrr<*p%qB{p
zief!n$R0)+2{a3UkpRmJvUX8;XQ`%M*!dqWs5^A~B{H%3M`fbN<ZY~+_FwsWAwQSh
zX8)~FeV`i*@~nlXF@%<E`ook|4>8^Cy#DrVe`Cmyvlixq!-BQp4|Hh45AZdsiA=kf
z^_<9^CUhTEhY?d%SY2hog&PTo#gWk&h@x|9meDxZP&)`vq`Iww^rkfl+XM6aA5Vwf
zQqSom7bhW6IL9|EZ28$Tg}pBzFC|Q@k`iZHb#lQuw*S?@<r12aTxTI#EW50S;}W(S
zuc-YRaw>7%WYG2Wr(y=UKMMqiHu{BgetC-Lx3(T3-lj1HFAL6bo|q5QtPDKgg@{nk
zRSxMr<wP!Mn$F5GO6L9I0~gr$sqQe}+a>%1NHhNgME^VR4~<`zPN)ig>_gEu2w8^=
zgxRW-jL+N&{G^@LWm}HER&%SdR$z4GCK@BEg{BRTa`)1^RB8aHb`n7Qp1*CY>3D!Q
zmTbeCmDRFjtx}VnlS<g`c=T~n`9qCR6U)-?uzfP^31;C%*Ly0-+@;A)#BM*M_b|;E
zJ~F|1pq++>+Q`(a3(4(+lL2lY5#aXi^Q`<+DQ(0L)Y1K+I(qpQz~aIj$KHOlSF80l
z^52q<IhhX`2OCwPq~qy<P=hDp>6o%NkiC<yuyP`BlnYj$06iek{6hACPE33kbxd?k
z){hDwKsI6v$YH#2RAtGsgC+b&b?$TuLTM#@81gztD77%Aeu|c4V2ILUwvHNnjGM|o
z*cQA4l2Rg$l!CEv#2A9rSMe7r7*}hL!?(|3ZLG9*gN*s{Q}Z`4`rE@Y-LF6p1Kn>B
z13LHz398`6KmM3vh7UCoARIS8G?`ttmFOmqD<k}$pfLK#6I}KPB3XY3dQI^m8!@<>
zS*DTAAicbpE@)_yk3>CEGF;wF>hTJI)Bj!ev;2jYL7(7aI7R-8t=%3UfBeXl?&|zz
zaiwV524c=-cCL#W;@>7S2dwjLnI|pyIo`bYptFBMC9Bypjjf>GI&sObV5MJSzu4}K
zf;9CaJrF_1*lvVj6%?=30A6V|HI<rObt>x)Xwcu_o@{J*{zLPf9x&fspbH>W(PFWt
z<b<uGKD^t9b^rHxZiC5Nt;=IPPd@UMmG?x_?bIo{vSQ-nRMaZ*zF-F)5!FLaGk4|z
zA`p^N3|`ABBN&xt_&)25!Bi}jrQ_=(r~CcB9E!S~RfuVwdC`Q%ghN<mFq|04?5X9c
zzcMr1A*(HV)wkfkg}-@oU5O1Wr9j}!l1S{U?5-=~4@ElJU*n0se2mha{NMC9e`EJ=
z$j%BQ23*_<CBoU8tjYQ$(}WUwUM}Xa%L_1<?aA(PaC+-Jim}|5B=2o4(`+prjrbxm
zJ?dnZ*{YMBZo@V!!TjL1Ynr(>@A_@d&m5&ZeduHZlCC_m1C8!RAh$`{xNt^M$tFX(
z=mvJ_)D5e?{oXJ<-%;LuW_Ej5#RjknCXXN_t4u2`QV_Pf%_)O0tUt7_OS{5M9ok6B
zv1ffPzrvg;siq~JC#?0S(JS&VsDL!lAr?|`+@QU1Uz)uJW&EemefM|`x%Dj3^w-;V
zELN;DPCv}=fb$&>=j5x1<=1gOck%RMxfy4_h>7mbU5<b3fKq8$n69=e{?0?nGCQ0w
z%Pz>Jwhfz?zmZDq%oc}G7{wtB5Sx6I4~~^=;8@v+JXV66ZtmGiE88#~Lk%5@tS7()
zobpzv4^1q#y>o2AN9lb<nsP;?H2%#~1T1fa4Dl?ZyiOrq#iN&_Z3xSp%kR^_-N47w
z8U9B{f;%{+zc>>pA}+{eXMwXUI1?N;15^iP`YR0jnmL)^ug($RiNG);e}-i47QG_6
zbp5UAKXSx@a-myC9G^+PW8-QneEL5eeIxKYE@}3G&U6c_^MRz3H-tXyFISfPo<jc&
zaEyEgKe{p&5+%PfWig6vA69KSC++KkeiG$R=Cr!-_wamwi53O)zTLiTyIR8j&(aX@
zYQN6vsXcWPgmbkt&ZLJ6zV-~~r4K|p)Rdxvkx2XzP=n)%=mr|V$<7i2fEv-}JMiW$
zbAZe#jc}w`kQoxZfm$BqlHncSI%kT;>M;gM3cb5KpUg+!&8Ho3kZW_dt|pr}*4B$N
zFVf0#*>aWt0=l`5uXO}U0$1>b|L~Z+74w}$J3|U5>>zC;x7g-KA8a#+vEUTas1xQG
z>(X*6jS7&rhuxT~mzR8Z^MYPrYtCaQC)&cq<D=tI?_PT=y`CtVhgNF*8znp(_5pnt
zkJLFbge}#bRwcP4CjaQ*X>^)PX){lWd&4@p%`hTVa37JfSJJOtEL+7GxJD<QDRwhQ
z)Pgk4GTo&7TnTZt@mxyafSrfFHAY|2CDj6Hnv0P)A7t~4DfG*)MBIc|Wk&3=mR@q#
zA<P&^CF@4q>-oU&LT>@H^clNAnc|-5wSnrkU|2I)0BSBg1Z<PZa=XQ0{$;`<jDESm
zY4mH?B-aV~3St3_U%h?3%Iz67yJ5p|Z_Zk)xryHTR9nJT$*Ih1Hhq(>fh$Qb4h@t`
zb!Q83cCcr}NH$x|PpWxOsaHQAz<5yH%bT0=j!~N}r*-hzRoDr8OV}s}FOFct1p!JY
zP)mfWi9h~9mqcfs(SdMZo~&TiF~ZzosG$wrD2HQu36H=C4UkCYr0II=G;nLtg`kqo
z^%@`|rMXhChrsy!?qH2)O+V%+dEYsZ_A=`a1)kl?MjEvZ&#oyup#x0-oOWgeovv>A
za-_t?>;tt!I5@_gp_(5hG3vF=`}I$Lyy~F;$w2>jIC1q<Vw}e6$Duhrb8eE}D{k;O
zeEopWdqItVN{nCP!lJ1kkh>O?wTw0PdouG93er9-%Kte~IjT~dvUav8f7KI=(02sp
z+xV!^4pFtIOm#U1K*W6DtKJg4nF>bD^R7TOmWX8N_Vapd0O56JQoe)8d*cR6h&W!p
z>0aj9a4*}xn|F<D`OL&n5qNAr=@lZA7(3S&ee(DFI}_=PD4XYAv3D@sVwsnM5z~Om
z&E*Y#w<q5(B2c@q4MZ8@<@im=9lKR1fVnZ{=RN;(xbngp{L%IpyOD-&fX~Qw#UnqJ
zmf8fvm;t6JRD)%0cY5-jmb`EHcRF~U-&79hep3!Gqc8|OPLR6q+ZHJFAiy9LURpG&
z`r*O7ui<Y28#>4j9iGe-Yg@m=9!M3M3LqD}mrzrIo?>C-es?T6d27FWui`M0YKcAQ
zSMyzH(O*RLqF-Brayv`s-%ax&mvk@NmrzT(!{0CI=H^-FZv(Pr)Yt~uaxY+*e1=YW
z&+!L3hUut1zpFqBs$kVZL>RMgL>Tm@$G}twL4&$ggy8z$54SfwV-3ww(*O`seW=1!
zmTi<;=P5OasCjB&(Zzpi13rjrY_6r?hzOY#u{Q`dUhUiX&Z!y-N9DJ5jZJR3#(H(r
zT;IGCfZc>#*P9$0mKDst9W`7dJGaRGvFShvw^>%EV+ZxmlD)sdlKc@Hx*29@_9V%9
z+^vX4(1D#8_hmhy=jqD5k7XkT%2@tIFS7H!KTDF<@_8V?fWM(Y$B9O@HFpua{t;2+
zxIz8d#iilM<|EtQ^w_dD^Up6|$I8)5V7OuiEI@KX*eK_~^hBr9c;F&?TBSRQ4wSKA
zcnoDM7OD3r7aQTg&_7|t51`_-r?%p>!pXN;9;Jq3!Tp;EM+Eh(eu&eO^=st?`wWj}
zCe*t5(<T6+oI!|%ain@od}PV{`2!Mt@(!my4(Aeg#R5}d8`MEY*vWV!JIGxG&MO_<
zPh{G;oYCij8=gB#9r1EVVaZawv-pQ8OKf~<!ki}`BRz1LCJPcSJf<xaSdfLb<hSw+
z@+>aY(5T8FJIHmQ%V!cHuvO+ew3Z@ut&bnKwAVz{wwv6_iM-41%p&Y53|2}|v8ude
zXE;q$o}GVfc&*Y=xjVX}ZjbY=TCd#EwHeb1HRG=^Y=BpnS$xDlS+ZVh&%0iDcHurZ
z-*Ia-2N<46vVkl*@kx+RdP*R(PEx%dJ9)&e*KW<<g8#reQKvOKSGn??9IGe6AyqR)
zHE0gk*27FxnlPSk=h=b&PR3XQNyyLxYYz`|otv=LYcMtiU)>H)DW1D*=`U(%n%bWH
zB2=jsuWxR6FJwZ;O{;-jxZl|h(>g5*sG8ys?clD@N&eu>3Mx|R`~P<fca)re5G%*j
z!$8%M=}yu7{Quiw$gSlVDhx@1I#QlgHaw(U@){u8VnJFHH`#pUZ<2Oumh?!)`WSh&
z1ZoGQO1*HjtVdp$@9ZFR>W9Zq@r6e^Xfr#*W^jYo;0gXmKwe%}OT!t|$)xsOLn^TU
zy@eBj&p_icL6-N*YzX;cJ6vs>{PS-th?Tr72OKY)n29n@mAePSt!P5GFD?5UD{nHe
z@-E(7dBct%H=cFhZ#+k*R@!6+LayAk#!Bo|5`Ixe@SqlOP}|{mzNaDG03gkD9vkCB
zKl6(g6TB-t634EpdDJakyfF1eY^5IXU)2gE=3lN76=dfv-H){iA0WwQi7d=Hl)0SR
zMnRuf6DvPx8-w6pQ#-V9Dejq@Q`A;r?$pygZae=$-zALw$J#?KtBp7u%*FpLOnk%O
zfB9IK|2>pnQrr6!srQw<i0YuILLC$x&_UrglA{N3Whu(qabi2JWfJlHFpn8MygGAC
zSbF!%hXeMM%#Ck6LYYkLE~v{dv8l32BI8SwFC6#gxWALdKZVg`&n0^+Pg&wrIkmqK
ziA=)U(cae`0W1M*qEcKa563T%{Zch4M`mz8$s-H|#h`ncweHw7UVM~qne%CjRrmz7
z#?lnl|2KMBP$Kc)rd&Sb7z(JrEKfYDQ%eg#cJ{m4C?YYJJ`m<IXY~gq5okOWQGIJX
zMV&UJw3>jB__A(ywkVfC1h#tK?pNP|^><}Gn<CsFQ1jD``%~jTqoi?gNE+wQ^8%z1
zUd{M<3GfvHKaCo^plLfB*lxR=&>3@F&bk7eF%Rva+WKCY@Z)!s`rZ5)B!jb=m<oCL
z74%1m;-ovbL2!bUVd(fuh18f#*ZbwbGf&YOu=3>jD-?|Ai>4;t`S*D1u4si!{8^Ta
zJOWoe*lkwzs0zjR!vrcO$N?rSFd8@3yQ$w3X1U#6|L!rKY87+JbC;c3RH!(f8b@X#
zGLGkIt4FO`M+_dUu_|=hN@anIvb^G{6Y|uQl!c)=eGul7W=s8t@_3hTnwkTiX(r4?
z+H+~0a2_LupOULhU&ACPo3KrK?bV>7NVBuZY%1rwdBi|V3N-0B&bl1sV8L+nD-2Jg
z(xR9mKfSOy`2^BE&FaZ}1GXPG^3Su5otJ_2-)J2><7bi~tmTcFm$kx|=xQZ6VJoxi
zyvT9Ph5+t6M~2`=ay<7>YCuyR1IUqa`A0_v->vH+&>XA~nuF=Fazk@4!xx=7r?D0P
zguPEb1m<87V8uowo$H1;v0h<Gl6-bg1modQ*Df!-Cr3C-ZvvW2=b^gBNlF7t(w6d6
zO`G`~O(OBXsamggF>{|JR&l4oJr-i|igCcB+#uHR^%~B85cA!6f)MLAoIl)EP3@`M
z;+J?|T?`4j=d<6ev;E>%6Zwm}i(<%A>|OUlG?$X>DjAb>vO>E{+F8@K5kQ|`C%c@n
zgeB>91ILS2?mF{<Sz^3ix1srKkf1ebCELy0xPL4JiX&}#WDEM?r>BHB@8M$WZ4a3<
zSNRDy*E$SpI3_V4QlYCG395TkerD_{Q;RbH$4KNrIt&^}YcL$a0VyIQ8{90ludtJa
zA+m1_iz)6TCw_39U9!4eUmqhnV^f_IIngw&=Att|ab{9TRnahc@CD6DJ(c6J8AqGe
z^lAxPxnhLL@JG7OB)n1g4rVO{!mxLe;_u9$JX4Sj-11wVNx5b|siCu1D(Nh#_tNuu
zaQSTX=v_fe?)cXzk6MU~bs>Gt^}>;A-UZUOBs-$~*oWpepNQB@;$BsWUUHB!3)*cJ
zXBYZzK2l*{@a`!dX3l+bj>aj7iZhAgYSm(?lHa)1#mp7#>a~+?oHIq$Pv>~^De0ds
zDUDm>#>U-L6pys$aT+O68jrpV?i6Gdi8cJ41@X~jMBgTx!-+qC(mHZ+(2>jOgNWvY
zJpA?i&NuX1QsI&BrNXmSNgiI2{<+llCiD%g+*plh<B|y<KD_y0de|WOLInL<<iyw*
zc#j|HTNs*c9dTX8Bz>-5*7CV@jF{TP@zkDw3Kv*7ok&9yr)_Y6n%r&HQ{v@psg>V&
zJ{QXhc?H=fZP3cfj7R>NvTA4d$uEinK^Humb2V0EhHm(foRF#!2HG$;SYp(;Tz=LA
zivZ?>$Dr(<4Q^zj|44qA=Qkc=$l94VoHHdIl)ZJ40|m+a3TY_eGJksRmrrH)A};+(
zaQ8EeR0aIpryR?kQp=%z-TX2_)8FDsnebI|?C^Q%B_DsgezE3tVfB}I5iT9-t%hr)
z^Boo`M14iwS&il0+@BT?*QZ4v8qU=R-_YDK)dnH160-RVWxeB5hUpYS_%|r%HYHns
zaZLhu^Jp=o6_8_A_BP|5bG21rt`v{R{d)OI<?X)x__|e8LX&(Bx@SIb!k_!}5`p^|
zPw4j1bAWBeCZU2s-W6?wL9p$wn_aT;xvS;^O=J3_7d${!#OwTKA)!7jA2?<8w{gn;
zYRWnP?G|gNR{w82vfnJR;6gc2QE4?c;+y5vZc8&jyQ%yoq~_W2TG1ZmD-`ZY=lvd*
za*|M^7xY5N3>!=HebJdYL?~sT2nQywXsXZr0h%{z3-J5!=&7YVAALi$BgZi%E(Va(
zZ_LWgH1sX%(*{RtO=1be24*&tE-QC}$<`vI9(>-1IO0{T?BH2z42tgvYP@+1mwAr)
zWDJRSjAaUj`p@$`zZqSi2p!Vs0!8SkHyuz1H;gV@VD~)aq)rj(;5>aR)FAH?wWqm+
zTI)h_;Hp?q&FEaJU!g3ZgXkoR#rI~2=@qUycQ0Sj0Aax)11Ky=2R&DL@La>db7hW_
z0)g-@sDwJh7P%L4{OaqQ6;ly0GeT%Y>Ia)>M3DS~y_}g<6B<pA8#{F!a$W65&Lpd$
zdPvp*+nZ7v+OVOjK24Kspc{>YdS9#j@s+6FSF}b^w%0KJ0!}1EsrJUvV=q3Rt4xIx
zAIkrVcHmG`XMA4<|0R2g!{p)NF~vzxcP9+h-FaUcFpa*3?C#tORaS`_!b>H+z5qJ@
zB?GD-=p8JDQ|UTW?0pm;IS^}WO4!?59l;TugusX)zmpbBI;J(x8w@SHVs_6s|DAh%
zh={h!d;OidUjH!_k&#+V4Dj}Ugqla3XGn|kNx$m0MN1fyI~SEhIE98k!;<!V^(8AE
zM3i$1FKmm3HZ6au5o~N)x>)ph<nY8_AvY~LVAI0E`xRE@%8H_~2k(AE@QXQo&kt)1
zVCy<k<9?j^ii3|c^wTiQNOtvw=k2_loIR(9^k{^+AKU?B_~|!P5z$t%TWu(8z$5&N
z9+&Pzzqbj{XCbaRWS54i?F9$hzTi-5lmSKXBfjYf;$F<nnTZEl;UV{pj)x2TvS{WC
zM}bfq#rfqG`DoSOd%YETuR^?vf9Lgv@?S?dQ~j(_7Z;RAZ(HhWHq^QNzES`{DdX8r
zBH=4n7*mApnY0LVya|t^M0(hR;ZOH)J7kjuM)dG>{kH*#Sh<}_k@<kTX~z1CA7?P{
zQ~U}$o}Km;re>-c{^bz7mEXNv=4F2I0#*p;4K2v^`WpuZ04ev8olMndU<@Bv8Nq^^
zeQBwu7E+Jrk=k37))84y;MmeC<{qT!2FyS;qi7=W36J6VvJDZm^C{8V#yW>-N8Nz-
z$@~?@p28maqNSjOKl)~qQI;KV{}i&-b@iVzWfJ`ifxSn~0ck)VkOm?~AkqNt?2^T(
z+uRygL(j_wZuqGH00E412bQ3h8zvrFl!*sTAXcqMNpGRi_&&<RdlKp4JzXo`e3$qx
zdwS7T6rb;?vG-Ny@y}RdicALDSGMX?iKzNiuhKZ%clD_k9I5O+8-<<Dt-=n~=+tB*
z{&TE?BB!Hn3vv%w{odOaJM8|U2W|?h9-)L)rBY|~V%&O;&x(>uBZtQgq5p6m4u0Xq
zKRxZQTT+i$63i?IRTlm^P0hKK^t`h_LBDf$ay2_R{XWa#PwNZV4$sGm)<2YxElpYJ
z@{2O!4p%1_CPOlHT9#K&3usG~-Vg%xNX1dN1F!vReX5Ee)9q0K9CUPy|ELX>3Hu*V
zhWZfj`|NydY)XHM`ppXO=8qQcf`8!q6Zu3&tcCDEf31D;<M75PBv`{A0&74Ve+aCR
zI<wM*MuTdJS$F=hW8yoV>0Di!Uz8OE$Hr3-k8^K2<eZ)K-8&@<-YI&T>bAqmh_lMR
zBf9;PX|)AWH=onu?7w!y_Jp?jlF8idUmSSq{_4Q<b|vPouBYI?@AAL{MVG#xHvXLh
z&u<<pj{FDRc%YLl2>{)AGn7tN`PQaRwoz&Wf^11C2O!9mC<yYYE+yANUu9N270s}S
zQ-OLl!et~R=uV;cMv*>Qg)Wp&R)N8Q@m4soi;%oPscqIQt_zEMFL4)Rpo?s*{PKxH
z>3pbLaki;X=KzkRq>Q@khZSojITS=|;+Pm>w(U-w^a>^-mMr=xgQLH~%j*fq>AA(z
z-vSLjH^v8`L0148RDH;&6o~{4o;qP_UkZNTIUvwrFR!PbMq7IEN@rT#PfrZy;o(r$
zz-@|@UnLo!<W~%Y?1X8O*q5FqKE27WwCqxzFEr`?S;*VPKc;s^4DtLN>fqo2i0S|}
zm@sF<N3oS=h%U7Jn1ym*B}N~`cJ!A$D3LH%{1rA;x$c;|^Bi_prrZnjJ56G@!|?x(
z>2>!K+J8|g1!x<LKq<9tChjk22}tI^<^>ByODLQmmHg4(w$zNJemxPa8KYKOee~_L
z&sb_`(;R!=7ixkwKe{>!>RH=2T-+$m+z;Vr9=-w7cCa`zy5x?t4>@)(Xr7=F_p^Cc
zMFZ^x9J%qlo?J^rRiWMhSj&Y&YU(*EXvaAj_xSjG@pX#a5<bW)VPl`N+P39I{)-^r
zZO|agt<hcLF+p2J^@DD;>uzx0MK-waMx@25<Vv!Ym4EmS374=uW*YXO8erf^SuD?y
zlMiN%3f2nOepZ!+Bsk8|P@D>pM|tLAe^q1JMGEpaJ50Z!T)LhyLielJZ~d=hY(d|@
z!+yhb<n}h8o<UjocZ^&=mm^O$2(J4Bxi?r`K)zn?9eur%Oiwei0|tP>lDO@DLFhip
zdE~_7Zs1I(vg%ovR)96Mn$~wcX+rRaQ4(4;vg>Aj#+!A>63kk-uHr-W$g-I3E?xQQ
z@!+Fa`WLDGz$aD|;5-j<Sjhl&Gk#8hX8bD*cib4RgcFva9FLsu2=O|_L6<Jg@?fD!
z4+Fxw@&p}lO)e7r;hGGA4VC^3HarEc$>$~5|1|SZd*(*}WU+D6f92suUtXfyi?<Xa
z(h8)i2pM^OCe_<BVxlU-DF+!b<GCHi3Y#!;n|;*CAbOtwMgLOZtp<@jVtX}Ks*IN~
zGp76o&bp>p8z&%e2MCH15YAQnZ1joJ4|2uP8#o%s&y}K+)a*4un_#(4D0PqlMmd)W
zlfFWs@R;}r4i4@mF@b&fZsYGeN(Y^sS^`;M(xQ151YPWCz87lnq5ts;@hfzm>rH4r
zAB8Mq%cv;TRs~wJiY~5qV992^9m;Sde#GIa>EV&~rJA8@dM;ILLwH-ft?+iC`i}Ew
z<mrJ)+Y3bL8fn0L7?U8FGA7Su9oU)=jJKfzf#IithP3SjIp(aT+Ty@%fJAEmNVG&G
z63uXS*#?p5VOU&ENk{N0nHWx&zGlPyAYI-1+B<P0P`e?>whO8TwHtsY>kX$bXib_(
z;vtIfR0QprN&x;oE3k+3-Or8xXuuy>rBD2P26#yKOK_etI`t!g;Q}9##}b)1n2FqP
zOr2CGXqBg(%6hv*B`H;!9&`z2iBWHx)jIO@(~%d2clZaqTW7z*_%7Ie*q;NBn4Yg-
z>%nRr$_wU@`DBtZMtthkB?rGARDa~O6@Sg&Op+whR7T9u!>Q&+**GU(v=(VD{pjy<
ziu!MQ(D>gLlfr5wO=sqEmg{pwt4ZZ}`{^tXMm%5Q^XHP0uLbin8>~|oTykP#0fK-w
ztGlShFf07DA4fNtGgUO-cVeffnE4XX@I(a>OHFED5{$j(7)>jrd1~<Z_8tW8bsKpR
z;fPB?B`q6!pH;Aj%6~KOPBw9#^#Bvi7kcmF*F<`AoIqSF0Lmo+k=X?~`PUXQDw%c2
zo*R|V#b*fY15=6MR2v&Tjk}9PZReDgV|xPh_MBN09k8dJ@OpQsyNfU&T}7pi91N5J
z?gf+g+6^bwmK3w)v&rW4^zIxEXL)4kpADv5;UE~9DVnUP2Y`(*0N9X60ycU_XI7@a
z$*9jEWz-f2wThQQpy>s01m6XC2{64Nwu>k);gPpUR_s(joR!JW04I4!N~Z=v`G0P^
ziB-N%TauUWk`>H+f*ooi)hGQE-0-&bE+oN%?I=h>E9&f0Y`T_c=nNV$&hiu3e#tS<
zY0I2xyFTMhbl&j_;;7>N_@Q?BD{R5PdX?@g>_9cu49S@H3fwRlmh}~O2+QzTC%jLc
zu#SH!_C2g|NObMY*#|uo(`O&q(%*@|TUCDkDza+fE9`dKix56tVk}>r1En%`Xuhkk
zefQgDqMzdW*is^JBp>yhn(hfp9nTNtqkLm-IV=WO+_Rif^NCq(&$&e6(rU7NkI=cV
zFd}oj{t)kds|;_lmd&nFib-`iC$f}je_{SCRghaX;r)h$_J;9ei;~-t%UIT8Utup+
zpMxJ~tLn99EFt>z1M&}64U@;j5|(bA`U=aKKcz&)s(D={(Z70)<;4A`6d%w-r*&M`
zNtWycBNUd*2VFCP$xm-O<c|~AWM0E#q<uZ2`U1*`iTC<E7OcFpDmC2`JwWADysuSF
zxsC!$Hc!EP`4ifaGxVR#W!jIV-!$Vue{5ukD6gQ}$C8RIX7@~YourkOPKx2v&1^YK
z@;U4RljnI2+>52jjd+M+Lzp%Hl=~f9%(RTDg{LFDN(}_Pl;J6t&@zZmI952V3V4au
z!%VV!ae4M7w|91}l1$%KQh#*I+?S5@K}eYJg%ROij%5=BdXt9^g~i0tL0;vT>qJxQ
zo;vZDV()Ye4RqLCVXMp$YCqv4DLY(U->yT<LDRr@UeTJu?`6igEpC3gOGUI#l4c9-
z&H7{1_L}5pN$4I@FE;GEyL!8cwJz1U7spyMc`P};`vVq!a8YQ2;*c<2;4B8OIKm*8
zK~lVFpK|KH)#EN#n00&jyxJ<7W0vMM64s|)7dqcOK-|D}bFhRYZuQcj`&SsL?W@HV
z&$;@wcs)x3|BCT@K@?gZ&VzY!DqlwO#lr6QF{qK)++yS|H$CTW_1;~G2)n(XO_O@m
zV{S53j>yX;><j%k5jgs|zei~GbYC*czUQ9FY?=1(u?O}MjA1v}n3<%bm4(OZY0c!N
z#vAWac|_tmM2wg&cgi)TJI0fu*%MJ~c4moAVN6i6E~%z0=~U(%?Jd`c<BJ$Q7^J1^
zr+H=)!vY(JbIhTm=zba14DFKLdKvf(uQ0-2*=aMeq>_%VJ0JWCn<IFNW4ItZ72wAj
z9y@^HWj8^8=sFyJ?`W<d)syE>8d_{kUs3wc)mwLeVO;VPZH`c%@_KGD1fTsVk9#Pw
zs<eiXwUMKgrq1Bt#H9L_rk~4OlmB3Oi-_=VCD(RfVPHJKpqpi#SM!dojMhU{UcLt#
zTxGvX$0v9-;X$mdmAc*1t*<aal^xL%Utw80r-e+*636fR{}U4Le<UmfmFCE#a@{~_
zj;sch=7<(`9}Ec(8F}PYb6^8@IGIMMWzk?gD9zm<N`qzt?lG!=|0eqybdy~TZnAM_
zL=jGAY~u)q#-_fb;I#*JheQ?|K1RBYx`Pp@JAmjlKCol+$auD$4iF@rjy?+qX`yx6
z&%Y?8{fR2^<FXE7hQEh)Fh0Klki?1ogQaRVtO`4RR@muFr3?8zwU~k(yE~9@ddF7Z
zcTs>JCu5sAJn^;PN;Gc2+S=Bz_FvQ4*=iGQyDQwe&gWa@TlMM;hA{CN5C)i~IC?eh
z5p8bv155cMENL{>W%!||!Lks87ycIR^gUaUX1L@2(dH{P<nt8Vo6_zNRUVx+21)jA
z6Z>V~+5pv-p_C1`dEcyR^26Gmc4bwOMLJJyDq)iKtctSVNZ~=e5K%lWVbR0Pc#a2Z
zju%!7w|t^eOwpexb%Vh5nz_d*ksHnRgnt$ey1lIKeP?uX6k+oonumlaued}TYFe?&
zyfFQ+Hw8j%uyn<<!s2jlGpR6Ktd##=oWx}VkpuD^qiFUJOy+_&fXU<y!Rv5}JONC`
z>(;<q5pv)6ac+a;pO!BdNL)O(S4(TNFni#drfI*;@sDB-H`8*n+z_`Y2$0&h*L+Bl
ztVy;SO0`(IrCOZYf`MlK1Op{dh!-BaJPO3uQcuQH#U?>9N&Qz`--VskoL&f48%Yyn
zB6s)sf#3LMA>D;gTnJ-BsJ16gZA^}6R{52UroJg-#-u~cBb5ro@Ief{KC@~;`47lu
z0Hz$5c0l-TI2JwpcG&7GkMn0NJ=aD`$-G-*a&th5*%{HI5OqDpi=7_*?O`jYF_~ro
z@`aGTM(i<*uQ2^L)|}7Z_l)vc6P|xQ`rf`KIWYv0a?=;f)mK4lpXv!R|0nZ%DNC+~
zyTL{?OZF{kiWYp43+k_@Y#F6q&R&qND)qjIK~N|xxX%PPna>S2GC0Nl@OTMP(~{<8
zv8GE5A2kje^+aY5QCYRWHT!by>igs@tSd|tp^i<Br1?j11gX7@-SUKBBG`J^qWD*i
zu%OA_xC!IJtxeSH(Dw>X_vEVTXu7(t<#@oKJtw}DgkPyZ2&RRLqg&ERBfnntrX<U5
zbn#Q#^58m#XE(=3r8#MVeeC$-r_gd2Vj(MA#6w1Wn7y&r+st462XN*Qf>D!eW-)@|
kIN(b(3mm{8o)81YGJuK9l-dI!;G1{-uL1HH)z|m`2Lbv=CIA2c

literal 0
HcmV?d00001

diff --git a/analysis/The Coalescent Process.Rmd b/analysis/The Coalescent Process.Rmd
index 95b1f07..ab38a26 100644
--- a/analysis/The Coalescent Process.Rmd	
+++ b/analysis/The Coalescent Process.Rmd	
@@ -2,12 +2,12 @@
 title: "The Coalescent Process"
 author: "Jennifer Blanc"
 date: "3/13/2019"
-output: html_document
+output: 
+  html_document:
+    self_contained: no
 ---
 
 ```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE)
-rm(list = ls())
 set.seed(12)
 library(ggplot2)
 library(reshape2)
@@ -16,12 +16,16 @@ library(reshape2)
 
 ### Pre-requisites 
 
-* Wright-Fisher Model 
-* Intro probability  
+* [Wright-Fisher Model](https://stephens999.github.io/fiveMinuteStats/wright_fisher_model.html)  
+* Intro probability    
 
 ### The Coalescent Process
 
-In population genetics, we are often interested in modeling the relationship between a sample of genetic variants. The coalescent works backwards in time to model the genelogical relationship between the alleles in a sample. Coalescent theory was pioneered by John Kingman in the 1980s and remains one of the most widely used population genetics models. To understand how the coalescent works, it is helpful to first think about the wright-fisher model where we have idealized population with 2N haploid individuals, a neutral genetic locus, a constant population size of 2N, and offspring being randomly assigned to parents. While these might sound like idealized and unrealistic conditions, we can still derive interesting and applicable results from this model. Furthermore, it is the foundation on which a variety of more complicated models are built.     
+In population genetics, we are often interested in modeling the relationship between a sample of genetic variants. The coalescent works backwards in time to model the genelogical relationship between the alleles in a sample. [Coalescent theory](https://en.wikipedia.org/wiki/Coalescent_theory) was pioneered by [John Kingman](https://en.wikipedia.org/wiki/John_Kingman) in the 1980s and remains one of the most widely used population genetics models. To understand how the coalescent works, it is helpful to first think about the [wright-fisher model](https://stephens999.github.io/fiveMinuteStats/wright_fisher_model.html) where we have idealized population with 2N haploid individuals, a neutral genetic locus, a constant population size of 2N, and offspring being randomly assigned to parents (see Figure 1). While these might sound like idealized and unrealistic conditions, we can still derive interesting and applicable results from this model. Furthermore, it is the foundation on which a variety of more complicated models are built.   
+
+```{r, out.width = "200px", echo=FALSE, out.height="300px", fig.align='center', fig.cap="   Figure 1"}
+knitr::include_graphics("Discrete_Coal_Pic.jpg", dpi = 50)
+```
 
 ### Time to Coalescence  
 
@@ -83,9 +87,16 @@ We have now written the probability of a coalescent event in a sample of size n
 
 $$E[T_n] = \frac{2N}{{n}\choose{2}}$$  
 
+We can now return to the wright-fisher popultion from Figure 1. There we had a population of size 2N = 10 and we were considering a sample of size n = 3. Looking at the discrete generations in Figure 1, we can see the geneological relationship between these three samples is highlighted in dark green. This geneology is reproduced on the left in figure 2. The right side of figure 2 shows the continous generation coalescent model where the time to first coalescent event for our sample, $T_3$, is shown by the dotted line in the middle and the last coalescent event, $T_2$, by the top dotted line.  
+
+```{r, out.width = "300px", echo=FALSE, out.height="400px", fig.align='center', fig.cap="   Figure 2"}
+knitr::include_graphics("Cont_Coal_Pic.jpg", dpi = 50)
+```
+
+
 ## Simulating Coalescent Times  
 
-Finally, we will simulate the coalescent times for an entire sample. Before we do that, we are going to re-scale our expression slightly. Since we are typically working with large populations, the expected time for two alleles to coalesce ($2N$) is usually a very large number, which can make looking at coalescent times on the scale of generations difficult. Instead, we are going to set $k = 2N \mathrm{ generations} \Longleftrightarrow t = 1$. So each $t$ time-step backwards, actually represents 2N generations. We can now re-write the probability of a coalescent event before $t +1$ as:   
+We will now simulate the coalescent times for an entire sample. Before we do that, we are going to re-scale our expression slightly. Since we are typically working with large populations, the expected time for two alleles to coalesce ($2N$) is usually a very large number, which can make looking at coalescent times on the scale of generations difficult. Instead, we are going to set $k = 2N \mathrm{ generations} \Longleftrightarrow t = 1$. So each $t$ time-step backwards, actually represents 2N generations. We can now re-write the probability of a coalescent event before $t +1$ as:   
 
 $$P_{C,t+1} = {{n}\choose{2}}e^{-{{n}\choose{2}}}t$$
 
@@ -186,8 +197,6 @@ dat <- melt(data, id = "n")
 ggplot(dat, aes(x = n, y = value, color = variable)) + geom_point() + geom_line() + theme_classic() + scale_color_manual(name = " ", labels = c("E[Ln]", "E[TMRCA]"), values  = c("red4", "navy")) + ylab("Time") + theme(legend.position="bottom")
 ```
 
-```{r}
-sessionInfo()
-```
+Figures from: The Coalescent and Measurably Evolving Populations Alexei Drummond Department of Computer Science University of Auckland, NZ  
 
 

From fc0216e268fa4d2d53bcab971d2d3b977796e4a3 Mon Sep 17 00:00:00 2001
From: Jennifer Blanc <jgblanc@uchicago.edu>
Date: Fri, 5 Apr 2019 11:31:04 -0500
Subject: [PATCH 3/3] Delete The Ancestral Process.Rmd

---
 analysis/The Ancestral Process.Rmd | 195 -----------------------------
 1 file changed, 195 deletions(-)
 delete mode 100644 analysis/The Ancestral Process.Rmd

diff --git a/analysis/The Ancestral Process.Rmd b/analysis/The Ancestral Process.Rmd
deleted file mode 100644
index e71d295..0000000
--- a/analysis/The Ancestral Process.Rmd	
+++ /dev/null
@@ -1,195 +0,0 @@
----
-title: "The Ancestral Process"
-author: "Jennifer Blanc"
-date: "3/14/2019"
-output: html_document
----
-
-```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE)
-rm(list = ls())
-library(reshape2)
-library(ggplot2)
-library(expm)
-library(ggpubr)
-set.seed(12)
-```
-
-### Pre-requisites 
-
-* The Coalescent Process
-* Continuous Time Markov Chains 
-
-### Ancestral Process  
-
-In coalescent theory, the ancestral process is the number of active ancestral lineages at time t. Recall that we can think of the coalescent process as modeling the relationship back in time between a sample of haploid individuals at a single neutral locus, in a population of constant size 2N. In the previous vignette, we showed the time until two individuals coalesce, $T_2$, is an exponential random variable with parameter 1. We then extended this idea to show that for a sample of size n, each coalescent coalescent time, $T_n, T_{n-1}, ... T_2$, is an independent random variable where $T_i \sim exp{{i}\choose{2}}$. So each time a coalescent event happens, we are transitioning from $i$ to $i-1$ lineages after waiting $T_i \sim exp{{i}\choose{2}}$ amount of time.  
-
-Keeping this in mind, we can define the ancestral process as a continuous time Markov chain where the states are the number of active lineages (1,2, .., n) in generation t:  
-
-$$A_n(t) := \text{number of distinct lineages in generation t of a sample of size n at time 0} $$
-
-### Deriving Kingman's Coalescent
-
-We will first re-derive the results from the previous vignette while explicitly considering the ancestral process as a continuous time Markov chain. Mainly, we will show that the coalescent process describes the ancestral process for a sample of fixed size n in the limit as N approaches infinity in the Wright-Fisher model. Remember that in the wright-fisher model we assume a fixed population size (in this case N) and that the j parents of i individuals are sampled with replacement from the N possible parents in the previous generation. This process is described by the expression below where i is the number of active lineages in the current generation and j is the number of ancestral lineages in the previous generation. 
-
-$$G_{i,j} := P(\text{i individuals have j separate parents)}= \frac{N(N-1)...(N-j+1)\delta_k^{(j)}}{N^i} \\ 1 \le j \le i$$
-
-This probability can be broken into three parts:  
-$\bullet \ N(N-1)...(N-j+1)$ is a descending factorial that describes the number of ways that you can choose j parents  
-$\bullet \ \delta_k^{(j)}$ is a Stirling number of the second kind, it describes the number of ways to assign k individuals to their j parents  
-$\bullet \ N^i$ is the total number of ways of assigning the number of i individuals to their parents (there are N possible parents) 
-
-Taken together, $G_{i,j}$ describes the probability of i individuals having j separate ancestors (remember that $j \le i$ because individuals cannot have more than one parent). Since our CTMC is defined as the number of lineages at a given time, $G_{ij}$ describes the transition probabilities for our CTMC:  
-
-$$P(A_n(t+1) = i|A_n(t) = j) = G_{ij}$$
-
-When we look closer at this transition probability, we can see that Kingman's coalescent does not apply when we have a fixed population size N. Right now, i individuals can have anywhere from i to 1 distinct parents. However, we know that Kingman's coalescent only allows only one coalescent event per generation, in other words, i individuals can have only i or i - 1 parents. We will now show that for large N, we recover Kingman's coalescent.    
-
-First let's consider the case where i individuals have i - 1 parents.  
-
-$$G_{i,i-1} = \delta_i^{(i-1)}\frac{N(N-1)...N(N-k+1)}{N^i} =  {{i}\choose{2}}\frac{1}{N}  + O(N^{-2})$$
-
-
-We can arrive at the above statement by first recalling that $\delta_k^{(k-1)} = {{k}\choose{2}}$ and recognizing that the other terms are proportional to $1/N^2$ and will decay quickly when N is is very large. 
-
-Next, we will consider the other two possibilities. First that i individuals have fewer than i - 1 parents and the case where they have exactly i parents:  
-
-$$G_{i,j} = \delta_i^{(j)}\frac{N(N-1)...N(N-j+1)}{N^i}= O(N^{-2}) \ \text{for} \ j <i-1$$
-
-$$G_{i,i} =  \delta_i^{(i)}N^{-i} N(N-1)..(N-i+1) = 1-{{i}\choose{2}}\frac{1}{N} + O(N^{-2})$$  
-
-Now we can show formally, that in the limit as N tends to infinity, the ancestral process converges to the continuous-time coalescent process. In order to do this we are again going to re-scale time so that one unit of t corresponds to N generations.   
-
-$$P(T_i^{(N)} > t) = G_{i,i}^{[Nt]} \rightarrow e^{-{{i}\choose{2}}t} \ \text{as} \ N \rightarrow \infty$$
-
-Here the notation $[Nt]$ means only the integer part of Nt. So while t could be any value greater than zero, the probability of $G_{i,i}^{[Nt]}$ only makes sense for whole generations. We have shown that the probability of a coalescent event, $T_i$ happens after time t happens with probability $e^{-{{i}\choose{2}}}$. This is that same coalescent time that we arrived at in the last vignette.   
-
-### Properties of the Ancestral Process  
-
-In the last section, we showed that as N tends to infinity, the ancestral process converges to the coalescent process. We can use the ancestral process transition probabilities, which we defined above to write the general transition matrix for this process:  
-
-$$G_{i,j}^{Nt} = (I + N^{-1}Q + O(N^{-2}))^{Nt} \rightarrow e^{Qt}, \ \text{as} \ N \rightarrow \infty$$
-
-Thus, the ancestral process is approximated by the Markov chain $A_n(t)$ whose behave is determined by the rate matrix Q and $A_n(0) = n$. Here Q is is an n x n matrix with non-zero entries given by:  
-
-$$q_{i,i} = -{{i}\choose{2}}, q_{i, i-1} = {{i}\choose{2}}, i = n, n-1, ..., 2 $$
-
-$$Q =
-  \left[ {\begin{array}{ccccc}
-   0 &  &  & \\
-   {{i}\choose{2}} & -{{i}\choose{2}} &  & \\
-   0 & {{i}\choose{2}} & -{{i}\choose{2}} &  \\
-   0 & 0 & {{i}\choose{2}} & -{{i}\choose{2}} \\
-   \vdots & & & & \ddots \\
-  \end{array} } \right]$$
-
-Intuitively, we can think about $q_{i,i-1}$ as the rate individuals coalesce, going from i to i-1 lineages. Therefore, $q_{i,i}$ is the rate that process is leaving the state corresponding to i=1. Since these are the only two options, we know the $A_n(t)$ is pure death process that starts at $A_n(0)$ and decreases by only one jump at a time. We can use Q to calculate the distribution of $A_n(t)$.  
-
-One way we can do this is write $Q = RDL$, where D is a diagonal matrix of eigenvalues $\lambda_i = - {{i}\choose{2}}$ and L and R are matrices of right and left eigenvalues of Q, normalized so that RL = LR = I. Using this method and some calculation, its possible to solve for the the distribution of $A_n(t)$:    
-
-$$G_{n,j} = P(A_n(t) = j) = \sum\limits_{i=j}^n e^{-i(i-1)t/2}\frac{(2i-1)(-1)^{i-j}j_{(i-1)}n_{[i]}}{j!(i-j)!n_{(i)}}  \\ s_{(n)} = s(s+1)...(s + n-1) \\ s_{[n]} = s(s-1) ... (s -n +1) \\ s_{(0)} = s_{[0]} = 1$$  
-
-The expectation of this distribution: 
-
-$$E[A_n(t)] = \sum\limits_{k=1}^n e^{-i(i-1)t/2}\frac{(2i-1)n_{[i]}}{n_{(i)}}$$
-
-Now we have a both a distribution and its expectation for the number of ancestral lineages for a sample of size n and different times, t. Lets plots the mean number of ancestors over time for different values of n:  
-
-
-```{r}
-up <- function(a,n) { # Function for s(n) calculations
-  x <- seq(a, (a+n-1))
-  return(prod(x))
-}
-
-down <- function(a,n) { #Function for s[n] calculations
-  x <- seq(a, (a-n+1))
-  return(prod(x))
-}
-
-
-mean_ancestral_lineages <- function(t, n) { 
-  tmp <- rep(NA, n)
-  for (i in 1:n) { # Evaluate espression for every value of i
-    tmp[i] <- exp(-i * (i-1) * (t /2)) * ((((2*i)-1) * down(n,i)) / up(n,i))
-  }
-  return(sum(tmp)) # Return sum over values of i 
-}
-```
-
-```{r}
-# Calculate the mean number of ancestral lineages for n = 5, 10, 50 
-x <- seq(0, 2, 0.01)
-df <- as.data.frame(cbind(x,(sapply(x,FUN = mean_ancestral_lineages, n = 5)), (sapply(x,FUN = mean_ancestral_lineages, n = 10)), (sapply(x,FUN = mean_ancestral_lineages, n = 50))))
-dat <- melt(df, id = "x")
-
-# Plot the data 
-ggplot(data = dat, aes(x = x, y = value, color = variable)) + geom_line() + theme_bw() + xlab("time") + ylab("Expecent Number Ancestral Lineages") + scale_color_manual(name = "n", labels = c("10", "50", "100"), values = c("red4", "navy", "darkgreen")) + geom_hline(yintercept = 1, linetype = 2)
-```
-
-
-Again, we can see that most coalescent events happen quickly as the number of active lineages quickly decreases and that most of the time is take up by the last coalescent event.  
-
-
-### Example  
-
-Let's do an example for a sample size of n = 5. Recalling that the ancestral process is a pure death CTMC with death rate of ${{i}\choose{2}}$, we can write down the Q matrix for our example:  
-
-$$Q =
-  \left[ {\begin{array}{ccccc}
-   0 &   \\
-   1 & -1 &  \\
-   0 & 2 & -2 \\
-   0 & 0 & 6 & -6 \\
-   0 & 0 & 0 & 10 & -10
-  \end{array} } \right]$$  
-
-
-Now, we can use the fact the the transition matrix, $G_{i,j}(t)$, is equal to $e^{Qt}$ when $N \rightarrow \infty$ to calculate the transition matrix for different time points.  
-
-```{r}
-# Enter Q matrix 
-Q <- matrix(c(0,0,0,0,0,1,-1,0,0,0,0,2,-2,0,0,0,0,6,-6,0,0,0,0,10,-10), nrow = 5, byrow = T)
-
-# Calculate transition matrix at different time points 
-P0.1 <- expm(Q * 0.1)
-P0.5 <- expm(Q * 0.5)
-P1 <- expm(Q * 1)
-P10 <- expm(Q * 10)
-```
-
-
-Let's visualize the transition matrix with heat maps for each time point:  
-```{r, fig.width=12}
-# Prepare data for plotting 
-dat_plot1 <- cbind(melt(P0.1, varnames = names(dimnames(P0.1)), na.rm = FALSE, as.is = FALSE, value.name = "value"), 0.1)
-dat_plot2 <- cbind(melt(P0.5, varnames = names(dimnames(P0.5)), na.rm = FALSE, as.is = FALSE, value.name = "value"), 0.5)
-dat_plot3 <- cbind(melt(P1, varnames = names(dimnames(P1)), na.rm = FALSE, as.is = FALSE, value.name = "value"), 1)
-dat_plot4 <- cbind(melt(P10, varnames = names(dimnames(P10)), na.rm = FALSE, as.is = FALSE, value.name = "value"), 10)
-
-# Make a heat map for each tranisition matrix
-p1 <- ggplot(data = dat_plot1, aes(x=Var2, y=Var1, fill=value)) + 
-  geom_tile() + theme_classic() + scale_fill_gradient2(low = "white", high = "navy") + geom_text(aes(Var2, Var1, label = round(value,2)), color = "black", size = 4) + scale_x_reverse() + scale_y_reverse() + xlab("j") + ylab("i")
-p2 <- ggplot(data = dat_plot2, aes(x=Var2, y=Var1, fill=value)) + 
-  geom_tile() + theme_classic() + scale_fill_gradient2(low = "white", high = "navy") + geom_text(aes(Var2, Var1, label = round(value,2)), color = "black", size = 4) + scale_x_reverse() + scale_y_reverse() + xlab("j") + ylab("i")
-p3 <- ggplot(data = dat_plot3, aes(x=Var2, y=Var1, fill=value)) + 
-  geom_tile() + theme_classic() + scale_fill_gradient2(low = "white", high = "navy") + geom_text(aes(Var2, Var1, label = round(value,2)), color = "black", size = 4) + scale_x_reverse() + scale_y_reverse() + xlab("j") + ylab("i")
-p4 <- ggplot(data = dat_plot4, aes(x=Var2, y=Var1, fill=value)) + 
-  geom_tile() + theme_classic() + scale_fill_gradient2(low = "white", high = "navy") + geom_text(aes(Var2, Var1, label = round(value,2)), color = "black", size = 4) + scale_x_reverse() + scale_y_reverse() + xlab("j") + ylab("i")
-
-# View plots 
-ggarrange(p1,p2,p3,p4, common.legend = T, ncol = 4, labels = c("0.1", "0.5", "1", "10"), hjust = -.3)
-```
-
-In each transition matrix, each cell in the lower diagonal matrix represents the probability of transitioning from i lineages to j lineages for the given time period. If we focus on the fifth row in each transition matrix, we can see that for t = 0.1 the highest probability is that one coalescent event will happen and the process will go from 5 lineages to 4. If we look at t = 0.5 and t = 1, the highest probability is to go from i = 5 to j = 3 and j = 2, respectively. Unsurprisingly, given that the coalescent is a pure death process, by t = 10 all of 5 of our samples will have coalesced and there will be only 1 lineage left. 
-
-```{r}
-sessionInfo()
-```
-
-
-
-
-
-
-