_results/figS5.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


# Takes ~ 1.6 min
t0 <- Sys.time()
ratio_sims <- crossing(eta = 0.5,
                      d1 = -1 * c(1.1, 0.6, 0.1),
                      d2 = c(1.1, 0.6, 0.1),
                      final_t = 40e3L,
                      save_every = 1L,
                      f = (formals(quant_gen)$f *
                               seq(0.1, 2, length.out = 11)) %>%
                          round(3),
                      a0 = formals(quant_gen)$a0 * c(0.1, 1, 10),
                      V0 = list(cbind(c(0, 0.01), c(0.01, 0)))) %>%
    mutate(save_every = ifelse(d2 == 1.1 & f == 0.01, 50e3L, save_every)) %>%
    split(1:nrow(.)) %>%
    mclapply(run_determ_1rep) %>%
    do.call(what = bind_rows)
t1 <- Sys.time()
t1 - t0


# ratio_sims[c(166, 167, 168, 265, 266, 267),] %>%
#     distinct(f, d1, d2)


for (i in 1:nrow(ratio_sims)) {
    if (!equil_check(ratio_sims[i,])) {
        cat(sprintf("Not at equilibrium on row %i\n", i))
    } else if ((max(Re(ratio_sims$eig[[i]])) - 1) > 1e-6) {
        cat(sprintf("Unstable equilibrium on row %i\n", i))
    }
}; rm(i)




# ------------------------------------------------*
# ------------------------------------------------*
# Abundances ----
# ------------------------------------------------*
# ------------------------------------------------*


nratio_sims <- ratio_sims %>%
    mutate(cp_N = map2_dbl(time, N,
                           function(.x, .y) {
                               z <- .y[.x == max(.x)]
                               return(z[2] / z[1])
                           })) %>%
    select(eta:a0, cp_N) %>%
    mutate(d1 = abs(d1)) %>%
    rename(dx = d1, dp = d2)

nratio_sims %>%
    {.[["cp_N"]] %>% min() %>% print(); .} %>%
    filter(cp_N == min(cp_N))

nratio_p <- nratio_sims %>%
    mutate(dx = factor(dx, levels = sort(unique(dx)),
                       labels = sprintf("d[x] == %.1f", sort(unique(dx)))),
           dp = factor(dp, levels = sort(unique(dp)),
                       labels = sprintf("d[p] == %.1f", sort(unique(dp))))) %>%
    ggplot(aes(f, log2(cp_N), color = factor(a0))) +
    geom_hline(yintercept = 0, color = "gray70", linetype = 1) +
    geom_point() +
    ylab(expression(log[2]( hat(N)[c] / hat(N)[p] ))) +
    facet_grid(dx ~ dp, labeller = label_parsed) +
    scale_color_viridis_d(expression(alpha[0]),
                          option = "A", begin = 0.1, end = 0.9) +
    guides(color = guide_legend(override.aes = list(size = 4))) +
    theme(legend.title = element_text(hjust = 0.5),
          strip.text.y = element_text(angle = 0))


save_plot(nratio_p, 6, 6, .prefix = "S5-")


# ------------------------------------------------*
# ------------------------------------------------*
# Investments ----
# ------------------------------------------------*
# ------------------------------------------------*

vratio_sims <- ratio_sims %>%
    mutate(V1 = map2_dbl(time, V1, ~ sum(.y[.x == max(.x)])),
           V2 = map2_dbl(time, V2, ~ sum(.y[.x == max(.x)])),
           cp_V = V1 / V2) %>%
    select(eta:a0, V1, V2, cp_V) %>%
    mutate(d1 = abs(d1)) %>%
    rename(dx = d1, dp = d2)

vratio_sims %>%
    {.[["cp_V"]] %>% max() %>% print(); .} %>%
    filter(cp_V == max(cp_V)) %>%
    {.[["V1"]] %>% range() %>% print(); .} %>%
    {.[["V2"]] %>% range() %>% print(); .} %>%
    identity()



# vratio_p <-
vratio_sims %>%
    filter(f < 0.2 | f > 0.25) %>%
    mutate(dx = factor(dx, levels = sort(unique(dx)),
                       labels = sprintf("d[x] == %.1f", sort(unique(dx)))),
           dp = factor(dp, levels = sort(unique(dp)),
                       labels = sprintf("d[p] == %.1f", sort(unique(dp))))) %>%
    ggplot(aes(f, log2(cp_V), color = factor(a0))) +
    geom_hline(yintercept = 0, color = "gray70", linetype = 1) +
    geom_point() +
    ylab(expression(log[2]( hat(x)[c] / hat(p)[p] ))) +
    facet_grid(dx ~ dp, labeller = label_parsed) +
    scale_color_viridis_d(expression(alpha[0]),
                          option = "A", begin = 0.1, end = 0.9) +
    guides(color = guide_legend(override.aes = list(size = 4))) +
    theme(legend.title = element_text(hjust = 0.5),
          strip.text.y = element_text(angle = 0))

save_plot(vratio_p, 6, 6, .prefix = "S6-")


vratio_sims %>%
    filter(f <= formals(quant_gen)$f) %>%
    .[["cp_V"]] %>% max()