library(demovuln)

A <- matrix(
  c(0.0, 0.4,
    2.0, 0.7),
  nrow = 2,
  byrow = FALSE
)

model <- matrix_population_model(
  A,
  adult_stages = 2,
  juvenile_stages = 1
)

grid <- perturbation_grid(
  magnitudes = seq(0, 1, length.out = 11),
  durations = 0:3,
  periods = c(1, 2, 3, 5, 10)
)

out <- run_grid(
  model,
  target = "adult_survival",
  grid = grid,
  t_max = 50,
  recovery_steps = 10
)

out$vulnerability
#> [1] 61.26872
head(out$table)
#>           target magnitude duration period feasible population_reduction
#> 1 adult_survival       0.0        0      1     TRUE                    0
#> 2 adult_survival       0.1        0      1     TRUE                    0
#> 3 adult_survival       0.2        0      1     TRUE                    0
#> 4 adult_survival       0.3        0      1     TRUE                    0
#> 5 adult_survival       0.4        0      1     TRUE                    0
#> 6 adult_survival       0.5        0      1     TRUE                    0
#>   final_population baseline_final_population
#> 1                1                         1
#> 2                1                         1
#> 3                1                         1
#> 4                1                         1
#> 5                1                         1
#> 6                1                         1

tab <- subset(out$table, duration == 1 & feasible)

mat <- tapply(
  tab$population_reduction,
  list(tab$magnitude, tab$period),
  mean
)

image(
  x = as.numeric(rownames(mat)),
  y = as.numeric(colnames(mat)),
  z = mat,
  xlab = "Magnitude",
  ylab = "Period",
  main = "Population reduction"
)

contour(
  x = as.numeric(rownames(mat)),
  y = as.numeric(colnames(mat)),
  z = mat,
  add = TRUE
)