Structure of a gadget3 model

R formula, or the tilde operator

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: 0x55c15c6765c8>

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 8

We 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 5

Note that:

R hasn’t at any point worked out cows + pigs, the code was stored for later use.
The environment (i.e. all variables defined at that point) is “remembered”

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:

g3_to_r(list(f))

## function (param) 
## {
##     cur_time <- -1L
##     cows <- 16
##     pigs <- 32
##     while (TRUE) {
##         cur_time <- cur_time + 1L
##         cows + pigs
##     }
## }
## <bytecode: 0x55c15ca8f5c0>
## <environment: 0x55c15bc9b5c0>

…or a TMB template:

g3_to_tmb(list(f))

## #include <TMB.hpp>
## 
## 
## 
## template<class Type>
## Type objective_function<Type>::operator() () {
##     DATA_SCALAR(reporting_enabled); DATA_UPDATE(reporting_enabled);
##     int cur_time = -1;
##     Type cows = (double)(16);
##     Type pigs = (double)(32);
## 
##     while (true) {
##         cur_time += 1;
##         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.

Actions

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)

## $`001`
## ~{
##     debug_label("g3a_time: Start of time period")
##     if (cur_time == 0 && assert_msg(g3_param("retro_years", value = 0, 
##         optimise = FALSE) >= 0, "retro_years must be >= 0")) 
##         return(NaN)
##     if (cur_time == 0 && assert_msg(project_years >= 0, "project_years must be >= 0")) 
##         return(NaN)
##     if (cur_time > total_steps) 
##         return(nll)
##     cur_year <- start_year + (cur_time%/%step_count)
##     cur_year_projection <- cur_year > end_year - g3_param("retro_years", 
##         value = 0, optimise = FALSE)
##     cur_step <- (cur_time%%step_count) + 1L
##     cur_step_final <- cur_step == step_count
## }
## <environment: 0x55c15c001dc8>

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.

g3_to_r(g3a_time(1990, 1999))

## function (param) 
## {
##     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"]]
##     step_lengths <- 12L
##     end_year <- 1999L
##     retro_years <- param[["retro_years"]]
##     start_year <- 1990L
##     total_steps <- length(step_lengths) * (end_year - retro_years - 
##         start_year + project_years) + length(step_lengths) - 
##         1L
##     nll <- 0
##     cur_year <- 0L
##     step_count <- length(step_lengths)
##     cur_year_projection <- FALSE
##     cur_step <- 0L
##     cur_step_final <- FALSE
##     while (TRUE) {
##         cur_time <- cur_time + 1L
##         {
##             comment("g3a_time: Start of time period")
##             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)
##             if (cur_time > total_steps) 
##                 return(nll)
##             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_final <- cur_step == step_count
##         }
##     }
## }
## <bytecode: 0x55c15d286bb8>
## <environment: 0x55c15d8ad528>

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:

g3_to_r(g3a_time(1990, ~start_year + 4 ))

## function (param) 
## {
##     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"]]
##     step_lengths <- 12L
##     start_year <- 1990L
##     end_year <- start_year + 4
##     retro_years <- param[["retro_years"]]
##     total_steps <- length(step_lengths) * (end_year - retro_years - 
##         start_year + project_years) + length(step_lengths) - 
##         1L
##     nll <- 0
##     cur_year <- 0L
##     step_count <- length(step_lengths)
##     cur_year_projection <- FALSE
##     cur_step <- 0L
##     cur_step_final <- FALSE
##     while (TRUE) {
##         cur_time <- cur_time + 1L
##         {
##             comment("g3a_time: Start of time period")
##             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)
##             if (cur_time > total_steps) 
##                 return(nll)
##             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_final <- cur_step == step_count
##         }
##     }
## }
## <bytecode: 0x55c15ae5b9d0>
## <environment: 0x55c15bdd7c68>

Stocks

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:

