---
title: "R Notebook"
output: html_notebook
---

## Illustrating the coalescent

```{r}
options(tidyverse.quiet = TRUE)
library(tidyverse)
library(phyclust)
library(ape)

n_sample <- 100
ms_out <- ms(nsam = 2*n_sample, opts = "-T")
ms_tree <- read.tree(text = ms_out)
plot(ms_tree, show.tip.label = FALSE)
```

## Comparing the efficiency of coalescent simulation with direct simulation

Here's a simple comparison of simulating drift directly and simulating it with the coalescent. I've set $N_e$ to 100 and $n_sample$ to 25.

```{r}
## These are the parameters you can play with
##
n_e <- 1000
n_sample <- 25

## There's no easy way to tell when the forward simulation has run long enough 
## to reach stationarity. I just fix the number of iterations to the expected
## fixation time and hope that that's close enough
##
n_gen <- 4*n_e

## The number of times to repeat each simulation
##
n_reps <- 10000

## This is the forward simulation.
##
start_time <- Sys.time()
for (rep in 1:n_reps) {
  p <- numeric(n_gen)
  p[1] <- 0.5
  for (i in 2:n_gen) {
    k <- rbinom(1, 2*n_e, p[i-1])
    p[i] <- k/(2*n_e)
  }
  sample <- rbinom(1, 2*n_sample, p[n_gen])
}
end_time <- Sys.time()
cat("Elapsed time for forward simulation:       ", round(end_time - start_time, 2))

## Here's the coalescent simulation.
##
start_time <- Sys.time()
for (rep in 1:n_reps) {
  ms_out <- ms(nsam = 2*n_sample, opts = "-T")
}
end_time <- Sys.time()
cat("\nElapsed time for coalescent simulation: ", round(end_time - start_time, 2))
```

## The coalescent with migration

As you've seen in lab this week, `ms()` also allows you to simulate the coalescent with migration. What you didn't see was the actual call to `ms()` that does it. Here's an example with Wright's island model: 5 islands with a sample size of 25 from each island and $N_em = 0.1$.

```{r}
ms_out <- ms(nsam = 125, # The total number of alleles in the sample
             nreps = 1,  # The number of coalescent samples to produce
             opts = "-T -I 5 25 25 25 25 25 0.1",
                         # -T produce a tree with coalescent times
                         # -I migration follows Wright's island model
                         # 5 populations
                         # 25 the number of samples from each population
                         # 0.1 Ne*m
            )
ms_tree <- read.tree(text = ms_out)
plot(ms_tree, show.tip.label = FALSE)
```