Writing G3 Actions

Global & native functions

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 logspace_add_vec, used in many actions to avoid div/0, has the following definition:

# vector<Type> form of logspace_add
g3_env$logspace_add_vec <- g3_native(r = function(a,b) {
    # https://github.com/kaskr/adcomp/issues/7#issuecomment-642559660
    pmax(a, b) + log1p(exp(pmin(a,b) - pmax(a, b)))
}, cpp = '[](vector<Type> a, Type b) -> vector<Type> {
    vector<Type> res(a.size());
    for(int i = 0; i < a.size(); i++) {
        res[i] = logspace_add(a[i], b);
    }
    return res;
}')

Stock steps

Most 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.

Writing R code destined for C++

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.

Sub-formulas and g3_global_formula

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:

Because the definition of E refers to age, gadget3 has automatically inserted it into the loop.
F however can be defined outside the loop, so is.
They are all defined in the step, not necessarily for the whole model function. In TMB this is enforced with scoping.

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.

g3_native

There are limits to what the Gadget3 transpilation can manage, particularly once you need to operate on matrices. To get around this, you can define functions using g3_native, which annotates functions with a C++ definition. For example, it is used to define logspace_add_vec:

# vector<Type> form of logspace_add
g3_env$logspace_add_vec <- g3_native(r = function(a,b) {
    pmax(a, b) + log1p(exp(pmin(a,b) - pmax(a, b)))
}, cpp = '[](vector<Type> a, Type b) -> vector<Type> {
    vector<Type> res(a.size());

    for(int i = 0; i < a.size(); i++) {
        res[i] = logspace_add(a[i], b);
    }

    // NB: This is a string within R code, so we have to escape slashes.
    //     For example, "\\n"

    return res;
}')