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



#'
#' Palette for 4-species community.
#' It shows 4 distinct colors with the first two being blues and latter two
#' being orange/reds.
#' It's also ordered so that nearest communities will be contrasting.
#'
four_spp_pal <- c(rev(turbo(2, begin = 0.1, end = 0.3)),
                  turbo(2, begin = 0.75, end = 0.9))



#' Get x,y coordinates from angle and radius.
#'
#' If statements used to avoid round errors when `.angle %in% c(0,90)`
#'
make_xy <- function(.angle, .radius) {
    if (.angle == 0) {
        .x <- .radius
        .y <- 0
    } else if (.angle == 90) {
        .x <- 0
        .y <- .radius
    } else {
        .x <- cos(.angle * pi / 180) * .radius
        .y <- sin(.angle * pi / 180) * .radius
    }
    return(c(.x, .y))
}

#'
#' This function removes time points after a community reaches equilibrium.
#' This is useful to keep plot sizes down.
#'
remove_redundant_times <- function(.df, max_t = NULL) {

    .df  <- .df %>%
        arrange(time, spp)

    .df_ref <- .df %>%
        filter(time == max(time))

    #' This calculates maximum average square differences between
    #' `.df_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 `.df_ref` and `.df[.df$time == t,]`.
    #' Returns `Inf` if they have a different number of species.
    max_asd <- function(.t) {
        .df_t <- filter(.df, time == .t)
        if (nrow(.df_t) != nrow(.df_ref)) return(Inf)
        .asd <- ((.df_t$V1 - .df_ref$V1)^2 + (.df_t$V2 - .df_ref$V2)^2) / 2
        return(max(.asd))
    }
    for (t in sort(unique(.df$time))) {
        masd <- max_asd(t)
        if (masd < 1e-4^2) break
    }
    if (t < max(.df$time)) {
        .df <- .df %>% filter(time <= t)
    }

    return(.df)

}



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

# 2 communities, 2 invader permutations ----

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

#' Make matrix inside a list that contains starting investments based on
#' their angle from the origin. This is used for `V0` arg below.
make_V0 <- function(.angles) {
    rbind(cos(.angles * pi / 180) * 1,
          sin(.angles * pi / 180) * 1) %>%
        list()
}


res_angles <- 20 + c(-1, 1) * 2
inv_d_mag <- 0.6

#'
#' Note that the code below prints lambda, which should be == 1 for a
#' neutrally stable community.
#'
equil_comms <- map_dfr(list(res_angles, 90 - res_angles),
        function(.a) {
            .above <- all(.a > 45)
            stopifnot(.above || all(.a < 45))
            .comm_type <- ifelse(.above, "part", "conf")
            .X <- tibble(eta = 0, d1 = -1 * inv_d_mag, d2 = 1 * inv_d_mag,
                   final_t = 20e3L, spp_gap_t = 0L, save_every = 0L,
                   V0 = make_V0(.a)) %>%
                run_determ_1rep()
            L <- max(Re(.X[["eig"]][[1]]))
            print(L)
            .X %>%
                select(spp, N, V1, V2) %>%
                unnest(c(spp, N, V1, V2)) %>%
                mutate(time = 1L) %>%
                mutate(Omega = calc_omega(., d_vec = c(-1, 1) * inv_d_mag)) %>%
                mutate(comm = .comm_type,
                       invest = sqrt(V1^2 + V2^2)) %>%
                select(comm, everything(), -time)
        })


