run_tmb: Gadget3 actions into TMB code

In gadget3: Globally-Applicable Area Disaggregated General Ecosystem Toolbox V3

Gadget3 actions into TMB code

Description

Turn g3 actions into CPP code that can be compiled using TMB

Usage

g3_to_tmb(actions, trace = FALSE, strict = FALSE)

g3_tmb_adfun(cpp_code, parameters = attr(cpp_code, "parameter_template"),
    compile_flags =
        if (.Platform$OS.type == "windows") c("-O1", "-march=native")
        else c("-O3", "-flto", "-march=native"),
    work_dir = tempdir(),
    output_script = FALSE, ...)

g3_tmb_par(parameters, include_random = TRUE)

g3_tmb_lower(parameters)

g3_tmb_upper(parameters)

g3_tmb_parscale(parameters)

g3_tmb_relist(parameters, par)

Arguments

actions	A list of actions (i.e. list of formula objects), as produced by g3a_* functions.
trace	If TRUE, turn all comments into print statements.
strict	If TRUE, enable extra sanity checking in actions. Any invalid conditions (e.g. more/less fish after growth) will result in a warning.
cpp_code	cpp_code as produced by g3_to_tmb.
parameters	Parameter table as produced by attr(g3_to_tmb(...), 'parameter_template'), modified to provide initial conditions, etc.
par	Parameter vector, as produced by one of nlminb(...)$par obj.fun$env$last.par g3_tmb_par() The first will not include random parameters by default, the others will.
include_random	Should random parameters assumed to be part of par? Should be TRUE if using obj.fun$fn, obj.fun$report directly, e.g. obj.fun$fn(g3_tmb_par(param_tbl)). In other cases, FALSE.
compile_flags	List of extra flags to compile with, use e.g. "-g" to enable debugging output.
work_dir	Directory to write and compile .cpp files in. Defaults to R's current temporary directory
output_script	If TRUE, create a temporary R script that runs MakeADFun, and return the location. This can then be directly used with gdbsource or callr::rscript.
...	Any other options handed directly to MakeADFun

Details

g3_tmb_adfun

g3_tmb_adfun will do both the compile and MakeADFun steps of making a model. If the code is identical to an already-loaded model then it won't be recompiled, so repeated calls to g3_tmb_adfun to change parameters are fast.

If MakeADFun is crashing your R session, then you can use output_script to run in a separate R session. Use this with gdbsource to debug your model.

Value

g3_to_tmb

A string of C++ code that can be used as an input to g3_tmb_adfun, with the following attributes:

actions

The original actions list given to the function

model_data

An environment containing data attached to the model

parameter_template

A data.frame to be filled in and used as parameters in the other g3_tmb_* functions

Use e.g. attr(cpp_code, 'parameter_template') to retrieve them.

g3_tmb_adfun

An ADFun as produced by TMB's MakeADFun, or location of temporary script if output_script is TRUE

g3_tmb_par

Values extracted from parameters table converted into a vector of values for obj$fn(par) or nlminb

g3_tmb_lower

Lower bounds extracted from parameters table converted into a vector of values for nlminb. Random parameters are always excluded

g3_tmb_upper

Lower bounds extracted from parameters table converted into a vector of values for nlminb. Random parameters are always excluded

g3_tmb_parscale

Parscale extracted from parameters table, converted into a vector of values for nlminb. Random parameters are always excluded

g3_tmb_relist

The parameters table value column, but with optimised values replaced with contents of par vector. i.e. the inverse operation to g3_tmb_par. par can either include or discount random variables.

Examples

ling_imm <- g3_stock(c(species = 'ling', 'imm'), seq(20, 156, 4)) %>% g3s_age(3, 10)

initialconditions_action <- g3a_initialconditions_normalparam(
    ling_imm,
    factor_f = g3a_renewal_initabund(by_stock_f = 'species'),
    by_stock = 'species',
    by_age = TRUE)

abundance_action <- g3l_abundancedistribution(
    'abundance',
    data.frame(year = 2000:2004, number = 100),
    stocks = list(ling_imm),
    function_f = g3l_distribution_sumofsquares())

# Timekeeping action
time_action <- g3a_time(
    start_year = 2000,
    end_year = 2004,
    c(3, 3, 3, 3))

