The following describes the structure of a gadget3 model, from the bottom up.
Crucial to gadget3 is the R
formula, created using the tilde operator (~). These
get used in several places in R for various things, but at their core
the tilde operator stores the R code on the left and right hand sides,
as well as the environment it was created in, which amounts to all
variables defined at the time.
For example, let’s declare a function that produces a formula, and make some:
get_formula <- function (size) {
# NB: The reason we make a function here is so we have an isolated environment
# to make examples cleaner.
cows <- size * 2
pigs <- size * 4
return(~cows + pigs)
}
f <- get_formula(8)
g <- get_formula(5)str() shows f contains the formula’s code
(not the result of cows + pigs), and has an environment attached:
str(f)## Class 'formula' language ~cows + pigs
## ..- attr(*, ".Environment")=<environment: 0x55c0aead6fb8>We can look into this environment, and see the values that got set for cows & pigs:
str(as.list(environment(f)))## List of 3
## $ pigs: num 32
## $ cows: num 16
## $ size: num 8We can do similarly for g, and see the results are
different:
str(as.list(environment(g)))## List of 3
## $ pigs: num 20
## $ cows: num 10
## $ size: num 5Note that:
cows + pigs, the code
was stored for later use.A g3 model at it’s core is a list of formula objects that make up the model. We can even use Gadget3 to compile our simple example above into an R function:
## function (param = parameter_template)
## {
## if (is.data.frame(param)) {
## param_lower <- structure(param$lower, names = param$switch)
## param_upper <- structure(param$upper, names = param$switch)
## param <- structure(param$value, names = param$switch)
## }
## else {
## param_lower <- lapply(param, function(x) NA)
## param_upper <- lapply(param, function(x) NA)
## }
## cows <- 16
## pigs <- 32
## while (TRUE) {
## cows + pigs
## }
## }
## <bytecode: 0x55c0ab650d10>
## <environment: 0x55c0b00368b8>…or a TMB template:
## #include <TMB.hpp>
##
##
## #ifndef TYPE_IS_SCALAR
## #ifdef TMBAD_FRAMEWORK
## #define TYPE_IS_SCALAR(TestT) typename = std::enable_if_t<std::is_same<TestT, int>::value || std::is_same<TestT, double>::value || std::is_same<TestT, TMBad::global::ad_aug>::value>
## #endif // TMBAD_FRAMEWORK
## #ifdef CPPAD_FRAMEWORK
## #define TYPE_IS_SCALAR(TestT) typename = std::enable_if_t<std::is_same<TestT, int>::value || std::is_same<TestT, double>::value || std::is_same<TestT, CppAD::AD<double>>::value || std::is_same<TestT, CppAD::AD<CppAD::AD<double>>>::value || std::is_same<TestT, CppAD::AD<CppAD::AD<CppAD::AD<double>>>>::value>
## #endif // CPPAD_FRAMEWORK
## #endif // TYPE_IS_SCALAR
##
##
##
## template<class Type>
## Type objective_function<Type>::operator() () {
## DATA_SCALAR(reporting_enabled); DATA_UPDATE(reporting_enabled);
## Type cows = (double)(16);
## Type pigs = (double)(32);
##
## while (true) {
## cows + pigs;
## }
## }Obviously this in itself isn’t a very useful function, but you can see that the environment for the formula has provided the initial values for the variables, and the code itself has been put into the main loop.
Also note that, despite appearances, we’re not generically converting R into C++. There’s a subset of R that gadget3 understands and knows how to convert. Using R libraries isn’t possible, for example. To simplify scoping rules, we assume variables are either global, and should have different names, or are iterators in loops, in which case they will be local to that loop.
In reality you will never be providing formulae to insert into models
directly, you’d be using the g3 action functions to generate these for
you. All action functions, prefixed with g3a_, produce a
list of formula objects—an action in gadget3 parlance. These are where
the gadget functionality is implemented.
One of the simplest is g3a_time, which produces code
that will count years/steps, and stop when the end of the time period is
reached. For example:
g3a_time(1990, 1999)## $`-01`
## ~{
## debug_label("g3a_time: Start of time period")
## cur_time <- cur_time + 1L
## if (cur_time == 0 && assert_msg(g3_param("retro_years", value = 0,
## optimise = FALSE, source = "g3a_time") >= 0, "retro_years must be >= 0"))
## return(NaN)
## if (cur_time == 0 && assert_msg(project_years >= 0, "project_years must be >= 0"))
## return(NaN)
## cur_year <- start_year + (cur_time%/%step_count)
## cur_year_projection <- cur_year > end_year - g3_param("retro_years",
## value = 0, optimise = FALSE, source = "g3a_time")
## cur_step <- (cur_time%%step_count) + 1L
## cur_step_size <- step_lengths[[cur_step]]/12
## cur_step_final <- cur_step == step_count
## }
## <environment: 0x55c0ad5f76c0>
##
## $`001`
## ~{
## if (cur_time > total_steps)
## return(nll)
## }
## <environment: 0x55c0aa35f678>Like in the example above, the definitions are part of the formula’s environment, and if we compile it we see our years ending up in the code.
## function (param = parameter_template)
## {
## if (is.data.frame(param)) {
## param_lower <- structure(param$lower, names = param$switch)
## param_upper <- structure(param$upper, names = param$switch)
## param <- structure(param$value, names = param$switch)
## }
## else {
## param_lower <- lapply(param, function(x) NA)
## param_upper <- lapply(param, function(x) NA)
## }
## stopifnot("retro_years" %in% names(param))
## assert_msg <- function(expr, message) {
## if (isFALSE(expr)) {
## warning(message)
## return(TRUE)
## }
## return(FALSE)
## }
## cur_time <- -1L
## stopifnot("project_years" %in% names(param))
## project_years <- param[["project_years"]]
## cur_year <- 0L
## start_year <- 1990L
## step_lengths <- 12L
## step_count <- length(step_lengths)
## cur_year_projection <- FALSE
## end_year <- 1999L
## cur_step <- 0L
## cur_step_size <- step_lengths[[1]]/12
## cur_step_final <- FALSE
## as_integer <- as.integer
## retro_years <- param[["retro_years"]]
## total_steps <- length(step_lengths) * (end_year - as_integer(retro_years) -
## start_year + as_integer(project_years)) + length(step_lengths) -
## 1L
## nll <- 0
## while (TRUE) {
## {
## comment("g3a_time: Start of time period")
## cur_time <- cur_time + 1L
## if (cur_time == 0 && assert_msg(param[["retro_years"]] >=
## 0, "retro_years must be >= 0"))
## return(NaN)
## if (cur_time == 0 && assert_msg(project_years >=
## 0, "project_years must be >= 0"))
## return(NaN)
## cur_year <- start_year + (cur_time%/%step_count)
## cur_year_projection <- cur_year > end_year - param[["retro_years"]]
## cur_step <- (cur_time%%step_count) + 1L
## cur_step_size <- step_lengths[[cur_step]]/12
## cur_step_final <- cur_step == step_count
## }
## {
## if (cur_time > total_steps)
## return(nll)
## }
## }
## }
## <bytecode: 0x55c0af5f0ed0>
## <environment: 0x55c0ad104ec8>In this case our years have been hard-coded, but their definitions could themselves be formula and the end result will be added to the model. For example:
## function (param = parameter_template)
## {
## if (is.data.frame(param)) {
## param_lower <- structure(param$lower, names = param$switch)
## param_upper <- structure(param$upper, names = param$switch)
## param <- structure(param$value, names = param$switch)
## }
## else {
## param_lower <- lapply(param, function(x) NA)
## param_upper <- lapply(param, function(x) NA)
## }
## stopifnot("retro_years" %in% names(param))
## assert_msg <- function(expr, message) {
## if (isFALSE(expr)) {
## warning(message)
## return(TRUE)
## }
## return(FALSE)
## }
## cur_time <- -1L
## stopifnot("project_years" %in% names(param))
## project_years <- param[["project_years"]]
## cur_year <- 0L
## start_year <- 1990L
## step_lengths <- 12L
## step_count <- length(step_lengths)
## cur_year_projection <- FALSE
## end_year <- start_year + 4
## cur_step <- 0L
## cur_step_size <- step_lengths[[1]]/12
## cur_step_final <- FALSE
## as_integer <- as.integer
## retro_years <- param[["retro_years"]]
## total_steps <- length(step_lengths) * (end_year - as_integer(retro_years) -
## start_year + as_integer(project_years)) + length(step_lengths) -
## 1L
## nll <- 0
## while (TRUE) {
## {
## comment("g3a_time: Start of time period")
## cur_time <- cur_time + 1L
## if (cur_time == 0 && assert_msg(param[["retro_years"]] >=
## 0, "retro_years must be >= 0"))
## return(NaN)
## if (cur_time == 0 && assert_msg(project_years >=
## 0, "project_years must be >= 0"))
## return(NaN)
## cur_year <- start_year + (cur_time%/%step_count)
## cur_year_projection <- cur_year > end_year - param[["retro_years"]]
## cur_step <- (cur_time%%step_count) + 1L
## cur_step_size <- step_lengths[[cur_step]]/12
## cur_step_final <- cur_step == step_count
## }
## {
## if (cur_time > total_steps)
## return(nll)
## }
## }
## }
## <bytecode: 0x55c0ae13cde0>
## <environment: 0x55c0b00a91a8>Beyond g3a_time(), pretty much any action will be
describing changes to a stock, possibly via. interacting with another
stock (or fleet). To keep track of this state, we use g3_stock
objects. These describe several things:
g3_stock objects can be created with either g3_stock()
or g3_fleet(), the former will store lengthgroups, which
fleets do not have.
Actions will store data about the stocks, for instance the current
number of individuals, in arrays called stock instances. We can see what
sort of arrays a stock will make by using
g3_stock_instance:
ling_imm <- g3_stock('ling_imm', seq(0, 50, 10))
g3_stock_instance(ling_imm)## length
## 0:10 10:20 20:30 30:40 40:50 50:Inf
## NA NA NA NA NA NATo add complexity to our model, we can use other g3s_
functions, such as g3s_livesonareas() or
g3s_age(), which adds area or age dimensions to a
stock:
ling_imm <- g3_stock('ling_imm', seq(0, 50, 10)) %>%
g3s_age(3, 10)
g3_stock_instance(ling_imm)## age
## length age3 age4 age5 age6 age7 age8 age9 age10
## 0:10 NA NA NA NA NA NA NA NA
## 10:20 NA NA NA NA NA NA NA NA
## 20:30 NA NA NA NA NA NA NA NA
## 30:40 NA NA NA NA NA NA NA NA
## 40:50 NA NA NA NA NA NA NA NA
## 50:Inf NA NA NA NA NA NA NA NA
ling_imm <- g3_stock('ling_imm', seq(0, 50, 10)) %>%
g3s_livesonareas(c(1,2)) %>%
g3s_age(3, 10)
g3_stock_instance(ling_imm)[,,'age3']## area
## length area1 area2
## 0:10 NA NA
## 10:20 NA NA
## 20:30 NA NA
## 30:40 NA NA
## 40:50 NA NA
## 50:Inf NA NAWhen you use an action such as g3a_growmature(), you
provide g3_stock(s) to act on, and formula objects to fill in gaps in
the code. g3a_growmature() will iterate over all areas/ages
the stock has, and apply growth for each length group it finds.
g3a_growmature() itself doesn’t care about areas or age,
it just does the same to each. However, the formula you supply can.
age and area variables will be set with the
current age/area, which you can use when writing formula. For
example:
fn <- g3_to_r(g3a_growmature(
ling_imm,
impl_f = g3a_grow_impl_bbinom(
delta_len_f = ~age * 10,
delta_wgt_f = ~area * 20,
beta_f = ~g3_param("ling.bbin"),
maxlengthgroupgrowth = 4),
transition_f = ~TRUE))
fn## function (param = parameter_template)
## {
## if (is.data.frame(param)) {
## param_lower <- structure(param$lower, names = param$switch)
## param_upper <- structure(param$upper, names = param$switch)
## param <- structure(param$value, names = param$switch)
## }
## else {
## param_lower <- lapply(param, function(x) NA)
## param_upper <- lapply(param, function(x) NA)
## }
## stopifnot("ling.bbin" %in% names(param))
## growth_bbinom <- function(delt_l, binn, beta) {
## alpha <- (beta * delt_l)/(binn - delt_l)
## x <- 0:binn
## na <- length(alpha)
## n <- length(x) - 1
## alpha <- rep(alpha, n + 1)
## x <- rep(x, each = na)
## val <- exp(lgamma(n + 1) + lgamma(alpha + beta) + lgamma(n -
## x + beta) + lgamma(x + alpha) - lgamma(n - x + 1) -
## lgamma(x + 1) - lgamma(n + alpha + beta) - lgamma(beta) -
## lgamma(alpha))
## dim(val) <- c(na, n + 1)
## return(val)
## }
## dif_pmax <- function(a, b, scale) {
## logspace_add <- function(a, b) pmax(a, b) + log1p(exp(pmin(a,
## b) - pmax(a, b)))
## b <- as.vector(b)
## logspace_add(a * scale, b * scale)/scale
## }
## avoid_zero <- function(a) {
## dif_pmax(a, 0, 1000)
## }
## g3a_grow_matrix_wgt <- function(delta_w) {
## na <- dim(delta_w)[[1]]
## n <- dim(delta_w)[[2]] - 1
## wgt.matrix <- array(0, c(na, na))
## for (lg in 1:na) {
## if (lg == na) {
## wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na -
## lg + 1)]
## }
## else if (lg + n > na) {
## wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na -
## lg + 1)]
## }
## else {
## wgt.matrix[lg, lg:(n + lg)] <- delta_w[lg, ]
## }
## }
## return(wgt.matrix)
## }
## g3a_grow_matrix_len <- function(delta_l) {
## na <- dim(delta_l)[[1]]
## n <- dim(delta_l)[[2]] - 1
## growth.matrix <- array(0, c(na, na))
## for (lg in 1:na) {
## if (lg == na) {
## growth.matrix[na, na] <- sum(delta_l[lg, ])
## }
## else if (lg + n > na) {
## growth.matrix[lg, lg:(na - 1)] <- delta_l[lg,
## 1:(na - lg)]
## growth.matrix[lg, na] <- sum(delta_l[lg, (na -
## lg + 1):(n + 1)])
## }
## else {
## growth.matrix[lg, lg:(n + lg)] <- delta_l[lg,
## ]
## }
## }
## return(growth.matrix)
## }
## g3a_grow_apply <- function(growth.matrix, wgt.matrix, input_num,
## input_wgt) {
## na <- dim(growth.matrix)[[1]]
## growth.matrix <- growth.matrix * as.vector(input_num)
## wgt.matrix <- growth.matrix * (wgt.matrix + as.vector(input_wgt))
## growth.matrix.sum <- colSums(growth.matrix)
## return(array(c(growth.matrix.sum, colSums(wgt.matrix)/avoid_zero(growth.matrix.sum)),
## dim = c(na, 2)))
## }
## assert_msg <- function(expr, message) {
## if (isFALSE(expr)) {
## warning(message)
## return(TRUE)
## }
## return(FALSE)
## }
## ling_imm__minage <- 3L
## ling_imm__maxage <- 10L
## ling_imm__growth_l <- array(NA, dim = structure(6:5, names = c("length",
## "delta")), dimnames = list(length = c("0:10", "10:20",
## "20:30", "30:40", "40:50", "50:Inf"), delta = c("0",
## "1", "2", "3", "4")))
## ling_imm__plusdl <- 10
## ling_imm__growth_w <- array(NA, dim = structure(6:5, names = c("length",
## "delta")), dimnames = list(length = c("0:10", "10:20",
## "20:30", "30:40", "40:50", "50:Inf"), delta = c("0",
## "1", "2", "3", "4")))
## ling_imm__num <- array(0, dim = c(length = 6L, area = 2L,
## age = 8L), dimnames = list(length = c("0:10", "10:20",
## "20:30", "30:40", "40:50", "50:Inf"), area = c("area1",
## "area2"), age = c("age3", "age4", "age5", "age6", "age7",
## "age8", "age9", "age10")))
## ling_imm__wgt <- array(1, dim = c(length = 6L, area = 2L,
## age = 8L), dimnames = list(length = c("0:10", "10:20",
## "20:30", "30:40", "40:50", "50:Inf"), area = c("area1",
## "area2"), age = c("age3", "age4", "age5", "age6", "age7",
## "age8", "age9", "age10")))
## ling_imm__prevtotal <- 0
## while (TRUE) {
## {
## comment("g3a_grow for ling_imm")
## for (age in seq(ling_imm__minage, ling_imm__maxage,
## by = 1)) {
## ling_imm__age_idx <- age - ling_imm__minage +
## 1L
## for (ling_imm__area_idx in seq_along(ling_imm__areas)) {
## area <- ling_imm__areas[[ling_imm__area_idx]]
## growth_delta_l <- (ling_imm__growth_l[] <- growth_bbinom(avoid_zero((age *
## 10)/ling_imm__plusdl), 4L, avoid_zero(param[["ling.bbin"]])))
## growth_delta_w <- (ling_imm__growth_w[] <- area *
## 20)
## growthmat_w <- g3a_grow_matrix_wgt(growth_delta_w)
## growthmat_l <- g3a_grow_matrix_len(growth_delta_l)
## growthresult <- g3a_grow_apply(growthmat_l,
## growthmat_w, ling_imm__num[, ling_imm__area_idx,
## ling_imm__age_idx], ling_imm__wgt[, ling_imm__area_idx,
## ling_imm__age_idx])
## {
## if (FALSE)
## ling_imm__prevtotal <- sum(ling_imm__num[,
## ling_imm__area_idx, ling_imm__age_idx])
## comment("Update ling_imm using delta matrices")
## ling_imm__num[, ling_imm__area_idx, ling_imm__age_idx] <- growthresult[,
## (1)]
## ling_imm__wgt[, ling_imm__area_idx, ling_imm__age_idx] <- growthresult[,
## (2)]
## if (FALSE)
## assert_msg(~abs(ling_imm__prevtotal - sum(ling_imm__num[,
## ling_imm__area_idx, ling_imm__age_idx])) <
## 1e-04, "g3a_growmature: ling_imm__num totals are not the same before and after growth")
## }
## }
## }
## }
## }
## }
## <bytecode: 0x55c0ae8905b8>
## <environment: 0x55c0b056a0b8>…you can see our provided formula have been used to calculate
ling_imm__growth_l and ling_imm__growth_w, and
age and area are available to us thanks to the
loops provided by the stock.
We can also define formulas and reference them. Gadget3 will add the definitions into the code where appropriate. For example:
custom_delta_l <- ~area * 99
custom_delta_w <- ~ling_imm__plusdl * 44
fn <- g3_to_r(g3a_growmature(
ling_imm,
impl_f = g3a_grow_impl_bbinom(
delta_len_f = ~age * custom_delta_l * 10,
delta_wgt_f = ~area * custom_delta_w * 20,
beta_f = ~g3_param("ling.bbin"),
maxlengthgroupgrowth = 4),
transition_f = ~TRUE))
fn## function (param = parameter_template)
## {
## if (is.data.frame(param)) {
## param_lower <- structure(param$lower, names = param$switch)
## param_upper <- structure(param$upper, names = param$switch)
## param <- structure(param$value, names = param$switch)
## }
## else {
## param_lower <- lapply(param, function(x) NA)
## param_upper <- lapply(param, function(x) NA)
## }
## stopifnot("ling.bbin" %in% names(param))
## growth_bbinom <- function(delt_l, binn, beta) {
## alpha <- (beta * delt_l)/(binn - delt_l)
## x <- 0:binn
## na <- length(alpha)
## n <- length(x) - 1
## alpha <- rep(alpha, n + 1)
## x <- rep(x, each = na)
## val <- exp(lgamma(n + 1) + lgamma(alpha + beta) + lgamma(n -
## x + beta) + lgamma(x + alpha) - lgamma(n - x + 1) -
## lgamma(x + 1) - lgamma(n + alpha + beta) - lgamma(beta) -
## lgamma(alpha))
## dim(val) <- c(na, n + 1)
## return(val)
## }
## dif_pmax <- function(a, b, scale) {
## logspace_add <- function(a, b) pmax(a, b) + log1p(exp(pmin(a,
## b) - pmax(a, b)))
## b <- as.vector(b)
## logspace_add(a * scale, b * scale)/scale
## }
## avoid_zero <- function(a) {
## dif_pmax(a, 0, 1000)
## }
## g3a_grow_matrix_wgt <- function(delta_w) {
## na <- dim(delta_w)[[1]]
## n <- dim(delta_w)[[2]] - 1
## wgt.matrix <- array(0, c(na, na))
## for (lg in 1:na) {
## if (lg == na) {
## wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na -
## lg + 1)]
## }
## else if (lg + n > na) {
## wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na -
## lg + 1)]
## }
## else {
## wgt.matrix[lg, lg:(n + lg)] <- delta_w[lg, ]
## }
## }
## return(wgt.matrix)
## }
## g3a_grow_matrix_len <- function(delta_l) {
## na <- dim(delta_l)[[1]]
## n <- dim(delta_l)[[2]] - 1
## growth.matrix <- array(0, c(na, na))
## for (lg in 1:na) {
## if (lg == na) {
## growth.matrix[na, na] <- sum(delta_l[lg, ])
## }
## else if (lg + n > na) {
## growth.matrix[lg, lg:(na - 1)] <- delta_l[lg,
## 1:(na - lg)]
## growth.matrix[lg, na] <- sum(delta_l[lg, (na -
## lg + 1):(n + 1)])
## }
## else {
## growth.matrix[lg, lg:(n + lg)] <- delta_l[lg,
## ]
## }
## }
## return(growth.matrix)
## }
## g3a_grow_apply <- function(growth.matrix, wgt.matrix, input_num,
## input_wgt) {
## na <- dim(growth.matrix)[[1]]
## growth.matrix <- growth.matrix * as.vector(input_num)
## wgt.matrix <- growth.matrix * (wgt.matrix + as.vector(input_wgt))
## growth.matrix.sum <- colSums(growth.matrix)
## return(array(c(growth.matrix.sum, colSums(wgt.matrix)/avoid_zero(growth.matrix.sum)),
## dim = c(na, 2)))
## }
## assert_msg <- function(expr, message) {
## if (isFALSE(expr)) {
## warning(message)
## return(TRUE)
## }
## return(FALSE)
## }
## ling_imm__plusdl <- 10
## ling_imm__minage <- 3L
## ling_imm__maxage <- 10L
## ling_imm__growth_w <- array(NA, dim = structure(6:5, names = c("length",
## "delta")), dimnames = list(length = c("0:10", "10:20",
## "20:30", "30:40", "40:50", "50:Inf"), delta = c("0",
## "1", "2", "3", "4")))
## ling_imm__growth_l <- array(NA, dim = structure(6:5, names = c("length",
## "delta")), dimnames = list(length = c("0:10", "10:20",
## "20:30", "30:40", "40:50", "50:Inf"), delta = c("0",
## "1", "2", "3", "4")))
## ling_imm__num <- array(0, dim = c(length = 6L, area = 2L,
## age = 8L), dimnames = list(length = c("0:10", "10:20",
## "20:30", "30:40", "40:50", "50:Inf"), area = c("area1",
## "area2"), age = c("age3", "age4", "age5", "age6", "age7",
## "age8", "age9", "age10")))
## ling_imm__wgt <- array(1, dim = c(length = 6L, area = 2L,
## age = 8L), dimnames = list(length = c("0:10", "10:20",
## "20:30", "30:40", "40:50", "50:Inf"), area = c("area1",
## "area2"), age = c("age3", "age4", "age5", "age6", "age7",
## "age8", "age9", "age10")))
## ling_imm__prevtotal <- 0
## while (TRUE) {
## {
## custom_delta_w <- (ling_imm__plusdl * 44)
## {
## comment("g3a_grow for ling_imm")
## for (age in seq(ling_imm__minage, ling_imm__maxage,
## by = 1)) {
## ling_imm__age_idx <- age - ling_imm__minage +
## 1L
## for (ling_imm__area_idx in seq_along(ling_imm__areas)) {
## area <- ling_imm__areas[[ling_imm__area_idx]]
## growth_delta_w <- (ling_imm__growth_w[] <- area *
## custom_delta_w * 20)
## custom_delta_l <- (area * 99)
## growth_delta_l <- (ling_imm__growth_l[] <- growth_bbinom(avoid_zero((age *
## custom_delta_l * 10)/ling_imm__plusdl),
## 4L, avoid_zero(param[["ling.bbin"]])))
## growthmat_w <- g3a_grow_matrix_wgt(growth_delta_w)
## growthmat_l <- g3a_grow_matrix_len(growth_delta_l)
## growthresult <- g3a_grow_apply(growthmat_l,
## growthmat_w, ling_imm__num[, ling_imm__area_idx,
## ling_imm__age_idx], ling_imm__wgt[, ling_imm__area_idx,
## ling_imm__age_idx])
## {
## if (FALSE)
## ling_imm__prevtotal <- sum(ling_imm__num[,
## ling_imm__area_idx, ling_imm__age_idx])
## comment("Update ling_imm using delta matrices")
## ling_imm__num[, ling_imm__area_idx, ling_imm__age_idx] <- growthresult[,
## (1)]
## ling_imm__wgt[, ling_imm__area_idx, ling_imm__age_idx] <- growthresult[,
## (2)]
## if (FALSE)
## assert_msg(~abs(ling_imm__prevtotal -
## sum(ling_imm__num[, ling_imm__area_idx,
## ling_imm__age_idx])) < 1e-04, "g3a_growmature: ling_imm__num totals are not the same before and after growth")
## }
## }
## }
## }
## }
## }
## }
## <bytecode: 0x55c0a87fbfd8>
## <environment: 0x55c0ad740498>Note that custom_delta_l is recalculated every loop,
since it uses area, whereas custom_delta_w is
calculated once for the step.
In the g3a_growmature function above, we see a reference
to g3_param("ling.bbin"). The g3_param()
function is pseudo-code that specifies the model should accept a
parameter at this point. In the R code, this has been converted to
param[["ling.bbin"]], so when we call our R function, we
can provide a value, e.g. fn(list(ling.bbin = 6)). We can
also use g3_param_vector() to provide a model with a vector
of values. See ?g3_param for more information on this and
other functions available to you.
When converting to TMB, there are a lot more options for using the optimisation features it offers.
Any useful model will have multiple actions, so the outputs from the
g3a_ functions need to be combined. To do this, you can
pass a list of actions to any of the g3_to_* functions, for
example:
ling_model <- g3_to_r(list(
g3a_age(ling_imm),
g3a_growmature(
ling_imm,
impl_f = g3a_grow_impl_bbinom(
delta_len_f = ~age * 10,
delta_wgt_f = ~area * 20,
beta_f = ~g3_param("ling.bbin"),
maxlengthgroupgrowth = 4)),
g3a_time(1990, 1999)))A useful technique is to break down the actions into separate lists, e.g.
ling_imm_actions <- list(
g3a_age(ling_imm),
g3a_growmature(
ling_imm,
impl_f = g3a_grow_impl_bbinom(
delta_len_f = ~age * 10,
delta_wgt_f = ~area * 20,
beta_f = ~g3_param("ling.bbin"),
maxlengthgroupgrowth = 4)))
time_actions <- list(
g3a_time(1990, 1999))
ling_model <- g3_to_r(c(ling_imm_actions, time_actions))Your actions can be combined in any order, the reasons for why are in the next section.
Gadget2 has a strict Order of Calculations, and re-ordering calculations may well have averse effects. Actions also have common steps to perform, for example overstocking for a prey has to be calculated, and calculated once, after each fleet has finished harvesting.
To manage this, formulas within an action are labelled with a string, for example:
lln <- g3_fleet('lln') %>% g3s_livesonareas(1)
action <- g3a_predate_fleet(
lln,
list(ling_imm),
suitabilities = list(
ling_imm = g3_suitability_exponentiall50(
~g3_param('ling.lln.alpha'),
~g3_param('ling.lln.l50')),
ling_mat = g3_suitability_exponentiall50(
~g3_param('ling.lln.alpha'),
~g3_param('ling.lln.l50'))),
catchability_f = g3a_predate_catchability_totalfleet(1))
names(action)## [1] "000:000:g3a_suitability_report:lln :ling_imm "
## [2] "003:g3a_predate :005:lln :ling_imm "
## [3] "003:g3a_predate :004:ling_imm "
## [4] "003:g3a_predate :000:000:ling_imm "
## [5] "003:g3a_predate :000:001:lln "
## [6] "003:g3a_predate :001:lln :ling_imm :be1c9cb11c841a862f53"
## [7] "003:g3a_predate :002:lln :ling_imm :be1c9cb11c841a862f53"
## [8] "003:g3a_predate :099:lln :ling_imm "We see each step has a colon-separated name, including the following parts:
run_at parameter of
g3a_predate_fleet(), or any other action.g3a_predate_fleet() call.When a G3 model is made, it will sort the combined list of steps by name, and remove duplicates by name. This means that:
003:001) before
overconsumption is calculated (003:004). There’s no explict
process needed to interleave multiple predation actions.