## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7.5,
  fig.height = 4.8,
  warning = FALSE,
  message = FALSE
)

## ----install, eval=FALSE------------------------------------------------------
#  install.packages("remotes")
#  remotes::install_github("agimenezromero/demovuln-r", build_vignettes = TRUE)

## ----packages-----------------------------------------------------------------
library(demovuln)

if (!requireNamespace("ggplot2", quietly = TRUE)) {
  stop("This tutorial uses ggplot2 for plotting. Install it with install.packages('ggplot2').")
}
library(ggplot2)

## ----define-matrix------------------------------------------------------------
stage_names <- c("Juvenile", "Subadult", "Adult")

A <- matrix(
  c(
    0.00, 0.00, 1.60,
    0.35, 0.55, 0.05,
    0.00, 0.25, 0.86
  ),
  nrow = 3,
  byrow = TRUE,
  dimnames = list(
    destination = stage_names,
    source = stage_names
  )
)

A

## ----create-model-------------------------------------------------------------
model <- matrix_population_model(
  A,
  fecundity_rows = 1,
  adult_stages = 3,
  juvenile_stages = c(1, 2),
  name = "Example three-stage model"
)

model$lambda_
stable_stage_distribution(model)

## ----target-masks-------------------------------------------------------------
targets <- c("adult_survival", "juvenile_survival", "fecundity", "all")
target_labels <- c(
  adult_survival = "Adult survival",
  juvenile_survival = "Juvenile survival",
  fecundity = "Fecundity",
  all = "All parameters"
)

mask_to_data <- function(mask, target_label) {
  idx <- expand.grid(
    destination = seq_len(nrow(mask)),
    source = seq_len(ncol(mask))
  )
  idx$selected <- as.vector(mask)
  idx$target <- target_label
  idx
}

mask_table <- do.call(
  rbind,
  lapply(targets, function(target) {
    mask <- build_target_mask(model, target = target)
    mask_to_data(mask, target_labels[[target]])
  })
)

mask_table$source_stage <- factor(mask_table$source, labels = stage_names)
mask_table$destination_stage <- factor(mask_table$destination, labels = stage_names)

p_masks <- ggplot(mask_table, aes(x = source_stage, y = destination_stage, fill = selected)) +
  geom_tile(color = "grey70") +
  facet_wrap(~ target, nrow = 1) +
  scale_y_discrete(limits = rev(stage_names)) +
  scale_fill_manual(values = c(`FALSE` = "white", `TRUE` = "steelblue")) +
  labs(
    x = "Source stage at time t",
    y = "Destination stage at time t + 1",
    fill = "Perturbed",
    title = "Matrix entries selected by each demographic target"
  ) +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p_masks

## ----survival-affects-fecundity-----------------------------------------------
B_with_fecundity <- apply_perturbation(
  model,
  target = "adult_survival",
  magnitude = 0.5,
  survival_affects_fecundity = TRUE
)

B_without_fecundity <- apply_perturbation(
  model,
  target = "adult_survival",
  magnitude = 0.5,
  survival_affects_fecundity = FALSE
)

B_with_fecundity
B_without_fecundity

## ----single-simulation--------------------------------------------------------
generation_time <- 20

sim_adult <- simulate_dynamics(
  model,
  target = "adult_survival",
  magnitude = 0.30,
  duration = 4,
  period = 10,
  t_max = 3 * generation_time,
  recovery_steps = 4 * generation_time,
  normalize_by_lambda = TRUE
)

sim_adult$reduction

## ----single-trajectory-plot---------------------------------------------------
trajectory_df <- data.frame(
  time = seq_along(sim_adult$abundance) - 1,
  time_generations = (seq_along(sim_adult$abundance) - 1) / generation_time,
  baseline = sim_adult$baseline_abundance,
  perturbed = sim_adult$abundance
)

trajectory_plot_df <- rbind(
  data.frame(
    time_generations = trajectory_df$time_generations,
    abundance = trajectory_df$baseline,
    scenario = "Unperturbed baseline"
  ),
  data.frame(
    time_generations = trajectory_df$time_generations,
    abundance = trajectory_df$perturbed,
    scenario = "Perturbed"
  )
)

p_single <- ggplot(trajectory_plot_df, aes(x = time_generations, y = abundance, linetype = scenario)) +
  geom_line(linewidth = 0.9) +
  labs(
    x = "Time (reference generations)",
    y = "Relative population size",
    linetype = NULL,
    title = "One temporally structured perturbation regime"
  ) +
  theme_minimal(base_size = 12)

p_single

