From b4901b1da8b31d01f4b2fcfb34bfeeee58b7afa1 Mon Sep 17 00:00:00 2001 From: Jim Ianelli Date: Wed, 14 Feb 2024 21:02:06 -0800 Subject: [PATCH] Had issue w/ remotes install (probably cuz I broke GOA pollock) --- content/AFSC-BSAI_AtkaMackerel.qmd | 203 +---------------------------- content/AFSC-GOA-pollock.qmd | 7 +- 2 files changed, 8 insertions(+), 202 deletions(-) diff --git a/content/AFSC-BSAI_AtkaMackerel.qmd b/content/AFSC-BSAI_AtkaMackerel.qmd index 9e3ad4c..9ba2951 100644 --- a/content/AFSC-BSAI_AtkaMackerel.qmd +++ b/content/AFSC-BSAI_AtkaMackerel.qmd @@ -333,8 +333,9 @@ obj <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE) # fitting the model -opt <- nlminb(start=obj$par, objective=obj$fn, gradient=obj$gr, - control = list(eval.max = 8000, iter.max = 800)) +#opt <- nlminb(start=obj$par, objective=obj$fn, gradient=obj$gr, +# control = list(eval.max = 8000, iter.max = 800)) + # method = "BFGS", # control = list(maxit=1000000, reltol = 1e-15)) #print(opt) @@ -346,200 +347,8 @@ opt <- nlminb(start=obj$par, objective=obj$fn, gradient=obj$gr, #max(abs(obj$gr())) -sdr <- TMB::sdreport(obj) -sdr_fixed <- summary(sdr, "fixed") -report <- obj$report() +#sdr <- TMB::sdreport(obj) +#sdr_fixed <- summary(sdr, "fixed") +#report <- obj$report() ### Plotting - -mycols <- c("FIMS" = "blue", "ASAP" = "red", "ASAP_orig" = "darkgreen") - -for (i in 1:rdat$parms$nindices){ - index_results <- data.frame( - survey = i, - year = years, - observed = rdat$index.obs[[i]], - FIMS = report$exp_index[[rdat$parms$nfleet+i]], - ASAP = rdat$index.pred[[i]] - ) - if (i==1){ - allinds_results <- index_results - }else{ - allinds_results <- rbind(allinds_results, index_results) - } -} -#print(allinds_results) - -comp_index <- ggplot(allinds_results, aes(x = year, y = observed)) + - geom_point() + - geom_line(aes(x = year, y = FIMS), color = "blue") + - geom_line(aes(x = year, y = ASAP), color = "red") + - facet_wrap(~survey, scales = "free_y", nrow = 2) + - xlab("Year") + - ylab("Index") + - ggtitle("Blue=FIMS, Red=ASAP") + - theme_bw() -#print(comp_index) - -catch_results <- data.frame( - observed = fishing_fleet_index$index_data, - FIMS = report$exp_index[[1]], - ASAP = as.numeric(rdat$catch.pred[1,]) -) -#print(catch_results) - -comp_catch <- ggplot(catch_results, aes(x = years, y = observed)) + - geom_point() + - xlab("Year") + - ylab("Catch (mt)") + - geom_line(aes(x = years, y = FIMS), color = "blue") + - geom_line(aes(x = years, y = ASAP), color = "red") + - ggtitle("Blue=FIMS, Red=ASAP") + - theme_bw() -#print(comp_catch) - -pop_results <- data.frame( - Year = c(years, max(years)+1, years, years, years, years, max(years)+1, years), - Metric = c(rep("SSB", 2*nyears+1), rep("F_mort", 2*nyears), rep("Recruitment", 2*nyears+1)), - Model = c(rep("FIMS", nyears+1), rep("ASAP", nyears), rep(c("FIMS", "ASAP"), each=nyears), - rep("FIMS", nyears+1), rep("ASAP", nyears)), - Value = c(report$ssb[[1]], rdat$SSB, report$F_mort[[1]], rdat$F.report, report$recruitment[[1]], as.numeric(rdat$N.age[,1])) -) -#print(pop_results) - -# ggplot(filter(pop_results, Year <=2019), aes(x=Year, y=Value, color=Model)) + -# geom_line() + -# facet_wrap(~Metric, ncol=1, scales = "free_y") + -# theme_bw() + -# scale_color_manual(values = mycols) - -orig_years <- seq(orig$parms$styr, orig$parms$endyr) -orig_pop_results <- data.frame( - Year = rep(orig_years, 3), - Metric = rep(c("SSB", "F_mort", "Recruitment"), each = length(orig_years)), - Model = "ASAP_orig", - Value = c(orig$SSB, orig$F.report, as.numeric(orig$N.age[,1])) -) - -pop_results_3 <- rbind(pop_results, orig_pop_results) -#print(pop_results_3) - -# ggplot(filter(pop_results_3, Year <=2019), aes(x=Year, y=Value, color=Model)) + -# geom_line() + -# facet_wrap(~Metric, ncol=1, scales = "free_y") + -# theme_bw() + -# scale_color_manual(values = mycols) - -comp_FRSSB3 <- ggplot(pop_results_3, aes(x=Year, y=Value, color=Model)) + - geom_line() + - facet_wrap(~Metric, ncol=1, scales = "free_y") + - theme_bw() + - scale_color_manual(values = mycols) -#print(comp_FRSSB3) - -FIMS_naa_results <- data.frame( - Year = rep(c(years, max(years)+1), each = nages), - Age = rep(ages, nyears+1), - Metric = "NAA", - Model = "FIMS", - Value = report$naa[[1]] -) - -ASAP_naa_results <- data.frame( - Year = rep(years, each = nages), - Age = rep(ages, nyears), - Metric = "NAA", - Model = "ASAP", - Value = as.numeric(t(rdat$N.age)) -) - -orig_naa_results <- data.frame( - Year = rep(orig_years, each = nages), - Age = rep(ages, length(orig_years)), - Metric = "NAA", - Model = "ASAP_orig", - Value = as.numeric(t(orig$N.age)) -) -naa_results <- rbind(FIMS_naa_results, ASAP_naa_results, orig_naa_results) -#print(naa_results) - -# ggplot(filter(naa_results, Year <= 2019), aes(x=Year, y=Value, color=Model)) + -# geom_line() + -# facet_wrap(~Age, ncol=1, scales = "free_y") + -# ylab("NAA") + -# theme_bw() + -# scale_color_manual(values = mycols) - -comp_naa2 <- ggplot(filter(naa_results, Year <= 2019, Model %in% c("ASAP", "FIMS")), aes(x=Year, y=Value, color=Model)) + - geom_line() + - facet_wrap(~Age, ncol=1, scales = "free_y") + - ylab("NAA") + - theme_bw() + - scale_color_manual(values = mycols) -#print(comp_naa2) - -# ggplot(filter(naa_results, Year == 1973, Model %in% c("ASAP", "FIMS")), aes(x=Age, y=Value, color=Model)) + -# geom_line() + -# ylab("NAA in Year 1") + -# theme_bw() + -# scale_color_manual(values = mycols) - - -saveplots <- TRUE -if(saveplots){ - ggsave(filename = "figures/NEFSC_YT_compare_index.png", plot = comp_index, width = 4, height = 4, units = "in") - ggsave(filename = "figures/NEFSC_YT_compare_catch.png", plot = comp_catch, width = 4, height = 4, units = "in") - ggsave(filename = "figures/NEFSC_YT_compare_FRSSB3.png", plot = comp_FRSSB3, width = 5, height = 6.5, units = "in") - ggsave(filename = "figures/NEFSC_YT_compare_NAA2.png", plot = comp_naa2, width = 5, height = 6.5, units = "in") -} - -``` -## Comparison figures -![Catch](figures/NEFSC_YT_compare_catch.png){width=4in} -![Indices](figures/NEFSC_YT_compare_index.png){width=4in} -![FRSSB](figures/NEFSC_YT_compare_FRSSB3.png){width=4in} -![NAA](figures/NEFSC_YT_compare_NAA2.png){width=4in} - -## Comparison table - -The likelihood components from FIMS and ASAP for the same data are shown in the table below. Note that the ASAP file had to turn on the use likelihood constants option to enable this comparison (this option should not be used when recruitment deviations are estimated). - -```{r} -jnlltab <- data.frame(Component=c("Total","Index","Age Comp", "Rec"), - FIMS = c(report$jnll, report$index_nll, report$age_comp_nll, report$rec_nll), - ASAP = c(rdat$like$lk.total, - (rdat$like$lk.catch.total + rdat$like$lk.index.fit.total), - (rdat$like$lk.catch.age.comp + rdat$like$lk.index.age.comp), - rdat$like$lk.Recruit.devs)) -print(jnlltab) -``` - -## What was your experience using FIMS? What could we do to improve usability? - -Relatively easy to use by following the vignette. Creating wrappers for data input would help so that each element did not need to be assigned directly. - -## List any issues that you ran into or found - -Please [open an issue](https://github.com/NOAA-FIMS/FIMS/issues/new/choose) if you found something new. - -* SSB calculations in FIMS assume 0.5 multiplier, which differs from ASAP [Issue #521](https://github.com/NOAA-FIMS/FIMS/issues/521). -* Output all derived values (this is mostly done) -* Fix recruitment estimation [Issue #364](https://github.com/NOAA-FIMS/FIMS/issues/364) -* Handle missing data, especially surveys [Issue #502](https://github.com/NOAA-FIMS/FIMS/issues/502) -* Weights at age that change over time -* Separate weights at age for catch, SSB, Jan-1 population, indices, etc. -* Fishery selectivity blocks or random effects -* Allow time-varying CVs and ESS (or alternative functions) -* Option for Index in numbers -* Timing of Index and SSB calculations within the year -* One-step-ahead residuals -* Reference points, projections, pushbutton retro - -## What features are most important to add based on this case study? - -* Missing values, would allow inclusion of the other 3 indices (too many missing years to fill for this example) - -```{r} -# Clear C++ objects from memory -clear() -``` diff --git a/content/AFSC-GOA-pollock.qmd b/content/AFSC-GOA-pollock.qmd index 374b1d2..a35a0b6 100644 --- a/content/AFSC-GOA-pollock.qmd +++ b/content/AFSC-GOA-pollock.qmd @@ -18,12 +18,9 @@ installed_packages <- packages %in% rownames(installed.packages()) if (any(installed_packages == FALSE)) { install.packages(packages[!installed_packages], repos = "http://cran.us.r-project.org") } -<<<<<<< HEAD -remotes::install_github("NOAA-FIMS/FIMS") -======= +#remotes::install_github("NOAA-FIMS/FIMS") #remotes::install_github("NOAA-FIMS/FIMS", ref='isNA') ->>>>>>> 84cf899 (Jim added an initial atka case) -remotes::install_github("kaskr/TMB_contrib_R/TMBhelper") +#remotes::install_github("kaskr/TMB_contrib_R/TMBhelper") library(ggplot2) library(tidyr)