Debugging a gadget3 model

Source: vignettes/model-debugging.Rmd

There are several techniques available to help debug problems in gadget3 models. The following lists common errors and ways you can investigate further.

Viewing & editing model code

gadget3 models are converted to R or TMB source code, which is often useful to study to see what the model will do at any point. It is also possible to edit the source code and e.g. add debugging breaks with recover().

An R model function can be viewed / edited directly using the R edit() function:

# Create model based on your actions
model_fn <- g3_to_r(actions)

# Edit source code & re-run model
model_fn <- edit(model_fn)
model_fn(params.in)

The TMB source code could also be inspected using edit():

tmb_model <- g3_to_tmb(actions)

# Edit source code & re-run model
tmb_model <- edit(tmb_model)
obj.fn <- g3_tmb_adfun(tmb_model, params.in, Type = "Fun")
r <- obj.fn$report()

As the R & TMB model will be equivalent, you can view the source of whichever you are more comfortable with.

NaN in likelihood / reports

A non-finite value, once generated, will spread to all related parts of the model, so trying to debug by looking at the final result is often useless.

To pinpoint where the problem is, you can add g3a_trace_var() to find where the problem originates. See the manpage for how to use it.

“Error in optim(…): initial value in ‘vmmin’ is not finite”

The optimiser is telling you that the initial values you provided do not produce a finite likelihood. Add g3a_trace_var() and run your model outside the optimiser, with e.g:

tmb_model <- g3_to_tmb(c(actions, list(
  g3a_report_detail(actions),
  g3a_trace_var(actions),
  gadget3::g3l_bounds_penalty(actions) )))
obj.fn <- g3_tmb_adfun(tmb_model, params.in, Type = "Fun")
r <- obj.fn$report()

TMB model crashes your R session

If the model crashes whilst forming the TMB ADFun object, then it takes your R session with it. If this is happening, wrap g3_tmb_adfun() with TMB::gdbsource() as follows:

# Model setup will look something like this
tmb_ling <- g3_to_tmb(...)
tmb_param <- attr(tmb_ling, 'parameter_template')

writeLines(TMB::gdbsource(g3_tmb_adfun(
    tmb_ling,
    tmb_param,
    compile_flags = "-g",
    output_script = TRUE)))

output_script = TRUE tells g3_tmb_adfun() to, after compilation, write a temporary R script that will build the TMB ADFun object (and presumably crash in the process). TMB::gdbsource() in turn runs a provided R script in a fresh R session wrapped in gdb. By default it will print a stacktrace and quit.

In the stacktrace you should see a line number for where the crash occured. You can then look at the TMB model source code, find the line number and thus the action that is currently failing.

Interactive debugging of TMB models

In theory you can use interactive = TRUE with TMB::gdbsource(), however as this eats error messages it’s better to do this by hand:

> g3_tmb_adfun(tmb_ling, tmb_param, compile_flags = "-g", output_script = TRUE)
[1] "/tmp/RtmpysTVvW/file3da4a6f13a80c.R"

R -d gdb
(gdb) run --vanilla < /tmp/RtmpysTVvW/file3da4a6f13a80c.R
 . . . Compilation, crash at some point . . .
(gdb) up
(gdb) up
(gdb) up
(gdb) call ling_imm__consratio.print()
Array dim: 35 1 8
Array val: -nan -nan -nan -nan
(gdb) call ling_imm__num.print()
Array dim: 35 1 8
Array val: -nan -nan -nan -nan -nan

(gdb) print cur_time
$5 = 0
  • Print the contents of single values with std::cout << ling__Linf << std::endl;
  • Use ().cols() and ().rows() to get the size of an array expression
  • Print the contents of array objects with ling_imm__consratio.print()

Note that for the .print() method to be available for arrays, it has to be referenced at least once in the model source, otherwise it won’t be compiled in. Use edit(tmb_ling) to add it somewhere first.

Random effects

Some notes on debugging errors with random effects models.

Tracing inner model

Control arguments for the inner TMB:::newton() model can be provided to g3_tmb_adfun(), e.g. to add tracing:

obj.fn <- g3_tmb_adfun(bounded_code, params.in, inner.control = list(trace = 3, maxit = 100))

Missing value for m

The error:

Error in if (m < 0) { : missing value where TRUE/FALSE needed

…comes from TMB’s newton optimiser, and essentially says there is NaN in the hessian matrix. m in this case is equivalent to:

local({ min(diag(spHess(random = TRUE, set_tail = random[1]))) }, envir = obj.fn$env)

Missing value par - parold 

Error in if (norm(par - parold) < step.tol) { :
  missing value where TRUE/FALSE needed

The par list of parameters to optimise became NaN, for various reasons, but likely solveCholesky() failed.

Separate netwton optimisation

Generally, TMB:::newton() is called from obj.fn$env$ff. You can extract the arguments provided by editing this function:

obj.fn$env$ff <- edit(obj.fn$env$ff)

Around line 16, add code to extract the arguments provided, e.g:

        assign("newt.args", c(list(par = eval(random.start),
            fn = f0, gr = function(x) f0(x, order = 1), he = H0,
            env = env), inner.control), env = globalenv())

Then you can perform a single run to extract arguments:

obj.fn$fn()
do.call(TMB:::newton, newt.args)

…or modify the newton function:

newt <- TMB:::newton
newt <- edit(newt)
do.call(newt, newt.args)