The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Structure of a gadget3 model

The following describes the structure of a gadget3 model, from the bottom up.

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:

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:

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:

When a G3 model is made, it will sort the combined list of steps by name, and remove duplicates by name. This means that:

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.