## ----compare-target-trajectories----------------------------------------------
simulations <- lapply(
  c("adult_survival", "juvenile_survival", "fecundity"),
  function(target) {
    simulate_dynamics(
      model,
      target = target,
      magnitude = 0.30,
      duration = 4,
      period = 10,
      t_max = 3 * generation_time,
      recovery_steps = 4 * generation_time,
      normalize_by_lambda = TRUE
    )
  }
)
names(simulations) <- c("adult_survival", "juvenile_survival", "fecundity")

trajectory_panel_data <- do.call(
  rbind,
  lapply(names(simulations), function(target) {
    sim <- simulations[[target]]
    time <- seq_along(sim$abundance) - 1
    rbind(
      data.frame(
        time_generations = time / generation_time,
        abundance = sim$baseline_abundance,
        scenario = "Unperturbed baseline",
        target = target_labels[[target]]
      ),
      data.frame(
        time_generations = time / generation_time,
        abundance = sim$abundance,
        scenario = "Perturbed",
        target = target_labels[[target]]
      )
    )
  })
)

p_trajectories <- ggplot(
  trajectory_panel_data,
  aes(x = time_generations, y = abundance, linetype = scenario)
) +
  geom_line(linewidth = 0.85) +
  facet_wrap(~ target, nrow = 1) +
  labs(
    x = "Time (reference generations)",
    y = "Relative population size",
    linetype = NULL,
    title = "The same perturbation regime applied to different demographic targets"
  ) +
  theme_minimal(base_size = 11)

p_trajectories

## ----compare-target-reductions------------------------------------------------
regime_reductions <- data.frame(
  target = target_labels[names(simulations)],
  population_reduction = vapply(simulations, function(x) x$reduction, numeric(1))
)

regime_reductions

## ----direct-grid--------------------------------------------------------------
grid_direct <- perturbation_grid(
  magnitudes = seq(0, 0.8, by = 0.1),
  durations = seq(0, generation_time, by = 2),
  periods = c(5, 10, 20, 40)
)

grid_scenarios(grid_direct)[1:10, ]

## ----frequency-grid-----------------------------------------------------------
frequencies <- c(0.25, 0.5, 1, 2, 4)

grid <- perturbation_grid_from_frequencies(
  magnitudes = seq(0, 0.8, by = 0.1),
  durations = seq(0, generation_time, by = 2),
  frequencies = frequencies,
  generation_time = generation_time,
  rounding = "nearest"
)

frequency_period_table <- data.frame(
  frequency_per_generation = frequencies,
  period = pmax(1L, as.integer(round(generation_time / frequencies)))
)

frequency_period_table

## ----run-grid-one-target------------------------------------------------------
out_adult <- run_grid(
  model,
  target = "adult_survival",
  grid = grid,
  t_max = 3 * generation_time,
  recovery_steps = 4 * generation_time,
  normalize_by_lambda = TRUE,
  skip_infeasible = TRUE
)

out_adult$vulnerability
head(out_adult$table)

## ----run-grid-all-targets-----------------------------------------------------
grid_results <- lapply(targets, function(target) {
  run_grid(
    model,
    target = target,
    grid = grid,
    t_max = 3 * generation_time,
    recovery_steps = 4 * generation_time,
    normalize_by_lambda = TRUE,
    skip_infeasible = TRUE
  )
})
names(grid_results) <- targets

vulnerability_table <- data.frame(
  target = unname(target_labels[targets]),
  Phi = vapply(grid_results, function(x) x$vulnerability, numeric(1))
)

vulnerability_table

## ----vulnerability-barplot----------------------------------------------------
p_phi <- ggplot(vulnerability_table, aes(x = target, y = Phi)) +
  geom_col(width = 0.7, fill = "grey40") +
  labs(
    x = "Perturbed demographic target",
    y = expression(Phi),
    title = "Integrated vulnerability across demographic targets"
  ) +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 25, hjust = 1))

p_phi

## ----grid-table-for-heatmaps--------------------------------------------------
grid_table <- do.call(
  rbind,
  lapply(targets, function(target) {
    tab <- grid_results[[target]]$table
    tab$target <- target_labels[[target]]
    tab
  })
)

heatmap_data <- subset(grid_table, feasible & period == generation_time)

p_heatmap <- ggplot(
  heatmap_data,
  aes(x = magnitude, y = duration, fill = population_reduction)
) +
  geom_tile() +
  facet_wrap(~ target, nrow = 2) +
  scale_fill_gradient(low = "white", high = "firebrick", limits = c(0, 100)) +
  labs(
    x = "Perturbation magnitude",
    y = "Perturbation duration (projection intervals)",
    fill = "Population\nreduction (%)",
    title = "Vulnerability surfaces at one event per reference generation"
  ) +
  theme_minimal(base_size = 11)

