## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
library(aae.pop)
library(scales)

## -----------------------------------------------------------------------------
popmat_simple <- matrix(
  c(
    0, 2,       # reproduction of 10 new individuals per adult (female)
    0.45, 0.15  # move from class 1 to class 2 with 0.8 probability,
                # survive in class 2 with 0.25 probability
  ),
  nrow = 2,
  byrow = TRUE
)

## -----------------------------------------------------------------------------
# setup a matrix with four classes and all elements equal to zero
metapop_simple <- matrix(0, nrow = 4, ncol = 4)

# fill the first population's vital rates
#   (top two rows, two left-hand columns)
metapop_simple[1:2, 1:2] <- popmat_simple

# repeat for the second population, this time
#   filling the bottom right square
metapop_simple[3:4, 3:4] <- popmat_simple

## ----echo = FALSE-------------------------------------------------------------
metapop_simple

## ----dev = "png", fig.alt = "Line plot showing 100 simulated population trajectories for a metapopulation (combined abundances of all sub-populations shown). Lines show a population increasing from below 50 individuals to above 100 individuals, with minimal variation among lines due to differences in initial conditions."----
# allow 30 % of surviving adults to move between populations
#   (columns move to rows, so population 1 to 2 is in [4, 2] and
#    population 2 to 1 is in [2, 4])
metapop_simple[4, 2] <- 0.3 * metapop_simple[2, 2]
metapop_simple[2, 4] <- 0.3 * metapop_simple[4, 4]

# this 30 % needs to be removed from the adults that survive
#   and remain within the same population
metapop_simple[2, 2] <- metapop_simple[2, 2] - metapop_simple[4, 2]
metapop_simple[4, 4] <- metapop_simple[4, 4] - metapop_simple[2, 4]

# what happens if we simulate from this?
sims <- simulate(dynamics(metapop_simple), nsim = 100)
plot(sims, col = alpha("#2171B5", 0.4))

## ----dev = "png", fig.alt = "Line plots showing 100 simulated population trajectories for a metapopulation (sub-populations 1 and 2 shown separately from sub-populations 3 and 4). Lines show populations increasing from below 20 individuals to above 50 individuals, with minimal variation among lines due to differences in initial conditions."----
pop1_sims <- subset(sims, 1:2)
pop2_sims <- subset(sims, 3:4)
plot(pop1_sims, col = alpha("#2171B5", 0.4))
plot(pop2_sims, col = alpha("#2171B5", 0.4))

## ----echo = FALSE-------------------------------------------------------------
matrix(c(0, 1, 1, 0), nrow = 2)

## -----------------------------------------------------------------------------
structure_simple <- matrix(c(0, 1, 1, 0), nrow = 2)

## ----echo = FALSE-------------------------------------------------------------
x <- matrix(1, nrow = 5, ncol = 5)
diag(x) <- 0
x

## ----echo = FALSE-------------------------------------------------------------
x <- matrix(0, nrow = 5, ncol = 5)
x[1, 2] <- x[3, 2] <- x[5, 4] <- x[4, 1] <- x[3, 1] <- x[2, 4] <- 1
x

## -----------------------------------------------------------------------------
kern <- matrix(c(0, 0, 0, 0.3), nrow = 2)
dispersal_simple <- dispersal(kernel = kern, proportion = TRUE)

## -----------------------------------------------------------------------------
dispersal_simple <- list(dispersal_simple, dispersal_simple)

## -----------------------------------------------------------------------------
# create a population dynamics object from the population matrix
dynamics_simple <- dynamics(popmat_simple)

# turn this into a metapopulation
metapopulation_simple <- metapopulation(
  structure = structure_simple,
  dynamics = list(dynamics_simple, dynamics_simple), # both populations have identical rates
  dispersal = dispersal_simple
)

## -----------------------------------------------------------------------------
all.equal(metapop_simple, metapopulation_simple$matrix)

## ----dev = "png", fig.alt = "Line plots showing 100 simulated population trajectories for a metapopulation (sub-populations 1 and 2 shown separately from sub-populations 3 and 4). Lines show populations increasing from below 20 individuals to above 50 individuals, with minimal variation among lines due to differences in initial conditions."----
sims <- simulate(metapopulation_simple, nsim = 100)
plot(subset(sims, 1:2), col = alpha("#2171B5", 0.4), main = "Population 1")
plot(subset(sims, 3:4), col = alpha("#2171B5", 0.4), main = "Population 2")

## ----dev = "png", fig.alt = "Line plots showing 100 simulated population trajectories for a metapopulation (summed abundance of all sub-populations). Lines show populations increasing from approximately 200 individuals with rapid fluctuations in early years, stabilising to a slightly increasing population near 700 individuals at 50 years."----
# define population matrix
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
)

# define metapopulation structure
mp_str <- matrix(0, nrow = 5, ncol = 5)
mp_str[1, 2] <- mp_str[3, 2] <- mp_str[5, 4] <- mp_str[4, 1] <- mp_str[3, 1] <- mp_str[2, 4] <- 1

# define dispersal, assume all dispersal links are the same
mp_kern <- matrix(0, nrow = 5, ncol = 5)
mp_kern[3, 4] <- mp_kern[4, 5] <- 0.2  # adults can move once per year with 20 % probability
mp_disp <- dispersal(kernel = mp_kern, proportion = TRUE)

# define metapopulation
mp <- metapopulation(
  structure = mp_str,
  dynamics = lapply(1:5, function(x) dynamics(popmat)),   # five populations, all identical rates
  dispersal = lapply(1:6, function(x) mp_disp)  # six transitions in mp_str, all identical
)

