Code snippets

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!!!)

template<class Type>
Type posfun(Type x, Type eps, Type &pen)
{
  if ( x >= eps ){
    return x;
  } else {
    pen += Type(0.01) * pow(x-eps,2);
    return eps/(Type(2.0)-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()

---

Next snippet...