---
title: "Debugging a gadget3 model"
output:
  html_document:
    toc: true
    theme: null
vignette: >
  %\VignetteIndexEntry{Model Debugging}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r, message=FALSE, echo=FALSE}
library(gadget3)
library(magrittr)
if (nzchar(Sys.getenv('G3_TEST_TMB'))) options(gadget3.tmb.work_dir = gadget3:::vignette_base_dir('work_dir'))
```

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:

```{r, eval=FALSE}
# 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()``:

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
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()
```

## Increasing objective function verbosity

TMB objective functions have several options to increase logging output

```{r, eval=FALSE}
obj.fn <- g3_tmb_adfun(
    model_cpp,
    params.in,
    # control allows you to tweak trace level for outer optimiser, see "?stats::optim"
    control = list(trace = 2),
    # inner.control allows you to increase trace level for the random effect optimiser, see "?TMB::newton"
    inner.control = list(trace = 3, maxit = 1000, tol = 1e-7) )
obj.fn$env$tracepar <- TRUE  # Print result of every likelihood evaluation
obj.fn$env$tracemgc <- TRUE  # Print maximum gradient at every gradient evaluation
obj.fn$env$silent <- FALSE  # Show/hide tracing messages
```

## 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:

```{r, eval=FALSE}
# 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:

```{r, eval=FALSE}
obj.fn <- g3_tmb_adfun(bounded_code, params.in, inner.control = list(trace = 3, maxit = 100))
```

#### Logging messages

```
Trying early exit - PD Hess?
```

Not much improvement to the solution has been found, so TMB checks to see if the current solution has a
[positive-definite hessian](https://en.wikipedia.org/wiki/Second_partial_derivative_test#Functions_of_many_variables).
If yes, we can stop here.

#### 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:

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
obj.fn$env$ff <- edit(obj.fn$env$ff)
```

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

```{r, eval=FALSE}
        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:

```{r, eval=FALSE}
obj.fn$fn()
do.call(TMB:::newton, newt.args)
```

...or modify the newton function:

```{r, eval=FALSE}
newt <- TMB:::newton
newt <- edit(newt)
do.call(newt, newt.args)
```