Skip to content

Commit

Permalink
Finished updating visualisation functions and reorganised code...
Browse files Browse the repository at this point in the history
  • Loading branch information
@cvisintin
cvisintin committed Apr 26, 2021
1 parent 9698462 commit 7006cbe
Show file tree
Hide file tree
Showing 13 changed files with 881 additions and 144 deletions.
5 changes: 5 additions & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,11 @@ export(kernel_dispersal)
export(landscape)
export(modified_transition)
export(mortality)
export(plot_hab_spatial)
export(plot_k_spatial)
export(plot_k_trend)
export(plot_pop_spatial)
export(plot_pop_trend)
export(population_dynamics)
export(set_proportion_dispersing)
export(simulation)
Expand Down
128 changes: 6 additions & 122 deletions R/simulation_results-class.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,7 +248,7 @@ plot.simulation_results <- function (x,
}


#' Extract objects from a 'simulation_results' object
#' Extract spatial object from a 'simulation_results' object
#'
#' The simulation results object is a list of lists containing spatial (and other) objects and
#' is organised by the following tree diagram:
Expand Down Expand Up @@ -278,7 +278,7 @@ plot.simulation_results <- function (x,
#' @param x a simulation_results object
#' @param replicate which replicate to extract from a \code{simulation_results}
#' object
#' @param timestep timestep(s) to extract from a \code{simulation_results}
#' @param timestep which timestep to extract from a \code{simulation_results}
#' @param landscape_object which landscape object to extract from a
#' \code{simulation_results} object - can be specified by name
#' (e.g. "suitability") or index number
Expand Down Expand Up @@ -308,7 +308,7 @@ plot.simulation_results <- function (x,
#'
#' # Extract the population raster for the second life-stage in the first
#' # replicate and ninth timestep
#' extract_spatial(sim, timestep = 9, stage = 2)
#' extract_spatial(sim, replicate = 1, timestep = 9, stage = 2)
#' }

extract_spatial <- function (x,
Expand All @@ -318,15 +318,15 @@ extract_spatial <- function (x,
stage = 1,
misc = 1) {
if (landscape_object == 1 | landscape_object == "population") {
x[[replicate]][[timestep]][[landscape_object]][[stage]]
return(x[[replicate]][[timestep]][[landscape_object]][[stage]])
}
if (landscape_object > 3 |
!landscape_object == "population" |
!landscape_object == "suitability" |
!landscape_object == "carrying_capacity") {
x[[replicate]][[timestep]][[landscape_object]][[misc]]
return(x[[replicate]][[timestep]][[landscape_object]][[misc]])
}

return(x[[replicate]][[timestep]][[landscape_object]])
}

# #' @rdname simulation
Expand All @@ -350,122 +350,6 @@ extract_spatial <- function (x,
#
# }

#' Compare minimum expected populations
#'
#' Compare minimum expected populations from two or more 'simulation_results' objects.
#'
#' @param x a simulation_results object
#' @param ... additional simulation results objects
#' @param interval the desired confidence interval representing the uncertainty around
#' the expected minimum population estimates from simulation comparisons; expressed as
#' a whole integer between 0 and 100 (default value is 95).
#' @param all_points should the expected minimum populations from all simulation
#' replicates be shown on the plot?
#' @param simulation_names an optional character vector of simulation names to override
#' the defaults
#'
#' @export
#'
#' @examples
#'
#' \dontrun{
#' ls <- landscape(population = egk_pop, suitability = egk_hab, carrying_capacity = egk_k)
#'
#' # Create populations dynamics with and without ceiling density dependence
#' pd1 <- population_dynamics(change = growth(egk_mat),
#' dispersal = kernel_dispersal(max_distance = 1000,
#' dispersal_kernel = exponential_dispersal_kernel(distance_decay = 500)),
#' density_dependence = ceiling_density())
#' pd2 <- population_dynamics(change = growth(egk_mat),
#' dispersal = kernel_dispersal(max_distance = 3000,
#' dispersal_kernel = exponential_dispersal_kernel(distance_decay = 1500)))
#'
#' # Run first simulation with ceiling density dependence and three replicates
#' sim1 <- simulation(landscape = ls,
#' population_dynamics = pd1,
#' habitat_dynamics = NULL,
#' timesteps = 20,
#' replicates = 3)
#'
#' # Run second simulation without ceiling density dependence and three replicates
#' sim2 <- simulation(landscape = ls,
#' population_dynamics = pd2,
#' habitat_dynamics = NULL,
#' timesteps = 20,
#' replicates = 3)
#'
#' compare_emp(sim1, sim2)
#' }

compare_emp <- function (x, ..., interval = 95, all_points = FALSE, simulation_names = NULL) {

# read in simulation objects to compare
sim_objects <- list(x, ...)
n_objects <- length(sim_objects)

# get names of simulations
sim_names <- as.character(substitute(list(x, ...)))[-1L]

interval_range <- c((100 - interval) / 2, 100 - (100 - interval) / 2) / 100

# initiate table of values
df <- data.frame("name" = sim_names,
"emp_mean" = NA,
"emp_lower" = NA,
"emp_upper" = NA)

if(is.null(simulation_names)) simulation_names <- sim_names

# populate table with emp mean and error values
for (i in seq_len(n_objects)){
pops <- get_pop_simulation(sim_objects[[i]])
min_total_pops <- apply(pops, 3, function(x) min(rowSums(x)))
emp_mean <- mean(min_total_pops)
emp_lower <- stats::quantile(min_total_pops, interval_range)[1]
emp_upper <- stats::quantile(min_total_pops, interval_range)[2]
df[i, -1] <- c(emp_mean, emp_lower, emp_upper)
}

graphics::par(mar=c(4, 4.5, 1.5, 1.5) + 0.1)

graphics::plot(NULL,
xlim = c(0.5, n_objects + 0.5),
xaxt = "n",
xlab = "Simulation Name",
ylim = range(c(df$emp_lower, df$emp_upper)),
yaxt = "n",
ylab = "",
main = "",
lwd = 0.5)
if (all_points == TRUE) {
for (i in seq_len(n_objects)){
pops <- get_pop_simulation(sim_objects[[i]])
min_total_pops <- apply(pops, 3, function(x) min(rowSums(x)))
graphics::points(jitter(rep(i, length(min_total_pops))),
min_total_pops,
col = "lightgrey",
pch = 19,
cex = 0.8)
}
}
graphics::points(seq_len(n_objects),
df$emp_mean,
pch = 19)
graphics::arrows(seq_len(n_objects),
df$emp_lower,
seq_len(n_objects),
df$emp_upper,
length=0.05,
angle=90,
code=3)
graphics::axis(1, at = seq_len(n_objects), labels = simulation_names)
graphics::axis(2, at = pretty(range(c(df$emp_lower, df$emp_upper)), 5))
graphics::mtext(paste0("Minimum Population (", interval, "% Interval)"),
side = 2,
line = 2.5)
}


##########################
### internal functions ###
##########################
Expand Down
38 changes: 20 additions & 18 deletions R/utils.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,24 +18,32 @@ round_pop <- function (population) {

}

get_pop_replicate <- function(x, ...) {
total_stages <- raster::nlayers(x[[1]]$population)
idx <- which(!is.na(raster::getValues(x[[1]]$population[[1]])))
pops <- lapply(x, function(x) raster::extract(x$population, idx))
get_pop_replicate <- function(replicate_result, init_pops, ...) {
total_stages <- raster::nlayers(replicate_result[[1]]$population)
idx <- which(!is.na(raster::getValues(replicate_result[[1]]$population[[1]])))

init_pops <- raster::extract(init_pops, idx)
init_pop_sums <- colSums(init_pops)

pops <- lapply(replicate_result, function(x) raster::extract(x$population, idx))
pop_sums <- lapply(pops, function(x) colSums(x))

pop_matrix <- matrix(unlist(pop_sums), ncol = total_stages, byrow = TRUE)
pop_matrix <- unname(rbind(init_pop_sums, pop_matrix))

return(pop_matrix)
}

get_pop_simulation <- function(x, ...) {
total_stages <- raster::nlayers(x[[1]][[1]]$population)
timesteps <- length(x[[1]])
sims <- length(x)

pop_array <- array(dim=c(timesteps, total_stages, sims))
get_pop_simulation <- function(sim_result, ...) {
total_stages <- raster::nlayers(sim_result[[1]][[1]]$population)
timesteps <- length(sim_result[[1]])
sims <- length(sim_result)

pop_array <- array(dim = c(timesteps + 1, total_stages, sims))

for(i in seq_len(sims)) {
pop_array[, , i] <- get_pop_replicate(x[[i]])
pop_array[, , i] <- get_pop_replicate(replicate_result = sim_result[[i]],
init_pops = attr(sim_result, "initial_population"))
}
return(pop_array)
}
Expand Down Expand Up @@ -116,12 +124,6 @@ rmultinom_large_int <- function (population) {

}

pretty_int <- function (...) {
at <- pretty(...)
at <- at[at %% 1 == 0]
at[at != 0]
}

int_or_proper_length_vector <- function (input, n_stages, parameter) {
if (length(input) != 1 & length(input) != n_stages) {
stop(paste0("Please provide either a single number or vector of ",
Expand Down Expand Up @@ -150,4 +152,4 @@ global_object_error <- function(error) {
}

stop (message, call. = FALSE)
}
}
Loading

0 comments on commit 7006cbe

Please sign in to comment.