Skip to content

Code snippets

Jump to bottom
Mollie Brooks edited this page Oct 3, 2021 · 39 revisions

Place your useful code snippets here for the benefit of others.

  • Your C++ tricks
  • Techniques that need highlighting even if they are in the docs already
  • Snippets may later migrate into documentation

Ensuring function remains positive (stolen from DF and ADMB manual with shame!!! - and modified such that branching is being taped, see also Things you should NOT do in TMB) Note, thas this function only works for models without random effects.

template<class Type>
Type posfun(Type x, Type eps, Type &pen){
  pen += CppAD::CondExpLt(x, eps, Type(0.01) * pow(x-eps,2), Type(0));
  return CppAD::CondExpGe(x, eps, x, eps/(Type(2)-x/eps));
}

Element-wise product of matrices. Remember that m1*m1 is a linear algebra matrix multiplication. To do element-wise add ".array()" to it.

matrix<double> m1(2,2);
m1.array()*m1.array();

Saving and loading objects created by MakeADFun

After saving an object obj created by MakeADFun the external pointers are lost, i.e. all functions in he list cannot be directly used. Running retape re-creates all external pointers.

## Creating an object, saving it and deleting it.
obj <- MakeADFun(data, parameters)
save(obj, "temp.RData")
rm(obj)

## After loading and object, `retape` re-creates the external pointers
load("temp.RData")
## obj$fn(), obj$report() fail to run because of NULL pointers
obj$retape()
obj$fn() ## Works!
obj$report() ## Also works.

Automatically modifying parameter bounds when using map

During model development it is convenient to turn parameters on/off with map. Here is how to automatically modify upper and lower bounds for nlminb() or optim():

# Set bounds for all parameter 
L = list(a=3,b=-1,d=0.001,sigma=.0001)
U = list(a=5,b=2, d=10, sigma=10)

# List parameters that should be fixed in nlminb() or optim()
map = list(a=factor(NA))
#map=list()  # All parameters active

# Remove inactive parameters from bounds
member <- function(x,y) !is.na(match(x,y))
L = unlist(L[!member(names(L),names(map))])
U = unlist(U[!member(names(U),names(map))])

obj <- MakeADFun(data=data,parameters = parameters,map=map)
opt <- nlminb(obj$par,obj$fn,obj$gr,lower=L,upper=U)

# Or more simply:
# Set bounds for all parameter
L = c(a=3,b=-1,d=0.001,sigma=.0001)
U = c(a=5,b=2, d=10, sigma=10)

# List parameters that should be fixed in nlminb() or optim()
map = list(a=factor(NA), d=factor(NA))
#map=list()  # All parameters active

# Remove inactive parameters from bounds
(L <- L[-match(names(map), names(L))])
(U <- U[-match(names(map), names(U))])

obj <- MakeADFun(data=data,parameters = parameters,map=map)
opt <- nlminb(obj$par,obj$fn,obj$gr,lower=L,upper=U)

Convert estimates and SE from tabular to list format

pl <- model$env$parList()
jointrep <- sdreport(model, getJointPrecision=TRUE)
allsd <- sqrt(diag(solve(jointrep$jointPrecision)))
plsd <- model$env$parList(par=allsd)

This becomes useful when we work with larger models (this snippet is shamelessly snatched from Anders Nielsen)

The above failed however, the following gave correlation matrix: cov2cor(jointrep$cov.fixed)

For windows machines, debugging has previously caused the R terminal to crash. This behavior is demonstrated in the file

https://github.com/James-Thorson/TMB_experiments/blob/master/windows_debugger/problem.R

However, the gdbsource function can be used to identify the line number of the CPP file, by loading the TMBdebug package:

devtools::install_github("kaskr/TMB_contrib_R/TMBdebug")
library( TMBdebug )

Visualize sparse Hessian

Visualize the sparseness of the Hessian (for a random-effects model):

model$env$spHess(random=TRUE)

Evaluate model (mceval)

Evaluate reported model predictions using any parameter values, e.g., current biomass from MCMC draws:

model$report(mcmc.out[1,])

See what DLLs are loaded

This is presumably helpful for the forgetful among us...

TMB:::getUserDLL()

RStudio: jump directly to line with the first compilation error

  • Setup: see "RStudio integration" below.
  • It is recommended to have run precompile() to speed up compilation.

In RStudio:

  • Load TMB: library(TMB) in the R console
  • Open your .cpp file
  • Mark the "Source on Save" field, and save the file. This will compile and jump to the first error.

RStudio integration

From TMB version 1.7.12 one can integrate TMB with Rstudio by TMB:::setupRStudio(). The integration is triggered per session by library(TMB). Wait 30-60 sec for features to load in the session.

What works:

  • No side effects if not running RStudio or not running TMB.
  • Static code analyzer gives C++ auto-completion (start typing then hit tab key) and on-the-fly error checking.
  • Snippet to quickly setup an initial model template (by typing 'tmb' in an empty R script).
  • Compile error messages are displayed nicely in separate window (however note that getwd() must be the containing folder of the cpp file).
  • 'sourceCpp' button compiles cpp file and displays errors.
  • Works in conjunction with precompile().
  • Also works in conjunction with add-on package 'TMBdebug'.
  • Tested in Linux and Windows with RStudio version 1.0.35; Mac tested with RStudio 1.1.442.

Next snippet...

Clone this wiki locally