_results/fig2.R

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

# start ----


# 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




# =============================================================================*
# =============================================================================*

# *--------------------------* ----
# Fig 2: Tradeoff effects on two-axis outcomes ----

# =============================================================================*
# =============================================================================*


# Takes ~20 sec
set.seed(34506729, "L'Ecuyer")
tradeoff_sims <- crossing(eta = c(-1, 0, 1) * 0.6,
                          d1 = -1e-6, d2 = 1e-6,
                          final_t = 200e3L,
                          V0 = map(rep(1:24, 10),
                                   function(.i) {
                                       .n <- 2
                                       m <- matrix(runif(.n * 2, 0, 0.5),
                                                   nrow = 2)
                                       return(m)
                                   })) %>%
    split(1:nrow(.)) %>%
    mclapply(run_determ_1rep) %>%
    do.call(what = bind_rows)



tradeoff_sim_df <- tradeoff_sims %>%
    select(eta, spp:V2) %>%
    group_by(eta) %>%
    mutate(rep = 1:n()) %>%
    ungroup() %>%
    unnest(spp:V2) %>%
    mutate(eta = factor(eta, levels = sort(unique(eta)),
                        labels = c("sub-additive", "additive",
                                   "super-additive")),
           rep = factor(rep))


# ## 1/5 have complex eigenvalues.
# ## This only happens when you have eta >= 0
# tradeoff_sims %>%
#     mutate(first_cmplx = map_int(eig,
#                                  function(.z) {
#                                      if (any(which(Im(.z) != 0))) {
#                                          min(which(Im(.z) != 0))
#                                          } else 0L}),
#            max_poss = map_int(eig, length)) %>%
#     select(eta, first_cmplx, max_poss) %>%
#     filter(first_cmplx > 0)
# ## All are either stable or neutrally stable:
# tradeoff_sims %>%
#     mutate(E = map_dbl(eig, ~ max(Re(.x)))) %>%
#     select(eta, E) %>%
#     mutate(E = E - 1) %>%
#     arrange(desc(E))


tradeoffs_outcomes <- list()
tradeoffs_outcomes[["points"]] <- tradeoff_sim_df %>%
    filter(eta != "additive") %>%
    group_by(eta) %>%
    filter(unq_spp_filter(V1, V2, .prec = 0.001)) %>%
    ungroup()
tradeoffs_outcomes[["lines"]] <- tradeoff_sim_df %>%
    filter(eta == "additive") %>%
    mutate(radius = sqrt(V1^2 + V2^2)) %>%
    group_by(eta) %>%  # this keeps the variable around for facetting
    summarize(radius = median(radius), .groups = "drop") %>%
    mutate(V1 = map(radius,
                    ~ .x * sin(seq(0, 0.5, length.out = 1001) *
                                   pi)),
           V2 = map(radius,
                    ~ .x * cos(seq(0, 0.5, length.out = 1001) *
                                   pi))) %>%
    unnest(cols = c(V1, V2))





#' Perturb one species' investments and output the change in investments
#' that would occur.
#' This is used to visualize stable points or neutrally stable arcs.
#' `.dd` should be tibble that looks like the following:
#     eta              N    V1    V2
#     <fct>        <dbl> <dbl> <dbl>
#   1 sub-additive 3709.  1.06  1.06
#   2 sub-additive 3709.  1.06  1.06
#'
#' `.add_var` is additive genetic variance. `NULL` or `NA` results in
#' the default for `quant_gen` being used.
#' Default to 1, which effectively just outputs the selection strengths.
#'
perturb_invests <- function(.dd, .add_var = 1) {

    if (is.null(.add_var) | is.na(.add_var)) {
        .add_var <- with(list(n = 2), eval(formals(quant_gen)$add_var))
    }

    .eta <- case_when(grepl("^sub", .dd$eta[1]) ~ -0.6,
                      grepl("^super", .dd$eta[1]) ~ 0.6,
                      TRUE ~ 0)
    C <- matrix(0, 2, 2)
    C[lower.tri(C)] <- .eta
    C <- C + t(C)
    diag(C) <- 1

    OUT <- tibble(v1 = c(0.5, 1.5, 1.5),
                  v2 = c(0.2, 0.5, 1.0)) %>%
        bind_rows(rename(., v1 = v2, v2 = v1)) %>%
        pmap_dfr(function(v1, v2) {
            # v1 = 0.25; v2 = 0.5
            # rm(v1, v2, .diffs, u, V, ss)

            # Change the species whose investments are most similar:
            .diffs <- ((.dd$V1 - v1)^2 + (.dd$V2 - v2)^2) / 2
            u <- which(.diffs == min(.diffs))[1]

            V <- select(.dd, starts_with("V"))
            V$V1[u] <- v1
            V$V2[u] <- v2

            ss <- sauron:::sel_str_cpp(V = t(V),
                                       N = .dd$N,
                                       f = formals(quant_gen)$f,
                                       a0 = formals(quant_gen)$a0,
                                       C = C,
                                       r0 = formals(quant_gen)$r0,
                                       D = diag(c(-1e-6, 1e-6)))
            dV <- .add_var * ss[,u]

            tibble(V1 = v1, V2 = v2, dV1 = dV[1], dV2 = dV[2])
        }) %>%
        mutate(eta = .dd$eta[[1]],
               newV1 = V1 + dV1,
               newV2 = V2 + dV2) %>%
        select(eta, V1, V2, newV1, newV2, everything())

    return(OUT)

}

