library(magrittr)
ling_imm <- g3_stock('ling_imm', seq(20, 156, 4)) %>% g3s_age(3, 10)
ling_mat <- g3_stock('ling_mat', seq(20, 156, 4)) %>% g3s_age(5, 15)
lln <- g3_fleet('lln')
# Invent a ldist.lln table for our tests
ldist.lln.raw <- data.frame(
year = c(1999, 2000),
age = sample(5:9, 100, replace = TRUE),
length = sample(10:70, 100, replace = TRUE),
number = 1,
stringsAsFactors = FALSE)
# Group length into 10-long bins
# NB: The last 2 bins will be empty, but gadget3 will use the factor levels, include them as zero
# NB: Generally one would use mfdb::mfdb_sample_count() source and group data for you
ldist.lln.raw |> dplyr::group_by(
year = year, age = age,
length = cut(length, breaks = seq(10, 100, by = 10), right = FALSE)
) |> dplyr::summarise(number = sum(number), .groups = 'keep') -> ldist.lln
# Turn age into a factor, indicating all ages we should be interested in
ldist.lln$age <- factor(ldist.lln$age, levels = 5:15)
# We can see the results of this being turned into an array:
g3_distribution_preview(ldist.lln)
likelihood_actions <- list(
g3l_catchdistribution(
'ldist_lln',
ldist.lln,
fleets = list(lln),
stocks = list(ling_imm, ling_mat),
g3l_distribution_sumofsquares()))
# Make an (incomplete) model using the action, extract the observation array
fn <- suppressWarnings(g3_to_r(likelihood_actions))
environment(fn)$cdist_sumofsquares_ldist_lln_obs__num
# Apply age-reading error matrix to model data
more_likelihood_actions <- list(
g3l_catchdistribution(
'ldist_lln_readerror',
ldist.lln,
fleets = list(lln),
stocks = list(ling_imm, ling_mat),
transform_fs = list(age = g3_formula(
g3_param_array('reader1matrix', value = diag(5))[g3_idx(preage), g3_idx(age)]
)),
g3l_distribution_sumofsquares()))
# Apply per-stock age-reading error matrix to model data
more_likelihood_actions <- list(
g3l_catchdistribution(
'ldist_lln_readerror',
ldist.lln,
fleets = list(lln),
stocks = list(ling_imm, ling_mat),
transform_fs = list(age = list(
ling_imm = quote( g3_param_array('imm_readermatrix',
value = diag(ling_imm__maxage - ling_imm__minage + 1)
)[ling_imm__preage_idx, ling_imm__age_idx] ),
ling_mat = quote( g3_param_array('mat_readermatrix',
value = diag(ling_mat__maxage - ling_mat__minage + 1)
)[ling_mat__preage_idx, ling_mat__age_idx] ),
unused = 0)),
g3l_distribution_sumofsquares()))
## Stomach content: predator-prey species preference
prey_a <- g3_stock('prey_a', seq(1, 10)) |> g3s_age(1,3)
prey_b <- g3_stock('prey_b', seq(1, 10)) |> g3s_age(1,3)
pred_a <- g3_stock('pred_a', seq(50, 80, by = 10)) |> g3s_age(0, 10)
otherfood <- g3_stock('otherfood', 0)
# Produce data.frame with columns:
# * predator_length or predator_age
# * stock
# * number or weight
pred_a_preypref_obs <- expand.grid(
year = 2000:2005,
predator_length = c(50,70),
stock = c('prey_a', 'prey_b', 'otherfood'),
number = 0 )
# Create catchdistribution likelihood component
actions <- list(
g3l_catchdistribution(
'pred_a_preypref',
pred_a_preypref_obs,
fleets = list(pred_a),
stocks = list(prey_a, prey_b, otherfood),
g3l_distribution_sumofsquares(),
nll_breakdown = TRUE,
report = TRUE ),
NULL)
## Stomach content: predator-prey size preference
# Produce data.frame with columns:
# * predator_length or predator_age
# * (prey) length
# * number or weight
pred_a_sizepref_obs <- expand.grid(
year = 2000:2005,
predator_length = c(50,70),
length = seq(1, 10),
number = 0 )
# Create catchdistribution likelihood component
actions <- list(
g3l_catchdistribution(
'pred_a_sizepref',
pred_a_sizepref_obs,
predators = list(pred_a),
# NB: Only referencing stocks included in observation data
stocks = list(prey_a),
function_f = g3l_distribution_sumofsquares(),
# Use transform_fs to apply digestioncoefficients
transform_fs = list(length = list(prey_a = g3_formula(
quote( diag(d0 + d1 * prey_a__midlen^d2) ),
d0 = g3_parameterized('d0', by_stock = TRUE),
d1 = g3_parameterized('d1', by_stock = TRUE),
d2 = g3_parameterized('d2', by_stock = TRUE) ))),
nll_breakdown = TRUE,
report = TRUE ),
NULL)
Run the code above in your browser using DataLab