The dimensions one would use if making an array to store data on that stock. For example if we want an array for number of individuals, what lengthgroups do we use? How many ages do we store? How many (and which) areas do we have?
What code do we need to iterate over the stock? Say I want to add 1 individual to each lengthgroup, how do I loop over all the other dimensions?
If looping over one stock, how do I find corresponding entries in another stock? For example, my fleet is only interested in prey that is in the same area.

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     NA

To 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    NA

When 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) 
## {
##     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)
##     }
##     avoid_zero_vec <- function(a) {
##         (pmax(a * 1000, 0) + log1p(exp(pmin(a * 1000, 0) - pmax(a * 
##             1000, 0))))/1000
##     }
##     avoid_zero <- function(a) {
##         (pmax(a * 1000, 0) + log1p(exp(pmin(a * 1000, 0) - pmax(a * 
##             1000, 0))))/1000
##     }
##     g3a_grow_apply <- function(delta_l, delta_w, input_num, input_wgt) {
##         na <- dim(delta_l)[[1]]
##         n <- dim(delta_l)[[2]] - 1
##         avoid_zero_vec <- function(a) {
##             (pmax(a * 1000, 0) + log1p(exp(pmin(a * 1000, 0) - 
##                 pmax(a * 1000, 0))))/1000
##         }
##         growth.matrix <- array(0, c(na, na))
##         wgt.matrix <- array(0, c(na, na))
##         for (lg in 1:na) {
##             if (lg == na) {
##                 growth.matrix[na, na] <- sum(delta_l[lg, ])
##                 wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na - 
##                   lg + 1)]
##             }
##             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)])
##                 wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na - 
##                   lg + 1)]
##             }
##             else {
##                 growth.matrix[lg, lg:(n + lg)] <- delta_l[lg, 
##                   ]
##                 wgt.matrix[lg, lg:(n + lg)] <- delta_w[lg, ]
##             }
##         }
##         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_vec(growth.matrix.sum)), 
##             dim = c(na, 2)))
##     }
##     assert_msg <- function(expr, message) {
##         if (isFALSE(expr)) {
##             warning(message)
##             return(TRUE)
##         }
##         return(FALSE)
##     }
##     cur_time <- -1L
##     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__prevtotal <- 0
##     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")))
##     while (TRUE) {
##         cur_time <- cur_time + 1L
##         {
##             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]]
##                   {
##                     comment("Calculate length/weight delta matrices for current lengthgroups")
##                     ling_imm__growth_l[] <- growth_bbinom(avoid_zero_vec((age * 
##                       10)/ling_imm__plusdl), 4L, avoid_zero(param[["ling.bbin"]]))
##                     ling_imm__growth_w[] <- area * 20
##                     if (FALSE) 
##                       ling_imm__prevtotal <- sum(ling_imm__num[, 
##                         ling_imm__area_idx, ling_imm__age_idx])
##                     comment("Update ling_imm using delta matrices")
##                     {
##                       growthresult <- g3a_grow_apply(ling_imm__growth_l, 
##                         ling_imm__growth_w, ling_imm__num[, ling_imm__area_idx, 
##                           ling_imm__age_idx], ling_imm__wgt[, 
##                           ling_imm__area_idx, ling_imm__age_idx])
##                       {
##                         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: 0x55c159e24788>
## <environment: 0x55c15bee0d10>