p_heatmap

## ----custom-target------------------------------------------------------------
custom_mask <- matrix(
  FALSE,
  nrow = nrow(A),
  ncol = ncol(A),
  dimnames = dimnames(A)
)
custom_mask["Adult", "Subadult"] <- TRUE
custom_mask["Adult", "Adult"] <- TRUE

custom_mask

A_custom_perturbed <- apply_perturbation(
  model,
  target = "custom",
  custom_mask = custom_mask,
  magnitude = 0.40
)

A_custom_perturbed

## ----custom-simulation--------------------------------------------------------
sim_custom <- simulate_dynamics(
  model,
  target = "custom",
  custom_mask = custom_mask,
  magnitude = 0.40,
  duration = 4,
  period = 10,
  t_max = 3 * generation_time,
  recovery_steps = 4 * generation_time,
  normalize_by_lambda = TRUE
)

sim_custom$reduction

## ----custom-grid--------------------------------------------------------------
out_custom <- run_grid(
  model,
  target = "custom",
  custom_mask = custom_mask,
  grid = grid,
  t_max = 3 * generation_time,
  recovery_steps = 4 * generation_time,
  normalize_by_lambda = TRUE
)

out_custom$vulnerability

## ----initial-state-example----------------------------------------------------
sim_initial_state <- simulate_dynamics(
  model,
  target = "adult_survival",
  magnitude = 0.30,
  duration = 4,
  period = 10,
  t_max = 3 * generation_time,
  recovery_steps = 4 * generation_time,
  initial_state = c(0.80, 0.15, 0.05),
  normalize_by_lambda = TRUE
)

sim_initial_state$reduction

## ----stage-vectors-example----------------------------------------------------
sim_stages <- simulate_dynamics(
  model,
  target = "juvenile_survival",
  magnitude = 0.30,
  duration = 4,
  period = 10,
  t_max = 3 * generation_time,
  recovery_steps = 4 * generation_time,
  return_stage_vectors = TRUE
)

head(sim_stages$stage_vectors)

## ----recovery-option-example--------------------------------------------------
sim_recovery_off <- simulate_dynamics(
  model,
  target = "adult_survival",
  magnitude = 0.30,
  duration = 4,
  period = 10,
  t_max = 3 * generation_time,
  recovery_steps = 4 * generation_time,
  force_during_recovery = FALSE
)

sim_recovery_on <- simulate_dynamics(
  model,
  target = "adult_survival",
  magnitude = 0.30,
  duration = 4,
  period = 10,
  t_max = 3 * generation_time,
  recovery_steps = 4 * generation_time,
  force_during_recovery = TRUE
)

c(
  perturbations_stop_before_recovery = sim_recovery_off$reduction,
  perturbations_continue_during_recovery = sim_recovery_on$reduction
)

## ----infeasible-example-------------------------------------------------------
grid_with_infeasible <- perturbation_grid(
  magnitudes = 0.5,
  durations = c(1, 4, 8),
  periods = c(2, 5)
)

out_with_infeasible <- run_grid(
  model,
  target = "adult_survival",
  grid = grid_with_infeasible,
  t_max = 30,
  skip_infeasible = FALSE
)

out_with_infeasible$table

## ----minimal-workflow, eval=FALSE---------------------------------------------
#  library(demovuln)
#  
#  A <- matrix(
#    c(
#      0.00, 0.00, 1.60,
#      0.35, 0.55, 0.05,
#      0.00, 0.25, 0.86
#    ),
#    nrow = 3,
#    byrow = TRUE
#  )
#  
#  model <- matrix_population_model(
#    A,
#    fecundity_rows = 1,
#    adult_stages = 3,
#    juvenile_stages = c(1, 2)
#  )
#  
#  generation_time <- 20
#  
#  grid <- perturbation_grid_from_frequencies(
#    magnitudes = seq(0, 0.8, by = 0.1),
#    durations = seq(0, generation_time, by = 2),
#    frequencies = c(0.25, 0.5, 1, 2, 4),
#    generation_time = generation_time
#  )
#  
#  out <- run_grid(
#    model,
#    target = "adult_survival",
#    grid = grid,
#    t_max = 3 * generation_time,
#    recovery_steps = 4 * generation_time
#  )
#  
#  out$vulnerability
#  head(out$table)

## ----session-info-------------------------------------------------------------
sessionInfo()