## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----packages-----------------------------------------------------------------
library(mpmaggregate)
library(kableExtra)
library(rphylopic)
library(expm)

## -----------------------------------------------------------------------------
# Leslie population projection matrix for Marmota flaviventris
# MatrixID = 249455

# Not run: requires Rcompadre and internet access

# Matrices can be retrieved with:
# library(Rcompadre)
# comadre <- cdb_fetch("comadre",version = "4.23.3.1")
# mpm <- comadre[comadre$MatrixID == 249455, ]
# matA <- matA(mpm)[[1]]
# matU <- matU(mpm)[[1]]
# matF <- matF(mpm)[[1]]


# Matrices are defined locally to avoid requiring the use of the internet and Rcompadre

# Leslie population projection matrix
matA <- matrix(c(
  0.000000, 0.258040, 0.627254, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037,
  0.545125, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
  0.000000, 0.498625, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
  0.000000, 0.000000, 0.592875, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
  0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
  0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000, 0.000000,
  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000, 0.000000,
  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000, 0.000000,
  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.643500, 0.000000
), nrow = 10, byrow = TRUE)


# Reproduction matrix
matF <- matrix(c(
  0, 0.25804, 0.627254, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037, 0.700037,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0,
  0, 0,       0,        0,        0,        0,        0,        0,        0,        0
), nrow = 10, byrow = TRUE)

# Survival/transition matrix
matU <- matrix(c(
  0,       0,       0,       0,      0,      0,      0,      0,      0,      0,
  0.545125,0,       0,       0,      0,      0,      0,      0,      0,      0,
  0,       0.498625,0,       0,      0,      0,      0,      0,      0,      0,
  0,       0,       0.592875,0,      0,      0,      0,      0,      0,      0,
  0,       0,       0,       0.6435, 0,      0,      0,      0,      0,      0,
  0,       0,       0,       0,      0.6435, 0,      0,      0,      0,      0,
  0,       0,       0,       0,      0,      0.6435, 0,      0,      0,      0,
  0,       0,       0,       0,      0,      0,      0.6435, 0,      0,      0,
  0,       0,       0,       0,      0,      0,      0,      0.6435, 0,      0,
  0,       0,       0,       0,      0,      0,      0,      0,      0.6435, 0
), nrow = 10, byrow = TRUE)

# Sanity check: A = U + F (after folding any C into F)
stopifnot(all.equal(unname(matA), unname(matU + matF)))
   

## ----helpers------------------------------------------------------------------
`%||%` <- function(x, y) if (!is.null(x)) x else y

# Return R0, the net reproductive rate
get_R0 <- function(U, F) {
  I <- diag(nrow(U))
  N <- solve(I - U)
  K <- F %*% N
  spectral_radius(K)
}

# Return generation time
get_Ta <- function(U, F) {
  generation_time(F, U, framework = "lambda")$generation_time
}

# Return cohort generation time
get_Tc <- function(U, F) {
  generation_time(F, U, framework = "R0")$generation_time
}

## ----aggregate----------------------------------------------------------------
# Aggregate into 5 x 5 matrices (m = 5)
ngroups <- 5

agg1 <- leslie_aggregate(
  matA = matA, ngroups = ngroups,
  framework = "lambda", criterion = "standard"
)

agg2 <- leslie_aggregate(
  matA = matA, ngroups = ngroups,
  framework = "lambda", criterion = "elasticity"
)

agg3 <- leslie_aggregate(
  matA = matA, ngroups = ngroups,
  framework = "R0", criterion = "standard"
)

agg4 <- leslie_aggregate(
  matA = matA, ngroups = ngroups,
  framework = "R0", criterion = "elasticity"
)

# Split a Leslie matrix into U and F components (Leslie form: fertility in row 1)
split_matA <- function(A) {
  n <- nrow(A)
  U <- A
  U[1, ] <- 0
  F <- matrix(0, n, n)
  F[1, ] <- A[1, ]
  list(U = U, F = F, A = A)
}

m0 <- list(U = matU, F = matF, A = matA)
m1 <- split_matA(agg1$matAk_agg)
m2 <- split_matA(agg2$matAk_agg)
m3 <- split_matA(agg3$matAk_agg)
m4 <- split_matA(agg4$matAk_agg)

