_results/_preamble.R

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

# start ----

# load packages
suppressPackageStartupMessages({
    library(sauron)
    library(tidyverse)
    library(grid)
    library(scales)
    library(parallel)
    library(viridis)
    library(patchwork)
})


# This sets plotting device on LAN computer:
if (file.exists(".Rprofile")) source(".Rprofile")


# number of threads to use for simulations
.N_THREADS <- max(1, parallel::detectCores()-2)
# `mclapply` automatically uses this many threads
options(mc.cores = .N_THREADS)


#'
#' Save plot using proper device.
#'
save_plot <- function(plot_obj, .width, .height,
                      .name = NULL, .prefix = NULL, .suffix = NULL,
                      .png = FALSE) {
    if (is.null(.name)) {
        fn <- gsub("_p$", "", paste(substitute(plot_obj)))
        if (!is.null(.prefix)) fn <- paste0(.prefix, fn)
        if (!is.null(.suffix)) fn <- paste0(fn, .suffix)
    } else fn <- gsub(".pdf$", "", .name)
    fn <- sprintf("_results/results-plots/%s.pdf", fn)
    message(fn)
    png_fn <- gsub(".pdf$", ".png", fn)
    if (.png) message("... and ", png_fn)
    if (inherits(plot_obj, "function")) {
        cairo_pdf(fn, width = .width, height = .height)
        plot_obj()
        dev.off()
        if (.png) {
            png(png_fn, .width, .height, "in", res = 300)
            plot_obj()
            dev.off()
        }
    } else {
        plot_fun <- ifelse(inherits(plot_obj, "egg"), print, plot)
        cairo_pdf(fn, width = .width, height = .height)
        plot_fun(plot_obj)
        dev.off()
        if (.png) {
            png(png_fn, .width, .height, "in", res = 300)
            plot_fun(plot_obj)
            dev.off()
        }
    }
    invisible(NULL)
}

unq_spp_filter <- function(..., .prec = 0.001) {
    V <- list(...)
    V <- do.call(cbind, V)
    V <- split(V, row(V))
    return(as.logical(sauron:::unq_spp_cpp(V, precision = .prec)))
}

group_spp <- function(..., .prec = 0.001) {
    V <- list(...)
    V <- do.call(cbind, V)
    V <- split(V, row(V))
    return(sauron:::group_spp_cpp(V, precision = .prec))
}


pivot <- function(.df) {
    if ("pheno" %in% colnames(.df)) {
        .df <- rename(.df, V = geno, Vp = pheno) %>%
            pivot_wider(names_from = axis, values_from = c(V, Vp),
                        names_sep = "")
    } else {
        .df <- mutate(.df, axis = paste0("V", axis)) %>%
            pivot_wider(names_from = axis, values_from = geno)
    }
    return(.df)
}


#'
#' Remove axis titles and/or text for plots arrayed into grids.
#' Assumes byrow = TRUE
#'
adjust_grid_axes <- function(x, i, rows, cols,
                             no_titles = "none",
                             keep_text = "none") {
    # i = 1:3; rows = 3; cols = 1; no_titles = "none"; keep_text = "none"
    # rm(i, rows, cols, no_titles, keep_text)
    no_titles <- match.arg(no_titles, c("none", "x", "y", "xy"))
    keep_text <- match.arg(keep_text, c("none", "x", "y", "xy"))
    if (i <= (rows-1) * cols || grepl("x", no_titles)) {
        x <- x + theme(axis.title.x = element_blank())
    }
    if ((i-1) %% cols != 0 || grepl("y", no_titles)) {
        x <- x + theme(axis.title.y = element_blank())
    }
    if (i <= (rows-1) * cols && !grepl("x", keep_text)) {
        x <- x + theme(axis.text.x = element_blank())
    }
    if ((i-1) %% cols != 0 && !grepl("y", keep_text)) {
        x <- x + theme(axis.text.y = element_blank())
    }
    return(x)
}


