Skip to content

Code snippets

Jump to bottom
Hans J. Skaug edited this page May 18, 2016 · 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)

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)

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 identify the line number of the CPP file, as demonstrated in the file:

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

In essence, the R debugger requires creating an R script, with the same name as the CPP file, where running the R script in an empty R terminal results in the TMB crash. By running such an R script in the "gdbsource" function, you can determine the line number of the CPP file that is responsible for the crash.

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

The following experimental code integrates TMB with Rstudio (fairly recent version required). Can be placed in ~/.Rprofile to run on startup:

setHook(packageEvent("TMB", "onLoad"),
        function(...) {
            if("tools:rstudio" %in% search()) {
                tmb.env <- asNamespace("TMB")
                compile.orig <- tmb.env$compile
                unlockBinding("compile", tmb.env)
                ## Rstudio handle compilation errors:
                rs.env <- as.environment("tools:rstudio")
                tmb.env$compile <- function(file,...) {
                    .Call(rs.env$.rs.routines$rs_sourceCppOnBuild,
                          file, FALSE, FALSE)
                    status <- try( compile.orig(file, ...) )
                    succeeded <- (status == 0)
                    .Call(rs.env$.rs.routines$rs_sourceCppOnBuildComplete,
                          succeeded, "")
                    if(!succeeded) stop("Compilation failed")
                    status
                }
                ## Bind 'sourceCpp' button to TMB compile
                rcpp.env <- asNamespace("Rcpp")
                unlockBinding("sourceCpp", rcpp.env)
                rcpp.env$sourceCpp <- tmb.env$compile
            }
        } )

What works:

  • Tested in Linux and Windows with RStudio version 0.99.902.
  • 'sourceCpp' button works and error messages are displayed nicely in separate window (however note that getwd() must be the containing folder of the cpp file).
  • Works in conjunction with precompile().
  • Also works in conjunction with add-on package 'TMBdebug'.
  • No side effects if not running RStudio or not running TMB.

Next snippet...

Clone this wiki locally