There are also special g3 functions that can be used in formula that affect the resulting code, rather than just being called when run.
See ?g3_param for more information.
All G3 models have the G3 global environment included, which contains
the definition of generally useful functions. For example
avoid_zero(), used to avoid div/0 conditions by adding a
small amount of error close to zero.
See R/aab_env.R for more information, and other existing
helpers.
R functions can be included in a formula as anything else, provided you use this through R. For example:
model_fn <- g3_to_r(list(g3a_time(1990, 1990), g3_formula(
nll <- nll + fn(9),
fn = function (x) x * 10 )))
model_fn()Your function is now part of the model environment
environment(model_fn)$fnWe do not yet translate functions from TMB to R,however
g3_native allows you to Some things aren’t easy to do with
code translation. g3_native allows you to define a function
with separate R and C++ definitions, for example:
normalize_vec <- g3_native(r = function (a) {
a / sum(a)
}, cpp = '[](vector<Type> a) -> vector<Type> {
return a / a.sum();
}')
model_fn <- g3_to_r(list(g3a_time(1990, 1990), g3_formula(
nll <- nll + normalize_vec(g3_param("a")),
normalize_vec = normalize_vec )))
model_fn(c(list(a = 10:20), attr(model_fn, "parameter_template")))
model_code <- g3_to_tmb(list(g3a_time(1990, 1990), g3_formula(
nll <- nll + normalize_vec(g3_param("a")),
normalize_vec = normalize_vec )))
model_codeMost actions interact with stocks, and fill out abstract formulae
with the stocks provided to the function. To do this you need to do a
series of substitutions, which are handled by g3_step().
This takes a formula, looks for stock_* named functions and
mangles the formula as appropriate. For example, a snippet from
action_mature.R.
out <- new.env(parent = emptyenv())
out[[step_id(run_at, 1, stock)]] <- g3_step(f_substitute(~{
debug_label("g3a_mature for ", stock)
# Matured stock will weigh the same
stock_with(stock, stock_with(matured, matured__wgt <- stock__wgt))
stock_iterate(stock, stock_intersect(matured, if (run_f) {
debug_label("Move matured ", stock, " into temporary storage")
stock_ss(matured__num) <- stock_ss(stock__num) * maturity_f
stock_ss(stock__num) <- stock_ss(stock__num) - stock_ss(matured__num)
}))
}, list(run_f = run_f, maturity_f = maturity_f)))Assume that stock has name “ling_imm” and
matured has name “ling_imm_maturing”.
The first line uses debug_label() to produce a
debug_label() function call,
debug_label("g3a_mature for ling_imm"). This will be used
as a code comment and a label for this block if producing diagrams.
Next stock_with() is used to to transform
matured__wgt <- stock__wgt to use the proper stock
names. We don’t care about dimensions since we’re copying over all the
data.
Finally, we use a combination of stock_iterate() and
stock_intersect(). stock_iterate() will create
a loop that loops over all of the stock’s dimensions, and
stock_ss() will subset stock__num, prividing
1-dimension lengthgroup vector. stock_intersect().
These iterators will then be available to the maturity_f
that the users provide, as demonstrated in the [Stocks] section.
For more information on the stock_* functions, see
?stock_ss.
Obviously R and C++’s type systems are quite different, and gadget3 attempts to bridge the gap.
In R, there is no difference between a scalar and a 1-element vector.
If you expect to treat the variable as a vector or array, then state
this explicitly with array, even if the result may be a
1-element vector. This means that methods that work on TMB array or
vector classes will be available.
One needs to be more careful with the type of constants in C++ than
R. In particular, x / 2 means integer division. As a
result, G3 will cast any numeric constant as a double outside of certain
situations, e.g. array indices, which will be integer values. However,
if you do want an integer it’s best to express this explictly,
i.e. 3L. Note that there is no difference in R code between
3 and 3.0.
R forumlas all you to include extra definitions when defining a formula. This allows you to break up a definition into more readable chunks. For example:
ling_imm <- g3_stock('ling_imm', seq(0, 50, 10)) %>%
g3s_age(3, 10)
nmort <- function() {
E <- ~stock__minlen * age
F <- ~stock__minage # TODO: Does this work now?
~E * F * 4
}
g3_to_r(list(g3a_naturalmortality(ling_imm, nmort())))Note that:
E refers to age,
gadget3 has automatically inserted it into the loop.F however can be defined outside the loop, so is.If you need to have something available to ther steps, it can be
defined using g3_global_formula and providing an
init_val:
ling_imm <- g3_stock('ling_imm', seq(0, 50, 10)) %>%
g3s_age(3, 10)
nmort <- function() {
# Define a counter
E <- g3_global_formula(
~E + 1, init_val = 0L)
# We can just give init_val, to define something global to the model
F <- g3_global_formula(
init_val = 99L)
~E * F * 4
}
g3_to_r(list(g3a_naturalmortality(ling_imm, nmort())))As well as making values available to other steps,
g3_global_formula() can also be used to ensure that the
value ends up in the model report, which will happen automatically for
any non-constant global in the model.
Additional steps can be attached to a formula, allowing setup for a variable at the beginning of a model, or implicit likelihood components. For example:
st_imm <- g3_stock(c("st", "imm"), 1:10)
g3_to_r(list(
g3a_naturalmortality(
st_imm,
g3_formula(
parrot**2,
parrot = 0,
"-01:ut:parrot" = g3_formula({
parrot <- runif(1)
}))),
NULL ))## 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)
## }
## parrot <- 0
## st_imm__num <- array(0, dim = c(length = 10L), dimnames = list(length = c("1:2",
## "2:3", "3:4", "4:5", "5:6", "6:7", "7:8", "8:9", "9:10",
## "10:Inf")))
## while (TRUE) {
## {
## parrot <- runif(1)
## }
## {
## comment("Natural mortality for st_imm")
## st_imm__num[] <- st_imm__num[] * (parrot^2)
## }
## }
## }
## <bytecode: 0x562b4ccb6300>
## <environment: 0x562b4bc1a210>The attached step gets inserted at the beginning of the model, and a
new random number is chosen for parrot at each model
timestep.