search

commfish / seak_sablefish

NSEI sablefish stock assessment
8 stars 5 forks source link

Fitting model with TMBhelper::fit_mod() and with random effects not working #51

Open jysullivan opened 4 years ago

jysullivan commented 4 years ago

Debugging TMBhelper's fit_mod() function, which currently doesn't work. Also can't figure getting NAN function eval for random effects even though one can estimate sigmaR when rec_devs/rinit_devs are treated as fixed effects. Tried to fix this before I left but didn't have time to finish it. Please reach out to ummjane@gmail.com for help with this feature. Here's some code to get started.

# # Using TMBhelper's fit_tmb():
# TMB::compile("mod.cpp")
# dyn.load(TMB::dynlib("mod"))
# 
# # work out the map for this phase if phases for parameters is less than the
# # current phase then map will contain a factor filled with NAs
# map_use <- list()
# map_use$dummy <- fill_vals(parameters$dummy, NA)
# 
# # if not using random effects, assign log_sigma_r an NA in the map so it's not estimated
# if (data$random_rec == FALSE) {
#   map_use$log_sigma_r <- fill_vals(parameters$log_sigma_r, NA)
# }
# 
# # if natural mortality is fixed, assign log_M an NA in the map so its not
# # estimated
# if (data$M_type == 0) {
#   map_use$log_M <- fill_vals(parameters$log_M, NA)
# }
# 
# # if not using the Dirichlet-multinonial, assign log_fsh_theta and
# # log_srv_theta NAs in the map so they're not estimated
# if (data$comp_type != 1) {
#   map_use$log_fsh_theta <- fill_vals(parameters$log_fsh_theta, NA)
#   map_use$log_srv_theta <- fill_vals(parameters$log_srv_theta, NA)
# }
# 
# # Temporary debug trying to figure out why I'm getting NA/NaN function
# # evaluation
# if (tmp_debug == TRUE) {
#   # map_use$log_spr_Fxx <- fill_vals(parameters$log_spr_Fxx, NA)
#   map_use$log_fsh_slx_pars <- fill_vals(parameters$log_fsh_slx_pars, NA)
#   map_use$log_srv_slx_pars <- fill_vals(parameters$log_srv_slx_pars, NA)
#   # map_use$fsh_logq <- fill_vals(parameters$fsh_logq, NA)
# }
# # 
# # Build upper and lower parameter bounds and remove any that are not
# # estimated (should be the inverse of the map_use)
# bounds <- build_bounds(param_list = parameters)
# bounds$upper <- bounds$upper[!names(bounds$upper) %in% names(map_use)]
# bounds$lower <- bounds$lower[!names(bounds$lower) %in% names(map_use)]
# 
# # Remove inactive parameters from bounds and vectorize
# lower <- unlist(bounds$lower)
# upper <- unlist(bounds$upper)
# 
# # if (data$random_rec == TRUE) {
# #   lower <- lower[-which(grepl(random_vars[1], names(lower)))]
# #   lower <- lower[-which(grepl(random_vars[2], names(lower)))]
# #   upper <- upper[-which(grepl(random_vars[1], names(upper)))]
# #   upper <- upper[-which(grepl(random_vars[2], names(upper)))]
# # }
# # 
# obj <- TMB::MakeADFun(data, parameters, random = NULL,#random_vars, 
#                       DLL = "mod", map = map_use)
# 
# TMBhelper::fit_tmb(obj = obj, fn=obj$fn, gr=obj$gr, startpar=NULL, lower=lower, upper=upper,
#                    getsd=TRUE, control=list(eval.max=1e4, iter.max=1e4, trace=0), bias.correct=FALSE,
#                    bias.correct.control=list(sd=FALSE, split=NULL, nsplit=NULL, vars_to_correct=NULL),
#                    savedir=NULL, loopnum=1, newtonsteps=0, n=Inf, getReportCovariance=FALSE, 
#                    getJointPrecision=FALSE,
#                    getHessian=FALSE, quiet=FALSE)