cole-monnahan-noaa / adnuts Goto Github PK

We've reorganized the Fisheries Integrated Toolbox GitHub Organization, and by doing so the .css file that makes your pkgdown site NOAA colors has been changed.

To update your pkgdown site, edit your extra.css to @import url("https://nmfs-fish-tools.github.io/nmfspalette/extra.css");

Here's an example from r4MAS.

Hi just noticed that sample_tmb is now not exported and the file is called sample_tmb_deprecated.r just wondering if there is an alternate place for accessing this now?
If so it may be worth adding a deprecated message using .deprecated for it telling people where to look.

Cheers

@colemonnahan, thank you for this excellent tool.

@okenk and I have been exploring {adnuts} for some Stock Synthesis models and she asked about converting between the parameter labels used by ADMB which appear in the adnuts output and the parameter labels internal to Stock Synthesis.

The following commands worked for me to get new labels into the plots, but there are surely better ways to go about this. Can you provide any guidance on the best way to update the parameter labels used by {adnuts}?

# get informative parameter labels used by Stock Synthesis
# model <- r4ss::SS_output() # read model output
est_par_labels <- model$parameters %>% 
  dplyr::filter(!is.na(Active_Cnt)) %>%
  dplyr::pull(Label)

# copy adnuts results to avoid messing things up
fit2 <- fit

# fit$samples is a 3D array where the third dimension is the parameter value
# The "lp__" column seems to be the log-posterior added by adnuts.
dimnames(fit2$samples)[[3]] <- c(est_par_labels, "lp__")

fit2$par_names = est_par_labels # not sure if this needs replacement or not

# for my model samples_unbounded = NULL, so I'm not sure about the following step
if (!is.null(fit2$samples_unbounded)) {
  dimnames(fit2$samples_unbounded)[[3]] <- c(est_par_labels, "lp__")
}

# make some plots
adnuts::pairs_admb(fit2, pars=1:6, order='slow', label.cex = 0.6)

# explore output using ShinyStan
adnuts::launch_shinyadmb(fit2)

The function .check_ADMB_version crashes if the model is compiled with the fast flag -f. This needs to be robust to that flag, and probably also try to catch if the model is ADMB-RE and gracefully exit early b/c it won't necessarily in ADMB.

check_identifiable() crashes when eigenvalues of a bad "covariance" matrix are complex. The following code handles the situation without blowing up. Might want to add an error message pointing out why the check failed.

check_identifiable <- function(model, path=getwd()){

Check eigendecomposition

fit <- adnuts:::.read_mle_fit(model, path)
hes <- adnuts:::.getADMBHessian(path)
ev <- eigen(hes);
if (any(is.complex(ev$values))){
WhichBad <- which( abs(Im(ev$values)) > .Machine$double.eps);
} else {
WhichBad <- which( ev$values < sqrt(.Machine$double.eps) );
}

if(length(WhichBad)==0){
message( "All parameters are identifiable" )
} else {
## Check for parameters
if(length(WhichBad)==1){
RowMax <- abs(ev$vectors[,WhichBad])
} else {
RowMax <- apply(ev$vectors[, WhichBad], MARGIN=1, FUN=function(vec){max(Mod(vec))} )
}
bad <- data.frame(ParNum=1:nrow(hes), Param=fit$par.names,
MLE=fit$est[1:nrow(hes)],
Param_check=ifelse(RowMax>0.1, "Bad","OK"))
row.names(bad) <- NULL
bad <- bad[bad$Param_check=='Bad',]
print(bad)
return(invisible(bad))
}
}

Running adnuts version v1.1.0 with Stock Synthesis (SS) version 3.30.16.03 I get the following error after merging chains from individual cores and trying to run the -mceval.

"SS is not configured to work with -noest; use -maxI instead which overrides maxphase in starter.ss"

Looks like for SS applications this is referring to the system call on Line 378 of sample_admb.R

It was easy enough to just run the -mceval via the command line to skirt this, but figured a warning may help direct folks or a code alteration may remedy this.

It would be helpful to have an easy way to check a fitted object for any issues, like divergences, exceeding max_treedepth, BFMI, Rhat, ESS, etc. Basically to adapt this code to work with adnuts.

Open to other ideas to check as well. Could look at unbounded variances (want these ~1) or differences in MLE and posterior estimates?

There are two entries for the adapt_delta element of the control object

• https://github.com/Cole-Monnahan-NOAA/adnuts/blob/main/R/sample_admb.R#L104
• https://github.com/Cole-Monnahan-NOAA/adnuts/blob/main/R/sample_admb.R#L108

which imply a difference in usage and class. Which one is correct?

The function .read_mle_fit searches for a par file using the line
ff <- list.files()[grep(x = list.files(), pattern = ".par")]

This does not actually search for files containing ".par" because of how R does regular expressions. Instead, it returns all files containing "par," including, say "plot_comparisons.png."

If you change the pattern argument from .par to \\.par it will actually search for ".par".

The user needs to use model='./model' instead of just model='model' when using a MacOS. I handled this for the MCMC part but it was not done consistently through with the checks and other aspects. Ideally the user specifies just "model" and the "/." part is added internally. This makes code OS-agnostic and cleaner on the user end.

I tried to fix this in a new branch and would like some help testing it (I don't have access to a MacOS machine). @wStockhausen would you please try this? Or maybe @jimianelli?

devtools::install_github('Cole-Monnahan-NOAA/adnuts', ref='fix_path')

Then test this with code like

fit <- sample_rwm(model=m, path=p, iter=100, chains=1)

Where the variables m and p vary. Basically try to break it. It should gracefully and informatively error out when either is wrong. Also what happens if you pass e.g., m='./simple' and it gets added on again internally? That works on Windows.

Recently @okenk and @jimianelli both mentioned they would like to have better options for a workflow. I'm open to that. Please use this issue to discuss what your preferred workflow is and how it could be implemented in adnuts.

Right now the user specifies the path to a directory that contains all the necessary folders, and the model name. It's assumed the executable is in that folder. Each is then copied to the working directory when running in parallel.

ff <- file.path(path, paste0(model, ".exe"))

windoze doesn't need to have this extension (I don't think) and it breaks models in bsd/linux so why have it...

Consider changing the line in the examples for sample from:

init <- function() rnorm(2)

to

init <- function() rnorm(numpars)

where numpars would be set based on the problem.

This is so that in case someone misses that the example is a 2-parameter regression, it would be more extensible. Also init= not really defined elsewhere (at least with clear instructions for alternatives from NULL.

Hi,

I have been trying the sample_tmb function on an 'obj' that includes random parameters (I see in your example that you use "random = NULL"). I get the following warning message:

In sample_tmb(obj = obj2, iter = 40, chains = 1, init = initran, :
Some parameters declated as random. Are you sure? For MCMC this is usually turned off

It takes about 9h for me to finish the sampling (for only 40 iterations).

My question for now is if sample_tmb can be used with the random parameters?

PS: I can upload an example if needed.

Thanks
Janne

Couldn't find in documentation and too dumb in R to figure out how to get default values

When running in a parallel, the dominate usage case, the output does not appear in the Rstudio console. A workaround is to run it in the R gui or terminal directly, or use a text editor where it shows.

According to this thread it's not unexpected and possibly switching to the future package could be a good solution.

https://stackoverflow.com/questions/62308162/why-dont-parallel-jobs-print-in-rstudio

If anyone has an experience with this package or getting console output to work for RStudio please let me know.

The warning comes from the low-level function which makes no sense. Should come from top-level sampling calls.