diff --git a/_quarto.yml b/_quarto.yml
index c0806ef..7194f76 100644
--- a/_quarto.yml
+++ b/_quarto.yml
@@ -43,6 +43,8 @@ website:
         text: R Markdown files
       - href: content/NEFSC-yellowtail.qmd
         text: NEFSC yellowtail flounder case study
+      - href: content/AFSC-GOA-pollock.qmd
+        text: AFSC GOA pollock case study
       - href: content/acknowledgements.qmd
         text: Acknowledgements
 
diff --git a/content/AFSC-GOA-pollock.qmd b/content/AFSC-GOA-pollock.qmd
new file mode 100644
index 0000000..33ac3c6
--- /dev/null
+++ b/content/AFSC-GOA-pollock.qmd
@@ -0,0 +1,277 @@
+---
+title: AFSC Case Study Gulf of Alaska Walleye Pollock
+format:
+  html:
+    code-fold: true
+---
+
+## The setup
+
+```{r}
+#| warning: false
+#| label: startup
+
+library(tidyverse)
+library(FIMS)
+library(TMB)
+library(TMBhelper)
+## devtools::install_github('afsc-assessments/GOApollock', ref='v0.1.2')
+library(GOApollock)
+theme_set(theme_bw())
+
+
+R_version <- version$version.string
+TMB_version <- packageDescription("TMB")$Version
+FIMS_commit <- substr(packageDescription("FIMS")$GithubSHA1, 1, 7)
+```
+
+-   R version: `r R_version`\
+-   TMB version: `r TMB_version`\
+-   FIMS commit: `r FIMS_commit`\
+-   Stock name: Gulf of Alaska (GOA) Walleye Pollock\
+-   Region: AFSC\
+-   Analyst: Cole Monnahan\
+
+## Simplifications to the original assessment
+
+The model presented in this case study was changed substantially from the operational version and should not be considered reflective of the pollock stock. This is intended as a demonstration and nothing more.
+
+To get the opertional model to more closely match FIMS I:
+
+* Droped surveys 1, 4, and 5
+* Removed ageing error
+* Removed length compositions
+* Removed bias correction in log-normal index likelihoods
+* Simplified catchability of survey 3 to be constant in time (removed random walk)
+* Updated maturity to be parametric rather than empirical
+* Used constant weight at age for all sources: spawning, fishery, surveys, and biomass calculations. The same matrix was used throughout.
+* Changee timing to be Jan 1 for spawning and all surveys
+* Removed prior on catchability for survey 2
+* Removed time-varying fisheries selectivity (constant double logistic)
+* Took off normalization of selectivity
+* Removed age accumulation for fishery age compositions
+
+
+## Script to prepare data for building FIMS object
+```{r}
+#| warning: false
+#| label: prepare-fims-data
+## define the dimensions and global variables
+years <- 1970:2023
+nyears <- length(years)
+nseasons <- 1
+nages <- 10
+ages <- 1:nages 
+## nfleets <- 1
+## This will fit the models bridging to FIMS (simplifying)
+## source("fit_bridge_models.R")
+## compare changes to model
+pkfitfinal <- readRDS("data_files/pkfitfinal.RDS")
+pkfit0 <- readRDS("data_files/pkfit0.RDS")
+parfinal <- pkfitfinal$obj$env$parList()
+pkinput0 <- readRDS('data_files/pkinput0.RDS')
+fimsdat <- pkdat0 <- pkinput0$dat
+pkinput <- readRDS('data_files/pkinput.RDS')
+```
+
+## Run FIMS model
+```{r}
+#| warning: false
+#| label: run-FIMS
+#| output: false
+#| eval: true
+
+## set up FIMS data objects
+clear()
+clear_logs()
+estimate_fish_selex <- TRUE
+estimate_survey_selex <- TRUE
+estimate_q2 <- TRUE
+estimate_q3 <- TRUE
+estimate_q6 <- TRUE
+estimate_F <- TRUE
+estimate_recdevs <- TRUE
+source("pk_prepare_FIMS_inputs.R")
+## make FIMS model
+success <- CreateTMBModel()
+parameters <- list(p = get_fixed())
+obj <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
+## report values for the two models
+rep0 <- pkfitfinal$rep
+rep1 <- obj$report() # FIMS initial values
+## try fitting the model
+opt <- TMBhelper::fit_tmb(obj, getsd=FALSE, newtonsteps=0, control=list(trace=100))
+## opt <- with(obj, nlminb(start=par, objective=fn, gradient=gr))
+max(abs(obj$gr())) # from Cole, can use TMBhelper::fit_tmb to get val to <1e-10
+rep2 <- obj$report() ## FIMS after estimation
+
+## Output plotting
+out1 <- get_long_outputs(rep1, rep0) %>%
+  mutate(platform=ifelse(platform=='FIMS', 'FIMS init', 'TMB'))
+out2 <- get_long_outputs(rep2, rep0) %>% filter(platform=='FIMS') %>%
+  mutate(platform='FIMS opt')
+out <- rbind(out1,out2)
+g <- ggplot(out, aes(year, value, color=platform)) + geom_line() +
+  facet_wrap('name', scales='free') + ylim(0,NA) +
+  labs(x=NULL, y=NULL)
+ggsave('figures/AFSC_PK_ts_comparison.png', g, width=9, height=5, units='in')
+g <- ggplot(filter(out, platform!='TMB'), aes(year, relerror, color=platform)) + geom_line() +
+  facet_wrap('name', scales='free') +
+  labs(x=NULL, y='Relative difference') + coord_cartesian(ylim=c(-.5,.5))
+ggsave('figures/AFSC_PK_ts_comparison_relerror.png', g, width=9, height=5, units='in')
+
+## Quick check on age comp fits
+p1 <- get_acomp_fits(rep0, rep1, rep2, fleet=1, years=pkdat0$fshyrs)
+g <- ggplot(p1, aes(age, paa, color=platform)) + facet_wrap('year') + geom_line()
+ggsave('figures/AFSC_PK_age_comp_fits_1.png', g, width=9, height=8, units='in')
+p2 <- get_acomp_fits(rep0, rep1, rep2, fleet=2, years=pkdat0$srv_acyrs2)
+g <- ggplot(p2, aes(age, paa, color=platform)) + facet_wrap('year') + geom_line()
+ggsave('figures/AFSC_PK_age_comp_fits_2.png', g, width=9, height=8, units='in')
+## p3 <- get_acomp_fits(rep0, rep1, rep2, fleet=3, years=pkdat0$srv_acyrs3)
+## g <- ggplot(p3, aes(age, paa, color=platform)) + facet_wrap('year') + geom_line()
+## p6 <- get_acomp_fits(rep0, rep1, rep2, fleet=4, years=pkdat0$srv_acyrs6)
+## g <- ggplot(p6, aes(age, paa, color=platform)) + facet_wrap('year') + geom_line()
+
+## index fits
+addsegs <- function(yrs, obs, CV){
+  getlwr <- function(obs, CV) qlnorm(p=.025, meanlog=log(obs), sdlog=sqrt(log(1+CV^2)))
+  getupr <- function(obs, CV) qlnorm(p=.975, meanlog=log(obs), sdlog=sqrt(log(1+CV^2)))
+  segments(yrs, y0=getlwr(obs,CV), y1=getupr(obs,CV))
+  points(yrs, obs, pch=22, bg='white')
+}
+png('figures/AFSC_PK_index_fits.png', res=300, width=6, height=7, units='in')
+par(mfrow=c(3,1), mar=c(3,3,.5,.5), mgp=c(1.5,.5,0), tck=-0.02)
+plot(years, rep0$Eindxsurv2, type='l',
+     ylim=c(0,2), lwd=5.5,
+      xlab=NA, ylab='Biomass (million t)')
+x1 <- out1 %>% filter(name=='Index2' & platform=='FIMS init')
+x2 <- out2 %>% filter(name=='Index2' & platform=='FIMS opt')
+lines(years,x1$value, col=2, lwd=1.5)
+lines(years,x2$value, col=3, lwd=1.5)
+addsegs(yrs=pkdat0$srvyrs2, obs=pkdat0$indxsurv2, CV=pkdat0$indxsurv_log_sd2)
+legend('topright', legend=c('TMB', 'FIMS init', 'FIMS opt'), lty=1, col=1:3)
+mtext('Survey 2', line=-1.5)
+plot(years, rep0$Eindxsurv3, type='l',
+     ylim=c(0,.6), lwd=5.5,
+      xlab=NA, ylab='Biomass (million t)')
+x1 <- out1 %>% filter(name=='Index3' & platform=='FIMS init')
+x2 <- out2 %>% filter(name=='Index3' & platform=='FIMS opt')
+lines(years,x1$value, col=2, lwd=1.5)
+lines(years,x2$value, col=3, lwd=1.5)
+addsegs(yrs=pkdat0$srvyrs3, obs=pkdat0$indxsurv3, CV=pkdat0$indxsurv_log_sd3)
+mtext('Survey 3', line=-1.5)
+legend('topright', legend=c('TMB', 'FIMS init', 'FIMS opt'), lty=1, col=1:3)
+plot(years, rep0$Eindxsurv6, type='l',
+     ylim=c(0,2.6), lwd=5.5,
+      xlab=NA, ylab='Biomass (million t)')
+x1 <- out1 %>% filter(name=='Index6' & platform=='FIMS init')
+x2 <- out2 %>% filter(name=='Index6' & platform=='FIMS opt')
+lines(years,x1$value, col=2, lwd=1.5)
+lines(years,x2$value, col=3, lwd=1.5)
+addsegs(yrs=pkdat0$srvyrs6, obs=pkdat0$indxsurv6, CV=pkdat0$indxsurv_log_sd6)
+mtext('Survey 6', line=-1.5)
+legend('topright', legend=c('TMB', 'FIMS init', 'FIMS opt'), lty=1, col=1:3)
+dev.off()
+```
+## Comparison figures for basic model
+![Time Series](figures/AFSC_PK_ts_comparison.png){width=7in}
+![Time Series (relative error)](figures/AFSC_PK_ts_comparison_relerror.png){width=7in}
+![Indices](figures/AFSC_PK_index_fits.png){width=7in}
+![Time Series](figures/AFSC_PK_ts_comparison.png){width=7in}
+![Fishery Age Composition Fits](figures/AFSC_PK_age_comp_fits_1.png){width=7in}
+![Survey 2 Age Composition Fits](figures/AFSC_PK_age_comp_fits_2.png){width=7in}
+
+## Extra analyses
+Two extra analyses are demonstrated. First is a likelihood profile over lnR0, showing component contributions and testing for data conflict (a Piner plot). The second is to run the model through the 'Stan' software using the 'tmbstan' R package. This samples from the posterior, which are put back into the model to get the posterior distribution for spawning stock biomass. Given its long run time the results are saved to a file and read in for post-hoc plotting.
+
+``` {r}
+#| label: extra-calcs
+#| warning: true
+#| output: false
+
+## Try a likelihood profile
+i <- 68 # this will break if model is changed at all
+map <- parameters
+map$p[i] <- NA # map off R0 specified below
+map$p <- as.factor(map$p)
+xseq <- as.numeric(c(opt$par[i],  seq(22,24, len=30)))
+res <- list()
+for(j in 1:length(xseq)){
+  print(j)
+  parameters$p[i] <- xseq[j]
+  obj2 <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE, map=map)
+  opt2 <- TMBhelper::fit_tmb(obj2, getsd=FALSE, newtonsteps=0, control=list(trace=0))
+  out <- obj2$report()
+  res[[j]] <- data.frame(j=j, lnR0=xseq[j], total=out$jnll, index=out$index_nll,
+                         age=out$age_comp_nll,recruit=out$rec_nll, maxgrad=max(abs(obj$gr())))
+}
+res <- bind_rows(res) %>%
+  pivot_longer( cols=c(total, index, age, recruit)) %>%
+  group_by(name) %>% mutate(deltaNLL=value-min(value))
+g <- ggplot(res, aes(lnR0, deltaNLL, color=name)) + geom_line()
+g <- g+ geom_point(data=filter(res, deltaNLL==0), size=2) +
+  labs(y='Delta NLL', color='NLL component')
+ggsave('figures/AFSC_PK_like_profile_R0.png', g, width=7, height=3.5)
+
+
+
+## ## Try Bayesian
+## library(tmbstan)
+## library(shinystan)
+## ## Some paraemters wandering off to Inf so fix those (need
+## ## priors). Needs a ton of work but is proof of concept. Major
+## ## problem is parallel fails.
+## map <- parameters
+## ## parameters$p[65:66]
+## map$p[65:66] <- NA # map off R0 specified below
+## map$p <- as.factor(map$p)
+## obj3 <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent=TRUE, map=map)
+## opt3 <- TMBhelper::fit_tmb(obj3, getsd=FALSE, newtonsteps=0, control=list(trace=0))
+## ## Fails when trying to do this in parallel unfortunately
+## fit <- tmbstan(obj3, chains=1, cores=1, open_progress=FALSE,
+##                init='last.par.best', control=list(max_treedepth=10))
+## launch_shinystan(fit)
+## df <- as.data.frame(fit)
+## df <- df[,-ncol(df)] # drop lp__
+## ## for each posterior draw, report to get SSB
+## postreps <- list()
+## for(ii in 1:nrow(df)){
+##   if(ii %% 10==0) print(ii)
+##   postreps[[ii]] <- obj3$rep(df[ii,])
+## }
+## ssbpost <- lapply(postreps, function(x) data.frame(year=years, ssb=x$ssb[[1]][-55]))%>%
+##   bind_rows %>% mutate(rep=rep(1:nrow(df), each=54))
+## saveRDS(ssbpost, file='pk_SSB_posteriors.RDS')
+ssbpost <- readRDS('data_files/pk_pollock_SSB_posteriors.RDS')
+g <- ggplot(ssbpost, aes(year, ssb/1e9, group=rep)) + geom_line(alpha=.1) +
+  ylim(0,NA) + labs(x=NULL, y='SSB (M t)', title='Posterior demo (unconverged!)')
+ggsave('figures/AFSC_PK_ssb_posterior.png', g, width=7, height=4, units='in')
+
+```
+
+## Plots for extra analyses
+![Likelihood Profile](figures/AFSC_PK_like_profile.png){width=7in}
+![SSB Posterior](figures/AFSC_PK_ssb_posterior.png){width=7in}
+
+## Comparison table
+
+The likelihood components from the TMB model do not include constants and thus are not directly comparable. To be fixed later. Relative differences between the modified TMB model and FIMS implementation are given in the figure above.
+
+## What was your experience using FIMS? What could we do to improve usability?
+
+To do
+
+## List any issues that you ran into or found
+
+* Output more derived quantities like selectivity, maturity, etc.
+* NLL components are not separated by fleet and need to be. So age comp NLL for fleets 1 and 2 need to be separate to make, e.g., the likelihood profile plot above.
+* Need more ADREPORTed values like SSB
+
+## What features are most important to add based on this case study?
+
+* More sophisticated control over selectivity so that ages 1 and 2 can be zeroed out for a double-logistic form, overriding the selectivity curve.
+
+```{r}
+# Clear C++ objects from memory
+clear()
+```
diff --git a/content/data_files/pk_pollock_SSB_posteriors.RDS b/content/data_files/pk_pollock_SSB_posteriors.RDS
new file mode 100644
index 0000000..95f124d
Binary files /dev/null and b/content/data_files/pk_pollock_SSB_posteriors.RDS differ
diff --git a/content/data_files/pkfit0.RDS b/content/data_files/pkfit0.RDS
new file mode 100644
index 0000000..fe74f62
Binary files /dev/null and b/content/data_files/pkfit0.RDS differ
diff --git a/content/data_files/pkfitfinal.RDS b/content/data_files/pkfitfinal.RDS
new file mode 100644
index 0000000..85a9c7f
Binary files /dev/null and b/content/data_files/pkfitfinal.RDS differ
diff --git a/content/data_files/pkinput.RDS b/content/data_files/pkinput.RDS
new file mode 100644
index 0000000..f1d4eb0
Binary files /dev/null and b/content/data_files/pkinput.RDS differ
diff --git a/content/data_files/pkinput0.RDS b/content/data_files/pkinput0.RDS
new file mode 100644
index 0000000..ec817c0
Binary files /dev/null and b/content/data_files/pkinput0.RDS differ
diff --git a/content/figures/AFSC_PK_age_comp_fits_1.png b/content/figures/AFSC_PK_age_comp_fits_1.png
new file mode 100644
index 0000000..1b9dbd6
Binary files /dev/null and b/content/figures/AFSC_PK_age_comp_fits_1.png differ
diff --git a/content/figures/AFSC_PK_age_comp_fits_2.png b/content/figures/AFSC_PK_age_comp_fits_2.png
new file mode 100644
index 0000000..286ea06
Binary files /dev/null and b/content/figures/AFSC_PK_age_comp_fits_2.png differ
diff --git a/content/figures/AFSC_PK_index_fits.png b/content/figures/AFSC_PK_index_fits.png
new file mode 100644
index 0000000..a764e46
Binary files /dev/null and b/content/figures/AFSC_PK_index_fits.png differ
diff --git a/content/figures/AFSC_PK_like_profile_R0.png b/content/figures/AFSC_PK_like_profile_R0.png
new file mode 100644
index 0000000..1397d07
Binary files /dev/null and b/content/figures/AFSC_PK_like_profile_R0.png differ
diff --git a/content/figures/AFSC_PK_ssb_posterior.png b/content/figures/AFSC_PK_ssb_posterior.png
new file mode 100644
index 0000000..ee9cb75
Binary files /dev/null and b/content/figures/AFSC_PK_ssb_posterior.png differ
diff --git a/content/figures/AFSC_PK_ts_comparison.png b/content/figures/AFSC_PK_ts_comparison.png
new file mode 100644
index 0000000..60af3e3
Binary files /dev/null and b/content/figures/AFSC_PK_ts_comparison.png differ
diff --git a/content/figures/AFSC_PK_ts_comparison_relerror.png b/content/figures/AFSC_PK_ts_comparison_relerror.png
new file mode 100644
index 0000000..936c40b
Binary files /dev/null and b/content/figures/AFSC_PK_ts_comparison_relerror.png differ
diff --git a/content/pk_prepare_FIMS_inputs.R b/content/pk_prepare_FIMS_inputs.R
new file mode 100644
index 0000000..a958a40
--- /dev/null
+++ b/content/pk_prepare_FIMS_inputs.R
@@ -0,0 +1,396 @@
+
+## build a FIMS and PK data set that match
+##  need to fill missing years with -999 so it's ignored in FIMS
+ind2 <- 0*pkfit0$rep$Eindxsurv2-999
+ind2[which(years %in% fimsdat$srvyrs2)] <- fimsdat$indxsurv2
+CV2 <- rep(1, length=nyears) # actually SE in log space
+CV2[which(years %in% fimsdat$srvyrs2)] <- fimsdat$indxsurv_log_sd2
+paa2 <- pkfit0$rep$Esrvp2*0-999
+paa2[which(years %in% fimsdat$srv_acyrs2),] <- fimsdat$srvp2
+Npaa2 <- rep(1,nyears)
+Npaa2[which(years %in% fimsdat$srv_acyrs2)] <- fimsdat$multN_srv2
+
+ind3 <- 0*pkfit0$rep$Eindxsurv3-999
+ind3[which(years %in% fimsdat$srvyrs3)] <- fimsdat$indxsurv3
+CV3 <- rep(1, length=nyears) # actually SE in log space
+CV3[which(years %in% fimsdat$srvyrs3)] <- fimsdat$indxsurv_log_sd3
+paa3 <- pkfit0$rep$Esrvp3*0-999
+paa3[which(years %in% fimsdat$srv_acyrs3),] <- fimsdat$srvp3
+Npaa3 <- rep(1,nyears)
+Npaa3[which(years %in% fimsdat$srv_acyrs3)] <- fimsdat$multN_srv3
+
+ind6 <- 0*pkfit0$rep$Eindxsurv6-999
+ind6[which(years %in% fimsdat$srvyrs6)] <- fimsdat$indxsurv6
+CV6 <- rep(1, length=nyears) # actually SE in log space
+CV6[which(years %in% fimsdat$srvyrs6)] <- fimsdat$indxsurv_log_sd6
+paa6 <- pkfit0$rep$Esrvp6*0-999
+paa6[which(years %in% fimsdat$srv_acyrs6),] <- fimsdat$srvp6
+Npaa6 <- rep(1,nyears)
+Npaa6[which(years %in% fimsdat$srv_acyrs6)] <- fimsdat$multN_srv6
+
+## repeat with fish catch at age, using expected in missing years
+caa <- pkfit0$rep$Ecatp*0-999
+caa[which(years %in% fimsdat$fshyrs),] <- fimsdat$catp
+Ncaa <- rep(1,nyears)
+Ncaa[which(years %in% fimsdat$fshyrs)] <- fimsdat$multN_fsh
+
+
+
+## put into fims friendly form
+res <- data.frame(type = character(),
+                  name = character(),
+                  age = integer(),
+                  datestart = character(),
+                  dateend = character(),
+                  value = double(),
+                  unit = character(),
+                  uncertainty = double())
+landings <- data.frame(type = "landings",
+                       name = "fleet1",
+                       age = NA,
+                       datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
+                       dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
+                       value = as.numeric(fimsdat$cattot)*1e3,
+                       unit = "mt",
+                       uncertainty = fimsdat$cattot_log_sd[1])
+index2 <- data.frame(type = "index",
+                    name = "survey2",
+                    age = NA,
+                    datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
+                    dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
+                    value = ifelse(ind2>0, ind2*1e9, ind2),
+                    unit = "",
+                    uncertainty = CV2)
+index3 <- data.frame(type = "index",
+                    name = "survey3",
+                    age = NA,
+                    datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
+                    dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
+                    value = ifelse(ind3>0, ind3*1e9, ind3),
+                    unit = "",
+                    uncertainty = CV3)
+index6 <- data.frame(type = "index",
+                    name = "survey6",
+                    age = NA,
+                    datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
+                    dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
+                    value = ifelse(ind6>0, ind6*1e9, ind6),
+                    unit = "",
+                    uncertainty = CV6)
+## these have -999 for missing data years
+catchage <- data.frame(type = "age",
+                       name = "fleet1",
+                       age = rep(seq(1,nages), nyears),
+                       datestart = rep(paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"), each=nages),
+                       dateend = rep(paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"), each=nages),
+                       value = as.numeric(t(caa)),
+                       unit = "",
+                       uncertainty = rep(Ncaa, each=nages))
+indexage2 <- data.frame(type = "age",
+                     name = "survey2",
+                     age = rep(seq(1,nages), nyears),
+                     datestart = rep(paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"), each=nages),
+                     dateend = rep(paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"), each=nages),
+                     value = as.numeric(t(paa2)),
+                     unit = "",
+                     uncertainty = rep(Npaa2, each=nages))
+indexage3 <- data.frame(type = "age",
+                     name = "survey3",
+                     age = rep(seq(1,nages), nyears),
+                     datestart = rep(paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"), each=nages),
+                     dateend = rep(paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"), each=nages),
+                     value = as.numeric(t(paa3)),
+                     unit = "",
+                     uncertainty = rep(Npaa3, each=nages))
+indexage6 <- data.frame(type = "age",
+                     name = "survey6",
+                     age = rep(seq(1,nages), nyears),
+                     datestart = rep(paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"), each=nages),
+                     dateend = rep(paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"), each=nages),
+                     value = as.numeric(t(paa6)),
+                     unit = "",
+                     uncertainty = rep(Npaa6, each=nages))
+indexage <- rbind(indexage2, indexage3, indexage6)
+index <- rbind(index2,index3, index6)
+## indexage=indexage2
+## index=index2
+res <- rbind(res, landings, index, catchage, indexage)
+## rm(landings, index, catchage, indexage)
+
+
+age_frame <- FIMS::FIMSFrameAge(res)
+fishery_catch <- FIMS::m_landings(age_frame)
+fishery_agecomp <- FIMS::m_agecomp(age_frame, "fleet1")
+survey_index2 <- FIMS::m_index(age_frame, "survey2")
+survey_agecomp2 <- FIMS::m_agecomp(age_frame, "survey2")
+survey_index3 <- FIMS::m_index(age_frame, "survey3")
+survey_agecomp3 <- FIMS::m_agecomp(age_frame, "survey3")
+survey_index6 <- FIMS::m_index(age_frame, "survey6")
+survey_agecomp6 <- FIMS::m_agecomp(age_frame, "survey6")
+# need to think about how to deal with multiple fleets - only using 1 fleeet for now
+fish_index <- methods::new(Index, nyears)
+fish_age_comp <- methods::new(AgeComp, nyears, nages)
+fish_index$index_data <- fishery_catch
+fish_age_comp$age_comp_data <- fishery_agecomp * catchage$uncertainty#rep(Ncaa, each=nages)
+
+
+### set up fishery
+## fleet selectivity: converted from time-varying ascending
+## slope/intercept to constant double-logistic
+## methods::show(DoubleLogisticSelectivity)
+fish_selex <- methods::new(DoubleLogisticSelectivity)
+fish_selex$inflection_point_asc$value <- parfinal$inf1_fsh_mean
+fish_selex$inflection_point_asc$is_random_effect <- FALSE
+fish_selex$inflection_point_asc$estimated <- estimate_fish_selex
+fish_selex$inflection_point_desc$value <- parfinal$inf2_fsh_mean
+fish_selex$inflection_point_desc$is_random_effect <- FALSE
+fish_selex$inflection_point_desc$estimated <- estimate_fish_selex
+fish_selex$slope_asc$value <- exp(parfinal$log_slp1_fsh_mean)
+fish_selex$slope_asc$is_random_effect <- FALSE
+fish_selex$slope_asc$estimated <- estimate_fish_selex
+fish_selex$slope_desc$value <- exp(parfinal$log_slp2_fsh_mean)
+fish_selex$slope_desc$is_random_effect <- FALSE
+fish_selex$slope_desc$estimated <- estimate_fish_selex
+## create fleet object
+fish_fleet <- methods::new(Fleet)
+fish_fleet$nages <- nages
+fish_fleet$nyears <- nyears
+fish_fleet$log_Fmort <- log(pkfitfinal$rep$F)
+fish_fleet$estimate_F <- estimate_F
+fish_fleet$random_F <- FALSE
+fish_fleet$log_q <- 0 # why is this length two in Chris' case study?
+fish_fleet$estimate_q <- FALSE
+fish_fleet$random_q <- FALSE
+fish_fleet$log_obs_error <- log(landings$uncertainty)
+## fish_fleet$log_obs_error$estimated <- FALSE
+# Next two lines not currently used by FIMS
+fish_fleet$SetAgeCompLikelihood(1)
+fish_fleet$SetIndexLikelihood(1)
+# Set Index, AgeComp, and Selectivity using the IDs from the modules defined above
+fish_fleet$SetObservedIndexData(fish_index$get_id())
+fish_fleet$SetObservedAgeCompData(fish_age_comp$get_id())
+fish_fleet$SetSelectivity(fish_selex$get_id())
+
+
+## Setup survey 2
+survey_fleet_index <- methods::new(Index, nyears)
+survey_age_comp <- methods::new(AgeComp, nyears, nages)
+survey_fleet_index$index_data <- survey_index2
+survey_age_comp$age_comp_data <- survey_agecomp2 * indexage2$uncertainty
+## survey selectivity: ascending logistic
+## methods::show(DoubleLogisticSelectivity)
+survey_selex <- methods::new(DoubleLogisticSelectivity)
+survey_selex$inflection_point_asc$value <- parfinal$inf1_srv2
+survey_selex$inflection_point_asc$is_random_effect <- FALSE
+survey_selex$inflection_point_asc$estimated <- estimate_survey_selex
+survey_selex$slope_asc$value <- exp(parfinal$log_slp1_srv2)
+survey_selex$slope_asc$is_random_effect <- FALSE
+survey_selex$slope_asc$estimated <- estimate_survey_selex
+## not estimated to make it ascending only, fix at input values
+survey_selex$inflection_point_desc$value <- parfinal$inf2_srv2
+survey_selex$inflection_point_desc$is_random_effect <- FALSE
+survey_selex$inflection_point_desc$estimated <- FALSE
+survey_selex$slope_desc$value <- exp(parfinal$log_slp2_srv2)
+survey_selex$slope_desc$is_random_effect <- FALSE
+survey_selex$slope_desc$estimated <- FALSE
+survey_fleet <- methods::new(Fleet)
+survey_fleet$is_survey <- TRUE
+survey_fleet$nages <- nages
+survey_fleet$nyears <- nyears
+survey_fleet$estimate_F <- FALSE
+survey_fleet$random_F <- FALSE
+survey_fleet$log_q <- parfinal$log_q2_mean
+survey_fleet$estimate_q <- estimate_q2
+survey_fleet$random_q <- FALSE
+# sd = sqrt(log(cv^2 + 1)), sd is log transformed
+survey_fleet$log_obs_error <- log(index2$uncertainty)
+## survey_fleet$log_obs_error$estimated <- FALSE
+survey_fleet$SetAgeCompLikelihood(1)
+survey_fleet$SetIndexLikelihood(1)
+survey_fleet$SetSelectivity(survey_selex$get_id())
+survey_fleet$SetObservedIndexData(survey_fleet_index$get_id())
+survey_fleet$SetObservedAgeCompData(survey_age_comp$get_id())
+
+## Setup survey 3
+survey_fleet_index <- methods::new(Index, nyears)
+survey_age_comp <- methods::new(AgeComp, nyears, nages)
+survey_fleet_index$index_data <- survey_index3
+survey_age_comp$age_comp_data <- survey_agecomp3 * indexage3$uncertainty
+## survey selectivity: ascending logistic
+## methods::show(LogisticSelectivity)
+survey_selex <- methods::new(LogisticSelectivity)
+survey_selex$inflection_point$value <- parfinal$inf1_srv3
+survey_selex$inflection_point$is_random_effect <- FALSE
+survey_selex$inflection_point$estimated <- estimate_survey_selex
+survey_selex$slope$value <- exp(parfinal$log_slp1_srv3)
+survey_selex$slope$is_random_effect <- FALSE
+survey_selex$slope$estimated <- estimate_survey_selex
+survey_fleet <- methods::new(Fleet)
+survey_fleet$is_survey <- TRUE
+survey_fleet$nages <- nages
+survey_fleet$nyears <- nyears
+survey_fleet$estimate_F <- FALSE
+survey_fleet$random_F <- FALSE
+survey_fleet$log_q <- parfinal$log_q3_mean
+survey_fleet$estimate_q <- estimate_q3
+survey_fleet$random_q <- FALSE
+# sd = sqrt(log(cv^2 + 1)), sd is log transformed
+survey_fleet$log_obs_error <- log(index3$uncertainty)
+## survey_fleet$log_obs_error$estimated <- FALSE
+survey_fleet$SetAgeCompLikelihood(2)
+survey_fleet$SetIndexLikelihood(2)
+survey_fleet$SetSelectivity(survey_selex$get_id())
+survey_fleet$SetObservedIndexData(survey_fleet_index$get_id())
+survey_fleet$SetObservedAgeCompData(survey_age_comp$get_id())
+
+## Setup survey 6
+survey_fleet_index <- methods::new(Index, nyears)
+survey_age_comp <- methods::new(AgeComp, nyears, nages)
+survey_fleet_index$index_data <- survey_index6
+survey_age_comp$age_comp_data <- survey_agecomp6 * indexage6$uncertainty
+## survey selectivity: ascending logistic
+## methods::show(DoubleLogisticSelectivity)
+survey_selex <- methods::new(DoubleLogisticSelectivity)
+survey_selex$inflection_point_asc$value <- parfinal$inf1_srv6
+survey_selex$inflection_point_asc$is_random_effect <- FALSE
+survey_selex$inflection_point_asc$estimated <- FALSE
+survey_selex$slope_asc$value <- exp(parfinal$log_slp1_srv6)
+survey_selex$slope_asc$is_random_effect <- FALSE
+survey_selex$slope_asc$estimated <- FALSE
+## not estimated to make it ascending only, fix at input values
+survey_selex$inflection_point_desc$value <- parfinal$inf2_srv6
+survey_selex$inflection_point_desc$is_random_effect <- FALSE
+survey_selex$inflection_point_desc$estimated <- estimate_survey_selex
+survey_selex$slope_desc$value <- exp(parfinal$log_slp2_srv6)
+survey_selex$slope_desc$is_random_effect <- FALSE
+survey_selex$slope_desc$estimated <- estimate_survey_selex
+survey_fleet <- methods::new(Fleet)
+survey_fleet$is_survey <- TRUE
+survey_fleet$nages <- nages
+survey_fleet$nyears <- nyears
+survey_fleet$estimate_F <- FALSE
+survey_fleet$random_F <- FALSE
+survey_fleet$log_q <- parfinal$log_q6
+survey_fleet$estimate_q <- estimate_q6
+survey_fleet$random_q <- FALSE
+# sd = sqrt(log(cv^2 + 1)), sd is log transformed
+survey_fleet$log_obs_error <- log(index6$uncertainty)
+## survey_fleet$log_obs_error$estimated <- FALSE
+survey_fleet$SetAgeCompLikelihood(3)
+survey_fleet$SetIndexLikelihood(3)
+survey_fleet$SetSelectivity(survey_selex$get_id())
+survey_fleet$SetObservedIndexData(survey_fleet_index$get_id())
+survey_fleet$SetObservedAgeCompData(survey_age_comp$get_id())
+
+
+
+# Population module
+# recruitment
+recruitment <- methods::new(BevertonHoltRecruitment)
+## methods::show(BevertonHoltRecruitment)
+recruitment$log_sigma_recruit$value <- log(parfinal$sigmaR)
+recruitment$log_rzero$value <- parfinal$mean_log_recruit + log(1e9)
+recruitment$log_rzero$is_random_effect <- FALSE
+recruitment$log_rzero$estimated <- TRUE
+## note: do not set steepness exactly equal to 1, use 0.99 instead in ASAP run
+recruitment$logit_steep$value <- -log(1.0 - .99999) + log(.99999 - 0.2)
+recruitment$logit_steep$is_random_effect <- FALSE
+recruitment$logit_steep$estimated <- FALSE
+recruitment$estimate_log_devs <- estimate_recdevs
+recruitment$log_devs <-  parfinal$dev_log_recruit[-1]
+
+## growth  -- assumes single WAA vector for everything, based on
+## Srv1 above
+waa <- pkinput$dat$wt_srv1
+ewaa_growth <- methods::new(EWAAgrowth)
+ewaa_growth$ages <- ages
+# NOTE: FIMS currently cannot use matrix of WAA, so have to ensure constant WAA over time in ASAP file for now
+ewaa_growth$weights <- waa[1,]
+## NOTE: FIMS assumes SSB calculated at the start of the year, so
+## need to adjust ASAP to do so as well for now, need to make
+## timing of SSB calculation part of FIMS later
+## maturity
+## NOTE: for now tricking FIMS into thinking age 0 is age 1, so need to adjust A50 for maturity because FIMS calculations use ages 0-5+ instead of 1-6
+maturity <- new(LogisticMaturity)
+maturity$inflection_point$value <- 4.5
+maturity$inflection_point$is_random_effect <- FALSE
+maturity$inflection_point$estimated <- FALSE
+maturity$slope$value <- 1.5
+maturity$slope$is_random_effect <- FALSE
+maturity$slope$estimated <- FALSE
+
+# population
+population <- new(Population)
+population$log_M <- log(as.numeric(t(matrix(rep(pkfitfinal$rep$M, each=nyears), nrow=nyears))))
+population$estimate_M <- FALSE
+population$log_init_naa <- c(log(pkfitfinal$rep$recruit[1]), log(pkfitfinal$rep$initN)) + log(1e9)
+population$estimate_init_naa <- FALSE # TRUE , NOTE: fixing at ASAP estimates to test SSB calculations
+population$nages <- nages
+population$ages <- ages
+population$nfleets <- 2 # 1 fleet and 1 survey
+population$nseasons <- nseasons
+population$nyears <- nyears
+## population$prop_female <- 1.0 # ASAP assumption
+population$SetMaturity(maturity$get_id())
+population$SetGrowth(ewaa_growth$get_id())
+population$SetRecruitment(recruitment$get_id())
+
+## Quick function to compare models
+get_long_outputs<- function(fims, tmb){
+  ## estimaed catch is pretty close I think
+  naa <- matrix(fims$naa[[1]], ncol=nages, byrow=TRUE)/1e9
+  catch <- data.frame(year=years, name='catch', FIMS=fims$exp_catch[[1]]/1e3, TMB=tmb$Ecattot)
+  ## spawnbio
+  ssb <- data.frame(year=years, name='SSB',
+                    FIMS=fims$ssb[[1]][-55]/1e9,
+                    TMB=tmb$Espawnbio)
+  bio <- data.frame(year=years, name='biomass', FIMS=fims$biomass[[1]][-55]/1e9, TMB=tmb$Etotalbio)
+  fmort <- data.frame(year=years, name='F', FIMS=fims$F_mort[[1]], TMB=tmb$F)
+  recruit <- data.frame(year=years, name='recruit', FIMS=naa[-55,1], TMB=tmb$recruit)
+  ## expected index survey 2
+  eind2 <- data.frame(year=years, name='Index2', FIMS=fims$exp_index[[2]]/1e9, TMB=tmb$Eindxsurv2)
+  eind3 <- data.frame(year=years, name='Index3', FIMS=fims$exp_index[[3]]/1e9, TMB=tmb$Eindxsurv3)
+  eind6 <- data.frame(year=years, name='Index6', FIMS=fims$exp_index[[4]]/1e9, TMB=tmb$Eindxsurv6)
+  xx <- rbind(catch,ssb,bio,fmort,eind2,eind3, eind6, recruit) %>%
+    pivot_longer(cols=-c(name, year), names_to='platform') %>%
+    group_by(year, name) %>%
+    mutate(relerror=(value-value[platform=='TMB'])/value[platform=='TMB']) %>%
+    ungroup()
+  return(xx)
+}
+
+get_acomp_fits <- function(tmb,fims1, fims2, fleet, years){
+  ind <- which(tmb$years %in% years)
+  if(fleet==1){
+    y <- tmb$res_fish[,11:20]
+    obs <- tmb$res_fish[,1:10]
+  } else if(fleet==2) {
+    y <- tmb$res_srv2[,11:20]
+    obs <- tmb$res_srv2[,1:10]
+  } else if(fleet==3) {
+    y <- tmb$res_srv3[,11:20]
+    obs <- tmb$res_srv3[,1:10]
+  } else if(fleet==4) {
+    y <- tmb$res_srv6[,11:20]
+    obs <- tmb$res_srv6[,1:10]
+  } else { stop("bad fleet") }
+
+  lab <- c('Fishery', 'Survey 2', 'Survey 3', 'Survey 6')[fleet]
+  x1 <- matrix(fims1$cnaa[[fleet]], ncol=10, byrow=TRUE)[ind,]
+  x1 <- x1/rowSums(x1)
+  x2 <- matrix(fims2$cnaa[[fleet]], ncol=10, byrow=TRUE)[ind,]
+  x2 <- x2/rowSums(x2)
+  dimnames(y) <-  dimnames(x1) <- dimnames(x2) <-
+    list(year=years, age=1:10)
+  dimnames(obs) <- list(year=years, age=1:10)
+  x1 <- reshape2::melt(x1,value.name='paa') %>%
+    cbind(platform='FIMS init', type=lab)
+  x2 <- reshape2::melt(x2,value.name='paa') %>%
+    cbind(platform='FIMS opt', type=lab)
+  y <- reshape2::melt(y,value.name='paa') %>%
+    cbind(platform='TMB', type=lab)
+  obs <- reshape2::melt(obs,value.name='paa') %>%
+    cbind(platform='obs', type=lab)
+  out <- rbind(x1,x2,y,obs)
+  out
+}