eff <- c(
  NA,
  agg1$effectiveness,
  agg2$effectiveness,
  agg3$effectiveness,
  agg4$effectiveness
)

models <- list(
  Original            = m0,
  Agg_lambda_std      = m1,
  Agg_lambda_elast    = m2,
  Agg_R0_std          = m3,
  Agg_R0_elast        = m4
)

model_labels <- c(
  Original         = "Original",
  Agg_lambda_std   = "Agg &lambda; + standard",
  Agg_lambda_elast = "Agg &lambda; + elasticity",
  Agg_R0_std       = "Agg R<sub>0</sub> + standard",
  Agg_R0_elast     = "Agg R<sub>0</sub> + elasticity"
)

# Basic checks
stopifnot(all(sapply(models, function(m) !is.null(m$A))))

# Show aggregated matrices
print("Agg λ + standard");   m1$A
print("Agg λ + elasticity"); m2$A
print("Agg R0 + standard");  m3$A
print("Agg R0 + elasticity");m4$A

## ----results-table alt--------------------------------------------------------
# Adjust lambdas so that they are all per capita growth rate over a year
# Adjust generation times so they all use units of 1 year
adj = c(1,2,2,2,2)

lambda = sapply(models, function(m) spectral_radius(m$A))
lambda = lambda^(1/adj)
R0     = sapply(models, function(m) get_R0(m$U, m$F))
Ta     = sapply(models, function(m) get_Ta(m$U, m$F))*adj
Tc     = sapply(models, function(m) get_Tc(m$U, m$F))*adj

results <- data.frame(
  Model  = unname(model_labels[names(models)]),
  lambda = lambda,
  R0     = R0,
  Ta     = Ta,
  Tc     = Tc,
  Effectiveness = eff,
  row.names = NULL,
  stringsAsFactors = FALSE
)

n_rows <- nrow(results)

knitr::kable(
  results,
  col.names = c(
    "Model",
    "<em>&lambda;</em>",
    "<em>R</em><sub>0</sub>",
    "<em>T</em><sub>a</sub> (years)",
    "<em>T</em><sub>c</sub> (years)",
    "<em>&rho;</em><sup>2</sup>"
  ),
  digits = 3,
  caption = paste0(
    "<strong>Table 1.</strong> Comparison of demographic properties for the ",
    "yellow-bellied marmot <em>Marmota flaviventris</em> (COMADRE MatrixID 249455) ",
    "using the original 10 × 10 Leslie matrix and four aggregated Leslie matrices ",
    "collapsed to 5 age groups. Metrics include population growth rate ",
    "(<em>&lambda;</em>), net reproductive rate (<em>R</em><sub>0</sub>), ",
    "generation time (years; <em>T</em><sub>a</sub>), cohort generation time ",
    "(years; <em>T</em><sub>c</sub>), and aggregation effectiveness ",
    "(<em>&rho;</em><sup>2</sup>). &ldquo;Standard&rdquo; denotes use of a standard aggregator ",
    "and &ldquo;elasticity&rdquo; denotes elasticity-consistent aggregation."
  ),
  format = "html",
  escape = FALSE
) |>
  kable_styling(full_width = TRUE) |>
  row_spec(0, extra_css = "border-bottom: 2px solid black;") |>
  row_spec(nrow(results), extra_css = "border-bottom: 2px solid black;") |>
  footnote(
    general_title = "",
    general = paste0(
      "<span style='font-size: 90%;'>",
      "<strong>Note.</strong> Population growth rates (<em>&lambda;</em>) are all transformed to ",
      "represent growth over 1 year. Generation times have all been transformed to ",
      "represent generation time measured in years.</span>"
    ),
    footnote_as_chunk = TRUE,
    escape = FALSE
  )

## ----elasticities-lambda^k----------------------------------------------------

# Age groups used for aggregation in example
groups <- list(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))

# Partitioning matrix induced by the age groups.
# This matrix is constructed internally by leslie_aggregate(),
# but we construct it explicitly here so it can be used to
# map quantities from the original age-class space to the
# aggregated space.
P <- mpm_partition(groups = groups, n = nrow(matA))

