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



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

# *--------------------------* ----
# Invest ~ Omega ----

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



#'
#' Total investment [`sqrt(x^2 + p^2)`] at equilibrium based on
#' eta and equilibrium Omega
#'
calc_invest <- function(Omega, eta) {
    # Below, rounding is to prevent R from making mistakes near boundaries
    # As an example, on my machine, if I run `(0.1 * 0.4 / 1e-4) - 400`,
    # I get `5.684342e-14`.
    stopifnot(length(eta) %in% c(1, length(Omega)))
    if (length(eta) == 1) eta <- rep(eta, length(Omega))
    a0 <- formals(quant_gen)$a0
    f <- formals(quant_gen)$f
    invests <- numeric(length(Omega))

    # Fill in non-zero values for sub-additive
    non_NA_sub <- eta < 0 & Omega >= round(f * (1 + eta) / a0, 12)
    logs_sub <- round(log(a0 * Omega[non_NA_sub] / (f * (1 + eta[non_NA_sub]))), 15)
    invests[non_NA_sub] <- sqrt(logs_sub)
    # ... and for additive / super-additive
    non_NA_add_super <- eta >= 0 & Omega >= round(f/a0, 12)
    logs_add_super <- round(log(a0 * Omega[non_NA_add_super] / f), 15)
    invests[non_NA_add_super] <- sqrt(logs_add_super)

    return(invests)
}


# -------------------*
# Heatmap
# -------------------*
omega_inv_hm_df <- crossing(Omega = seq(0, 10e3, length.out = 101),
                         eta = seq(-0.8, 0.8, length.out = 101) %>% round(5)) %>%
    mutate(inv = calc_invest(Omega, eta))
omega_inv_hm_p <- omega_inv_hm_df %>%
    ggplot(aes(Omega, eta)) +
    geom_raster(aes(fill = inv)) +
    geom_path(data = omega_inv_hm_df %>%
                  filter(inv == 0) %>%
                  mutate(O_incr = Omega %>% unique() %>% sort() %>% diff() %>% .[[1]],
                         e_incr = eta %>% unique() %>% sort() %>% diff() %>% .[[1]],
                         Omega = Omega + O_incr / 2) %>%
                  select(-O_incr, -inv) %>%
                  group_by(eta) %>%
                  filter(Omega == max(Omega)) %>%
                  ungroup() %>%
                  group_by(Omega) %>%
                  filter(eta %in% range(eta)) %>%
                  mutate(eta = ifelse(eta == min(eta),
                                      eta - e_incr / 2,
                                      eta + e_incr / 2)) %>%
                  ungroup(),
              color = "white", size = 0.75) +
    geom_text(data = tibble(Omega = 480, eta = 0.35, lab = "zero investment"),
              aes(label = lab), angle = 90, size = 10 / 2.83465, vjust = 0.5,
              color = "white") +
    scale_y_continuous(expression("Cost non-additivity (" * eta * ")")) +
    scale_x_continuous(expression("Equilibrium effective competitive neighborhood (" *
                                      hat(Omega)[i] * ")")) +
    scale_fill_viridis_c("Total\ninvestment", option = "A") +
    coord_fixed(ratio = diff(range(omega_inv_hm_df$Omega)) /
                    diff(range(omega_inv_hm_df$eta)),
                xlim = range(omega_inv_hm_df$Omega),
                ylim = range(omega_inv_hm_df$eta)) +
    NULL


# -------------------*
# Lines
# -------------------*

oil_etas <- seq(-0.8, 0, length.out = 5) %>% round(1)
omega_inv_line_df <- crossing(Omega = seq(0, 10e3, length.out = 1001),
                         eta = oil_etas) %>%
    mutate(inv = calc_invest(Omega, eta),
           eta_fct = ifelse(eta < 0, paste("eta ==", eta), "eta >= 0") %>%
               factor(levels = c(paste("eta ==", head(oil_etas, -1)),
                                 "eta >= 0")))


omega_inv_line_p <- omega_inv_line_df %>%
    ggplot(aes(Omega, inv, color = eta_fct)) +
    geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
    geom_line() +
    scale_y_continuous("Total investment") +
    scale_x_continuous(expression("Equilibrium effective competitive neighborhood (" *
                                      hat(Omega)[i] * ")")) +
    scale_color_viridis_d(NULL, option = "D", end = 0.8, direction = -1,
                          labels = parse_format()) +
    theme(legend.text.align = 0) +
    NULL

omega_inv_p <- omega_inv_hm_p + omega_inv_line_p +
    plot_annotation(tag_levels = "A") +
    plot_layout(ncol = 1, heights = c(1, 0.9))


save_plot(omega_inv_p, 6, 6, .prefix = "S2-")