# Generate a model from the above 2 actions
# NB: Obviously in reality we'd need more actions
cpp <- g3_to_tmb(list(initialconditions_action, abundance_action, time_action))

if (interactive()) {
  # Edit the resulting code
  cpp <- edit(cpp)
}

# Set initial conditions for parameters
tmb_param <- attr(cpp, 'parameter_template')
tmb_param$value$project_years <- 0
tmb_param$value$ling.init.F <- 0.4
tmb_param$value$ling.Linf <- 160
tmb_param$value$ling.K <- 90
tmb_param$value$ling.t0 <- 0
tmb_param[grepl('^ling.init.sd.', rownames(tmb_param)), 'value'] <- 50.527220
tmb_param[grepl('^ling_imm.init.\\d+', rownames(tmb_param)), 'value'] <- 1
tmb_param$value$ling_imm.init.scalar <- 200
tmb_param$value$ling_imm.walpha <- 2.27567436711055e-06
tmb_param$value$ling_imm.wbeta <- 3.20200445996187
tmb_param[grepl('\\.M$', rownames(tmb_param)), 'value'] <- 0.15

# We can set lower/upper bounds for multiple properties at once with grepl()
tmb_param[grepl('.', rownames(tmb_param)), 'lower'] <- -1000
tmb_param[grepl('.', rownames(tmb_param)), 'upper'] <- 1000

# parscale gives optim() a relative scale of parameters
tmb_param['parscale'] <- 1




if (!( nzchar(Sys.getenv('GITHUB_CI')) && .Platform$OS.type == "windows" )) {
  # Compile to a TMB ADFun
  tmb <- g3_tmb_adfun(cpp, tmb_param)
}

# NB: TMB::gdbsource() requires both "R" and "gdb" to be available
# NB: gdbsource hangs on windows - https://github.com/kaskr/adcomp/issues/385
if (all(nzchar(Sys.which(c('gdb', 'R')))) && .Platform$OS.type !="windows") {

  cpp_broken <- g3_to_tmb(list(
    initialconditions_action,
    abundance_action,
    g3_formula(quote( stop("This model is broken") )),
    time_action))

  # Build the model in an isolated R session w/debugger
  writeLines(TMB::gdbsource(g3_tmb_adfun(
      cpp_broken,
      compile_flags = "-g",
      output_script = TRUE)))

}

if (!( nzchar(Sys.getenv('GITHUB_CI')) && .Platform$OS.type == "windows" )) {
  # Perform a single run, using values in table
  result <- tmb$fn(g3_tmb_par(tmb_param))
}

if (!( nzchar(Sys.getenv('GITHUB_CI')) && .Platform$OS.type == "windows" )) {
  # perform optimisation using upper/lower/parscale from table
  fit <- optim(tmb$par, tmb$fn, tmb$gr,
      method = "L-BFGS-B",
      upper = g3_tmb_upper(tmb_param),
      lower = g3_tmb_lower(tmb_param),
      control = list(maxit=10, parscale=g3_tmb_parscale(tmb_param)))
}

if (!( nzchar(Sys.getenv('GITHUB_CI')) && .Platform$OS.type == "windows" )) {
  # perform optimisation without bounds
  fit <- optim(tmb$par, tmb$fn, tmb$gr)
}

if (!( nzchar(Sys.getenv('GITHUB_CI')) && .Platform$OS.type == "windows" )) {
  # Go back to a list of parameters, suitable for the R version
  # NB: This will not set the values for random parameters
  param_list <- g3_tmb_relist(tmb_param, fit$par)
}

if (!( nzchar(Sys.getenv('GITHUB_CI')) && .Platform$OS.type == "windows" )) {
  # Update parameters with values from last run, *including* random parameters.
  param_list <- g3_tmb_relist(tmb_param, tmb$env$last.par)
}

if (!( nzchar(Sys.getenv('GITHUB_CI')) && .Platform$OS.type == "windows" )) {
  # Rebuild, only including "Fun" (i.e. without auto-differentiation)
  # Result will only work for tmb$report
  tmb <- g3_tmb_adfun(cpp, tmb_param, type = "Fun")
  result <- tmb$report(g3_tmb_par(tmb_param))
}

gadget3 documentation built on July 3, 2024, 9:07 a.m.