# Dimensions of the original and aggregated matrices
m <- ngroups
n <- nrow(matA)
k <- n / m

matAk <- matA %^% k

# Elasticity of lambda^k with respect to entries of A^k
E_A <- mpm_elasticity(matA = matAk, framework = "lambda")$elasticity

E_list <- list(
  "Original"            = P %*% E_A %*% t(P),
  "Agg λ + standard"    = mpm_elasticity(matA = m1$A, framework = "lambda")$elasticity,
  "Agg λ + elasticity"  = mpm_elasticity(matA = m2$A, framework = "lambda")$elasticity,
  "Agg R0 + standard"   = mpm_elasticity(matA = m3$A, framework = "lambda")$elasticity,
  "Agg R0 + elasticity" = mpm_elasticity(matA = m4$A, framework = "lambda")$elasticity
)

to_long <- function(M, model_name) {
  idx <- which(M != 0, arr.ind = TRUE)
  data.frame(
    model = model_name,
    row = idx[, 1],
    col = idx[, 2],
    entry = paste(idx[, 1], idx[, 2], sep = ","),
    elasticity = M[idx],
    stringsAsFactors = FALSE
  )
}

elast_df <- do.call(rbind, Map(to_long, E_list, names(E_list)))
elast_df <- elast_df[elast_df$elasticity > 0, ]
elast_df <- elast_df[order(elast_df$row, elast_df$col), ]

models_order <- names(E_list)
entries <- unique(elast_df$entry)

elast_mat <- sapply(models_order, function(m) elast_df$elasticity[elast_df$model == m])
elast_mat <- t(elast_mat)
rownames(elast_mat) <- models_order
colnames(elast_mat) <- entries

## ----elasticity-plot, fig.width=12, fig.height=6------------------------------
# Not run: requires internet access

# Marmot 2012 Feb 2 by T. Michael Keesey, CC0 1.0, Marmota monax (Linnaeus 1758)
# uuid <- "eee50efb-40dc-47d0-b2cb-52a14a5e0e51"
# img <- get_phylopic(uuid = uuid, format = "raster", height = 869)
# Write image 
# png::writePNG(unclass(img), "inst/extdata/marmot.png")

# Load image locally
img_file <- system.file("extdata", "marmot.png", package = "mpmaggregate")

if (!nzchar(img_file)) {
  img_file <- "../inst/extdata/marmot.png"
}

img <- png::readPNG(img_file)

rc <- do.call(rbind, strsplit(colnames(elast_mat), ","))
r_vec <- as.integer(rc[, 1])
c_vec <- as.integer(rc[, 2])

make_entry_label <- function(r, c) {
  if (r == 1) {
    bquote(F[1, .(c)])
  } else {
    bquote(U[.(r), .(c)])
  }
}
entry_labels <- mapply(make_entry_label, r_vec, c_vec, SIMPLIFY = FALSE)

cols <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e")
op <- par(mar = c(7, 6, 4, 2))

bp <- barplot(
  elast_mat,
  beside = TRUE,
  ylim = c(0, .35),
  col = cols,
  border = NA,
  ylab = expression("Elasticity of " * lambda^k),
  xaxt = "n",
  space = c(0.2, 1)
)

axis(side = 1, at = colMeans(bp), labels = entry_labels, las = 2, cex.axis = 0.9)
legend("topleft", legend = rownames(elast_mat), fill = cols, bty = "n", inset = 0.02)

lims <- par("usr")
mydiff <- lims[4] - lims[3]
add_phylopic_base(
  img = img,
  x = NULL,
  y = mydiff * .9,
  height = .2 * mydiff * dim(img)[1] / 1024
)

par(op)

## ----elasticities-r0----------------------------------------------------------
# Age groupings used for aggregation (reintroduced here so this section
# can be read independently)
groups <- list(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))

# Partitioning matrix induced by the age groups.
# Constructed internally by leslie_aggregate(), but reproduced here
# because it is used to map quantities to aggregated space.
P <- mpm_partition(groups = groups, n = nrow(matA))