#'
#' Run 1 rep of deterministic simulations (with given starting trait values)
#' based on a single row of a dataframe.
#' `.d` must contain the columns `eta`, `d1`, `d2`, and `V0`.
#' All other columns must be within the arguments to `quant_gen`.
#' Arguments `n` and `n_reps` are not allowed to be changed.
#'
#' This is used in a few of the functions below.
#'
run_determ_1rep <- function(.d) {

    # .d <- tibble(eta = 0.6, d1 = -0.5, d2 = 0.5, V0 = list(cbind(c(1,0), c(0,1))))
    # rm(.d, .dd, args, other_args, sim0, n, .E, X, out)
    stopifnot(is.data.frame(.d) && nrow(.d) == 1)
    stopifnot(all(c("eta", "d1", "d2", "V0") %in% colnames(.d)))
    stopifnot(is.numeric(.d$d1) && .d$d1 < 0)
    stopifnot(is.numeric(.d$d2) && .d$d2 > 0)
    if ("n_reps" %in% colnames(.d)) {
        stop("\nThis function not designed for > 1 rep")
    }
    if ("n" %in% colnames(.d)) {
        stop("\nYou should adjust `n` using `V0`.")
    }

    .dd <- as.list(.d)
    if (is.list(.dd$V0)) .dd$V0 <- .dd$V0[[1]]
    stopifnot(is.matrix(.dd$V0) && nrow(.dd$V0) == 2)

    args <- list(q = 2, eta = .dd$eta, V0 = .dd$V0,
                 n = ncol(.dd$V0), d = c(.dd$d1, .dd$d2), n_reps = 1,
                 spp_gap_t = 0L,
                 final_t = 100e3L,
                 save_every = 0L,
                 sigma_V0 = 0,
                 show_progress = FALSE)
    other_args <- as.list(select(.d, -eta, -d1, -d2, -V0))
    if (length(other_args) > 0) {
        stopifnot(!is.null(names(other_args)) && all(names(other_args) != ""))
        stopifnot(all(names(other_args) %in% names(formals(quant_gen))))
        for (n in names(other_args)) args[[n]] <- unlist(other_args[[n]])
    }
    sim0 <- do.call(quant_gen, args)

    # stable?
    .E <- sim0 %>%
        jacobians() %>%
        .[[1]] %>%
        eigen(only.values = TRUE) %>%
        .[["values"]]

    X <- sim0 %>%
        .[["nv"]] %>%
        pivot() %>%
        select(-rep)

    out <- .d %>%
        mutate(n_spp = ncol(.dd$V0),
               spp = list(X$spp),
               N = list(X$N),
               V1 = list(X$V1),
               V2 = list(X$V2),
               eig = list(.E))
    if ("time" %in% colnames(X)) {
        out <- mutate(out, time = list(X$time)) %>%
            relocate(time, .after = spp)
    }

    return(out)

}




#'
#' This function returns a logical indicating whether a community has
#' reached equilibrium.
#' It checks the last `n_times` time points to verify that they all
#' have the same community as the final time point.
#'
#' It takes as input a data frame with columns for species, time, N, V1, and V2.
#'
equil_check <- function(.output, n_times = 100) {

    stopifnot(all(c("spp", "time", "N", "V1", "V2") %in% colnames(.output)))

    if (inherits(.output$spp, "list")) {
        .output <- .output %>%
            unnest(c( spp, time, N, V1, V2))
    }
    if (!"rep" %in% colnames(.output)) {
        .output$rep <- 1L
    }

    #' This calculates maximum average square differences between
    #' `.ref` (community at `time == max(time)`) and the community at time `.t`.
    #' Because `.df` is sorted by species, rows should correspond to
    #' the same species between `.ref` and `.df[.df$time == t,]`.
    #' Returns `Inf` if they have a different number of species.
    max_asd <- function(.t, .ref) {
        .df_t <- filter(.df, time == .t)
        if (any(is.na(.df_t$V1))) return(Inf)
        if (nrow(.df_t) != nrow(.ref)) return(Inf)
        .asd <- ((.df_t$V1 - .ref$V1)^2 + (.df_t$V2 - .ref$V2)^2) / 2
        return(max(.asd))
    }

    .df_list <- .output %>%
        split(.$rep) %>%
        map(~ select(.x, spp, time, N, V1, V2) %>%
                arrange(time, spp))

    .equil_reached <- logical(length(.df_list))
    for (i in 1:length(.df_list)) {
        .df  <- .df_list[[i]]
        .df_ref <- .df %>% filter(time == max(time))
        .check_times <- sort(unique(.df$time))
        .check_times <- .check_times[.check_times < max(.df$time)]
        .check_times <- tail(.check_times,  n_times)
        # Comparisons for each time:
        .comparisons <- map_dbl(.check_times, max_asd, .ref = .df_ref)
        # Are all below the threshold?
        .equil_reached[i] <- all(.comparisons < 1e-4^2)
    }

    return(.equil_reached)
}