The Moran game simulator

The following compiles and loads the Moran Genealogy Process (MGP) simulator, encoded in the file mgp.cc. The contents of this file are displayed below.

stopifnot(system("R CMD SHLIB mgp.cc")==0)
dyn.load("mgp.so")

The following functions run simulators of the various genealogy processes. The first simulates a realization of the MGP (without sampling). It returns snapshots of the MGP at the specified times.

playMGP <- function (n, rate = n, times, t0 = 0, tree = TRUE) {
  if (is(rate,"data.frame")) rate <- NULL
  if (is(n,"data.frame")) n <- attr(n,"state")
  x <- .Call("playMGP",n,rate,times,t0,tree,PACKAGE="mgp")
  state <- x$state
  x$state <- NULL
  x <- as_tibble(x)
  attr(x,"state") <- state
  x
}

The second simulates one realization of the Sampled MGP. If drop = TRUE, the black balls will be dropped.

playSMGP <- function (n, rate = n, times, t0 = 0, drop = TRUE, tree = TRUE) {
  if (is(rate,"data.frame")) rate <- NULL
  if (is(n,"data.frame")) n <- attr(n,"state")
  x <- .Call("playSMGP",n,rate,times,t0,drop,tree,PACKAGE="mgp")
  state <- x$state
  x$state <- NULL
  x <- tail(as_tibble(x),-1)
  attr(x,"state") <- state
  x
}

The third simulates a realization of the Observed Chain of the stationary Sampled Moran Genealogy Process.

playObsSMGP <- function (n, rate = n, times, t0 = 0, tree = TRUE) {
  if (is(rate,"data.frame")) rate <- NULL
  if (is(n,"data.frame")) n <- attr(n,"state")
  x <- .Call("playObsSMGP",n,rate,times,t0,tree,PACKAGE="mgp")
  state <- x$state
  x$state <- NULL
  x <- tail(as_tibble(x),-1)
  attr(x,"state") <- state
  x
}

The following function extracts information about an MGP state. In addition to the sample times and event counts, it returns the final tree, plus the times of the epochs (the interval boundaries) and the number of lineages on each interval. Finally, it computes the conditional log likelihoods and cumulative hazards, of each successive sample, conditional on the previous ones.

getMGPinfo <- function (x, drop = TRUE) {
  .Call("get_MGPS_info",attr(x,"state"),drop,PACKAGE="mgp")
}

Plotting utilities.

The treeplot function uses ggtree [1] to generate a plot of the genealogy from its Newick representation.

library(ggtree)
library(tidyverse)
library(doParallel)
library(doRNG)

theme_set(theme_bw())

treeplot <- function (dat, times = dat$times, ladderize = FALSE, points = FALSE) {
  
  read.tree(text=dat$tree[1]) %>%
    fortify() %>%
    summarize(tmin=times[1]+min(x)-max(x)) %>%
    pull(tmin) -> tmin
  
  foreach (k=seq_along(dat$tree)) %dopar% {
    
    read.tree(text=dat$tree[k]) %>%
      fortify(ladderize=ladderize) %>%
      mutate(
        tfin=times[k],
        x=x-max(x)+tfin,
        nodecol=factor(
          case_when(label=="o"~"black",label=="g"~"green",label=="i"~"invisible",TRUE~"blue")
        ),
        tipcol=factor(
          case_when(label=="g"~"green",TRUE~"red")
        )
      ) %>%
      ggplot(aes(x=x,y=y))+
      geom_tree(layout="rectangular")+
      expand_limits(x=c(tmin,dat$times))+
      scale_x_continuous()+
      theme_tree2() -> pl
    
    if (points) {
      pl+
        geom_nodepoint(aes(color=nodecol))+
        geom_tippoint(aes(color=tipcol))+
        scale_color_manual(
          values=c(green="darkgreen",blue="blue",black="black",red="red",
                   invisible=alpha("black",0))
        )+
        guides(color=FALSE) -> pl
    }
    
    pl
    
  }
}

One simulation of the MGP

burnin <- 0
tout <- seq(0,400,by=1)
playMGP(n=100,times=tout,t0=-burnin) -> x
registerDoParallel()
x %>% treeplot() -> pls

png_path <- file.path(tempdir(),"frame%05d.png")
png(png_path,height=3.7,width=6,res=100,units="in")
for (i in seq_along(pls)) print(pls[[i]])
dev.off()

library(gifski)
png_files <- sprintf(png_path,seq_along(pls))
gif_file <- "mgp1.gif"
gifski(png_files,gif_file,delay=0.02,loop=TRUE)
unlink(png_files)