do_invasion <- function(.inv_angle, .res_comm_type,
                        .equil_comms = equil_comms, .radius = 1) {

    # .inv_angle = 4; .res_comm_type = "conf"; .radius = NULL;
    # .equil_comms = equil_comms
    # rm(.inv_angle, .res_comm_type, .inv_V, .radius, .V0, .N0, i_sim,
    #    equil_row, .equil_comms)

    .res_comm_type <- match.arg(.res_comm_type, c("conf", "part"))

    stopifnot(length(.inv_angle) == 1 && is.numeric(.inv_angle))
    stopifnot(.inv_angle <= 90 && .inv_angle >= 0)

    if (is.null(.radius)) {
        .radius <- .equil_comms %>%
            filter(comm == .res_comm_type) %>%
            .[["invest"]] %>%
            mean()
    }

    .inv_V <- cbind(make_xy(.inv_angle, .radius))

    .V0 <- .equil_comms %>%
        filter(comm == .res_comm_type) %>%
        .[,c("V1", "V2")] %>%
        as.matrix() %>%
        unname() %>%
        t() %>%
        cbind(.inv_V)
    .N0 <- c(filter(.equil_comms, comm == .res_comm_type)$N, 1)

    i_sim <- tibble(eta = 0, d1 = -1 * inv_d_mag, d2 = 1 * inv_d_mag,
                     final_t = 20e3L, spp_gap_t = 0L, save_every = 1L,
                     V0 = list(.V0), N0 = list(.N0)) %>%
        run_determ_1rep() %>%
        select(spp, time, N, V1, V2) %>%
        unnest(c(spp, time, N, V1, V2)) %>%
        arrange(time, spp) %>%
        mutate(Omega = calc_omega(., d_vec = c(-1, 1) * inv_d_mag)) %>%
        mutate(inv_angle = .inv_angle, res_comm = .res_comm_type,
               invest = sqrt(V1^2 + V2^2))

    # Add row of equilibrium community for plotting:
    equil_row <- .equil_comms %>%
        filter(comm == .res_comm_type) %>%
        mutate(time = -1L, inv_angle = .inv_angle) %>%
        rename(res_comm = comm)

    i_sim <- bind_rows(equil_row, i_sim) %>%
        select(res_comm, inv_angle, spp, everything())

    return(i_sim)

}


inv_sims <- crossing(inv_angle = c(5, 90 - 5), res_comm = c("part", "conf")) %>%
    as.list() %>%
    unname() %>%
    c(list(FUN = do_invasion, SIMPLIFY = FALSE)) %>%
    do.call(what = mapply) %>%
    do.call(what = bind_rows) %>%
    mutate(spp_int = spp %>% paste() %>% as.integer(),
           spp = factor(spp_int, levels = 1:max(spp_int),
                        labels = c(paste("resident", 1:(max(spp_int)-1)),
                                   "invader")),
           inv_type = ifelse(inv_angle < 45, "conf", "part"),
           sim_grp = sprintf("res=%s--inv=%s", res_comm, inv_type) %>%
               factor(levels = c("res=conf--inv=conf", "res=conf--inv=part",
                                 "res=part--inv=conf", "res=part--inv=part"))) %>%
    select(sim_grp, spp, time, N, V1, V2, Omega, invest)

# Should be 12:
# inv_sims %>% filter(time == max(time)) %>% nrow()






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

# Z-transform and plot ----

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


#'
#' Z-transforms of beginning and ending points plotted together:
#'
z_inv_sims <- inv_sims %>%
    filter(time == 0 | time == max(time)) %>%
    mutate(time = as.integer(time == 0) %>%
               factor(levels = 1:0, labels = c("begin", "end"))) %>%
    select(sim_grp, spp, time, N, Omega, invest) %>%
    #' Z-scoring values so they can be plotted next to each other.
    #' I'm removing N for invaders at beginning bc they all start at 1,
    #' and that would mess up the rest of the points for the `N` values.
    mutate(N = ifelse(time == "begin" & spp == "invader", NA, N),
           N = (N - mean(N, na.rm = TRUE)) / sd(N, na.rm = TRUE),
           invest = (invest - mean(invest)) / sd(invest),
           Omega = (Omega - mean(Omega)) / sd(Omega)) %>%
    pivot_longer(N:invest, names_to = "var", values_to = "val") %>%
    mutate(var = factor(var, levels = c("invest", "Omega", "N"))) %>%
    #' Replacing the N values we turned to `NA` previously to a very low value
    #' that we'll cut off in the plots.
    mutate(val = ifelse(is.na(val), -3, val)) %>%
    identity()




