Nothing
R/generics.R
In ipmr: Integral Projection Models
Defines functions
.all_params_in_mod
print.proto_ipm
Documented in
print.proto_ipm
# Print----------------
#' @title Print proto_ipms or *_ipm objects
#' @rdname print_star
#'
#' @param x An object of class \code{proto_ipm} or produced by \code{make_ipm()}.
#' @param ... Ignored
#'
#' @details For printing \code{proto_ipm} objects, indices are wrapped in
#' \code{<index>} to assist with debugging. These are not carried into the model,
#' just a visual aid.
#'
#' @export
print.proto_ipm <- function(x, ...) {
n_kerns <- dim(x)[1]
cls_switch <- class(x)[1]
pretty_class <- .pretty_class(cls_switch)
msg <- paste(ifelse(.uses_age(x),
"An age structured",
"A"),
pretty_class,
"proto_ipm with",
n_kerns,
"kernels defined:\n")
msg <- c(msg, paste(x$kernel_id, collapse = ', '))
msg <- c(msg, '\n\nKernel formulae:\n\n')
cat(msg)
forms <- kernel_formulae(x)
print(forms)
cat("\n\nVital rates:\n\n")
vrs <- vital_rate_exprs(x)
print(vrs)
cat("\n\nParameter names:\n\n")
pars <- parameters(x)
print(names(pars))
# Check if all parameters are in vital rate expressions, and if all
# vital rate expression args are found in parameters.
in_mod <- .all_params_in_mod(x)
in_mod_tst <- vapply(in_mod, function(x) x$out, logical(1L)) %>%
all()
cat("\nAll parameters in vital rate expressions found in 'data_list': ",
in_mod_tst)
if(!in_mod_tst) {
stoch_params <- c(paste("simple_d", c('i', 'd'), "_stoch_param", sep = ""),
paste("general_d", c('i', 'd'), "_stoch_param", sep = ""))
if(inherits(x, stoch_params)) {
cat("\nDid not test parameters in 'define_env_state()'.\n")
}
cat("\nItems in vital rate expressions not found in 'data_list':\n")
for(i in seq_along(in_mod)) {
if(!in_mod[[i]]$out) {
nm <- paste("\n*", names(in_mod)[i], "*\n", sep = "")
miss_msg <- paste(nm,
paste(in_mod[[i]]$missing, collapse = "\n"),
sep = "")
cat(miss_msg, '\n')
}
}
cat("\nThis may be because 'define_domains()' has not yet been called.\n",
"If it has been called, then double check the model definition and",
"'data_list'.\n")
}
cat("\n\nDomains for state variables:\n\n")
print(domains(x))
cat("\n\nPopulation states defined:\n\n")
print(pop_state(x))
cat("\n\nInternally generated model iteration procedure:\n\n")
if(!all(is.na(x$domain))){
k_row <- .init_iteration(x, TRUE, direction = "right")
print(kernel_formulae(k_row))
} else {
cat("Not enough information to generate an iteration procedure.\n",
"define_impl() and define_pop_state() must be called first to generate one.")
}
# Add more later -----------
cat("\n\n")
invisible(x)
}
#' @noRd
.all_params_in_mod <- function(proto) {
kern_list <- list()
for(row in seq_len(nrow(proto))) {
use_proto <- proto[row, ]
kern_nm <- proto$kernel_id[row]
vr_exprs <- vital_rate_exprs(use_proto)
params <- parameters(proto)
# Age x size models and/or par_sets models need expansion before we
# we can test whether all parameters are in the model
if(.uses_age(proto) | any(proto$uses_par_sets)) {
all_levs <- .flatten_to_depth(proto$par_set_indices[proto$uses_par_sets],
1L) %>%
.[!duplicated(names(.))]
par_set_nms <- names(all_levs)
if(length(all_levs) > 1) {
lev_df <- expand.grid(all_levs,
stringsAsFactors = FALSE)
} else {
lev_df <- as.data.frame(all_levs)
}
# Hold the results
all_vr_exprs <- list()
it <- 1
# Replace suffixes in expressions
for(i in seq_along(vr_exprs)) {
base_expr <- rlang::expr_text(vr_exprs[[i]])
for(j in seq_len(dim(lev_df)[2])) {
par_set_nm <- names(lev_df)[j]
for(k in seq_len(dim(lev_df)[1])) {
# Replace suffixes in expression and name
all_vr_exprs[[it]] <- gsub(par_set_nm,
lev_df[k, j],
base_expr)
names(all_vr_exprs)[it] <- gsub(par_set_nm,
lev_df[k, j],
names(vr_exprs)[i])
it <- it + 1
}
}
}
vr_exprs <- lapply(all_vr_exprs, rlang::parse_expr)
vr_args <- .vr_args(vr_exprs)
} else {
vr_args <- .vr_args(vr_exprs)
}
# If a vital rate expression creates another parameter, then remove that
# from this test, as it shouldn't appear in the data_list anyway.
to_exclude <- names(vr_exprs)
domains <- names(domains(proto)) %>%
vapply(function(x) paste(x, c("_1", "_2"), sep = ""),
character(2L)) %>%
as.vector()
vr_args <- Filter(.can_be_number, vr_args)
out <- vapply(vr_args[!vr_args %in% to_exclude],
function(nm, pars) nm %in% pars,
logical(1L),
pars = c(names(params), domains))
if(!all(out)) {
vr_args <- vr_args[!vr_args %in% to_exclude]
ind <- which(!vr_args %in% c(names(params), domains))
miss_args <- vr_args[ind]
out <- FALSE
} else {
out <- TRUE
miss_args <- NA_character_
}
kern_list <- c(kern_list,
list(list(out = out,
missing = miss_args)))
names(kern_list)[row] <- kern_nm
}
return(kern_list)
}
.can_be_number <- function(x) {
suppressWarnings(is.na(as.numeric(x)))
}
#' @noRd
.vr_args <- function(vr_exprs) {
all_args <- lapply(vr_exprs, function(y) {
temp <- rlang::expr_text(y)
.args_from_txt(temp)
}) %>%
unlist() %>%
unique()
return(all_args)
}
#' @noRd
.uses_age <- function(x) {
inherits(x, "age_x_size")
}
#' @noRd
.pretty_class <- function(cls_switch) {
out <-
switch(
cls_switch,
# IPM classes -----------------------------
'simple_di_det_ipm' = "simple, density independent, deterministic",
'simple_di_stoch_kern_ipm' = "simple, density independent, stochastic, kernel-resampled",
'simple_di_stoch_param_ipm' = "simple, density independent, stochastic, parameter-resampled",
'simple_dd_det_ipm' = "simple, density dependent, deterministic",
'simple_dd_stoch_kern_ipm' = "simple, density dependent, stochastic, kernel-resampled",
'simple_dd_stoch_param_ipm' = "simple, density dependent, stochastic, parameter-resampled",
'general_di_det_ipm' = "general, density independent, deterministic",
'general_di_stoch_kern_ipm' = "general, density independent, stochastic, kernel-resampled",
'general_di_stoch_param_ipm' = "general, density independent, stochastic, parameter-resampled",
'general_dd_det_ipm' = "general, density dependent, deterministic",
'general_dd_stoch_kern_ipm' = "general, density dependent, stochastic, kernel-resampled",
'general_dd_stoch_param_ipm' = "general, density dependent, stochastic, parameter-resampled",
# proto_ipm classes -----------------------
'simple_di_det' = "simple, density independent, deterministic",
'simple_di_stoch_kern' = "simple, density independent, stochastic, kernel-resampled",
'simple_di_stoch_param' = "simple, density independent, stochastic, parameter-resampled",
'simple_dd_det' = "simple, density dependent, deterministic",
'simple_dd_stoch_kern' = "simple, density dependent, stochastic, kernel-resampled",
'simple_dd_stoch_param' = "simple, density dependent, stochastic, parameter-resampled",
'general_di_det' = "general, density independent, deterministic",
'general_di_stoch_kern' = "general, density independent, stochastic, kernel-resampled",
'general_di_stoch_param' = "general, density independent, stochastic, parameter-resampled",
'general_dd_det' = "general, density dependent, deterministic",
'general_dd_stoch_kern' = "general, density dependent, stochastic, kernel-resampled",
'general_dd_stoch_param' = "general, density dependent, stochastic, parameter-resampled"
)
return(out)
}
#' @rdname print_star
#' @title Generics for IPM classes
#'
#' @param comp_lambda A logical indicating whether or not to calculate lambdas
#' for the iteration kernels and display them.
#' @param sig_digits The number of significant digits to round to if \code{
#' comp_lambda = TRUE}.
#' @param type_lambda Either \code{'all'} or \code{'stochastic'}. See
#' \code{\link{lambda}} for more details.
#' @param check_conv A logical: for deterministic models, check if population state
#' has converged to asymptotic dynamics? If \code{TRUE} and the model has not
#' converged, a message will be printed.
#' @param ... Ignored
#'
#' @return \code{x} invisibly.
#'
#' @export
print.simple_di_det_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'last',
sig_digits = 3,
check_conv = TRUE,
...) {
pretty_cls <- .pretty_class(class(x)[1])
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) defined.', sep = "")
if(comp_lambda) {
lambdas <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nDeterministic lambda = ',
round(lambdas, sig_digits),
sep = "")
msg <- c(msg, l_msg)
if(check_conv &&
!all(is_conv_to_asymptotic(x))) {
# Captures the name of the model that the user gave rather than
# just print "x isn't converged"
mod_nm <- deparse(substitute(x))
message(
paste(mod_nm,
' has has not converged to asymptotic dynamics!',
sep = "")
)
}
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.simple_dd_det_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'last',
sig_digits = 3,
...) {
mod_nm <- deparse(substitute(x))
pretty_cls <- .pretty_class(class(x)[1])
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) defined.', sep = "")
if(comp_lambda) {
lambdas <- lambda(x, type_lambda = type_lambda)
l_msg <- paste('\nLambda for the final time step of the model is: ',
round(lambdas, sig_digits),
sep = "")
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.simple_di_stoch_kern_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'stochastic',
sig_digits = 3,
...) {
mod_nm <- deparse(substitute(x))
pretty_cls <- .pretty_class(class(x)[1])
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$iterators),
' iteration kernel(s) and ',
length(x$sub_kernels),
' sub-kernel(s) defined.', sep = "")
if(comp_lambda) {
all_lams <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nStochastic lambda for ',
mod_nm,
' = ',
round(all_lams, sig_digits),
sep = "" )
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.simple_dd_stoch_kern_ipm <- print.simple_di_stoch_kern_ipm
#' @rdname print_star
#' @export
print.simple_di_stoch_param_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'stochastic',
sig_digits = 3,
...) {
pretty_cls <- .pretty_class(class(x)[1])
mod_nm <- deparse(substitute(x))
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$iterators),
' iteration kernel(s) and ',
length(x$sub_kernels),
' sub-kernel(s) defined.', sep = "")
if(comp_lambda) {
all_lams <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nStochastic lambda for ',
mod_nm,
' = ',
round(all_lams, sig_digits),
sep = "" )
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.simple_dd_stoch_param_ipm <- print.simple_di_stoch_param_ipm
#' @rdname print_star
#' @export
print.general_di_det_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'last',
sig_digits = 3,
check_conv = TRUE,
...) {
pretty_cls <- .pretty_class(class(x)[1])
pops <- x$pop_state[!grepl("lambda", names(x$pop_state))]
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) and ',
length(pops),
' population vectors defined.', sep = "")
mod_nm <- deparse(substitute(x))
if(comp_lambda) {
ret_lam <- lambda(x, type_lambda = type_lambda)
l_msg <- paste('\nLambda for the final time step of the model is: ',
round(ret_lam, sig_digits),
'\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
}
if(check_conv && !all(is_conv_to_asymptotic(x))) {
message(
paste(mod_nm,
' has has not converged to asymptotic dynamics!',
sep = "")
)
}
msg <- c(msg, l_msg)
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.general_dd_det_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'last',
sig_digits = 3,
...) {
pretty_cls <- .pretty_class(class(x)[1])
pops <- x$pop_state[!grepl("lambda", names(x$pop_state))]
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) and ',
length(pops),
' population vectors defined.', sep = "")
mod_nm <- deparse(substitute(x))
if(comp_lambda) {
ret_lam <- lambda(x, type_lambda = type_lambda)
l_msg <- paste('\nLambda for the final time step of the model is: ',
round(ret_lam, sig_digits),
'\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg)
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.general_di_stoch_kern_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'stochastic',
sig_digits = 3,
...) {
pretty_cls <- .pretty_class(class(x)[1])
mod_nm <- deparse(substitute(x))
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) and ',
length(x$pop_state[!grepl("lambda", names(x$pop_state))]),
' population vectors defined.',
sep = "")
if(comp_lambda) {
all_lams <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nStochastic lambda for ',
mod_nm,
' = ',
round(all_lams, sig_digits),
sep = "" )
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
# No check for convergence here - there may quite reasonably be no
# convergence for
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.general_dd_stoch_kern_ipm <- print.general_di_stoch_kern_ipm
#' @rdname print_star
#' @export
print.general_di_stoch_param_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'stochastic',
sig_digits = 3,
...) {
pretty_cls <- .pretty_class(class(x)[1])
mod_nm <- deparse(substitute(x))
msg <- paste0(
'A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s), ',
length(x$pop_state[!grepl("lambda", names(x$pop_state))]),
' population vectors, and ',
dim(x$env_seq)[2],
' environmental parameters defined.',
sep = ""
)
env_nms <- vapply(dimnames(x$env_seq)[[2]],
FUN = function(x) strsplit(x, '\\.')[[1]][2],
FUN.VALUE = character(1L))
msg <- c(msg,
paste(
'\nEnvironmental variables are: ',
paste(
paste(env_nms, collapse = ', ')
)
))
if(comp_lambda) {
all_lams <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nStochastic lambda for ',
mod_nm,
' = ',
round(all_lams, sig_digits),
sep = "" )
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
# No check for convergence here - there may quite reasonably be no
# convergence for
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.general_dd_stoch_param_ipm <- print.general_di_stoch_param_ipm
#' @export
print.ipmr_vital_rate_exprs <- function(x, ...) {
proto <- attr(x, "proto")
if(any(proto$uses_par_sets) | .uses_age(proto)) {
out <- .pretty_print_par_sets(x, proto)
} else {
out <- lapply(x, rlang::expr_text)
}
attr(out, "proto") <- NULL
out <- paste(names(out), out, sep = ": ") %>%
paste(., collapse = "\n")
class(out) <- "character"
cat(out)
invisible(x)
}
.pretty_print_par_sets <- function(x, proto) {
nm_par_sets <- proto$par_set_indices %>%
.flatten_to_depth(1L) %>%
names() %>%
unique()
if(.uses_age(proto)) {
ages <- proto$age_indices %>%
.flatten_to_depth(1L) %>%
names() %>%
unique()
nm_par_sets <- c(nm_par_sets, ages)
}
nm_par_sets <- nm_par_sets[!is.na(nm_par_sets)]
x <- lapply(x, rlang::expr_text)
for(i in seq_along(nm_par_sets)) {
to_print <- paste("<", nm_par_sets[i], ">", sep = "")
x <- lapply(x,
function(y, nm, to_print) {
gsub(nm, to_print, y)
},
nm = nm_par_sets[i],
to_print = to_print)
names(x) <- gsub(nm_par_sets[i], to_print, names(x))
}
return(x)
}
#' @export
print.ipmr_kernel_exprs <- function(x, ...) {
proto <- attr(x, "proto")
if(any(proto$uses_par_sets) | .uses_age(proto)) {
out <- .pretty_print_par_sets(x, proto)
} else {
out <- lapply(x, rlang::expr_text)
}
attr(out, "proto") <- NULL
out <- paste(names(out), out, sep = ": ") %>%
paste(., collapse = "\n")
class(out) <- "character"
cat(out)
invisible(x)
}
#' @export
print.ipmr_parameters <- function(x, ...) {
out <- as.numeric(x)
names(out) <- names(x)
print(out)
invisible(x)
}
#' @export
# Will need to update w/ other names for other integration rules as those
# become available
print.ipmr_domains <- function(x, ...) {
# Strips class and proto attributes
out <- lapply(x, function(y) {
y
})
out <- paste(names(out), out, sep = ": ") %>%
paste(., collapse = "\n")
out <- gsub("c\\(", "", out)
out <- gsub("\\)", "", out)
cat(out)
invisible(x)
}
#' @export
print.ipmr_pop_state <- function(x, ...) {
out <- lapply(x, function(y) y)
if(out[[1]] == "No population state defined.") {
cat(out[[1]])
} else {
if(rlang::is_named(x) && !grepl("pop_state_", names(out))) {
message("Names of population states don't appear to be prefixed with 'n_'.")
}
names(out) <- gsub("pop_state_", "n_", names(out))
out <- paste(names(out), out, sep = ": ") %>%
paste(., collapse = "\n")
cat(out)
}
invisible(x)
}
#' @export
print.ipmr_mp_mesh <- function(x, ...) {
uses <- Filter(Negate(is.null), x)
msg <- vapply(uses, function(z) {
if(length(z) > 1) {
paste("Lower bound: ", range(z)[1], ", Upper bound: ", range(z)[2],
", # of Meshpoints: ", sqrt(length(z)),
sep = "")
} else {
as.character(z)
}
},
character(1L)) %>%
paste(names(.), " - ", ., sep = "") %>%
paste(., collapse = "\n")
cat(msg)
invisible(x)
}
#' @export
print.ipmr_matrix <- function(x, ...) {
rng <- range(x)
msg <- paste("Minimum value: ",
round(rng[1], 5),
", maximum value: ",
round(rng[2], 5), sep = "")
pp_tst <- all(x >= 0)
msg <- c(msg,
paste("\nAll entries greater than or equal to 0: ",
pp_tst ,
sep = ""))
cat("\n", msg, "\n")
invisible(x)
}
#' @export
print.ipmr_vital_rate_funs <- function(x, ...) {
if(length(x) == 1 && x == "No vital rate functions specified") {
cat("No vital rates specified\n")
invisible(x)
return()
}
vr_nms <- names(x)
cls <- vapply(x,
function(x) .mat_cls(x),
character(1L))
rngs <- vapply(x,
function(y) round(range(y), 4),
numeric(2L))
out <- vapply(seq_along(cls),
function(it, vr_nm, cls, rngs) {
paste(
paste(vr_nms[it],
" (not yet discretized)",
sep = ""),
cls[it],
sep = ": "
) %>%
paste(" with minimum value: ", rngs[1, it],
" and maximum value: ", rngs[2, it],
sep = "")
},
vr_nm = vr_nms, cls = cls, rngs = rngs,
FUN.VALUE = character(1L)) %>%
paste(., collapse = "\n")
out <- paste(out, '\n', sep = "")
cat(out)
invisible(x)
}
.pretty_dim <- function(x) {
col <- ncol(x)
row <- nrow(x)
cls <- switch(class(x)[1],
"CC" = "kernel",
"DC" = "column vector",
"CD" = "row vector",
"DD" = "matrix")
paste("A", row, "x", col, cls, sep = " ")
}
.mat_cls <- function(x) {
.pretty_dim(x)
}
# Lambda------------
#' @title Compute the per-capita growth rate for an IPM object
#' @rdname lambda
#'
#' @description Compute the per-capita growth rate for a given model. Can handle
#' stochastic and deterministic models, and has the option to discard burn in for
#' stochastic models.
#'
#' @param ipm An object returned by \code{make_ipm()}.
#' @param type_lambda Either \code{'all'}, \code{'last'},
#' or \code{'stochastic'}. \code{'all'}
#' returns a vector of lambda values for each time step of the simulation (equal
#' in length to the \code{iterations} argument of \code{make_ipm()}).
#' \code{'last'} returns the lambda value for the final timestep.
#' \code{'stochastic'} returns a single value, which by default is
#' \code{mean(log(lambda(ipm, type_lambda = "all")))}, with the proportion of
#' \code{burn_in} iterations removed from the beginning of the simulation. Set
#' \code{log} to \code{FALSE} to get \code{lambda} on the linear scale for
#' stochastic models (i.e. \code{exp(mean(log(lambdas)))}).
#' @param ... other arguments passed to methods.
#' @param burn_in The proportion of iterations to discard. Default is 0.1
#' (i.e. first 10\% of iterations in the simulation).
#' @param log Return lambda on the log scale? This is \code{TRUE} by default for
#' stochastic models, and \code{FALSE} for deterministic models.
#'
#' @return When \code{type_lambda = "all"}, an array. Rows correspond to time
#' steps, and columns correspond to parameter sets (if any). For other types,
#' a numeric vector.
#'
#' @export
lambda <- function(ipm, ...) {
UseMethod('lambda')
}
#' @rdname lambda
#' @export
lambda.simple_di_det_ipm <- function(ipm,
type_lambda = 'last',
log = FALSE,
...) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
"all" = TRUE,
'last' = FALSE)
out <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(log) out <- log(out)
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_di_stoch_kern_ipm <- function(ipm,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL,
...) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
"all" = TRUE,
"stochastic" = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_di_stoch_param_ipm <- function(ipm,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL,
...) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
"all" = TRUE,
"stochastic" = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @importFrom methods hasArg
#' @export
lambda.general_di_det_ipm <- function(ipm,
type_lambda = 'last',
log = FALSE,
...) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'last' = FALSE,
'stochastic' = stop("Cannot compute stochastic lambda for deterministic IPM",
call. = FALSE))
out <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(log) out <- log(out)
return(out)
}
#' @rdname lambda
#' @export
lambda.general_di_stoch_kern_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.general_di_stoch_param_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_dd_det_ipm <- function(ipm,
type_lambda = "all",
...,
log = FALSE) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(
type_lambda,
'all' = TRUE,
'last' = FALSE,
'stochastic' = stop("Cannot compute stochastic lambda for deterministic IPM",
call. = FALSE))
out <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(log) out <- log(out)
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_dd_stoch_kern_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_dd_stoch_param_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.general_dd_det_ipm <- function(ipm,
type_lambda = 'last',
...,
log = FALSE) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'last' = FALSE,
'stochastic' = stop("Cannot compute stochastic lambda for deterministic IPM",
call. = FALSE))
out <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(log) out <- log(out)
return(out)
}
#' @rdname lambda
#' @export
lambda.general_dd_stoch_kern_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.general_dd_stoch_param_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
# Plot ---------------
# Authors - whoever wrote the code from the IPM book
#' @title Plot a matrix or an *_ipm object
#' @rdname plot_star
#'
#' @param x,y Either the values of the meshpoints or \code{NULL}. If \code{NULL},
#' then a sequence is generated so that meshpoints are given sequential bin numbers.
#' @param A,ipm A matrix or a result from \code{make_ipm}, or \code{NULL} if \code{x}
#' is specified as the matrix or IPM object.
#' @param col A vector of colors to use for plotting
#' @param bw A logical indicating whether to use a greyscale palette for plotting
#' @param do_contour A logical indicating whether or not draw contour lines
#' on the plot
#' @param contour_cex A numeric specifying how large to make labels for the
#' contour lines.
#' @param do_legend A logical indicating whether to draw a legend for the plot
#' @param ... further arguments passed to legend
#'
#' @return \code{A} or \code{ipm} invisibly
#'
#' @details
#' If an IPM kernel is overwhelmed by information in say, a fecundity sub-kernel,
#' use the \code{exponent} argument in \code{plot.*_ipm} to make it more visually
#' appealing.
#'
#' @importFrom grDevices grey rainbow
#' @importFrom graphics abline axis contour image layout par
#' @export
plot.ipmr_matrix <- function(x = NULL, y = NULL,
A,
col = grDevices::rainbow(100, start=0.67, end=0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
contour_cex = 1,
...) {
if(missing(A)) A <- x
x <- seq_len(ncol(A))
y <- seq_len(nrow(A))
nx = length(x)
ny = length(y)
x1 = c(1.5 * x[1] - 0.5 * x[2], 1.5 * x[nx] - 0.5 * x[nx - 1])
y1 = c(1.5 * y[1] - 0.5 * y[2], 1.5 * y[ny] - 0.5 * y[ny - 1])
if(bw) col = grDevices::grey( (200:50) / 200 )
graphics::image(list(x = x,
y = y,
z = t(A)),
xlim = x1,
ylim = rev(y1),
col = col,
bty = "u",
...)
graphics::abline(v = range(x1))
graphics::abline(h = range(y1))
if(do_contour) graphics::contour(x,
y,
t(A),
nlevels = 5,
labcex = contour_cex,
add = TRUE)
if(do_legend) {
l.y = seq(min(A), max(A),length = 100)
graphics::image(list(x = 1:2,
y = l.y,
z = rbind(l.y, l.y)),
col = col,
bty = "o",
xaxt = "n",
yaxt = "n")
graphics::axis(side = 2,
cex.axis = 1.5,
at = pretty(seq(min(A), max(A), length=10)))
}
invisible(A)
}
#' @rdname plot_star
#' @param exponent The exponent to raise each kernel to. Setting this to a low
#' number can help visualize kernels that are overwhelmed by a few very large
#' numbers.
#' @param n_row,n_col If plotting multiple (sub-)kernels, how many rows and
#' columns to arrange them in.
#'
#' @export
plot.simple_di_det_ipm <- function(x = NULL, y = NULL,
ipm = NULL,
col = rainbow(100, start=0.67, end=0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
exponent = 1,
n_row = 1,
n_col = 1,
...) {
dots <- list(...)
old_par <- graphics::par(no.readonly = TRUE)
on.exit(par(old_par))
# This is used so that users can just say plot(my_model) instead of
# plot(ipm = my_model). ipmr_matrix expects x and y to both be NULL
if(!is.null(x) && is.null(ipm)){
ipm <- x
x <- NULL
}
plot_list <- ipm$sub_kernels
par(mar = c(6, 2, 3, 1),
mfrow = c(n_row, n_col))
if(do_legend) {
graphics::layout(mat = cbind(matrix(1, 5, 5), rep(2, 5)))
}
for(i in seq_along(plot_list)){
use_kern <- plot_list[[i]]
nm <- names(plot_list)[i]
plot.ipmr_matrix(x = x,
y = y,
A = use_kern ^ exponent,
col = col,
bw = bw,
do_contour = do_contour,
do_legend = do_legend,
dots,
main = nm)
}
invisible(ipm)
}
#' @rdname plot_star
#' @export
plot.simple_di_stoch_param_ipm <- function(x = NULL, y = NULL,
ipm = NULL,
col = rainbow(100, start=0.67, end=0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
exponent = 1,
n_row = 1,
n_col = 1,
...) {
dots <- list(...)
old_par <- graphics::par(no.readonly = TRUE)
on.exit(par(old_par))
# This is used so that users can just say plot(my_model) instead of
# plot(ipm = my_model). ipmr_matrix expects x and y to both be NULL
if(!is.null(x) && is.null(ipm)){
ipm <- x
x <- NULL
}
plot_list <- ipm$sub_kernels
if(all(is.na(plot_list[[1]]))) {
stop("Cannot plot sub-kernels unless 'make_ipm(return_sub_kernels = TRUE)'!")
}
plt_seq <- seq_along(plot_list)
par(mar = c(6, 2, 3, 1),
mfrow = c(n_row, n_col))
if(do_legend) {
graphics::layout(mat = cbind(matrix(1, 5, 5), rep(2, 5)))
}
for(i in seq_along(plot_list)){
use_kern <- plot_list[[i]]
nm <- names(plot_list)[i]
plot.ipmr_matrix(x = x,
y = y,
A = use_kern ^ exponent,
col = col,
bw = bw,
do_contour = do_contour,
do_legend = do_legend,
dots,
main = nm)
}
invisible(ipm)
}
#' @rdname plot_star
#' @export
plot.simple_di_stoch_kern_ipm <- function(x = NULL, y = NULL,
ipm = NULL,
col = rainbow(100, start=0.67, end=0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
exponent = 1,
n_row = 1,
n_col = 1,
...) {
dots <- list(...)
old_par <- graphics::par(no.readonly = TRUE)
on.exit(par(old_par))
# This is used so that users can just say plot(my_model) instead of
# plot(ipm = my_model). ipmr_matrix expects x and y to both be NULL
if(!is.null(x) && is.null(ipm)){
ipm <- x
x <- NULL
}
plot_list <- ipm$sub_kernels
par(mar = c(6, 2, 3, 1),
mfrow = c(n_row, n_col))
if(do_legend) {
graphics::layout(mat = cbind(matrix(1, 5, 5), rep(2, 5)))
}
for(i in seq_along(plot_list)){
use_kern <- plot_list[[i]]
nm <- names(plot_list)[i]
plot.ipmr_matrix(x = x,
y = y,
A = use_kern ^ exponent,
col = col,
bw = bw,
do_contour = do_contour,
do_legend = do_legend,
dots,
main = nm)
}
invisible(ipm)
}
#' @rdname plot_star
#' @inheritParams format_mega_kernel
#'
#' @export
plot.general_di_det_ipm <- function(x = NULL, y = NULL,
ipm = NULL,
mega_mat = NA_character_,
col = rainbow(100,
start = 0.67,
end = 0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
exponent = 1,
n_row = 1,
n_col = 1,
...) {
dots <- list(...)
old_par <- graphics::par(no.readonly = TRUE)
on.exit(par(old_par))
# This is used so that users can just say plot(my_model) instead of
# plot(ipm = my_model). ipmr_matrix expects x and y to both be NULL
if(!is.null(x) && is.null(ipm)){
ipm <- x
x <- NULL
}
mega_mat <- rlang::enquo(mega_mat)
if(rlang::quo_text(mega_mat) == "NA") {
stop("Plotting general IPMs requires building a 'mega_mat'.\n",
"Please specify an expression for the 'mega_mat' argument.")
}
graphics::par(mar = c(6, 5, 3, 2),
mfrow = c(n_row, n_col))
plot_list <- format_mega_kernel(ipm, mega_mat = !! mega_mat)
if(do_legend) {
graphics::layout(mat = cbind(matrix(1, 5, 5), rep(2, 5)))
}
lapply(plot_list, function(ipm) plot.ipmr_matrix(x = x,
y = y,
A = ipm ^ exponent,
col = col,
bw = bw,
do_contour = do_contour,
do_legend = do_legend,
dots))
invisible(ipm)
}
.ncol_nrow <- function(plt_seq) {
out <- list(ncol = NA,
nrow = NA)
if(length(plt_seq) > 25) {
out$ncol <- 5
out$nrow <- 5
} else {
out$ncol <- out$nrow <- ceiling(sqrt(length(plt_seq)))
}
return(out)
}
# right_ev ----------------
#' @rdname eigenvectors
#'
#' @title Compute the standardized left and right eigenvectors via iteration
#'
#' @param ipm Output from \code{make_ipm()}.
#' @param ... Other arguments passed to methods
#'
#' @return A list of named numeric vector(s) corresponding to the stable trait distribution
#' function (\code{right_ev}) or the reproductive values for each trait (\code{left_ev}).
#'
#' @export
right_ev <- function(ipm, ... ) {
UseMethod('right_ev')
}
#' @rdname eigenvectors
#' @param iterations The number of times to iterate the model to reach
#' convergence. Default is 100.
#' @param tolerance Tolerance to evaluate convergence to asymptotic dynamics.
#'
#' @export
right_ev.simple_di_det_ipm <- function(ipm,
iterations = 100,
tolerance = 1e-10,
...) {
mod_nm <- deparse(substitute(ipm))
# if it's already been iterated to convergence, we don't have much work to do.
if(.already_iterated(ipm)) {
# get index for population vector of final iteration
final_it <- dim(ipm$pop_state[[1]])[2]
if(all(is_conv_to_asymptotic(ipm,
tolerance = tolerance))) {
out <- .extract_conv_ev_simple(ipm$pop_state)
} else {
# If it's iterated, but not converged, we can at least use the final
# iteration which is hopefully closer to convergence than the user-defined
# initial population vector
message(
paste(
"'",
mod_nm,
"'",
' did not converge to asymptotic dynamics after ',
final_it,
' iterations.\n',
'Will re-iterate the model ',
iterations,
' times and check for convergence.',
sep = ""
)
)
init_pop_vec <- lapply(ipm$pop_state[!grepl("lambda", names(ipm$pop_state))],
function(x, final_it) {
x[ , final_it]
},
final_it = ncol(ipm$pop_state[[1]]))
#
# pop_nm <- ipm$proto_ipm$state_var %>%
# unlist() %>%
# unique() %>%
# .[1] %>%
# paste('n_', ., sep = "")
#
# pop_nm <- rlang::ensym(pop_nm)
test_conv <- ipm$proto_ipm %>%
define_pop_state(!!! init_pop_vec) %>%
make_ipm(iterate = TRUE,
iterations = iterations)
if(all(is_conv_to_asymptotic(test_conv,
tolerance = tolerance))) {
out <- .extract_conv_ev_simple(test_conv$pop_state)
message('model is now converged :)')
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
final_it + iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
}
} else {
message(
paste(
"'",
mod_nm,
"'",
' has not been iterated yet. ',
'Generating a population vector using runif() and\niterating the model ',
iterations,
' times to check for convergence to asymptotic dynamics',
sep = ""
)
)
# A model that hasn't been iterated yet
# Create variables for internal usage
pop_nm <- .get_pop_nm_simple(ipm)
pop_states <- paste('n', pop_nm, sep = '_')
# Final step is to generate an initial pop_state. This is always just
# vector drawn from a random uniform distribution
len_pop_state <- dim(ipm$sub_kernels[[1]])[1]
init_pop <- stats::runif(len_pop_state)
# Drop _t for defining the initial population vector so define_pop_state
# doesn't complain
test_conv <- ipm$proto_ipm %>%
define_pop_state(!! pop_states[1] := init_pop) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
normalize_pop_size = TRUE)
if(all(is_conv_to_asymptotic(test_conv,
tolerance = tolerance))) {
out <- .extract_conv_ev_simple(test_conv$pop_state)
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
}
# Stuff into a list and standardize
# out <- rlang::list2(!! out_nm := (out / sum(out)))
out <- Filter(function(x) !any(is.na(x)), out)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'w', sep = '_')
class(out) <- "ipmr_w"
return(out)
}
#' @rdname eigenvectors
#' @param burn_in The proportion of early iterations to discard from the
#' stochastic simulation
#' @export
right_ev.simple_di_stoch_kern_ipm <- function(ipm,
burn_in = 0.25,
...) {
mod_nm <- deparse(substitute(ipm))
# Identify state variable name
pop_nm <- .get_pop_nm_simple(ipm)
out <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq],
burn_seq = burn_in_seq)
# Normalize just in case. easy for simple IPMs where we know there's only
# one trait contributing to total population size. Tricker in general IPMs.
out[[1]] <- apply(out[[1]], 2, function(x) x / sum(x))
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'w', sep = '_')
class(out) <- "ipmr_w"
return(out)
}
#' @rdname eigenvectors
#' @export
right_ev.simple_di_stoch_param_ipm <- right_ev.simple_di_stoch_kern_ipm
#' @rdname eigenvectors
#' @export
right_ev.general_di_det_ipm <- function(ipm,
iterations = 100,
tolerance = 1e-10,
...) {
mod_nm <- deparse(substitute(ipm))
final_it <- dim(ipm$pop_state[[1]])[2]
if(all(is_conv_to_asymptotic(ipm,
tolerance = tolerance))) {
out <- .extract_conv_ev_general(ipm$pop_state, ipm$proto_ipm)
} else {
# Not yet converged, so re-iterate the model and see what we get
message(
paste(
"'",
mod_nm,
"'",
' did not converge to asymptotic dynamics after ',
final_it,
' iterations.\n',
'Will re-iterate the model ',
iterations,
' times and check for convergence.',
sep = ""
)
)
use_pop_state <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
init_pop_vec <- lapply(use_pop_state,
function(x, final_it) x[ , final_it],
final_it = final_it)
# Remove NAs generated by the unsubstituted pop_state place holders
init_pop_vec <- Filter(f = function(y) !any(is.na(y)), x = init_pop_vec)
names(init_pop_vec) <- gsub('pop_state', 'n', names(init_pop_vec))
test_conv <- ipm$proto_ipm %>%
define_pop_state(
pop_vectors = init_pop_vec
) %>%
make_ipm(iterate = TRUE,
iterations = iterations)
if(all(is_conv_to_asymptotic(test_conv,
tolerance = tolerance))) {
out <- .extract_conv_ev_general(test_conv$pop_state, test_conv$proto_ipm)
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
}
out <- Filter(function(x) !any(is.na(x)), out)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'w', sep = '_')
class(out) <- "ipmr_w"
return(out)
}
#' @rdname eigenvectors
#' @export
right_ev.general_di_stoch_kern_ipm <- function(ipm,
burn_in = 0.25,
...) {
mod_nm <- deparse(substitute(ipm))
out <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq, drop = FALSE],
burn_seq = burn_in_seq)
# Normalize everything just in case
pop_sizes <- lapply(out, colSums) %>%
do.call(what = ".add", args = .)
for(i in seq_along(out)) {
for(j in seq_len(ncol(out[[1]]))) {
out[[i]][ , j] <- out[[i]][ ,j] / pop_sizes[j]
}
}
out <- Filter(function(x) !any(is.na(x)), out)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'w', sep = '_')
class(out) <- "ipmr_w"
return(out)
}
#' @rdname eigenvectors
#' @export
right_ev.general_di_stoch_param_ipm <- right_ev.general_di_stoch_kern_ipm
# left_ev -----------------
#' @rdname eigenvectors
#' @export
left_ev <- function(ipm, ...) {
UseMethod('left_ev')
}
#' @export
#' @rdname eigenvectors
#' @importFrom stats runif
left_ev.simple_di_det_ipm <- function(ipm,
iterations = 100,
tolerance = 1e-10,
...) {
mod_nm <- deparse(substitute(ipm))
# Identify state variable name
pop_nm <- .get_pop_nm_simple(ipm)
if(any(ipm$proto_ipm$uses_par_sets)) {
ps_nms <- .flatten_to_depth(ipm$proto_ipm$par_set_indices, 1L) %>%
names() %>%
unique()
ps_nms <- ps_nms[ps_nms != "levels"]
ps_suff <- paste(ps_nms, collapse = "_")
pop_nm <- paste(pop_nm, ps_suff, sep = "_")
}
# Create variables for internal usage
pop_states <- paste('n', pop_nm, sep = '_')
# Final step is to generate an initial pop_state. This is always just
# vector drawn from a random uniform distribution
len_pop_state <- dim(ipm$sub_kernels[[1]])[1]
init_pop <- stats::runif(len_pop_state)
test_conv <- ipm$proto_ipm %>%
define_pop_state(!! pop_states[1] := init_pop) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
iteration_direction = "left",
normalize_pop_size = TRUE)
if(all(is_conv_to_asymptotic(test_conv,
tolerance = tolerance))) {
out <- .extract_conv_ev_simple(test_conv$pop_state)
# out_nm <- paste(pop_nm, 'v', sep = "_")
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
# Stuff into a list and standardize
out <- Filter(function(x) !any(is.na(x)), out)
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'v', sep = '_')
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#' @param kernel_seq The sequece of parameter set indices used to select kernels
#' during the iteration procedure. If \code{NULL}, will use the sequence stored
#' in the \code{ipm} object. Should usually be left as \code{NULL}.
#' @export
left_ev.simple_di_stoch_kern_ipm <- function(ipm,
iterations = 10000,
burn_in = 0.25,
kernel_seq = NULL,
...) {
mod_nm <- deparse(substitute(ipm))
# Identify state variable name
pop_nm <- .get_pop_nm_simple(ipm)
# Create variables for internal usage
pop_states <- paste('n', pop_nm, c('t', 't_1'), sep = '_')
# Final step is to generate an initial pop_state. This is always just
# vector drawn from a random uniform distribution
len_pop_state <- nrow(ipm$sub_kernels[[1]])
init_pop <- stats::runif(len_pop_state)
if(is.null(kernel_seq)) kernel_seq <- ipm$env_seq
# Drop _t for defining the initial population vector so define_pop_state
# doesn't complain
pop_states[1] <- gsub("_t$", "", pop_states[1])
test_conv <- ipm$proto_ipm %>%
define_pop_state(!! pop_states[1] := init_pop) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
iteration_direction = "left",
normalize_pop_size = TRUE,
kernel_seq = kernel_seq)
out <- test_conv$pop_state
out$lambda <- NULL
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq],
burn_seq = burn_in_seq)
out_nm <- paste(pop_nm, 'v', sep = "_")
names(out) <- out_nm
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#'
#' @section \strong{Deterministic eigenvectors}:
#' For \code{right_ev}, if the model has already been iterated and has
#' converged to asymptotic dynamics, then it will just extract the final
#' population state and return that in a named list. Each element of the list
#' is a vector with length \code{>= 1} and corresponds each state variable's
#' portion of the eigenvector.
#' If the model has been iterated, but has not yet converged to asymptotic dynamics,
#' \code{right_ev} will try to iterate it further using the final population state
#' as the starting point. The default number of iterations is 100, and can be
#' adjusted using the \code{iterations} argument.
#' If the model hasn't been iterated, then \code{right_ev} will try iterating it
#' for \code{iterations} number of time steps and check for convergence. In the
#' latter two cases, if the model still has not converged to asymptotic dynamics,
#' it will return \code{NA} with a warning.
#'
#' For \code{left_ev}, the transpose iteration (\emph{sensu} Ellner & Rees 2006,
#' Appendix A) is worked out based on the \code{state_start} and \code{state_end}
#' in the model's \code{proto_ipm} object. The model is then iterated for
#' \code{iterations} times to produce a standardized left eigenvector.
#'
#' @section \strong{Stochastic eigenvectors}:
#' \code{left_ev} and \code{right_ev} return different things for stochastic models.
#' \code{right_ev} returns the trait distribution through time from the stochastic
#' simulation (i.e. \code{ipm$pop_state}), and normalizes it such that the
#' distribution at each time step integrates to 1 (if it is not already).
#' It then discards the first \code{burn_in * iterations} time steps of the
#' simulation to eliminate transient dynamics. See Ellner, Childs, & Rees 2016,
#' Chapter 7.5 for more details.
#'
#' \code{left_ev} returns a similar result as \code{right_ev}, except the trait
#' distributions are the result of left multiplying the kernel and trait
#' distribution. See Ellner, Childs, & Rees 2016, Chapter 7.5 for more
#' details.
#'
#' @export
left_ev.general_di_det_ipm <- function(ipm,
iterations = 100,
tolerance = 1e-10,
...) {
mod_nm <- deparse(substitute(ipm))
use_pop_state <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
init_pop_vec <- lapply(use_pop_state, function(x) x[ , 1])
# Remove NAs generated by the unsubstituted pop_state place holders
init_pop_vec <- Filter(f = function(y) !any(is.na(y)), x = init_pop_vec)
names(init_pop_vec) <- gsub('pop_state', 'n', names(init_pop_vec))
test_conv <- ipm$proto_ipm %>%
define_pop_state(
pop_vectors = init_pop_vec
) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
iteration_direction = "left",
normalize_pop_size = TRUE)
if(all(is_conv_to_asymptotic(test_conv, tolerance = tolerance))) {
out <- .extract_conv_ev_general(test_conv$pop_state, test_conv$proto_ipm)
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
out <- Filter(function(x) !any(is.na(x)), out)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'v', sep = '_')
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#' @export
left_ev.general_di_stoch_kern_ipm <- function(ipm,
iterations = 10000,
burn_in = 0.25,
kernel_seq = NULL,
...) {
mod_nm <- deparse(substitute(ipm))
use_pop_state <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
init_pop_vec <- lapply(use_pop_state, function(x) x[ , 1])
names(init_pop_vec) <- gsub('pop_state', 'n', names(init_pop_vec))
if(is.null(kernel_seq)) kernel_seq <- ipm$env_seq
test_conv <- ipm$proto_ipm %>%
define_pop_state(
pop_vectors = init_pop_vec
) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
iteration_direction = "left",
normalize_pop_size = TRUE,
kernel_seq = kernel_seq)
out <- test_conv$pop_state[!grepl('lambda', names(test_conv$pop_state))]
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq, drop = FALSE],
burn_seq = burn_in_seq)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'v', sep = '_')
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#' @export
# Needs to re-iterate the model, but use the same kernel sequence that is used
# in right_ev/make_ipm. In effect, we want a deterministic model using the
# sub-kernels slot in lieu of parameters.
left_ev.general_di_stoch_param_ipm <- function(ipm,
iterations = 10000,
burn_in = 0.25,
kernel_seq = NULL,
...) {
mod_nm <- deparse(substitute(ipm))
sub_kernels <- ipm$sub_kernels
proto_ipm <- ipm$proto_ipm
usr_its <- ncol(ipm$pop_state[[1]]) - 1
if(usr_its < iterations) {
stop("'iterations' should be less than or equal to the number of iterations in 'make_ipm()'!")
}
# Set up the model iteration expressions and initial population state.
# This is a leaner version of make_ipm.general_di_stoch_param, because
# we don't need all the user functions or other stuff - just the basic
# infrastructure to house model iteration evaluation.
proto_list <- .initialize_kernels(proto_ipm,
iterate = TRUE,
iter_dir = "left")
others <- proto_list$others
k_row <- proto_list$k_row
main_env <- .make_main_env(others$domain,
proto_ipm$usr_funs[[1]],
age_size = .uses_age(others))
temp <- .prep_di_output(others, k_row, proto_ipm,
iterations, normal= TRUE)
main_env <- .add_pop_state_to_main_env(temp$pop_state, main_env)
# Throw this in just in case we have parameter set indices in addition
# to continuous environmental variation.
kern_seq <- .make_kern_seq(others, iterations, kernel_seq)
# Create an index to select sub-kernels based on which iteration they were
# generated in.
kern_it_ind <- vapply(names(sub_kernels), function(x) {
as.integer(strsplit(x, "(_it_)")[[1]][2])
}, integer(1L))
for(i in seq_len(iterations)) {
# This selects sub-kernels by matching the _it_XXX suffix to the current
# iteration.
use_kerns <- sub_kernels[kern_it_ind == i]
# Now, drop the _it_XXX suffix and pretend they're just regular kernels
# generated by .make_sub_kernel_general_lazy.
rm_regex <- paste("_it_", i, '$', sep = "")
names(use_kerns) <- gsub(rm_regex, "", names(use_kerns))
# We don't care about any of the other output from model iteration this,
# so skip all the naming/stashing other outputs.
pop_state <- .iterate_model(proto_ipm = proto_ipm,
k_row = k_row,
sub_kern_list = use_kerns,
current_iteration = i,
kern_seq = kern_seq,
pop_state = temp$pop_state,
main_env = main_env,
normal = TRUE)
temp$pop_state <- pop_state
}
out <- temp$pop_state
out[grepl("lambda", names(out))] <- NULL
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq, drop = FALSE],
burn_seq = burn_in_seq)
names(out) <- gsub("pop_state_", "", names(out))
names(out) <- paste(names(out), 'v', sep = '_')
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#' @export
left_ev.simple_di_stoch_param_ipm <- left_ev.general_di_stoch_param_ipm
#' @noRd
# Helper to be used in do.call(".add", X) where X is a list with potentially more
# than two elements. "+" won't work in this case, and sum() will return a single
# number.
.add <- function(...) {
dots <- list(...)
Reduce("+", x = dots, init = 0)
}
#' @rdname as_matrix
#' @title Convert to bare matrices
#' @description Converts objects to \code{c("matrix", "array")}.
#' @param x An object of class \code{ipmr_matrix}, or the output from
#' \code{make_ipm}.
#' @param ... ignored.
#'
#' @return A matrix.
#' @export
as.matrix.ipmr_matrix <- function(x, ...) {
class(x) <- c("matrix", "array")
return(x)
}
#' @rdname as_matrix
#' @export
as.matrix.ipmr_ipm <- function(x, ...) {
lapply(x$sub_kernels, as.matrix)
}
Try the ipmr package in your browser
Any scripts or data that you put into this service are public.
ipmr documentation built on Feb. 16, 2023, 10:04 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 .all_params_in_mod print.proto_ipm
Documented in print.proto_ipm
# Print----------------
#' @title Print proto_ipms or *_ipm objects
#' @rdname print_star
#'
#' @param x An object of class \code{proto_ipm} or produced by \code{make_ipm()}.
#' @param ... Ignored
#'
#' @details For printing \code{proto_ipm} objects, indices are wrapped in
#' \code{<index>} to assist with debugging. These are not carried into the model,
#' just a visual aid.
#'
#' @export
print.proto_ipm <- function(x, ...) {
n_kerns <- dim(x)[1]
cls_switch <- class(x)[1]
pretty_class <- .pretty_class(cls_switch)
msg <- paste(ifelse(.uses_age(x),
"An age structured",
"A"),
pretty_class,
"proto_ipm with",
n_kerns,
"kernels defined:\n")
msg <- c(msg, paste(x$kernel_id, collapse = ', '))
msg <- c(msg, '\n\nKernel formulae:\n\n')
cat(msg)
forms <- kernel_formulae(x)
print(forms)
cat("\n\nVital rates:\n\n")
vrs <- vital_rate_exprs(x)
print(vrs)
cat("\n\nParameter names:\n\n")
pars <- parameters(x)
print(names(pars))
# Check if all parameters are in vital rate expressions, and if all
# vital rate expression args are found in parameters.
in_mod <- .all_params_in_mod(x)
in_mod_tst <- vapply(in_mod, function(x) x$out, logical(1L)) %>%
all()
cat("\nAll parameters in vital rate expressions found in 'data_list': ",
in_mod_tst)
if(!in_mod_tst) {
stoch_params <- c(paste("simple_d", c('i', 'd'), "_stoch_param", sep = ""),
paste("general_d", c('i', 'd'), "_stoch_param", sep = ""))
if(inherits(x, stoch_params)) {
cat("\nDid not test parameters in 'define_env_state()'.\n")
}
cat("\nItems in vital rate expressions not found in 'data_list':\n")
for(i in seq_along(in_mod)) {
if(!in_mod[[i]]$out) {
nm <- paste("\n*", names(in_mod)[i], "*\n", sep = "")
miss_msg <- paste(nm,
paste(in_mod[[i]]$missing, collapse = "\n"),
sep = "")
cat(miss_msg, '\n')
}
}
cat("\nThis may be because 'define_domains()' has not yet been called.\n",
"If it has been called, then double check the model definition and",
"'data_list'.\n")
}
cat("\n\nDomains for state variables:\n\n")
print(domains(x))
cat("\n\nPopulation states defined:\n\n")
print(pop_state(x))
cat("\n\nInternally generated model iteration procedure:\n\n")
if(!all(is.na(x$domain))){
k_row <- .init_iteration(x, TRUE, direction = "right")
print(kernel_formulae(k_row))
} else {
cat("Not enough information to generate an iteration procedure.\n",
"define_impl() and define_pop_state() must be called first to generate one.")
}
# Add more later -----------
cat("\n\n")
invisible(x)
}
#' @noRd
.all_params_in_mod <- function(proto) {
kern_list <- list()
for(row in seq_len(nrow(proto))) {
use_proto <- proto[row, ]
kern_nm <- proto$kernel_id[row]
vr_exprs <- vital_rate_exprs(use_proto)
params <- parameters(proto)
# Age x size models and/or par_sets models need expansion before we
# we can test whether all parameters are in the model
if(.uses_age(proto) | any(proto$uses_par_sets)) {
all_levs <- .flatten_to_depth(proto$par_set_indices[proto$uses_par_sets],
1L) %>%
.[!duplicated(names(.))]
par_set_nms <- names(all_levs)
if(length(all_levs) > 1) {
lev_df <- expand.grid(all_levs,
stringsAsFactors = FALSE)
} else {
lev_df <- as.data.frame(all_levs)
}
# Hold the results
all_vr_exprs <- list()
it <- 1
# Replace suffixes in expressions
for(i in seq_along(vr_exprs)) {
base_expr <- rlang::expr_text(vr_exprs[[i]])
for(j in seq_len(dim(lev_df)[2])) {
par_set_nm <- names(lev_df)[j]
for(k in seq_len(dim(lev_df)[1])) {
# Replace suffixes in expression and name
all_vr_exprs[[it]] <- gsub(par_set_nm,
lev_df[k, j],
base_expr)
names(all_vr_exprs)[it] <- gsub(par_set_nm,
lev_df[k, j],
names(vr_exprs)[i])
it <- it + 1
}
}
}
vr_exprs <- lapply(all_vr_exprs, rlang::parse_expr)
vr_args <- .vr_args(vr_exprs)
} else {
vr_args <- .vr_args(vr_exprs)
}
# If a vital rate expression creates another parameter, then remove that
# from this test, as it shouldn't appear in the data_list anyway.
to_exclude <- names(vr_exprs)
domains <- names(domains(proto)) %>%
vapply(function(x) paste(x, c("_1", "_2"), sep = ""),
character(2L)) %>%
as.vector()
vr_args <- Filter(.can_be_number, vr_args)
out <- vapply(vr_args[!vr_args %in% to_exclude],
function(nm, pars) nm %in% pars,
logical(1L),
pars = c(names(params), domains))
if(!all(out)) {
vr_args <- vr_args[!vr_args %in% to_exclude]
ind <- which(!vr_args %in% c(names(params), domains))
miss_args <- vr_args[ind]
out <- FALSE
} else {
out <- TRUE
miss_args <- NA_character_
}
kern_list <- c(kern_list,
list(list(out = out,
missing = miss_args)))
names(kern_list)[row] <- kern_nm
}
return(kern_list)
}
.can_be_number <- function(x) {
suppressWarnings(is.na(as.numeric(x)))
}
#' @noRd
.vr_args <- function(vr_exprs) {
all_args <- lapply(vr_exprs, function(y) {
temp <- rlang::expr_text(y)
.args_from_txt(temp)
}) %>%
unlist() %>%
unique()
return(all_args)
}
#' @noRd
.uses_age <- function(x) {
inherits(x, "age_x_size")
}
#' @noRd
.pretty_class <- function(cls_switch) {
out <-
switch(
cls_switch,
# IPM classes -----------------------------
'simple_di_det_ipm' = "simple, density independent, deterministic",
'simple_di_stoch_kern_ipm' = "simple, density independent, stochastic, kernel-resampled",
'simple_di_stoch_param_ipm' = "simple, density independent, stochastic, parameter-resampled",
'simple_dd_det_ipm' = "simple, density dependent, deterministic",
'simple_dd_stoch_kern_ipm' = "simple, density dependent, stochastic, kernel-resampled",
'simple_dd_stoch_param_ipm' = "simple, density dependent, stochastic, parameter-resampled",
'general_di_det_ipm' = "general, density independent, deterministic",
'general_di_stoch_kern_ipm' = "general, density independent, stochastic, kernel-resampled",
'general_di_stoch_param_ipm' = "general, density independent, stochastic, parameter-resampled",
'general_dd_det_ipm' = "general, density dependent, deterministic",
'general_dd_stoch_kern_ipm' = "general, density dependent, stochastic, kernel-resampled",
'general_dd_stoch_param_ipm' = "general, density dependent, stochastic, parameter-resampled",
# proto_ipm classes -----------------------
'simple_di_det' = "simple, density independent, deterministic",
'simple_di_stoch_kern' = "simple, density independent, stochastic, kernel-resampled",
'simple_di_stoch_param' = "simple, density independent, stochastic, parameter-resampled",
'simple_dd_det' = "simple, density dependent, deterministic",
'simple_dd_stoch_kern' = "simple, density dependent, stochastic, kernel-resampled",
'simple_dd_stoch_param' = "simple, density dependent, stochastic, parameter-resampled",
'general_di_det' = "general, density independent, deterministic",
'general_di_stoch_kern' = "general, density independent, stochastic, kernel-resampled",
'general_di_stoch_param' = "general, density independent, stochastic, parameter-resampled",
'general_dd_det' = "general, density dependent, deterministic",
'general_dd_stoch_kern' = "general, density dependent, stochastic, kernel-resampled",
'general_dd_stoch_param' = "general, density dependent, stochastic, parameter-resampled"
)
return(out)
}
#' @rdname print_star
#' @title Generics for IPM classes
#'
#' @param comp_lambda A logical indicating whether or not to calculate lambdas
#' for the iteration kernels and display them.
#' @param sig_digits The number of significant digits to round to if \code{
#' comp_lambda = TRUE}.
#' @param type_lambda Either \code{'all'} or \code{'stochastic'}. See
#' \code{\link{lambda}} for more details.
#' @param check_conv A logical: for deterministic models, check if population state
#' has converged to asymptotic dynamics? If \code{TRUE} and the model has not
#' converged, a message will be printed.
#' @param ... Ignored
#'
#' @return \code{x} invisibly.
#'
#' @export
print.simple_di_det_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'last',
sig_digits = 3,
check_conv = TRUE,
...) {
pretty_cls <- .pretty_class(class(x)[1])
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) defined.', sep = "")
if(comp_lambda) {
lambdas <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nDeterministic lambda = ',
round(lambdas, sig_digits),
sep = "")
msg <- c(msg, l_msg)
if(check_conv &&
!all(is_conv_to_asymptotic(x))) {
# Captures the name of the model that the user gave rather than
# just print "x isn't converged"
mod_nm <- deparse(substitute(x))
message(
paste(mod_nm,
' has has not converged to asymptotic dynamics!',
sep = "")
)
}
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.simple_dd_det_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'last',
sig_digits = 3,
...) {
mod_nm <- deparse(substitute(x))
pretty_cls <- .pretty_class(class(x)[1])
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) defined.', sep = "")
if(comp_lambda) {
lambdas <- lambda(x, type_lambda = type_lambda)
l_msg <- paste('\nLambda for the final time step of the model is: ',
round(lambdas, sig_digits),
sep = "")
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.simple_di_stoch_kern_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'stochastic',
sig_digits = 3,
...) {
mod_nm <- deparse(substitute(x))
pretty_cls <- .pretty_class(class(x)[1])
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$iterators),
' iteration kernel(s) and ',
length(x$sub_kernels),
' sub-kernel(s) defined.', sep = "")
if(comp_lambda) {
all_lams <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nStochastic lambda for ',
mod_nm,
' = ',
round(all_lams, sig_digits),
sep = "" )
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.simple_dd_stoch_kern_ipm <- print.simple_di_stoch_kern_ipm
#' @rdname print_star
#' @export
print.simple_di_stoch_param_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'stochastic',
sig_digits = 3,
...) {
pretty_cls <- .pretty_class(class(x)[1])
mod_nm <- deparse(substitute(x))
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$iterators),
' iteration kernel(s) and ',
length(x$sub_kernels),
' sub-kernel(s) defined.', sep = "")
if(comp_lambda) {
all_lams <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nStochastic lambda for ',
mod_nm,
' = ',
round(all_lams, sig_digits),
sep = "" )
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.simple_dd_stoch_param_ipm <- print.simple_di_stoch_param_ipm
#' @rdname print_star
#' @export
print.general_di_det_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'last',
sig_digits = 3,
check_conv = TRUE,
...) {
pretty_cls <- .pretty_class(class(x)[1])
pops <- x$pop_state[!grepl("lambda", names(x$pop_state))]
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) and ',
length(pops),
' population vectors defined.', sep = "")
mod_nm <- deparse(substitute(x))
if(comp_lambda) {
ret_lam <- lambda(x, type_lambda = type_lambda)
l_msg <- paste('\nLambda for the final time step of the model is: ',
round(ret_lam, sig_digits),
'\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
}
if(check_conv && !all(is_conv_to_asymptotic(x))) {
message(
paste(mod_nm,
' has has not converged to asymptotic dynamics!',
sep = "")
)
}
msg <- c(msg, l_msg)
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.general_dd_det_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'last',
sig_digits = 3,
...) {
pretty_cls <- .pretty_class(class(x)[1])
pops <- x$pop_state[!grepl("lambda", names(x$pop_state))]
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) and ',
length(pops),
' population vectors defined.', sep = "")
mod_nm <- deparse(substitute(x))
if(comp_lambda) {
ret_lam <- lambda(x, type_lambda = type_lambda)
l_msg <- paste('\nLambda for the final time step of the model is: ',
round(ret_lam, sig_digits),
'\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg)
}
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.general_di_stoch_kern_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'stochastic',
sig_digits = 3,
...) {
pretty_cls <- .pretty_class(class(x)[1])
mod_nm <- deparse(substitute(x))
msg <- paste0('A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s) and ',
length(x$pop_state[!grepl("lambda", names(x$pop_state))]),
' population vectors defined.',
sep = "")
if(comp_lambda) {
all_lams <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nStochastic lambda for ',
mod_nm,
' = ',
round(all_lams, sig_digits),
sep = "" )
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
# No check for convergence here - there may quite reasonably be no
# convergence for
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.general_dd_stoch_kern_ipm <- print.general_di_stoch_kern_ipm
#' @rdname print_star
#' @export
print.general_di_stoch_param_ipm <- function(x,
comp_lambda = TRUE,
type_lambda = 'stochastic',
sig_digits = 3,
...) {
pretty_cls <- .pretty_class(class(x)[1])
mod_nm <- deparse(substitute(x))
msg <- paste0(
'A ',
pretty_cls,
' IPM with ',
length(x$sub_kernels),
' sub-kernel(s), ',
length(x$pop_state[!grepl("lambda", names(x$pop_state))]),
' population vectors, and ',
dim(x$env_seq)[2],
' environmental parameters defined.',
sep = ""
)
env_nms <- vapply(dimnames(x$env_seq)[[2]],
FUN = function(x) strsplit(x, '\\.')[[1]][2],
FUN.VALUE = character(1L))
msg <- c(msg,
paste(
'\nEnvironmental variables are: ',
paste(
paste(env_nms, collapse = ', ')
)
))
if(comp_lambda) {
all_lams <- lambda(x, type_lambda = type_lambda)
l_msg <- paste0('\nStochastic lambda for ',
mod_nm,
' = ',
round(all_lams, sig_digits),
sep = "" )
det_lam_msg <- paste('\nCall lambda(',
mod_nm,
', type_lambda = "all") for deterministic lambdas\n',
'from each iteration.',
sep = "")
msg <- c(msg, l_msg, det_lam_msg)
}
# No check for convergence here - there may quite reasonably be no
# convergence for
cat(msg)
invisible(x)
}
#' @rdname print_star
#' @export
print.general_dd_stoch_param_ipm <- print.general_di_stoch_param_ipm
#' @export
print.ipmr_vital_rate_exprs <- function(x, ...) {
proto <- attr(x, "proto")
if(any(proto$uses_par_sets) | .uses_age(proto)) {
out <- .pretty_print_par_sets(x, proto)
} else {
out <- lapply(x, rlang::expr_text)
}
attr(out, "proto") <- NULL
out <- paste(names(out), out, sep = ": ") %>%
paste(., collapse = "\n")
class(out) <- "character"
cat(out)
invisible(x)
}
.pretty_print_par_sets <- function(x, proto) {
nm_par_sets <- proto$par_set_indices %>%
.flatten_to_depth(1L) %>%
names() %>%
unique()
if(.uses_age(proto)) {
ages <- proto$age_indices %>%
.flatten_to_depth(1L) %>%
names() %>%
unique()
nm_par_sets <- c(nm_par_sets, ages)
}
nm_par_sets <- nm_par_sets[!is.na(nm_par_sets)]
x <- lapply(x, rlang::expr_text)
for(i in seq_along(nm_par_sets)) {
to_print <- paste("<", nm_par_sets[i], ">", sep = "")
x <- lapply(x,
function(y, nm, to_print) {
gsub(nm, to_print, y)
},
nm = nm_par_sets[i],
to_print = to_print)
names(x) <- gsub(nm_par_sets[i], to_print, names(x))
}
return(x)
}
#' @export
print.ipmr_kernel_exprs <- function(x, ...) {
proto <- attr(x, "proto")
if(any(proto$uses_par_sets) | .uses_age(proto)) {
out <- .pretty_print_par_sets(x, proto)
} else {
out <- lapply(x, rlang::expr_text)
}
attr(out, "proto") <- NULL
out <- paste(names(out), out, sep = ": ") %>%
paste(., collapse = "\n")
class(out) <- "character"
cat(out)
invisible(x)
}
#' @export
print.ipmr_parameters <- function(x, ...) {
out <- as.numeric(x)
names(out) <- names(x)
print(out)
invisible(x)
}
#' @export
# Will need to update w/ other names for other integration rules as those
# become available
print.ipmr_domains <- function(x, ...) {
# Strips class and proto attributes
out <- lapply(x, function(y) {
y
})
out <- paste(names(out), out, sep = ": ") %>%
paste(., collapse = "\n")
out <- gsub("c\\(", "", out)
out <- gsub("\\)", "", out)
cat(out)
invisible(x)
}
#' @export
print.ipmr_pop_state <- function(x, ...) {
out <- lapply(x, function(y) y)
if(out[[1]] == "No population state defined.") {
cat(out[[1]])
} else {
if(rlang::is_named(x) && !grepl("pop_state_", names(out))) {
message("Names of population states don't appear to be prefixed with 'n_'.")
}
names(out) <- gsub("pop_state_", "n_", names(out))
out <- paste(names(out), out, sep = ": ") %>%
paste(., collapse = "\n")
cat(out)
}
invisible(x)
}
#' @export
print.ipmr_mp_mesh <- function(x, ...) {
uses <- Filter(Negate(is.null), x)
msg <- vapply(uses, function(z) {
if(length(z) > 1) {
paste("Lower bound: ", range(z)[1], ", Upper bound: ", range(z)[2],
", # of Meshpoints: ", sqrt(length(z)),
sep = "")
} else {
as.character(z)
}
},
character(1L)) %>%
paste(names(.), " - ", ., sep = "") %>%
paste(., collapse = "\n")
cat(msg)
invisible(x)
}
#' @export
print.ipmr_matrix <- function(x, ...) {
rng <- range(x)
msg <- paste("Minimum value: ",
round(rng[1], 5),
", maximum value: ",
round(rng[2], 5), sep = "")
pp_tst <- all(x >= 0)
msg <- c(msg,
paste("\nAll entries greater than or equal to 0: ",
pp_tst ,
sep = ""))
cat("\n", msg, "\n")
invisible(x)
}
#' @export
print.ipmr_vital_rate_funs <- function(x, ...) {
if(length(x) == 1 && x == "No vital rate functions specified") {
cat("No vital rates specified\n")
invisible(x)
return()
}
vr_nms <- names(x)
cls <- vapply(x,
function(x) .mat_cls(x),
character(1L))
rngs <- vapply(x,
function(y) round(range(y), 4),
numeric(2L))
out <- vapply(seq_along(cls),
function(it, vr_nm, cls, rngs) {
paste(
paste(vr_nms[it],
" (not yet discretized)",
sep = ""),
cls[it],
sep = ": "
) %>%
paste(" with minimum value: ", rngs[1, it],
" and maximum value: ", rngs[2, it],
sep = "")
},
vr_nm = vr_nms, cls = cls, rngs = rngs,
FUN.VALUE = character(1L)) %>%
paste(., collapse = "\n")
out <- paste(out, '\n', sep = "")
cat(out)
invisible(x)
}
.pretty_dim <- function(x) {
col <- ncol(x)
row <- nrow(x)
cls <- switch(class(x)[1],
"CC" = "kernel",
"DC" = "column vector",
"CD" = "row vector",
"DD" = "matrix")
paste("A", row, "x", col, cls, sep = " ")
}
.mat_cls <- function(x) {
.pretty_dim(x)
}
# Lambda------------
#' @title Compute the per-capita growth rate for an IPM object
#' @rdname lambda
#'
#' @description Compute the per-capita growth rate for a given model. Can handle
#' stochastic and deterministic models, and has the option to discard burn in for
#' stochastic models.
#'
#' @param ipm An object returned by \code{make_ipm()}.
#' @param type_lambda Either \code{'all'}, \code{'last'},
#' or \code{'stochastic'}. \code{'all'}
#' returns a vector of lambda values for each time step of the simulation (equal
#' in length to the \code{iterations} argument of \code{make_ipm()}).
#' \code{'last'} returns the lambda value for the final timestep.
#' \code{'stochastic'} returns a single value, which by default is
#' \code{mean(log(lambda(ipm, type_lambda = "all")))}, with the proportion of
#' \code{burn_in} iterations removed from the beginning of the simulation. Set
#' \code{log} to \code{FALSE} to get \code{lambda} on the linear scale for
#' stochastic models (i.e. \code{exp(mean(log(lambdas)))}).
#' @param ... other arguments passed to methods.
#' @param burn_in The proportion of iterations to discard. Default is 0.1
#' (i.e. first 10\% of iterations in the simulation).
#' @param log Return lambda on the log scale? This is \code{TRUE} by default for
#' stochastic models, and \code{FALSE} for deterministic models.
#'
#' @return When \code{type_lambda = "all"}, an array. Rows correspond to time
#' steps, and columns correspond to parameter sets (if any). For other types,
#' a numeric vector.
#'
#' @export
lambda <- function(ipm, ...) {
UseMethod('lambda')
}
#' @rdname lambda
#' @export
lambda.simple_di_det_ipm <- function(ipm,
type_lambda = 'last',
log = FALSE,
...) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
"all" = TRUE,
'last' = FALSE)
out <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(log) out <- log(out)
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_di_stoch_kern_ipm <- function(ipm,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL,
...) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
"all" = TRUE,
"stochastic" = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_di_stoch_param_ipm <- function(ipm,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL,
...) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
"all" = TRUE,
"stochastic" = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @importFrom methods hasArg
#' @export
lambda.general_di_det_ipm <- function(ipm,
type_lambda = 'last',
log = FALSE,
...) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'last' = FALSE,
'stochastic' = stop("Cannot compute stochastic lambda for deterministic IPM",
call. = FALSE))
out <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(log) out <- log(out)
return(out)
}
#' @rdname lambda
#' @export
lambda.general_di_stoch_kern_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.general_di_stoch_param_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_dd_det_ipm <- function(ipm,
type_lambda = "all",
...,
log = FALSE) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(
type_lambda,
'all' = TRUE,
'last' = FALSE,
'stochastic' = stop("Cannot compute stochastic lambda for deterministic IPM",
call. = FALSE))
out <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(log) out <- log(out)
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_dd_stoch_kern_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.simple_dd_stoch_param_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.general_dd_det_ipm <- function(ipm,
type_lambda = 'last',
...,
log = FALSE) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'last' = FALSE,
'stochastic' = stop("Cannot compute stochastic lambda for deterministic IPM",
call. = FALSE))
out <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(log) out <- log(out)
return(out)
}
#' @rdname lambda
#' @export
lambda.general_dd_stoch_kern_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
#' @rdname lambda
#' @export
lambda.general_dd_stoch_param_ipm <- function(ipm,
...,
type_lambda = 'stochastic',
burn_in = 0.1,
log = NULL) {
.check_lambda_args(ipm, type_lambda)
all_lams <- switch(type_lambda,
'all' = TRUE,
'stochastic' = TRUE,
'last' = FALSE)
temp <- .lambda_pop_size(ipm, all_lambdas = all_lams)
if(all_lams) {
burn_ind <- seq_len(round(length(temp) * burn_in))
}
if(is.null(log)) {
if(type_lambda == "stochastic") {
message("log(lambda) is returned by default for stochastic models. Set ",
"'log = FALSE' for lambda on linear scale.")
}
log <- switch(type_lambda,
'all' = FALSE,
'last' = FALSE,
'stochastic' = TRUE)
}
out <- switch(type_lambda,
'all' = temp,
'last' = temp,
'stochastic' = .thin_stoch_lambda(temp, burn_ind, log))
if(log & type_lambda %in% c("all", "last")) {
out <- log(out)
}
return(out)
}
# Plot ---------------
# Authors - whoever wrote the code from the IPM book
#' @title Plot a matrix or an *_ipm object
#' @rdname plot_star
#'
#' @param x,y Either the values of the meshpoints or \code{NULL}. If \code{NULL},
#' then a sequence is generated so that meshpoints are given sequential bin numbers.
#' @param A,ipm A matrix or a result from \code{make_ipm}, or \code{NULL} if \code{x}
#' is specified as the matrix or IPM object.
#' @param col A vector of colors to use for plotting
#' @param bw A logical indicating whether to use a greyscale palette for plotting
#' @param do_contour A logical indicating whether or not draw contour lines
#' on the plot
#' @param contour_cex A numeric specifying how large to make labels for the
#' contour lines.
#' @param do_legend A logical indicating whether to draw a legend for the plot
#' @param ... further arguments passed to legend
#'
#' @return \code{A} or \code{ipm} invisibly
#'
#' @details
#' If an IPM kernel is overwhelmed by information in say, a fecundity sub-kernel,
#' use the \code{exponent} argument in \code{plot.*_ipm} to make it more visually
#' appealing.
#'
#' @importFrom grDevices grey rainbow
#' @importFrom graphics abline axis contour image layout par
#' @export
plot.ipmr_matrix <- function(x = NULL, y = NULL,
A,
col = grDevices::rainbow(100, start=0.67, end=0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
contour_cex = 1,
...) {
if(missing(A)) A <- x
x <- seq_len(ncol(A))
y <- seq_len(nrow(A))
nx = length(x)
ny = length(y)
x1 = c(1.5 * x[1] - 0.5 * x[2], 1.5 * x[nx] - 0.5 * x[nx - 1])
y1 = c(1.5 * y[1] - 0.5 * y[2], 1.5 * y[ny] - 0.5 * y[ny - 1])
if(bw) col = grDevices::grey( (200:50) / 200 )
graphics::image(list(x = x,
y = y,
z = t(A)),
xlim = x1,
ylim = rev(y1),
col = col,
bty = "u",
...)
graphics::abline(v = range(x1))
graphics::abline(h = range(y1))
if(do_contour) graphics::contour(x,
y,
t(A),
nlevels = 5,
labcex = contour_cex,
add = TRUE)
if(do_legend) {
l.y = seq(min(A), max(A),length = 100)
graphics::image(list(x = 1:2,
y = l.y,
z = rbind(l.y, l.y)),
col = col,
bty = "o",
xaxt = "n",
yaxt = "n")
graphics::axis(side = 2,
cex.axis = 1.5,
at = pretty(seq(min(A), max(A), length=10)))
}
invisible(A)
}
#' @rdname plot_star
#' @param exponent The exponent to raise each kernel to. Setting this to a low
#' number can help visualize kernels that are overwhelmed by a few very large
#' numbers.
#' @param n_row,n_col If plotting multiple (sub-)kernels, how many rows and
#' columns to arrange them in.
#'
#' @export
plot.simple_di_det_ipm <- function(x = NULL, y = NULL,
ipm = NULL,
col = rainbow(100, start=0.67, end=0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
exponent = 1,
n_row = 1,
n_col = 1,
...) {
dots <- list(...)
old_par <- graphics::par(no.readonly = TRUE)
on.exit(par(old_par))
# This is used so that users can just say plot(my_model) instead of
# plot(ipm = my_model). ipmr_matrix expects x and y to both be NULL
if(!is.null(x) && is.null(ipm)){
ipm <- x
x <- NULL
}
plot_list <- ipm$sub_kernels
par(mar = c(6, 2, 3, 1),
mfrow = c(n_row, n_col))
if(do_legend) {
graphics::layout(mat = cbind(matrix(1, 5, 5), rep(2, 5)))
}
for(i in seq_along(plot_list)){
use_kern <- plot_list[[i]]
nm <- names(plot_list)[i]
plot.ipmr_matrix(x = x,
y = y,
A = use_kern ^ exponent,
col = col,
bw = bw,
do_contour = do_contour,
do_legend = do_legend,
dots,
main = nm)
}
invisible(ipm)
}
#' @rdname plot_star
#' @export
plot.simple_di_stoch_param_ipm <- function(x = NULL, y = NULL,
ipm = NULL,
col = rainbow(100, start=0.67, end=0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
exponent = 1,
n_row = 1,
n_col = 1,
...) {
dots <- list(...)
old_par <- graphics::par(no.readonly = TRUE)
on.exit(par(old_par))
# This is used so that users can just say plot(my_model) instead of
# plot(ipm = my_model). ipmr_matrix expects x and y to both be NULL
if(!is.null(x) && is.null(ipm)){
ipm <- x
x <- NULL
}
plot_list <- ipm$sub_kernels
if(all(is.na(plot_list[[1]]))) {
stop("Cannot plot sub-kernels unless 'make_ipm(return_sub_kernels = TRUE)'!")
}
plt_seq <- seq_along(plot_list)
par(mar = c(6, 2, 3, 1),
mfrow = c(n_row, n_col))
if(do_legend) {
graphics::layout(mat = cbind(matrix(1, 5, 5), rep(2, 5)))
}
for(i in seq_along(plot_list)){
use_kern <- plot_list[[i]]
nm <- names(plot_list)[i]
plot.ipmr_matrix(x = x,
y = y,
A = use_kern ^ exponent,
col = col,
bw = bw,
do_contour = do_contour,
do_legend = do_legend,
dots,
main = nm)
}
invisible(ipm)
}
#' @rdname plot_star
#' @export
plot.simple_di_stoch_kern_ipm <- function(x = NULL, y = NULL,
ipm = NULL,
col = rainbow(100, start=0.67, end=0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
exponent = 1,
n_row = 1,
n_col = 1,
...) {
dots <- list(...)
old_par <- graphics::par(no.readonly = TRUE)
on.exit(par(old_par))
# This is used so that users can just say plot(my_model) instead of
# plot(ipm = my_model). ipmr_matrix expects x and y to both be NULL
if(!is.null(x) && is.null(ipm)){
ipm <- x
x <- NULL
}
plot_list <- ipm$sub_kernels
par(mar = c(6, 2, 3, 1),
mfrow = c(n_row, n_col))
if(do_legend) {
graphics::layout(mat = cbind(matrix(1, 5, 5), rep(2, 5)))
}
for(i in seq_along(plot_list)){
use_kern <- plot_list[[i]]
nm <- names(plot_list)[i]
plot.ipmr_matrix(x = x,
y = y,
A = use_kern ^ exponent,
col = col,
bw = bw,
do_contour = do_contour,
do_legend = do_legend,
dots,
main = nm)
}
invisible(ipm)
}
#' @rdname plot_star
#' @inheritParams format_mega_kernel
#'
#' @export
plot.general_di_det_ipm <- function(x = NULL, y = NULL,
ipm = NULL,
mega_mat = NA_character_,
col = rainbow(100,
start = 0.67,
end = 0),
bw = FALSE,
do_contour = FALSE,
do_legend = FALSE,
exponent = 1,
n_row = 1,
n_col = 1,
...) {
dots <- list(...)
old_par <- graphics::par(no.readonly = TRUE)
on.exit(par(old_par))
# This is used so that users can just say plot(my_model) instead of
# plot(ipm = my_model). ipmr_matrix expects x and y to both be NULL
if(!is.null(x) && is.null(ipm)){
ipm <- x
x <- NULL
}
mega_mat <- rlang::enquo(mega_mat)
if(rlang::quo_text(mega_mat) == "NA") {
stop("Plotting general IPMs requires building a 'mega_mat'.\n",
"Please specify an expression for the 'mega_mat' argument.")
}
graphics::par(mar = c(6, 5, 3, 2),
mfrow = c(n_row, n_col))
plot_list <- format_mega_kernel(ipm, mega_mat = !! mega_mat)
if(do_legend) {
graphics::layout(mat = cbind(matrix(1, 5, 5), rep(2, 5)))
}
lapply(plot_list, function(ipm) plot.ipmr_matrix(x = x,
y = y,
A = ipm ^ exponent,
col = col,
bw = bw,
do_contour = do_contour,
do_legend = do_legend,
dots))
invisible(ipm)
}
.ncol_nrow <- function(plt_seq) {
out <- list(ncol = NA,
nrow = NA)
if(length(plt_seq) > 25) {
out$ncol <- 5
out$nrow <- 5
} else {
out$ncol <- out$nrow <- ceiling(sqrt(length(plt_seq)))
}
return(out)
}
# right_ev ----------------
#' @rdname eigenvectors
#'
#' @title Compute the standardized left and right eigenvectors via iteration
#'
#' @param ipm Output from \code{make_ipm()}.
#' @param ... Other arguments passed to methods
#'
#' @return A list of named numeric vector(s) corresponding to the stable trait distribution
#' function (\code{right_ev}) or the reproductive values for each trait (\code{left_ev}).
#'
#' @export
right_ev <- function(ipm, ... ) {
UseMethod('right_ev')
}
#' @rdname eigenvectors
#' @param iterations The number of times to iterate the model to reach
#' convergence. Default is 100.
#' @param tolerance Tolerance to evaluate convergence to asymptotic dynamics.
#'
#' @export
right_ev.simple_di_det_ipm <- function(ipm,
iterations = 100,
tolerance = 1e-10,
...) {
mod_nm <- deparse(substitute(ipm))
# if it's already been iterated to convergence, we don't have much work to do.
if(.already_iterated(ipm)) {
# get index for population vector of final iteration
final_it <- dim(ipm$pop_state[[1]])[2]
if(all(is_conv_to_asymptotic(ipm,
tolerance = tolerance))) {
out <- .extract_conv_ev_simple(ipm$pop_state)
} else {
# If it's iterated, but not converged, we can at least use the final
# iteration which is hopefully closer to convergence than the user-defined
# initial population vector
message(
paste(
"'",
mod_nm,
"'",
' did not converge to asymptotic dynamics after ',
final_it,
' iterations.\n',
'Will re-iterate the model ',
iterations,
' times and check for convergence.',
sep = ""
)
)
init_pop_vec <- lapply(ipm$pop_state[!grepl("lambda", names(ipm$pop_state))],
function(x, final_it) {
x[ , final_it]
},
final_it = ncol(ipm$pop_state[[1]]))
#
# pop_nm <- ipm$proto_ipm$state_var %>%
# unlist() %>%
# unique() %>%
# .[1] %>%
# paste('n_', ., sep = "")
#
# pop_nm <- rlang::ensym(pop_nm)
test_conv <- ipm$proto_ipm %>%
define_pop_state(!!! init_pop_vec) %>%
make_ipm(iterate = TRUE,
iterations = iterations)
if(all(is_conv_to_asymptotic(test_conv,
tolerance = tolerance))) {
out <- .extract_conv_ev_simple(test_conv$pop_state)
message('model is now converged :)')
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
final_it + iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
}
} else {
message(
paste(
"'",
mod_nm,
"'",
' has not been iterated yet. ',
'Generating a population vector using runif() and\niterating the model ',
iterations,
' times to check for convergence to asymptotic dynamics',
sep = ""
)
)
# A model that hasn't been iterated yet
# Create variables for internal usage
pop_nm <- .get_pop_nm_simple(ipm)
pop_states <- paste('n', pop_nm, sep = '_')
# Final step is to generate an initial pop_state. This is always just
# vector drawn from a random uniform distribution
len_pop_state <- dim(ipm$sub_kernels[[1]])[1]
init_pop <- stats::runif(len_pop_state)
# Drop _t for defining the initial population vector so define_pop_state
# doesn't complain
test_conv <- ipm$proto_ipm %>%
define_pop_state(!! pop_states[1] := init_pop) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
normalize_pop_size = TRUE)
if(all(is_conv_to_asymptotic(test_conv,
tolerance = tolerance))) {
out <- .extract_conv_ev_simple(test_conv$pop_state)
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
}
# Stuff into a list and standardize
# out <- rlang::list2(!! out_nm := (out / sum(out)))
out <- Filter(function(x) !any(is.na(x)), out)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'w', sep = '_')
class(out) <- "ipmr_w"
return(out)
}
#' @rdname eigenvectors
#' @param burn_in The proportion of early iterations to discard from the
#' stochastic simulation
#' @export
right_ev.simple_di_stoch_kern_ipm <- function(ipm,
burn_in = 0.25,
...) {
mod_nm <- deparse(substitute(ipm))
# Identify state variable name
pop_nm <- .get_pop_nm_simple(ipm)
out <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq],
burn_seq = burn_in_seq)
# Normalize just in case. easy for simple IPMs where we know there's only
# one trait contributing to total population size. Tricker in general IPMs.
out[[1]] <- apply(out[[1]], 2, function(x) x / sum(x))
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'w', sep = '_')
class(out) <- "ipmr_w"
return(out)
}
#' @rdname eigenvectors
#' @export
right_ev.simple_di_stoch_param_ipm <- right_ev.simple_di_stoch_kern_ipm
#' @rdname eigenvectors
#' @export
right_ev.general_di_det_ipm <- function(ipm,
iterations = 100,
tolerance = 1e-10,
...) {
mod_nm <- deparse(substitute(ipm))
final_it <- dim(ipm$pop_state[[1]])[2]
if(all(is_conv_to_asymptotic(ipm,
tolerance = tolerance))) {
out <- .extract_conv_ev_general(ipm$pop_state, ipm$proto_ipm)
} else {
# Not yet converged, so re-iterate the model and see what we get
message(
paste(
"'",
mod_nm,
"'",
' did not converge to asymptotic dynamics after ',
final_it,
' iterations.\n',
'Will re-iterate the model ',
iterations,
' times and check for convergence.',
sep = ""
)
)
use_pop_state <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
init_pop_vec <- lapply(use_pop_state,
function(x, final_it) x[ , final_it],
final_it = final_it)
# Remove NAs generated by the unsubstituted pop_state place holders
init_pop_vec <- Filter(f = function(y) !any(is.na(y)), x = init_pop_vec)
names(init_pop_vec) <- gsub('pop_state', 'n', names(init_pop_vec))
test_conv <- ipm$proto_ipm %>%
define_pop_state(
pop_vectors = init_pop_vec
) %>%
make_ipm(iterate = TRUE,
iterations = iterations)
if(all(is_conv_to_asymptotic(test_conv,
tolerance = tolerance))) {
out <- .extract_conv_ev_general(test_conv$pop_state, test_conv$proto_ipm)
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
}
out <- Filter(function(x) !any(is.na(x)), out)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'w', sep = '_')
class(out) <- "ipmr_w"
return(out)
}
#' @rdname eigenvectors
#' @export
right_ev.general_di_stoch_kern_ipm <- function(ipm,
burn_in = 0.25,
...) {
mod_nm <- deparse(substitute(ipm))
out <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq, drop = FALSE],
burn_seq = burn_in_seq)
# Normalize everything just in case
pop_sizes <- lapply(out, colSums) %>%
do.call(what = ".add", args = .)
for(i in seq_along(out)) {
for(j in seq_len(ncol(out[[1]]))) {
out[[i]][ , j] <- out[[i]][ ,j] / pop_sizes[j]
}
}
out <- Filter(function(x) !any(is.na(x)), out)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'w', sep = '_')
class(out) <- "ipmr_w"
return(out)
}
#' @rdname eigenvectors
#' @export
right_ev.general_di_stoch_param_ipm <- right_ev.general_di_stoch_kern_ipm
# left_ev -----------------
#' @rdname eigenvectors
#' @export
left_ev <- function(ipm, ...) {
UseMethod('left_ev')
}
#' @export
#' @rdname eigenvectors
#' @importFrom stats runif
left_ev.simple_di_det_ipm <- function(ipm,
iterations = 100,
tolerance = 1e-10,
...) {
mod_nm <- deparse(substitute(ipm))
# Identify state variable name
pop_nm <- .get_pop_nm_simple(ipm)
if(any(ipm$proto_ipm$uses_par_sets)) {
ps_nms <- .flatten_to_depth(ipm$proto_ipm$par_set_indices, 1L) %>%
names() %>%
unique()
ps_nms <- ps_nms[ps_nms != "levels"]
ps_suff <- paste(ps_nms, collapse = "_")
pop_nm <- paste(pop_nm, ps_suff, sep = "_")
}
# Create variables for internal usage
pop_states <- paste('n', pop_nm, sep = '_')
# Final step is to generate an initial pop_state. This is always just
# vector drawn from a random uniform distribution
len_pop_state <- dim(ipm$sub_kernels[[1]])[1]
init_pop <- stats::runif(len_pop_state)
test_conv <- ipm$proto_ipm %>%
define_pop_state(!! pop_states[1] := init_pop) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
iteration_direction = "left",
normalize_pop_size = TRUE)
if(all(is_conv_to_asymptotic(test_conv,
tolerance = tolerance))) {
out <- .extract_conv_ev_simple(test_conv$pop_state)
# out_nm <- paste(pop_nm, 'v', sep = "_")
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
# Stuff into a list and standardize
out <- Filter(function(x) !any(is.na(x)), out)
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'v', sep = '_')
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#' @param kernel_seq The sequece of parameter set indices used to select kernels
#' during the iteration procedure. If \code{NULL}, will use the sequence stored
#' in the \code{ipm} object. Should usually be left as \code{NULL}.
#' @export
left_ev.simple_di_stoch_kern_ipm <- function(ipm,
iterations = 10000,
burn_in = 0.25,
kernel_seq = NULL,
...) {
mod_nm <- deparse(substitute(ipm))
# Identify state variable name
pop_nm <- .get_pop_nm_simple(ipm)
# Create variables for internal usage
pop_states <- paste('n', pop_nm, c('t', 't_1'), sep = '_')
# Final step is to generate an initial pop_state. This is always just
# vector drawn from a random uniform distribution
len_pop_state <- nrow(ipm$sub_kernels[[1]])
init_pop <- stats::runif(len_pop_state)
if(is.null(kernel_seq)) kernel_seq <- ipm$env_seq
# Drop _t for defining the initial population vector so define_pop_state
# doesn't complain
pop_states[1] <- gsub("_t$", "", pop_states[1])
test_conv <- ipm$proto_ipm %>%
define_pop_state(!! pop_states[1] := init_pop) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
iteration_direction = "left",
normalize_pop_size = TRUE,
kernel_seq = kernel_seq)
out <- test_conv$pop_state
out$lambda <- NULL
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq],
burn_seq = burn_in_seq)
out_nm <- paste(pop_nm, 'v', sep = "_")
names(out) <- out_nm
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#'
#' @section \strong{Deterministic eigenvectors}:
#' For \code{right_ev}, if the model has already been iterated and has
#' converged to asymptotic dynamics, then it will just extract the final
#' population state and return that in a named list. Each element of the list
#' is a vector with length \code{>= 1} and corresponds each state variable's
#' portion of the eigenvector.
#' If the model has been iterated, but has not yet converged to asymptotic dynamics,
#' \code{right_ev} will try to iterate it further using the final population state
#' as the starting point. The default number of iterations is 100, and can be
#' adjusted using the \code{iterations} argument.
#' If the model hasn't been iterated, then \code{right_ev} will try iterating it
#' for \code{iterations} number of time steps and check for convergence. In the
#' latter two cases, if the model still has not converged to asymptotic dynamics,
#' it will return \code{NA} with a warning.
#'
#' For \code{left_ev}, the transpose iteration (\emph{sensu} Ellner & Rees 2006,
#' Appendix A) is worked out based on the \code{state_start} and \code{state_end}
#' in the model's \code{proto_ipm} object. The model is then iterated for
#' \code{iterations} times to produce a standardized left eigenvector.
#'
#' @section \strong{Stochastic eigenvectors}:
#' \code{left_ev} and \code{right_ev} return different things for stochastic models.
#' \code{right_ev} returns the trait distribution through time from the stochastic
#' simulation (i.e. \code{ipm$pop_state}), and normalizes it such that the
#' distribution at each time step integrates to 1 (if it is not already).
#' It then discards the first \code{burn_in * iterations} time steps of the
#' simulation to eliminate transient dynamics. See Ellner, Childs, & Rees 2016,
#' Chapter 7.5 for more details.
#'
#' \code{left_ev} returns a similar result as \code{right_ev}, except the trait
#' distributions are the result of left multiplying the kernel and trait
#' distribution. See Ellner, Childs, & Rees 2016, Chapter 7.5 for more
#' details.
#'
#' @export
left_ev.general_di_det_ipm <- function(ipm,
iterations = 100,
tolerance = 1e-10,
...) {
mod_nm <- deparse(substitute(ipm))
use_pop_state <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
init_pop_vec <- lapply(use_pop_state, function(x) x[ , 1])
# Remove NAs generated by the unsubstituted pop_state place holders
init_pop_vec <- Filter(f = function(y) !any(is.na(y)), x = init_pop_vec)
names(init_pop_vec) <- gsub('pop_state', 'n', names(init_pop_vec))
test_conv <- ipm$proto_ipm %>%
define_pop_state(
pop_vectors = init_pop_vec
) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
iteration_direction = "left",
normalize_pop_size = TRUE)
if(all(is_conv_to_asymptotic(test_conv, tolerance = tolerance))) {
out <- .extract_conv_ev_general(test_conv$pop_state, test_conv$proto_ipm)
} else {
warning(
paste(
"'",
mod_nm,
"'",
' did not converge after ',
iterations,
' iterations. Returning NA, please try again with more iterations.',
sep = ""
),
call. = FALSE
)
return(NA_real_)
}
out <- Filter(function(x) !any(is.na(x)), out)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'v', sep = '_')
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#' @export
left_ev.general_di_stoch_kern_ipm <- function(ipm,
iterations = 10000,
burn_in = 0.25,
kernel_seq = NULL,
...) {
mod_nm <- deparse(substitute(ipm))
use_pop_state <- ipm$pop_state[!grepl("lambda", names(ipm$pop_state))]
init_pop_vec <- lapply(use_pop_state, function(x) x[ , 1])
names(init_pop_vec) <- gsub('pop_state', 'n', names(init_pop_vec))
if(is.null(kernel_seq)) kernel_seq <- ipm$env_seq
test_conv <- ipm$proto_ipm %>%
define_pop_state(
pop_vectors = init_pop_vec
) %>%
make_ipm(iterate = TRUE,
iterations = iterations,
iteration_direction = "left",
normalize_pop_size = TRUE,
kernel_seq = kernel_seq)
out <- test_conv$pop_state[!grepl('lambda', names(test_conv$pop_state))]
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq, drop = FALSE],
burn_seq = burn_in_seq)
# Replaces the leading n_ with a trailing _v
names(out) <- substr(names(out), 3, nchar(names(out)))
names(out) <- paste(names(out), 'v', sep = '_')
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#' @export
# Needs to re-iterate the model, but use the same kernel sequence that is used
# in right_ev/make_ipm. In effect, we want a deterministic model using the
# sub-kernels slot in lieu of parameters.
left_ev.general_di_stoch_param_ipm <- function(ipm,
iterations = 10000,
burn_in = 0.25,
kernel_seq = NULL,
...) {
mod_nm <- deparse(substitute(ipm))
sub_kernels <- ipm$sub_kernels
proto_ipm <- ipm$proto_ipm
usr_its <- ncol(ipm$pop_state[[1]]) - 1
if(usr_its < iterations) {
stop("'iterations' should be less than or equal to the number of iterations in 'make_ipm()'!")
}
# Set up the model iteration expressions and initial population state.
# This is a leaner version of make_ipm.general_di_stoch_param, because
# we don't need all the user functions or other stuff - just the basic
# infrastructure to house model iteration evaluation.
proto_list <- .initialize_kernels(proto_ipm,
iterate = TRUE,
iter_dir = "left")
others <- proto_list$others
k_row <- proto_list$k_row
main_env <- .make_main_env(others$domain,
proto_ipm$usr_funs[[1]],
age_size = .uses_age(others))
temp <- .prep_di_output(others, k_row, proto_ipm,
iterations, normal= TRUE)
main_env <- .add_pop_state_to_main_env(temp$pop_state, main_env)
# Throw this in just in case we have parameter set indices in addition
# to continuous environmental variation.
kern_seq <- .make_kern_seq(others, iterations, kernel_seq)
# Create an index to select sub-kernels based on which iteration they were
# generated in.
kern_it_ind <- vapply(names(sub_kernels), function(x) {
as.integer(strsplit(x, "(_it_)")[[1]][2])
}, integer(1L))
for(i in seq_len(iterations)) {
# This selects sub-kernels by matching the _it_XXX suffix to the current
# iteration.
use_kerns <- sub_kernels[kern_it_ind == i]
# Now, drop the _it_XXX suffix and pretend they're just regular kernels
# generated by .make_sub_kernel_general_lazy.
rm_regex <- paste("_it_", i, '$', sep = "")
names(use_kerns) <- gsub(rm_regex, "", names(use_kerns))
# We don't care about any of the other output from model iteration this,
# so skip all the naming/stashing other outputs.
pop_state <- .iterate_model(proto_ipm = proto_ipm,
k_row = k_row,
sub_kern_list = use_kerns,
current_iteration = i,
kern_seq = kern_seq,
pop_state = temp$pop_state,
main_env = main_env,
normal = TRUE)
temp$pop_state <- pop_state
}
out <- temp$pop_state
out[grepl("lambda", names(out))] <- NULL
burn_in_seq <- seq_len(ncol(out[[1]])) %>%
.[seq.int(from = ceiling(ncol(out[[1]]) * burn_in),
to = ncol(out[[1]]),
by = 1)]
out <- lapply(out,
function(x, burn_seq) x[ , burn_seq, drop = FALSE],
burn_seq = burn_in_seq)
names(out) <- gsub("pop_state_", "", names(out))
names(out) <- paste(names(out), 'v', sep = '_')
class(out) <- "ipmr_v"
return(out)
}
#' @rdname eigenvectors
#' @export
left_ev.simple_di_stoch_param_ipm <- left_ev.general_di_stoch_param_ipm
#' @noRd
# Helper to be used in do.call(".add", X) where X is a list with potentially more
# than two elements. "+" won't work in this case, and sum() will return a single
# number.
.add <- function(...) {
dots <- list(...)
Reduce("+", x = dots, init = 0)
}
#' @rdname as_matrix
#' @title Convert to bare matrices
#' @description Converts objects to \code{c("matrix", "array")}.
#' @param x An object of class \code{ipmr_matrix}, or the output from
#' \code{make_ipm}.
#' @param ... ignored.
#'
#' @return A matrix.
#' @export
as.matrix.ipmr_matrix <- function(x, ...) {
class(x) <- c("matrix", "array")
return(x)
}
#' @rdname as_matrix
#' @export
as.matrix.ipmr_ipm <- function(x, ...) {
lapply(x$sub_kernels, as.matrix)
}
Try the ipmr 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.