# simulate
sims <- simulate(mp, nsim = 100)

# plot
plot(sims, col = alpha("#2171B5", 0.4), main = "All populations")

## ----dev = "png", fig.alt = "Line plots showing 100 simulated population trajectories for sub-populations 1 and 5 of the modelled metapopulation. Lines show population 1 rapidly increasing but then gradually declining back to its initial values at 50 years. Population 5 fluctuates strongly in early years but then gradually increases to 200 individuals at year 50."----
plot(subset(sims, 1:5), col = alpha("#2171B5", 0.4), main = "Population 1")
plot(subset(sims, 21:25), col = alpha("#2171B5", 0.4), main = "Population 5")

## ----dev = "png", fig.alt = "Line plots showing 100 simulated population trajectories for a metapopulation with covariate-driven dispersal. Lines show the total population fluctuating in early years but then decreasing to below 20 individuals at 50 years."----
# define two covariate sequences and combine them into a data.frame
covar_seq1 <- runif(20, min = 0.5, max = 1.0)
covar_seq2 <- runif(20, min = 0.75, max = 1.0)
covar_seq <- data.frame(pop1 = covar_seq1, pop2 = covar_seq2)

# define covariate responses for each population, pulling out
#   the appropriate covariate sequence in each
covar_fn1 <- function(mat, x) {
  mat * x$pop1
}
covar_fn2 <- function(mat, x) {
  mat * x$pop2
}

# define the covariates object for each population
covar1 <- covariates(masks = transition(popmat_simple), funs = covar_fn1)
covar2 <- covariates(masks = transition(popmat_simple), funs = covar_fn2)

# create a population dynamics object for each population
dynamics_simple1 <- dynamics(popmat_simple, covar1)
dynamics_simple2 <- dynamics(popmat_simple, covar2)

# turn this into a metapopulation
metapopulation_simple <- metapopulation(
  structure = structure_simple,
  dynamics = list(dynamics_simple1, dynamics_simple2),
  dispersal = dispersal_simple
)

# simulate
sims <- simulate(
  metapopulation_simple,
  nsim = 100,
  args = list(covariates = format_covariates(covar_seq))
)

# plot
plot(sims, col = alpha("#2171B5", 0.4), main = "Population-specific covariates")

## ----dev = "png", fig.alt = "Line plots showing 100 simulated population trajectories for a metapopulation with stochastic dispersal. Lines show the total population increasing steadily from 50 to approximately 150 individuals at 50 years."----
# define a mask that is TRUE whenever the dispersal kernel is non-zero
stoch_mask <- kern > 0

# define a stochastic dispersal function (small change in probability)
stoch_fun <- function(x) {
  rmultiunit(length(x), mean = x, sd = 0.05)
}

# using the kernel defined earlier (30 % probability of adults moving)
dispersal_simple <- dispersal(
  kernel = kern,
  stochasticity_masks = stoch_mask,
  stochasticity_funs = stoch_fun,
  proportion = TRUE
)

# dispersal is one-directional, so need one dispersal object for each population
dispersal_simple <- list(dispersal_simple, dispersal_simple)

# define metapopulation using the dynamics and structure objects defined above
metapopulation_simple <- metapopulation(
  structure = structure_simple,
  dynamics = list(dynamics_simple, dynamics_simple),
  dispersal = dispersal_simple
)

# simulate
sims <- simulate(metapopulation_simple, nsim = 100)

# plot
plot(sims, col = alpha("#2171B5", 0.4), main = "Stochastic dispersal")

## ----dev = "png", fig.alt = "Line plots showing 100 simulated population trajectories for a metapopulation with density-dependent dispersal. Lines show the total population fluctuating for approximately 20 years but then gradually increasing to 50-70 individuals at 50 years."----
# define a mask that is TRUE whenever the dispersal kernel is non-zero
#   (this is the same as for stochasticity and this mask 
#    could be used for both)
dd_mask <- kern > 0

# define Ricker type density dependence based on the abundance in the
#    destination population (recall n is a vector of two classes per
#    population, so elements 1 and 2 are population 1 and elements
#    3 and 4 are population 2)
dd_fun1 <- function(x, n) {
  # dispersing from population 1 to 2
  x * exp(1 - sum(n[3:4]) / 20) / exp(1) 
}
dd_fun2 <- function(x, n) {
  # dispersing from population 2 to 1
  x * exp(1 - sum(n[1:2]) / 20) / exp(1)
}

# using the kernel defined earlier (30 % of surviving adults moving),
#    add density masks and functions for each population
dispersal_simple1 <- dispersal(
  kernel = kern,
  density_masks = dd_mask,
  density_funs = dd_fun1,
  proportion = TRUE
)
dispersal_simple2 <- dispersal(
  kernel = kern,
  density_masks = dd_mask,
  density_funs = dd_fun2,
  proportion = TRUE
)

# combine the dispersal objects for both populations
dispersal_simple <- list(dispersal_simple1, dispersal_simple2)

# define metapopulation using the dynamics and structure objects defined above
metapopulation_simple <- metapopulation(
  structure = structure_simple,
  dynamics = list(dynamics_simple, dynamics_simple),
  dispersal = dispersal_simple
)

# simulate
sims <- simulate(metapopulation_simple, nsim = 100)

# plot
plot(sims, col = alpha("#2171B5", 0.4), main = "Density-dependent dispersal")

