_results/fig1.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 1: Partitioning vs conflict ----

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

#'
#' This part has simulations showing the difference between partitioning
#' vs agonistic traits.
#' They include just 1 suite of traits at a time.
#'


one_axis_one_fixed_spp <- function(.d, .vj, .Nj,
                                   .vi_0 = 0.05, .Ni_0 = 2500) {

    .max_t <- length(.vj)
    stopifnot(.max_t == length(.Nj))

    add_var <- with(list(n = 1), eval(formals(quant_gen)$add_var))

    # Changes through time:
    V = numeric(.max_t + 1)
    N = numeric(.max_t + 1)
    V[1] <- .vi_0
    N[1]  <- .Ni_0


    for (t in 1:.max_t) {
        ss_mat <- sauron:::sel_str_cpp(V = cbind(V[t], .vj[t]),
                                       N = c(N[t], .Nj[t]),
                                       f = formals(quant_gen)$f,
                                       a0 = formals(quant_gen)$a0,
                                       r0 = formals(quant_gen)$r0,
                                       C = matrix(1), D = matrix(.d))
        V[t+1] <- V[t] + add_var * ss_mat[1,1]

        fit <- sauron:::F_it_cpp(i = 0,
                                 V = cbind(V[t], .vj[t]),
                                 N = c(N[t], .Nj[t]),
                                 f = formals(quant_gen)$f,
                                 a0 = formals(quant_gen)$a0,
                                 r0 = formals(quant_gen)$r0,
                                 C = matrix(1), D = matrix(.d))
        N[t+1] <- N[t] * fit
    }

    tibble(d = .d, time = rep(0:.max_t, 2),
           spp = factor(rep(1:2, each = .max_t + 1),
                        levels = 1:2,
                        labels = sprintf("%s species",
                                         c("evolving", "fixed"))),
           N = c(N, c(.Nj, tail(.Nj, 1))),
           V = c(V, c(.vj, tail(.vj, 1)))) %>%
        arrange(time, spp)

}






# The fixed investments and abundances (from species "j") is based on the
# equilibriums for two species with these d values.
one_axis_sims <- tibble(.d = c(-1, 1),
                        .vj = list(c(rep(0, 10e3), rep(1, 3000))),
                        .Nj = list(rep(2500, 13e3))) %>%
    pmap_dfr(one_axis_one_fixed_spp) %>%
    filter(!(grepl("fixed", spp) & d > 0)) %>%
    mutate(plot_fct = case_when(grepl("fixed", spp) ~ 0L, d < 0 ~ 1L,
                                d > 0 ~ 2L, TRUE ~ NA_integer_) %>%
               factor(levels = 0:2,
                      labels = c("species A",
                                 "species B - conflict traits",
                                 "species B - partitioning traits")))

if (any(is.na(one_axis_sims$plot_fct))) {
    stop("\n~~~~~~~~~~~~~~\nERROR IN CREATING `plot_fct`\n~~~~~~~~~~~~~~\n")
}

vmod <- 3000

# Colors for the two species:
one_axis_pal <- c("dodgerblue3", "goldenrod")



one_axis_plot_list <- one_axis_sims %>%
    filter(time > 9.5e3) %>%
    # filter(time < 10e3 + 1250) %>%
    mutate(V = V * vmod) %>%
    pivot_longer(N:V) %>%
    mutate(name = factor(name, levels = c("N", "V"))) %>%
    split(interaction(.$plot_fct, .$name, sep = "___", drop = TRUE)) %>%
    map(function(.d) {
        # The code below makes the two lines not connect for species A
        # when it switches its investment value.
        if (grepl("A$", .d$plot_fct[1]) && .d$name[1] == "V") {
            .d <- bind_rows(.d,
                            filter(.d, value > 0, name == "V") %>%
                                filter(time == min(time)) %>%
                                mutate(value = NA, time = time - 0.5))
        }
        .color <- one_axis_pal[ifelse(grepl("A$", .d$plot_fct[1]), 1, 2)]
        .p <- ggplot(.d, aes(time - 10e3, value)) +
            ggtitle(paste(.d$plot_fct[1], "--", .d$name[1])) +
            # geom_hline(yintercept = 0, linetype = 1, color = "gray70", size = 0.5) +
            geom_vline(xintercept = 0, linetype = 2, color = "gray70", size = 0.5) +
            geom_line(size = 1, color = .color) +
            scale_x_continuous("Time", limits = c(-500, 3000)) +
            scale_y_continuous("Abundance", limits = c(0, 4200)) +
            coord_cartesian(expand = FALSE, clip = "off") +
            theme(strip.text = element_blank(),
                  axis.ticks = element_blank(),
                  axis.text = element_blank(),
                  axis.title = element_blank(),
                  plot.title = element_blank(),
                  plot.margin = margin(0, 0, b=1, l=1, unit = "mm"),
                  strip.background = element_blank()) +
            NULL
        return(.p)
    })

# wrap_plots(one_axis_plot_list, nrow = 3)




#'
#' Save separate plots in a folder, then edit the overall figure in Illustrator:
#'
if (.RESAVE_PLOTS) {

    if (!dir.exists("_results/results-plots/fig1_sub-plots")) {
        dir.create("_results/results-plots/fig1_sub-plots")
    }
    for (i in 1:length(one_axis_plot_list)) {
        fn <- names(one_axis_plot_list)[i] %>%
            str_replace("species ", "spp") %>%
            str_replace(" - ", "_") %>%
            str_replace("___", "-")
        if (nchar(fn) > 9) {
            fn <- fn %>%
                str_sub(start = rbind(c(1, 9), c(-2, -1))) %>%
                str_flatten()
        }
        fn <- paste0("fig1_sub-plots/1-1axis-", fn)
        gt <- ggplot_gtable(ggplot_build(one_axis_plot_list[[i]]))
        gt$layout$z[grepl("^panel", gt$layout$name)] <- max(gt$layout$z) + 1
        save_plot(function() {grid.newpage(); grid.draw(gt)},
                  2.5, 1.25, .name = fn)
    }; rm(i, fn)

}