## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
library(aae.pop)
library(scales)

# set up a basic population dynamics object to demonstrate examples below
popmat <- rbind(
  c(0,    0,    2,    4,    7),  # reproduction from 3-5 year olds
  c(0.25, 0,    0,    0,    0),  # survival from age 1 to 2
  c(0,    0.45, 0,    0,    0),  # survival from age 2 to 3
  c(0,    0,    0.70, 0,    0),  # survival from age 3 to 4
  c(0,    0,    0,    0.85, 0)   # survival from age 4 to 5
)
demostoch_mask <- all_classes(popmat) # affects all classes
demostoch_fn <- function(x) {
  rpois(length(x), lambda = x)
}
demostoch <- demographic_stochasticity(
  masks = demostoch_mask,
  funs = demostoch_fn
)
popdyn <- dynamics(popmat, demostoch)

## -----------------------------------------------------------------------------
sims <- simulate(
  popdyn,
  nsim = 100,
  seed = 123,
  options = list(ntime = 30, keep_slices = FALSE, tidy_abundances = floor)
)

## ----eval = FALSE-------------------------------------------------------------
#  # ntime
#  options(aae.pop_ntime = 30)
#  
#  # keep_slices
#  options(aae.pop_keep_slices = FALSE)
#  
#  # tidy_abundances
#  options(aae.pop_tidy_abundances = floor)

## -----------------------------------------------------------------------------
sims <- simulate(
  popdyn,
  options = list(initialise_args = 20)
)

## ----eval = FALSE-------------------------------------------------------------
#  options(aae.pop_lambda = 20)

## ----eval = FALSE-------------------------------------------------------------
#  # define your own initialisation function
#  my_initials_function <- function(n, other_arguments) {
#    simulate_n_values_somehow(n = n, other_arguments)  # e.g. rnbinom, rlnorm, ...
#  }
#  
#  # this neeeds wrapping up so aae.pop can use it
#  initials_function_wrapper <- function(n, args) {
#    do.call(my_initials_function, c(list(n), args))
#  }
#  
#  # set this at the global level (i.e. for an entire R session)
#  options(aae.pop_initialisation = initials_function_wrapper)

## ----dev = "png", fig.alt = "Line plot showing simulated trajectories through time with three different initial conditions. Larger initial population sizes correspond with higher population trajectories."----
# need a matrix/array with one row per replicate and one column
#    per class (5 age classes in this example)
my_initials <- matrix(
  c(
    100, 50, 25, 15, 5,   # first replicate
    200, 100, 50, 30, 10, # second replicate
    50, 25, 12, 8, 2      # third replicate
  ),
  nrow = 3,    # will set 3 replicates in simulate
  ncol = 5,    # and 5 age classes
  byrow = TRUE # this just makes formatting easier
)

# simulate with these values
sims <- simulate(
  popdyn,
  nsim = 3,
  seed = 123,
  init = my_initials,
  options = list(tidy_abundances = floor)
)

# plot this
plot(
  sims,
  col = alpha("#2171B5", 0.9),
  main = "Three different initial conditions"
)

## ----dev = "png", fig.alt = "Line plot showing 100 simulated trajectories initialised with the same initial conditions. Lines are variable around a central average."----
# need a vector with one value per class (5 age classes, here)
my_initials <- c(100, 50, 25, 15, 5)

# simulate with these values
sims <- simulate(
  popdyn,
  nsim = 100,
  seed = 123,
  init = my_initials,
  options = list(tidy_abundances = floor)
)

# plot this
plot(
  sims,
  col = alpha("#2171B5", 0.4),
  main = "Many replicates with the same initial condition"
)

## ----dev = "png", fig.alt = "Line plot showing 100 simulated trajectories initialised with stochastic initial conditions. Lines are variable around a central average with additional variation due to different initial conditions."----
# generate a matrix with one row per replicate and one column
#   per class (5 classes, 100 replicates), setting a different 
#   mean (lambda) in each column
my_initials <- matrix(
  rpois(5 * 100, lambda = c(100, 50, 25, 15, 5)),
  nrow = 100,
  ncol = 5,
  byrow = TRUE
)

# simulate with these values
sims <- simulate(
  popdyn,
  nsim = 100,
  seed = 123,
  init = my_initials,
  options = list(tidy_abundances = floor)
)

# plot this
plot(
  sims,
  col = alpha("#2171B5", 0.4),
  main = "Many replicates with stochastic, structured initial conditions"
)

## ----eval = FALSE-------------------------------------------------------------
#  # change updater
#  options(aae.pop_update = update_my_way)

## ----echo = FALSE-------------------------------------------------------------
update_crossprod

## ----echo = FALSE-------------------------------------------------------------
update_binomial_leslie