R/wham_plots_tables.R

Defines functions plot.annual.SPR.targets plot.SPR.table get_YPR get_SPR plot.maturity plot.waa plot.index.age.comp.bubbles plot.index.input plot.catch.age.comp.bubbles plot.catch.by.fleet plot.M plot.cv plot.fleet.F plot.SARC.R.SSB plot.recr.ssb.yr plot.recruitment.devs plot.NAA plot.SSB.AA plot.SSB.F.trend plot.q.trend plot.index.sel.blocks plot.fleet.sel.blocks plot.index.age.comp.resids plot.catch.age.comp.resids multinomial.pearson.fn plot.index.age.comp plot.catch.age.comp plot.NAA.res plot.NAA.4.panel plot.index.4.panel plot.catch.4.panel plot.fleet.stdresids.fn plot.index.stdresids.fn plot.ecov.stdresids.fn plot.all.stdresids.fn get.wham.results.fn RMSE.table.fn get.RMSEs.fn hi.cor.fn residual.bubbleplot.fn plot.ll.table.fn fit.summary.text.plot.fn mypalette plot.osa.residuals plot.ecov

plot.ecov <- function(mod, plot.pad = FALSE, do.tex=FALSE, do.png=FALSE, fontfam="", res=72, od){
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  years <- seq(from=dat$year1_Ecov, by=1, length.out=dat$n_years_Ecov)
  years_full <- seq(from=dat$year1_Ecov, by=1, length.out=dat$n_years_Ecov+dat$n_years_proj_Ecov)

  ecov.pred = mod$rep$Ecov_x
  ecov.obs = dat$Ecov_obs[1:dat$n_years_Ecov,,drop=F]
  # ecov.obs.sig = mod$rep$Ecov_obs_sigma # Ecov_obs_sigma now a derived quantity in sdrep
  if(class(mod$sdrep)[1] == "sdreport"){
    sdrep = summary(mod$sdrep)
  } else {
    sdrep = mod$sdrep
  }

  ecov.obs.sig = mod$rep$Ecov_obs_sigma # Ecov_obs_sigma is filled with fixed, or estimated values (fe or re) for each covariate depending on the respective options
  # if("Ecov_obs_logsigma" %in% names(mod$env$par)){
  #   ecov.obs.sig = matrix(exp(sdrep[rownames(sdrep) %in% "Ecov_obs_logsigma",1]), ncol=dat$n_Ecov) # all the same bc obs_sig_var --> 0
  #   if(dim(ecov.obs.sig)[1] == 1) ecov.obs.sig = matrix(rep(ecov.obs.sig, dim(ecov.obs)[1]), ncol=dat$n_Ecov, byrow=T)
  # } else {
  #   ecov.obs.sig = exp(mod$input$par$Ecov_obs_logsigma)
  # }
  ecov.use = dat$Ecov_use_obs[1:dat$n_years_Ecov,,drop=F]
  ecov.obs.sig = ecov.obs.sig[1:dat$n_years_Ecov,,drop=F]
  ecov.obs.sig[ecov.use == 0] <- NA
  ecov.pred.se = matrix(sdrep[rownames(sdrep) %in% "Ecov_x",2], ncol=dat$n_Ecov)

  # default: don't plot the padded entries that weren't used in ecov likelihood
  if(!plot.pad) ecov.obs[ecov.use == 0] <- NA

  ecov.res = (ecov.obs - ecov.pred[1:dat$n_years_Ecov,]) / ecov.obs.sig # standard residual (obs - pred)

  ecovs <- 1:dat$n_Ecov
  plot.colors = mypalette(dat$n_Ecov)
  for (i in ecovs)
  {
    if(do.tex) cairo_pdf(file.path(od, paste0("Ecov_",i,".pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("Ecov_",i,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)

    ecov.pred.low <- ecov.pred[,i] - 1.96 * ecov.pred.se[,i]
    ecov.pred.high <- ecov.pred[,i] + 1.96 * ecov.pred.se[,i]
    ecov.low <- ecov.obs[,i] - 1.96 * ecov.obs.sig[,i]
    ecov.high <- ecov.obs[,i] + 1.96 * ecov.obs.sig[,i]
    y.min <- ifelse(min(ecov.low,na.rm=T) < 0, 1.1*min(ecov.low,na.rm=T), 0.9*min(ecov.low,na.rm=T))
    y.max <- ifelse(max(ecov.high,na.rm=T) < 0, 0.9*max(ecov.high,na.rm=T), 1.1*max(ecov.high,na.rm=T))
    if(max(ecov.pred[,i],na.rm=T) > y.max) y.max <- max(ecov.pred[,i],na.rm=T)
    if(min(ecov.pred[,i],na.rm=T) < y.min) y.min <- min(ecov.pred[,i],na.rm=T)
    plot(years_full, ecov.pred[,i], type='n', xlab="Year", ylab=unlist(dat$Ecov_label)[i],
         ylim=c(y.min, y.max))
    polygon(c(years_full,rev(years_full)), c(ecov.pred.low, rev(ecov.pred.high)), col=adjustcolor(plot.colors[i], alpha.f=0.4), border = "transparent")
    arrows(years, ecov.low, years, ecov.high, length=0)
    points(years, ecov.obs[,i], pch=19)
    lines(years_full, ecov.pred[,i], col=plot.colors[i], lwd=3)
    title (paste0("Ecov ",i, ": ",unlist(dat$Ecov_label)[i]), outer=T, line=-1)
    if(dat$n_years_proj_Ecov > 0) abline(v=tail(years,1), lty=2)

    if(do.tex | do.png) dev.off() else par(origpar)
  }
  #plot standard deviations if estimated as random effects
  ecov.obs.sigma.cv = as.list(mod$sdrep, "Std")$Ecov_obs_logsigma_re
  for (i in ecovs) if(dat$Ecov_obs_sigma_opt[i]==4)
  {
    if(do.tex) cairo_pdf(file.path(od, paste0("Ecov_",i,"_sig_re.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("Ecov_",i,'_sig_re.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    sigmas = ecov.obs.sig[,i]
    sigmas.cv = ecov.obs.sigma.cv[,i]
    sigmas.low <- exp(log(sigmas) - 1.96 * sigmas.cv)
    sigmas.high <- exp(log(sigmas) + 1.96 * sigmas.cv)
    y.min = min(sigmas.low, na.rm =T)
    y.max = max(sigmas.high, na.rm =T)
    plot(years, sigmas, type='n', xlab="Year", ylab=unlist(dat$Ecov_label)[i],
         ylim=c(y.min, y.max))
    polygon(c(years,rev(years)), c(sigmas.low, rev(sigmas.high)), col=adjustcolor(plot.colors[i], alpha.f=0.4), border = "transparent")
    points(years, sigmas, pch=19)
    lines(years, sigmas, col=plot.colors[i], lwd=3)
    title (paste0("Ecov ",i, ": ",unlist(dat$Ecov_label)[i], " SD"), outer=T, line=-1)
    if(dat$n_years_proj_Ecov > 0) abline(v=tail(years,1), lty=2)

    if(do.tex | do.png) dev.off() else par(origpar)
  }
}

plot.osa.residuals <- function(mod, do.tex=FALSE, do.png=FALSE, fontfam="", res=72, od){
  origpar <- par(no.readonly = TRUE)
  years <- mod$years
  if("logcatch" %in% mod$osa$type){
    dat <- subset(mod$osa, type == "logcatch")
    dat$fleet <- factor(as.character(dat$fleet))
    n.fleets <- length(table(dat$fleet))
    plot.colors = mypalette(n.fleets)
    for(f in 1:n.fleets){
      tmp <- subset(dat, fleet==names(table(dat$fleet))[f])
      if(do.tex) cairo_pdf(file.path(od, paste0("OSA_resid_catch_4panel_fleet_",f,".pdf")), family = fontfam, height = 10, width = 10)
      if(do.png) png(filename = file.path(od, paste0("OSA_resid_catch_4panel_fleet_",f,'.png')), width = 10*res, height = 10*res, res = res, pointsize = 12, family = fontfam)
      par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(2,2))

      # set plot lims using max residual for any component (easier to compare if all the same)
      ylim.max <- max(abs(range(dat$residual, na.rm=TRUE)))
      if(is.infinite(ylim.max)) {
        cat("Infinite osa residuals for aggregate catch in fleet ", f, ", so using +/-10 for range of y axis \n")
        ylim.max = 10
      }
      ylims <- c(-ylim.max, ylim.max)

      # 1. trend vs. year
      plot(years, tmp$residual, type='p', col=plot.colors[f], pch=19, xlab="Year", ylab="OSA Residuals",
           ylim=ylims)
      abline(h=0, col=plot.colors[f], lwd=2)

      # 2. trend vs. fitted val
      plot(mod$rep$pred_log_catch[1:length(mod$years),f], tmp$residual, type='p', col=plot.colors[f], pch=19, xlab="Log(Predicted Catch)", ylab="OSA Residuals",
           ylim=ylims)
      abline(h=0, col=plot.colors[f], lwd=2)

      # 3. histogram
      xfit<-seq(-ylim.max, ylim.max, length=100)
      yfit<-dnorm(xfit)
      hist(tmp$residual, ylim=c(0,1.05*max(yfit)), xlim=ylims, plot=T, xlab="OSA Residuals", ylab="Probability Density", col=plot.colors[f], freq=F, main=NULL, breaks="scott")
      lines(xfit, yfit)

      # 4. QQ plot modified from car:::qqPlot.default
      notNA <- tmp$residual[!is.na(tmp$residual)]
      ord.x <- notNA[order(notNA)]
      n <- length(ord.x)
      P <- ppoints(n)
      z <- qnorm(P, mean=0, sd=1)
      plot(z, ord.x, xlab="Std Normal Quantiles", ylab="OSA Residual Quantiles", main="", type = "n")
      grid(lty = 1, equilogs = FALSE)
      box()
      points(z, ord.x, col=plot.colors[f], pch=19)
      abline(0,1, col=plot.colors[f])
      conf = 0.95
      zz <- qnorm(1 - (1 - conf)/2)
      SE <- (1/dnorm(z)) * sqrt(P * (1 - P)/n)
      upper <- z + zz * SE
      lower <- z - zz * SE
      lines(z, upper, lty=2, col=plot.colors[f])
      lines(z, lower, lty=2, col=plot.colors[f])

      title (paste0("OSA residual diagnostics: Fleet ",f), outer=T, line=-1)
      if(do.tex | do.png) dev.off() else par(origpar)
    }
  }
  if("catchpaa" %in% mod$osa$type) {
    dat <- subset(mod$osa, type == "catchpaa")
    n.fleets <- mod$input$data$n_fleets
    plot.colors = mypalette(n.fleets)
    for(f in 1:n.fleets) if(any(mod$input$data$use_catch_paa[,f]==1)){
      if(do.tex) cairo_pdf(file.path(od, paste0("OSA_resid_paa_6panel_fleet_",f,".pdf")), family = fontfam, height = 10, width = 10)
      if(do.png) png(filename = file.path(od, paste0("OSA_resid_paa_6panel_fleet_",f,'.png')), width = 10*res, height = 10*res, res = res, pointsize = 12, family = fontfam)
      
      par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(3,2))
      yind = which(mod$input$data$use_catch_paa[,f] ==1)
      resids = vals = ages = pyears = cohorts = matrix(NA, nrow = mod$input$data$n_years_model, 
        ncol = mod$input$data$n_ages)
      for(j in yind){
        tmp = subset(dat, year == j & fleet == paste0("fleet_",f))
        resids[j,tmp$age] = tmp$residual
        vals[j,tmp$age] = tmp$val
        ages[j,tmp$age] = tmp$age
        pyears[j,tmp$age] = mod$years[tmp$year]
        cohorts[j,tmp$age] = tmp$cohort
      }

      # set plot lims using max residual for any component (easier to compare if all the same)
      ylim.max <- max(abs(range(resids, na.rm=TRUE)))
      if(is.infinite(ylim.max)) {
        cat("Infinite osa residuals for catch proportions at age in fleet ", f, ", so using +/-10 for range \n")
        ylim.max = 10
      }
      ylims <- c(-ylim.max, ylim.max)

      #plot_years = rep(years, NCOL(resids))#years[as.integer(tmp$year)]
      # 1. trend vs. year
      plot(as.integer(pyears), resids, type='p', col=plot.colors[f], pch=19, xlab="Year", ylab="OSA Residuals",
           ylim=ylims)
      grid(lty = 2, col = "grey")
      abline(h=0, col=plot.colors[f], lwd=2)
      abline(h=-qnorm(0.975), col=plot.colors[f], lwd=2, lty = 2)
      abline(h=qnorm(0.975), col=plot.colors[f], lwd=2, lty = 2)
      # 2. trend vs. age
      plot(as.integer(ages), resids, type='p', col=plot.colors[f], pch=19, xlab="Age", ylab="OSA Residuals",
           ylim=ylims, xaxt ="n")
      axis(1, at = 1:mod$input$data$n_ages, labels = mod$ages.lab)
      grid(lty = 2, col = "grey")
      abline(h=0, col=plot.colors[f], lwd=2)
      abline(h=-qnorm(0.975), col=plot.colors[f], lwd=2, lty = 2)
      abline(h=qnorm(0.975), col=plot.colors[f], lwd=2, lty = 2)
      # 3. trend vs. cohort
      plot(min(years) + as.integer(cohorts), resids, type='p', col=plot.colors[f], pch=19, xlab="Cohort", ylab="OSA Residuals",
           ylim=ylims)
      grid(lty = 2, col = "grey")
      abline(h=0, col=plot.colors[f], lwd=2)
      abline(h=-qnorm(0.975), col=plot.colors[f], lwd=2, lty = 2)
      abline(h=qnorm(0.975), col=plot.colors[f], lwd=2, lty = 2)
      # 4. trend vs. observed val
      plot(vals, resids, type='p', col=plot.colors[f], pch=19, xlab="Observed proportion", ylab="OSA Residuals",
           ylim=ylims)
      grid(lty = 2, col = "grey")
      abline(h=0, col=plot.colors[f], lwd=2)
      abline(h=-qnorm(0.975), col=plot.colors[f], lwd=2, lty = 2)
      abline(h=qnorm(0.975), col=plot.colors[f], lwd=2, lty = 2)
      # 5. histogram
      xfit<-seq(-ylim.max, ylim.max, length=100)
      yfit<-dnorm(xfit)
      hist(resids, ylim=c(0,1.05*max(yfit)), xlim=ylims, plot=T, xlab="OSA Residuals", ylab="Probability Density", col=plot.colors[f], freq=F, main=NULL, breaks="scott")
      lines(xfit, yfit)
      # 6. QQ plot modified from car:::qqPlot.default
      notNA <- resids[!is.na(resids)]
      ord.x <- notNA[order(notNA)]
      n <- length(ord.x)
      P <- ppoints(n)
      z <- qnorm(P, mean=0, sd=1)
      plot(z, ord.x, xlab="Std Normal Quantiles", ylab="OSA Residual Quantiles", main="", type = "n")
      grid(lty = 1, equilogs = FALSE)
      box()
      points(z, ord.x, col=plot.colors[f], pch=19)
      abline(0,1, col=plot.colors[f])
      conf = 0.95
      zz <- qnorm(1 - (1 - conf)/2)
      SE <- (1/dnorm(z)) * sqrt(P * (1 - P)/n)
      upper <- z + zz * SE
      lower <- z - zz * SE
      lines(z, upper, lty=2, col=plot.colors[f])
      lines(z, lower, lty=2, col=plot.colors[f])

      title (paste0("age composition OSA residual diagnostics: Fleet ",f), outer=T, line=-1)
      if(do.tex | do.png) dev.off() else par(origpar)
    }
  }

  if("logindex" %in% mod$osa$type){
    dat <- subset(mod$osa, type == "logindex")
    dat$fleet <- factor(as.character(dat$fleet))
    n.fleets <- length(table(dat$fleet))
    plot.colors = mypalette(n.fleets)
    for(f in 1:n.fleets){
      tmp <- subset(dat, fleet==names(table(dat$fleet))[f])
      if(do.tex) cairo_pdf(file.path(od, paste0("OSA_resid_catch_4panel_index_",f,".pdf")), family = fontfam, height = 10, width = 10)
      if(do.png) png(filename = file.path(od, paste0("OSA_resid_catch_4panel_index_",f,'.png')), width = 10*res, height = 10*res, res = res, pointsize = 12, family = fontfam)
      par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(2,2))

      # set plot lims using max residual for any component (easier to compare if all the same)
      ylim.max <- max(abs(range(dat$residual, na.rm=TRUE)))
      if(is.infinite(ylim.max)) {
        cat("Infinite osa residuals for aggregate indices in  index ", f, ", so using +/-10 for range of y axis \n")
        ylim.max = 10
      }
      ylims <- c(-ylim.max, ylim.max)

      # make robust to missing years
      tmp$year <- years[tmp$year]
      tmp <- rbind(tmp[c("year","residual")], data.frame(year=years[!(years %in% tmp$year)], residual=rep(NA, length(years[!(years %in% tmp$year)]))))
      tmp <- tmp[order(tmp$year),]
      # 1. trend vs. year
      plot(tmp$year, tmp$residual, type='p', col=plot.colors[f], pch=19, xlab="Year", ylab="OSA Residuals",
           ylim=ylims)
      abline(h=0, col=plot.colors[f], lwd=2)

      # 2. trend vs. fitted val
      plot(mod$rep$pred_log_indices[1:length(mod$years),f], tmp$residual, type='p', col=plot.colors[f], pch=19, xlab="Log(Predicted Index)", ylab="OSA Residuals",
           ylim=ylims)
      abline(h=0, col=plot.colors[f], lwd=2)

      # 3. histogram
      xfit<-seq(-ylim.max, ylim.max, length=100)
      yfit<-dnorm(xfit)
      hist(tmp$residual, ylim=c(0,1.05*max(yfit)), xlim=ylims, plot=T, xlab="OSA Residuals", ylab="Probability Density", col=plot.colors[f], freq=F, main=NULL, breaks="scott")
      lines(xfit, yfit)

      # 4. QQ plot modified from car:::qqPlot.default
      notNA <- tmp$residual[!is.na(tmp$residual)]
      ord.x <- notNA[order(notNA)]
      n <- length(ord.x)
      P <- ppoints(n)
      z <- qnorm(P, mean=0, sd=1)
      plot(z, ord.x, xlab="Std Normal Quantiles", ylab="OSA Residual Quantiles", main="", type = "n")
      grid(lty = 1, equilogs = FALSE)
      box()
      points(z, ord.x, col=plot.colors[f], pch=19)
      abline(0,1, col=plot.colors[f])
      conf = 0.95
      zz <- qnorm(1 - (1 - conf)/2)
      SE <- (1/dnorm(z)) * sqrt(P * (1 - P)/n)
      upper <- z + zz * SE
      lower <- z - zz * SE
      lines(z, upper, lty=2, col=plot.colors[f])
      lines(z, lower, lty=2, col=plot.colors[f])

      title (paste0("OSA residual diagnostics: Index ",f), outer=T, line=-1)
      if(do.tex | do.png) dev.off() else par(origpar)
    }
  }
  if("indexpaa" %in% mod$osa$type) {
    dat <- subset(mod$osa, type == "indexpaa")
    dat$fleet <- factor(as.character(dat$fleet))
    dat$residual[which(is.infinite(dat$residual))] = NA #some happen for zeros or last age class
    #dat$residual[which(as.integer(dat$age) == mod$input$data$n_ages)] = NA #remove last age class
    
    n.indices <- mod$input$data$n_indices
    plot.colors = mypalette(n.indices)

    for(i in 1:n.indices) if(any(mod$input$data$use_index_paa[,i]==1)){
      if(do.tex) cairo_pdf(file.path(od, paste0("OSA_resid_paa_6panel_index_",i,".pdf")), family = fontfam, height = 10, width = 10)
      if(do.png) png(filename = file.path(od, paste0("OSA_resid_paa_6panel_index_",i,'.png')), width = 10*res, height = 10*res, res = res, pointsize = 12, family = fontfam)

      par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(3,2))
      yind = which(mod$input$data$use_index_paa[,i] ==1)
      resids = vals = ages = pyears = cohorts = matrix(NA, nrow = mod$input$data$n_years_model, 
        ncol = mod$input$data$n_ages)
      for(j in yind){
        tmp = subset(dat, year == j & fleet == paste0("index_",i))
        resids[j,tmp$age] = tmp$residual
        vals[j,tmp$age] = tmp$val
        ages[j,tmp$age] = tmp$age
        pyears[j,tmp$age] = mod$years[tmp$year]
        cohorts[j,tmp$age] = tmp$cohort
      }

      # set plot lims using max residual for any component (easier to compare if all the same)
      ylim.max <- max(abs(range(resids, na.rm=TRUE)))
      if(is.infinite(ylim.max)) {
        cat("Infinite osa residuals for proportions at age in  index ", f, ", so using +/-10 for range \n")
        ylim.max = 10
      }
      ylims <- c(-ylim.max, ylim.max)

      #plot_years = rep(years, NCOL(resids))#years[as.integer(tmp$year)]
      # 1. trend vs. year
      plot(as.integer(pyears), resids, type='p', col=plot.colors[i], pch=19, xlab="Year", ylab="OSA Residuals",
           ylim=ylims)
      grid(lty = 2, col = "grey")
      abline(h=0, col=plot.colors[i], lwd=2)
      abline(h=-qnorm(0.975), col=plot.colors[i], lwd=2, lty = 2)
      abline(h=qnorm(0.975), col=plot.colors[i], lwd=2, lty = 2)
      # 2. trend vs. age
      plot(as.integer(ages), resids, type='p', col=plot.colors[i], pch=19, xlab="Age", ylab="OSA Residuals",
           ylim=ylims, xaxt ="n")
      axis(1, at = 1:mod$input$data$n_ages, labels = mod$ages.lab)
      grid(lty = 2, col = "grey")
      abline(h=0, col=plot.colors[i], lwd=2)
      abline(h=-qnorm(0.975), col=plot.colors[i], lwd=2, lty = 2)
      abline(h=qnorm(0.975), col=plot.colors[i], lwd=2, lty = 2)
      # 3. trend vs. cohort
      plot(min(years) + as.integer(cohorts), resids, type='p', col=plot.colors[i], pch=19, xlab="Cohort", ylab="OSA Residuals",
           ylim=ylims)
      grid(lty = 2, col = "grey")
      abline(h=0, col=plot.colors[i], lwd=2)
      abline(h=-qnorm(0.975), col=plot.colors[i], lwd=2, lty = 2)
      abline(h=qnorm(0.975), col=plot.colors[i], lwd=2, lty = 2)
      # 4. trend vs. observed val
      plot(vals, resids, type='p', col=plot.colors[i], pch=19, xlab="Observed proportion", ylab="OSA Residuals",
           ylim=ylims)
      grid(lty = 2, col = "grey")
      abline(h=0, col=plot.colors[i], lwd=2)
      abline(h=-qnorm(0.975), col=plot.colors[i], lwd=2, lty = 2)
      abline(h=qnorm(0.975), col=plot.colors[i], lwd=2, lty = 2)
      # 5. histogram
      xfit<-seq(-ylim.max, ylim.max, length=100)
      yfit<-dnorm(xfit)
      hist(resids, ylim=c(0,1.05*max(yfit)), xlim=ylims, plot=T, xlab="OSA Residuals", ylab="Probability Density", col=plot.colors[i], freq=F, main=NULL, breaks="scott")
      lines(xfit, yfit)
      # 6. QQ plot modified from car:::qqPlot.default
      notNA <- resids[!is.na(resids)]
      ord.x <- notNA[order(notNA)]
      n <- length(ord.x)
      P <- ppoints(n)
      z <- qnorm(P, mean=0, sd=1)
      plot(z, ord.x, xlab="Std Normal Quantiles", ylab="OSA Residual Quantiles", main="", type = "n")
      grid(lty = 1, equilogs = FALSE)
      box()
      points(z, ord.x, col=plot.colors[i], pch=19)
      abline(0,1, col=plot.colors[i])
      conf = 0.95
      zz <- qnorm(1 - (1 - conf)/2)
      SE <- (1/dnorm(z)) * sqrt(P * (1 - P)/n)
      upper <- z + zz * SE
      lower <- z - zz * SE
      lines(z, upper, lty=2, col=plot.colors[i])
      lines(z, lower, lty=2, col=plot.colors[i])

      title (paste0("age composition OSA residual diagnostics: Index ",i), outer=T, line=-1)
      if(do.tex | do.png) dev.off() else par(origpar)
    }
  }

  if(!all(mod$env$data$Ecov_model == 0)){
    dat <- subset(mod$osa, type == "Ecov")
    dat$fleet <- factor(as.character(dat$fleet))
    n.fleets <- length(table(dat$fleet))
    plot.colors = mypalette(n.fleets)
    for(f in 1:n.fleets){
      tmp <- subset(dat, fleet==names(table(dat$fleet))[f])
      ecov_years = mod$env$data$year1_Ecov + 0:(mod$env$data$n_years_Ecov-1)
      tmp$year <- ecov_years[tmp$year+1] # year in osa is MODEL year, not Ecov year
      # tmp$year <- seq(mod$env$data$year1_Ecov, by=1, length.out=mod$env$data$n_years_Ecov)
#      if(mod$env$data$year1_Ecov < mod$env$data$year1_model){
#        tmp <- rbind(data.frame(year=seq(mod$env$data$year1_Ecov, length.out=mod$env$data$year1_model-mod$env$data$year1_Ecov),
#                                fleet=names(table(dat$fleet))[f],
#                                age=NA, type="Ecov", val=NA, ind=NA, residual=NA), tmp)
#      }
#      end_Ecov <- mod$env$data$year1_Ecov + mod$env$data$n_years_Ecov - 1
#      end_model <- mod$env$data$year1_model + mod$env$data$n_years_model - 1
#      if(end_Ecov > end_model){
#        tmp <- rbind(tmp, data.frame(year=seq(end_model+1, length.out=end_Ecov-end_model),
#                                fleet=names(table(dat$fleet))[f],
#                                age=NA, type="Ecov", val=NA, ind=NA, residual=NA))
#      }
      tmp$pred <- mod$rep$Ecov_x[1:dim(tmp)[1],f]
      tmp <- subset(tmp, !is.nan(tmp$residual))
      if(do.tex) cairo_pdf(file.path(od, paste0("OSA_resid_ecov_4panel_",f,".pdf")), family = fontfam, height = 10, width = 10)
      if(do.png) png(filename = file.path(od, paste0("OSA_resid_ecov_4panel_",f,'.png')), width = 10*res, height = 10*res, res = res, pointsize = 12, family = fontfam)
      par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(2,2))

      # set plot lims using max residual for any component (easier to compare if all the same)
      ylim.max <- max(abs(range(dat$residual, na.rm=TRUE)))
      if(is.infinite(ylim.max)) {
        cat("Infinite osa residuals for Environmental observations in series ", f, ", so using +/-10 for range of y axis \n")
        ylim.max = 10
      }
      ylims <- c(-ylim.max, ylim.max)

      # 1. trend vs. year
      plot(tmp$year, tmp$residual, type='p', col=plot.colors[f], pch=19, xlab="Year", ylab="OSA Residuals",
           ylim=ylims)
      abline(h=0, col=plot.colors[f], lwd=2)

      # 2. trend vs. fitted val
      plot(tmp$pred, tmp$residual, type='p', col=plot.colors[f], pch=19, xlab=paste0("Predicted ", mod$env$data$Ecov_label[[1]][f]), ylab="OSA Residuals",
           ylim=ylims)
      abline(h=0, col=plot.colors[f], lwd=2)

      # 3. histogram
      xfit<-seq(-ylim.max, ylim.max, length=100)
      yfit<-dnorm(xfit)
      hist(tmp$residual, ylim=c(0,1.05*max(yfit)), xlim=ylims, plot=T, xlab="OSA Residuals", ylab="Probability Density", col=plot.colors[f], freq=F, main=NULL, breaks="scott")
      lines(xfit, yfit)

      # 4. QQ plot modified from car:::qqPlot.default
      notNA <- tmp$residual[!is.na(tmp$residual)]
      ord.x <- notNA[order(notNA)]
      n <- length(ord.x)
      P <- ppoints(n)
      z <- qnorm(P, mean=0, sd=1)
      plot(z, ord.x, xlab="Std Normal Quantiles", ylab="OSA Residual Quantiles", main="", type = "n")
      grid(lty = 1, equilogs = FALSE)
      box()
      points(z, ord.x, col=plot.colors[f], pch=19)
      abline(0,1, col=plot.colors[f])
      conf = 0.95
      zz <- qnorm(1 - (1 - conf)/2)
      SE <- (1/dnorm(z)) * sqrt(P * (1 - P)/n)
      upper <- z + zz * SE
      lower <- z - zz * SE
      lines(z, upper, lty=2, col=plot.colors[f])
      lines(z, lower, lty=2, col=plot.colors[f])

      title (paste0("OSA residual diagnostics: Ecov ",f," (",mod$env$data$Ecov_label[[1]][f],")"), outer=T, line=-1)
      if(do.tex | do.png) dev.off() else par(origpar)
    }
  }
}

mypalette = function(n){
  palette.fn <- colorRampPalette(c("dodgerblue","green","red"), space = "Lab")
  palette.fn(n)
}

fit.summary.text.plot.fn <- function(mod){
  acm = c("Multinomial", "Dirichlet-multinomial", "Dirichlet (miss0)", "Dirichlet (pool0)","Logistic normal (miss0)",
    "Logistic normal AR1 corr (miss0)", "Logistic normal (pool0)", "ZI-logistic normal(1)","ZI-logistic normal(2)", "MV Tweedie")
  selmods = c("Age-specific", "Logistic(+)", "Double-Logistic", "Logistic(-)")
  recs <- c("Random walk","Random about mean","Bev-Holt","Ricker")
  env.mod <- c("RW", "AR1")
  # env.where <- c('Recruitment','Growth','Mortality')
  env.where <- c('Recruitment','Mortality',paste0("q for index ", 1:mod$env$data$n_indices))  
  env.how <- c("controlling", "limiting", "lethal", "masking", "directive")
  fleet_selblocks = lapply(1:mod$env$data$n_fleets, function(x) unique(mod$env$data$selblock_pointer_fleets[,x]))
  index_selblocks = lapply(1:mod$env$data$n_indices, function(x) unique(mod$env$data$selblock_pointer_indices[,x]))
  plot(1:10,1:10,type='n',axes=F,xlab="",ylab="")
  nl = 10
  text(5,nl <- nl-0.5,mod$model_name)
  text(5,nl <- nl-0.5,paste0("Model years: ", min(mod$years), "-", max(mod$years)))
  proj.yrs <- tail(mod$years_full, mod$env$data$n_years_proj)
  proj.txt <- ifelse(mod$env$data$do_proj == 0, "none", paste0(min(proj.yrs), "-", max(proj.yrs)))
  text(5,nl <- nl-0.5,paste0("Projection years: ", proj.txt))
  text(5,nl <- nl-0.5,paste0("Number of fleets: ", mod$env$data$n_fleets))
  text(5,nl <- nl-0.5, paste0("Fleet Age Comp Models: ", paste(acm[mod$env$data$age_comp_model_fleets], collapse = ", ")))
  text(5,nl <- nl-0.5,paste0("Number of indices: ", mod$env$data$n_indices))
  text(5,nl <- nl-0.5, paste0("Index Age Comp Models: ", paste(acm[mod$env$data$age_comp_model_indices], collapse = ", ")))
  text(5,nl <- nl-0.5,paste0("Recruitment model: ", recs[mod$env$data$recruit_model]))
  if(!all(mod$env$data$Ecov_model == 0)){
    for(ec in 1:mod$env$data$n_Ecov) {
      ec.where = env.where[which(mod$env$data$Ecov_where[ec,]==1)]
      ec.mod = env.mod[mod$env$data$Ecov_model[ec]]
      ec.label = mod$env$data$Ecov_label[[1]][ec]
      if(length(ec.mod)) {
        out = paste0("Environmental covariate ", ec,": ", ec.label," modeled as ",ec.mod,".")
        if(length(ec.where)) out = paste0(out, paste0(" Effects on ", paste(ec.where,collapse = ','), " estimated."))
        else out = paste(out, " No effects estimated.")
        text(5,nl <- nl-0.5, out)
      }
    }
  } else {
    text(5,nl <- nl-0.5, "No Environmental covariates.")
  }
  text(5,nl <- nl-0.5,paste0("Number of Selectivity blocks: ", mod$env$data$n_selblocks))
  text(5,nl <- nl-0.5, paste0("Selectivity Block Types: ", paste(selmods[mod$env$data$selblock_models], collapse = ", ")))
  for(i in 1:length(fleet_selblocks)) text(5,nl <- nl-0.5, paste0("Fleet ", i, " Selectivity Blocks: ", paste(fleet_selblocks[i], collapse = ", ")))
  for(i in 1:length(index_selblocks)) text(5,nl <- nl-0.5, paste0("Index ", i, " Selectivity Blocks: ", paste(index_selblocks[i], collapse = ", ")))
  text(5,nl <- nl-0.5,paste0("Run date: ",format(mod$date, usetz = TRUE)))
  text(5,nl <- nl-0.5,paste0("Run directory: ",mod$dir))
  text(5,nl <- nl-0.5,paste0("WHAM version: ",mod$wham_version)) 
  text(5,nl <- nl-0.5,paste0("TMB version: ",mod$TMB_version)) 

  if(mod$is_sdrep)
  {
    text(5,nl <- nl-0.5,paste0("sdreport() performed",
      ifelse(mod$na_sdrep, ", but with NAs for some variance estimates.", " with all variance estimates provided.")))
  }
  else
  {
    text(5,nl <- nl-0.5,"Warning: run did not provide pos-def Hessian or sdreport() not performed", col="red")
  }
  mgind = which(abs(mod$final_gradient) == max(abs(mod$final_gradient)))
  text(5,nl <- nl-0.5,paste0("Maximum absolute gradient: ", names(mod$opt$par)[mgind]," ", format(mod$final_gradient[mgind],digits =4)))
  text(5,nl <- nl-0.5,paste0("Number of fixed effects = ",length(mod$opt$par),", Number of random effects = ", length(mod$env$random)))

  return()
}

plot.ll.table.fn <- function(mod,plot.colors){
  par(mfrow=c(1,1) )

  npar <- length(mod$opt$par)
  lls = mod$rep[c(grep("nll",names(mod$rep)), grep("lprior_b", names(mod$rep)))]
  ll.names = names(lls)
  #n.like <- length(lls)
  n_fleets = mod$env$data$n_fleets
  n_indices = mod$env$data$n_indices
  obs.lls = lls[names(lls) %in% c("nll_agg_catch", "nll_catch_acomp", "nll_agg_indices", "nll_index_acomp")]
  obs.lls = lapply(obs.lls, function(x) apply(x,2,sum))
  names(obs.lls$nll_agg_catch) = paste0("Fleet ", 1:n_fleets, " Catch")
  names(obs.lls$nll_catch_acomp) = paste0("Fleet ", 1:n_fleets, " Age Comp")
  names(obs.lls$nll_agg_indices) = paste0("Index ", 1:n_indices, " Catch")
  names(obs.lls$nll_index_acomp) = paste0("Index ", 1:n_indices, " Age Comp")
  names(obs.lls) = NULL
  obs.lls = unlist(obs.lls)
  n.obs.ll = length(obs.lls)
  obs.dists = character(n.obs.ll)
  names(obs.dists) = names(obs.lls)
  obs.dists[grep("Catch", names(obs.lls))] = "log(x) ~ Gaussian"
  acm = c("Multinomial", "Dirichlet-multinomial", "Dirichlet (miss0)", "Dirichlet (pool0)","Logistic normal (miss0)",
    "Logistic normal AR1 corr (miss0)", "Logistic normal (pool0)", "ZI-logistic normal(1)","ZI-logistic normal(2)", "MV Tweedie")
  obs.dists[paste0("Fleet ", 1:n_fleets, " Age Comp")] = paste0("x ~ ", acm[mod$env$data$age_comp_model_fleets])
  obs.dists[paste0("Index ", 1:n_indices, " Age Comp")] = paste0("x ~ ", acm[mod$env$data$age_comp_model_indices])

  proc.lls = lls[names(lls) %in% c("nll_M", "nll_NAA", "nll_recruit", "lprior_b")]
  names(proc.lls) = c("M", "NAA", "recruit", "W_b_M")[match(names(proc.lls),c("nll_M", "nll_NAA", "nll_recruit", "lprior_b"))]
  proc.lls = unlist(lapply(proc.lls, sum))
  n.proc.ll = length(proc.lls)
  proc.dists = rep("log(x) ~ Gaussian", n.proc.ll)
  likes = -c(obs.lls, proc.lls)
  n.likes = length(likes)
  my.range <- 1.2*range(likes)#c(min(likes), 1.2*max(likes))
  par(mar=c(5,10,1,1), oma=c(1,0,0,0))
  if(missing(plot.colors)) plot.colors = mypalette(n.likes)
  barplot(horiz=TRUE, likes, beside=FALSE, col=plot.colors, xlab="Conditional log-likelihood components",  axisnames=FALSE,  axes=FALSE,  space=0,
    xlim=my.range)
  axis(side=1, at=pretty(seq(my.range[1],my.range[2]), n=10), labels=pretty(seq(my.range[1],my.range[2]), n=10), cex=.75 )
  axis(side=2, at=seq(0.5,(n.likes-0.5)), labels= names(likes), las=2)
  axis(side=2, at=seq(0.5,(n.likes-0.5))-0.2, labels = c(obs.dists,proc.dists), las = 2, tick = FALSE)
  text(x= likes, y=seq(0.5,(n.likes-0.5)), labels=round(likes,0), cex=0.8, pos=ifelse(likes>0, 4, 2))
  box()
  #title(paste0("Components of Obj. Function (", round(as.numeric(asap$like[1]),0), "), npar=", npar), cex=0.9 )
  title(sub=paste0("Model: ", mod$model_name, "     ", mod$date))
}

residual.bubbleplot.fn = function(x, x.cats, y.cats, scale.bubble = 0.25*25, ylab = "Pearson Residuals", xlab = "Age"){
  resids <- x  # NOTE obs-pred
  n_y = NCOL(x)
  n_x = NROW(x)
  pos.col = "#ffffffaa"
  neg.col = "#8c8c8caa"
  range.resids<-range(abs((as.vector(x))), na.rm=T)
  #scale.resid.bubble <- 25

  z <- x * scale.bubble #* scale.resid.bubble
  resid.col = ifelse(z > 0.0, pos.col, neg.col)

  plot(1:n_x, 1:n_y,  xlim = c(1,n_x), ylim = c(1,n_y), xlab = x_lab, ylab = ylab, type = "n", axes=F)
  axis(1, lab = x.cats, lwd = 2)
  axis(2, lab = y.cats, cex.axis=0.75, las=1, lwd = 2)
  box(lwd = 2)
  abline(h = 1:n_y, col="lightgray")
  segments(x0 = 1:n_x, y0 = rep(1,n_x), x1= 1:n_x, y1 = rep(n_y,n_x), col = "lightgray", lty = 1)
  for (j in 1:n_y)
  {
    points(1:n_x, rep(j, n_x), cex = abs(z[j,]), col="black", bg = resid.col[j,],  pch = 21)
  }

  bubble.legend1 = c(0.1, 1, 2)
  #bubble.legend1 <- c(0.01,0.05,0.1)
  bubble.legend2 <- bubble.legend1 * scale.bubble#scale.resid.bubble*scale.bubble
  legend("topright", xpd = TRUE, legend = bubble.legend1, pch = rep(1, 3), pt.cex = bubble.legend2, horiz = TRUE, col = 'black')
  legend("topleft", xpd = TRUE, legend = c("Neg.", "Pos."), pch = rep(21, 2), pt.cex = 3, horiz = TRUE,
    pt.bg = c(neg.resid.col, pos.resid.col), col = "black")
  text(x = trunc(n_x/2), y = 0, cex = 0.8, label=paste0("Max(abs(resid)) = ",round(range.resids[2],2)))
}

hi.cor.fn <- function(mod, out.dir = "", do.tex = FALSE, do.csv = FALSE, cor.limit = 0.99)
{
  cor = mod$sdrep$cov.fixed
  cor = cor/sqrt(diag(cor) %*% t(diag(cor)))
  npar = NROW(cor)
	which.hi <- matrix(nrow = 0, ncol = 2)
	for(i in 2:npar) for(j in 1:(i-1))  if(cor[i,j]>cor.limit) which.hi = rbind(which.hi, c(i, j))

  hi.cor = cor[which.hi[,1],which.hi[,2]]
  if(length(which.hi))
  {
    out = cbind.data.frame(
      "Parameter 1" = rownames(cor)[which.hi[,1]], "Parameter 1 value" = round(mod$opt$par[which.hi[,1]],3),
      "Parameter 2" = colnames(cor)[which.hi[,2]], "Parameter 2 value" = round(mod$opt$par[which.hi[,2]],3),
      "Correlation" = round(cor[which.hi],3))

    if(do.tex) x = latex(out, file = paste0(out.dir,"hi_cor.tex"), rowlabel = '', table.env = FALSE)
    if(do.csv) write.csv(out, file = paste0(out.dir,"hi_cor.csv"))
    return(out)
  }
  else cat(paste0("No fixed effects had correlations estimated greater than ", cor.limit, " \n"))
}

get.RMSEs.fn <- function(model)
{
  out <- list(RMSE=list(), RMSE_n = list())
  if(class(mod$sdrep)[1] == "sdreport"){
    sdrep = summary(model$sdrep)
  } else {
    sdrep = model$sdrep
  }  
  temp = sdrep[rownames(sdrep) %in% "log_catch_resid",]
  catch_stdresid <- matrix(temp[,1]/temp[,2], model$env$data$n_years_model, model$env$data$n_fleets)
  temp = sdrep[rownames(sdrep) %in% "log_index_resid",]
  index_stdresid <- matrix(temp[,1]/temp[,2], model$env$data$n_years_model, model$env$data$n_indices)
  # temp = model$env$data$catch_paa - aperm(model$rep$pred_catch_paa[1:model$env$data$n_years_model,,,drop=FALSE],c(2,1,3))
  #temp = model$env$data$catch_paa - aperm(model$rep$pred_catch_paa,c(2,1,3))

  out$RMSE$catch <- sqrt(apply(catch_stdresid^2,2,mean, na.rm = TRUE))
  out$RMSE_n$catch = apply(catch_stdresid^2,2,function(x) sum(!is.na(x)))
  out$RMSE$catch_total = sqrt(mean(catch_stdresid^2, na.rm = TRUE))
  out$RMSE_n$catch_total <- sum(!is.na(catch_stdresid^2))
  out$RMSE$index <- sqrt(apply(index_stdresid^2,2,mean, na.rm = TRUE))
  out$RMSE_n$index = apply(index_stdresid^2,2,function(x) sum(!is.na(x)))
  out$RMSE$index_total = sqrt(mean(index_stdresid^2, na.rm = TRUE))
  out$RMSE_n$index_total <- sum(!is.na(index_stdresid^2))
  return(out)
}
RMSE.table.fn <- function(model)
{
  origpar <- par(no.readonly = TRUE)
  RMSEs = get.RMSEs.fn(model)
	par(mfrow=c(1,1), oma=rep(2,4), mar=c(0,0,1,0))
	max.txt<-16
	rmses <- which(unlist(RMSEs$RMSE)>0)
	n.rmses<-length(rmses)
	plot(seq(1,15), seq(1,15), type='n', axes=F, xlab="", ylab="", xlim=c(1,max.txt+2+ 6+2+ 8+2), ylim=c(n.rmses+4, 1))
	text(rep(1, n.rmses), seq(3,n.rmses+2), labels=names(unlist(RMSEs$RMSE)[rmses]), pos=4)
	text(rep(max.txt+2, n.rmses), seq(3,n.rmses+2), labels=unlist(RMSEs$RMSE_n)[rmses], pos=4)
	text(rep(max.txt+2+ 6+2, n.rmses), seq(3,n.rmses+2), labels=signif(as.numeric(unlist(RMSEs$RMSE)[rmses]),3), pos=4)
	text(c(1, max.txt+2, max.txt+2+ 6+2), rep( 2, 3), labels=c("Component","# resids","RMSE"), font=2, pos=4)
	title(main="Root Mean Square Error computed from Standardized Residuals", outer=T, cex=0.85)
  par(origpar)
	#if (save.plots) savePlot(paste(od, "RMSE_Comp_Table.",plotf, sep=""), type=plotf)
}
#RMSE.table.fn(mod)
#plot.RMSE.table(x)

get.wham.results.fn = function(mod, out.dir, do.tex = FALSE, do.png = FALSE)
{
  years = mod$years #not used in cpp
  ages = mod$ages #not used in cpp, assumes last age is a plus group
  ny = length(years)
  na = length(ages)
  ni = mod$env$data$n_indices
  nf = mod$env$data$n_fleets

  res = list()
  res$ll = -mod$opt$obj
  res$np = length(mod$opt$par)
  res$aic = 2*(mod$opt$obj + res$np)
  #rho = mohns_rho(mod) #if no peels, then this will be NULL
  tcol <- col2rgb('black')
  black.poly <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
  tcol <- col2rgb('red')
  red.poly <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
  tcol <- col2rgb('blue')
  blue.poly <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')

  if(do.tex | do.png)
  {
    use_outer = TRUE
    x_line = y_line = 2.5
  }
  else
  {
    use_outer = FALSE
    x_line = -1
    y_line = 3
  }
  origpar <- par(no.readonly = TRUE)
  par(mar = c(5,5,1,1), oma = c(1,1,1,1), mfrow = c(3,1))
  if(do.tex | do.png)
  {
    fn = paste0(out.dir, '/SSB')
    if(do.tex) cairo_pdf(paste0(fn,'.pdf'), family = fontfam, height = 10, width = 10)
    else png(filename = paste0(fn, '.png'), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar = c(0,0,0,0), oma = c(4,4,1,1), mfrow = c(1,1))
  }
  if(class(mod$sdrep)[1] == "sdreport"){
    temp = summary(mod$sdrep)
  } else {
    temp = mod$sdrep
  }
  temp = temp[rownames(temp) == "log_SSB",]
  temp = exp(cbind(temp, temp[,1] + qnorm(0.975)*cbind(-temp[,2],temp[,2])))/1000
  max.y <- max(temp[,4])
  na.lim <- is.na(max.y)
  if(!na.lim){
    plot(years,temp[,1], type = 'n', ylim = c(0,max.y), xlab = "", ylab = '', axes = FALSE)
    axis(1, lwd = 2, cex.axis = 1.5)
    axis(2, lwd = 2, cex.axis = 1.5)
    grid(col = gray(0.7))
    lines(years,temp[,1], lwd = 2)
    polygon(c(years,rev(years)), c(temp[,3],rev(temp[,4])), col = black.poly, border = "transparent")
    box(lwd = 2)
    mtext(side = 1, "Year", cex = 2, outer = TRUE, line = x_line)
    mtext(side = 2, "SSB (kmt)", cex = 2, outer = use_outer, line = y_line)
  } else {
    max.y <- max(temp[,1])
    plot(years,temp[,1], type = 'n', ylim = c(0,max.y), xlab = "", ylab = '', axes = FALSE)
    axis(1, lwd = 2, cex.axis = 1.5)
    axis(2, lwd = 2, cex.axis = 1.5)
    grid(col = gray(0.7))
    lines(years,temp[,1], lwd = 2)
    # polygon(c(years,rev(years)), c(temp[,3],rev(temp[,4])), col = black.poly, border = "transparent")
    box(lwd = 2)
    mtext(side = 1, "Year", cex = 2, outer = TRUE, line = x_line)
    mtext(side = 2, "SSB (kmt)", cex = 2, outer = use_outer, line = y_line)
  }
  if(do.tex | do.png) dev.off() else par(origpar)

  if(do.tex | do.png)
  {
    fn = paste0(out.dir, '/recruits')
    if(do.tex) cairo_pdf(paste0(fn,'.pdf'), family = fontfam, height = 10, width = 10)
    else png(filename = paste0(fn, '.png'), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar = c(0,0,0,0), oma = c(4,4,1,1), mfrow = c(1,1))
  }
  if(class(mod$sdrep)[1] == "sdreport"){
    temp = summary(mod$sdrep)
  } else {
    temp = mod$sdrep
  }
  ind = rownames(temp) == "log_NAA_rep"
  templo = exp(array(temp[ind,1] - qnorm(0.975)*temp[ind,2], dim = c(ny, na)))
  temphi = exp(array(temp[ind,1] + qnorm(0.975)*temp[ind,2], dim = c(ny, na)))
  temp = exp(array(temp[ind,1], dim = c(ny, 8)))
  max.y <- max(temphi[,1])
  na.lim <- is.na(max.y)
  if(!na.lim){
    plot(years,temp[,1], type = 'n', ylim = c(0,max.y), xlab = "", ylab = '', axes = FALSE)
    axis(1, lwd = 2, cex.axis = 1.5)
    axis(2, lwd = 2, cex.axis = 1.5)
    grid(col = gray(0.7))
    lines(years,temp[,1], lwd = 2)
    polygon(c(years,rev(years)), c(templo[,1],rev(temphi[,1])), col = black.poly, border = "transparent")
    box(lwd = 2)
    if(use_outer) mtext(side = 1, "Year", cex = 2, outer = TRUE, line = x_line)
    mtext(side = 2, "Recruits (1000s)", cex = 2, outer = use_outer, line = y_line)
  } else {
    max.y <- max(temp[,1])
    plot(years,temp[,1], type = 'n', ylim = c(0,max.y), xlab = "", ylab = '', axes = FALSE)
    axis(1, lwd = 2, cex.axis = 1.5)
    axis(2, lwd = 2, cex.axis = 1.5)
    grid(col = gray(0.7))
    lines(years,temp[,1], lwd = 2)
    # polygon(c(years,rev(years)), c(templo[,1],rev(temphi[,1])), col = black.poly, border = "transparent")
    box(lwd = 2)
    if(use_outer) mtext(side = 1, "Year", cex = 2, outer = TRUE, line = x_line)
    mtext(side = 2, "Recruits (1000s)", cex = 2, outer = use_outer, line = y_line)
  }
  if(do.tex | do.png) dev.off() else par(origpar)

  if(do.tex | do.png)
  {
    fn = paste0(out.dir, '/Fbar')
    if(do.tex) cairo_pdf(paste0(fn,'.pdf'), family = fontfam, height = 10, width = 10)
    else png(filename = paste0(fn, '.png'), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar = c(0,0,0,0), oma = c(4,4,1,1), mfrow = c(1,1))
  }
  if(class(mod$sdrep)[1] == "sdreport"){
    temp = summary(mod$sdrep)
  } else {
    temp = mod$sdrep
  }
  temp = temp[rownames(temp) == "log_Fbar",]
  temp = exp(cbind(temp, temp[,1] + qnorm(0.975)*cbind(-temp[,2],temp[,2])))
  max.y <- max(temp[,4])
  na.lim <- is.na(max.y)
  if(!na.lim){
    plot(years,temp[,1], type = 'n', ylim = c(0,max.y), xlab = "", ylab = '', axes = FALSE)
    axis(1, lwd = 2, cex.axis = 1.5)
    axis(2, lwd = 2, cex.axis = 1.5)
    grid(col = gray(0.7))
    lines(years,temp[,1], lwd = 2)
    polygon(c(years,rev(years)), c(temp[,3],rev(temp[,4])), col = black.poly, border = "transparent")
    box(lwd = 2)
    if(use_outer) mtext(side = 1, "Year", cex = 2, outer = TRUE,line = x_line)
    ar = c(min(mod$env$data$Fbar_ages), max(mod$env$data$Fbar_ages))
    mtext(side = 2, paste0("Average F (",ar[1],"-",ar[2],")"), cex = 2, outer = use_outer, line = y_line)
  } else {
    max.y <- max(temp[,1])
    plot(years,temp[,1], type = 'n', ylim = c(0,max.y), xlab = "", ylab = '', axes = FALSE)
    axis(1, lwd = 2, cex.axis = 1.5)
    axis(2, lwd = 2, cex.axis = 1.5)
    grid(col = gray(0.7))
    lines(years,temp[,1], lwd = 2)
    # polygon(c(years,rev(years)), c(temp[,3],rev(temp[,4])), col = black.poly, border = "transparent")
    box(lwd = 2)
    if(use_outer) mtext(side = 1, "Year", cex = 2, outer = TRUE,line = x_line)
    ar = c(min(mod$env$data$Fbar_ages), max(mod$env$data$Fbar_ages))
    mtext(side = 2, paste0("Average F (",ar[1],"-",ar[2],")"), cex = 2, outer = use_outer, line = y_line)
  }
  if(do.tex | do.png) dev.off() else par(origpar)
  # par(origpar)
}

plot.all.stdresids.fn = function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
  years = mod$years
  # load Ecov residuals
  xe <- NULL
  if(!all(mod$env$data$Ecov_model == 0)){
    ny = mod$env$data$n_years_Ecov
    ni = mod$env$data$n_Ecov
    if(class(mod$sdrep)[1] == "sdreport"){
      temp = summary(mod$sdrep)
    } else {
      temp = mod$sdrep
    }
    ind = rownames(temp) == "Ecov_resid"
    templo = matrix(temp[ind,1] - qnorm(0.975)*temp[ind,2], ny, ni)
    temphi = matrix(temp[ind,1] + qnorm(0.975)*temp[ind,2], ny, ni)
    temp = matrix(temp[ind,1], ny, ni)
    xe = data.frame(Label = integer(),
      Year = numeric(),
      Stdres = numeric(),
      lo = numeric(),
      hi = numeric())
    for(i in 1:ni)
    {
      ind = which(mod$env$data$Ecov_use_obs[,i] == 1)
      td = data.frame(Label = rep(i,length(ind)),
        Year = mod$input$data$Ecov_year[ind],
        Stdres = temp[ind,i],
        lo = templo[ind,i],
        hi = temphi[ind,i])
      xe <- rbind(xe, td)
    }
    xe$row = xe$Label
    xe$Label = factor(xe$Label)
    levels(xe$Label) = mod$env$data$Ecov_label[[1]]
    xe$type = "Ecov"
  }

  # load Index residuals
  ny = mod$env$data$n_years_model
  ni = mod$env$data$n_indices
  if(class(mod$sdrep)[1] == "sdreport"){
    temp = summary(mod$sdrep)
  } else {
    temp = mod$sdrep
  }
  ind = rownames(temp) == "log_index_resid"
  templo = matrix(temp[ind,1] - qnorm(0.975)*temp[ind,2], ny, ni)
  temphi = matrix(temp[ind,1] + qnorm(0.975)*temp[ind,2], ny, ni)
  temp = matrix(temp[ind,1], ny, ni)
  xi = data.frame(Label = integer(),
    Year = numeric(),
    Stdres = numeric(),
    lo = numeric(),
    hi = numeric())
  for(i in 1:ni)
  {
    ind = which(mod$env$data$use_indices[,i] == 1)
    td = data.frame(Label = rep(i,length(ind)),
      Year = years[ind],
      Stdres = temp[ind,i],
      lo = templo[ind,i],
      hi = temphi[ind,i])
    xi <- rbind(xi, td)
  }
  xi$row = xi$Label
  xi$Label = factor(xi$Label)
  levels(xi$Label) = paste0("Index ",1:length(table(xi$Label)))
  xi$type = "Index"
  # if(!is.null(index.names)) levels(x$Index) = index.names

  # load catch data (fleet)
  ny = mod$env$data$n_years_model
  ni = mod$env$data$n_fleets
  if(missing(years)) years = mod$years
  if(class(mod$sdrep)[1] == "sdreport"){
    temp = summary(mod$sdrep)
  } else {
    temp = mod$sdrep
  }
  ind = rownames(temp) == "log_catch_resid"
  templo = matrix(temp[ind,1] - qnorm(0.975)*temp[ind,2], ny, ni)
  temphi = matrix(temp[ind,1] + qnorm(0.975)*temp[ind,2], ny, ni)
  temp = matrix(temp[ind,1], ny, ni)
  xc = data.frame(Label = integer(),
    Year = numeric(),
    Stdres = numeric(),
    lo = numeric(),
    hi = numeric())
  for(i in 1:ni)
  {
    td = data.frame(Label = rep(i,ny),
      Year = years,
      Stdres = temp[,i],
      lo = templo[,i],
      hi = temphi[,i])
    xc <- rbind(xc, td)
  }
  xc$row = xc$Label
  xc$Label = factor(xc$Label)
  levels(xc$Label) = paste0("Fleet ",1:length(table(xc$Label)))
  xc$type = "Catch"
  # if(!is.null(fleet.names)) levels(xc$Fleet) = fleet.names

  x <- rbind(xe, xi, xc)
  x$row = factor(x$row)
  x$type = factor(x$type)

  ggp = ggplot2::ggplot(x, ggplot2::aes(x=Year, y = Stdres, color=type, fill=type)) +
    # ggplot2::geom_ribbon(ggplot2::aes(ymin=lo, ymax=hi, fill=type), alpha=0.3, linetype = 0) +
    # ggplot2::geom_line(size=1.1) +
    #ggplot2::geom_smooth(method = "lm", alpha=0.2) +
    ggplot2::geom_smooth(formula = y ~ x, method = "lm", alpha=0.2) +
    ggplot2::geom_point(size=0.8) +
    ggplot2::ylab("Standardized residual") +
    # expand_limits(y=0) +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "none") +
    # ggplot2::scale_color_manual(values=plot.colors) +
    # ggplot2::scale_fill_manual(values=plot.colors) +
    ggplot2::facet_wrap(~Label)
    # ggplot2::facet_grid(type ~ row)
    # ggplot2::facet_grid(row ~ type)
  if(do.tex) cairo_pdf(file.path(od, paste0("Residuals_time.pdf")), family = fontfam, height = 10, width = 10)
  if(do.png) png(filename = file.path(od, paste0("Residuals_time.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
  print(ggp)
  if(do.tex | do.png) dev.off()
  # return(ggp)
}

plot.ecov.stdresids.fn = function(mod, years, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, plot.colors, od)
{
  ny = mod$env$data$n_years_Ecov
  ni = mod$env$data$n_Ecov
  if(missing(years)) years = mod$years
  if(missing(plot.colors)) plot.colors = mypalette(ni)
  if(class(mod$sdrep)[1] == "sdreport"){
    temp = summary(mod$sdrep)
  } else {
    temp = mod$sdrep
  }
  ind = rownames(temp) == "Ecov_resid"
  templo = matrix(temp[ind,1] - qnorm(0.975)*temp[ind,2], ny, ni)
  temphi = matrix(temp[ind,1] + qnorm(0.975)*temp[ind,2], ny, ni)
  temp = matrix(temp[ind,1], ny, ni)
  x = data.frame(Ecov = integer(),
    Year = numeric(),
    Stdres = numeric(),
    lo = numeric(),
    hi = numeric())
  for(i in 1:ni)
  {
    ind = which(mod$env$data$Ecov_use_obs[,i] == 1)
    td = data.frame(Ecov = rep(i,length(ind)),
      Year = years[ind],
      Stdres = temp[ind,i],
      lo = templo[ind,i],
      hi = temphi[ind,i])
    x <- rbind(x, td)
  }
  x$Ecov = factor(x$Ecov)
  levels(x$Ecov) = mod$env$data$Ecov_label[[1]]
  names(plot.colors) = levels(x$Ecov)
  ggp = ggplot2::ggplot(x, ggplot2::aes(x=Year, y = Stdres, color=Ecov)) +
    ggplot2::geom_line(size=1.1) +
    ggplot2::geom_ribbon(ggplot2::aes(ymin=lo, ymax=hi, fill=Ecov), alpha=0.3, linetype = 0) +
    ggplot2::ylab("Standardized residual") +
    # expand_limits(y=0) +
    ggplot2::theme_bw() +
    ggplot2::scale_color_manual(values=plot.colors) +
    ggplot2::scale_fill_manual(values=plot.colors) +
    ggplot2::facet_wrap(~Ecov, ncol=1)
  if(do.tex) cairo_pdf(file.path(od, paste0("Residuals_ecov_time.pdf")), family = fontfam, height = 10, width = 10)
  if(do.png) png(filename = file.path(od, paste0("Residuals_ecov_time.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
  print(ggp)
  if(do.tex | do.png) dev.off()
  # return(ggp)
}

plot.index.stdresids.fn = function(mod, years, index.names = NULL, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, plot.colors, od)
{
  ny = mod$env$data$n_years_model
  ni = mod$env$data$n_indices
  if(missing(years)) years = mod$years
  if(missing(plot.colors)) plot.colors = mypalette(ni)
  if(class(mod$sdrep)[1] == "sdreport"){
    temp = summary(mod$sdrep)
  } else {
    temp = mod$sdrep
  }
  ind = rownames(temp) == "log_index_resid"
  templo = matrix(temp[ind,1] - qnorm(0.975)*temp[ind,2], ny, ni)
  temphi = matrix(temp[ind,1] + qnorm(0.975)*temp[ind,2], ny, ni)
  temp = matrix(temp[ind,1], ny, ni)
  x = data.frame(Index = integer(),
    Year = numeric(),
    Stdres = numeric(),
    lo = numeric(),
    hi = numeric())
  for(i in 1:ni)
  {
    ind = which(mod$env$data$use_indices[,i] == 1)
    td = data.frame(Index = rep(i,length(ind)),
      Year = years[ind],
      Stdres = temp[ind,i],
      lo = templo[ind,i],
      hi = temphi[ind,i])
    x <- rbind(x, td)
  }
  x$Index = factor(x$Index)
  names(plot.colors)= levels(x$Index)
  if(!is.null(index.names)) levels(x$Index) = index.names
  ggp = ggplot2::ggplot(x, ggplot2::aes(x=Year, y = Stdres, color=Index)) +
    ggplot2::geom_line(size=1.1) +
    ggplot2::geom_ribbon(ggplot2::aes(ymin=lo, ymax=hi, fill=Index), alpha=0.3, linetype = 0) +
    ggplot2::ylab("Standardized residual") +
    ggplot2::expand_limits(y=0) +
    ggplot2::theme_bw() +
    ggplot2::scale_color_manual(values=plot.colors) +
    ggplot2::scale_fill_manual(values=plot.colors) +
    ggplot2::facet_wrap(~Index, ncol=1)
  if(do.tex) cairo_pdf(file.path(od, paste0("Residuals_log_index_time.pdf")), family = fontfam, height = 10, width = 10)
  if(do.png) png(filename = file.path(od, paste0("Residuals_log_index_time.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
  print(ggp)
  if(do.tex | do.png) dev.off()
  # return(ggp)
}

plot.fleet.stdresids.fn = function(mod, years, fleet.names = NULL, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, plot.colors, od)
{
  ny = mod$env$data$n_years_model
  ni = mod$env$data$n_fleets
  if(missing(years)) years = mod$years
  if(missing(plot.colors)) plot.colors = mypalette(ni)
  if(class(mod$sdrep)[1] == "sdreport"){
    temp = summary(mod$sdrep)
  } else {
    temp = mod$sdrep
  }
  ind = rownames(temp) == "log_catch_resid"
  templo = matrix(temp[ind,1] - qnorm(0.975)*temp[ind,2], ny, ni)
  temphi = matrix(temp[ind,1] + qnorm(0.975)*temp[ind,2], ny, ni)
  temp = matrix(temp[ind,1], ny, ni)
  x = data.frame(Fleet = integer(),
    Year = numeric(),
    Stdres = numeric(),
    lo = numeric(),
    hi = numeric())
  for(i in 1:ni)
  {
    td = data.frame(Fleet = rep(i,ny),
      Year = years,
      Stdres = temp[,i],
      lo = templo[,i],
      hi = temphi[,i])
    x <- rbind(x, td)
  }
  x$Fleet = factor(x$Fleet)
  if(!is.null(fleet.names)) levels(x$Fleet) = fleet.names
  names(plot.colors)= levels(x$Fleet)
  ggp = ggplot2::ggplot(x, ggplot2::aes(x=Year, y = Stdres, color=Fleet)) +
    ggplot2::geom_line(size=1.1) +
    ggplot2::scale_color_manual(values=plot.colors) +
    ggplot2::scale_fill_manual(values=plot.colors) +
    ggplot2::geom_ribbon(ggplot2::aes(ymin=lo, ymax=hi, fill=Fleet), alpha=0.3, linetype = 0) +
    ggplot2::ylab("Standardized residual") +
    ggplot2::expand_limits(y=0) +
    ggplot2::theme_bw() +
    ggplot2::facet_wrap(~Fleet, ncol=1)
  if(do.tex) cairo_pdf(file.path(od, paste0("Residuals_log_catch_time.pdf")), family = fontfam, height = 10, width = 10)
  if(do.png) png(filename = file.path(od, paste0("Residuals_log_catch_time.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
  print(ggp)
  if(do.tex | do.png) dev.off()
  # return(ggp)
}

plot.catch.4.panel <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, use.i, plot.colors, od)
{
  origpar <- par(no.readonly = TRUE)
  years <- mod$years
  dat = mod$env$data
  years_full = mod$years_full
  pred_log_catch = mod$rep$pred_log_catch
  pred_catch = exp(pred_log_catch)
  sigma = dat$agg_catch_sigma %*% diag(exp(mod$parList$log_catch_sig_scale)) # dims: [ny,nf] x [nf]
  catch = dat$agg_catch
  log_stdres = (log(catch) - pred_log_catch[1:length(years),])/sigma # cpp already bias-corrects if bias_correct_oe = 1
  if(!missing(use.i)) fleets <- use.i
  else fleets <- 1:dat$n_fleets
  if(missing(plot.colors)) plot.colors = mypalette(dat$n_fleets)
	for (i in fleets)
	{
		if(do.tex) cairo_pdf(file.path(od, paste0("Catch_4panel_fleet_",i,".pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("Catch_4panel_fleet_",i,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(2,2))
		plot(years_full, pred_catch[,i], col=plot.colors[i], lwd=2, type='l', xlab="Year", ylab="Total Catch",
			ylim=c(0, 1.1*max(c(catch[,i],pred_catch[,i]))))
    points(years, catch[,i], col=plot.colors[i], pch=1)
    if(mod$env$data$n_fleets == 1){
      if(length(years_full) > length(years)){
        abline(v=tail(years,1), lty=2, lwd=1)
      }
    }
		log.ob.min <- log(catch[,i])-1.96*sigma[,i]
		log.ob.max <- log(catch[,i])+1.96*sigma[,i]
		plot(years_full, log(pred_catch[,i]), col=plot.colors[i], lwd=2, type='l', xlab="Year", ylab="Ln(Total Catch)",
			ylim=c(min(log.ob.min,log(pred_catch[,i]), na.rm=T), 1.1*max(log.ob.max,log(pred_catch[,i]), na.rm=T)))
		points(years, log(catch[,i]), pch=1, col=plot.colors[i])
    if(mod$env$data$n_fleets == 1) abline(v=tail(years,1), lty=2, lwd=1)
		arrows(years, log.ob.min, years, log.ob.max, length=0)
		title (paste0("Fleet ",i, " Catch"), outer=T, line=-1)
		plot(years, log_stdres[,i], type='h', lwd=2, col=plot.colors[i], xlab="Year", ylab="Log-scale Std. Residual")
		abline(h=0)
		hist(log_stdres[,i], plot=T, xlab="Std. Residual", ylab="Probability Density", freq=F, main=NULL)
		if(do.tex | do.png) dev.off() else par(origpar)
	}
}

plot.index.4.panel <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, use.i, plot.colors, od)
{
  origpar <- par(no.readonly = TRUE)
  years <- mod$years
  dat = mod$env$data
  pred_index = exp(mod$rep$pred_log_indices[1:length(years),,drop=F])
  index = dat$agg_indices
  # index[index < 0] = NA # robustify to missing values entered as negative
  index[dat$use_indices == 0] = NA # don't plot unused values
  sigma = dat$agg_index_sigma %*% diag(exp(mod$parList$log_index_sig_scale)) # dims: [ny,ni] x [ni]
  log_stdres = (log(index)-log(pred_index))/sigma
  if(!missing(use.i)) indices <- use.i
  else indices <- 1:dat$n_indices
  if(missing(plot.colors)) plot.colors = mypalette(dat$n_indices)
	for (i in indices)
	{
		if(do.tex) cairo_pdf(file.path(od, paste0("Index_4panel_",i,".pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("Index_4panel_",i,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(2,2))
		plot(years, index[,i], type='p', col=plot.colors[i], pch=1, xlab="Year", ylab="Index",
			ylim=c(0, 1.1*max(index[,i], na.rm=T)))
		lines(years, pred_index[,i], col=plot.colors[i], lwd=2)
		log.ob.min <- log(index[,i])-1.96*sigma[,i]
		log.ob.max <- log(index[,i])+1.96*sigma[,i]
    y.min <- min(c(log.ob.min,log(pred_index[,i]))[is.finite(c(log.ob.min,log(pred_index[,i])))], na.rm=T)
    y.max <- 1.1*max(c(log.ob.max,log(pred_index[,i]))[is.finite(c(log.ob.max,log(pred_index[,i])))], na.rm=T)
		plot(years, log(index[,i]), type='p', col=plot.colors[i], pch=1, xlab="Year", ylab="Ln(Index)", ylim=c(y.min, y.max))
		lines(years, log(pred_index[,i]), col=plot.colors[i], lwd=2)
		arrows(years, log.ob.min, years, log.ob.max, length=0)
		title (paste0("Index ",i), outer=T, line=-1)
		plot(years, log_stdres[,i], type='h', lwd=2, col=plot.colors[i], xlab="Year", ylab="Log-scale Std. Residual")
		abline(h=0)
		hist(log_stdres[,i], plot=T, xlab="Std. Residual", ylab="Probability Density", freq=F, main=NULL)
		if(do.tex | do.png) dev.off() else par(origpar)
	}
  # par(origpar)
}

plot.NAA.4.panel <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, use.i, plot.colors, od)
{
  origpar <- par(no.readonly = TRUE)
  par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(2,2))
  years <- mod$years
  years_full = mod$years_full
  dat = mod$env$data
  pred_NAA = mod$rep$pred_NAA
  NAA = mod$rep$NAA
  sigma = exp(mod$parList$log_NAA_sigma[mod$env$data$NAA_sigma_pointers])
  sigma = matrix(sigma, length(years_full), dat$n_ages, byrow = TRUE)
  log_stdres = (log(NAA)-log(pred_NAA))/sigma
  if(!missing(use.i)) ages <- use.i
  else ages <- 1:dat$n_ages
  if(missing(plot.colors)) plot.colors = mypalette(dat$n_ages)
	for (i in ages)
	{
		if(do.tex) cairo_pdf(file.path(od, paste0("NAA_4panel_",i,".pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("NAA_4panel_",i,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(2,2))
		plot(years_full, NAA[,i], type='p', col=plot.colors[i], pch=1, xlab="Year", ylab="Abundance (1000s)",
			ylim=c(0, 1.1*max(NAA[,i])))
		lines(years_full, pred_NAA[,i], col=plot.colors[i], lwd=2)
		if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
    log.ob.min <- log(NAA[,i])-1.96*sigma[,i]
		log.ob.max <- log(NAA[,i])+1.96*sigma[,i]
		plot(years_full, log(NAA[,i]), type='p', col=plot.colors[i], pch=1, xlab="Year", ylab="Ln(Abundance)",
			ylim=c(min(log.ob.min,log(pred_NAA[,i])), 1.1*max(log.ob.max,log(pred_NAA[,i]))))
		lines(years_full, log(pred_NAA[,i]), col=plot.colors[i], lwd=2)
		arrows(years_full, log.ob.min, years_full, log.ob.max, length=0)
    if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
		title (paste0("Conditional Expected and Posterior Estimates of Age ",i, " Abundance "), outer=T, line=-1)
		plot(years_full, log_stdres[,i], type='h', lwd=2, col=plot.colors[i], xlab="Year", ylab="Log-scale (Conditional) Std. Residual")
		abline(h=0)
		hist(log_stdres[,i], plot=T, xlab="(Conditional) Std. Residual", ylab="Probability Density", freq=F, main=NULL)
		if(do.tex | do.png) dev.off() else par(origpar)
	}
  # par(origpar)
}

plot.NAA.res <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, plot.colors, od)
{
  origpar <- par(no.readonly = TRUE)
  ages = mod$ages
	n.ages <- length(ages)
  if(missing(plot.colors)) plot.colors = mypalette(n.ages)
  years = mod$years_full
	n.yrs <- length(years)
  pred_NAA = mod$rep$pred_NAA
  NAA = mod$rep$NAA
  sigma = exp(mod$parList$log_NAA_sigma[mod$env$data$NAA_sigma_pointers])
  sigma = matrix(sigma, length(years), n.ages, byrow = TRUE)
  log_stdres = (log(NAA)-log(pred_NAA))/sigma

  ymin <- min(apply(log_stdres,1, function(x) sum(x[which(x<0)])))
  ymax <- max(apply(log_stdres,1, function(x) sum(x[which(x>0)])))
  dat <- as.data.frame.table(log_stdres)
  x.at <- seq(1,n.yrs,5)
  x.lab <- years[x.at]
  if(do.tex) cairo_pdf(file.path(od, paste0("NAA_res_barplot_stacked.pdf")), family = fontfam, height = 10, width = 10)
  if(do.png) png(filename = file.path(od, paste0("NAA_res_barplot_stacked.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
	par(mfrow=c(1,1), mar=c(5,5,1,1), oma = c(0,0,0,0))
	naa_res_plot <- lattice::barchart(Freq ~ Var1, data = dat, groups = Var2, stack = TRUE, col = plot.colors, xlab = "Year", ylab = "Std. Abundance Residuals", box.ratio = 10, reference = TRUE,
	  scales = list(x=list(at = x.at, labels = x.lab), alternating = FALSE),
	  key = list(text = list(lab = as.character(ages)), rectangles = list(col = plot.colors), columns = n.ages, title = "Age"))
  print(naa_res_plot)
  if(do.tex | do.png) dev.off() else par(origpar)
  # par(origpar)
}

plot.catch.age.comp <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, use.i, plot.colors, od)
{
  origpar <- par(no.readonly = TRUE)
  years = mod$years
  ages = 1:mod$env$data$n_ages
  ages.lab = mod$ages.lab
  if(!missing(use.i)) fleets <- use.i
  else fleets <- 1:mod$env$data$n_fleets
  if(missing(plot.colors)) plot.colors = mypalette(mod$env$data$n_fleets)

	for (i in fleets)
	{
    acomp.obs = mod$env$data$catch_paa[i,,]
    acomp.pred = aperm(mod$rep$pred_catch_paa[1:mod$env$data$n_years_model,,,drop=FALSE], c(2,1,3))[i,,]
    if(do.tex) cairo_pdf(file.path(od, paste0("Catch_age_comp_fleet_",i,".pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("Catch_age_comp_fleet_",i,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar=c(1,1,2,1), oma=c(4,4,2,1), mfcol=c(5,3))
    my.title <- paste0("Fleet ", i)
    for (j in 1:length(years))
    {
      plot(1:length(ages), acomp.obs[j,], type='p', col=plot.colors[i], pch=1, xlab="", ylab="",
        ylim=c(0, 1), axes = FALSE)
      if(j %% 15 == 1)
      {
        title(my.title, outer=TRUE, line=0)
        mtext(side = 2, "Proportion", line = 2, outer = TRUE)
        mtext(side = 1, "Age", line = 2, outer = TRUE)
      }
      if(j %% 15 %in% 1:5) axis(2)
      else axis(2, labels = FALSE)
      if(j %in% seq(5,length(years)+1,5)) axis(1, at = ages, labels = ages.lab)
      else axis(1, labels = FALSE)
      grid()
      box()
      lines(1:length(ages), acomp.pred[j,], col=plot.colors[i],  lwd=2)
      title(paste("Year = ", years[j], sep=""), outer = FALSE)

      # if 5x3 multipanel is full, save png and open new one
      if((j %% 15 == 0) & (do.tex | do.png) & (j < length(years))){
        dev.off()
        if(do.tex) cairo_pdf(file.path(od, paste0("Catch_age_comp_fleet_",i,"_",letters[j/15],".pdf")), family = fontfam, height = 10, width = 10)
        if(do.png) png(filename = file.path(od, paste0("Catch_age_comp_fleet_",i,"_",letters[j/15],".png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
        par(mar=c(1,1,2,1), oma=c(4,4,2,1), mfcol=c(5,3))
      }
    }  #end loop on n_years
    if(length(years) %% 15 != 0) frame()
    if(do.tex | do.png) dev.off() else par(origpar)
	}  #end loop on n_fleets
  # par(origpar)
}

plot.index.age.comp <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, use.i, plot.colors, od)
{
  origpar <- par(no.readonly = TRUE)
  years = mod$years
  ages = 1:mod$env$data$n_ages
  ages.lab = mod$ages.lab
  if(!missing(use.i)) indices <- use.i
  else indices <- 1:mod$env$data$n_indices
  if(missing(plot.colors)) plot.colors = mypalette(mod$env$data$n_indices)

	for (i in indices)
	{
    acomp.obs = mod$env$data$index_paa[i,,]
    acomp.pred = aperm(mod$rep$pred_IAA[1:length(years),,,drop=FALSE], c(2,1,3))[i,,] #biomass is accounted for on the cpp side
    acomp.pred = acomp.pred/apply(acomp.pred,1,sum)
    if(do.tex) cairo_pdf(file.path(od, paste0("Catch_age_comp_index_",i,".pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("Catch_age_comp_index_",i,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar=c(1,1,2,1), oma=c(4,4,2,1), mfcol=c(5,3))
    my.title <- paste0("Index ", i)
    for (j in 1:length(years))
    {
      plot(1:length(ages), acomp.obs[j,], type='p', col=plot.colors[i], pch=1, xlab="", ylab="",
        ylim=c(0, 1), axes = FALSE)
      if(j %% 15 == 1)
      {
        title(my.title, outer=TRUE, line=0)
        mtext(side = 2, "Proportion", line = 2, outer = TRUE)
        mtext(side = 1, "Age", line = 2, outer = TRUE)
      }
      if(j %% 15 %in% 1:5) axis(2)
      else axis(2, labels = FALSE)
      if(j %in% seq(5,length(years)+1,5)) axis(1, at = ages, labels = ages.lab)
      else axis(1, labels = FALSE)
      grid()
      box()
      lines(1:length(ages), acomp.pred[j,], col=plot.colors[i],  lwd=2)
      title(paste("Year = ", years[j], sep=""), outer = FALSE)

      # if 5x3 multipanel is full, save png and open new one
      if((j %% 15 == 0) & (do.tex | do.png) & (j < length(years))){
        dev.off()
        if(do.tex) cairo_pdf(file.path(od, paste0("Catch_age_comp_index_",i,"_",letters[j/15],".pdf")), family = fontfam, height = 10, width = 10)
        if(do.png) png(filename = file.path(od, paste0("Catch_age_comp_index_",i,"_",letters[j/15],".png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
        par(mar=c(1,1,2,1), oma=c(4,4,2,1), mfcol=c(5,3))
      }
    }  #end loop on n_years
    if(length(years) %% 15 != 0) frame()
		if(do.tex | do.png) dev.off() else par(origpar)
	}  #end loop on n_indices
  # par(origpar)
}

multinomial.pearson.fn = function(mod, ind = 1)
{
  dat = mod$env$data
  rep = mod$rep
  x = dat$index_paa[ind,,] - rep$pred_index_paa[1:dat$n_years_model,ind,]
  temp = dat$index_Neff[,ind]
  temp = rep$pred_index_paa[1:dat$n_years_model,ind,]*(1-rep$pred_index_paa[1:dat$n_years_model,ind,])/temp
  x = x/sqrt(temp)
  x[which(dat$use_index_paa[,ind] == 0),] = NA
  return(x)
}
#mean(x < 0, na.rm = TRUE)

plot.catch.age.comp.resids <- function(mod, ages, ages.lab, scale.catch.bubble2 = 2, pos.resid.col = "#ffffffaa", neg.resid.col = "#8c8c8caa",
                                       do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, osa = FALSE, use.i, od)
{
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = mod$ages.lab
	n_ages <- dat$n_ages
  years = mod$years
	nyrs <- length(years)
	if(!missing(use.i)) fleets <- use.i
	else fleets <- 1:mod$env$data$n_fleets
  tylab <- "Year"

	for (i in fleets)
	{
    yind = which(dat$use_catch_paa[,i] ==1)
    if(length(yind)){
      if(osa & "catchpaa" %in% mod$osa$type){
        df = subset(mod$osa, type == "catchpaa")
        my.title <- "Age Comp OSA Quantile Residuals for Fleet "
        fname = paste0("Catch_age_comp_osa_resids_fleet_",i)
        resids = matrix(NA, nrow = dat$n_years_model, ncol = dat$n_ages)
        vals = resids
        for(j in yind){
          tmp = subset(df, year == j & fleet == paste0("fleet_",i))
          resids[j,tmp$age] = tmp$residual
          vals[j,tmp$age] = tmp$val
          if(dat$age_comp_model_fleets[i] %in% c(1:2,10)) vals[j,tmp$age]/sum(vals[j,tmp$age]) #obs are numbers not proportions
        }

        scale.resid.bubble.catch <- 2
      } else{
        acomp.obs = dat$catch_paa[i,,]
        acomp.pred = mod$rep$pred_catch_paa[1:length(years),i,]
        #acomp.pred = aperm(mod$rep$pred_catch_paa[1:length(years),,,drop=FALSE], c(2,1,3))[i,,] #biomass is accounted for on the cpp side
        #acomp.pred = acomp.pred/apply(acomp.pred,1,sum)
        my.title <- "Age Comp Residuals (Observed-Predicted) for Fleet "
        resids <- acomp.obs - acomp.pred  # NOTE obs-pred
        resids[dat$use_catch_paa[,i]==0,] = NA # don't plot residuals for catch paa not fit in model
        fname = paste0("Catch_age_comp_resids_fleet_",i)
        scale.resid.bubble.catch <- 25
      }
      if(do.tex) cairo_pdf(file.path(od, paste0(fname,".pdf")), family = fontfam, height = 10, width = 10)
      if(do.png) png(filename = file.path(od, paste0(fname,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
      par(mar=c(4,4,2,2), oma=c(1,1,1,1), mfrow=c(1,1))
      z1 <- resids
      if(any(!is.na(resids))) range.resids<-range(abs((as.vector(z1))), na.rm=T)
      else range.resids = c(0,0)

      z3 <- z1 * scale.resid.bubble.catch * scale.catch.bubble2
      resid.col <- matrix(NA, nrow=nyrs, ncol=n_ages)   # set color for residual bubbles
      resid.col <- ifelse(z3 > 0.0, pos.resid.col, neg.resid.col)
      plot(ages, rev(ages),  xlim = c(1, n_ages+1), ylim = c(years[nyrs],(years[1]-2)), xlab = "Age", ylab = tylab,
        type = "n", axes=F)
      axis(1, at= ages, lab=ages.lab)
      axis(2, at = rev(years), lab = rev(years), cex.axis=0.75, las=1)
      box()
      abline(h=years, col="lightgray")
      segments(x0=ages, y0=rep(years[1],n_ages), x1=ages, y1=rep(years[nyrs],n_ages), col = "lightgray", lty = 1)

      for (j in 1:nyrs) points(ages, rep(years[j], n_ages), cex=abs(z3[j,]), col="black", bg = resid.col[j,],  pch = 21)

      bubble.legend1 <- round(quantile(abs(resids), probs = c(0.05,0.5,0.95), na.rm = TRUE),3)
      bubble.legend2 <- bubble.legend1 * scale.resid.bubble.catch*scale.catch.bubble2
      legend("topright", xpd=T, legend=bubble.legend1, pch=rep(1, 3), pt.cex=bubble.legend2, horiz=T , col='black')
      legend("topleft", xpd=T, legend=c("Neg.", "Pos."), pch=rep(21, 2), pt.cex=3, horiz=T, pt.bg=c(neg.resid.col, pos.resid.col), col="black")
      legend("top", xpd = TRUE, legend = paste("Max(resid)=",round(range.resids[2],2), sep=""), horiz = TRUE)
      title (paste0(my.title,i), outer=T, line=-1)
      if(do.tex | do.png) dev.off() else par(origpar)
    } #end if any age comp for fleet
	}   #end loop n_fleets
  # par(origpar)
}

plot.index.age.comp.resids <- function(mod, ages, ages.lab, scale.catch.bubble2 = 2, pos.resid.col = "#ffffffaa", neg.resid.col = "#8c8c8caa",
                                       do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, osa=FALSE, use.i, od)
{
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = mod$ages.lab
	n_ages <- dat$n_ages
  years = mod$years
	nyrs <- length(years)
	if(!missing(use.i)) indices <- use.i
	else indices <- 1:dat$n_indices

	for (i in indices)
	{
    yind = which(dat$use_index_paa[,i] ==1)
    if(length(yind)){
      if(osa & "indexpaa" %in% mod$osa$type){
        df = subset(mod$osa, type == "indexpaa")
        my.title <- "Age Comp OSA Quantile Residuals for Index "
        fname = paste0("Catch_age_comp_osa_resids_index_",i)
        resids = matrix(NA, nrow = dat$n_years_model, ncol = dat$n_ages)
        vals = resids
        for(j in yind){
          tmp = subset(df, year == j & fleet == paste0("index_",i))
          resids[j,tmp$age] = tmp$residual
          vals[j,tmp$age] = tmp$val
          if(dat$age_comp_model_indices[i] %in% c(1:2,10)) vals[j,tmp$age]/sum(vals[j,tmp$age]) #obs are numbers not proportions
        }
        scale.resid.bubble.catch <- 2
      } else {
        acomp.obs = dat$index_paa[i,,]
        acomp.pred = mod$rep$pred_index_paa[1:length(years),i,]
        #acomp.pred = aperm(mod$rep$pred_IAA[1:length(years),,,drop=FALSE], c(2,1,3))[i,,] #biomass is accounted for on the cpp side
        #acomp.pred = acomp.pred/apply(acomp.pred,1,sum)
        my.title <- "Age Comp Residuals (Observed-Predicted) for Index "
        resids <- acomp.obs - acomp.pred  # NOTE obs-pred
        resids[dat$use_index_paa[,i]==0,] = NA # don't plot residuals for index paa not fit in model
        fname = paste0("Catch_age_comp_resids_index_",i)
        scale.resid.bubble.catch <- 25
      }
      tylab <- "Year"
      
      if(do.tex) cairo_pdf(file.path(od, paste0(fname,".pdf")), family = fontfam, height = 10, width = 10)
      if(do.png) png(filename = file.path(od, paste0(fname,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
      par(mar=c(4,4,2,2), oma=c(1,1,1,1), mfrow=c(1,1))
      z1 <- resids
      if(any(!is.na(resids))) range.resids<-range(abs((as.vector(z1))), na.rm=T)
      else range.resids <- c(0,0)

      z3 <- z1 * scale.resid.bubble.catch * scale.catch.bubble2
      resid.col <- matrix(NA, nrow=nyrs, ncol=n_ages)   # set color for residual bubbles
      resid.col <- ifelse(z3 > 0.0, pos.resid.col, neg.resid.col)
      plot(ages, rev(ages),  xlim = c(1, n_ages+1), ylim = c(years[nyrs],(years[1]-2)), xlab = "Age", ylab = tylab,
        type = "n", axes=F)
      axis(1, at= ages, lab=ages.lab)
      axis(2, at = rev(years), lab = rev(years), cex.axis=0.75, las=1)
      box()
      abline(h=years, col="lightgray")
      segments(x0=ages, y0=rep(years[1],n_ages), x1=ages, y1=rep(years[nyrs],n_ages), col = "lightgray", lty = 1)

      for (j in 1:nyrs) if(dat$use_index_paa[j,i] == 1)
        points(ages, rep(years[j], n_ages), cex=abs(z3[j,]), col="black", bg = resid.col[j,],  pch = 21)

      bubble.legend1 <- round(quantile(abs(resids), probs = c(0.05,0.5,0.95), na.rm = TRUE),3)
      bubble.legend2 <- bubble.legend1 * scale.resid.bubble.catch*scale.catch.bubble2
      legend("topright", xpd=T, legend=bubble.legend1, pch=rep(1, 3), pt.cex=bubble.legend2, horiz=T , col='black')
      legend("topleft", xpd=T, legend=c("Neg.", "Pos."), pch=rep(21, 2), pt.cex=3, horiz=T, pt.bg=c(neg.resid.col, pos.resid.col), col="black")
      legend("top", xpd = TRUE, legend = paste("Max(resid)=",round(range.resids[2],2), sep=""), horiz = TRUE)
      title (paste0(my.title,i), outer=T, line=-1)
      if(do.tex | do.png) dev.off() else par(origpar)
    } #end if any age comp for index
	}   #end loop n_fleets
  # par(origpar)
}

plot.fleet.sel.blocks <- function(mod, ages, ages.lab, plot.colors, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, use.i, od)
{
  origpar <- par(no.readonly = TRUE)
  par(mfrow=c(1,1))
	cc<-0
  dat = mod$env$data
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = ages
  sb_p = dat$selblock_pointer_fleets #selblock pointer by year and fleet
  if(missing(plot.colors)) plot.colors = mypalette(length(unique(sb_p)))
	years <- mod$years
	if(!missing(use.i)) fleets <- use.i
	else fleets <- 1:mod$env$data$n_fleets

	for (i in fleets)
	{
	  if(do.tex) cairo_pdf(file.path(od, paste0("Selectivity_fleet",i,".pdf")), family = fontfam, height = 10, width = 10)
	  if(do.png) png(filename = file.path(od, paste0("Selectivity_fleet",i,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
	  blocks = unique(sb_p[,i])
		n.blocks <- length(blocks)
    # sel = rbind(mod$rep$selblocks[blocks,])
    sel = do.call(rbind, lapply(mod$rep$selAA, function(x) apply(x,2,mean)))[blocks,,drop=FALSE]
		minyr <- rep(NA, n.blocks)
		maxyr <- rep(NA, n.blocks)
		my.col <- rep(NA, n.blocks)
		for (j in 1:n.blocks)
		{
			cc<-cc+1
			my.col[j] <- plot.colors[cc]
			minyr[j] <- min(years[sb_p[,i] == blocks[j]])
			maxyr[j] <- max(years[sb_p[,i] == blocks[j]])
			if (j==1)
			{
				plot(ages, sel[j,], type='l', col=my.col[j], xlim=c(min(ages),max(ages)), ylim=c(0,1.1), xlab="Age", ylab="Selectivity",
					lwd=2, axes = FALSE)
				grid(col = gray(0.7), lwd = 2)
				axis(1, at = ages, labels = ages.lab, lwd = 2)
				axis(2, lwd = 2)
				box(lwd=2)
			}
			if (j>1) lines(ages, sel[j,], type='l', col=my.col[j], lwd=2)
		}
		title(paste0("Fleet ",i))
		legend("topright", col=my.col, legend=paste0(minyr, " - ", maxyr), lwd=2, bg = "white")
		if(do.tex | do.png) dev.off() else par(origpar)
	}
  # par(origpar)
}

plot.index.sel.blocks <- function(mod, ages, ages.lab, plot.colors, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, use.i, od)
{
  origpar <- par(no.readonly = TRUE)
  par(mfrow=c(1,1))
	cc<-0
  dat = mod$env$data
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = ages
  sb_p = dat$selblock_pointer_indices #selblock pointer by year and index
  if(missing(plot.colors)) plot.colors = mypalette(length(unique(sb_p)))
	years <- mod$years
	if(!missing(use.i)) indices <- use.i
	else indices <- 1:mod$env$data$n_indices

	for (i in indices)
	{
	  if(do.tex) cairo_pdf(file.path(od, paste0("Selectivity_index",i,".pdf")), family = fontfam, height = 10, width = 10)
	  if(do.png) png(filename = file.path(od, paste0("Selectivity_index",i,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
	  blocks = unique(sb_p[,i])
		n.blocks <- length(blocks)
    # sel = rbind(mod$rep$selblocks[blocks,])
    sel = do.call(rbind, lapply(mod$rep$selAA, function(x) apply(x,2,mean)))[blocks,,drop=FALSE]
		minyr <- rep(NA, n.blocks)
		maxyr <- rep(NA, n.blocks)
		my.col <- rep(NA, n.blocks)
		for (j in 1:n.blocks)
		{
			cc<-cc+1
			my.col[j] <- plot.colors[cc]
			minyr[j] <- min(years[sb_p[,i] == blocks[j]])
			maxyr[j] <- max(years[sb_p[,i] == blocks[j]])
			if (j==1)
			{
				plot(ages, sel[j,], type='l', col=my.col[j], xlim=c(min(ages),max(ages)), ylim=c(0,1.1), xlab="Age", ylab="Selectivity",
					lwd=2, axes = FALSE)
				grid(col = gray(0.7), lwd = 2)
				axis(1, at = ages, labels = ages.lab, lwd = 2)
				axis(2, lwd = 2)
				box(lwd=2)
			}
			if (j>1) lines(ages, sel[j,], type='l', col=my.col[j], lwd=2)
		}
		title(paste0("Index ",i))
		legend("topright", col=my.col, legend=paste0(minyr, " - ", maxyr), lwd=2, bg = "white")
		if(do.tex | do.png) dev.off() else par(origpar)
	}
  # par(origpar)
}

plot.q.trend<-function(mod, alpha = 0.05)
{
  origpar <- par(no.readonly = TRUE)
  years_full <- mod$years_full
  years <- mod$years
  tcol <- col2rgb('black')
  tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
  if(class(mod$sdrep)[1] == "sdreport"){
    std = summary(mod$sdrep)
  } else {
    std = mod$sdrep
  }
  par(mfrow=c(2,1), mar=c(1,1,1,1), oma = c(4,4,0,0))

  ssb.ind <- which(rownames(std) == "log_SSB")
  log.ssb <- std[ssb.ind,1]
  ssb = exp(log.ssb)/1000
  ssb.cv <- std[ssb.ind,2]
  log.ssb.ci <- log.ssb + cbind(qnorm(1-alpha/2)*ssb.cv, -qnorm(1-alpha/2)*ssb.cv)
  ssb.ci = exp(log.ssb.ci)/1000
  no.ssb.ci <- all(is.na(ssb.ci))
  if(!no.ssb.ci){ # have CI
    plot(years_full, ssb, type='l', lwd=2, xlab="", ylab="", ylim=c(0,max(ssb.ci)), axes = FALSE)
    axis(1, labels = FALSE)
    axis(2)
    box()
    mtext(side = 2, "SSB (kmt)", outer = FALSE, line = 3)
    grid(col = gray(0.7))
    polygon(c(years_full,rev(years_full)), c(ssb.ci[,1],rev(ssb.ci[,2])), col = tcol, border = tcol, lwd = 1)
  } else { # no CI but plot SSB trend
    plot(years_full, ssb, type='l', lwd=2, xlab="", ylab="", ylim=c(0,max(ssb)), axes = FALSE)
    axis(1, labels = FALSE)
    axis(2)
    box()
    mtext(side = 2, "SSB (kmt)", outer = FALSE, line = 3)
    grid(col = gray(0.7))
    # polygon(c(years,rev(years)), c(ssb.ci[,1],rev(ssb.ci[,2])), col = tcol, border = tcol, lwd = 1)
  }
  if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
  par(origpar)
}

plot.SSB.F.trend<-function(mod, alpha = 0.05)
{
  origpar <- par(no.readonly = TRUE)
  years_full <- mod$years_full
  years <- mod$years
  tcol <- col2rgb('black')
  tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
  if(class(mod$sdrep)[1] == "sdreport"){
    std = summary(mod$sdrep)
  } else {
    std = mod$sdrep
  }
	par(mfrow=c(2,1), mar=c(1,1,1,1), oma = c(4,4,0,0))

	ssb.ind <- which(rownames(std) == "log_SSB")
	log.ssb <- std[ssb.ind,1]
  ssb = exp(log.ssb)/1000
	ssb.cv <- std[ssb.ind,2]
  log.ssb.ci <- log.ssb + cbind(qnorm(alpha/2)*ssb.cv, qnorm(1-alpha/2)*ssb.cv)
  ssb.ci = exp(log.ssb.ci)/1000
  no.ssb.ci <- all(is.na(ssb.ci))
  if(!no.ssb.ci){ # have CI
    ymax <- max(ssb.ci[,2][which(!is.na(ssb.ci[,2]) & !is.infinite(ssb.ci[,2]))])
    ssb.ci[which(is.infinite(ssb.ci))] <-1e10
  	plot(years_full, ssb, type='l', lwd=2, xlab="", ylab="", ylim=c(0,ymax), axes = FALSE)
  	axis(1, labels = FALSE)
  	axis(2)
  	box()
  	mtext(side = 2, "SSB (kmt)", outer = FALSE, line = 3)
  	grid(col = gray(0.7))
  	polygon(c(years_full,rev(years_full)), c(ssb.ci[,1],rev(ssb.ci[,2])), col = tcol, border = tcol, lwd = 1)
	} else { # no CI but plot SSB trend
    ymax <- max(ssb[which(!is.na(ssb) & !is.infinite(ssb))])
    plot(years_full, ssb, type='l', lwd=2, xlab="", ylab="", ylim=c(0,ymax), axes = FALSE)
    axis(1, labels = FALSE)
    axis(2)
    box()
    mtext(side = 2, "SSB (kmt)", outer = FALSE, line = 3)
    grid(col = gray(0.7))
    # polygon(c(years,rev(years)), c(ssb.ci[,1],rev(ssb.ci[,2])), col = tcol, border = tcol, lwd = 1)
  }
  if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
  # F trend
  n_ages = mod$env$data$n_ages
	faa.ind <- which(rownames(std) == "log_FAA_tot")
	log.faa <- matrix(std[faa.ind,1], length(years_full), n_ages)
	faa.cv <- matrix(std[faa.ind,2], length(years_full), n_ages)
	age.full.f <- apply(log.faa,1, function(x) max(which(x == max(x))))
  full.f.ind = cbind(1:length(years_full), age.full.f)
	#full.f.ind <- c(age.full.f[1], n_ages + cumsum(age.full.f[-1]))
  log.full.f <- log.faa[full.f.ind]
  full.f.cv <- faa.cv[full.f.ind]
  log.f.ci <- log.full.f + cbind(qnorm(alpha/2)*full.f.cv, qnorm(1-alpha/2)*full.f.cv)
  full.f = exp(log.full.f)
  no.f.ci <- all(is.na(log.f.ci))
  if(!no.f.ci){ # have CI
    f.ci <- exp(log.f.ci)
    ymax <- max(f.ci[,2][which(!is.na(f.ci[,2]) & !is.infinite(f.ci[,2]))])
    f.ci[which(is.infinite(f.ci))] <-1e10
  	plot(years_full, full.f, type='l', lwd=2, col='black', xlab="", ylab="", ylim=c(0,ymax), axes = FALSE)
  	axis(1)
  	axis(2)
  	box()
  	mtext(side = 1, "Year", outer = FALSE, line = 3)
  	mtext(side = 2, "Fully-selected F", outer = FALSE, line = 3)
  	grid(col = gray(0.7))
  	polygon(c(years_full,rev(years_full)), c(f.ci[,1],rev(f.ci[,2])), col = tcol, border = tcol, lwd = 1)
  } else { # CI all NA
    ymax <- max(full.f[,1][which(!is.na(full.f[,1]) & !is.infinite(full.f[,1]))])
    plot(years_full, full.f, type='l', lwd=2, col='black', xlab="", ylab="", ylim=c(0,ymax), axes = FALSE)
    axis(1)
    axis(2)
    box()
    mtext(side = 1, "Year", outer = FALSE, line = 3)
    mtext(side = 2, "Fully-selected F", outer = FALSE, line = 3)
    grid(col = gray(0.7))
    # polygon(c(years,rev(years)), exp(c(log.f.ci[,1],rev(log.f.ci[,2]))), col = tcol, border = tcol, lwd = 1)
  }
  if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
  par(origpar)
}  #end function

plot.SSB.AA <- function(mod, ages, ages.lab, plot.colors, prop=FALSE)
{
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = mod$ages.lab
	n.ages <- length(ages)
  if(missing(plot.colors)) plot.colors = mypalette(n.ages)
	years<-  mod$years
  years_full <-  mod$years_full
	n.yrs <- length(years_full)
  ssbfrac = dat$fracyr_SSB
	ssb.aa <- (mod$rep$NAA * exp(-ssbfrac * (mod$rep$FAA_tot + mod$rep$MAA)) * dat$waa[dat$waa_pointer_ssb,,] * dat$mature)/1000
	ssb.max <- max(apply(ssb.aa,1,sum))

	par(mfrow=c(1,1), mar=c(5,5,1,1), oma = c(0,0,0,0))
	if(!prop){ # plot SSB at age
  	res <- barplot(t(ssb.aa), beside=F, cex.names=0.75, width=1, space=rep(0,n.yrs), xlab = 'Year', ylab =paste('SSB at age (', "kmt", ')', sep = ''),
  		ylim = c(0,1.15*ssb.max), xlim=c(0.5,n.yrs+1-0.5), col=plot.colors)
    if(length(years_full) > length(years)) abline(v=length(years), lty=2, lwd=2)
  	legend('top', horiz=TRUE, legend=ages.lab, pch=15, col=plot.colors, cex=0.8)
    axis(1, at = seq(5,n.yrs,5)-0.5, labels = years_full[seq(5,n.yrs,5)])
    box()
	}
	if(prop){ # plot *proportion* SSB at age
  	barplot(t(ssb.aa/apply(ssb.aa,1,sum)), beside=F, cex.names=0.75, width=1, space=rep(0,n.yrs), xlab = 'Year', ylab ='Proportion SSB at age',
  		ylim = c(0,1.1), xlim=c(0.5,n.yrs+1-0.5), col=plot.colors)
    if(length(years_full) > length(years)) abline(v=length(years), lty=2, lwd=2)
  	legend('top', horiz=TRUE, legend=ages.lab, pch=15, col=plot.colors, cex=0.8)
    axis(1, at = seq(5,n.yrs,5)-0.5, labels = years_full[seq(5,n.yrs,5)])
    box()
	}
  par(origpar)
}  #end funciton
#plot.SSB.AA(ssm)

#------------------------------------
plot.NAA <- function(mod, ages, ages.lab, plot.colors, scale = 1000, units = expression(10^6), prop=FALSE)
{
  origpar <- par(no.readonly = TRUE)
  par(mfrow=c(1,1))
  dat = mod$env$data
	## stacked barplot of NAA
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = mod$ages.lab
	n.ages <- length(ages)
  if(missing(plot.colors)) plot.colors = mypalette(n.ages)
  years<-  mod$years
  years_full <-  mod$years_full
  n.yrs <- length(years_full)
	NAA <- mod$rep$NAA
	N.max=max(apply(NAA,1,sum))/scale

	par(mfrow=c(1,1), mar=c(5,5,1,1), oma = c(0,0,0,0))
	if(!prop){ # plot numbers at age
  	barplot(t(NAA)/scale, beside=F, cex.names=0.75, width=1, space=rep(0,n.yrs), xlab = 'Year',
  		ylab =as.expression(substitute(paste("January 1 numbers at age (", units, ")", sep = ''), list(units = units[[1]]))),
  		ylim = c(0,1.15*N.max), xlim=c(0.5,n.yrs+1-0.5), col=plot.colors)
    if(length(years_full) > length(years)) abline(v=length(years), lty=2, lwd=2)
  	legend('top', horiz=TRUE, legend=ages.lab, pch=15, col=plot.colors, cex=0.8)
    axis(1, at = seq(5,n.yrs,5)-0.5, labels = years_full[seq(5,n.yrs,5)])
    box()
	}
	if(prop){ # plot *proportion* of numbers at age
  	barplot(t(NAA/apply(NAA,1,sum)), beside=F, cex.names=0.75, width=1, space=rep(0,n.yrs), xlab = 'Year',
  		ylab ='January 1 proportions at age', ylim = c(0,1.1), xlim=c(0.5,n.yrs+1-0.5), col=plot.colors)
    if(length(years_full) > length(years)) abline(v=length(years), lty=2, lwd=2)
  	legend('top', horiz=TRUE, legend=ages.lab, pch=15, col=plot.colors, cex=0.8 )
    axis(1, at = seq(5,n.yrs,5)-0.5, labels = years_full[seq(5,n.yrs,5)])
    box()
	}
  par(origpar)
} # end function
#plot.NAA(ssm)
#------------------------------------
#__annual estimates of Recruitment deviations (2 panel plot)____
plot.recruitment.devs <- function(mod, age.recruit = 1, units = expression(10^3), save.plots = FALSE)
{
  origpar <- par(no.readonly = TRUE)
	par(mfrow=c(2,1), mar=c(1,5,1,1), oma = c(4,1,1,0))
  dat = mod$env$data
	years <- mod$years
  nyrs = length(years)
	R <- mod$rep$NAA[,1]
	R.pred <- mod$rep$pred_NAA[,1]
	R.resids <- (log(R)-log(R.pred))/exp(mod$parList$log_NAA_sigma)[dat$NAA_sigma_pointers[1]]

	plot(years, R.pred, type='l', col='#114466', lty=1, lwd=2, xlab="",
	ylab= as.expression(substitute(paste("Age-", age.recruit, " Recruits (", units, ")", sep = ''),
		list(age.recruit = age.recruit, units = units[[1]]))), xlim = range(years), ylim=c(0, 1.1*max(R, R.pred)), axes = FALSE)
	lines(years, R, col="grey35", lwd=2, pch = 19, type = 'b')
	axis(1, labels = FALSE)
	axis(2)
	box()

	plot(years, R.resids, type='h', col='black', xlab="", ylab="Standardized (Log) Residuals", lwd=2,
	  xlim = range(years), ylim=1.1*range(R.resids))
	abline(h=0, lwd=1)
	mtext(side = 1, outer = TRUE, "Year", line = 2)
	par(origpar)
} # end function
#plot.recruitment.devs(ssm)

#scatter plot of SSB, R with 2-digit year as symbol (lag by 1 year)
plot.recr.ssb.yr <- function(mod, ssb.units = "kmt", recruits.units = expression(10^6), alpha = 0.05,
  scale.ssb = 1000, scale.recruits = 1000, age.recruit = 1, plot.colors, loglog=FALSE)
{
  origpar <- par(no.readonly = TRUE)
  par(mfrow=c(1,1), mar = c(4,5,1,1), oma = c(1,1,1,1))

  std <- summary(mod$sdrep, "report")
  cov <- mod$sdrep$cov
  dat = mod$env$data
  years = mod$years
  nyrs <- length(years)
  nages <- dat$n_ages
  #nstates <- model$dimensions$n_states
  #nprojyrs <- model$dimensions$n_proj_years
	ssb.ind <- which(rownames(std) == "log_SSB")
	log.ssb <- std[ssb.ind,1]
	ssb.cv <- std[ssb.ind,2]
  log.ssb.ci <- log.ssb + cbind(qnorm(1-alpha/2)*ssb.cv, -qnorm(1-alpha/2)*ssb.cv)
  R.ind = which(rownames(std) == "log_NAA_rep")[1:nyrs]
  log.R = std[R.ind,1]
  R.cv = std[R.ind,2]
	ssb.R.cov <- cov[c(ssb.ind,R.ind),c(ssb.ind,R.ind)]

  if(missing(plot.colors)) plot.colors = mypalette(nyrs-age.recruit)
	SR <- matrix(NA, (nyrs-age.recruit), 3)
	SR[,1] <- years[1:(nyrs-age.recruit)]
	SR[,2] <- exp(log.ssb[1:(nyrs-age.recruit)])/scale.ssb
	SR[,3] <- exp(log.R[age.recruit +1:(nyrs-age.recruit)])/scale.recruits
	yr.text <- substr(years[1:(nyrs-age.recruit)],3,4)
	npts <- nyrs-age.recruit
	ci.regs <- lapply(1:(nyrs-age.recruit), function(x)
	{
	  tcov <- cov[c(ssb.ind[x],R.ind[x+age.recruit]),c(ssb.ind[x],R.ind[x+age.recruit])]
	  tsd <- std[c(ssb.ind[x],R.ind[x+age.recruit]),2]
    tcor = tcov/(tsd %*% t(tsd))
	  el <- exp(ellipse::ellipse(tcor, level = 1-alpha, scale = tsd, centre = c(log(SR[x,2]),log(SR[x,3]))))
	  return(el)
	})
  # if(missing(max.y)) max.y <- max(sapply(ci.regs, function(x) max(x[,2])))
  # if(missing(max.x)) max.x <- max(sapply(ci.regs, function(x) max(x[,1])))
  max.y <- max(sapply(ci.regs, function(x) max(x[,2])))
  max.x <- max(sapply(ci.regs, function(x) max(x[,1])))
  na.lims <- any(is.na(c(max.y,max.x))) # check for na lims

	if(!loglog){ # untransformed SSB and Rec
    if(!na.lims){
    	plot(SR[,2], SR[,3], type='n', col='black',
    		xlab=as.expression(substitute(paste("SSB (", ssb.units, ")", sep = ""), list(ssb.units = ssb.units[[1]]))),
    		ylab= as.expression(substitute(paste("Age-", age.recruit, " Recruits (", units, ")", sep = ''),
    			list(age.recruit = age.recruit[[1]], units = recruits.units[[1]]))), ylim=c(0, max.y), xlim=c(0,max.x), axes = FALSE)
    	grid(col = gray(0.7), lty = 2)
    	axis(1)
    	axis(2)
    	box()
    	lines(SR[,2], SR[,3], col = gray(0.7), lwd =2)
    	for(i in 1:length(ci.regs))
      {
        tcol <- col2rgb(plot.colors[i])
        poly <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
        polygon(ci.regs[[i]][,1],ci.regs[[i]][,2], border = poly)
      }
      points(SR[npts,2], SR[npts,3], pch=19, col="#ffaa22", cex=2.5)
      text(SR[,2], SR[,3], yr.text, cex=0.9, col=plot.colors)
    } else { # NA lims
      max.y <- max(SR[,3])
      max.x <- max(SR[,2])
      plot(SR[,2], SR[,3], type='n', col='black',
        xlab=as.expression(substitute(paste("SSB (", ssb.units, ")", sep = ""), list(ssb.units = ssb.units[[1]]))),
        ylab= as.expression(substitute(paste("Age-", age.recruit, " Recruits (", units, ")", sep = ''),
          list(age.recruit = age.recruit[[1]], units = recruits.units[[1]]))), ylim=c(0, max.y), xlim=c(0,max.x), axes = FALSE)
      grid(col = gray(0.7), lty = 2)
      axis(1)
      axis(2)
      box()
      lines(SR[,2], SR[,3], col = gray(0.7), lwd =2)
      points(SR[npts,2], SR[npts,3], pch=19, col="#ffaa22", cex=2.5)
      text(SR[,2], SR[,3], yr.text, cex=0.9, col=plot.colors)
    }
	}

	if(loglog){ # log(SSB) and log(Rec)
    if(!na.lims){
    	plot(log(SR[,2]), log(SR[,3]), type='n', col='black',
    		xlab=as.expression(substitute(paste("Log-SSB (", ssb.units, ")", sep = ""), list(ssb.units = ssb.units[[1]]))),
    		ylab= as.expression(substitute(paste("Age-", age.recruit, " Log-Recruits (", units, ")", sep = ''),
    			list(age.recruit = age.recruit[[1]], units = recruits.units[[1]]))), ylim=c(log(min(SR[,3])), log(max.y)), xlim=c(log(min(SR[,2])),log(max.x)), axes = FALSE)
    	grid(col = gray(0.7), lty = 2)
    	axis(1)
    	axis(2)
    	box()
    	lines(log(SR[,2]), log(SR[,3]), col = gray(0.7), lwd =2)
    	for(i in 1:length(ci.regs))
      {
        tcol <- col2rgb(plot.colors[i])
        poly <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
        polygon(log(ci.regs[[i]][,1]),log(ci.regs[[i]][,2]), border = poly)
      }
      points(log(SR[npts,2]), log(SR[npts,3]), pch=19, col="#ffaa22", cex=2.5)
      text(log(SR[,2]), log(SR[,3]), yr.text, cex=0.9, col=plot.colors)
    } else { # na lims but correct max.x and max.y already calculated for untransformed plot SSB-Rec
      max.y <- max(SR[,3])
      max.x <- max(SR[,2])
      plot(log(SR[,2]), log(SR[,3]), type='n', col='black',
        xlab=as.expression(substitute(paste("Log-SSB (", ssb.units, ")", sep = ""), list(ssb.units = ssb.units[[1]]))),
        ylab= as.expression(substitute(paste("Age-", age.recruit, " Log-Recruits (", units, ")", sep = ''),
          list(age.recruit = age.recruit[[1]], units = recruits.units[[1]]))), ylim=c(log(min(SR[,3])), log(max.y)), xlim=c(log(min(SR[,2])),log(max.x)), axes = FALSE)
      grid(col = gray(0.7), lty = 2)
      axis(1)
      axis(2)
      box()
      lines(log(SR[,2]), log(SR[,3]), col = gray(0.7), lwd =2)
      points(log(SR[npts,2]), log(SR[npts,3]), pch=19, col="#ffaa22", cex=2.5)
      text(log(SR[,2]), log(SR[,3]), yr.text, cex=0.9, col=plot.colors)
    }
	}
  par(origpar)
}  #end function

#------------------------------------
plot.SARC.R.SSB <- function(mod, scale.ssb=1, scale.recruits=1, age.recruit = 1, ssb.units = 'mt', recruits.units = expression(10^3))
{
  origpar <- par(no.readonly = TRUE)
  par(mar = c(5,5,1,5), oma = c(0,0,0,1), family='serif')
  years = mod$years
  years_full = mod$years_full
  nyrs <- length(years_full)
  if(class(mod$sdrep)[1] == "sdreport"){
    std = summary(mod$sdrep)
  } else {
    std = mod$sdrep
  }
	ssb.ind <- which(rownames(std) == "log_SSB")
	log.ssb <- std[ssb.ind,1]
	R <- mod$rep$NAA[,1]
	ssb.plot <- exp(log.ssb[1:(nyrs-age.recruit)])/scale.ssb
	recr.plot <- R[age.recruit + 1:(nyrs-age.recruit)]/scale.recruits
	yr.text <- substr(years_full,3,4)

	max.r <- max(recr.plot)
	max.ssb <- max(ssb.plot)
	scale.r <- max(ssb.plot)/max(recr.plot)
	ylimr <- c(0,1.1*max(recr.plot))
	barplot(recr.plot/scale.recruits, axisnames=FALSE, width=1, space=rep(0,nyrs-age.recruit), offset=rep(-0.5,nyrs-age.recruit), axes=FALSE, xpd=FALSE,
		xlab = '', ylab ='', ylim = ylimr, xlim=c(0.5,nyrs-age.recruit - 0.5), col="lightcyan2")
	xr <-pretty(c(0,recr.plot/scale.recruits))
	axis(2, at = xr, lab = xr )
	axis(side=1, las=2, at=seq(0.5,nyrs-age.recruit-0.5, by=2),
	labels=as.character(seq(years_full[1],years_full[nyrs-age.recruit], by=2)), cex=0.75, las=2)

	y.ssb <- (ssb.plot)*max.r/max.ssb
	lines(seq(0.5,nyrs-age.recruit-0.5, by=1), y.ssb, lwd=2, col = 'navyblue')
	x <- pretty(c(0,ssb.plot))
	axis(4, at = c(0,x*max.r/max.ssb), lab = c(0,x), col='navyblue', col.axis="navyblue")
	box()
	mtext(side = 1, 'Year', line = 3)
	mtext(side = 4, as.expression(substitute(paste("SSB (", ssb.units, ")", sep = ""), list(ssb.units = ssb.units[[1]]))), line = 3, col='navyblue')
	mtext(side = 2, as.expression(substitute(paste("Age-", age.recruit, " Recruits (", units, ")", sep = ''),
		list(age.recruit = age.recruit[[1]], units = recruits.units[[1]]))), line = 3)
  if(length(years_full) > length(years)) abline(v=length(years)-age.recruit, lty=2, lwd=2)
	par(origpar)
}  # end function
#plot.SARC.R.SSB(ssm, ssm.aux)

plot.fleet.F <- function(mod, plot.colors)
{
  origpar <- par(no.readonly = TRUE)
  par(mfrow=c(1,1))
  years = mod$years
	n_fleets <- mod$env$data$n_fleets
  if(missing(plot.colors)) plot.colors = mypalette(n_fleets)
	if(n_fleets == 1){ # can calculate full F from FAA_tot in projection years if only one fleet
    years_full = mod$years_full
    nyrs <- length(years_full)
    if(class(mod$sdrep)[1] == "sdreport"){
      std = summary(mod$sdrep)
    } else {
      std = mod$sdrep
    }
    faa.ind <- which(rownames(std) == "log_FAA_tot")
    log.faa <- matrix(std[faa.ind,1], nyrs, mod$env$data$n_ages)
    full.f.ind <- cbind(1:nyrs, apply(log.faa,1, function(x) max(which(x == max(x)))))
    log.full.f <- log.faa[full.f.ind]
    plot(years_full, exp(log.full.f), xlab="Year", ylab="Full F", ylim=c(0,max(exp(log.full.f))), type='l', lty=1, lwd=2, col = plot.colors[1])
    if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
    grid(col = gray(0.7))
  } else { # multiple fleets, can't partition F in projection years to fleets
    for(i in 1:n_fleets){
  		if(i==1){
  		  plot(years, mod$rep$F[,i], xlab="Year", ylab="Full F", ylim=c(0,max(mod$rep$F)),	type='l', lty=1, lwd=2, col = plot.colors[i])
  		  grid(col = gray(0.7))
  		} else {
  		  lines(years, mod$rep$F[,i],lty=i, lwd=2, col=plot.colors[i])
      }
  	}
  }
	leg.names <- paste0("Fleet ",1:n_fleets)
	legend('topleft', legend=leg.names, col=plot.colors,lwd=rep(2, n_fleets), lty=seq(1, n_fleets), horiz=TRUE, bty='n')
	par(origpar)
}   # end function
#plot.fleet.F(ssm,ssm.aux)

#------------------------------------
plot.cv <- function(mod)
{
  origpar <- par(no.readonly = TRUE)
	par(mfrow=c(1,1), mar=c(4,4,2,2))
  years = mod$years
  years_full = mod$years_full
  nyrs <- length(years_full)
	if(class(mod$sdrep)[1] == "sdreport"){
    std = summary(mod$sdrep)
  } else {
    std = mod$sdrep
  }
	ssb.ind <- which(rownames(std) == "log_SSB")
	log.ssb <- std[ssb.ind,1]
	ssb.cv <- std[ssb.ind,2]
	R.ind <- which(rownames(std) == "log_NAA_rep")[1:nyrs]
	log.R <- std[R.ind,1]
  R.cv <- std[R.ind,2]
	faa.ind <- which(rownames(std) == "log_FAA_tot")
	log.faa <- matrix(std[faa.ind,1], nyrs, mod$env$data$n_ages)
	faa.cv <- matrix(std[faa.ind,2], nyrs, mod$env$data$n_ages)

	full.f.ind <- cbind(1:nyrs, apply(log.faa,1, function(x) max(which(x == max(x)))))
  log.full.f <- log.faa[full.f.ind]
  full.f.cv <- faa.cv[full.f.ind]

  any.na <- any(is.na(c(R.cv, faa.cv, full.f.cv)))

  if(!any.na){
  	plot(years_full, R.cv, type='l', lwd=2, col='black', xlab="Year", ylab="CV", ylim=c(0, 1.1*max(R.cv, ssb.cv, full.f.cv)))
  	lines(years_full, ssb.cv, lwd=2, col="blue")
  	lines(years_full, full.f.cv, lwd=2, col="green3")
    if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
  	legend('top', legend=c("Recruits", "SSB", "Full F"), col=c("black", "blue", "green3"), lty=rep(1,3), lwd=rep(2,3), horiz=T)
  }
  par(origpar)
}  # end function
#------------------------------------

#------------------------------------
plot.M <- function(mod, ages, ages.lab, alpha = 0.05, plot.colors)
{
  dat = mod$env$data
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = mod$ages.lab
	n_ages <- length(ages)
	meanMAA <- apply(mod$rep$MAA,2,mean)
  years = mod$years
  years_full = mod$years_full
  n_years = length(years_full)
	if(missing(plot.colors)) plot.colors = mypalette(n_years)
	n.M.by.age <- lapply(1:n_ages, function(x) table(mod$rep$MAA[,x]))
	plot(ages,meanMAA,lwd=2,xlab="Age",ylab="Natural Mortality Rate", ylim=c(0,1.1*(max(mod$rep$MAA))), type= 'n', xlim = c(min(ages)-0.5,max(ages)+1),
		axes=FALSE)
	grid(col = gray(0.7), lwd = 2)
	axis(1, at = ages, labels = ages.lab, lwd = 2)
	axis(2, lwd = 2)
	box(lwd = 2)
	sapply(1:n_ages, function(x)
	{
		y <- sort(unique(mod$rep$MAA[,x]))
		# ind1 <- sapply(y, function(z) which(mod$rep$MAA[,x] == y)[1])
		segments(x-0.2, y, x+0.2, y, lwd = 2,col = plot.colors[1:length(y)])
		text(x+0.2, y, paste('n =', table(mod$rep$MAA[,x])), pos = 4)
	})
}

#------------------------------------
#--------Data Plots------------------
plot.catch.by.fleet <- function(mod, units = "mt", plot.colors)
{
  origpar <- par(no.readonly = TRUE)
  par(mfrow = c(1,1))
  dat = mod$env$data
  years = mod$years
  nyrs = length(years)
	catch.obs <- dat$agg_catch
	n_fleets <- dat$n_fleets
  if(missing(plot.colors)) plot.colors = mypalette(n_fleets)
	barplot(t(catch.obs), xlab="Year", ylab= paste0("Catch (", units, ")"), ylim=c(0,1.1*max(apply(catch.obs,1,sum))), col=plot.colors,space=0)
	axis(side=1, at = seq(2,nyrs,2)-0.5, labels = years[seq(2,nyrs,2)], cex=0.75)
	box(lwd = 2)
	if (n_fleets > 1)
  {
    legend('top', legend=paste0("Fleet ",1:n_fleets), horiz=TRUE, pch=15, col=plot.colors)

    # do proportions only if n_fleets > 1
		catch.prop <- catch.obs/apply(catch.obs,1,sum)
		barplot(t(catch.prop), xlab="Year", ylab="Proportion of Catch", ylim=c(0,1.1), col=plot.colors, space=0)
    axis(side=1, las=2, at = seq(2,nyrs,2)-0.5, labels = years[seq(2,nyrs,2)], cex=0.75, las=2)
    box(lwd = 2)
		legend('top', legend=paste0("Fleet ",1:n_fleets), horiz=TRUE, pch=15, col=plot.colors)
	}
	par(origpar)
}

# Bubble plots of catch age comps (set is.catch.flag to False to plot Discard age comps)
plot.catch.age.comp.bubbles <- function(mod, ages, ages.lab, bubble.col = "#8c8c8caa", i=1)
{
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  years = mod$years
  nyrs = length(years)
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = mod$ages.lab
  n_ages = length(ages)
	par(mar=c(4,4,2,2), oma=c(1,1,1,1), mfrow=c(1,1))
  n_fleets = dat$n_fleets
	# for (i in 1:n_fleets)
	# {
		acomp.obs <- dat$catch_paa[i,,]
		catch.yrs <- which(dat$use_catch_paa[,i] == 1)
		my.title <- "Age Comps for Catch for Fleet "
		if (length(catch.yrs)>0)
		{
			scale.catch.obs <- 5
			z3 <- as.matrix(acomp.obs) * scale.catch.obs

			plot(ages, rev(ages),  xlim = range(ages), ylim = c(years[nyrs],(years[1]-2)), xlab = "Age", ylab = "", type = "n", axes=FALSE)
			axis(1, at= ages, lab = ages.lab)
			axis(2, at = rev(years), lab = rev(years), cex.axis=0.75, las=1)
			box()
			abline(h=years, col="lightgray")
			segments(x0=ages, y0=rep(years[1],n_ages), x1=ages, y1=rep(years[nyrs],n_ages), col = "lightgray", lty = 1)
			for (j in 1:nyrs) points(ages, rep(years[j], n_ages), cex=z3[j,], col="black", bg = bubble.col, pch = 21)

			bubble.legend1 <- c(0.05,0.2,0.4)
			bubble.legend2 <- bubble.legend1 * scale.catch.obs
			legend("topright", xpd=TRUE, legend=bubble.legend1, pch=rep(21, 3), pt.cex=bubble.legend2, horiz=T , col='black', pt.bg = bubble.col)
			title (paste0(my.title,i), outer=TRUE, line=-1)
		} # end catch.yrs test
	# }   #end loop n_fleets
  par(origpar)
}

#------------------------------------
plot.index.input <- function(mod, plot.colors)
{
  origpar <- par(no.readonly = TRUE)
  par(mfrow=c(2,1), mar = c(1,1,1,1), oma = c(4,4,2,0))
  years = mod$years
  nyrs = length(years)
  dat = mod$env$data
	indvals <- dat$agg_indices
	indvals[which(dat$use_indices!=1)] <- NA
	n_indices = dat$n_indices
	# rescale to mean 1 and stdev 1
	rescaled <- indvals
	my.mean <- apply(indvals,2,mean, na.rm=TRUE)
	my.std <- apply(indvals,2,sd, na.rm=TRUE)
	for (i in 1:n_indices) rescaled[,i] <- (indvals[,i] - my.mean[i]) / my.std[i]

	my.range <- range(rescaled, na.rm=TRUE)
  if(missing(plot.colors)) plot.colors = mypalette(n_indices)
	plot(years,rescaled[,1],xlab="",ylab="", axes = FALSE, xlim = range(years), ylim=my.range,col=plot.colors[1],type='n')
	grid(col = gray(0.7))
	axis(1, labels = FALSE)
	axis(2)
	box()
	mtext(side = 2, "Rescaled Indices", outer = FALSE, line = 3)
	for (i in 1:n_indices) lines(years,rescaled[,i],col=plot.colors[i])
  legend("top", legend = paste0("Index " , 1:n_indices), col = plot.colors, lty = 1, horiz = TRUE, xpd = NA, inset = c(0,-0.1), bty = "n")

	# now repeat on log scale
	log.indvals <- log(indvals)
	log.rescaled <- log.indvals
	my.log.mean <- apply(log.indvals,2,mean, na.rm=T)
	my.log.std <- apply(log.indvals,2,sd, na.rm=T)
	for (i in 1:n_indices) log.rescaled[,i] <- (log.indvals[,i] - my.log.mean[i]) / my.log.std[i]

	my.log.range <- range(log.rescaled, na.rm=T)
	plot(years,log.rescaled[,1], xlab="",ylab="",ylim=my.log.range,col=plot.colors[1],type='n')
	grid(col = gray(0.7))
	for (i in 1:n_indices) lines(years,log.rescaled[,i],col=plot.colors[i])
	mtext(side = 2, "Rescaled Log(Indices)", outer = FALSE, line = 3)
	mtext(side = 1, "Year", outer = FALSE, line = 3)
	par(origpar)
}
#plot.index.input(ssm)

#------------------------------------
# Bubble plots of index age comps
plot.index.age.comp.bubbles <- function(mod, ages, ages.lab, bubble.col = "#8c8c8caa", i=1)
{
  origpar <- par(no.readonly = TRUE)
  par(mar=c(4,4,2,2), oma=c(1,1,1,1), mfrow=c(1,1))
  years = mod$years
  nyrs = length(years)
  dat = mod$env$data
  n_indices = dat$n_indices
  if(missing(ages)) ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = mod$ages.lab
	n_ages <- length(ages)

	# for (i in 1:n_indices)
	# {
		acomp.obs <- dat$index_paa[i,,]
		index.yrs <- which(dat$use_index_paa[,i] == 1)
		my.title <- "Age Comps for Index "
		if (length(index.yrs)>0)
		{
			scale.index.obs <- 5
			z3 <- as.matrix(acomp.obs) * scale.index.obs

			plot(ages, rev(ages),  xlim = range(ages), ylim = c(years[nyrs],(years[1]-2)),
			xlab = "Age", ylab = "", type = "n", axes=FALSE)
			axis(1, at= ages, lab=ages.lab)
			axis(2, at = rev(years), lab = rev(years), cex.axis=0.75, las=1)
			box()
			abline(h=years, col="lightgray")
			segments(x0=ages, y0=rep(years[1],n_ages), x1=ages, y1=rep(years[nyrs],n_ages), col = "lightgray", lty = 1)
			for (j in 1:nyrs) points(ages, rep(years[j], n_ages), cex=z3[j,], col="black", bg = bubble.col, pch = 21)

			bubble.legend1 <- c(0.05,0.2,0.4)
			bubble.legend2 <- bubble.legend1 * scale.index.obs
			legend("topright", xpd=TRUE, legend=bubble.legend1, pch=rep(21, 3), pt.cex=bubble.legend2, horiz=TRUE, col='black', pt.bg = bubble.col)
			title (paste0(my.title,i), outer=T, line=-1)
		} # end index.yrs test
	# }   #end loop n_fleets
	par(origpar)
}

#------------------------------------
plot.waa <- function(mod,type="ssb",plot.colors,ind=1)
{
  origpar <- par(no.readonly = TRUE)
  years = mod$years_full
  nyrs = length(years)
  dat = mod$env$data
  if(missing(plot.colors)) plot.colors = mypalette(dat$n_ages)
  point = switch(type,
    ssb = dat$waa_pointer_ssb,
    jan1 = dat$waa_pointer_jan1,
    fleets = dat$waa_pointer_fleets[ind],
    indices = dat$waa_pointer_indices[ind],
    totcatch = dat$waa_pointer_totcatch
  )
  labs = switch(type,
    ssb = "SSB",
    jan1 = "January 1 Biomass",
    fleets = paste0("Fleet ", ind),
    indices = paste0("Index ", ind),
    # fleets = paste0("Fleet ", 1:dat$n_fleets),
    # indices = paste0("Index ", 1:dat$n_indices),
    totcatch = "Total Catch"
  )
  waa = dat$waa[point,,]
  n = ifelse(length(dim(waa)) == 2, 1, dim(waa)[1])
	for(i in 1:n)
	{
		if(n>1) WAA.plot <- dat$waa[i,,]
    else WAA.plot = waa
		plot(years,years,xlab="Year",ylab="Weight",ylim=c(0,max(WAA.plot)),type='n')
		for (a in 1:dat$n_ages)
		{
			lines(years,WAA.plot[,a],col=plot.colors[a],lwd=2)
			lines(years,rep(mean(WAA.plot[,a]),length(years)),lty=2,col=plot.colors[a])
		}
		title(main = paste0("Annual Weight-at-Age for ", labs[i]))
	}  # end k-loop
	par(origpar)
}  # end function

#------------------------------------
plot.maturity <- function(mod, ages.lab, plot.colors)
{
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  years = mod$years
  n_years = length(years)
  ages = 1:dat$n_ages
  if(missing(ages.lab)) ages.lab = mod$ages.lab
	meanmaturity <- apply(dat$mature,2,mean)
	if(missing(plot.colors)) plot.colors <- mypalette(n_years)

	plot(ages,meanmaturity,type='l',lwd=2,xlab="Age",ylab="Maturity",ylim=c(0,max(dat$mature)), axes = FALSE)
  axis(1, at = ages, labels = ages.lab, lwd = 2)
  axis(2, lwd = 2)
  box(lwd = 2)
	if (length(unique(dat$mature)) > length(ages))
	{
		for (i in 1:n_years) points(jitter(ages, factor=0.4), dat$mature[i,],col=plot.colors[i])
    midi <- floor(n_years/2)
		legend('topleft', horiz=FALSE, legend=c(years[1],years[midi],years[n_years]), pch=c(1,1,1), col=c(plot.colors[1], plot.colors[midi], plot.colors[n_years]))
	}
	title(main="Maturity", outer=FALSE)
	par(origpar)
}

#-------SSB/R -----------------------------
get_SPR = function(F, M, sel, mat, waassb, fracyrssb, at.age = FALSE)
{
  n_ages = length(sel)
  SPR = numeric()
  n = 1
  F = F * sel
  Z = F + M
  for(a in 1:(n_ages-1))
  {
    SPR[a] = n[a] * mat[a] * waassb[a] * exp(-fracyrssb * Z[a])
    n[a+1] = n[a] * exp(-Z[a])
  }
  n[n_ages] = n[n_ages]/(1-exp(-Z[n_ages]))
  SPR[n_ages] = n[n_ages] * mat[n_ages] * waassb[n_ages] * exp(-fracyrssb * Z[n_ages])
  if(at.age) return(SPR)
  else return(sum(SPR))
}

#-------Y/R -----------------------------
get_YPR = function(F, M, sel, waacatch, at.age = FALSE)
{
  n_ages = length(sel)
  YPR = numeric()
  n = 1
  F = F * sel
  Z = F + M
  for(a in 1:(n_ages-1))
  {
    YPR[a] = n[a] * F[a] * waacatch[a] * (1.0 - exp(-Z[a]))/Z[a]
    n[a+1] = n[a] * exp(-Z[a])
  }
  n[n_ages] = n[n_ages]/(1 - exp(-Z[n_ages]))
  YPR[n_ages] = n[n_ages] * F[n_ages] * waacatch[n_ages] * (1.0 - exp(-Z[n_ages]))/Z[n_ages]
  if(at.age) return(YPR)
  else return(sum(YPR))
}

plot.SPR.table <- function(mod, nyrs.ave = 5, plot=TRUE)
{
  origpar <- par(no.readonly = TRUE)
  spr.targ.values <- seq(0.2, 0.8, 0.05)
	n.spr <- length(spr.targ.values)
  dat = mod$env$data
	n_ages<- dat$n_ages
	years <- mod$years
	n_years <- length(years)
  avg.ind = (n_years-nyrs.ave+1):n_years
	#fec.age <- apply(dat$waa[dat$waa_pointer_ssb,,][avg.ind,],2,mean)
	mat.age <- apply(dat$mature[avg.ind,],2,mean)
	ssb.waa <- apply(dat$waa[dat$waa_pointer_ssb,,][avg.ind,],2,mean)
	catch.waa <- apply(dat$waa[dat$waa_pointer_totcatch,,][avg.ind,],2,mean)
	M.age <- apply(mod$rep$MAA[avg.ind,],2,mean)
  sel = apply(apply(mod$rep$FAA[avg.ind,,,drop=FALSE], 2:3, mean),2,sum) #sum across fleet averages 
  #sel = apply(mod$rep$FAA_tot[avg.ind,],2,mean) #average FAA, then do selectivity
	sel <- sel/max(sel)
	spawn.time <- mean(dat$fracyr_SSB[avg.ind])
  spr0 = get_SPR(F=0, M=M.age, sel=sel, mat=mat.age, waassb=ssb.waa, fracyrssb = spawn.time)
	F.start <- 0.11  # starting guess for optimization routine to find F_SPR%

	f.spr.vals <- rep(NA, n.spr)
	ypr.spr.vals <- rep(NA, n.spr)
	conv.vals <- rep(NA, n.spr)

	for (i in 1:n.spr)
	{
		t.spr <- spr.targ.values[i]

		spr.f <- function(F.start)
		{
      spr = get_SPR(F=F.start, M=M.age, sel=sel, mat=mat.age, waassb=ssb.waa, fracyrssb = spawn.time)
			abs(spr/spr0 - t.spr)
		}
		yyy <- nlminb(start=F.start, objective=spr.f, lower=0, upper=3)
		f.spr.vals[i] <- yyy$par
    ypr.spr.vals[i] = get_YPR(F = f.spr.vals[i], M=M.age, sel = sel, waacatch= catch.waa)
	}  #end i-loop over SPR values

	spr.target.table<- as.data.frame(cbind(spr.targ.values, f.spr.vals, ypr.spr.vals))
	colnames(spr.target.table) <- c("%SPR", "F(%SPR)", "YPR")
	par(mfrow=c(1,1), mar=c(4,4,2,4))

	if(plot){ # plot, not table
  	plot(spr.targ.values, ypr.spr.vals, type='n', xlab="% SPR Target", ylab="", lwd=2, col="blue3",
      ylim=c(0,1.2*max(ypr.spr.vals)), axes = FALSE)
    box(lwd = 2)
    axis(1, lwd = 2)
  	abline(v=seq(0.2,0.8, by=0.05), col="grey85")
  	lines(spr.targ.values, ypr.spr.vals, lwd=2, col="blue3" )
  	points(spr.targ.values, ypr.spr.vals, pch=19, col="blue3" )

    scale.f.spr <- max(f.spr.vals)/max(ypr.spr.vals)
  	axis(side=2, #at=seq(0,1,by=0.1)/scale.f.spr, lab=seq(0,1,by=0.1),
      las=2, col='blue3', col.axis="blue3", lwd = 2)
    mtext(side = 2, "Yield per Recruit", line = 3, col = "blue3")


  	lines(spr.targ.values, f.spr.vals/scale.f.spr, lwd = 2, col="red")
  	points(spr.targ.values, f.spr.vals/scale.f.spr, pch = 19, col="red")
  	axis(side=4, at=seq(0,1,by=0.1)/scale.f.spr, lab=seq(0,1,by=0.1), las=2, col='red', col.axis="red", lwd = 2)
  	mtext(side=4, "F (%SPR)", line=3, col="red")

  	title (paste("SPR Target Reference Points (Years Avg = ", nyrs.ave,")", sep=""), outer=T, line=-1)
	}

	if(!plot){ # table, not plot
  	par(mfrow=c(1,1), mar=c(2,2,2,2))
  	plot(seq(1,15), seq(1,15), type='n', axes=F, bty='n',xlab="",ylab="")

  	text(x=2,y=14, labels="% SPR", font=2, pos=4)
  	text(x=5, y=14, labels="F(%SPR)" , font=2, pos=4)
  	text(x=9, y=14, labels="YPR", font=2, pos=4 )
  	for (i in 1:n.spr)
  	{
  		text(x=2, y=seq(n.spr,1, by=-1), labels=round(spr.targ.values,2), cex=1.0, pos=4, font=1)
  		text(x=5, y=seq(n.spr,1, by=-1), labels=round(f.spr.vals,4), cex=1.0, pos=4, font=1)
  		text(x=9, y=seq(n.spr,1, by=-1), labels=round(ypr.spr.vals,4), cex=1.0, pos=4, font=1)
  	}
  	title (paste("SPR Target Reference Points (Years Avg = ", nyrs.ave,")", sep=""), outer=TRUE, line=-1)
	}
  par(origpar)
	#write.csv( spr.target.table, file=paste(od,"SPR_Target_Table.csv", sep=""),  row.names=F)
} # end function

#------------------------------------
plot.annual.SPR.targets <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, plot.colors, od)
{
  origpar <- par(no.readonly = TRUE)
  spr.targ.values <- seq(0.2, 0.5, by=0.1)
	n.spr <- length(spr.targ.values)
  dat = mod$env$data
  n_ages = dat$n_ages
  years = mod$years
  n_years = length(years)
  years_full = mod$years_full
  n_years_full = length(years_full)

  fec.age <- dat$waa[dat$waa_pointer_ssb,,]
	mat.age <- dat$mature
	wgt.age <- dat$waa[dat$waa_pointer_totcatch,,]
	M.age <- mod$rep$MAA
	sel.age <- mod$rep$FAA_tot/apply(mod$rep$FAA_tot,1,max)

	spawn.time <- dat$fracyr_SSB #vector length n_years

	spr0 <- numeric()
	F.start <- 0.11  # starting guess for optimization routine to find F_SPR%

	f.spr <- matrix(NA, n_years_full, n.spr)
	ypr.spr <- matrix(NA, n_years_full, n.spr)
	conv <- matrix(NA, n_years_full, n.spr)

	for (j in 1:n_years_full)
	{
    spr0[j] = get_SPR(F=0, M=M.age[j,], sel=sel.age[j,], mat=mat.age[j,], waassb=fec.age[j,], fracyrssb = spawn.time[j])
    #spr0.vals[j] <- s.per.recr(n_ages=n_ages, fec.age=fec.age[j,], mat.age=mat.age[j,], M.age= M.age[j,], F.mult=0, sel.age=sel.age[j,],
    #  spawn.time=spawn.time)
	  for (i in 1:n.spr)
	  {
			t.spr <- spr.targ.values[i]

			spr.f <- function(F.start)
			{
        spr = get_SPR(F=F.start, M=M.age[j,], sel=sel.age[j,], mat=mat.age[j,], waassb=fec.age[j,], fracyrssb = spawn.time[j])
        abs(spr/spr0[j] - t.spr)
				#abs(s.per.recr(n_ages=n_ages, fec.age=fec.age[j,], mat.age=mat.age[j,], M.age= M.age[j,], F.mult=F.start, sel.age=sel.age[j,],
				#	spawn.time=spawn.time)/spr0.vals[j] - t.spr)
			}
			yyy <- nlminb(start=F.start, objective=spr.f, lower=0, upper=3)
			f.spr[j,i] <- yyy$par
			ypr.spr[j,i] = get_YPR(F = f.spr[j,i], M=M.age[j,], sel = sel.age[j,], waacatch= wgt.age[j,])
      #ypr(n_ages, wgt.age=wgt.age[j,], M.age=M.age[j,],  F.mult=f.spr.vals[j,i], sel.age=sel.age[j,])
		}  # end j-loop over n_years
	}  #end i-loop over SPR values

  plot.colors = mypalette(n.spr)
  if(do.tex) cairo_pdf(file.path(od, paste0("FSPR_annual_time.pdf")), family = fontfam, height = 10, width = 10)
  if(do.png) png(filename = file.path(od, paste0("FSPR_annual_time.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
	par(mfrow=c(2,1), mar=c(4,5,2,1))
	plot(years_full, f.spr[,1], type='n', xlab="Years", ylab=expression(paste("Full ",italic(F)["%SPR"])), lwd=2, ylim=c(0,1.2*max(f.spr)))
	for (i in 1:n.spr) lines(years_full, f.spr[,i], lwd=2, col=plot.colors[i])
  if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
	legend('top', legend= sapply(spr.targ.values, function(x) as.expression(bquote(italic(F)[paste(.(x*100),"%")]))), col=plot.colors, horiz=TRUE, lwd=2, cex=0.9)
	title (main=expression(paste("Annual ", italic(F)["%SPR"], " Reference Points")), line=1)
	#if (save.plots) savePlot(paste(od, "Annual_FSPR.", plotf, sep=''), type=plotf)
	plot(years_full, ypr.spr[,1], type='n', xlab="Years", ylab=expression(paste("YPR(",italic(F)["%SPR"],")")), lwd=2, ylim=c(0,1.2*max(ypr.spr)))
	for (i in 1:n.spr) lines(years_full, ypr.spr[,i], lwd=2, col=plot.colors[i])
	if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
  #legend('top', legend=c("YPR20%", "YPR30%", "YPR40%", "YPR50%"), col=plot.colors, horiz=TRUE, lwd=2, cex=0.9)
	title (main=expression(paste("Annual YPR(",italic(F)["%SPR"], ") Reference Points")), line=1)
	if(do.tex | do.png) dev.off() else par(origpar)

	if(do.tex) cairo_pdf(file.path(od, paste0("FSPR_freq_annual_F.pdf")), family = fontfam, height = 10, width = 10)
	if(do.png) png(filename = file.path(od, paste0("FSPR_freq_annual_F.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
	par(mfrow=c(2,2))
  f.hist = lapply(1:length(spr.targ.values), function(x) hist(f.spr[1:n_years,x], plot = FALSE))
  par(mfrow=c(2,2), mar=c(4,1,2,2), oma=c(2,3,2,0))
  for(i in 1:length(spr.targ.values))
  {
    plot(f.hist[[i]]$mids, f.hist[[i]]$counts, xlab=bquote("Full " ~ italic(F)[paste(.(spr.targ.values[i]*100),"%")]), ylab="", type='h', lwd=2, ylim=c(0, max(f.hist[[i]]$counts)), col=plot.colors[i])
    lines(f.hist[[i]]$mids, f.hist[[i]]$counts, lwd=2, col=plot.colors[i])
  }
  mtext(side=2, outer = TRUE, "Frequency", line = 1)
	title (main=expression(paste("Frequencies of Annual ", italic(F)["%SPR"], " Reference Points")), outer=TRUE, line=0, cex.main = 2)
	if(do.tex | do.png) dev.off() else par(origpar)

	if(do.tex) cairo_pdf(file.path(od, paste0("FSPR_freq_annual_YPR.pdf")), family = fontfam, height = 10, width = 10)
	if(do.png) png(filename = file.path(od, paste0("FSPR_freq_annual_YPR.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
	ypr.hist = lapply(1:length(spr.targ.values), function(x) hist(ypr.spr[1:n_years,x], plot = FALSE))
  par(mfrow=c(2,2), mar=c(4,1,2,2), oma=c(2,3,2,0))
  for(i in 1:length(spr.targ.values))
  {
    plot(ypr.hist[[i]]$mids, ypr.hist[[i]]$counts, xlab=bquote(paste("YPR(", italic(F)[paste(.(spr.targ.values[i]*100),"%")],")")), ylab="", type='h', lwd=2, ylim=c(0, max(ypr.hist[[i]]$counts)), col=plot.colors[i])
    lines(ypr.hist[[i]]$mids, ypr.hist[[i]]$counts, lwd=2, col=plot.colors[i])
  }
  mtext(side=2, outer = TRUE, "Frequency", line = 1)
	title (main=expression(paste("Frequencies of Annual YPR(",italic(F)["%SPR"], ") Reference Points")), outer = TRUE, line=0, cex.main = 2)
	if(do.tex | do.png) dev.off() else par(origpar)

	# par(origpar)
} # end function

plot.SR.pred.line <- function(mod, ssb.units = "mt", SR.par.year, recruits.units = "thousands", scale.ssb = 1,
	scale.recruits = 1, age.recruit = 1, plot.colors)
{
  if(mod$env$data$recruit_model == 3 & mod$env$data$use_steepness != 1) #B-H stock recruit function with alpha/beta
  {
    if(class(mod$sdrep)[1] == "sdreport"){
      std = summary(mod$sdrep)
    } else {
      std = mod$sdrep
    }
    ssb.ind = which(rownames(std) == "log_SSB")
    years = mod$years
    nyrs = length(years)
    log.ssb <- std[ssb.ind,1]
    R <- mod$rep$NAA[,1]
    SR <- matrix(NA, (nyrs-age.recruit), 3)
    SR[,1] <- years[1:(nyrs-age.recruit)]/scale.ssb
    SR[,2] <- exp(log.ssb[1:(nyrs-age.recruit)])
    SR[,3] <- R[age.recruit +1:(nyrs-age.recruit)]/scale.recruits
    log_a <- mod$parList$mean_rec_pars[1]
    log_b <- mod$parList$mean_rec_pars[2]
    if(missing(SR.par.year)) SR.par.year = nyrs
    a.b.ind = which(rownames(std) == "mean_rec_pars")
    l.ab <- TMB:::as.list.sdreport(mod$sdrep, what = "Estimate")$mean_rec_pars
    #l.ab = std[a.b.ind,1]
    a.b.cov = mod$sdrep$cov.fixed[a.b.ind,a.b.ind]
    if(length(a.b.ind) == 1 & is.null(mod$input$map$mean_rec_pars)) {
      warning("The number of estimated stock-recruit parameters is less than 2, but input$map$mean_rec_pars is not provided.")
      return()
    }
    not.na.ind <- 1:2
    if(!is.null(mod$input$map$mean_rec_pars)) not.na.ind <- which(!is.na(mod$input$map$mean_rec_pars))
    if(length(not.na.ind)<2) {
      a.b.cov <- matrix(0, 2,2)
      if(length(not.na.ind)==1) a.b.cov[not.na.ind,not.na.ind] <- mod$sdrep$cov.fixed[a.b.ind,a.b.ind]
    }
    lR.fn = function(la, lb, S) la  + log(S) - log(1 + exp(lb)*S)
    dlR.dp = Deriv::Deriv(lR.fn, x = c("la","lb"))
    seq.ssb <- seq(0, max(SR[,2]), length.out=300)
    sd.pred.lR = sapply(seq.ssb, function(x)
    {
      d = dlR.dp(l.ab[1], l.ab[2], S= x)
      return(c(sqrt(t(d) %*% a.b.cov %*% d)))
    })
    pred.lR <- lR.fn(l.ab[1], l.ab[2], seq.ssb)
    #exp(log_a) *seq.ssb/(1 + exp(log_b)*seq.ssb)
    ci.pred.lR = pred.lR + qnorm(0.975)*cbind(-sd.pred.lR, sd.pred.lR) - log(scale.recruits)

    if(missing(plot.colors)) plot.colors = mypalette(nyrs-age.recruit)
    tcol <- col2rgb(plot.colors[1])
    tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')

    plot(SR[,2], SR[,3], type='p', col='black', pch=19,
      xlab=as.expression(substitute(paste("SSB (", ssb.units, ")", sep = ""), list(ssb.units = ssb.units[[1]]))),
      ylab= as.expression(substitute(paste("Age-", age.recruit, " Recruits (", units, ")", sep = ''),
        list(age.recruit = age.recruit[[1]], units = recruits.units[[1]]))), ylim=c(0, max(SR[,3])), xlim=c(0,1.1*max(SR[,2])))
    lines(seq.ssb, exp(pred.lR), col=plot.colors[1], lwd=2)
    polygon(c(seq.ssb,rev(seq.ssb)), exp(c(ci.pred.lR[,1],rev(ci.pred.lR[,2]))), col = tcol, border = "transparent")
  }
}

plot.FXSPR.annual <- function(mod, alpha = 0.05, status.years, max.x, max.y, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
  origpar <- par(no.readonly = TRUE)
  percentSPR = mod$env$data$percentSPR
	n_ages = mod$env$data$n_ages
  years_full = mod$years_full
	n_years_full = length(years_full)
  years = mod$years
  n_years = length(years)
  if(missing(status.years)){
    status.years = n_years
    status.lwd = 1
    if(mod$env$data$do_proj==1){
      status.years <- c(status.years, n_years_full)
      status.lwd <- c(2,1)
    }
  } else {
    status.lwd <- rep(1, length(status.years))
  }
  tcol <- col2rgb('black')
  tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
  std <- summary(mod$sdrep, "report")
	inds <- list(Y.t = which(rownames(std) == "log_Y_FXSPR"))
	inds$F.t <- which(rownames(std) == "log_FXSPR")
	inds$SSB.t <- which(rownames(std) == "log_SSB_FXSPR")
	inds$ssb <- which(rownames(std) == "log_SSB")
	inds$faa <- which(rownames(std) == "log_FAA_tot")
  log.faa <- matrix(std[inds$faa,1], n_years_full, n_ages)
  age.full.f <- apply(log.faa,1, function(x) max(which(x == max(x))))
  inds$full.f <- (age.full.f-1)*n_years_full + 1:n_years_full  + min(inds$faa) - 1 #cbind(1:n_years, age.full.f)
  na.sd <- sapply(inds, function(x) any(is.na(std[x,2])))
  ylabs <- c(
    bquote(paste('Yield(',italic(F)[paste(.(percentSPR), "%")], ')')),
    bquote(italic(F)[paste(.(percentSPR), "%")]),
    bquote(paste('SSB(', italic(F)[paste(.(percentSPR), "%")],')')))
	cov <- mod$sdrep$cov
  log.rel.ssb.rel.F.cov <- lapply(1:n_years_full, function(x)
  {
    K <- cbind(c(1,-1,0,0),c(0,0,1,-1))
    ind <- c(inds$ssb[x],inds$SSB.t[x],inds$full.f[x],inds$F.t[x])
    tcov <- cov[ind,ind]
    return(t(K) %*% tcov %*% K)
  })

  # FSPR absolute --------------------------------------------------
  if(do.tex) cairo_pdf(file.path(od, paste0("FSPR_absolute.pdf")), family = fontfam, height = 10, width = 10)
  if(do.png) png(filename = file.path(od, paste0("FSPR_absolute.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
  par(mfrow = c(3,1), mar = c(2,5,1,1), oma = c(4,2,1,1))
  for(i in 1:3)
  {
    t.ind <- inds[[i]]
    t.ylab <- ylabs[i]
    vals <- std[t.ind,1][1:n_years_full]
    cv <- std[t.ind,2][1:n_years_full]
    ci <-  vals + cbind(qnorm(1-alpha/2)*cv, -qnorm(1-alpha/2)*cv)
    if(all(!is.nan(unique(vals)))){
  	  if(!na.sd[i]) plot(years_full, exp(vals), xlab = '', ylab = t.ylab, ylim = c(0,max(exp(ci),na.rm= TRUE)), type = 'l', cex.lab = 2)
      if(na.sd[i]) plot(years_full, exp(vals), xlab = '', ylab = t.ylab, ylim = c(0,max(exp(vals),na.rm= TRUE)), type = 'l', cex.lab = 2)
  	  grid(col = gray(0.7))
      polyy = exp(c(ci[,1],rev(ci[,2])))
      polyx = c(years_full,rev(years_full))
      polyx = polyx[!is.na(polyy)]
      polyy = polyy[!is.na(polyy)]
      polygon(polyx, polyy, col = tcol, border = tcol, lwd = 1)
      if(mod$env$data$do_proj==1) abline(v=tail(years,1), lty=2, lwd=1)
      mtext(side = 1, outer = TRUE, "Year", cex = 2, line = 2)
    } else { # all nan, print error message
      plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
      text(x = 0.5, y = 0.55, paste("Error in plot, all values are NaN"), cex = 1.6, col = "black")
      text(x = 0.5, y = 0.45, ylabs[[i]], cex = 1.6, col = "black")
    }
	}
  if(do.tex | do.png) dev.off() else par(origpar)

  # FSPR relative --------------------------------------------------
  if(do.tex) cairo_pdf(file.path(od, paste0("FSPR_relative.pdf")), family = fontfam, height = 10, width = 10)
  if(do.png) png(filename = file.path(od, paste0("FSPR_relative.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
  par(mfrow=c(2,1))
  log.rel.ssb.vals <- std[inds$ssb,1][1:n_years_full] - std[inds$SSB.t,1][1:n_years_full]
  cv <- sapply(log.rel.ssb.rel.F.cov, function(x) return(sqrt(x[1,1])))
  ci <-  log.rel.ssb.vals + cbind(-qnorm(1-alpha/2)*cv, qnorm(1-alpha/2)*cv)
  if(all(!is.nan(unique(log.rel.ssb.vals)))){
    plot(years_full, exp(log.rel.ssb.vals), xlab = '', ylab = bquote(paste("SSB/", SSB[paste(.(percentSPR),"%")])), ylim = c(0,5), type = 'l')
    grid(col = gray(0.7))
    polyy = exp(c(ci[,1],rev(ci[,2])))
    polyy[which(is.infinite(polyy))] <-1e10
    polyx = c(years_full,rev(years_full))
    polyx = polyx[!is.na(polyy)]
    polyy = polyy[!is.na(polyy)]
    polygon(polyx, polyy, col = tcol, border = tcol, lwd = 1)
    abline(h=1, lty = 2)
    abline(h=0.5, lty = 2, col = 'red')
    if(mod$env$data$do_proj==1) abline(v=tail(mod$years,1), lty=2, lwd=1)
  } else {
    plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
    text(x = 0.5, y = 0.58, paste("Error in plot, all values are NaN"), cex = 1.6, col = "black")
    text(x = 0.5, y = 0.42, bquote(paste("SSB/", SSB[paste(.(percentSPR),"%")])), cex = 1.6, col = "black")    
  }

  log.rel.f.vals <- std[inds$full.f,1][1:n_years_full] - std[inds$F.t,1][1:n_years_full]
  cv <- sapply(log.rel.ssb.rel.F.cov, function(x) return(sqrt(x[2,2])))
  ci <-  log.rel.f.vals + cbind(-qnorm(1-alpha/2)*cv, qnorm(1-alpha/2)*cv)
  if(all(!is.nan(unique(log.rel.f.vals)))){
    polyy = exp(c(ci[,1],rev(ci[,2])))
    polyy[which(is.infinite(polyy))] <-1e10
    if(!na.sd["full.f"]) plot(years_full, exp(log.rel.f.vals), xlab = '', ylab = bquote(paste(italic(F),"/", italic(F)[paste(.(percentSPR),"%")])),
      ylim = c(0,max(polyy,1, na.rm = TRUE)), type = 'l')
    if(na.sd["full.f"]) plot(years_full, exp(log.rel.f.vals), xlab = '', ylab = bquote(paste(italic(F),"/", italic(F)[paste(.(percentSPR),"%")])),
      ylim = c(0,max(exp(log.rel.f.vals),1, na.rm = TRUE)), type = 'l')
    grid(col = gray(0.7))
    polyx = c(years_full,rev(years_full))
    polyx = polyx[!is.na(polyy)]
    polyy = polyy[!is.na(polyy)]
    polygon(polyx, polyy, col = tcol, border = tcol, lwd = 1)
    abline(h=1, lty = 2, col = 'red')
    if(mod$env$data$do_proj==1) abline(v=tail(years,1), lty=2, lwd=1)
    mtext(side =1, "Year", outer = TRUE, line = 2, cex = 1.5)
  } else {
    plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
    text(x = 0.5, y = 0.58, paste("Error in plot, all values are NaN"), cex = 1.6, col = "black")
    text(x = 0.5, y = 0.42, bquote(paste(italic(F),"/", italic(F)[paste(.(percentSPR),"%")])), cex = 1.6, col = "black")      
  }
  if(do.tex | do.png) dev.off() else par(origpar)

  # Kobe plot - only if sdreport was successful ---------------------------
  do.kobe <- sapply(log.rel.ssb.rel.F.cov[status.years], function(x) !all(!is.finite(x))) # only if some non-infinite values for at least some status years
  # if(!mod$na_sdrep){
  if(any(do.kobe)){  
    rel.ssb.rel.F.ci.regs <- lapply(status.years, function(x){
      if(is.na(log.rel.f.vals[x]) | any(diag(log.rel.ssb.rel.F.cov[[x]])<0)) return(matrix(NA,100,2))
      else return(exp(ellipse::ellipse(log.rel.ssb.rel.F.cov[[x]], centre = c(log.rel.ssb.vals[x],log.rel.f.vals[x]), level = 1-alpha)))
      })
    p.ssb.lo.f.lo <- p.ssb.lo.f.hi <- p.ssb.hi.f.lo <- p.ssb.hi.f.hi <- rep(NA,length(status.years))
    for(i in 1:length(status.years)){
      check.zero.sd <- diag(log.rel.ssb.rel.F.cov[[status.years[i]]])==0 | diag(log.rel.ssb.rel.F.cov[[status.years[i]]]) < 0
      if(!any(is.na(check.zero.sd))) if(!any(check.zero.sd)){
        p.ssb.lo.f.lo[i] <- mnormt::sadmvn(lower = c(-Inf,-Inf), upper = c(-log(2), 0), mean = c(log.rel.ssb.vals[status.years[i]],log.rel.f.vals[status.years[i]]), varcov = log.rel.ssb.rel.F.cov[[status.years[i]]])
        p.ssb.lo.f.hi[i] <- mnormt::sadmvn(lower = c(-Inf,0), upper = c(-log(2), Inf), mean = c(log.rel.ssb.vals[status.years[i]],log.rel.f.vals[status.years[i]]), varcov = log.rel.ssb.rel.F.cov[[status.years[i]]])
        p.ssb.hi.f.lo[i] <- mnormt::sadmvn(lower = c(-log(2),-Inf), upper = c(Inf, 0), mean = c(log.rel.ssb.vals[status.years[i]],log.rel.f.vals[status.years[i]]), varcov = log.rel.ssb.rel.F.cov[[status.years[i]]])
        p.ssb.hi.f.hi[i] <- mnormt::sadmvn(lower = c(-log(2),0), upper = c(Inf, Inf), mean = c(log.rel.ssb.vals[status.years[i]],log.rel.f.vals[status.years[i]]), varcov = log.rel.ssb.rel.F.cov[[status.years[i]]])
      }
    }

    vals <- exp(cbind(log.rel.ssb.vals, log.rel.f.vals))
    if(missing(max.x)) max.x <- max(sapply(rel.ssb.rel.F.ci.regs, function(x) max(x[,1],na.rm = TRUE)),1.5)
    if(missing(max.y)) max.y <- max(sapply(rel.ssb.rel.F.ci.regs, function(x) max(x[,2],na.rm = TRUE)),1.5)
    if(is.infinite(max.y)) max.y <- 10
    if(is.infinite(max.x)) max.x <- 10
    if(do.tex) cairo_pdf(file.path(od, paste0("Kobe_status.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("Kobe_status.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mfrow = c(1,1))
    plot(vals[status.years,1],vals[status.years,2], ylim = c(0,max.y), xlim = c(0,max.x), xlab = bquote(paste("SSB/", SSB[paste(.(percentSPR),"%")])),
      ylab = bquote(paste(italic(F),"/", italic(F)[paste(.(percentSPR),"%")])),type = 'n')
    lims = par("usr")
    tcol <- col2rgb('red')
    tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
    polygon(c(lims[1],0.5,0.5,lims[1]),c(1,1,lims[4],lims[4]), border = tcol, col = tcol)
    tcol <- col2rgb('green')
    tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
    polygon(c(0.5,lims[2],lims[2],0.5),c(lims[3],lims[3],1,1), border = tcol, col = tcol)
    tcol <- col2rgb('yellow')
    tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
    polygon(c(lims[1],0.5,0.5,lims[1]),c(lims[3],lims[3],1,1), border = tcol, col = tcol)
    polygon(c(0.5,lims[2],lims[2],0.5),c(1,1,lims[4],lims[4]), border = tcol, col = tcol)
    legend("topleft", legend = paste0("Prob = ", round(p.ssb.lo.f.hi,2)), bty = "n", text.font=status.lwd)
    legend("topright", legend = paste0("Prob = ", round(p.ssb.hi.f.hi,2)), bty = "n", text.font=status.lwd)
    legend("bottomleft", legend = paste0("Prob = ", round(p.ssb.lo.f.lo,2)), bty = "n", text.font=status.lwd)
    legend("bottomright", legend = paste0("Prob = ", round(p.ssb.hi.f.lo,2)), bty = "n", text.font=status.lwd)
    text(vals[status.years,1],vals[status.years,2], substr(years_full[status.years],3,4), font=status.lwd)
    for(i in 1:length(status.years)) if(all(!is.infinite(rel.ssb.rel.F.ci.regs[[i]]))) { 
      polygon(rel.ssb.rel.F.ci.regs[[i]][,1], rel.ssb.rel.F.ci.regs[[i]][,2], lwd=status.lwd[i])#, border = gray(0.7))
    }
    if(do.tex | do.png) dev.off() else par(origpar)
    return(list(p.ssb.lo.f.lo = p.ssb.lo.f.lo, p.ssb.hi.f.lo = p.ssb.hi.f.lo, p.ssb.hi.f.hi = p.ssb.hi.f.hi, p.ssb.lo.f.hi = p.ssb.lo.f.hi))
  } else { return(NULL) }
}  # end function

plot.MSY.annual <- function(mod, alpha = 0.05, max.x, max.y, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  n_ages = dat$n_ages
  years = mod$years
  n_years = length(years)
  years_full = mod$years_full
  n_years_full = length(years_full)
  std = summary(mod$sdrep, "report")
  cov <- mod$sdrep$cov
	if(dat$recruit_model == 3) #Beverton-Holt assumed in model fit
	{ # test to make sure steepness was estimated
    tcol <- col2rgb('black')
    tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
		inds <- list(MSY = which(rownames(std) == "log_MSY"))
		inds$FMSY <- which(rownames(std) == "log_FMSY")
		inds$SSBMSY <- which(rownames(std) == "log_SSB_MSY")
		inds$RMSY <- which(rownames(std) == "log_R_MSY")
		inds$ssb <- which(rownames(std) == "log_SSB")
  	inds$faa <- which(rownames(std) == "log_FAA_tot")
	  log.faa <- matrix(std[inds$faa,1], n_years_full, n_ages)
	  age.full.f <- apply(log.faa,1, function(x) max(which(x == max(x))))
	  inds$full.f <- (age.full.f-1)*n_years_full + 1:n_years_full  + min(inds$faa) - 1 #cbind(1:n_years, age.full.f)
	  ylabs <- c(expression(MSY),expression(italic(F)[MSY]), expression(SSB[MSY]), expression(italic(R)[MSY]))
	  log.rel.ssb.rel.F.cov <- lapply(1:n_years_full, function(x)
	  {
	    K <- cbind(c(1,-1,0,0),c(0,0,1,-1))
	    ind <- c(inds$ssb[x],inds$SSBMSY[x],inds$full.f[x],inds$FMSY[x])
	    tcov <- cov[ind,ind]
	    return(t(K) %*% tcov %*% K)
	  })

    # 4-panel MSY plot
    if(do.tex) cairo_pdf(file.path(od, paste0("MSY_4panel_F_SSB_R.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("MSY_4panel_F_SSB_R.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mfrow=c(2,2))
    for(i in 1:4)
    {
      t.ind <- inds[[i]]
      t.ylab <- ylabs[i]
      vals <- std[t.ind,1][1:n_years_full]
  	  cv <- std[t.ind,2][1:n_years_full]
      ci <-  vals + cbind(qnorm(1-alpha/2)*cv, -qnorm(1-alpha/2)*cv)
      na.ci <- all(is.na(ci))
      # remove NaN and Inf
      # rm.rows <- which(vals < 0)
      vals[!is.finite(vals)] = NA
      # vals[rm.rows] = NA
      # rm.rows <- which(ci[,2] < 0)
      ci[!is.finite(exp(ci))] = NA
      # ci[rm.rows,] = NA
		  if(!na.ci) plot(years_full, exp(vals), xlab = 'Year', ylab = t.ylab, ylim = c(0,max(exp(ci),na.rm=TRUE)), type = 'l')
      if(na.ci) plot(years_full, exp(vals), xlab = 'Year', ylab = t.ylab, ylim = c(0,max(exp(vals),na.rm=TRUE)), type = 'l')
		  grid(col = gray(0.7))
      not.na.ci <- !is.na(ci[,1])
		  polygon(c(years_full[not.na.ci],rev(years_full[not.na.ci])), exp(c(ci[,1][not.na.ci],rev(ci[,2][not.na.ci]))), col = tcol, border = tcol, lwd = 1)
      # polygon(c(years_full,rev(years_full)), exp(c(ci[,1],rev(ci[,2]))), col = tcol, border = tcol, lwd = 1)
      if(mod$env$data$do_proj==1) abline(v=tail(years,1), lty=2, lwd=1)
		}
    if(do.tex | do.png) dev.off() else par(origpar)

    # 2-panel SSB_MSY and F_MSY
    if(do.tex) cairo_pdf(file.path(od, paste0("MSY_relative_status.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("MSY_relative_status.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mfrow=c(2,1))
    rel.ssb.vals <- exp(std[inds$ssb,1][1:n_years_full] - std[inds$SSBMSY,1][1:n_years_full])
    cv <- sapply(log.rel.ssb.rel.F.cov, function(x) return(sqrt(x[1,1])))
    ci <-  rel.ssb.vals * exp(cbind(-qnorm(1-alpha/2)*cv, qnorm(1-alpha/2)*cv))
    na.ci <- all(is.na(ci))
      # remove NaN and Inf
      # rm.rows <- which(rel.ssb.vals < 0)
      rel.ssb.vals[!is.finite(rel.ssb.vals)] = NA
      # rel.ssb.vals[rm.rows] = NA
      # rm.rows <- which(ci[,2] < 0 | ci[,1] < 0)
      ci[!is.finite(ci)] = NA
      # ci[rm.rows,] = NA    
		if(!na.ci) plot(years_full, rel.ssb.vals, xlab = 'Year', ylab = expression(paste("SSB/", SSB[MSY])), ylim = c(0,max(ci,na.rm=TRUE)), type = 'l')
    if(na.ci) plot(years_full, rel.ssb.vals, xlab = 'Year', ylab = expression(paste("SSB/", SSB[MSY])), ylim = c(0,max(rel.ssb.vals,na.rm=TRUE)), type = 'l')
	  grid(col = gray(0.7))
	  polygon(c(years_full,rev(years_full)), exp(c(ci[,1],rev(ci[,2]))), col = tcol, border = "transparent", lwd = 1)
	  abline(h=1, lty = 2)
	  abline(h=0.5, lty = 2, col = 'red')
    if(mod$env$data$do_proj==1) abline(v=tail(years,1), lty=2, lwd=1)

    rel.f.vals <- exp(std[inds$full.f,1][1:n_years_full] - std[inds$FMSY,1][1:n_years_full])
    cv <- sapply(log.rel.ssb.rel.F.cov, function(x) return(sqrt(x[2,2])))
    ci <-  rel.f.vals * exp(cbind(-qnorm(1-alpha/2)*cv, qnorm(1-alpha/2)*cv))
    na.ci <- all(is.na(ci))
      # remove NaN and Inf
      # rm.rows <- which(rel.f.vals < 0)
      rel.f.vals[!is.finite(rel.f.vals)] = NA
      # rel.f.vals[rm.rows] = NA
      # rm.rows <- which(ci[,2] < 0 | ci[,1] < 0)
      ci[!is.finite(ci)] = NA
      # ci[rm.rows,] = NA    
		if(!na.ci) plot(years_full, rel.f.vals, xlab = 'Year', ylab = expression(paste(italic(F),"/", italic(F)[MSY])),
      ylim = c(0,max(ci,na.rm=TRUE)), type = 'l')
    if(na.ci) plot(years_full, rel.f.vals, xlab = 'Year', ylab = expression(paste(italic(F),"/", italic(F)[MSY])),
      ylim = c(0,max(rel.f.vals,na.rm=TRUE)), type = 'l')
	  grid(col = gray(0.7))
	  polygon(c(years_full,rev(years_full)), c(ci[,1],rev(ci[,2])), col = tcol, border = tcol, lwd = 1)
	  abline(h=1, lty = 2, col = 'red')
    if(mod$env$data$do_proj==1) abline(v=tail(years,1), lty=2, lwd=1)
    if(do.tex | do.png) dev.off() else par(origpar)
	}
}  # end function

#------------------------------------
plot.yield.curves <- function(mod, nyrs.ave = 5, plot=TRUE, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
	n_ages = dat$n_ages
  years = mod$years
  n_years = length(years)
  avg.ind = (n_years-nyrs.ave+1):n_years
	mat.age <- apply(dat$mature[avg.ind,],2,mean)
	ssb.waa <- apply(dat$waa[dat$waa_pointer_ssb,,][avg.ind,],2,mean)
	catch.waa <- apply(dat$waa[dat$waa_pointer_totcatch,,][avg.ind,],2,mean)
	M.age <- apply(mod$rep$MAA[avg.ind,],2,mean)
  sel = apply(apply(mod$rep$FAA[avg.ind,,,drop=FALSE], 2:3, mean),2,sum) #sum across fleet averages 
  #sel = apply(mod$rep$FAA_tot[avg.ind,],2,mean) #average FAA, then do selectivity
	sel <- sel/max(sel)
	spawn.time <- mean(dat$fracyr_SSB[avg.ind])
  spr0 = get_SPR(F=0, M=M.age, sel=sel, mat=mat.age, waassb=ssb.waa, fracyrssb = spawn.time)

	F.range <- seq(0,2.0, by=0.01)
	nF <- length(F.range)

  spr = sapply(F.range, function(x) get_SPR(F=x, M=M.age, sel=sel, mat=mat.age, waassb=ssb.waa, fracyrssb = spawn.time))
  ypr = sapply(F.range, function(x) get_YPR(F=x, M=M.age, sel=sel, waacatch=catch.waa))
  pr = spr/spr0

  if(plot){
    if(do.tex) cairo_pdf(file.path(od, paste0("YPR_F_curve_plot.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("YPR_F_curve_plot.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
  	par(mfrow=c(1,1), mar=c(4,4,2,4))
  	plot(F.range, ypr, type='n', xlab="Full F", ylab="", lwd=2, col="blue3", ylim=c(0,max(ypr)), axes = FALSE)
    box()
    axis(1)
    axis(2, col='blue3', col.axis="blue3")
  	abline(v=seq(0.1,2.0, by=0.1), col="grey85")
  	lines(F.range, ypr, lwd=2, col="blue3" )
  	points(F.range, ypr, pch=19, col="blue3")
    mtext(side =2, "Yield per Recruit", col = "blue3", line = 2)
  	scale.pr <- max(pr)/max(ypr)

  	lines(F.range, pr/scale.pr, col="red", lwd=2)
  	points(F.range, pr/scale.pr, pch=19, col="red")
  	axis(side=4, at=seq(0,1,by=0.1)/scale.pr, lab=seq(0,1,by=0.1), las=2, col='red', col.axis="red")
  	mtext(side=4, "% SPR", line=3, col="red")
  	title (paste("YPR-SPR Reference Points (Years Avg = ", nyrs.ave,")", sep=""), outer=T, line=-1)
  	if(do.tex | do.png) dev.off() else par(origpar)
  }

  if(!plot){
    if(do.tex) cairo_pdf(file.path(od, paste0("YPR_F_curve_table.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("YPR_F_curve_table.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)

  	par(mfrow=c(1,1), mar=c(2,2,2,2))
  	plot(seq(1,42), seq(-3,38), type='n', axes=F, bty='n',xlab="",ylab="")
  	text(x=0,y=36, labels="F", font=2, pos=4)
  	text(x=4, y=36, labels="YPR" , font=2, pos=4)
  	text(x=9, y=36, labels="SPR", font=2, pos=4)
  	text(x=15,y=36, labels="F", font=2, pos=4)
  	text(x=19, y=36, labels="YPR" , font=2, pos=4)
  	text(x=24, y=36, labels="SPR", font=2, pos=4)
  	text(x=30,y=36, labels="F", font=2, pos=4)
  	text(x=34, y=36, labels="YPR" , font=2, pos=4)
  	text(x=39, y=36, labels="SPR", font=2, pos=4)

    text(x=0, y=35:1, labels=F.range[1:35], cex=0.82, pos=4, font=1)
    text(x=4, y=35:1, labels=round(ypr[1:35],4), cex=0.82, pos=4, font=1)
    text(x=9, y=35:1, labels=round(pr[1:35],4), cex=0.82, pos=4, font=1)

    text(x=15, y=35:1, labels=F.range[36:70], cex=0.82, pos=4, font=1)
    text(x=19, y=35:1, labels=round(ypr[36:70],4), cex=0.82, pos=4, font=1)
    text(x=24, y=35:1, labels=round(pr[36:70],4), cex=0.82, pos=4, font=1)

    text(x=30, y=35:1, labels=F.range[71:105], cex=0.82, pos=4, font=1)
    text(x=34, y=35:1, labels=round(ypr[71:105],4), cex=0.82, pos=4, font=1)
    text(x=39, y=35:1, labels=round(pr[71:105],4), cex=0.82, pos=4, font=1)

  	title (paste("YPR-SPR Reference Points (Years Avg = ", nyrs.ave,")", sep=""), outer=T, line=-1)
  	if(do.tex | do.png) dev.off() else par(origpar)
  }
  ypr.table = cbind.data.frame(F.range, ypr, pr)
  colnames(ypr.table) <- c("Full.F", "YPR", "SPR")
	#write.csv(ypr.table, file=paste(od,"YPR_Table.csv", sep=""), row.names=F)
  # par(origpar)
} # end function

#------------------------------------
#------------------------------------
plot.exp.spawn <- function(mod, nyrs.ave = 5)
{
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  n_ages = dat$n_ages
  years = mod$years
  n_years = length(years)
  avg.ind = (n_years-nyrs.ave+1):n_years
	mat.age <- apply(dat$mature[avg.ind,],2,mean)
	ssb.waa <- apply(dat$waa[dat$waa_pointer_ssb,,][avg.ind,],2,mean)
	catch.waa <- apply(dat$waa[dat$waa_pointer_totcatch,,][avg.ind,],2,mean)
	M.age <- apply(mod$rep$MAA[avg.ind,],2,mean)
  sel = apply(apply(mod$rep$FAA[avg.ind,,,drop=FALSE], 2:3, mean),2,sum) #sum across fleet averages 
  #sel = apply(mod$rep$FAA_tot[avg.ind,],2,mean) #average FAA, then do selectivity
	sel <- sel/max(sel)
	spawn.time <- mean(dat$fracyr_SSB[avg.ind])
  spr0 = get_SPR(F=0, M=M.age, sel=sel, mat=mat.age, waassb=ssb.waa, fracyrssb = spawn.time)

	F.range <- seq(0,2.0, by=0.01)
	nF <- length(F.range)
	# plot maturity vs selectivity based on "nyrs.ave" average

	par(mfrow=c(1,1), mar=c(4,4,2,4))
	plot( seq(1,n_ages), mat.age, type='b', pch = 16, col='black', lwd=2, ylim=c(0,1.1), xlab="Age", ylab="Selectivity or Maturity at age")
	lines( seq(1,n_ages), sel, col='orange2', lwd=1)
	points( seq(1,n_ages), sel, col='orange2', pch=16)
	legend('topleft', legend=c("Maturity", "Selectivity"), col=c("black", "orange2"), lwd=c(2,1), pch=c(NA, 16))
  spr = sapply(F.range, function(x) get_SPR(F=x, M=M.age, sel=sel, mat=mat.age, waassb=ssb.waa, fracyrssb = spawn.time))
  expected.spawners = sapply(F.range, function(x) get_SPR(F=x, M=M.age, sel=sel, mat=mat.age, waassb=rep(1,n_ages), fracyrssb = spawn.time))
  ypr = sapply(F.range, function(x) get_YPR(F=x, M=M.age, sel=sel, waacatch=catch.waa))
  pr = spr/spr0

	par(mfrow=c(1,1), mar=c(4,4,2,4))

	plot(F.range, expected.spawners, type='n', xlab="Full F", ylab="", lwd=2, col="skyblue3", ylim=c(0,max(expected.spawners)), axes = FALSE)
  box()
  axis(1)
  axis(2, col='skyblue3', col.axis="skyblue3")
  mtext(side = 2, col = "skyblue3", "Expected Spawnings", line = 2)
	abline(v=seq(0.1,2.0, by=0.1), col="grey85")
	abline(h=c(1,2,3), col="skyblue1")
	lines(F.range, expected.spawners, lwd=2, col="skyblue3")
	points(F.range, expected.spawners, pch=19, col="skyblue3")
	scale.pr <- max(pr)/max(expected.spawners)

	lines(F.range, pr/scale.pr, col="red", lwd=2)
	points(F.range, pr/scale.pr, col="red", pch=19)
	axis(side=4, at=seq(0,1,by=0.1)/scale.pr, lab=seq(0,1,by=0.1), las=2, col='red', col.axis="red")
	mtext(side=4, "% SPR", line=3, col="red")
	title (paste("Expected Spawnings and SPR Reference Points (Years Avg = ", nyrs.ave,")",	sep=""), cex=0.9, outer=T, line=-1)

	par(mfrow=c(1,1), mar=c(2,2,2,2))
	plot(seq(1,42), seq(-3,38), type='n', axes=F, bty='n',xlab="",ylab="")
	text(x=0,y=36, labels="F", font=2, pos=4)
	text(x=4, y=36, labels="E[Sp]" , font=2, pos=4,cex=0.9)
	text(x=9, y=36, labels="SPR", font=2, pos=4)
	text(x=15,y=36, labels="F", font=2, pos=4)
	text(x=19, y=36, labels="E[Sp]" , font=2, pos=4,cex=0.9)
	text(x=24, y=36, labels="SPR", font=2, pos=4)
	text(x=30,y=36, labels="F", font=2, pos=4)
	text(x=34, y=36, labels="E[Sp]" , font=2, pos=4,cex=0.9)
	text(x=39, y=36, labels="SPR", font=2, pos=4)

  text(x=0, y=seq(35,1, by=-1), labels=F.range[1:35], cex=0.82, pos=4, font=1)
  text(x=4, y=seq(35,1, by=-1), labels=round(expected.spawners[1:35],4), cex=0.82, pos=4, font=1)
  text(x=9, y=seq(35,1, by=-1), labels=round(pr[1:35],4), cex=0.82, pos=4, font=1)

  text(x=15, y=seq(35,1, by=-1), labels=F.range[36:70], cex=0.82, pos=4, font=1)
  text(x=19, y=seq(35,1, by=-1), labels=round(expected.spawners[36:70],4), cex=0.82, pos=4, font=1)
  text(x=24, y=seq(35,1, by=-1), labels=round(pr[36:70],4), cex=0.82, pos=4, font=1)

  text(x=30, y=seq(35,1, by=-1), labels=F.range[71:105], cex=0.82, pos=4, font=1)
  text(x=34, y=seq(35,1, by=-1), labels=round(expected.spawners[71:105],4), cex=0.82, pos=4, font=1)
  text(x=39, y=seq(35,1, by=-1), labels=round(pr[71:105],4), cex=0.82, pos=4, font=1)

	title (paste("Expected Spawnings & SPR Reference Points (Years Avg = ", nyrs.ave,")", sep=""), cex=0.9, outer=T, line=-1)


	exp.spawn.table<- as.data.frame(cbind(F.range, expected.spawners, pr))
	colnames(exp.spawn.table) <- c("Full.F", "Exp.Spawn", "SPR")
	#write.csv(exp.spawn.table, file=paste(od,"Exp.Spawn.Table.csv", sep=""), row.names=F)
  par(origpar)
} # end function

plot.retro <- function(mod,y.lab,y.range1,y.range2, alpha = 0.05, what = "SSB", age=NULL, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
  origpar <- par(no.readonly = TRUE)
  years = mod$years
	n_years <- length(years) # don't use projections
  npeels = length(mod$peels)
  if(npeels)
  {
    if(what == "NAA_age") {
      age.ind <- (1:mod$input$data$n_ages)[age] 
      #age.ind <- match(age, mod$ages.lab)
      what.print <- paste0(what, age)
    } else {
      what.print <- what
    }

    # standard retro plot
    if(do.tex) cairo_pdf(file.path(od, paste0(what.print,"_retro.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0(what.print,"_retro.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    plot.colors = mypalette(npeels+1)
    tcol = col2rgb(plot.colors)
    tcol = rgb(tcol[1,],tcol[2,],tcol[3,], maxColorValue = 255, alpha = 200)
    if(what == "NAA_age"){
      res = list(head(mod$rep[["NAA"]],n_years))
      res[2:(npeels+1)] = lapply(mod$peels, function(x) x$rep[["NAA"]])
    } else {
      res = list(head(mod$rep[[what]],n_years))
      res[2:(npeels+1)] = lapply(mod$peels, function(x) x$rep[[what]])
    }
    if(what == "NAA")
    {
      par(mfcol = c(ceiling(NCOL(res[[1]])/2),2))
      for(i in 1:NCOL(res[[1]]))
      {
        y.range1 <- range(sapply(res, function(x) range(x[,i])))
        plot(years,res[[1]][,i],lwd=1,col=plot.colors[1],type='l',xlab="Year",ylab=paste0("Numbers at age ", mod$ages.lab[i]),ylim=y.range1)
        grid(col = gray(0.7), lty = 2)
        for (j in 1:npeels)
        {
          lines(years[1:(n_years-j)],res[[j+1]][,i], col = tcol[j+1])
          points(years[n_years-j],res[[j+1]][n_years-j,i],pch=16,col=plot.colors[j+1])
        }
      }
    }
    if(what == "NAA_age") # only specified age
    {
      par(mfrow = c(1,1))
      i = age.ind
      y.range1 <- range(sapply(res, function(x) range(x[,i])))
      plot(years,res[[1]][,i],lwd=1,col=plot.colors[1],type='l',xlab="Year",ylab=paste0("Numbers at age ", mod$ages.lab[i]),ylim=y.range1)
      grid(col = gray(0.7), lty = 2)
      for (j in 1:npeels)
      {
        lines(years[1:(n_years-j)],res[[j+1]][,i], col = tcol[j+1])
        points(years[n_years-j],res[[j+1]][n_years-j,i],pch=16,col=plot.colors[j+1])
      }
    }
    if(what %in% c("SSB","Fbar"))
    {
      if(missing(y.range1)) y.range1 <- range(sapply(res, function(x) range(x)))
      par(mfrow = c(1,1))
      plot(years,res[[1]],lwd=1,col=plot.colors[1],type='l',xlab="Year",ylab=what,ylim=y.range1)
      grid(col = gray(0.7), lty = 2)
      for (i in 1:npeels)
      {
        lines(years[1:(n_years-i)],res[[i+1]], col = tcol[i+1])
        points(years[n_years-i],res[[i+1]][n_years-i],pch=16,col=plot.colors[i+1])
      }
    }
    if(do.tex | do.png) dev.off() else par(origpar)

    # relative retro plot
    if(do.tex) cairo_pdf(file.path(od, paste0(what.print,"_retro_relative.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0(what.print,"_retro_relative.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    if(missing(y.lab)) y.lab = bquote(paste("Mohn's ", rho, "(",.(what),")"))
    if(what %in% c("NAA","NAA_age")) rel.res = lapply(1:length(res), function(x) res[[x]]/res[[1]][1:(n_years - x + 1),] - 1)
    if(what %in% c("SSB","Fbar")) rel.res = lapply(1:length(res), function(x) res[[x]]/res[[1]][1:(n_years - x + 1)] - 1)
    rho.vals = mohns_rho(mod)

    if(what == "NAA")
    {
      par(mfcol = c(ceiling(NCOL(rel.res[[1]])/2),2))
      for(i in 1:NCOL(res[[1]]))
      {
        y.range2 <- c(-1,max(sapply(rel.res, function(x) range(x[,i]))))
        plot(years,rel.res[[1]][,i],lwd=1,col=plot.colors[1],type='l',xlab="Year",ylab=bquote(paste("Mohn's ", rho, "(Numbers at age ", .(mod$ages.lab[i]), ")")),ylim = y.range2)
        grid(col = gray(0.7), lty = 2)
        for (j in 1:npeels)
        {
          lines(years[1:(n_years-j)],rel.res[[j+1]][,i], col = tcol[j+1])
          points(years[n_years-j],rel.res[[j+1]][n_years-j,i],pch=16,col=plot.colors[j+1])
        }
        rho.nm = ifelse(i==1, "R",paste0("N",mod$ages.lab[i]))
        rho.plot <- round(rho.vals[rho.nm],3)
        legend("bottomleft", legend = bquote(rho == .(rho.plot)), bty = "n")
      }
    }
    if(what == "NAA_age")
    {
      par(mfrow = c(1,1))
      i = age.ind
      y.range2 <- c(-1,max(sapply(rel.res, function(x) range(x[,i]))))
      plot(years,rel.res[[1]][,i],lwd=1,col=plot.colors[1],type='l',xlab="Year",ylab=bquote(paste("Mohn's ", rho, "(Numbers at age ", .(mod$ages.lab[i]), ")")),ylim = y.range2)
      grid(col = gray(0.7), lty = 2)
      for (j in 1:npeels)
      {
        lines(years[1:(n_years-j)],rel.res[[j+1]][,i], col = tcol[j+1])
        points(years[n_years-j],rel.res[[j+1]][n_years-j,i],pch=16,col=plot.colors[j+1])
      }
      rho.nm = ifelse(i==1, "R",paste0("N",mod$ages.lab[i]))
      rho.plot <- round(rho.vals[rho.nm],3)
      legend("bottomleft", legend = bquote(rho == .(rho.plot)), bty = "n")
    }
    if(what %in% c("SSB","Fbar"))
    {
      if(missing(y.range2)) y.range2 <- c(-1,max(sapply(rel.res, function(x) range(x))))
      par(mfrow = c(1,1))
      plot(years,rel.res[[1]],lwd=1,col="black",type='l',xlab="Year",ylab=y.lab,ylim=y.range2)
      grid(col = gray(0.7), lty = 2)
      for (i in 1:npeels)
      {
        lines(years[1:(n_years-i)],rel.res[[i+1]], col = tcol[i+1])
        points(years[n_years-i],rel.res[[i+1]][n_years-i],pch=16,col=plot.colors[i+1])
      }
      rho.plot <- round(rho.vals[what],3)
      legend("bottomleft", legend = bquote(rho == .(rho.plot)), bty = "n")
    }
    if(do.tex | do.png) dev.off() else par(origpar)
    # par(origpar)
    return(rho.vals)
  }
}

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-----Catch curve and cohort correspondence plots-------------------------------
#-------------------------------------------------------------------------------
# function to compute catch at age (matrix)
# from proportions at age (matrix), weight at age (matrix) and total weight (vector)
wtprop2caa <- function(fleet)#totwt,waa,props)
{

  return(paa * totwt/apply(paa*waa,1,sum))
	caa <- mod$rep$pred_caa[,fleet,] * mod$env$data
	for (i in 1:length(totwt))
	{
		if (sum(props[i,]) == 0) caa[i,] <- 0
		if (sum(props[i,]) > 0) caa[i,] <- props[i,] * (totwt[i] / sum(props[i,] * waa[i,]))
	}
	return(caa)
}

#-------------------------------------------------------------------------------
# function replace zeros and take log
rep0log <- function(mat)
{
  mat[mat==0] = NA
	return(log(mat))
}
#-------------------------------------------------------------------------------

# function make cohorts
makecohorts <- function(mat)
{
	NAmat <- matrix(NA, NCOL(mat), NCOL(mat))
	matcoh1 <- rbind(NAmat,mat,NAmat)
	nr <- NROW(matcoh1)
	nc <- NCOL(mat)
	matcoh <- matrix(NA, nrow=nr, ncol=nc)
	for (i in 1:nc) matcoh[1:(nr-i),i] <- matcoh1[i:(nr-1),i]
	return(matcoh)
}
#-------------------------------------------------------------------------------

# function to plot all ages by all ages of data by cohorts
# assumes matrix has ages as columns and that cohorts are in rows
plotcoh <- function(matcoh,mytitle="",mylabels=NA, save.plots = FALSE, mod)
{
	origpar <- par(no.readonly = TRUE)
  nc <- NCOL(matcoh)
	my.cor <- cor(matcoh,use="pairwise.complete.obs")
	my.cor.round <- round(my.cor,2)
	par(mfcol=c(nc,nc))
	par(oma=c(0,0,3,1),mar=c(1,1,0,0))
	for (i in 1:nc)
	{
		for (j in nc:1)
		{
			if (i == j)
			{
				plot(1:10,1:10,type='n',axes=FALSE)
				text(5,5,paste0("age-",mod$ages.lab[i]),cex=1.4)
			}
			if (i < j)
			{
				if (!is.na(my.cor[i,j]))
				{
					plot(matcoh[,i],matcoh[,j],axes=FALSE) # make sure have some data to plot
					xx <- matcoh[,i]
					yy <- matcoh[,j]
					my.fit <- lm(yy~xx)
					if (!is.na(my.fit$coefficients[2])) abline(my.fit,col="red")
					xrng <- data.frame(xx = seq(min(xx,na.rm=T),max(xx,na.rm=T),length.out=100))
					zz <- predict(my.fit,xrng,interval="confidence")
					lines(xrng[,1],zz[,2],col="blue")
					lines(xrng[,1],zz[,3],col="blue")
				}
				if (is.na(my.cor[i,j]))
				{  # if not data, just make empty box
					plot(1:10,1:10,type='n',axes=FALSE)
				}
        box()
			}
			if (i > j)
			{
				plot(1:10,1:10,type='n',axes=FALSE)
				txt <- format(my.cor.round[i,j], nsmall=2)
				text(5,5,txt)
				box()
			}
		}
	}
	title(mytitle, outer=TRUE)
	# par(origpar)
	return(my.cor)
}
#-------------------------------------------------------------------------------

plot_catch_at_age_consistency <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
	# create plots of ages vs each other (correctly lagged) on log scale for catch by fleet
	cat.corr <- list()
  dat = mod$env$data
  rep = mod$rep
  n_years = length(mod$years)
	for (i in 1:dat$n_fleets)
	{
		if (dat$n_fleets == 1) title1 = "Catch"
		if (dat$n_fleets >= 2) title1 = paste("Catch for Fleet ",i, sep="")

		# get catch at age
    catchob = dat$catch_paa[i,,] * dat$agg_catch[,i]/apply(dat$catch_paa[i,,] * dat$waa[dat$waa_pointer_fleets[i],1:n_years,],1,sum)
    catchpr = rep$pred_catch_paa[1:n_years,i,] * exp(rep$pred_log_catch[1:n_years,i])/apply(rep$pred_catch_paa[1:n_years,i,] * dat$waa[dat$waa_pointer_fleets[i],1:n_years,],1,sum)
		# replace zeros with NA and take logs
		cob <- rep0log(catchob)
		cpr <- rep0log(catchpr)

		# make cohorts
		cob.coh <- makecohorts(cob)
		cpr.coh <- makecohorts(cpr)

		# make the plots
		if(do.tex) cairo_pdf(file.path(od, paste0("catch_at_age_consistency_obs_fleet",i,".pdf")), family = fontfam, height = 10, width = 10)
		if(do.png) png(filename = file.path(od, paste0("catch_at_age_consistency_obs_fleet",i,".png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
		cob.cor <- plotcoh(cob.coh,mytitle=paste(title1," Observed", sep=""),mod=mod)
		if(do.tex | do.png) dev.off()

		if(do.tex) cairo_pdf(file.path(od, paste0("catch_at_age_consistency_pred_fleet",i,".pdf")), family = fontfam, height = 10, width = 10)
		if(do.png) png(filename = file.path(od, paste0("catch_at_age_consistency_pred_fleet",i,".png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
		cpr.cor <- plotcoh(cpr.coh,mytitle=paste(title1," Predicted", sep=""),mod=mod)
		if(do.tex | do.png) dev.off()

		cat.corr[[i]] <- list(cob.cor,cpr.cor)
	}
	return(cat.corr)
}

#-------------------------------------------------------------------------------

convert_survey_to_at_age <- function(mod)
{
	# takes West Coast style surveys and converts them to catch at age matrices
  dat = mod$env$data
  rep = mod$rep
	index.mats <- list(ob = list(), pr = list())
	weight.mat.counter <- 0
	for (i in 1:dat$n_indices)
	{
		if (sum(dat$use_index_paa[,i])>0)
		{  # used age composition for the index
			# get the aggregate index observed and predicted time series
			agg.ob <- dat$agg_indices[which(dat$use_index_paa[,i]==1),i]
			agg.pr <- exp(rep$pred_log_indices[which(dat$use_index_paa[,i]==1),i]) # bias corrected

			# get proportions for correct years and ages only
			props.ob <- dat$index_paa[i,which(dat$use_index_paa[,i]==1),]
			props.pr <- rep$pred_index_paa[which(dat$use_index_paa[,i]==1),i,]

      # figure out units for aggregate and proportions
			agg.units <- dat$units_indices[i]
			prp.units <- dat$units_index_paa[i]

      waa <- dat$waa[dat$waa_pointer_indices[i],which(dat$use_index_paa[,i]==1),]
			# get weight (matrix if necessary)
			if (agg.units==1 | prp.units==1)
			{  # either in weight
			}

			# create index.obs and pred based on which of the four possible combinations of units is used for this index
			if (agg.units==1 && prp.units==1)
			{  # both in weight
				index.ob <- agg.ob * props.ob / waa
				index.pr <- agg.pr * props.pr / waa
			}
			if (agg.units==1 && prp.units==2)
			{  # agg in weight, props in numbers
        index.ob = props.ob * agg.ob/apply(props.ob * waa,1,sum)
        index.pr = props.pr * agg.pr/apply(props.pr * waa,1,sum)
				#index.ob <- wtprop2caa(agg.ob,waa,props.ob)  # use catch function
				#index.pr <- wtprop2caa(agg.pr,waa,props.pr)
			}
			if (agg.units==2 && prp.units==1)
			{  # agg in numbers, props in weight
				# need to search for correct agg total in weight to result in observed agg total in number
				# for now just use simple approximation that agg.wt = sum(waa*prop) *ctot and then solve using both in weight approach
				agg.wt.ob <- apply((waa * props.ob),1,sum) * agg.ob
				agg.wt.pr <- apply((waa * props.pr),1,sum) * agg.pr
				index.ob <- agg.wt.ob * props.ob / waa
				index.pr <- agg.wt.pr * props.pr / waa
			}
			if (agg.units==2 && prp.units==2)
			{  # both in numbers
				index.ob <- agg.ob * props.ob
				index.pr <- agg.pr * props.pr
			}

			# put matrices into full year matrix (ages only for selected ages though)
			# need to do this to account for missing years of data interspersed in time series
			index.ob.full <- index.pr.full <- matrix(NA, nrow=length(mod$years), ncol=dat$n_ages)
			index.ob.full[dat$use_index_paa[,i]==1,] <- index.ob
			index.pr.full[dat$use_index_paa[,i]==1,] <- index.pr

			# save the results for this index
			index.mats$ob[[i]] <- index.ob.full
			index.mats$pr[[i]] <- index.pr.full
		}
		else
		{  # cannot use this index
			index.mats$ob[[i]] <- NA
			index.mats$pr[[i]] <- NA
		}
	}
	return(index.mats)
}
#-------------------------------------------------------------------------------

plot_index_at_age_consistency <- function(mod, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
  dat = mod$env$data
  n_years = length(mod$years)
	# now loop through indices, check to make sure actually estimating proportions at age
	# also need to use only the age range selected in the proportions
	index.corr <- list()

	# convert the west coast style indices to catch at age matrices
	index.mats <- convert_survey_to_at_age(mod)

	# loop through all the indices
	for (ind in 1:dat$n_indices)
	{
		if (sum(dat$use_index_paa[,ind]) >0)
		{  # used age composition for the index
			title1 <- paste("Index ",ind, sep="")

			# replace zeros with NA and take logs
			iob <- rep0log(index.mats$ob[[ind]])
			ipr <- rep0log(index.mats$pr[[ind]])

			# make cohorts
			iob.coh <- makecohorts(iob)
			ipr.coh <- makecohorts(ipr)


			# create labels for plot (if necessary)
			mylabels <- paste0("age-",mod$ages.lab, sep="")

			# make the plots
			if(do.tex) cairo_pdf(file.path(od, paste0("catch_at_age_consistency_obs_index",ind,".pdf")), family = fontfam, height = 10, width = 10)
			if(do.png) png(filename = file.path(od, paste0("catch_at_age_consistency_obs_index",ind,".png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
			iob.cor <- plotcoh(iob.coh,mytitle=paste(title1," Observed", sep=""),mylabels,mod=mod)
			if(do.tex | do.png) dev.off()

			if(do.tex) cairo_pdf(file.path(od, paste0("catch_at_age_consistency_pred_index",ind,".pdf")), family = fontfam, height = 10, width = 10)
			if(do.png) png(filename = file.path(od, paste0("catch_at_age_consistency_pred_index",ind,".png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
			ipr.cor <- plotcoh(ipr.coh,mytitle=paste(title1," Predicted", sep=""),mylabels,mod=mod)
			if(do.tex | do.png) dev.off()

			index.corr[[ind]] <- list(iob.cor,ipr.cor)
		}
	}
	return(index.corr)
}
#plot_index_at_age_consistency(ssm)

#-------------------------------------------------------------------------------

find_peak_age <- function(cohmat)
{
	# determines peak age within each cohort and replaces younger ages with NA
	ages <- seq(1,NCOL(cohmat))
	temp.sum <- apply(cohmat,1,sum, na.rm = TRUE)
	for (i in 1:NROW(cohmat))
	{
		if (temp.sum[i] > 0)
		{  # necessary to avoid cohorts that are all NA
			age.peak <- max(ages[cohmat[i,]==max(cohmat[i,],na.rm=TRUE)], na.rm=TRUE)  # find the peak age
			if (age.peak > 1) cohmat[i,1:(age.peak-1)] <- NA  # replace ages < peak.age with NA
		}
	}
	return(cohmat)
}
#-------------------------------------------------------------------------------

calc_Z_cohort <- function(cohmat)
{
	# calculate Z along cohort using linear regression and return point estimate and 80% confidence interval
	ages <- seq(1,NCOL(cohmat))
	z <- matrix(NA, nrow=NROW(cohmat), ncol=3)
	for (i in 1:NROW(cohmat))
	{
		if (length(cohmat[i,which(!is.na(cohmat[i,]))]) >= 2)
		{  # make sure there are at least 2 data points for point estimate
			z.fit <- lm(cohmat[i,]~ages)  # linear regression of cohort abundance vs age
			z[i,1] <- -1 * z.fit$coefficients[2]  # note change in sign for Z
			if (length(cohmat[i,!is.na(cohmat[i,])]) >= 3)
			{  # need at least 3 data points for CI
				z[i,2:3] <- -1 * rev(confint(z.fit, "ages", level=0.80))  # note change in sign and order for Z CI
			}
		}
	}
	return(z)
}
#-------------------------------------------------------------------------------

plot_catch_curves_for_catch <- function(mod, first.age=-999, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
	# create catch curve plots for catch by fleet
  origpar <- par(no.readonly = TRUE)
	lastyr <- max(mod$years)
  dat = mod$env$data
  rep = mod$rep
  n_years = length(mod$years)
	cohort <- (min(mod$years) - dat$n_ages-1):(lastyr+dat$n_ages-1)
	ages <- 1:dat$n_ages
  my.col = rep(mypalette(5),50)
	#my.col <- rep(c("blue","red","green","orange","gray50"),50)
	for (i in 1:dat$n_fleets)
	{
		if (dat$n_fleets == 1) title1 = "Catch"
		if (dat$n_fleets >= 2) title1 = paste0("Catch for Fleet ",i)

		# get catch at age
    catchob = dat$catch_paa[i,,] * dat$agg_catch[,i]/apply(dat$catch_paa[i,,] * dat$waa[dat$waa_pointer_fleets[i],1:n_years,],1,sum)
    catchpr = rep$pred_catch_paa[1:n_years,i,] * exp(rep$pred_log_catch[1:n_years,i])/apply(rep$pred_catch_paa[1:n_years,i,] * dat$waa[dat$waa_pointer_fleets[i],1:n_years,],1,sum)

		# replace zeros with NA and take logs
		cob <- rep0log(catchob)
		cpr <- rep0log(catchpr)

		# make cohorts
		cob.coh <- makecohorts(cob)
		cpr.coh <- makecohorts(cpr)

		# drop plus group
		cob.coh[,dat$n_ages] <- NA
		cpr.coh[,dat$n_ages] <- NA

		first.age.label <- 1
		if (first.age==1) title1 <- paste(title1," First Age = 1", sep="")

		# determine which ages to use for each cohort (default)
		if (first.age == -999)
    {
      cob.coh <- find_peak_age(cob.coh)
      cpr.coh <- find_peak_age(cpr.coh)
      first.age.label <- "find_peak"
      title1 <- paste0(title1," (Peak Age)")
		}

		# or drop youngest ages based on user control
		if (first.age > 1)
    {
      cob.coh[,1:(first.age-1)] <- NA
      cpr.coh[,1:(first.age-1)] <- NA
      title1 <- paste0(title1," First Age = ",first.age)
      first.age.label <- first.age
		}

		# compute Z by cohort
		z.ob <- calc_Z_cohort(cob.coh)
		z.pr <- calc_Z_cohort(cpr.coh)

		# make the plots
		if(do.tex) cairo_pdf(file.path(od, paste0("catch_curves_fleet",i,"_obs.pdf")), family = fontfam, height = 10, width = 10)
		if(do.png) png(filename = file.path(od, paste0("catch_curves_fleet",i,"_obs.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
		par(oma=c(1,1,1,1),mar=c(4,4,1,0.5),mfrow=c(2,1))
		plot(cohort,cohort,type='n',ylim=range(c(cob.coh,cpr.coh),na.rm=TRUE),xlab="",ylab="Log(Catch)", main=paste0(title1," Observed"))
		grid(col = gray(0.7))
		for (j in 1:NROW(cob.coh))
		{
			lines(cohort[j]:(cohort[j]+dat$n_ages-1),cob.coh[j,],type='p',lty=1,pch=1:dat$n_ages,col="gray50")
			lines(cohort[j]:(cohort[j]+dat$n_ages-1),cob.coh[j,],type='l',lty=1,col=my.col[j])
		}
		Hmisc::errbar(cohort,z.ob[,1],z.ob[,3],z.ob[,2],xlab="Year Class",ylab="Z",ylim=range(c(z.ob,z.pr),na.rm=T))
		grid(col = gray(0.7))
		if(do.tex | do.png) dev.off()

		if(do.tex) cairo_pdf(file.path(od, paste0("catch_curves_fleet",i,"_pred.pdf")), family = fontfam, height = 10, width = 10)
		if(do.png) png(filename = file.path(od, paste0("catch_curves_fleet",i,"_pred.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
		par(oma=c(1,1,1,1),mar=c(4,4,1,0.5),mfrow=c(2,1))
		plot(cohort,cohort,type='n',ylim=range(c(cob.coh,cpr.coh),na.rm=TRUE),xlab="",ylab="Log(Catch)", main=paste0(title1," Predicted"))
		grid(col = gray(0.7))
		for (j in 1:length(cob.coh[,1]))
		{
			lines(cohort[j]:(cohort[j]+dat$n_ages-1),cpr.coh[j,],type='p',lty=1,pch=1:dat$n_ages,col="gray50")
			lines(cohort[j]:(cohort[j]+dat$n_ages-1),cpr.coh[j,],type='l',lty=1,col=my.col[j])
		}
		Hmisc::errbar(cohort,z.pr[,1],z.pr[,3],z.pr[,2],xlab="Year Class",ylab="Z",ylim=range(c(z.ob,z.pr),na.rm=T))
		grid(col = gray(0.7))
		if(do.tex | do.png) dev.off()

    # write out .csv files for Z, one file for each fleet
		colnames(z.ob) <-c("Z.obs","low.80%", "high.80%")
		#write.csv(z.ob, file=paste(od,"Z.Ob.Fleet.",i,".csv", sep=""), row.names=cohort)

		colnames(z.pr) <-c("Z.pred","low.80%", "high.80%")
		#write.csv(z.pr, file=paste(od,"Z.Pr.Fleet.",i,".csv", sep=""), row.names=cohort)
	}
  if(!(do.tex | do.png)) par(origpar)
}
#plot_catch_curves_for_catch(ssm)

#-------------------------------------------------------------------------------

plot_catch_curves_for_index <- function(mod, first.age=-999, do.tex = FALSE, do.png = FALSE, fontfam="", res = 72, od)
{
	# create catch curve plots for each west coast style index
  origpar <- par(no.readonly = TRUE)
  lastyr <- max(mod$years)
  dat = mod$env$data
  rep = mod$rep
  ages = 1:dat$n_ages
  n_ages = dat$n_ages
  cohort <- (min(mod$years)-n_ages-min(ages)):(lastyr+n_ages-min(ages))
	my.col <- rep(c("blue","red","green","orange","gray50"),50)

	# convert the west coast style indices to catch at age matrices
	index.mats <- convert_survey_to_at_age(mod)

	# loop through all the indices
	for (ind in 1:dat$n_indices)
	{
		if (sum(dat$use_index_paa[,ind]) > 0)
		{  # used age composition for the index
			title1 <- paste0("Index ",ind)
			# replace zeros with NA and take logs
			iob <- rep0log(index.mats$ob[[ind]])
			ipr <- rep0log(index.mats$pr[[ind]])

			# make cohorts
			iob.coh <- makecohorts(iob)
			ipr.coh <- makecohorts(ipr)

			# drop plus group
			iob.coh[,NCOL(iob.coh)] <- NA
			ipr.coh[,NCOL(iob.coh)] <- NA

			first.age.label <- 1
			if (first.age==1) title1 <- paste0(title1," First Age = 1")

			# determine which ages to use for each cohort (default)
			if (first.age == -999)
			{
				iob.coh <- find_peak_age(iob.coh)
				ipr.coh <- find_peak_age(ipr.coh)
				first.age.label <- "find_peak"
				title1 <- paste0(title1," (Peak Age)")
			}

			# or drop youngest ages based on user control
			if (first.age > min(ages))
			{
				iob.coh[,1:(first.age-min(ages))] <- NA
				ipr.coh[,1:(first.age-min(ages))] <- NA
				title1 <- paste0(title1," First Age = ",first.age)
				first.age.label <- first.age
			}

			# compute Z by cohort
			z.ob <- calc_Z_cohort(iob.coh)
			z.pr <- calc_Z_cohort(ipr.coh)

			# make the plots
			if(!(all(is.na(iob.coh)) & all(is.na(ipr.coh))))
			{
			  if(do.tex) cairo_pdf(file.path(od, paste0("catch_curves_index",ind,"_obs.pdf")), family = fontfam, height = 10, width = 10)
			  if(do.png) png(filename = file.path(od, paste0("catch_curves_index",ind,"_obs.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
			  par(oma=c(1,1,1,1),mar=c(4,4,1,0.5),mfrow=c(2,1))
			  plot(cohort,cohort,type='n',ylim=range(c(iob.coh,ipr.coh),na.rm=T),xlab="",ylab="Log(Index)",
          main=paste0(title1," Observed"))
				grid(col = gray(0.7))
				for (i in 1:length(iob.coh[,1]))
				{
					lines(cohort[i]:(cohort[i]+n_ages-1),iob.coh[i,],type='p',lty=1,pch=1:n_ages,col="gray50")
					lines(cohort[i]:(cohort[i]+n_ages-1),iob.coh[i,],type='l',lty=1,col=my.col[i])
				}
				Hmisc::errbar(cohort,z.ob[,1],z.ob[,3],z.ob[,2],xlab="Year Class",ylab="Z",ylim=range(c(z.ob,z.pr),na.rm=TRUE))
				grid(col = gray(0.7))
				if(do.tex | do.png) dev.off()

				if(do.tex) cairo_pdf(file.path(od, paste0("catch_curves_index",ind,"_pred.pdf")), family = fontfam, height = 10, width = 10)
				if(do.png) png(filename = file.path(od, paste0("catch_curves_index",ind,"_pred.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
				par(oma=c(1,1,1,1),mar=c(4,4,1,0.5),mfrow=c(2,1))
				plot(cohort,cohort,type='n',ylim=range(c(iob.coh,ipr.coh),na.rm=T),xlab="",ylab="Log(Index)", main=paste0(title1," Predicted"))
				grid(col = gray(0.7))
				for (i in 1:NROW(iob.coh))
				{
					lines(cohort[i]:(cohort[i]+n_ages-1),ipr.coh[i,],type='p',lty=1,pch=1:n_ages,col="gray50")
					lines(cohort[i]:(cohort[i]+n_ages-1),ipr.coh[i,],type='l',lty=1,col=my.col[i])
				}
				Hmisc::errbar(cohort,z.pr[,1],z.pr[,3],z.pr[,2],xlab="Year Class",ylab="Z",ylim=range(c(z.ob,z.pr),na.rm=TRUE))
				grid(col = gray(0.7))
				if(do.tex | do.png) dev.off()
			}

			# write out .csv files for Z, one file for each fleet
			colnames(z.ob) <-c("Z.obs","low.80%", "high.80%")
			#write.csv(z.ob, file=paste(od,"Z.Ob.Index.",ind,".csv", sep=""), row.names=cohort)

			colnames(z.pr) <-c("Z.pred","low.80%", "high.80%")
			#write.csv(z.pr, file=paste(od,"Z.Pr.Index.",ind,".csv", sep=""), row.names=cohort)
		}
	}   # end loop over n_indices
	if(!(do.tex | do.png)) par(origpar)
}
#plot_catch_curves_for_index(ssm)

#-------------------------------------------------------------------------------
plot.ecov.diagnostic <- function(mod, use.i, plot.pad = FALSE, do.tex = FALSE, do.png = FALSE, fontfam="", od){
  origpar <- par(no.readonly = TRUE)
  dat = mod$env$data
  years <- seq(from=dat$year1_Ecov, by=1, length.out=dat$n_years_Ecov)
  years_full <- seq(from=dat$year1_Ecov, by=1, length.out=dat$n_years_Ecov+dat$n_years_proj_Ecov)

  ecov.pred = mod$rep$Ecov_x
  ecov.obs = dat$Ecov_obs[1:dat$n_years_Ecov,,drop=F]
  # ecov.obs.sig = mod$rep$Ecov_obs_sigma # Ecov_obs_sigma now a derived quantity in sdrep
  if(class(mod$sdrep)[1] == "sdreport"){
    sdrep = summary(mod$sdrep)
  } else {
    sdrep = mod$sdrep
  }
  ecov.obs.sig = mod$rep$Ecov_obs_sigma # Ecov_obs_sigma is filled with fixed, or estimated values (fe or re) for each covariate depending on the respective options
  # if("Ecov_obs_logsigma" %in% names(mod$env$par)){
  #   ecov.obs.sig = matrix(exp(sdrep[rownames(sdrep) %in% "Ecov_obs_logsigma",1]), ncol=dat$n_Ecov) # all the same bc obs_sig_var --> 0
  #   if(dim(ecov.obs.sig)[1] == 1) ecov.obs.sig = matrix(rep(ecov.obs.sig, dim(ecov.obs)[1]), ncol=dat$n_Ecov, byrow=T)
  # } else {
  #   ecov.obs.sig = exp(mod$input$par$Ecov_obs_logsigma)
  # }
  ecov.use = dat$Ecov_use_obs[1:dat$n_years_Ecov,,drop=F]
  ecov.obs.sig = ecov.obs.sig[1:dat$n_years_Ecov,,drop=F]
  ecov.obs.sig[ecov.use == 0] <- NA
  ecov.pred.se = matrix(sdrep[rownames(sdrep) %in% "Ecov_x",2], ncol=dat$n_Ecov)

  # default: don't plot the padded entries that weren't used in ecov likelihood
  if(!plot.pad) ecov.obs[ecov.use == 0] <- NA

  ecov.res = (ecov.obs - ecov.pred[1:dat$n_years_Ecov,]) / ecov.obs.sig # standard residual (obs - pred)

  if(!missing(use.i)) ecovs <- use.i
  else ecovs <- 1:dat$n_Ecov
  plot.colors = mypalette(dat$n_Ecov)
  for (i in ecovs)
  {
    if(do.tex) cairo_pdf(file.path(od, paste0("Ecov_",i,"_diagnostic.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("Ecov_",i,'_diagnostic.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)

    m <- rbind(c(1,1), c(2,3))
    layout(m)
    par(mar=c(4,4,2,0), oma=c(0,0,0.5,0.5))

    ecov.pred.low <- ecov.pred[,i] - 1.96 * ecov.pred.se[,i]
    ecov.pred.high <- ecov.pred[,i] + 1.96 * ecov.pred.se[,i]
    ecov.low <- ecov.obs[,i] - 1.96 * ecov.obs.sig[,i]
    ecov.high <- ecov.obs[,i] + 1.96 * ecov.obs.sig[,i]
    y.min <- ifelse(min(ecov.low,na.rm=T) < 0, 1.1*min(ecov.low,na.rm=T), 0.9*min(ecov.low,na.rm=T))
    y.max <- ifelse(max(ecov.high,na.rm=T) < 0, 0.9*max(ecov.high,na.rm=T), 1.1*max(ecov.high,na.rm=T))
    if(max(ecov.pred[,i],na.rm=T) > y.max) y.max <- max(ecov.pred[,i],na.rm=T)
    if(min(ecov.pred[,i],na.rm=T) < y.min) y.min <- min(ecov.pred[,i],na.rm=T)
    plot(years_full, ecov.pred[,i], type='n', xlab="Year", ylab=unlist(dat$Ecov_label)[i],
         ylim=c(y.min, y.max))
    polygon(c(years_full,rev(years_full)), c(ecov.pred.low, rev(ecov.pred.high)), col=adjustcolor(plot.colors[i], alpha.f=0.4), border = "transparent")
    arrows(years, ecov.low, years, ecov.high, length=0)
    points(years, ecov.obs[,i], pch=19)
    lines(years_full, ecov.pred[,i], col=plot.colors[i], lwd=3)
    title (paste0("Ecov ",i, ": ",unlist(dat$Ecov_label)[i]), outer=T, line=-1)
    if(dat$n_years_proj_Ecov > 0) abline(v=tail(years,1), lty=2)

    plot(years, ecov.res[,i], type='h', lwd=2, col=plot.colors[i], xlab="Year", ylab="Std. Residual")
    abline(h=0)

    hist(ecov.res[,i], breaks=10, plot=T, xlab="Std. Residual", ylab="Probability Density", freq=F, main=NULL)
    if(do.tex | do.png) dev.off() else par(origpar)
  }
}

#-------------------------------------------------------------------------------
# 2D tile plot by age and year (e.g. selAA, MAA)
plot.tile.age.year <- function(mod, type="selAA", do.tex = FALSE, do.png = FALSE, fontfam="", od){
  dat = mod$env$data
  rep = mod$rep
  years = mod$years
  n_years = length(years)
  n_ages = dat$n_ages
  ages <- 1:n_ages
  ages.lab = 1:n_ages
  if(!is.null(mod$ages.lab)) ages.lab = mod$ages.lab

  # selAA for all blocks using facet_wrap
  if(type=="selAA"){ 
    n_selblocks <- length(rep$selAA)
    sel_mod <- c("age-specific","logistic","double-logistic","decreasing-logistic")[dat$selblock_models]
    sel_re <- c("no","IID","AR1","AR1_y","2D AR1")[dat$selblock_models_re]
    df.selAA <- data.frame(matrix(NA, nrow=0, ncol=n_ages+2))
    colnames(df.selAA) <- c(paste0("Age_",1:n_ages),"Year","Block")
    for(i in 1:n_selblocks){
      tmp = as.data.frame(rep$selAA[[i]])
      tmp$Year <- years
      colnames(tmp) <- c(paste0("Age_",1:n_ages),"Year")
      tmp$Block = paste0("Block ",i,": ", sel_mod[i]," (",sel_re[i]," random effects)")
      df.selAA <- rbind(df.selAA, tmp)
    }
    df.plot <- df.selAA %>% tidyr::pivot_longer(-c(Year,Block),
              names_to = "Age", 
              names_prefix = "Age_",
              names_ptypes = list(Age = character()),
              values_to = "Selectivity")
    df.plot$Age <- as.factor(as.integer(df.plot$Age))
    levels(df.plot$Age) = ages.lab
    df.plot$Block <- factor(as.character(df.plot$Block), levels=names(table(df.plot$Block)))

    if(do.tex) cairo_pdf(file.path(od, paste0("SelAA_tile.pdf")), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, paste0("SelAA_tile.png")), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
      print(ggplot2::ggplot(df.plot, ggplot2::aes(x=Year, y=Age, fill=Selectivity)) + 
        ggplot2::geom_tile() +
        ggplot2::scale_x_continuous(expand=c(0,0)) +
        ggplot2::scale_y_discrete(expand=c(0,0)) + #, breaks = function(x) unique(floor(pretty(seq(0, (max(x) + 1) * 1.1))))) +        
#        ggplot2::scale_y_continuous(expand=c(0,0), breaks = function(x) unique(floor(pretty(seq(0, (max(x) + 1) * 1.1))))) +        
        ggplot2::theme_bw() + 
        ggplot2::facet_wrap(~Block, dir="v") +
        viridis::scale_fill_viridis())
    if(do.tex | do.png) dev.off()
  }

  # MAA
  if(type=="MAA"){ 
    if(mod$env$data$do_proj == 1){
      years_full = mod$years_full
    } else {
      years_full = years
    }
    df.MAA <- as.data.frame(rep$MAA)
    df.MAA$Year <- years_full
    colnames(df.MAA) <- c(paste0("Age_",1:n_ages),"Year")
    df.plot <- df.MAA %>% tidyr::pivot_longer(-Year,
              names_to = "Age", 
              names_prefix = "Age_",
              names_ptypes = list(Age = character()),
              values_to = "M")
    df.plot$Age <- as.factor(as.integer(df.plot$Age))
    levels(df.plot$Age) = ages.lab

    if(do.tex) cairo_pdf(file.path(od, paste0("MAA_tile.pdf")), family = fontfam, height = 5, width = 10)
    if(do.png) png(filename = file.path(od, paste0("MAA_tile.png")), width = 10*144, height = 5*144, res = 144, pointsize = 12, family = fontfam)
      print(ggplot2::ggplot(df.plot, ggplot2::aes(x=Year, y=Age, fill=M)) + 
        ggplot2::geom_tile() +
        ggplot2::scale_x_continuous(expand=c(0,0)) +
        ggplot2::scale_y_discrete(expand=c(0,0)) + #, breaks = function(x) unique(floor(pretty(seq(0, (max(x) + 1) * 1.1))))) +        
#        ggplot2::scale_y_continuous(expand=c(0,0), breaks = function(x) unique(floor(pretty(seq(0, (max(x) + 1) * 1.1))))) +
        ggplot2::theme_bw() + 
        viridis::scale_fill_viridis())
    if(do.tex | do.png) dev.off()
  }
}  

#pdf of a univariate logit-normal with any min and max
dlogitnorm = function(p,mu,sd,min,max) {
  logitp = log((p-min)/(max-p))
  (exp(-(logitp - mu)^2/(2 * sd^2))/(sd * sqrt(2*pi))) * (max-p)/((p-min) * (max-p))
}

plot_q_prior_post = function(mod, do.tex = F, do.png = F, fontfam="", od){
  origpar <- par(no.readonly = TRUE)
  ind = which(mod$input$data$use_q_prior == 1)
  if(length(ind) & "sdrep" %in% names(mod)){
    logit_q = cbind(as.list(mod$sdrep, "Est")$q_prior_re, as.list(mod$sdrep, "Std")$q_prior_re)
    ht = 10
    wd = 10*length(ind)
    priorq = approx_postq = list()
    if(do.tex) cairo_pdf(file.path(od, "q_prior_post.pdf"), family = fontfam, height = ht, width = wd)
    if(do.png) png(filename = file.path(od, "q_prior_post.png"), width = wd*144, height = ht*144, res = 144, pointsize = 12, family = fontfam)
    par(mar=c(4,4,3,2), oma=c(1,1,1,1), mfrow=c(1,length(ind)))
    pal = viridisLite::viridis(n=2)
    for(i in ind) {
      qmax = mod$input$data$q_upper[i]
      x = seq(0.001,qmax,0.001)
      y = log(x-0) - log((qmax-x))
      approx_postq[[i]] = dnorm(y, logit_q[i,1], logit_q[i,2])
      priorq[[i]] = dlogitnorm(x, mod$parList$logit_q[i], 0.3, 0, 1000) 
      maxx = max(x[which(approx_postq[[i]] > 1e-5)], x[which(priorq[[i]] > 1e-5)], na.rm = T)
      plot(x,priorq[[i]], type = 'l', xlab = "q", ylab = "pdf", col = pal[1], lwd = 2, ylim = c(0,max(approx_postq[[i]],priorq[[i]], na.rm =T)), xlim = c(0,maxx))
      lines(x,approx_postq[[i]], col = pal[2], lwd = 2)
      ci = qmax/(1+ exp(-(logit_q[i,1] + c(-1,1)*qnorm(0.975) * logit_q[i,2])))
      abline(v= ci, lty = 2)
      legend("topright", legend = c("prior", "approx. posterior", "95% CI"), col = c(pal, "black"), lty = c(1,1,2), lwd = 2)
      mtext(paste0("Index ", ind), side = 3, line = 1, outer = F, cex = 1.5)
    }
    if(do.tex | do.png) dev.off() else par(origpar)
  }  
}

plot_q = function(mod, do.tex = F, do.png = F, fontfam = '', od){
  origpar <- par(no.readonly = TRUE)
  if("sdrep" %in% names(mod)){
    if("q_re" %in% mod$input$random) se = as.list(mod$sdrep, "Std. Error", report=TRUE)$logit_q_mat
    else se = t(matrix(as.list(mod$sdrep, "Std. Error")$logit_q, nrow = NCOL(mod$rep$logit_q_mat), 
      ncol = NROW(mod$rep$logit_q_mat)))
    logit_q_lo = mod$rep$logit_q_mat - qnorm(0.975)*se
    logit_q_hi = mod$rep$logit_q_mat + qnorm(0.975)*se
    q = t(mod$input$data$q_lower + (mod$input$data$q_upper - mod$input$data$q_lower)/(1+exp(-t(mod$rep$logit_q_mat))))
    q_lo = t(mod$input$data$q_lower + (mod$input$data$q_upper - mod$input$data$q_lower)/(1+exp(-t(logit_q_lo))))
    q_hi = t(mod$input$data$q_lower + (mod$input$data$q_upper - mod$input$data$q_lower)/(1+exp(-t(logit_q_hi))))
    if(do.tex) cairo_pdf(file.path(od, "q_time_series.pdf"), family = fontfam, height = 10, width = 10)
    if(do.png) png(filename = file.path(od, "q_time_series.png"), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
    par(mar=c(4,4,3,2), oma=c(1,1,1,1))
    pal = viridisLite::viridis(n=mod$input$data$n_indices)
    ymax = max(c(q,q_hi), na.rm = TRUE)
    plot(mod$years_full, q[,1], type = 'n', lwd = 2, col = pal[1], ylim = c(0,ymax), ylab = "q", xlab = "Year")
    for( i in 1:mod$input$data$n_indices){
      lines(mod$years_full, q[,i], lwd = 2, col = pal[i])
      polygon(c(mod$years_full,rev(mod$years_full)), c(q_lo[,i],rev(q_hi[,i])), col=adjustcolor(pal[i], alpha.f=0.4), border = "transparent")
    }
    legend("topright", legend = paste0("Index ", 1:mod$input$data$n_indices), lwd = 2, col = pal, lty = 1)
  }
  if(do.tex | do.png) dev.off() else par(origpar)
}
timjmiller/wham documentation built on April 28, 2024, 5:39 a.m.