#' I've already checked that different reps result in the same community,
#' which is why I'm just perturbing rep # 1.
tradeoffs_outcomes[["arrows"]] <- tradeoff_sim_df %>%
    filter(rep == 1) %>%
    select(-rep) %>%
    split(.$eta) %>%
    map_dfr(perturb_invests)





tradeoffs_outcomes_p <- crossing(v1 = round(seq(0, 2.5, 0.01), 2), v2 = v1,
                                   eta = c(-1,0,1)*0.6) %>%
    # `f` is fitness without competition
    mutate(f = exp(formals(quant_gen)$r0 - formals(quant_gen)$f *
                       (v1^2 + 2 * eta * v1 * v2 + v2^2)),
           # `f0` is fitness without competition or cost to investing
           f0 = exp(formals(quant_gen)$r0),
           # proportional reduction from the maximum fitness without competition
           M = (f0 - f) / f0) %>%
    mutate(eta = factor(eta, levels = sort(unique(eta)),
                        labels = paste0(c("sub-", "", "super-"), "additive"))) %>%
    # group_by(eta) %>%
    # summarize(M = mean(M))
    ggplot(aes(v1, v2)) +
    geom_raster(aes(fill = M), interpolate = FALSE) +
    geom_abline(slope = 1, intercept = 0, linetype = 2, color = "gray70") +
    geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
    geom_vline(xintercept = 0, color = "gray70", size = 0.5) +
    facet_wrap(~ eta) +
    scale_x_continuous("Conflict investment", breaks = 0:2) +
    scale_y_continuous("Partitioning investment", breaks = 0:2) +
    scale_fill_viridis("Proportional\nfitness cost", option = "mako",
                       begin = 0.1, direction = -1) +
    coord_equal(xlim = c(-0.05, 2.5), ylim = c(-0.05, 2.5)) +
    theme(strip.text.x = element_text(size = 10, margin = margin(0,0,0,b=6)),
          strip.text.y = element_text(size = 10, margin = margin(0,0,0,l=6)),
          legend.text = element_text(size = 8),
          legend.title = element_text(size = 9, lineheight = 0.75),
          legend.key.height = unit(1, "lines"),
          legend.position = "right") +
    geom_path(data = tradeoffs_outcomes[["lines"]],
              aes(V1, V2), color = "black", size = 1, linetype = "23") +
    geom_point(data = tradeoffs_outcomes[["points"]],
               aes(V1, V2), color = "black", size = 4, shape = 19) +
    geom_point(data = tradeoffs_outcomes[["points"]],
               aes(V1, V2), size = 4, shape = 1) +
    geom_segment(data = tradeoffs_outcomes[["arrows"]],
                 aes(V1, V2, xend = newV1, yend = newV2),
                 arrow = arrow(length = unit(0.1, "lines"), type = "closed"),
                 linejoin = "mitre", size = 1) +
    NULL


# If you want to add letters:
# tradeoffs_outcomes_p <- function() {
#     grid.newpage()
#     grid.draw(tradeoffs_outcomes_p_main)
#     grid.draw(textGrob("A", x = 0.08, y = 0.99, just = c("left", "top")))
#     grid.draw(textGrob("B", x = 0.337, y = 0.99, just = c("left", "top")))
#     grid.draw(textGrob("C", x = 0.595, y = 0.99, just = c("left", "top")))
#     invisible(NULL)
# }
#
# tradeoffs_outcomes_p()


if (.RESAVE_PLOTS) save_plot(tradeoffs_outcomes_p, 6, 2.3, .prefix = "2-")