_results/figS3.R

In lucasnell/evo_alt_states: Simulations for alternative stable states and coexistence in evolution

# load packages and get helper functions
source("_results/_preamble.R")

# whether to re-do simulations (use rds files otherwise)
.REDO_SIMS <- FALSE
# whether to re-save plots
.RESAVE_PLOTS <- TRUE



dist_ss <- function(.eta, .Omega) {
    # .eta = 0; .Omega = 10e3
    # rm(.eta, .Omega, .f, .a0, .r0, .av, .rad, .ring_dists, out)

    .f <- formals(quant_gen)$f
    .a0 <- formals(quant_gen)$a0
    .r0 <- formals(quant_gen)$r0
    .av <- with(list(n = 1), eval(formals(quant_gen)$add_var))
    if (.eta >= 0) {
        .rad <- sqrt(log(.a0 / .f * .Omega))
    } else {
        .rad <- sqrt(log(.a0 / (.f * (1 + .eta)) * .Omega))
    }

    # movement along and to the ring:
    .ring_dists <- crossing(.angle = seq(0, 45, 15),
                            .dist = seq(0, 0.5, 0.1)) %>%
        pmap_dfr(function(.angle, .dist) {
            # .angle = 45; .dist = 0.4
            # rm(.angle, .dist, .rad_i, .m, .ss, ..out)
            .rad_i <- .rad + .dist
            .m <- rbind(ifelse(.angle == 90, 0, cos(.angle * pi / 180)) * .rad_i,
                        ifelse(.angle == 0, 0, sin(.angle * pi / 180)) * .rad_i)
            .ss <- sauron:::sel_str_cpp(V = cbind(.m, rbind(1, 1)),
                                        N = c(1, .Omega),
                                        f = .f,
                                        a0 = .a0,
                                        C = rbind(c(1, .eta), c(.eta, 1)),
                                        r0 = .r0,
                                        D = diag(c(0, 0)))[,1]
            .ss[(.m[,1] + .av * .ss < 0)] <- 0
            .ss <- sqrt(sum((.ss)^2))
            ..out <- tibble(angle = .angle, dist = .dist, ss = .ss)
            return(..out)
        })
    out <- .ring_dists %>%
        mutate(eta = .eta, Omega = .Omega) %>%
        select(eta, Omega, dist, angle, ss)
    return(out)
}

dist_ss_df <- crossing(.eta = seq(-0.6, 0.6, 0.2), .Omega = 10e3) %>%
    pmap_dfr(dist_ss)

dist_ss_p <- dist_ss_df %>%
    mutate(angle = factor(angle),
           eta = factor(eta, levels = sort(unique(eta)),
                        labels = sprintf("eta == %.1f", sort(unique(eta))))) %>%
    ggplot(aes(dist, ss)) +
    geom_hline(yintercept = 0, color = "gray70") +
    geom_line(aes(color = angle)) +
    geom_point(aes(color = angle)) +
    scale_color_viridis_d(begin = 0.3, guide = "none") +
    scale_y_continuous("Strength of selection") +
    scale_x_continuous("Distance from ring", breaks = c(0, 0.2, 0.4)) +
    coord_cartesian(ylim = c(0, NA)) +
    facet_wrap( ~ eta, nrow = 1, labeller = label_parsed) +
    theme(legend.position = "top",
          axis.text.x = element_text(size = 8)) +
    NULL


ring_df <- tibble(angle = seq(0, 90, length.out = 91)) %>%
    mutate(x = ifelse(angle == 90, 0, cos(angle * pi / 180)),
           y = ifelse(angle == 0, 0, sin(angle * pi / 180)),
           x = x * with(formals(quant_gen), sqrt(log(a0 / f * 10e3))),
           y = y * with(formals(quant_gen), sqrt(log(a0 / f * 10e3))))

ring_p <- ring_df %>%
    ggplot(aes(x, y)) +
    geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
    geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
    geom_path(linetype = "23", size = 1, color = "gray70") +
    geom_point(data = ring_df %>% filter(angle %in% seq(0, 45, 15)),
               aes(color = factor(angle)), size = 3) +
    scale_color_viridis_d(begin = 0.3, guide = "none") +
    scale_x_continuous("Conflict investment", breaks = 0:1) +
    scale_y_continuous("Partitioning investment", breaks = 0:1) +
    coord_equal(xlim = c(0, 1.75), ylim = c(0, 1.75))


dist_ss_ring_p <- ring_p +
    dist_ss_p +
    plot_layout(ncol = 1, heights = c(0.6, 1))


save_plot(dist_ss_ring_p, 6.5, 5, .prefix = "S3-")