# Dimensions of the original and aggregated matrices
# Reintroduced here for clarity
m <- ngroups
n <- nrow(matA)
k <- n / m

R0 <- spectral_radius(matF %*% solve(diag(n) - matU))

# Reference operator for aggregation eigenstructure in the R0 framework
# Divide by R0 so that it has Leslie matrix form
A_R0ref <- (matF + R0 * matU) / R0

Uk <- expm::`%^%`(matU, k)
Rk <- (expm::`%^%`(A_R0ref, k) - Uk) * R0
Ak <- Uk + Rk

# Elasticity of R0 with respect to entries of Ak (closest to what aggregated matrices yield)
E_A <- mpm_elasticity(matF = Rk, matU = Uk, framework = "R0", normalize = TRUE)$elasticity

E_list2 <- list(
  "Original"            = P %*% E_A %*% t(P),
  "Agg λ + standard"    = mpm_elasticity(matF = m1$F, matU = m1$U, framework = "R0", normalize = TRUE)$elasticity,
  "Agg λ + elasticity"  = mpm_elasticity(matF = m2$F, matU = m2$U, framework = "R0", normalize = TRUE)$elasticity,
  "Agg R0 + standard"   = mpm_elasticity(matF = m3$F, matU = m3$U, framework = "R0", normalize = TRUE)$elasticity,
  "Agg R0 + elasticity" = mpm_elasticity(matF = m4$F, matU = m4$U, framework = "R0", normalize = TRUE)$elasticity
)

elast_df2 <- do.call(rbind, Map(to_long, E_list2, names(E_list2)))
elast_df2 <- elast_df2[elast_df2$elasticity > 0, ]
elast_df2 <- elast_df2[order(elast_df2$row, elast_df2$col), ]

models_order2 <- names(E_list2)
entries2 <- unique(elast_df2$entry)

elast_mat2 <- sapply(models_order2, function(m) elast_df2$elasticity[elast_df2$model == m])
elast_mat2 <- t(elast_mat2)
rownames(elast_mat2) <- models_order2
colnames(elast_mat2) <- entries2

## ----elasticity-r0-plot, fig.width=12, fig.height=6---------------------------

# Not run: requires internet access

# Marmot 2012 Feb 2 by T. Michael Keesey, CC0 1.0, Marmota monax (Linnaeus 1758)
# uuid <- "eee50efb-40dc-47d0-b2cb-52a14a5e0e51"
# img <- get_phylopic(uuid = uuid, format = "raster", height = 869)
# Write image 
# png::writePNG(unclass(img), "inst/extdata/marmot.png")

# Load the image file locally
img_file <- system.file("extdata", "marmot.png", package = "mpmaggregate")

#For development
if (!nzchar(img_file)) {
  img_file <- "../inst/extdata/marmot.png"
}

img <- png::readPNG(img_file)

rc <- do.call(rbind, strsplit(colnames(elast_mat2), ","))
r_vec <- as.integer(rc[, 1])
c_vec <- as.integer(rc[, 2])

entry_labels <- mapply(make_entry_label, r_vec, c_vec, SIMPLIFY = FALSE)

cols <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e")
op <- par(mar = c(7, 6, 4, 2))

bp <- barplot(
  elast_mat2,
  beside = TRUE,
  ylim = c(0, .35),
  col = cols,
  border = NA,
  ylab = expression("Normalized elasticity of " * R[0]),
  xaxt = "n",
  space = c(0.2, 1)
)

axis(side = 1, at = colMeans(bp), labels = entry_labels, las = 2, cex.axis = 0.9)
legend("topleft", legend = rownames(elast_mat2), fill = cols, bty = "n", inset = 0.02)

lims <- par("usr")
mydiff <- lims[4] - lims[3]
add_phylopic_base(
  img = img,
  x = NULL,
  y = mydiff * .9,
  height = .2 * mydiff * dim(img)[1] / 1024
)

par(op)

## ----fertility-elasticity-mass------------------------------------------------
# Mass for elasticity of R0 w.r.t. fertilities
sum(elast_mat2[, 1:5]) / 5

# Mass for elasticity of lambda^k w.r.t. fertilities
sum(elast_mat[, 1:5]) / 5