z_inv_p_list <-
    levels(z_inv_sims$sim_grp) %>%
    as.list() %>%
    set_names() %>%
    map(function(.sg) {
        .spp_adds <- c(0, 0.1, 0.35) %>%
            {. - mean(.)}
        .dd <- z_inv_sims %>%
            filter(sim_grp == .sg) %>%
            mutate(x = as.integer(var) + .spp_adds[as.integer(spp)]) %>%
            pivot_wider(names_from = time, values_from = val)
        .dd %>%
            ggplot(aes(x, begin, color = spp)) +
            geom_hline(yintercept = 0, color = "gray70", size = 0.5) +
            geom_segment(aes(xend = x, yend = end),
                         arrow = arrow(length = unit(0.2, "lines"), type = "closed"),
                         linejoin = "mitre") +
            # geom_segment(aes(x = x - 0.05, xend = x + 0.05, y = end, yend = end)) +
            geom_point(shape = 19, size = 2, alpha = 0.75) +
            geom_point(shape = 1, size = 2) +
            scale_color_manual(values = four_spp_pal[-3], guide = "none") +
            scale_shape_manual(values = c(1, 19), guide = "none") +
            scale_y_continuous("Z-score", breaks = -2:2) +
            scale_x_continuous(NULL, breaks = 1:3, limits = c(0.5, 3.5),
                               labels = parse(text = c("I[i]", "Omega[i]", "N[i]"))) +
            coord_cartesian(ylim = c(-2, 2)) +
            theme(panel.background = element_blank(),
                  axis.title = element_blank(),
                  axis.text.x = element_blank(),
                  plot.margin = margin(0,0,0,0))
    })
# wrap_plots(z_inv_p_list, nrow = 2)




#'
#' Plots of starting communities.
#'
inv_v_p_list <-
    levels(inv_sims$sim_grp) %>%
    as.list() %>%
    set_names() %>%
    map(function(.sg) {
        .dd <- inv_sims %>%
            filter(sim_grp == .sg) %>%
            filter(time == 0)
        .dd %>%
            ggplot(aes(V1, V2, color = spp))+
            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) +
            # geom_path() +
            # geom_point(aes(shape = time), size = 2) +
            geom_point(shape = 19, size = 2, alpha = 0.75) +
            geom_point(shape = 1, size = 2) +
            scale_color_manual(values = four_spp_pal[-3], guide = "none") +
            scale_shape_manual(values = c(1, 19), guide = "none") +
            scale_y_continuous("Partitioning investment", breaks = 0:2) +
            scale_x_continuous("Conflict investment", breaks = 0:2) +
            coord_equal(xlim = c(0, 2), ylim = c(0, 2)) +
            theme(panel.background = element_blank(),
                  axis.title = element_blank(),
                  axis.ticks = element_blank(),
                  axis.text.x = element_blank(),
                  axis.text.y = element_blank(),
                  #' bottom margin here makes these plots' x-axes
                  #' align with those from plots in `inv_ts_p_list`.
                  plot.margin = margin(0,0,0,b=5.5))
    })
# wrap_plots(inv_v_p_list, nrow = 2)




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

    if (!dir.exists("_results/results-plots/fig4_sub-plots")) {
        dir.create("_results/results-plots/fig4_sub-plots")
    }
    for (n in levels(z_inv_sims$sim_grp)) {
        save_plot(z_inv_p_list[[n]], 1.8, 1.5,
                  .name = sprintf("fig4_sub-plots/Z-%s", n))
        save_plot(inv_v_p_list[[n]], 0.75, 0.75,
                  .name = sprintf("fig4_sub-plots/V-%s", n))
    }; rm(n)

}