Nothing
R/population_dispersal.R
In poems: Pattern-Oriented Ensemble Modeling System
Defines functions
population_dispersal
Documented in
population_dispersal
#' Nested functions for population dispersal.
#'
#' Modular functions for the population simulator for performing dispersal of stage
#' abundance at a specified time step via dispersal rates provided.
#'
#' @examples
#' # User-defined dispersal: one-quarter of dispersing stages move one population over
#' simulator <- SimulatorReference$new()
#' example_function <- function(params) {
#' params$simulator$attached$params <- params # attach to reference object
#' emigrants <- round(params$stage_abundance * params$dispersal_stages * 0.25)
#' return(params$stage_abundance - emigrants + emigrants[, c(7, 1:6)])
#' }
#' dispersal_function <- population_dispersal(
#' replicates = 4,
#' time_steps = 10,
#' years_per_step = 1,
#' populations = 7,
#' demographic_stochasticity = TRUE,
#' density_stages = c(0, 1, 1),
#' dispersal = example_function,
#' dispersal_stages = c(0, 1, 0.5),
#' dispersal_source_n_k = list(cutoff = -0.5, threshold = 1.5),
#' dispersal_target_k = 5,
#' dispersal_target_n = list(threshold = 10, cutoff = 15),
#' simulator = simulator
#' )
#' carrying_capacity <- rep(10, 7)
#' stage_abundance <- matrix(
#' c(
#' 7, 13, 0, 26, 0, 39, 47,
#' 2, 0, 6, 8, 0, 12, 13,
#' 0, 3, 4, 6, 0, 9, 10
#' ),
#' nrow = 3,
#' ncol = 7,
#' byrow = TRUE
#' )
#' occupied_indices <- (1:7)[-5]
#' dispersal_function(
#' r = 2, tm = 6, carrying_capacity, stage_abundance,
#' occupied_indices
#' )
#'
#' @param replicates Number of replicate simulation runs.
#' @param time_steps Number of simulation time steps.
#' @param years_per_step Number of years per time step.
#' @param populations Number of populations.
#' @param demographic_stochasticity Boolean for optionally choosing demographic stochasticity for the transformation.
#' @param density_stages Array of booleans or numeric (0,1) for each stage to indicate which stages are affected by density.
#' @param dispersal Either a matrix of dispersal rates between populations (source columns to target rows) or a list of data frames of non-zero dispersal rates and indices for constructing a compact dispersal matrix, and optional changing rates over time (as per class \code{\link{DispersalGenerator}} \emph{dispersal_data} attribute). Alternatively a user-defined function (optionally nested in a list with additional attributes) may be used: \code{function(params)}, where \emph{params} is a list passed to the function containing:
#' \describe{
#' \item{\code{replicates}}{Number of replicate simulation runs.}
#' \item{\code{time_steps}}{Number of simulation time steps.}
#' \item{\code{years_per_step}}{Number of years per time step.}
#' \item{\code{populations}}{Number of populations.}
#' \item{\code{stages}}{Number of life cycle stages.}
#' \item{\code{demographic_stochasticity}}{Boolean for optionally choosing
#' demographic stochasticity for the transformation.}
#' \item{\code{density_stages}}{Array of booleans or numeric (0,1) for each
#' stage to indicate which stages are affected by density.}
#' \item{\code{dispersal_stages}}{Array of relative dispersal (0-1) for
#' each stage to indicate the degree to which each stage participates in
#' dispersal. This factor modifies dispersal proportion, not dispersal rate.}
#' \item{\code{dispersal_source_n_k}}{Dispersal proportion (p) density dependence via source population abundance divided by carrying capacity (n/k), where p is reduced via a linear slope (defined by two list items) from n/k <= \emph{cutoff} (p = 0) to n/k >= \emph{threshold}.}
#' \item{\code{dispersal_target_k}}{Dispersal rate (r) density dependence via target population carrying capacity (k), where r is reduced via a linear slope (through the origin) when k <= \emph{threshold}.}
#' \item{\code{dispersal_target_n}}{Dispersal rate (r) density dependence via target population abundance (n), where r is reduced via a linear slope (defined by two list items) from n >= \emph{threshold} to n <= \emph{cutoff} (r = 0) or vice versa.}
#' \item{\code{dispersal_target_n_k}}{Dispersal rate (r) density dependence via target population abundance divided by carrying capacity (n/k), where r is reduced via a linear slope (defined by two list items) from n/k >= \emph{threshold} to n/k <= \emph{cutoff} (r = 0) or vice versa.}
#' \item{\code{r}}{Simulation replicate.}
#' \item{\code{tm}}{Simulation time step.}
#' \item{\code{carrying_capacity}}{Array of carrying capacity values for each population at time step.}
#' \item{\code{stage_abundance}}{Matrix of abundance for each stage (rows) and population (columns) at time step.}
#' \item{\code{occupied_indices}}{Array of indices for populations occupied at time step.}
#' \item{\code{simulator}}{\code{\link{SimulatorReference}} object with dynamically accessible \emph{attached} and \emph{results} lists.}
#' \item{\code{additional attributes}}{Additional attributes when the transformation is optionally nested in a list.}
#' }
#' returns the post-dispersal abundance matrix
#' @param dispersal_stages Array of relative dispersal (0-1) for each stage to
#' indicate the degree to which each stage participates in dispersal (default
#' is 1 for all stages). This factor modifies dispersal proportion, not
#' dispersal rate.
#' @param dispersal_source_n_k Dispersal proportion (p) density dependence via
#' source population abundance divided by carrying capacity (n/k), where p is
#' reduced via a linear slope (defined by two list items) from n/k <= \emph{cutoff} (p = 0) to n/k >= \emph{threshold} or vice versa.
#' @param dispersal_target_k Dispersal rate (r) density dependence via target population carrying capacity (k), where r is reduced via a linear slope (through the origin) when k <= \emph{threshold}.
#' @param dispersal_target_n Dispersal rate (r) density dependence via target population abundance (n), where r is reduced via a linear slope (defined by two list items) from n >= \emph{threshold} to n <= \emph{cutoff} (r = 0) or visa-versa.
#' @param dispersal_target_n_k Dispersal rate (r) density dependence via target population abundance divided by carrying capacity (n/k), where r is reduced via a linear slope (defined by two list items) from n/k >= \emph{threshold} to n/k <= \emph{cutoff} (r = 0) or vice versa.
#' @param simulator \code{\link{SimulatorReference}} object with dynamically accessible \emph{attached} and \emph{results} lists.
#' @return Dispersal function: \code{function(r, tm, carrying_capacity, stage_abundance, occupied_indices)}, where:
#' \describe{
#' \item{\code{r}}{Simulation replicate.}
#' \item{\code{tm}}{Simulation time step.}
#' \item{\code{carrying_capacity}}{Array of carrying capacity values for each population at time step.}
#' \item{\code{stage_abundance}}{Matrix of abundance for each stage (rows) and population (columns) at time step.}
#' \item{\code{occupied_indices}}{Array of indices for populations occupied at time step.}
#' \item{\code{returns}}{New stage abundance matrix with dispersal applied.}
#' }
#' @export population_dispersal
population_dispersal <- function(
replicates,
time_steps,
years_per_step,
populations,
demographic_stochasticity,
density_stages,
dispersal,
dispersal_stages,
dispersal_source_n_k = NULL,
dispersal_target_k = NULL,
dispersal_target_n = NULL,
dispersal_target_n_k = NULL,
simulator
) {
if (is.null(dispersal)) {
# no dispersal
return(NULL)
}
# User-defined function?
if (
is.function(dispersal) ||
(is.list(dispersal) &&
length(which(unlist(lapply(dispersal, is.function)))))
) {
# Derive stages
stages <- length(dispersal_stages)
# Unpack transformation function and additional attributes from a list
additional_attributes <- list()
if (is.list(dispersal)) {
function_index <- which(unlist(lapply(dispersal, is.function)))
additional_attributes <- dispersal[-function_index]
dispersal <- dispersal[[function_index]]
}
if (is.function(dispersal)) {
# user-defined function
# List of parameters to pass to the user-defined function
params <- c(
list(
replicates = replicates,
time_steps = time_steps,
years_per_step = years_per_step,
populations = populations,
stages = stages,
demographic_stochasticity = demographic_stochasticity,
density_stages = density_stages,
dispersal_stages = dispersal_stages,
dispersal_source_n_k = dispersal_source_n_k,
dispersal_target_k = dispersal_target_k,
dispersal_target_n = dispersal_target_n,
dispersal_target_n_k = dispersal_target_n_k,
simulator = simulator
),
additional_attributes
)
## Create a nested function for applying user-defined dispersal of stage abundance ##
user_defined_function <- function(
r,
tm,
carrying_capacity,
stage_abundance,
occupied_indices
) {
# Add attributes to be made available to the user-defined function
params$r <- r
params$tm <- tm
params$carrying_capacity <- carrying_capacity
params$stage_abundance <- stage_abundance
params$occupied_indices <- occupied_indices
# Run user-defined dispersal function
tryCatch(
{
stage_abundance <- dispersal(params)
},
error = function(e) {
stop(
paste(
"Error produced within user-defined dispersal function:",
as.character(e)
),
call. = FALSE
)
}
)
# Warn if any negative or non-finite
if (any(!is.finite(stage_abundance))) {
warning(
"Non-finite abundances returned by user-defined dispersal function",
call. = FALSE
)
}
if (any(stage_abundance[which(is.finite(stage_abundance))] < 0)) {
warning(
"Negative abundances returned by user-defined dispersal function",
call. = FALSE
)
}
return(stage_abundance)
}
return(user_defined_function)
}
}
# Assumed to be matrix or generated dispersal data
# Initialize reusable dispersal attributes
if (is.matrix(dispersal)) {
# create compact matrix data
# Calculate the indices of non-zero dispersals
dispersal_data <- data.frame(which(dispersal > 0, arr.ind = TRUE))
dispersal_data <- dispersal_data[
order(dispersal_data[, 2], dispersal_data[, 1]),
]
names(dispersal_data) <- c("target_pop", "source_pop")
# Calculate indices for constructing compacted dispersal matrices for emigrants and immigrants
dispersal_rows <- tabulate(dispersal_data$source_pop, nbins = populations)
dispersal_cols <- tabulate(dispersal_data$target_pop, nbins = populations)
dispersal_compact_rows <- max(dispersal_rows, dispersal_cols)
compact_emigrant_matrix <- array(
1:dispersal_compact_rows,
c(dispersal_compact_rows, populations)
)
compact_immigrant_matrix <- compact_emigrant_matrix *
(compact_emigrant_matrix <=
matrix(
dispersal_cols,
nrow = dispersal_compact_rows,
ncol = populations,
byrow = TRUE
))
compact_emigrant_matrix <- compact_emigrant_matrix *
(compact_emigrant_matrix <=
matrix(
dispersal_rows,
nrow = dispersal_compact_rows,
ncol = populations,
byrow = TRUE
))
# Map the row of each compact matrix to the original target (for emigrants) or source (for immigrants) populations
dispersal_data$emigrant_row <- which(
compact_emigrant_matrix > 0,
arr.ind = TRUE,
useNames = FALSE
)[, 1]
dispersal_data$immigrant_row <- which(
compact_immigrant_matrix > 0,
arr.ind = TRUE,
useNames = FALSE
)[, 1]
target_sorted_indices <- dispersal_data[
order(dispersal_data$target_pop, dispersal_data$source_pop),
c("target_pop", "source_pop")
]
dispersal_data$immigrant_row <- dispersal_data$immigrant_row[order(
target_sorted_indices$source_pop,
target_sorted_indices$target_pop
)]
# Add dispersal rates
dispersal_data$dispersal_rate <- dispersal[as.matrix(dispersal_data[c(
"target_pop",
"source_pop"
)])]
dispersals_change_over_time <- FALSE
# Release variables from memory
dispersal_rows <- NULL
dispersal_cols <- NULL
compact_emigrant_matrix <- NULL
compact_immigrant_matrix <- NULL
target_sorted_indices <- NULL
} else if (
is.list(dispersal) && is.data.frame(dispersal[[1]]) && nrow(dispersal[[1]])
) {
# compact matrix data in list (as per DispersalModel class)
# Unpack dispersal data and determine compact matrix dimensions
dispersal_data <- dispersal[[1]]
dispersal_compact_rows <- max(dispersal_data[, c(
"emigrant_row",
"immigrant_row"
)])
# Are dispersals changing over time?
dispersals_change_over_time <- (length(dispersal) > 1)
if (dispersals_change_over_time) {
dispersal_data_changes <- dispersal
dispersal_data_changes[[1]] <- dispersal_data_changes[[1]][NULL, ]
}
} else {
return(NULL)
}
# Release dispersal from memory
dispersal <- NULL
# Create a compact matrix of dispersal rates
dispersal_compact_matrix <- array(0, c(dispersal_compact_rows, populations))
dispersal_compact_matrix[as.matrix(dispersal_data[, c(
"emigrant_row",
"source_pop"
)])] <- dispersal_data$dispersal_rate
# Does dispersal depend on source population abundance N divided by carrying capacity K?
dispersal_depends_on_source_pop_n_k <- (is.list(dispersal_source_n_k) &&
(is.numeric(dispersal_source_n_k$cutoff) ||
is.numeric(dispersal_source_n_k$threshold)))
# Does dispersal depend on target population carrying capacity K, abundance N, or N/K?
dispersal_depends_on_target_pop_k <- is.numeric(dispersal_target_k)
dispersal_depends_on_target_pop_n <- (is.list(dispersal_target_n) &&
(is.numeric(dispersal_target_n$threshold) ||
is.numeric(dispersal_target_n$cutoff)))
dispersal_depends_on_target_pop_n_k <- (is.list(dispersal_target_n_k) &&
(is.numeric(dispersal_target_n_k$threshold) ||
is.numeric(dispersal_target_n_k$cutoff)))
# Setup density dependence dispersal parameters
if (dispersal_depends_on_source_pop_n_k) {
# Convert NULL to zero in source N/K cutoff or one in threshold
if (dispersal_depends_on_source_pop_n_k) {
if (is.null(dispersal_source_n_k$cutoff)) dispersal_source_n_k$cutoff <- 0
if (is.null(dispersal_source_n_k$threshold))
dispersal_source_n_k$threshold <- 1
}
# Check threshold > cutoff
if (dispersal_source_n_k$threshold <= dispersal_source_n_k$cutoff) {
dispersal_depends_on_source_pop_n_k <- FALSE
warning(
"Dispersal density dependence for source N/K threshold must be greater than cutoff => not used",
call. = FALSE
)
}
}
if (
dispersal_depends_on_target_pop_k ||
dispersal_depends_on_target_pop_n ||
dispersal_depends_on_target_pop_n_k
) {
if (dispersal_depends_on_target_pop_n) {
# Convert NULL to zero in target N threshold or cutoff
if (is.null(dispersal_target_n$threshold))
dispersal_target_n$threshold <- 0
if (is.null(dispersal_target_n$cutoff)) dispersal_target_n$cutoff <- 0
}
if (dispersal_depends_on_target_pop_n_k) {
# Convert NULL to zero in target N/K threshold or cutoff
if (is.null(dispersal_target_n_k$threshold))
dispersal_target_n_k$threshold <- 0
if (is.null(dispersal_target_n_k$cutoff)) dispersal_target_n_k$cutoff <- 0
}
# Create a map of compact array indices for mapping dispersers (emigrants) to target populations
dispersal_target_pop_map <- array(0, c(dispersal_compact_rows, populations))
dispersal_target_pop_map[as.matrix(dispersal_data[, c(
"emigrant_row",
"source_pop"
)])] <- dispersal_data$target_pop
}
# Create a map of compact array indices for mapping dispersers (emigrants) to immigrants
dispersal_compact_indices <- array(
1:(dispersal_compact_rows * populations),
c(dispersal_compact_rows, populations)
)
dispersal_immigrant_map <- array(0, c(dispersal_compact_rows, populations))
dispersal_immigrant_map[as.matrix(dispersal_data[, c(
"emigrant_row",
"source_pop"
)])] <- dispersal_compact_indices[as.matrix(dispersal_data[, c(
"immigrant_row",
"target_pop"
)])]
# Release variables from memory
dispersal_data <- NULL
dispersal_compact_indices <- NULL
## Create a nested function for performing dispersal ##
dispersal_function <- function(
r,
tm,
carrying_capacity,
stage_abundance,
occupied_indices
) {
# Calculate occupied population number
occupied_populations <- length(occupied_indices)
# Apply any spatio-temporal dispersal changes
dispersal_compact_matrix_tm <- simulator$attached$dispersal_compact_matrix_tm
if (tm == 1 || !dispersals_change_over_time) {
dispersal_compact_matrix_tm <- dispersal_compact_matrix
} else if (
dispersals_change_over_time && nrow(dispersal_data_changes[[tm]])
) {
# and tm > 1
dispersal_compact_matrix_tm[as.matrix(dispersal_data_changes[[tm]][, c(
"emigrant_row",
"source_pop"
)])] <- dispersal_data_changes[[tm]]$dispersal_rate
}
simulator$attached$dispersal_compact_matrix_tm <- dispersal_compact_matrix_tm
# Select dispersals for occupied populations
occupied_dispersals <- dispersal_compact_matrix_tm[, occupied_indices]
# Calculate density abundance
if (
dispersal_depends_on_source_pop_n_k ||
dispersal_depends_on_target_pop_n ||
dispersal_depends_on_target_pop_n_k
) {
density_abundance <- .colSums(
stage_abundance * as.numeric(density_stages),
m = length(density_stages),
n = populations
)
}
# Modify dispersal rates when dispersal depends on source population N/K
if (dispersal_depends_on_source_pop_n_k) {
# Density dependent multipliers
dd_multipliers <- array(1, populations)
# Calculate the source N/K multipliers
abundance_on_capacity <- density_abundance / carrying_capacity
dd_multipliers[which(
abundance_on_capacity <= dispersal_source_n_k$cutoff
)] <- 0
modify_pop_indices <- which(
carrying_capacity > 0 &
dd_multipliers > 0 &
abundance_on_capacity < dispersal_source_n_k$threshold
)
dd_multipliers[modify_pop_indices] <- ((abundance_on_capacity[
modify_pop_indices
] -
array(dispersal_source_n_k$cutoff, populations)[modify_pop_indices]) /
array(
dispersal_source_n_k$threshold - dispersal_source_n_k$cutoff,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
# Apply modifying multipliers to dispersals
occupied_dispersals <- (occupied_dispersals *
matrix(
dd_multipliers[occupied_indices],
nrow = dispersal_compact_rows,
ncol = occupied_populations,
byrow = TRUE
))
} # dispersal depends on source pop N/K?
# Select occupied dispersal non-zero indices
occupied_dispersal_indices <- which(as.logical(occupied_dispersals)) # > 0
# Modify dispersal rates when dispersal depends on target population K, N, or N/K
if (
dispersal_depends_on_target_pop_k ||
dispersal_depends_on_target_pop_n ||
dispersal_depends_on_target_pop_n_k
) {
# Density dependent multipliers
dd_multipliers <- array(1, populations)
# Calculate the (below-threshold) target K multipliers
if (dispersal_depends_on_target_pop_k) {
modify_pop_indices <- which(carrying_capacity < dispersal_target_k)
dd_multipliers[modify_pop_indices] <- (carrying_capacity[
modify_pop_indices
] /
array(dispersal_target_k, populations)[modify_pop_indices])
}
# Calculate the target N multipliers
if (dispersal_depends_on_target_pop_n) {
if (all(dispersal_target_n$threshold < dispersal_target_n$cutoff)) {
# overcrowded cell avoidance \
dd_multipliers[which(
density_abundance >= dispersal_target_n$cutoff
)] <- 0
modify_pop_indices <- which(
density_abundance > dispersal_target_n$threshold &
dd_multipliers > 0
)
dd_multipliers[modify_pop_indices] <- ((array(
dispersal_target_n$cutoff,
populations
)[modify_pop_indices] -
density_abundance[modify_pop_indices]) /
array(
dispersal_target_n$cutoff - dispersal_target_n$threshold,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
} else if (
all(dispersal_target_n$threshold > dispersal_target_n$cutoff)
) {
# seek company /
dd_multipliers[which(
density_abundance <= dispersal_target_n$cutoff
)] <- 0
modify_pop_indices <- which(
density_abundance < dispersal_target_n$threshold &
dd_multipliers > 0
)
dd_multipliers[modify_pop_indices] <- ((density_abundance[
modify_pop_indices
] -
array(dispersal_target_n$cutoff, populations)[modify_pop_indices]) /
array(
dispersal_target_n$threshold - dispersal_target_n$cutoff,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
}
}
# Calculate the target N/K multipliers
if (dispersal_depends_on_target_pop_n_k) {
dd_multipliers[which(carrying_capacity <= 0)] <- 0
abundance_on_capacity <- density_abundance / carrying_capacity
if (all(dispersal_target_n_k$threshold < dispersal_target_n_k$cutoff)) {
# overcrowded cell avoidance \
dd_multipliers[which(
abundance_on_capacity >= dispersal_target_n_k$cutoff
)] <- 0
modify_pop_indices <- which(
abundance_on_capacity > dispersal_target_n_k$threshold &
dd_multipliers > 0
)
dd_multipliers[modify_pop_indices] <- ((array(
dispersal_target_n_k$cutoff,
populations
)[modify_pop_indices] -
abundance_on_capacity[modify_pop_indices]) /
array(
dispersal_target_n_k$cutoff - dispersal_target_n_k$threshold,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
} else if (
all(dispersal_target_n_k$threshold > dispersal_target_n_k$cutoff)
) {
# seek company /
dd_multipliers[which(
abundance_on_capacity <= dispersal_target_n_k$cutoff
)] <- 0
modify_pop_indices <- which(
abundance_on_capacity < dispersal_target_n_k$threshold &
dd_multipliers > 0
)
dd_multipliers[modify_pop_indices] <- ((abundance_on_capacity[
modify_pop_indices
] -
array(dispersal_target_n_k$cutoff, populations)[
modify_pop_indices
]) /
array(
dispersal_target_n_k$threshold - dispersal_target_n_k$cutoff,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
}
}
# Select multipliers via target populations for non-zero occupied dispersals
selected_dd_multipliers <- dd_multipliers[dispersal_target_pop_map[,
occupied_indices
][occupied_dispersal_indices]]
# Apply modifying multipliers to dispersals
modify_indices <- which(selected_dd_multipliers < 1)
if (length(modify_indices)) {
modify_dipersal_indices <- occupied_dispersal_indices[modify_indices]
occupied_dispersals[modify_dipersal_indices] <- occupied_dispersals[
modify_dipersal_indices
] *
selected_dd_multipliers[modify_indices]
occupied_dispersal_indices <- which(as.logical(occupied_dispersals)) # > 0
}
} # dispersal depends on target pop N, K or N/K?
# Perform dispersal for each participating stage
for (stage in which(dispersal_stages > 0)) {
# Disperser generation via abundance and corresponding dispersal rates
occupied_abundance <- stage_abundance[stage, occupied_indices]
occupied_abundance_rep <- stage_abundance[
rep(stage, dispersal_compact_rows),
occupied_indices
]
dispersers <- array(0, c(dispersal_compact_rows, occupied_populations))
# Generate dispersers
if (demographic_stochasticity) {
# via binomial distribution
dispersers[occupied_dispersal_indices] <- stats::rbinom(
length(occupied_dispersal_indices),
occupied_abundance_rep[occupied_dispersal_indices],
occupied_dispersals[occupied_dispersal_indices] *
dispersal_stages[stage]
)
} else {
# deterministic
dispersers[occupied_dispersal_indices] <- round(
occupied_abundance_rep[occupied_dispersal_indices] *
occupied_dispersals[occupied_dispersal_indices] *
dispersal_stages[stage]
)
}
# Calculate emigrants
emigrants <- array(.colSums(
dispersers,
m = dispersal_compact_rows,
n = occupied_populations
))
# Check consistency of emigrants (not to exceed abundances)
excessive_indices <- which(emigrants > occupied_abundance)
if (length(excessive_indices) > 0) {
# reduce emigrants to equal abundance via random sampling
for (excessive_index in excessive_indices) {
excessive_rows <- which(as.logical(dispersers[, excessive_index])) # > 0
excessive_dispersers <- dispersers[excessive_rows, excessive_index]
disperser_reduction <- emigrants[excessive_index] -
occupied_abundance[excessive_index]
for (remove_row_index in rep(
excessive_rows,
times = excessive_dispersers
)[sample(sum(excessive_dispersers), size = disperser_reduction)]) {
dispersers[remove_row_index, excessive_index] <- dispersers[
remove_row_index,
excessive_index
] -
1
}
}
emigrants[excessive_indices] <- occupied_abundance[excessive_indices]
}
# Update occupied stage abundance
stage_abundance[stage, occupied_indices] <- stage_abundance[
stage,
occupied_indices
] -
emigrants
# Calculate immigrants via dispersal immigrant map
disperser_indices <- which(as.logical(dispersers)) # > 0
immigrant_array <- array(0, c(dispersal_compact_rows, populations))
immigrant_array[dispersal_immigrant_map[, occupied_indices][
disperser_indices
]] <- dispersers[disperser_indices]
immigrants <- .colSums(
immigrant_array,
m = dispersal_compact_rows,
n = populations
)
# Update population abundances
stage_abundance[stage, ] <- stage_abundance[stage, ] + immigrants
}
# Perform additional dispersal for overcrowded cells (only to cells with room)
if (
(dispersal_depends_on_target_pop_n &&
all(dispersal_target_n$threshold < dispersal_target_n$cutoff)) ||
(dispersal_depends_on_target_pop_n_k &&
all(dispersal_target_n_k$threshold < dispersal_target_n_k$cutoff))
) {
# Flags for dependencies
depends_on_target_pop_n <- (dispersal_depends_on_target_pop_n &&
all(dispersal_target_n$threshold < dispersal_target_n$cutoff))
depends_on_target_pop_n_k <- (dispersal_depends_on_target_pop_n_k &&
all(dispersal_target_n_k$threshold < dispersal_target_n_k$cutoff))
# Get all updated dispersal rates
dispersals <- dispersal_compact_matrix_tm
dispersals[, occupied_indices] <- occupied_dispersals
# Identify overcrowded cells based on stages affected by density
density_abundance <- .colSums(
stage_abundance * as.numeric(density_stages),
m = length(density_stages),
n = populations
)
stage_indices <- which(density_stages > 0 & dispersal_stages > 0)
excessive_indices <- c()
if (depends_on_target_pop_n) {
excessive_indices <- which(
density_abundance > dispersal_target_n$cutoff
)
}
if (depends_on_target_pop_n_k) {
excessive_indices <- unique(c(
excessive_indices,
which(
density_abundance / carrying_capacity > dispersal_target_n_k$cutoff
)
))
}
# Disperse excess from each overcrowded cell (in random order)
for (excessive_index in excessive_indices[sample(length(
excessive_indices
))]) {
# Determine dispersal targets and rates with enough room for excess from overcrowded cell
dispersal_indices <- which(dispersals[, excessive_index] > 0)
target_indices <- dispersal_target_pop_map[, excessive_index][
dispersal_indices
]
if (depends_on_target_pop_n && depends_on_target_pop_n_k) {
indices_with_room <- which((density_abundance <
dispersal_target_n$cutoff &
(density_abundance + 1) / carrying_capacity <=
dispersal_target_n_k$cutoff)[target_indices])
} else if (depends_on_target_pop_n) {
indices_with_room <- which((density_abundance <
dispersal_target_n$cutoff)[target_indices])
} else if (depends_on_target_pop_n_k) {
indices_with_room <- which(((density_abundance + 1) /
carrying_capacity <=
dispersal_target_n_k$cutoff)[target_indices])
}
dispersal_indices <- dispersal_indices[indices_with_room]
target_indices <- target_indices[indices_with_room]
# Disperse excess one at a time sampled via the cell stage abundance distribution
rep_stage_indices <- rep(
stage_indices,
times = stage_abundance[stage_indices, excessive_index]
)
abundance_excess <- 0
if (depends_on_target_pop_n) {
abundance_excess <- density_abundance[excessive_index] -
dispersal_target_n$cutoff
}
if (depends_on_target_pop_n_k) {
abundance_excess <- max(
abundance_excess,
density_abundance[excessive_index] -
floor(
dispersal_target_n_k$cutoff * carrying_capacity[excessive_index]
)
)
}
for (stage_i in rep_stage_indices[sample(
1:length(rep_stage_indices),
size = abundance_excess
)]) {
if (length(target_indices)) {
# Sample target cell
target_i <- target_indices[sample(
length(target_indices),
size = 1,
prob = dispersals[dispersal_indices, excessive_index]
)]
# Perform dispersal
stage_abundance[stage_i, excessive_index] <- stage_abundance[
stage_i,
excessive_index
] -
1 # emigrant
stage_abundance[stage_i, target_i] <- stage_abundance[
stage_i,
target_i
] +
1 # immigrant
# Update target density abundance and potential targets if it becomes full
density_abundance[target_i] <- density_abundance[target_i] + 1
if (
(depends_on_target_pop_n &&
density_abundance[target_i] >= dispersal_target_n$cutoff) ||
(depends_on_target_pop_n_k &&
density_abundance[target_i] / carrying_capacity[target_i] >=
dispersal_target_n_k$cutoff)
) {
# remove from potential targets
full_index <- which(target_indices == target_i)
target_indices <- target_indices[-full_index]
dispersal_indices <- dispersal_indices[-full_index]
}
}
}
}
}
return(stage_abundance)
}
return(dispersal_function)
}
Try the poems package in your browser
Any scripts or data that you put into this service are public.
poems documentation built on Nov. 12, 2025, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions population_dispersal
Documented in population_dispersal
#' Nested functions for population dispersal.
#'
#' Modular functions for the population simulator for performing dispersal of stage
#' abundance at a specified time step via dispersal rates provided.
#'
#' @examples
#' # User-defined dispersal: one-quarter of dispersing stages move one population over
#' simulator <- SimulatorReference$new()
#' example_function <- function(params) {
#' params$simulator$attached$params <- params # attach to reference object
#' emigrants <- round(params$stage_abundance * params$dispersal_stages * 0.25)
#' return(params$stage_abundance - emigrants + emigrants[, c(7, 1:6)])
#' }
#' dispersal_function <- population_dispersal(
#' replicates = 4,
#' time_steps = 10,
#' years_per_step = 1,
#' populations = 7,
#' demographic_stochasticity = TRUE,
#' density_stages = c(0, 1, 1),
#' dispersal = example_function,
#' dispersal_stages = c(0, 1, 0.5),
#' dispersal_source_n_k = list(cutoff = -0.5, threshold = 1.5),
#' dispersal_target_k = 5,
#' dispersal_target_n = list(threshold = 10, cutoff = 15),
#' simulator = simulator
#' )
#' carrying_capacity <- rep(10, 7)
#' stage_abundance <- matrix(
#' c(
#' 7, 13, 0, 26, 0, 39, 47,
#' 2, 0, 6, 8, 0, 12, 13,
#' 0, 3, 4, 6, 0, 9, 10
#' ),
#' nrow = 3,
#' ncol = 7,
#' byrow = TRUE
#' )
#' occupied_indices <- (1:7)[-5]
#' dispersal_function(
#' r = 2, tm = 6, carrying_capacity, stage_abundance,
#' occupied_indices
#' )
#'
#' @param replicates Number of replicate simulation runs.
#' @param time_steps Number of simulation time steps.
#' @param years_per_step Number of years per time step.
#' @param populations Number of populations.
#' @param demographic_stochasticity Boolean for optionally choosing demographic stochasticity for the transformation.
#' @param density_stages Array of booleans or numeric (0,1) for each stage to indicate which stages are affected by density.
#' @param dispersal Either a matrix of dispersal rates between populations (source columns to target rows) or a list of data frames of non-zero dispersal rates and indices for constructing a compact dispersal matrix, and optional changing rates over time (as per class \code{\link{DispersalGenerator}} \emph{dispersal_data} attribute). Alternatively a user-defined function (optionally nested in a list with additional attributes) may be used: \code{function(params)}, where \emph{params} is a list passed to the function containing:
#' \describe{
#' \item{\code{replicates}}{Number of replicate simulation runs.}
#' \item{\code{time_steps}}{Number of simulation time steps.}
#' \item{\code{years_per_step}}{Number of years per time step.}
#' \item{\code{populations}}{Number of populations.}
#' \item{\code{stages}}{Number of life cycle stages.}
#' \item{\code{demographic_stochasticity}}{Boolean for optionally choosing
#' demographic stochasticity for the transformation.}
#' \item{\code{density_stages}}{Array of booleans or numeric (0,1) for each
#' stage to indicate which stages are affected by density.}
#' \item{\code{dispersal_stages}}{Array of relative dispersal (0-1) for
#' each stage to indicate the degree to which each stage participates in
#' dispersal. This factor modifies dispersal proportion, not dispersal rate.}
#' \item{\code{dispersal_source_n_k}}{Dispersal proportion (p) density dependence via source population abundance divided by carrying capacity (n/k), where p is reduced via a linear slope (defined by two list items) from n/k <= \emph{cutoff} (p = 0) to n/k >= \emph{threshold}.}
#' \item{\code{dispersal_target_k}}{Dispersal rate (r) density dependence via target population carrying capacity (k), where r is reduced via a linear slope (through the origin) when k <= \emph{threshold}.}
#' \item{\code{dispersal_target_n}}{Dispersal rate (r) density dependence via target population abundance (n), where r is reduced via a linear slope (defined by two list items) from n >= \emph{threshold} to n <= \emph{cutoff} (r = 0) or vice versa.}
#' \item{\code{dispersal_target_n_k}}{Dispersal rate (r) density dependence via target population abundance divided by carrying capacity (n/k), where r is reduced via a linear slope (defined by two list items) from n/k >= \emph{threshold} to n/k <= \emph{cutoff} (r = 0) or vice versa.}
#' \item{\code{r}}{Simulation replicate.}
#' \item{\code{tm}}{Simulation time step.}
#' \item{\code{carrying_capacity}}{Array of carrying capacity values for each population at time step.}
#' \item{\code{stage_abundance}}{Matrix of abundance for each stage (rows) and population (columns) at time step.}
#' \item{\code{occupied_indices}}{Array of indices for populations occupied at time step.}
#' \item{\code{simulator}}{\code{\link{SimulatorReference}} object with dynamically accessible \emph{attached} and \emph{results} lists.}
#' \item{\code{additional attributes}}{Additional attributes when the transformation is optionally nested in a list.}
#' }
#' returns the post-dispersal abundance matrix
#' @param dispersal_stages Array of relative dispersal (0-1) for each stage to
#' indicate the degree to which each stage participates in dispersal (default
#' is 1 for all stages). This factor modifies dispersal proportion, not
#' dispersal rate.
#' @param dispersal_source_n_k Dispersal proportion (p) density dependence via
#' source population abundance divided by carrying capacity (n/k), where p is
#' reduced via a linear slope (defined by two list items) from n/k <= \emph{cutoff} (p = 0) to n/k >= \emph{threshold} or vice versa.
#' @param dispersal_target_k Dispersal rate (r) density dependence via target population carrying capacity (k), where r is reduced via a linear slope (through the origin) when k <= \emph{threshold}.
#' @param dispersal_target_n Dispersal rate (r) density dependence via target population abundance (n), where r is reduced via a linear slope (defined by two list items) from n >= \emph{threshold} to n <= \emph{cutoff} (r = 0) or visa-versa.
#' @param dispersal_target_n_k Dispersal rate (r) density dependence via target population abundance divided by carrying capacity (n/k), where r is reduced via a linear slope (defined by two list items) from n/k >= \emph{threshold} to n/k <= \emph{cutoff} (r = 0) or vice versa.
#' @param simulator \code{\link{SimulatorReference}} object with dynamically accessible \emph{attached} and \emph{results} lists.
#' @return Dispersal function: \code{function(r, tm, carrying_capacity, stage_abundance, occupied_indices)}, where:
#' \describe{
#' \item{\code{r}}{Simulation replicate.}
#' \item{\code{tm}}{Simulation time step.}
#' \item{\code{carrying_capacity}}{Array of carrying capacity values for each population at time step.}
#' \item{\code{stage_abundance}}{Matrix of abundance for each stage (rows) and population (columns) at time step.}
#' \item{\code{occupied_indices}}{Array of indices for populations occupied at time step.}
#' \item{\code{returns}}{New stage abundance matrix with dispersal applied.}
#' }
#' @export population_dispersal
population_dispersal <- function(
replicates,
time_steps,
years_per_step,
populations,
demographic_stochasticity,
density_stages,
dispersal,
dispersal_stages,
dispersal_source_n_k = NULL,
dispersal_target_k = NULL,
dispersal_target_n = NULL,
dispersal_target_n_k = NULL,
simulator
) {
if (is.null(dispersal)) {
# no dispersal
return(NULL)
}
# User-defined function?
if (
is.function(dispersal) ||
(is.list(dispersal) &&
length(which(unlist(lapply(dispersal, is.function)))))
) {
# Derive stages
stages <- length(dispersal_stages)
# Unpack transformation function and additional attributes from a list
additional_attributes <- list()
if (is.list(dispersal)) {
function_index <- which(unlist(lapply(dispersal, is.function)))
additional_attributes <- dispersal[-function_index]
dispersal <- dispersal[[function_index]]
}
if (is.function(dispersal)) {
# user-defined function
# List of parameters to pass to the user-defined function
params <- c(
list(
replicates = replicates,
time_steps = time_steps,
years_per_step = years_per_step,
populations = populations,
stages = stages,
demographic_stochasticity = demographic_stochasticity,
density_stages = density_stages,
dispersal_stages = dispersal_stages,
dispersal_source_n_k = dispersal_source_n_k,
dispersal_target_k = dispersal_target_k,
dispersal_target_n = dispersal_target_n,
dispersal_target_n_k = dispersal_target_n_k,
simulator = simulator
),
additional_attributes
)
## Create a nested function for applying user-defined dispersal of stage abundance ##
user_defined_function <- function(
r,
tm,
carrying_capacity,
stage_abundance,
occupied_indices
) {
# Add attributes to be made available to the user-defined function
params$r <- r
params$tm <- tm
params$carrying_capacity <- carrying_capacity
params$stage_abundance <- stage_abundance
params$occupied_indices <- occupied_indices
# Run user-defined dispersal function
tryCatch(
{
stage_abundance <- dispersal(params)
},
error = function(e) {
stop(
paste(
"Error produced within user-defined dispersal function:",
as.character(e)
),
call. = FALSE
)
}
)
# Warn if any negative or non-finite
if (any(!is.finite(stage_abundance))) {
warning(
"Non-finite abundances returned by user-defined dispersal function",
call. = FALSE
)
}
if (any(stage_abundance[which(is.finite(stage_abundance))] < 0)) {
warning(
"Negative abundances returned by user-defined dispersal function",
call. = FALSE
)
}
return(stage_abundance)
}
return(user_defined_function)
}
}
# Assumed to be matrix or generated dispersal data
# Initialize reusable dispersal attributes
if (is.matrix(dispersal)) {
# create compact matrix data
# Calculate the indices of non-zero dispersals
dispersal_data <- data.frame(which(dispersal > 0, arr.ind = TRUE))
dispersal_data <- dispersal_data[
order(dispersal_data[, 2], dispersal_data[, 1]),
]
names(dispersal_data) <- c("target_pop", "source_pop")
# Calculate indices for constructing compacted dispersal matrices for emigrants and immigrants
dispersal_rows <- tabulate(dispersal_data$source_pop, nbins = populations)
dispersal_cols <- tabulate(dispersal_data$target_pop, nbins = populations)
dispersal_compact_rows <- max(dispersal_rows, dispersal_cols)
compact_emigrant_matrix <- array(
1:dispersal_compact_rows,
c(dispersal_compact_rows, populations)
)
compact_immigrant_matrix <- compact_emigrant_matrix *
(compact_emigrant_matrix <=
matrix(
dispersal_cols,
nrow = dispersal_compact_rows,
ncol = populations,
byrow = TRUE
))
compact_emigrant_matrix <- compact_emigrant_matrix *
(compact_emigrant_matrix <=
matrix(
dispersal_rows,
nrow = dispersal_compact_rows,
ncol = populations,
byrow = TRUE
))
# Map the row of each compact matrix to the original target (for emigrants) or source (for immigrants) populations
dispersal_data$emigrant_row <- which(
compact_emigrant_matrix > 0,
arr.ind = TRUE,
useNames = FALSE
)[, 1]
dispersal_data$immigrant_row <- which(
compact_immigrant_matrix > 0,
arr.ind = TRUE,
useNames = FALSE
)[, 1]
target_sorted_indices <- dispersal_data[
order(dispersal_data$target_pop, dispersal_data$source_pop),
c("target_pop", "source_pop")
]
dispersal_data$immigrant_row <- dispersal_data$immigrant_row[order(
target_sorted_indices$source_pop,
target_sorted_indices$target_pop
)]
# Add dispersal rates
dispersal_data$dispersal_rate <- dispersal[as.matrix(dispersal_data[c(
"target_pop",
"source_pop"
)])]
dispersals_change_over_time <- FALSE
# Release variables from memory
dispersal_rows <- NULL
dispersal_cols <- NULL
compact_emigrant_matrix <- NULL
compact_immigrant_matrix <- NULL
target_sorted_indices <- NULL
} else if (
is.list(dispersal) && is.data.frame(dispersal[[1]]) && nrow(dispersal[[1]])
) {
# compact matrix data in list (as per DispersalModel class)
# Unpack dispersal data and determine compact matrix dimensions
dispersal_data <- dispersal[[1]]
dispersal_compact_rows <- max(dispersal_data[, c(
"emigrant_row",
"immigrant_row"
)])
# Are dispersals changing over time?
dispersals_change_over_time <- (length(dispersal) > 1)
if (dispersals_change_over_time) {
dispersal_data_changes <- dispersal
dispersal_data_changes[[1]] <- dispersal_data_changes[[1]][NULL, ]
}
} else {
return(NULL)
}
# Release dispersal from memory
dispersal <- NULL
# Create a compact matrix of dispersal rates
dispersal_compact_matrix <- array(0, c(dispersal_compact_rows, populations))
dispersal_compact_matrix[as.matrix(dispersal_data[, c(
"emigrant_row",
"source_pop"
)])] <- dispersal_data$dispersal_rate
# Does dispersal depend on source population abundance N divided by carrying capacity K?
dispersal_depends_on_source_pop_n_k <- (is.list(dispersal_source_n_k) &&
(is.numeric(dispersal_source_n_k$cutoff) ||
is.numeric(dispersal_source_n_k$threshold)))
# Does dispersal depend on target population carrying capacity K, abundance N, or N/K?
dispersal_depends_on_target_pop_k <- is.numeric(dispersal_target_k)
dispersal_depends_on_target_pop_n <- (is.list(dispersal_target_n) &&
(is.numeric(dispersal_target_n$threshold) ||
is.numeric(dispersal_target_n$cutoff)))
dispersal_depends_on_target_pop_n_k <- (is.list(dispersal_target_n_k) &&
(is.numeric(dispersal_target_n_k$threshold) ||
is.numeric(dispersal_target_n_k$cutoff)))
# Setup density dependence dispersal parameters
if (dispersal_depends_on_source_pop_n_k) {
# Convert NULL to zero in source N/K cutoff or one in threshold
if (dispersal_depends_on_source_pop_n_k) {
if (is.null(dispersal_source_n_k$cutoff)) dispersal_source_n_k$cutoff <- 0
if (is.null(dispersal_source_n_k$threshold))
dispersal_source_n_k$threshold <- 1
}
# Check threshold > cutoff
if (dispersal_source_n_k$threshold <= dispersal_source_n_k$cutoff) {
dispersal_depends_on_source_pop_n_k <- FALSE
warning(
"Dispersal density dependence for source N/K threshold must be greater than cutoff => not used",
call. = FALSE
)
}
}
if (
dispersal_depends_on_target_pop_k ||
dispersal_depends_on_target_pop_n ||
dispersal_depends_on_target_pop_n_k
) {
if (dispersal_depends_on_target_pop_n) {
# Convert NULL to zero in target N threshold or cutoff
if (is.null(dispersal_target_n$threshold))
dispersal_target_n$threshold <- 0
if (is.null(dispersal_target_n$cutoff)) dispersal_target_n$cutoff <- 0
}
if (dispersal_depends_on_target_pop_n_k) {
# Convert NULL to zero in target N/K threshold or cutoff
if (is.null(dispersal_target_n_k$threshold))
dispersal_target_n_k$threshold <- 0
if (is.null(dispersal_target_n_k$cutoff)) dispersal_target_n_k$cutoff <- 0
}
# Create a map of compact array indices for mapping dispersers (emigrants) to target populations
dispersal_target_pop_map <- array(0, c(dispersal_compact_rows, populations))
dispersal_target_pop_map[as.matrix(dispersal_data[, c(
"emigrant_row",
"source_pop"
)])] <- dispersal_data$target_pop
}
# Create a map of compact array indices for mapping dispersers (emigrants) to immigrants
dispersal_compact_indices <- array(
1:(dispersal_compact_rows * populations),
c(dispersal_compact_rows, populations)
)
dispersal_immigrant_map <- array(0, c(dispersal_compact_rows, populations))
dispersal_immigrant_map[as.matrix(dispersal_data[, c(
"emigrant_row",
"source_pop"
)])] <- dispersal_compact_indices[as.matrix(dispersal_data[, c(
"immigrant_row",
"target_pop"
)])]
# Release variables from memory
dispersal_data <- NULL
dispersal_compact_indices <- NULL
## Create a nested function for performing dispersal ##
dispersal_function <- function(
r,
tm,
carrying_capacity,
stage_abundance,
occupied_indices
) {
# Calculate occupied population number
occupied_populations <- length(occupied_indices)
# Apply any spatio-temporal dispersal changes
dispersal_compact_matrix_tm <- simulator$attached$dispersal_compact_matrix_tm
if (tm == 1 || !dispersals_change_over_time) {
dispersal_compact_matrix_tm <- dispersal_compact_matrix
} else if (
dispersals_change_over_time && nrow(dispersal_data_changes[[tm]])
) {
# and tm > 1
dispersal_compact_matrix_tm[as.matrix(dispersal_data_changes[[tm]][, c(
"emigrant_row",
"source_pop"
)])] <- dispersal_data_changes[[tm]]$dispersal_rate
}
simulator$attached$dispersal_compact_matrix_tm <- dispersal_compact_matrix_tm
# Select dispersals for occupied populations
occupied_dispersals <- dispersal_compact_matrix_tm[, occupied_indices]
# Calculate density abundance
if (
dispersal_depends_on_source_pop_n_k ||
dispersal_depends_on_target_pop_n ||
dispersal_depends_on_target_pop_n_k
) {
density_abundance <- .colSums(
stage_abundance * as.numeric(density_stages),
m = length(density_stages),
n = populations
)
}
# Modify dispersal rates when dispersal depends on source population N/K
if (dispersal_depends_on_source_pop_n_k) {
# Density dependent multipliers
dd_multipliers <- array(1, populations)
# Calculate the source N/K multipliers
abundance_on_capacity <- density_abundance / carrying_capacity
dd_multipliers[which(
abundance_on_capacity <= dispersal_source_n_k$cutoff
)] <- 0
modify_pop_indices <- which(
carrying_capacity > 0 &
dd_multipliers > 0 &
abundance_on_capacity < dispersal_source_n_k$threshold
)
dd_multipliers[modify_pop_indices] <- ((abundance_on_capacity[
modify_pop_indices
] -
array(dispersal_source_n_k$cutoff, populations)[modify_pop_indices]) /
array(
dispersal_source_n_k$threshold - dispersal_source_n_k$cutoff,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
# Apply modifying multipliers to dispersals
occupied_dispersals <- (occupied_dispersals *
matrix(
dd_multipliers[occupied_indices],
nrow = dispersal_compact_rows,
ncol = occupied_populations,
byrow = TRUE
))
} # dispersal depends on source pop N/K?
# Select occupied dispersal non-zero indices
occupied_dispersal_indices <- which(as.logical(occupied_dispersals)) # > 0
# Modify dispersal rates when dispersal depends on target population K, N, or N/K
if (
dispersal_depends_on_target_pop_k ||
dispersal_depends_on_target_pop_n ||
dispersal_depends_on_target_pop_n_k
) {
# Density dependent multipliers
dd_multipliers <- array(1, populations)
# Calculate the (below-threshold) target K multipliers
if (dispersal_depends_on_target_pop_k) {
modify_pop_indices <- which(carrying_capacity < dispersal_target_k)
dd_multipliers[modify_pop_indices] <- (carrying_capacity[
modify_pop_indices
] /
array(dispersal_target_k, populations)[modify_pop_indices])
}
# Calculate the target N multipliers
if (dispersal_depends_on_target_pop_n) {
if (all(dispersal_target_n$threshold < dispersal_target_n$cutoff)) {
# overcrowded cell avoidance \
dd_multipliers[which(
density_abundance >= dispersal_target_n$cutoff
)] <- 0
modify_pop_indices <- which(
density_abundance > dispersal_target_n$threshold &
dd_multipliers > 0
)
dd_multipliers[modify_pop_indices] <- ((array(
dispersal_target_n$cutoff,
populations
)[modify_pop_indices] -
density_abundance[modify_pop_indices]) /
array(
dispersal_target_n$cutoff - dispersal_target_n$threshold,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
} else if (
all(dispersal_target_n$threshold > dispersal_target_n$cutoff)
) {
# seek company /
dd_multipliers[which(
density_abundance <= dispersal_target_n$cutoff
)] <- 0
modify_pop_indices <- which(
density_abundance < dispersal_target_n$threshold &
dd_multipliers > 0
)
dd_multipliers[modify_pop_indices] <- ((density_abundance[
modify_pop_indices
] -
array(dispersal_target_n$cutoff, populations)[modify_pop_indices]) /
array(
dispersal_target_n$threshold - dispersal_target_n$cutoff,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
}
}
# Calculate the target N/K multipliers
if (dispersal_depends_on_target_pop_n_k) {
dd_multipliers[which(carrying_capacity <= 0)] <- 0
abundance_on_capacity <- density_abundance / carrying_capacity
if (all(dispersal_target_n_k$threshold < dispersal_target_n_k$cutoff)) {
# overcrowded cell avoidance \
dd_multipliers[which(
abundance_on_capacity >= dispersal_target_n_k$cutoff
)] <- 0
modify_pop_indices <- which(
abundance_on_capacity > dispersal_target_n_k$threshold &
dd_multipliers > 0
)
dd_multipliers[modify_pop_indices] <- ((array(
dispersal_target_n_k$cutoff,
populations
)[modify_pop_indices] -
abundance_on_capacity[modify_pop_indices]) /
array(
dispersal_target_n_k$cutoff - dispersal_target_n_k$threshold,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
} else if (
all(dispersal_target_n_k$threshold > dispersal_target_n_k$cutoff)
) {
# seek company /
dd_multipliers[which(
abundance_on_capacity <= dispersal_target_n_k$cutoff
)] <- 0
modify_pop_indices <- which(
abundance_on_capacity < dispersal_target_n_k$threshold &
dd_multipliers > 0
)
dd_multipliers[modify_pop_indices] <- ((abundance_on_capacity[
modify_pop_indices
] -
array(dispersal_target_n_k$cutoff, populations)[
modify_pop_indices
]) /
array(
dispersal_target_n_k$threshold - dispersal_target_n_k$cutoff,
populations
)[modify_pop_indices] *
dd_multipliers[modify_pop_indices])
}
}
# Select multipliers via target populations for non-zero occupied dispersals
selected_dd_multipliers <- dd_multipliers[dispersal_target_pop_map[,
occupied_indices
][occupied_dispersal_indices]]
# Apply modifying multipliers to dispersals
modify_indices <- which(selected_dd_multipliers < 1)
if (length(modify_indices)) {
modify_dipersal_indices <- occupied_dispersal_indices[modify_indices]
occupied_dispersals[modify_dipersal_indices] <- occupied_dispersals[
modify_dipersal_indices
] *
selected_dd_multipliers[modify_indices]
occupied_dispersal_indices <- which(as.logical(occupied_dispersals)) # > 0
}
} # dispersal depends on target pop N, K or N/K?
# Perform dispersal for each participating stage
for (stage in which(dispersal_stages > 0)) {
# Disperser generation via abundance and corresponding dispersal rates
occupied_abundance <- stage_abundance[stage, occupied_indices]
occupied_abundance_rep <- stage_abundance[
rep(stage, dispersal_compact_rows),
occupied_indices
]
dispersers <- array(0, c(dispersal_compact_rows, occupied_populations))
# Generate dispersers
if (demographic_stochasticity) {
# via binomial distribution
dispersers[occupied_dispersal_indices] <- stats::rbinom(
length(occupied_dispersal_indices),
occupied_abundance_rep[occupied_dispersal_indices],
occupied_dispersals[occupied_dispersal_indices] *
dispersal_stages[stage]
)
} else {
# deterministic
dispersers[occupied_dispersal_indices] <- round(
occupied_abundance_rep[occupied_dispersal_indices] *
occupied_dispersals[occupied_dispersal_indices] *
dispersal_stages[stage]
)
}
# Calculate emigrants
emigrants <- array(.colSums(
dispersers,
m = dispersal_compact_rows,
n = occupied_populations
))
# Check consistency of emigrants (not to exceed abundances)
excessive_indices <- which(emigrants > occupied_abundance)
if (length(excessive_indices) > 0) {
# reduce emigrants to equal abundance via random sampling
for (excessive_index in excessive_indices) {
excessive_rows <- which(as.logical(dispersers[, excessive_index])) # > 0
excessive_dispersers <- dispersers[excessive_rows, excessive_index]
disperser_reduction <- emigrants[excessive_index] -
occupied_abundance[excessive_index]
for (remove_row_index in rep(
excessive_rows,
times = excessive_dispersers
)[sample(sum(excessive_dispersers), size = disperser_reduction)]) {
dispersers[remove_row_index, excessive_index] <- dispersers[
remove_row_index,
excessive_index
] -
1
}
}
emigrants[excessive_indices] <- occupied_abundance[excessive_indices]
}
# Update occupied stage abundance
stage_abundance[stage, occupied_indices] <- stage_abundance[
stage,
occupied_indices
] -
emigrants
# Calculate immigrants via dispersal immigrant map
disperser_indices <- which(as.logical(dispersers)) # > 0
immigrant_array <- array(0, c(dispersal_compact_rows, populations))
immigrant_array[dispersal_immigrant_map[, occupied_indices][
disperser_indices
]] <- dispersers[disperser_indices]
immigrants <- .colSums(
immigrant_array,
m = dispersal_compact_rows,
n = populations
)
# Update population abundances
stage_abundance[stage, ] <- stage_abundance[stage, ] + immigrants
}
# Perform additional dispersal for overcrowded cells (only to cells with room)
if (
(dispersal_depends_on_target_pop_n &&
all(dispersal_target_n$threshold < dispersal_target_n$cutoff)) ||
(dispersal_depends_on_target_pop_n_k &&
all(dispersal_target_n_k$threshold < dispersal_target_n_k$cutoff))
) {
# Flags for dependencies
depends_on_target_pop_n <- (dispersal_depends_on_target_pop_n &&
all(dispersal_target_n$threshold < dispersal_target_n$cutoff))
depends_on_target_pop_n_k <- (dispersal_depends_on_target_pop_n_k &&
all(dispersal_target_n_k$threshold < dispersal_target_n_k$cutoff))
# Get all updated dispersal rates
dispersals <- dispersal_compact_matrix_tm
dispersals[, occupied_indices] <- occupied_dispersals
# Identify overcrowded cells based on stages affected by density
density_abundance <- .colSums(
stage_abundance * as.numeric(density_stages),
m = length(density_stages),
n = populations
)
stage_indices <- which(density_stages > 0 & dispersal_stages > 0)
excessive_indices <- c()
if (depends_on_target_pop_n) {
excessive_indices <- which(
density_abundance > dispersal_target_n$cutoff
)
}
if (depends_on_target_pop_n_k) {
excessive_indices <- unique(c(
excessive_indices,
which(
density_abundance / carrying_capacity > dispersal_target_n_k$cutoff
)
))
}
# Disperse excess from each overcrowded cell (in random order)
for (excessive_index in excessive_indices[sample(length(
excessive_indices
))]) {
# Determine dispersal targets and rates with enough room for excess from overcrowded cell
dispersal_indices <- which(dispersals[, excessive_index] > 0)
target_indices <- dispersal_target_pop_map[, excessive_index][
dispersal_indices
]
if (depends_on_target_pop_n && depends_on_target_pop_n_k) {
indices_with_room <- which((density_abundance <
dispersal_target_n$cutoff &
(density_abundance + 1) / carrying_capacity <=
dispersal_target_n_k$cutoff)[target_indices])
} else if (depends_on_target_pop_n) {
indices_with_room <- which((density_abundance <
dispersal_target_n$cutoff)[target_indices])
} else if (depends_on_target_pop_n_k) {
indices_with_room <- which(((density_abundance + 1) /
carrying_capacity <=
dispersal_target_n_k$cutoff)[target_indices])
}
dispersal_indices <- dispersal_indices[indices_with_room]
target_indices <- target_indices[indices_with_room]
# Disperse excess one at a time sampled via the cell stage abundance distribution
rep_stage_indices <- rep(
stage_indices,
times = stage_abundance[stage_indices, excessive_index]
)
abundance_excess <- 0
if (depends_on_target_pop_n) {
abundance_excess <- density_abundance[excessive_index] -
dispersal_target_n$cutoff
}
if (depends_on_target_pop_n_k) {
abundance_excess <- max(
abundance_excess,
density_abundance[excessive_index] -
floor(
dispersal_target_n_k$cutoff * carrying_capacity[excessive_index]
)
)
}
for (stage_i in rep_stage_indices[sample(
1:length(rep_stage_indices),
size = abundance_excess
)]) {
if (length(target_indices)) {
# Sample target cell
target_i <- target_indices[sample(
length(target_indices),
size = 1,
prob = dispersals[dispersal_indices, excessive_index]
)]
# Perform dispersal
stage_abundance[stage_i, excessive_index] <- stage_abundance[
stage_i,
excessive_index
] -
1 # emigrant
stage_abundance[stage_i, target_i] <- stage_abundance[
stage_i,
target_i
] +
1 # immigrant
# Update target density abundance and potential targets if it becomes full
density_abundance[target_i] <- density_abundance[target_i] + 1
if (
(depends_on_target_pop_n &&
density_abundance[target_i] >= dispersal_target_n$cutoff) ||
(depends_on_target_pop_n_k &&
density_abundance[target_i] / carrying_capacity[target_i] >=
dispersal_target_n_k$cutoff)
) {
# remove from potential targets
full_index <- which(target_indices == target_i)
target_indices <- target_indices[-full_index]
dispersal_indices <- dispersal_indices[-full_index]
}
}
}
}
}
return(stage_abundance)
}
return(dispersal_function)
}
Try the poems package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.