…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) 
## {
##     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)
##     }
##     avoid_zero_vec <- function(a) {
##         (pmax(a * 1000, 0) + log1p(exp(pmin(a * 1000, 0) - pmax(a * 
##             1000, 0))))/1000
##     }
##     avoid_zero <- function(a) {
##         (pmax(a * 1000, 0) + log1p(exp(pmin(a * 1000, 0) - pmax(a * 
##             1000, 0))))/1000
##     }
##     g3a_grow_apply <- function(delta_l, delta_w, input_num, input_wgt) {
##         na <- dim(delta_l)[[1]]
##         n <- dim(delta_l)[[2]] - 1
##         avoid_zero_vec <- function(a) {
##             (pmax(a * 1000, 0) + log1p(exp(pmin(a * 1000, 0) - 
##                 pmax(a * 1000, 0))))/1000
##         }
##         growth.matrix <- array(0, c(na, na))
##         wgt.matrix <- array(0, c(na, na))
##         for (lg in 1:na) {
##             if (lg == na) {
##                 growth.matrix[na, na] <- sum(delta_l[lg, ])
##                 wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na - 
##                   lg + 1)]
##             }
##             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)])
##                 wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na - 
##                   lg + 1)]
##             }
##             else {
##                 growth.matrix[lg, lg:(n + lg)] <- delta_l[lg, 
##                   ]
##                 wgt.matrix[lg, lg:(n + lg)] <- delta_w[lg, ]
##             }
##         }
##         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_vec(growth.matrix.sum)), 
##             dim = c(na, 2)))
##     }
##     assert_msg <- function(expr, message) {
##         if (isFALSE(expr)) {
##             warning(message)
##             return(TRUE)
##         }
##         return(FALSE)
##     }
##     cur_time <- -1L
##     ling_imm__plusdl <- 10
##     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__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__prevtotal <- 0
##     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")))
##     while (TRUE) {
##         cur_time <- cur_time + 1L
##         {
##             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]]
##                     custom_delta_l <- (area * 99)
##                     {
##                       comment("Calculate length/weight delta matrices for current lengthgroups")
##                       ling_imm__growth_l[] <- growth_bbinom(avoid_zero_vec((age * 
##                         custom_delta_l * 10)/ling_imm__plusdl), 
##                         4L, avoid_zero(param[["ling.bbin"]]))
##                       ling_imm__growth_w[] <- area * custom_delta_w * 
##                         20
##                       if (FALSE) 
##                         ling_imm__prevtotal <- sum(ling_imm__num[, 
##                           ling_imm__area_idx, ling_imm__age_idx])
##                       comment("Update ling_imm using delta matrices")
##                       {
##                         growthresult <- g3a_grow_apply(ling_imm__growth_l, 
##                           ling_imm__growth_w, ling_imm__num[, 
##                             ling_imm__area_idx, ling_imm__age_idx], 
##                           ling_imm__wgt[, ling_imm__area_idx, 
##                             ling_imm__age_idx])
##                         {
##                           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: 0x55c15b44cfd8>
## <environment: 0x55c15bf38230>

Note that custom_delta_l is recalculated every loop, since it uses area, whereas custom_delta_w is calculated once for the step.

Model parameterization

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.

Combining actions

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.

Ordering of actions

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] "003:000:lln                 "                                          
## [2] "003:000:ling_imm            "                                          
## [3] "003:004:ling_imm            "                                          
## [4] "003:002:lln                 :ling_imm            :3e19f1f5997cd960b5a2"
## [5] "003:001:lln                 :ling_imm            :3e19f1f5997cd960b5a2"
## [6] "003:005:lln                 :ling_imm            "                     
## [7] "003:005:lln                 "                                          
## [8] "003:003:lln                 :ling_imm            "

We see each step has a colon-separated name, including the following parts:

A number corresponding to the order of calculations. If you look at the gadget2 user-guide, consumption is step 3. This can be controlled using the run_at parameter of g3a_predate_fleet(), or any other action.
Another number specifying the order within an action this step should be performed.
The prey and where relevant the fleet name
A hash, that corresponds to the parameters of the 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:

Predation will happen at the correct part of the model cycle.
Predation steps themselves will be ordered. So all defined predation steps will collect their catch (003:001) before overconsumption is calculated (003:004). There’s no explict process needed to interleave multiple predation actions.
The overconsumption co-efficient for ling_imm will only be calculated once, regardless of the number of predation actions, since we will remove the duplicate versions of this step.
Conversely the hash ensures that there are never duplicate versions of the main predation step, so could be repated. This is more relevant for renewal, where there can be multiple renewal actions working on the same stock.