From 468d943e1fcbaed8bbb5d62339f5a6d922e10df6 Mon Sep 17 00:00:00 2001 From: Jim Ianelli Date: Thu, 15 Jan 2015 02:45:15 -0800 Subject: [PATCH 1/4] added echoes --- src/gmacs.tpl | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/src/gmacs.tpl b/src/gmacs.tpl index 12ee8882..0eb69c1d 100644 --- a/src/gmacs.tpl +++ b/src/gmacs.tpl @@ -1,4 +1,4 @@ -// ==================================================================================== // +/// ==================================================================================== // // Gmacs: A generalized size-structured stock assessment modeling framework. // // Authors: Athol Whitten and Jim Ianelli @@ -136,7 +136,7 @@ DATA_SECTION init_vector size_breaks(1,nclass+1); vector mid_points(1,nclass); !! mid_points = size_breaks(1,nclass) + 0.5 * first_difference(size_breaks); - !! ECHO(syr); ECHO(nyr); ECHO(nfleet); ECHO(nsex); ECHO(nshell);ECHO(nmature); ECHO(nclass); + !! ECHO(syr); ECHO(nyr); ECHO(nfleet); ECHO(nsex); ECHO(nshell);ECHO(nmature); ECHO(nclass); ECHO(size_breaks); // |-----------| // | ALLOMETRY | @@ -157,6 +157,7 @@ DATA_SECTION // |-------------------------------| init_vector fecundity(1,nclass); init_matrix maturity(1,nsex,1,nclass); + !! ECHO(fecundity); ECHO(maturity); // |-------------| // | FLEET NAMES | @@ -182,6 +183,7 @@ DATA_SECTION catch_cv(k) = column(dCatchData(k),6); catch_dm(k) = column(dCatchData(k),11); } + ECHO(nCatchDF); ECHO(nCatchRows); ECHO(dCatchData); END_CALCS //!! ECHO(obs_catch); ECHO(catch_cv); @@ -244,7 +246,7 @@ DATA_SECTION obs_cpue(k) = column(dSurveyData(k),5); cpue_cv(k) = column(dSurveyData(k),6); } - ECHO(obs_cpue); ECHO(cpue_cv); + ECHO(nSurveys);ECHO(nSurveyRows),ECHO(dSurveyData); ECHO(obs_cpue); ECHO(cpue_cv); END_CALCS @@ -263,12 +265,12 @@ DATA_SECTION { dmatrix tmp = trans(d3_SizeComps(k)).sub(1,nSizeCompCols(k)); d3_obs_size_comps(k) = trans(tmp); + // NOTE This normalizes all observations by row--may be incorrect if shell condition for (int i=1;i<=nSizeCompRows(k);i++) d3_obs_size_comps(k,i) /= sum(d3_obs_size_comps(k,i)); size_comp_sample_size(k) = column(d3_SizeComps(k),0); } - ECHO(nSizeComps); - ECHO(d3_obs_size_comps); + ECHO(nSizeComps);ECHO(nSizeCompRows); ECHO(nSizeCompCols); ECHO(d3_SizeComps); ECHO(d3_obs_size_comps); END_CALCS ivector ilike_vector(1,nlikes) LOC_CALCS @@ -338,11 +340,7 @@ DATA_SECTION } } - - ECHO(dPreMoltSize); - ECHO(iMoltIncSex); - ECHO(dMoltInc); - ECHO(dMoltIncCV); + ECHO(nGrowthObs); ECHO(dGrowthData); ECHO(dPreMoltSize); ECHO(iMoltIncSex); ECHO(dMoltInc); ECHO(dMoltIncCV); END_CALCS // |------------------| @@ -402,6 +400,7 @@ DATA_SECTION Grwth_lb = column(Grwth_control,2); Grwth_ub = column(Grwth_control,3); Grwth_phz = ivector(column(Grwth_control,4)); + ECHO(theta_control); ECHO(Grwth_control); END_CALCS @@ -419,6 +418,7 @@ DATA_SECTION init_imatrix slx_nret(1,nsex,1,nfleet); init_matrix slx_control(1,nslx,1,nc); + !! ECHO(slx_nsel_blocks); ECHO(slx_nret); ECHO(slx_control); ivector slx_indx(1,nslx); ivector slx_type(1,nslx); @@ -493,9 +493,7 @@ DATA_SECTION prior_qtype = column(q_controls,1); prior_qbar = column(q_controls,2); prior_qsd = column(q_controls,3); - ECHO(prior_qtype); - ECHO(prior_qbar); - ECHO(prior_qsd); + ECHO(q_controls); ECHO(prior_qtype); ECHO(prior_qbar); ECHO(prior_qsd); END_CALCS @@ -526,6 +524,7 @@ DATA_SECTION } } } + ECHO(f_controls); ECHO(f_phz); END_CALCS @@ -535,6 +534,7 @@ DATA_SECTION init_ivector nAgeCompType(1,nSizeComps); init_ivector bTailCompression(1,nSizeComps); init_ivector nvn_phz(1,nSizeComps); + !! ECHO(nAgeCompType); ECHO(bTailCompression); ECHO(nvn_phz); @@ -564,6 +564,7 @@ DATA_SECTION nMdev = m_nNodes; break; } + ECHO(m_type); ECHO(Mdev_phz); ECHO(m_stdev); ECHO(m_nNodes); ECHO(m_nodeyear); END_CALCS @@ -590,11 +591,11 @@ DATA_SECTION spr_fleet = int(model_controls(7)); spr_lambda = model_controls(8); bUseEmpiricalGrowth = int(model_controls(9)); + ECHO(model_controls); END_CALCS init_int eof_ctl; - !! if(eof_ctl!=9999){cout<<"Error reading control file"< Date: Thu, 15 Jan 2015 09:59:48 -0800 Subject: [PATCH 2/4] Added a check in get-recruitment.R for when Hessian does not converge. Checks to see if sd_log_recruits exist within replist$fit, otherwise one cannot get recruitment. --- Rsrc/R/get-recruitment.R | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Rsrc/R/get-recruitment.R b/Rsrc/R/get-recruitment.R index 8c9191fb..9e6f4196 100644 --- a/Rsrc/R/get-recruitment.R +++ b/Rsrc/R/get-recruitment.R @@ -6,6 +6,11 @@ #' @export get_recruitment <- function(replist){ A <- replist + if(length(A$fit$est[A$fit$names == "sd_log_recruits"]) == 0) { + stop("Appears that the Hessian was not positive definite\n + thus estimates of recruitment do not exist.\n + See this in replist$fit.") + } dfpar <- data.frame(par=A$fit$names,log_rec=A$fit$est,log_sd=A$fit$std) df <- subset(dfpar,par=="sd_log_recruits")[,-1] df$year <- A$mod_yrs From 947ae8ed7a891bfede822695ca01a207f7d2e269 Mon Sep 17 00:00:00 2001 From: kellijohnson Date: Thu, 15 Jan 2015 10:56:29 -0800 Subject: [PATCH 3/4] change check in get-recuruitment --- Rsrc/R/get-recruitment.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Rsrc/R/get-recruitment.R b/Rsrc/R/get-recruitment.R index 9e6f4196..35a66733 100644 --- a/Rsrc/R/get-recruitment.R +++ b/Rsrc/R/get-recruitment.R @@ -6,7 +6,7 @@ #' @export get_recruitment <- function(replist){ A <- replist - if(length(A$fit$est[A$fit$names == "sd_log_recruits"]) == 0) { + if(is.null(A$fit$logDetHess)) { stop("Appears that the Hessian was not positive definite\n thus estimates of recruitment do not exist.\n See this in replist$fit.") From f5b29560ecebd7c999038185a9f14247f2380afd Mon Sep 17 00:00:00 2001 From: kellijohnson Date: Thu, 15 Jan 2015 11:52:02 -0800 Subject: [PATCH 4/4] Add function to plot multiple recruitments across models --- Rsrc/R/plot-recruitment.R | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/Rsrc/R/plot-recruitment.R b/Rsrc/R/plot-recruitment.R index f1bb80a1..8a2fc494 100644 --- a/Rsrc/R/plot-recruitment.R +++ b/Rsrc/R/plot-recruitment.R @@ -14,3 +14,23 @@ plot_recruitment <- function(replist){ pRecruitment <- p + ggtheme return(pRecruitment) } + +#' Plot predicted recruitment across model runs +#' +#' +#' @param data A list of multiple objects created by read_admb function +#' @param modnames A vector of model names included in \code{data} +#' @return Plot of predicted recruitment compared across models +#' @author Cole Monnahan Kelli Johnson +#' @export +plot_models_recruitment <- function(data, modnames){ + recs <- lapply(data, get_recruitment) + df <- do.call("rbind", Map(cbind, recs, modname = modnames)) + + p <- ggplot(df,aes(x=factor(year),y=exp(log_rec), group=modname, colour=modname)) + p <- p + geom_line(stat = "identity", alpha=0.4) + p <- p + geom_pointrange(aes(factor(year),exp(log_rec),ymax=ub,ymin=lb)) + p <- p + labs(x="Year", y="Recruitment") + pRecruitment <- p + ggtheme + return(pRecruitment) +} \ No newline at end of file