Nothing
R/paleopop_simulator.R
In paleopop: Pattern-Oriented Modeling Framework for Coupled Niche-Population Paleo-Climatic Models
Defines functions
paleopop_simulator
Documented in
paleopop_simulator
#' Runs a customized population model simulation.
#'
#' Simulates a population model customized for paleontological time-scales, optimized
#' for single-generation transitions and large populations, across multiple generations
#' and returns simulation results. Each generational time-step includes:
#' \enumerate{
#' \item Density dependence calculations
#' \item Environmental stochasticity calculations
#' \item Generational transition calculations
#' \item Harvest calculations
#' \item Dispersal calculations
#' \item Results collection
#' }
#'
#' @param inputs Nested list/object with named elements:
#' \describe{
#' \item{\code{random_seed}}{Number to seed the random number generation for stochasticity.}
#' \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{initial_abundance}}{Array of initial abundances for each population.}
#' \item{\code{transition_rate}}{Rate of transition (or fecundity) between generations.}
#' \item{\code{standard_deviation}}{Standard deviation applied to transition rates.}
#' \item{\code{compact_decomposition}}{List containing a compact transposed (Cholesky) decomposition \emph{matrix} (t_decomposition_compact_matrix) and a corresponding \emph{map} of population indices (t_decomposition_compact_map), as per \code{\link{SpatialCorrelation}} class attributes.}
#' \item{\code{carrying_capacity}}{Matrix of carrying capacities (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{density_dependence}}{Density dependence type ("competition", "logistic", or "ceiling").}
#' \item{\code{growth_rate_max}}{Maximum growth rate (for "competition" or "logistic" density dependence).}
#' \item{\code{harvest}}{Boolean for utilizing harvesting.}
#' \item{\code{harvest_max}}{Proportion harvested per year (note: annual time scale - not generational).}
#' \item{\code{harvest_g}}{The \emph{G} parameter in the harvest function.}
#' \item{\code{harvest_z}}{The \emph{Z} parameter in the harvest function.}
#' \item{\code{harvest_max_n}}{Maximum density per grid cell.}
#' \item{\code{human_density}}{Matrix of human density (fraction) (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{dispersal_data}}{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.}
#' \item{\code{dispersal_target_k}}{Target population carrying capacity threshold for density dependent dispersal.}
#' \item{\code{abundance_threshold}}{Abundance threshold (that needs to be exceeded) for each population to persist.}
#' \item{\code{occupancy_threshold}}{Threshold for the number of populations occupied (that needs to be exceeded) for all populations to persist.}
#' \item{\code{results_selection}}{List of results selection from: "abundance", "ema", "extirpation", "harvested", "occupancy", "human_density".}
#' }
#' @return Simulation results as a nested list (as selected):
#' \describe{
#' \item{\code{abundance}}{Matrix of simulation abundances (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{ema}}{Matrix of expected minimum abundances (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{extirpation}}{Array of extirpation times for each population.}
#' \item{\code{harvested}}{Matrix of estimated individuals harvested (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{occupancy}}{Array of number of populations occupied at each time-step.}
#' \item{\code{human_density}}{Matrix of human densities, (\emph{populations} rows by \emph{time_steps} columns).}
#' }
#'
#' @examples
#' library(raster)
#' library(poems)
#' # Ring Island example region
#' coordinates <- data.frame(x = rep(seq(-178.02, -178.06, -0.01), 5),
#' y = rep(seq(19.02, 19.06, 0.01), each = 5))
#' template_raster <- Region$new(coordinates = coordinates)$region_raster # full extent
#' sealevel_raster <- template_raster
#' template_raster[][c(7:9, 12:14, 17:19)] <- NA # make Ring Island
#' sealevel_raster[][c(7:9, 12:14, 17:18)] <- NA
#' raster_stack <- raster::stack(x = append(replicate(9, template_raster), sealevel_raster))
#' region <- PaleoRegion$new(template_raster = raster_stack)
#'
#' # Model template
#' model_template <- PaleoPopModel$new(
#' region = region,
#' time_steps = 10,
#' years_per_step = 12, # years per generational time-step
#' standard_deviation = 0.1,
#' growth_rate_max = 0.6,
#' harvest = FALSE,
#' populations = region$region_cells,
#' initial_abundance = seq(9000, 0, -1000),
#' transition_rate = 1.0,
#' carrying_capacity = rep(1000, 17),
#' dispersal = (!diag(nrow = 17, ncol = 17))*0.05,
#' density_dependence = "logistic",
#' dispersal_target_k = 10,
#' occupancy_threshold = 1,
#' abundance_threshold = 10,
#' results_selection = c("abundance")
#' )
#'
#' # Simulations
#' results <- paleopop_simulator(model_template) # input as PaleoPopModel object
#' inputs <- model_template$get_attributes()
#' paleopop_simulator(inputs) # input as list of attributes
#'
#' @export paleopop_simulator
paleopop_simulator <- function(inputs) {
## Unpack inputs and calculate re-usable variables ##
# Stop if minimal inputs are not present
if (is.null(inputs$time_steps) || is.null(inputs$populations) || is.null(inputs$initial_abundance) ||
is.null(inputs$transition_rate) || is.null(inputs$carrying_capacity)) {
incomplete_inputs <- if (is.null(inputs$time_steps)) "time_steps"
incomplete_inputs <- c(incomplete_inputs, if (is.null(inputs$populations)) "populations")
incomplete_inputs <- c(incomplete_inputs, if (is.null(inputs$initial_abundance)) "initial_abundance")
incomplete_inputs <- c(incomplete_inputs, if (is.null(inputs$transition_rate)) "transition_rate")
incomplete_inputs <- c(incomplete_inputs, if (is.null(inputs$carrying_capacity)) "carrying_capacity")
stop(paste("Minimal inputs required to run simulation should include:",
paste(incomplete_inputs, collapse = ", ")), call. = FALSE)
}
# General model settings
time_steps <- inputs$time_steps
years_per_step <- ifelse(is.null(inputs$years_per_step), 1, inputs$years_per_step)
random_seed <- inputs$random_seed
if (!is.null(random_seed)) {
set.seed(random_seed)
}
# Transition rate and standard deviation
transition_rate <- inputs$transition_rate
standard_deviation <- ifelse(is.null(inputs$standard_deviation), 0, inputs$standard_deviation)
environmental_stochasticity <- (standard_deviation*transition_rate > 0)
# Population settings
populations <- inputs$populations
population_abundances <- array(inputs$initial_abundance, c(1, populations)) # single column matrix for faster dispersal calculations
density_dependence <- inputs$density_dependence
if (is.null(density_dependence)) {
density_dependence <- "none"
}
growth_rate_max <- ifelse(is.null(inputs$growth_rate_max), log(transition_rate), inputs$growth_rate_max)
abundance_threshold <- inputs$abundance_threshold
occupancy_threshold <- inputs$occupancy_threshold
dispersal_target_k <- inputs$dispersal_target_k
carrying_capacities <- matrix(inputs$carrying_capacity, nrow = populations)
if (any(is.na(carrying_capacities))) {
carrying_capacities[which(is.na(carrying_capacities))] <- 0
}
# Dispersal and correlation
# Create re-usable dispersal structures
dispersal_present <- (!is.null(inputs$dispersal_data) && nrow(inputs$dispersal_data[[1]]))
if (dispersal_present) {
# Unpack dispersal data and determine compact matrix dimensions
dispersal_data <- inputs$dispersal_data[[1]]
dispersal_compact_rows <- max(dispersal_data[, c("emigrant_row", "immigrant_row")])
# Create compact array of zeros for quick initialization
dispersal_zero_array <- array(0, c(dispersal_compact_rows, populations))
# Create a compact matrix of dispersal rates
dispersal_compact_matrix <- dispersal_zero_array
dispersal_compact_matrix[as.matrix(dispersal_data[, c("emigrant_row", "source_pop")])] <- dispersal_data$dispersal_rate
# Are dispersals changing over time?
dispersals_change_over_time <- (length(inputs$dispersal_data) > 1)
if (dispersals_change_over_time) {
dispersal_data_changes <- inputs$dispersal_data
dispersal_data_changes[[1]] <- dispersal_data_changes[[1]][NULL,]
}
# Does dispersal dependent on target population carrying capacity K?
dispersal_depends_on_target_pop_k <- (!is.null(dispersal_target_k) && dispersal_target_k > 0)
if (dispersal_depends_on_target_pop_k) {
# Create a map of compact array indices for mapping dispersers (emigrants) to target populations
dispersal_target_pop_map <- dispersal_zero_array
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 <- dispersal_zero_array
dispersal_immigrant_map[as.matrix(dispersal_data[, c("emigrant_row", "source_pop")])] <- dispersal_compact_indices[as.matrix(dispersal_data[, c("immigrant_row", "target_pop")])]
# Create indices for replicating abundances
dispersal_abundance_rep_indices <- rep(1, dispersal_compact_rows)
# Release variables from memory
dispersal_data <- NULL; dispersal_compact_indices <- NULL
} # dispersal present?
# Environmental correlation?
if (environmental_stochasticity) {
# Flag for environmental correlation
use_env_correlation <- FALSE
# Unpack compact decompostion
if (!is.null(inputs$compact_decomposition)) {
use_env_correlation <- TRUE
t_decomposition_compact_matrix <- inputs$compact_decomposition$matrix
t_decomposition_compact_map <- inputs$compact_decomposition$map
t_decomposition_compact_rows <- nrow(t_decomposition_compact_matrix)
}
} # environmental stochasticity?
# Human impact (harvest) settings
harvest <- ifelse(is.null(inputs$harvest), FALSE, inputs$harvest)
if (harvest) {
harvest_max <- inputs$harvest_max
harvest_g <- inputs$harvest_g
harvest_z <- inputs$harvest_z
harvest_max_n <- inputs$harvest_max_n
human_densities <- inputs$human_density
}
# Results required selection
results_selection <- inputs$results_selection
if (is.null(results_selection)) results_selection <- "abundance"
## Initialization ##
# Initialize carrying capacity
carrying_capacity_t_max <- ncol(carrying_capacities)
if (carrying_capacity_t_max == 1) { # no temporal trend in K
carrying_capacity <- carrying_capacities[, 1]
}
# Intialize results collection list components and other variables
results <- list()
if ("abundance" %in% results_selection) {
results$abundance <- array(0, c(populations, time_steps))
}
if ("ema" %in% results_selection) {
results$ema <- array(0, c(populations, time_steps))
min_abundances <- population_abundances[1,]
}
if ("extirpation" %in% results_selection) {
results$extirpation <- array(NA, populations)
results$extirpation[which(population_abundances == 0)] <- 0
}
if ("harvested" %in% results_selection) {
results$harvested <- array(0, c(populations, time_steps))
}
if ("occupancy" %in% results_selection) {
results$occupancy <- array(0, time_steps)
}
if ("human_density" %in% results_selection) {
results$human_density <- inputs$human_density
}
# Dispersal tracking (for testing/debug purposes)
if ("dispersal_tracking" %in% names(inputs) ||
("attached" %in% names(inputs) && "dispersal_tracking" %in% names(inputs$attached))) {
results$emigrants <- array(0, c(populations, time_steps))
results$immigrants <- array(0, c(populations, time_steps))
dispersal_tracking <- TRUE
} else {
dispersal_tracking <- FALSE
}
## Simulation time steps ##
for (tm in 1:time_steps) {
# Set transitions for each population
transitions <- array(transition_rate, populations)
# Apply occupancy threshold when present
if (!is.null(occupancy_threshold)) {
if (length(which(as.logical(population_abundances))) <= occupancy_threshold) {
population_abundances[] <- 0
}
}
# Load carrying capacity for each population for time if temporal trend in K
if (carrying_capacity_t_max > 1) {
carrying_capacity <- carrying_capacities[, min(tm, carrying_capacity_t_max)]
}
## Density dependence calculations ##
# Selective calculations based on occupied populations and carrying capacity
occupied_indices <- which(as.logical(carrying_capacity*population_abundances[1,])) # > 0
occupied_populations <- length(occupied_indices)
zero_indices <- which(carrying_capacity <= 0 & as.logical(population_abundances))
if (occupied_populations && !is.null(density_dependence) && density_dependence %in% c("competition", "logistic")) {
# Focus on occupied populations
selected_carrying_capacity <- carrying_capacity[occupied_indices]
selected_population_abundances <- population_abundances[occupied_indices]
# Calculate growth rate
if (density_dependence == "competition") { # Beverton-Holt
growth_rate <- growth_rate_max - log(exp(growth_rate_max)*selected_population_abundances/selected_carrying_capacity -
selected_population_abundances/selected_carrying_capacity + 1)
} else if (density_dependence == "logistic") { # Ricker
growth_rate <- growth_rate_max*(1 - selected_population_abundances/selected_carrying_capacity)
}
# Calculate density dependent multipliers and apply to transitions
density_dependence_multipliers <- exp(growth_rate)/transition_rate
transitions[occupied_indices] <- transitions[occupied_indices]*density_dependence_multipliers
# Set any negative values to zero
negative_indices <- occupied_indices[which(transitions[occupied_indices] < 0)]
if (length(negative_indices)) {
transitions[negative_indices] <- 0
}
} # occupied populations?
## Environmental stochasticity calculations ##
if (occupied_populations && environmental_stochasticity) {
# Generate correlated normal deviates for each occupied population (as per Burgman, Ferson & Akcakaya, 1993)
if (use_env_correlation) {
occupied_correlated_deviates <- .colSums(t_decomposition_compact_matrix[, occupied_indices]*stats::rnorm(populations)[t_decomposition_compact_map[, occupied_indices]],
m = t_decomposition_compact_rows, n = occupied_populations, na.rm = TRUE)
} else {
occupied_correlated_deviates <- stats::rnorm(occupied_populations)
}
# Sample from lognormal distribution (as per Burgman, Ferson & Akcakaya, 1993)
if (length(occupied_indices)) {
log_common <- log((standard_deviation/transitions[occupied_indices])^2 + 1)
log_common[which(transitions[occupied_indices] == 0)] <- 0
log_common[which(is.infinite(log_common))] <- 0
transitions[occupied_indices] <- transitions[occupied_indices]*exp(sqrt(log_common)*occupied_correlated_deviates - 0.5*log_common)
}
# Set any negative values to zero
negative_indices <- occupied_indices[which(transitions[occupied_indices] < 0)]
if (length(negative_indices)) {
transitions[negative_indices] <- 0
}
} # occupied populations & environmental stochasticity?
## Generational transition calculations ##
# Remove populations no longer having carrying capacity (determined above)
population_abundances[zero_indices] <- 0
# Sample the next generation's population abundances from a Poisson distribution
if (occupied_populations) {
population_abundances[occupied_indices] <- stats::rpois(occupied_populations, transitions[occupied_indices]*population_abundances[occupied_indices])
}
# Limit abundances to carrying capacities when "ceiling" density dependence
if (occupied_populations && density_dependence == "ceiling") {
above_capacity_indices <- occupied_indices[which(population_abundances[occupied_indices] > carrying_capacity[occupied_indices])]
if (length(above_capacity_indices)) {
population_abundances[above_capacity_indices] <- carrying_capacity[above_capacity_indices]
}
}
## Harvest calculations ##
if (harvest || "harvested" %in% results_selection) {
occupied_indices <- occupied_indices[which(as.logical(population_abundances[occupied_indices]))] # > 0
occupied_populations <- length(occupied_indices)
harvested <- array(0, populations)
}
if (occupied_populations && harvest) {
# Focus on human presence in occupied cells
harvest_rate <- array(0, occupied_populations)
human_presence_indices <- which(as.logical(human_densities[occupied_indices, tm]))
human_presence_occupied_indices <- occupied_indices[human_presence_indices]
# Calculate annual harvest rate
if (length(human_presence_indices)) {
prey_density <- population_abundances[human_presence_occupied_indices]/harvest_max_n
prey_z <- prey_density^harvest_z
max_functional_response <- (harvest_max*prey_z)/(harvest_g + prey_z) # at max human density
functional_response <- max_functional_response*human_densities[human_presence_occupied_indices, tm] # at current human density
harvest_rate[human_presence_indices] <- functional_response/prey_density
}
# Convert to generational time scale
harvest_rate <- 1 - (1 - harvest_rate)^years_per_step
# Calculate numbers harvested
harvested[occupied_indices] <- stats::rbinom(occupied_populations, population_abundances[occupied_indices], harvest_rate)
# Remove harvested from population
population_abundances[occupied_indices] <- population_abundances[occupied_indices] - harvested[occupied_indices]
} # occupied populations & harvest?
## Dispersal calculations ##
if (occupied_populations && dispersal_present) {
# Apply any spatio-temporal dispersal changes
if (tm == 1) {
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
}
# Select dispersals for occupied populations and their non-zero indices
occupied_dispersals <- dispersal_compact_matrix_tm[, occupied_indices]
occupied_dispersal_indices <- which(as.logical(occupied_dispersals)) # > 0
# Modify dispersal rates when dispersal depends on target population carrying capacity K
if (dispersal_depends_on_target_pop_k) {
# Calculate the (below-threshold) multipliers
dd_multipliers <- array(1, populations)
modify_pop_indices <- which(carrying_capacity < dispersal_target_k)
dd_multipliers[modify_pop_indices] <- carrying_capacity[modify_pop_indices]/dispersal_target_k
# 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
}
# Release variables from memory
modify_pop_indices <- NULL; dd_multipliers <- NULL; selected_dd_multipliers <- NULL; modify_indices <- NULL; modify_dipersal_indices <- NULL
} # dispersal depends on target pop k?
# Disperser generation via abundance and corresponding dispersal rates
occupied_abundances <- population_abundances[occupied_indices]
occupied_abundances_rep <- population_abundances[dispersal_abundance_rep_indices, occupied_indices]
dispersers <- array(0, c(dispersal_compact_rows, occupied_populations))
# Generate dispersers using a binomial distribution
dispersers[occupied_dispersal_indices] <- stats::rbinom(length(occupied_dispersal_indices), occupied_abundances_rep[occupied_dispersal_indices], occupied_dispersals[occupied_dispersal_indices])
# Release variables from memory
occupied_dispersals <- NULL; occupied_abundances_rep <- NULL; occupied_dispersal_indices <- NULL
# Calculate emigrants
emigrants <- array(0, occupied_populations)
emigrants[] <- .colSums(dispersers, m = dispersal_compact_rows, n = occupied_populations)
# Check consistency of emigrants (not to exceed abundances)
excessive_indices <- which(emigrants > occupied_abundances)
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_abundances[excessive_index]
# for (remove_row_index in sample(rep(excessive_rows, times = excessive_dispersers), disperser_reduction)) {
# dispersers[remove_row_index, excessive_index] <- dispersers[remove_row_index, excessive_index] - 1
# }
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_abundances[excessive_indices]
}
# Update population abundances
population_abundances[occupied_indices] <- population_abundances[occupied_indices] - emigrants
# Dispersal tracking
if (dispersal_tracking) {
results$emigrants[occupied_indices, tm] <- emigrants
}
# Release variables from memory
occupied_abundances <- NULL; emigrants <- NULL; excessive_indices <- NULL
# Calculate immigrants via dispersal immigrant map
disperser_indices <- which(as.logical(dispersers)) # >0
immigrant_array <- dispersal_zero_array
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
population_abundances[1,] <- population_abundances[1,] + immigrants
# Dispersal tracking
if (dispersal_tracking) {
results$immigrants[, tm] <- immigrants
}
# Release variables from memory
dispersers <- NULL; immigrant_array <- NULL; immigrants <- NULL
} # occupied populations & dispersal present?
## Apply threshold to population abundances ##
if (!is.null(abundance_threshold)) {
below_threshold_indices <- which(as.logical(population_abundances) & population_abundances <= abundance_threshold)
if (length(below_threshold_indices)) {
population_abundances[below_threshold_indices] <- 0
}
}
## Results collection ##
if ("abundance" %in% results_selection) {
results$abundance[, tm] <- population_abundances[1,]
}
if ("ema" %in% results_selection) {
min_abundances <- pmin(min_abundances, population_abundances[1,])
results$ema[,tm] <- min_abundances
}
if ("extirpation" %in% results_selection) {
results$extirpation <- pmin(results$extirpation, rep(tm, populations), na.rm = TRUE)
results$extirpation[which(as.logical(population_abundances))] <- NA
}
if ("harvested" %in% results_selection) {
results$harvested[, tm] <- harvested
}
if ("occupancy" %in% results_selection) {
results$occupancy[tm] <- sum(as.logical(population_abundances))
}
} ### End time steps ###
return(results)
} ### End paleopop ###
Try the paleopop package in your browser
Any scripts or data that you put into this service are public.
paleopop documentation built on Oct. 11, 2023, 5:14 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 paleopop_simulator
Documented in paleopop_simulator
#' Runs a customized population model simulation.
#'
#' Simulates a population model customized for paleontological time-scales, optimized
#' for single-generation transitions and large populations, across multiple generations
#' and returns simulation results. Each generational time-step includes:
#' \enumerate{
#' \item Density dependence calculations
#' \item Environmental stochasticity calculations
#' \item Generational transition calculations
#' \item Harvest calculations
#' \item Dispersal calculations
#' \item Results collection
#' }
#'
#' @param inputs Nested list/object with named elements:
#' \describe{
#' \item{\code{random_seed}}{Number to seed the random number generation for stochasticity.}
#' \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{initial_abundance}}{Array of initial abundances for each population.}
#' \item{\code{transition_rate}}{Rate of transition (or fecundity) between generations.}
#' \item{\code{standard_deviation}}{Standard deviation applied to transition rates.}
#' \item{\code{compact_decomposition}}{List containing a compact transposed (Cholesky) decomposition \emph{matrix} (t_decomposition_compact_matrix) and a corresponding \emph{map} of population indices (t_decomposition_compact_map), as per \code{\link{SpatialCorrelation}} class attributes.}
#' \item{\code{carrying_capacity}}{Matrix of carrying capacities (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{density_dependence}}{Density dependence type ("competition", "logistic", or "ceiling").}
#' \item{\code{growth_rate_max}}{Maximum growth rate (for "competition" or "logistic" density dependence).}
#' \item{\code{harvest}}{Boolean for utilizing harvesting.}
#' \item{\code{harvest_max}}{Proportion harvested per year (note: annual time scale - not generational).}
#' \item{\code{harvest_g}}{The \emph{G} parameter in the harvest function.}
#' \item{\code{harvest_z}}{The \emph{Z} parameter in the harvest function.}
#' \item{\code{harvest_max_n}}{Maximum density per grid cell.}
#' \item{\code{human_density}}{Matrix of human density (fraction) (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{dispersal_data}}{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.}
#' \item{\code{dispersal_target_k}}{Target population carrying capacity threshold for density dependent dispersal.}
#' \item{\code{abundance_threshold}}{Abundance threshold (that needs to be exceeded) for each population to persist.}
#' \item{\code{occupancy_threshold}}{Threshold for the number of populations occupied (that needs to be exceeded) for all populations to persist.}
#' \item{\code{results_selection}}{List of results selection from: "abundance", "ema", "extirpation", "harvested", "occupancy", "human_density".}
#' }
#' @return Simulation results as a nested list (as selected):
#' \describe{
#' \item{\code{abundance}}{Matrix of simulation abundances (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{ema}}{Matrix of expected minimum abundances (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{extirpation}}{Array of extirpation times for each population.}
#' \item{\code{harvested}}{Matrix of estimated individuals harvested (\emph{populations} rows by \emph{time_steps} columns).}
#' \item{\code{occupancy}}{Array of number of populations occupied at each time-step.}
#' \item{\code{human_density}}{Matrix of human densities, (\emph{populations} rows by \emph{time_steps} columns).}
#' }
#'
#' @examples
#' library(raster)
#' library(poems)
#' # Ring Island example region
#' coordinates <- data.frame(x = rep(seq(-178.02, -178.06, -0.01), 5),
#' y = rep(seq(19.02, 19.06, 0.01), each = 5))
#' template_raster <- Region$new(coordinates = coordinates)$region_raster # full extent
#' sealevel_raster <- template_raster
#' template_raster[][c(7:9, 12:14, 17:19)] <- NA # make Ring Island
#' sealevel_raster[][c(7:9, 12:14, 17:18)] <- NA
#' raster_stack <- raster::stack(x = append(replicate(9, template_raster), sealevel_raster))
#' region <- PaleoRegion$new(template_raster = raster_stack)
#'
#' # Model template
#' model_template <- PaleoPopModel$new(
#' region = region,
#' time_steps = 10,
#' years_per_step = 12, # years per generational time-step
#' standard_deviation = 0.1,
#' growth_rate_max = 0.6,
#' harvest = FALSE,
#' populations = region$region_cells,
#' initial_abundance = seq(9000, 0, -1000),
#' transition_rate = 1.0,
#' carrying_capacity = rep(1000, 17),
#' dispersal = (!diag(nrow = 17, ncol = 17))*0.05,
#' density_dependence = "logistic",
#' dispersal_target_k = 10,
#' occupancy_threshold = 1,
#' abundance_threshold = 10,
#' results_selection = c("abundance")
#' )
#'
#' # Simulations
#' results <- paleopop_simulator(model_template) # input as PaleoPopModel object
#' inputs <- model_template$get_attributes()
#' paleopop_simulator(inputs) # input as list of attributes
#'
#' @export paleopop_simulator
paleopop_simulator <- function(inputs) {
## Unpack inputs and calculate re-usable variables ##
# Stop if minimal inputs are not present
if (is.null(inputs$time_steps) || is.null(inputs$populations) || is.null(inputs$initial_abundance) ||
is.null(inputs$transition_rate) || is.null(inputs$carrying_capacity)) {
incomplete_inputs <- if (is.null(inputs$time_steps)) "time_steps"
incomplete_inputs <- c(incomplete_inputs, if (is.null(inputs$populations)) "populations")
incomplete_inputs <- c(incomplete_inputs, if (is.null(inputs$initial_abundance)) "initial_abundance")
incomplete_inputs <- c(incomplete_inputs, if (is.null(inputs$transition_rate)) "transition_rate")
incomplete_inputs <- c(incomplete_inputs, if (is.null(inputs$carrying_capacity)) "carrying_capacity")
stop(paste("Minimal inputs required to run simulation should include:",
paste(incomplete_inputs, collapse = ", ")), call. = FALSE)
}
# General model settings
time_steps <- inputs$time_steps
years_per_step <- ifelse(is.null(inputs$years_per_step), 1, inputs$years_per_step)
random_seed <- inputs$random_seed
if (!is.null(random_seed)) {
set.seed(random_seed)
}
# Transition rate and standard deviation
transition_rate <- inputs$transition_rate
standard_deviation <- ifelse(is.null(inputs$standard_deviation), 0, inputs$standard_deviation)
environmental_stochasticity <- (standard_deviation*transition_rate > 0)
# Population settings
populations <- inputs$populations
population_abundances <- array(inputs$initial_abundance, c(1, populations)) # single column matrix for faster dispersal calculations
density_dependence <- inputs$density_dependence
if (is.null(density_dependence)) {
density_dependence <- "none"
}
growth_rate_max <- ifelse(is.null(inputs$growth_rate_max), log(transition_rate), inputs$growth_rate_max)
abundance_threshold <- inputs$abundance_threshold
occupancy_threshold <- inputs$occupancy_threshold
dispersal_target_k <- inputs$dispersal_target_k
carrying_capacities <- matrix(inputs$carrying_capacity, nrow = populations)
if (any(is.na(carrying_capacities))) {
carrying_capacities[which(is.na(carrying_capacities))] <- 0
}
# Dispersal and correlation
# Create re-usable dispersal structures
dispersal_present <- (!is.null(inputs$dispersal_data) && nrow(inputs$dispersal_data[[1]]))
if (dispersal_present) {
# Unpack dispersal data and determine compact matrix dimensions
dispersal_data <- inputs$dispersal_data[[1]]
dispersal_compact_rows <- max(dispersal_data[, c("emigrant_row", "immigrant_row")])
# Create compact array of zeros for quick initialization
dispersal_zero_array <- array(0, c(dispersal_compact_rows, populations))
# Create a compact matrix of dispersal rates
dispersal_compact_matrix <- dispersal_zero_array
dispersal_compact_matrix[as.matrix(dispersal_data[, c("emigrant_row", "source_pop")])] <- dispersal_data$dispersal_rate
# Are dispersals changing over time?
dispersals_change_over_time <- (length(inputs$dispersal_data) > 1)
if (dispersals_change_over_time) {
dispersal_data_changes <- inputs$dispersal_data
dispersal_data_changes[[1]] <- dispersal_data_changes[[1]][NULL,]
}
# Does dispersal dependent on target population carrying capacity K?
dispersal_depends_on_target_pop_k <- (!is.null(dispersal_target_k) && dispersal_target_k > 0)
if (dispersal_depends_on_target_pop_k) {
# Create a map of compact array indices for mapping dispersers (emigrants) to target populations
dispersal_target_pop_map <- dispersal_zero_array
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 <- dispersal_zero_array
dispersal_immigrant_map[as.matrix(dispersal_data[, c("emigrant_row", "source_pop")])] <- dispersal_compact_indices[as.matrix(dispersal_data[, c("immigrant_row", "target_pop")])]
# Create indices for replicating abundances
dispersal_abundance_rep_indices <- rep(1, dispersal_compact_rows)
# Release variables from memory
dispersal_data <- NULL; dispersal_compact_indices <- NULL
} # dispersal present?
# Environmental correlation?
if (environmental_stochasticity) {
# Flag for environmental correlation
use_env_correlation <- FALSE
# Unpack compact decompostion
if (!is.null(inputs$compact_decomposition)) {
use_env_correlation <- TRUE
t_decomposition_compact_matrix <- inputs$compact_decomposition$matrix
t_decomposition_compact_map <- inputs$compact_decomposition$map
t_decomposition_compact_rows <- nrow(t_decomposition_compact_matrix)
}
} # environmental stochasticity?
# Human impact (harvest) settings
harvest <- ifelse(is.null(inputs$harvest), FALSE, inputs$harvest)
if (harvest) {
harvest_max <- inputs$harvest_max
harvest_g <- inputs$harvest_g
harvest_z <- inputs$harvest_z
harvest_max_n <- inputs$harvest_max_n
human_densities <- inputs$human_density
}
# Results required selection
results_selection <- inputs$results_selection
if (is.null(results_selection)) results_selection <- "abundance"
## Initialization ##
# Initialize carrying capacity
carrying_capacity_t_max <- ncol(carrying_capacities)
if (carrying_capacity_t_max == 1) { # no temporal trend in K
carrying_capacity <- carrying_capacities[, 1]
}
# Intialize results collection list components and other variables
results <- list()
if ("abundance" %in% results_selection) {
results$abundance <- array(0, c(populations, time_steps))
}
if ("ema" %in% results_selection) {
results$ema <- array(0, c(populations, time_steps))
min_abundances <- population_abundances[1,]
}
if ("extirpation" %in% results_selection) {
results$extirpation <- array(NA, populations)
results$extirpation[which(population_abundances == 0)] <- 0
}
if ("harvested" %in% results_selection) {
results$harvested <- array(0, c(populations, time_steps))
}
if ("occupancy" %in% results_selection) {
results$occupancy <- array(0, time_steps)
}
if ("human_density" %in% results_selection) {
results$human_density <- inputs$human_density
}
# Dispersal tracking (for testing/debug purposes)
if ("dispersal_tracking" %in% names(inputs) ||
("attached" %in% names(inputs) && "dispersal_tracking" %in% names(inputs$attached))) {
results$emigrants <- array(0, c(populations, time_steps))
results$immigrants <- array(0, c(populations, time_steps))
dispersal_tracking <- TRUE
} else {
dispersal_tracking <- FALSE
}
## Simulation time steps ##
for (tm in 1:time_steps) {
# Set transitions for each population
transitions <- array(transition_rate, populations)
# Apply occupancy threshold when present
if (!is.null(occupancy_threshold)) {
if (length(which(as.logical(population_abundances))) <= occupancy_threshold) {
population_abundances[] <- 0
}
}
# Load carrying capacity for each population for time if temporal trend in K
if (carrying_capacity_t_max > 1) {
carrying_capacity <- carrying_capacities[, min(tm, carrying_capacity_t_max)]
}
## Density dependence calculations ##
# Selective calculations based on occupied populations and carrying capacity
occupied_indices <- which(as.logical(carrying_capacity*population_abundances[1,])) # > 0
occupied_populations <- length(occupied_indices)
zero_indices <- which(carrying_capacity <= 0 & as.logical(population_abundances))
if (occupied_populations && !is.null(density_dependence) && density_dependence %in% c("competition", "logistic")) {
# Focus on occupied populations
selected_carrying_capacity <- carrying_capacity[occupied_indices]
selected_population_abundances <- population_abundances[occupied_indices]
# Calculate growth rate
if (density_dependence == "competition") { # Beverton-Holt
growth_rate <- growth_rate_max - log(exp(growth_rate_max)*selected_population_abundances/selected_carrying_capacity -
selected_population_abundances/selected_carrying_capacity + 1)
} else if (density_dependence == "logistic") { # Ricker
growth_rate <- growth_rate_max*(1 - selected_population_abundances/selected_carrying_capacity)
}
# Calculate density dependent multipliers and apply to transitions
density_dependence_multipliers <- exp(growth_rate)/transition_rate
transitions[occupied_indices] <- transitions[occupied_indices]*density_dependence_multipliers
# Set any negative values to zero
negative_indices <- occupied_indices[which(transitions[occupied_indices] < 0)]
if (length(negative_indices)) {
transitions[negative_indices] <- 0
}
} # occupied populations?
## Environmental stochasticity calculations ##
if (occupied_populations && environmental_stochasticity) {
# Generate correlated normal deviates for each occupied population (as per Burgman, Ferson & Akcakaya, 1993)
if (use_env_correlation) {
occupied_correlated_deviates <- .colSums(t_decomposition_compact_matrix[, occupied_indices]*stats::rnorm(populations)[t_decomposition_compact_map[, occupied_indices]],
m = t_decomposition_compact_rows, n = occupied_populations, na.rm = TRUE)
} else {
occupied_correlated_deviates <- stats::rnorm(occupied_populations)
}
# Sample from lognormal distribution (as per Burgman, Ferson & Akcakaya, 1993)
if (length(occupied_indices)) {
log_common <- log((standard_deviation/transitions[occupied_indices])^2 + 1)
log_common[which(transitions[occupied_indices] == 0)] <- 0
log_common[which(is.infinite(log_common))] <- 0
transitions[occupied_indices] <- transitions[occupied_indices]*exp(sqrt(log_common)*occupied_correlated_deviates - 0.5*log_common)
}
# Set any negative values to zero
negative_indices <- occupied_indices[which(transitions[occupied_indices] < 0)]
if (length(negative_indices)) {
transitions[negative_indices] <- 0
}
} # occupied populations & environmental stochasticity?
## Generational transition calculations ##
# Remove populations no longer having carrying capacity (determined above)
population_abundances[zero_indices] <- 0
# Sample the next generation's population abundances from a Poisson distribution
if (occupied_populations) {
population_abundances[occupied_indices] <- stats::rpois(occupied_populations, transitions[occupied_indices]*population_abundances[occupied_indices])
}
# Limit abundances to carrying capacities when "ceiling" density dependence
if (occupied_populations && density_dependence == "ceiling") {
above_capacity_indices <- occupied_indices[which(population_abundances[occupied_indices] > carrying_capacity[occupied_indices])]
if (length(above_capacity_indices)) {
population_abundances[above_capacity_indices] <- carrying_capacity[above_capacity_indices]
}
}
## Harvest calculations ##
if (harvest || "harvested" %in% results_selection) {
occupied_indices <- occupied_indices[which(as.logical(population_abundances[occupied_indices]))] # > 0
occupied_populations <- length(occupied_indices)
harvested <- array(0, populations)
}
if (occupied_populations && harvest) {
# Focus on human presence in occupied cells
harvest_rate <- array(0, occupied_populations)
human_presence_indices <- which(as.logical(human_densities[occupied_indices, tm]))
human_presence_occupied_indices <- occupied_indices[human_presence_indices]
# Calculate annual harvest rate
if (length(human_presence_indices)) {
prey_density <- population_abundances[human_presence_occupied_indices]/harvest_max_n
prey_z <- prey_density^harvest_z
max_functional_response <- (harvest_max*prey_z)/(harvest_g + prey_z) # at max human density
functional_response <- max_functional_response*human_densities[human_presence_occupied_indices, tm] # at current human density
harvest_rate[human_presence_indices] <- functional_response/prey_density
}
# Convert to generational time scale
harvest_rate <- 1 - (1 - harvest_rate)^years_per_step
# Calculate numbers harvested
harvested[occupied_indices] <- stats::rbinom(occupied_populations, population_abundances[occupied_indices], harvest_rate)
# Remove harvested from population
population_abundances[occupied_indices] <- population_abundances[occupied_indices] - harvested[occupied_indices]
} # occupied populations & harvest?
## Dispersal calculations ##
if (occupied_populations && dispersal_present) {
# Apply any spatio-temporal dispersal changes
if (tm == 1) {
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
}
# Select dispersals for occupied populations and their non-zero indices
occupied_dispersals <- dispersal_compact_matrix_tm[, occupied_indices]
occupied_dispersal_indices <- which(as.logical(occupied_dispersals)) # > 0
# Modify dispersal rates when dispersal depends on target population carrying capacity K
if (dispersal_depends_on_target_pop_k) {
# Calculate the (below-threshold) multipliers
dd_multipliers <- array(1, populations)
modify_pop_indices <- which(carrying_capacity < dispersal_target_k)
dd_multipliers[modify_pop_indices] <- carrying_capacity[modify_pop_indices]/dispersal_target_k
# 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
}
# Release variables from memory
modify_pop_indices <- NULL; dd_multipliers <- NULL; selected_dd_multipliers <- NULL; modify_indices <- NULL; modify_dipersal_indices <- NULL
} # dispersal depends on target pop k?
# Disperser generation via abundance and corresponding dispersal rates
occupied_abundances <- population_abundances[occupied_indices]
occupied_abundances_rep <- population_abundances[dispersal_abundance_rep_indices, occupied_indices]
dispersers <- array(0, c(dispersal_compact_rows, occupied_populations))
# Generate dispersers using a binomial distribution
dispersers[occupied_dispersal_indices] <- stats::rbinom(length(occupied_dispersal_indices), occupied_abundances_rep[occupied_dispersal_indices], occupied_dispersals[occupied_dispersal_indices])
# Release variables from memory
occupied_dispersals <- NULL; occupied_abundances_rep <- NULL; occupied_dispersal_indices <- NULL
# Calculate emigrants
emigrants <- array(0, occupied_populations)
emigrants[] <- .colSums(dispersers, m = dispersal_compact_rows, n = occupied_populations)
# Check consistency of emigrants (not to exceed abundances)
excessive_indices <- which(emigrants > occupied_abundances)
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_abundances[excessive_index]
# for (remove_row_index in sample(rep(excessive_rows, times = excessive_dispersers), disperser_reduction)) {
# dispersers[remove_row_index, excessive_index] <- dispersers[remove_row_index, excessive_index] - 1
# }
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_abundances[excessive_indices]
}
# Update population abundances
population_abundances[occupied_indices] <- population_abundances[occupied_indices] - emigrants
# Dispersal tracking
if (dispersal_tracking) {
results$emigrants[occupied_indices, tm] <- emigrants
}
# Release variables from memory
occupied_abundances <- NULL; emigrants <- NULL; excessive_indices <- NULL
# Calculate immigrants via dispersal immigrant map
disperser_indices <- which(as.logical(dispersers)) # >0
immigrant_array <- dispersal_zero_array
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
population_abundances[1,] <- population_abundances[1,] + immigrants
# Dispersal tracking
if (dispersal_tracking) {
results$immigrants[, tm] <- immigrants
}
# Release variables from memory
dispersers <- NULL; immigrant_array <- NULL; immigrants <- NULL
} # occupied populations & dispersal present?
## Apply threshold to population abundances ##
if (!is.null(abundance_threshold)) {
below_threshold_indices <- which(as.logical(population_abundances) & population_abundances <= abundance_threshold)
if (length(below_threshold_indices)) {
population_abundances[below_threshold_indices] <- 0
}
}
## Results collection ##
if ("abundance" %in% results_selection) {
results$abundance[, tm] <- population_abundances[1,]
}
if ("ema" %in% results_selection) {
min_abundances <- pmin(min_abundances, population_abundances[1,])
results$ema[,tm] <- min_abundances
}
if ("extirpation" %in% results_selection) {
results$extirpation <- pmin(results$extirpation, rep(tm, populations), na.rm = TRUE)
results$extirpation[which(as.logical(population_abundances))] <- NA
}
if ("harvested" %in% results_selection) {
results$harvested[, tm] <- harvested
}
if ("occupancy" %in% results_selection) {
results$occupancy[tm] <- sum(as.logical(population_abundances))
}
} ### End time steps ###
return(results)
} ### End paleopop ###
Try the paleopop 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.