R/SSplotSelex.R

In r4ss: R Code for Stock Synthesis

#' Plot selectivity
#'
#' Plot selectivity, including retention and other quantities, with additional
#' plots for time-varying selectivity.
#'
#'
#' @template replist
#' @param fleets Optional vector to subset fleets for which to make plots
#' @param infotable Optional table of information controlling appearance of
#' plot and legend. Is produced as output and can be modified and entered as
#' input.
#' @param fleetnames Optional replacement for fleenames used in data file
#' @param sizefactors Which elements of the factors column of SIZE_SELEX should
#' be included in plot of selectivity across multiple fleets?
#' @param agefactors Which elements of the factors column of AGE_SELEX should
#' be included in plot of selectivity across multiple fleets?
#' @param years Which years for selectivity are shown in multi-line plot
#' (default = last year of model).
#' @param minyr optional input for minimum year to show in plots
#' @param maxyr optional input for maximum year to show in plots
#' @param season Which season (if seasonal model) for selectivity shown in
#' multi-line plot (default = 1).
#' @param sexes Optional vector to subset genders for which to make plots
#' (1=females, 2=males)
#' @param selexlines Vector to select which lines get plotted. values are 1.
#' Selectivity, 2. Retention, 3. Discard mortality, 4. Keep.
#' @param subplot Vector controlling which subplots to create.
#' Numbering of subplots is as follows,
#'
#' *Plots with all fleets grouped together*
#' \itemize{
#'   \item 1 selectivity at length in end year for all fleets shown together
#'   \item 2 selectivity at age in end year for all fleets shown together
#'           (this includes both age-based selectivity "Asel" and age values derived
#'           from length-based, "Asel2". You can choose only one using
#'           "agefactors" if needed.)
#' }
#'
#' *Plots of time-varying length-based selectivity*
#' \itemize{
#'   \item 3 selectivity at length time-varying surface
#'   \item 4 selectivity at length time-varying contour
#'   \item 5 retention at length time-varying surface
#'   \item 6 retention at length time-varying surface
#'   \item 7 discard mortality time-varying surface
#'   \item 8 discard mortality time-varying contour
#' }
#'
#' *Selectivity at length in end year by fleet*
#' \itemize{
#'   \item 9 selectivity, retention, and discard mortality at length in ending year
#' }
#'
#' *Plots of time-varying age-based selectivity*
#' \itemize{
#'   \item 11 selectivity at age time-varying surface
#'   \item 12 selectivity at age time-varying contour
#' }
#'
#' *Selectivity at age in end year by fleet*
#' \itemize{
#'   \item 13 selectivity at age in ending year if time-varying
#'   \item 14 selectivity at age in ending year if NOT time-varying
#'   \item 15 matrix of selectivity deviations for semi-parametric selectivity
#' }
#'
#' *Selectivity for both/either age or length*
#' \itemize{
#'   \item 21 selectivity at age and length contour with overlaid growth curve
#'   \item 22 selectivity with uncertainty if requested at end of control file
#' }
#' @param skipAgeSelex10 Exclude plots for age selectivity type 10 (selectivity
#' = 1.0 for all ages beginning at age 1)?
#' @param lwd Line widths for plots
#' @param spacepoints number of years between points shown on top of lines (for
#' long timeseries, points every year get mashed together)
#' @param staggerpoints number of years to stagger the first point (if
#' `spacepoints > 1`) for each line (so that adjacent lines have points in
#' different years)
#' @param legendloc location of legend. See ?legend for more info.
#' @template pwidth
#' @template pheight
#' @template punits
#' @template ptsize
#' @template res
#' @param plot Plot to active plot device?
#' @param print Print to PNG files?
#' @param add Add to existing plot (not yet implemented)
#' @param labels vector of labels for plots (titles and axis labels)
#' @param col1 color for female growth curve
#' @param col2 color for male growth curve
#' @param cex.main character expansion for plot titles
#' @template mainTitle
#' @template mar
#' @param plotdir Directory where PNG files will be written. By default it will
#' be the directory where the model was run.
#' @template verbose
#' @param showmain Deprecated, use mainTitle instead.
#' @author Ian Stewart, Ian Taylor
#' @export
#' @seealso [SS_plots()], [SS_output()]
SSplotSelex <-
  function(replist, infotable = NULL,
           fleets = "all", fleetnames = "default",
           sizefactors = c("Lsel"),
           agefactors = c("Asel", "Asel2"),
           years = "endyr",
           minyr = -Inf,
           maxyr = Inf,
           season = 1,
           sexes = "all",
           selexlines = 1:6,
           subplot = 1:25,
           skipAgeSelex10 = TRUE,
           plot = TRUE, print = FALSE, add = FALSE,
           labels = c(
             "Length (cm)", # 1
             "Age (yr)", # 2
             "Year", # 3
             "Selectivity", # 4
             "Retention", # 5
             "Discard mortality"
           ), # 6
           col1 = "red", col2 = "blue", lwd = 2,
           spacepoints = 5,
           staggerpoints = 1,
           legendloc = "bottomright",
           pwidth = 6.5, pheight = 5.0, punits = "in",
           res = 300, ptsize = 10,
           cex.main = 1,
           mainTitle = TRUE,
           showmain = lifecycle::deprecated(),
           mar = NULL,
           plotdir = "default",
           verbose = TRUE) {

    # empty table into which information on line types etc. might be copied
    infotable2 <- NULL

    nsexes <- replist[["nsexes"]]
    nseasons <- replist[["nseasons"]]
    nfleets <- replist[["nfleets"]]
    lbinspop <- replist[["lbinspop"]]
    nlbinspop <- replist[["nlbinspop"]]
    sizeselex <- replist[["sizeselex"]]
    ageselex <- replist[["ageselex"]]
    accuage <- replist[["accuage"]]
    startyr <- replist[["startyr"]]
    endyr <- replist[["endyr"]]
    FleetNames <- replist[["FleetNames"]]
    growdat <- replist[["endgrowth"]]
    growthCVtype <- replist[["growthCVtype"]]
    mainmorphs <- replist[["mainmorphs"]]
    nareas <- replist[["nareas"]]
    ngpatterns <- replist[["ngpatterns"]]
    derived_quants <- replist[["derived_quants"]]

    if (lifecycle::is_present(showmain)) {
      lifecycle::deprecate_warn(
        when = "1.41.1",
        what = "SSplotSelex(showmain)",
        details = "Please use mainTitle instead. Ability to use showmain will be dropped in next release."
      )
      mainTitle <- showmain
    }

    # message about skipping plots
    if (is.null(ageselex)) {
      message("Skipping age-based selectivity plots: no output available")
    }
    if (is.null(sizeselex)) {
      message("Skipping length-based selectivity plots: no output available")
    }

    # table to store information on each plot
    plotinfo <- NULL

    if (plotdir == "default") {
      plotdir <- replist[["inputs"]][["dir"]]
    }

    ians_blues <- c(
      "white", "grey", "lightblue", "skyblue", "steelblue1",
      "slateblue", topo.colors(6), "blue", "blue2", "blue3",
      "blue4", "black"
    )

    if (fleets[1] == "all") {
      fleets <- 1:nfleets
    } else {
      if (length(intersect(fleets, 1:nfleets)) != length(fleets)) {
        return("Input 'fleets' should be 'all' or a vector of values between 1 and nfleets.")
      }
    }
    # note lower-case value is the one used below (either equal to vector from replist, or input by user)
    if (fleetnames[1] == "default") fleetnames <- FleetNames

    if (sexes[1] == "all") {
      sexes <- 1:nsexes
    } else {
      if (length(intersect(sexes, 1:nsexes)) != length(sexes)) {
        return("Input 'sexes' should be 'all' or a vector of values between 1 and nsexes.")
      }
    }
    if (years[1] == "endyr") years <- endyr

    # set default plot margins
    if (is.null(mar)) {
      if (mainTitle) {
        mar <- c(5, 4, 4, 2) + 0.1
      } else {
        mar <- c(5, 4, 2, 2) + 0.1
      }
    }

    ################################################################################
    ### make plot of selectivity at length for all fleets together
    # make plots
    plotAllSel <- function(factor = "Lsel") {
      # first subset values
      if (factor %in% unique(sizeselex[["Factor"]])) {
        agebased <- FALSE
        allselex <- sizeselex[sizeselex[["Factor"]] == factor &
          sizeselex[["Fleet"]] %in% fleets &
          sizeselex[["Sex"]] %in% sexes, ]
      }
      if (factor %in% unique(ageselex[["Factor"]])) {
        agebased <- TRUE
        allselex <- ageselex[ageselex[["Factor"]] == factor &
          ageselex[["Seas"]] == season &
          ageselex[["Fleet"]] %in% fleets &
          ageselex[["Sex"]] %in% sexes, ]
      }
      if (!factor %in% unique(c(sizeselex[["Factor"]], ageselex[["Factor"]]))) {
        stop(
          "Factor '", factor, "' not found in age- or length-based selectivity. ",
          "This may be due to having 'detailed age-structured reports' ",
          "turned off in the starter file."
        )
      }
      if (nrow(allselex) == 0) {
        stop("Combination of season, fleets, & sexes didn't produce any results")
      }
      # figure out which fleets have time-varying qunatities
      time <- rep(FALSE, nfleets)
      for (ifleet in fleets) {
        time[ifleet] <- any(apply(
          allselex[allselex[["Fleet"]] == ifleet &
            allselex[["Yr"]] %in% (startyr:endyr), ],
          2, function(x) {
            any(x != x[1])
          }
        ))
      }
      if (any(time)) {
        if (length(years) > 1 & length(fleets) > 1) {
          message("plot not yet configured to work well with multiple years and multiple fleets")
        }
        # do a bunch of tedious filtering to get unique year ranges
        inputyears <- years
        years <- NULL
        years2 <- NULL
        year_ranges <- NULL
        for (i in 1:length(inputyears)) {
          if (inputyears[i] >= startyr) {
            newyear <- min(endyr, allselex[["Yr"]][allselex[["Yr"]] >= inputyears[i]])
            newyear2 <- max(startyr, allselex[["Yr"]][allselex[["Yr"]] <= inputyears[i]])
            if (newyear2 <= newyear) {
              newyear_range <- paste(newyear2, "-", newyear, sep = "")
              if (newyear == newyear2 & newyear > startyr - 3) newyear_range <- newyear
              if (!newyear_range %in% year_ranges) {
                years <- c(years, newyear)
                years2 <- c(years2, newyear2)
                year_ranges <- c(year_ranges, newyear_range)
              }
            }
          }
        }
        if (all(years2 == startyr & years == endyr)) {
          years <- endyr
          years2 <- startyr
          year_ranges <- paste(startyr, "-", endyr, sep = "")
        }
        bad <- rep(FALSE, length(years))
        for (i in 1:length(years)) {
          y <- years[i]
          y2 <- years2[i]
          if (sum(years == y) > 1) bad[years == y & years2 == y] <- TRUE
          if (sum(years2 == y2) > 1) bad[years == y2 & years2 == y2] <- TRUE
        }
        years <- years[!bad]
        years2 <- years2[!bad]
        year_ranges <- year_ranges[!bad]

        if ((startyr - 3) %in% inputyears) {
          years <- c(years, startyr - 3)
          year_ranges <- c(year_ranges, "Benchmarks")
        }
      } else {
        year_ranges <- ""
      }
      allselex <- allselex2 <- allselex[allselex[["Yr"]] %in% years, ]
      if (nrow(allselex) == 0) {
        stop("No values found for this combination of years and factor")
      }

      # do some processing
      Sex <- allselex[["Sex"]]
      if (!agebased) {
        allselex <- allselex[, -(1:5)]
        xlab <- labels[1]
      }
      if (agebased) {
        allselex <- allselex[, -(1:7)]
        xlab <- labels[2]
      }
      if (!is.null(infotable)) {
        infotable2 <- infotable
        good <- Sex %in% infotable[["Sex"]]
        allselex <- allselex[good, ]
        allselex2 <- allselex2[good, ]
        if (nrow(infotable2) != nrow(allselex)) {
          stop("Problem with input 'infotable'. Number of rows doesn't match.")
        }
      } else {
        # make table of info for each row (unless it is supplied already)
        infotable2 <- allselex2[c("Fleet", "Sex", "Yr")]
        infotable2[["ifleet"]] <- NA
        infotable2[["FleetName"]] <- fleetnames[infotable2[["Fleet"]]]
        infotable2[["longname"]] <- infotable2[["FleetName"]]
        for (i in 1:nrow(infotable2)) {
          infotable2[["Yr_range"]][i] <- year_ranges[years == infotable2[["Yr"]][i]]
        }

        if (length(unique(infotable2[["Yr"]])) > 1) {
          infotable2[["longname"]] <- paste(infotable2[["FleetName"]], infotable2[["Yr_range"]])
        }
        # check for whether there are differences between males and females
        twosex <- all(1:2 %in% infotable2[["Sex"]]) &&
          any(allselex[infotable2[["Sex"]] == 1, ] != allselex[infotable2[["Sex"]] == 2, ])
        if (!twosex) { # show only sex with lowest number if no differences between sexes
          good <- infotable2[["Sex"]] == min(infotable2[["Sex"]])
          allselex <- allselex[good, ]
          allselex2 <- allselex2[good, ]
          infotable2 <- infotable2[good, ]
        } else {
          infotable2[["longname"]] <- paste(infotable2[["longname"]], c("(f)", "(m)")[infotable2[["Sex"]]])
        }

        # add index from 1 up to number of fleets plotted
        allfleets <- sort(unique(infotable2[["Fleet"]]))
        for (ifleet in 1:length(allfleets)) {
          infotable2[["ifleet"]][infotable2[["Fleet"]] == allfleets[ifleet]] <- ifleet
        }
        # choose colors
        colvec <- rich.colors.short(length(allfleets))
        infotable2[["col"]] <- colvec[infotable2[["ifleet"]]]
        # choose line types
        infotable2[["lty"]] <- 1
        # either line by Sex
        infotable2[["lwd"]] <- lwd
        if (twosex) infotable2[["lty"]] <- infotable2[["Sex"]]
        # or line by year (with line width by Sex)
        allyears <- sort(unique(infotable2[["Yr"]]))
        if (length(allyears) > 1) {
          for (iyear in 1:length(allyears)) {
            infotable2[["lty"]][infotable2[["Yr"]] == allyears[iyear]] <- iyear
          }
          if (twosex) infotable2[["lwd"]][infotable2[["Sex"]] == 2] <- lwd / 2
        }
        # choose plot characters
        infotable2[["pch"]] <- infotable2[["ifleet"]] %% 25
      }

      main <- factor
      if (factor == "Lsel") main <- paste("Length-based selectivity")
      if (factor == "Asel") main <- paste("Age-based selectivity")
      if (factor == "Asel2") main <- paste("Derived age-based from length-based selectivity")
      if (factor == "Ret") main <- paste("Retention")
      if (length(fleets) > 1) main <- paste(main, "by fleet")
      if (length(fleets) == 1) main <- paste(main, "for", fleetnames[fleets])

      if (length(unique(infotable2[["Yr"]])) == 1) {
        main <- paste(main, "in", unique(infotable2[["Yr"]]))
      }

      bins <- as.numeric(names(allselex))

      # make empty plot
      if (!add) {
        par(mar = mar)
        plot(0,
          xlim = range(bins), ylim = c(0, 1), type = "n",
          main = ifelse(mainTitle, main, ""), cex.main = cex.main,
          xlab = xlab, ylab = labels[4]
        )
      }
      # add grey lines
      abline(h = 0, col = "grey")
      abline(h = 1, col = "grey")
      # add lines for selectivities
      matplot(
        x = bins, y = t(allselex), col = infotable2[["col"]],
        lty = infotable2[["lty"]], lwd = infotable2[["lwd"]], type = "l", add = TRUE
      )
      # add points on top of lines (first subsetting to optionally show fewer points)
      allselex2 <- allselex
      if (spacepoints > 0) {
        for (iline in 1:nrow(allselex)) {
          allselex2[iline, (1:ncol(allselex)) %%
            spacepoints != (staggerpoints * iline) %% spacepoints] <- NA
        }
        matplot(
          x = bins, y = t(allselex2), col = infotable2[["col"]],
          lwd = infotable2[["lwd"]], pch = infotable2[["pch"]], type = "p", add = TRUE
        )
      } else {
        # if 'spacepoints' is less than or equal to 0, don't show points
        infotable2[["pch"]] <- NA
      }
      # add legend
      if (nrow(infotable2) > 1) {
        legend(legendloc,
          inset = c(0, 0.05),
          legend = infotable2[["longname"]],
          col = infotable2[["col"]],
          seg.len = 4,
          lty = infotable2[["lty"]],
          pch = infotable2[["pch"]],
          lwd = infotable2[["lwd"]],
          bty = "n"
        )
      }
      return(infotable2)
    }

    if (1 %in% subplot & !is.null(sizeselex)) {
      for (ifactor in 1:length(sizefactors)) {
        if (plot) {
          infotable2 <- plotAllSel(factor = sizefactors[ifactor])
        }
        if (print) {
          file <- paste("sel01_multiple_fleets_length", ifactor, ".png", sep = "")
          caption <- "Selectivity at length for multiple fleets."
          plotinfo <- save_png(
            plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
            pheight = pheight, punits = punits, res = res, ptsize = ptsize,
            caption = caption
          )
          infotable2 <- plotAllSel(factor = "Lsel")
          dev.off()
        }
      }
    }
    if (2 %in% subplot & !is.null(ageselex)) {
      # remove factor == "Asel" if all age-based selectivity == 1
      if ("Asel" %in% agefactors &&
        all(ageselex[
          ageselex[["Factor"]] == "Asel",
          paste(0:replist[["accuage"]])
        ] == 1)) {
        agefactors <- setdiff(agefactors, "Asel")
        message("Skipping plot of age-based selectivity as all values = 1.0")
      }

      if (length(agefactors) > 0) {
        for (ifactor in 1:length(agefactors)) {
          factor <- agefactors[ifactor]
          if (plot) {
            infotable2 <- plotAllSel(factor = factor)
          }
          if (print) {
            file <- paste("sel02_multiple_fleets_age", ifactor, ".png", sep = "")
            caption <- "Selectivity at age for multiple fleets."
            if (factor == "Asel2") {
              caption <- paste(
                "Selectivity at age derived from selectivity at",
                "length for multiple fleets."
              )
            }
            plotinfo <- save_png(
              plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
              pheight = pheight, punits = punits, res = res, ptsize = ptsize,
              caption = caption
            )
            infotable2 <- plotAllSel(factor = factor)
            dev.off()
          }
        } # end loop over age factors
      } # end check for any agefactors
    }

    ################################################################################
    ### loop over fleets and sexes to make individual plot of length-based patterns

    # first check if any of these plots are requested
    if (any(3:9 %in% subplot) & !is.null(sizeselex)) {
      # selex and retention
      for (i in fleets)
      {
        for (m in sexes)
        {
          if (m == 1 & nsexes == 1) sextitle1 <- "Time-"
          if (m == 1 & nsexes == 2) sextitle1 <- "Female time-"
          if (m == 2) sextitle1 <- "Male time-"
          if (m == 1 & nsexes == 1) sextitle2 <- "Ending"
          if (m == 1 & nsexes == 2) sextitle2 <- "Female ending"
          if (m == 2) sextitle2 <- "Male ending"
          intret <- sizeselex[sizeselex[["Factor"]] == "Ret" &
            sizeselex[["Yr"]] != startyr - 3 &
            sizeselex[["Sex"]] == m, ]
          intmort <- sizeselex[sizeselex[["Factor"]] == "Mort" &
            sizeselex[["Yr"]] != startyr - 3 &
            sizeselex[["Sex"]] == m, ]
          intkeep <- sizeselex[sizeselex[["Factor"]] == "Keep" &
            sizeselex[["Yr"]] != startyr - 3 &
            sizeselex[["Sex"]] == m, ]
          intdead <- sizeselex[sizeselex[["Factor"]] == "Dead" &
            sizeselex[["Yr"]] != startyr - 3 &
            sizeselex[["Sex"]] == m, ]
          intselex <- sizeselex[sizeselex[["Factor"]] == "Lsel" &
            sizeselex[["Yr"]] != startyr - 3 &
            sizeselex[["Sex"]] == m, ]
          plotselex <- intselex[intselex[["Fleet"]] == i, ]
          plotret <- intret[intret[["Fleet"]] == i, ]
          plotmort <- intmort[intmort[["Fleet"]] == i, ]

          # test for time-varying length selectivity
          time <- any(apply(plotselex[-c(1, nrow(plotselex)), -(1:5)], 2, function(x) {
            any(x != x[1])
          }))
          if (time) {
            x <- lbinspop
            subset <- plotselex[["Yr"]] >= minyr & plotselex[["Yr"]] <= maxyr
            y <- plotselex[["Yr"]][subset]
            z <- plotselex[subset, -(1:5)]
            z <- matrix(as.numeric(as.matrix(z)), ncol = ncol(z))
            z <- t(z)
            main <- paste(sextitle1, "varying selectivity for ", fleetnames[i], sep = "")
            if (plot) {
              if (3 %in% subplot) {
                persp(x, y, z,
                  col = "white", xlab = labels[1],
                  ylab = labels[3], zlab = labels[4], expand = 0.5,
                  box = TRUE, main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, ticktype = "detailed",
                  phi = 35, theta = -10
                )
              }
              if (4 %in% subplot) {
                contour(x, y, z,
                  nlevels = 5, xlab = labels[1],
                  ylab = labels[3], main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, col = ians_blues, lwd = lwd
                )
              }
            }
            if (print) {
              if (3 %in% subplot) {
                file <- paste("sel03_len_timevary_surf_flt", i, "sex", m, ".png", sep = "")
                caption <- paste("Surface plot of", main)
                plotinfo <- save_png(
                  plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                  pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                  caption = caption
                )
                persp(x, y, z,
                  col = "white", xlab = labels[1], ylab = labels[3],
                  zlab = labels[4], expand = 0.5, box = TRUE,
                  main = ifelse(mainTitle, main, ""), cex.main = cex.main,
                  ticktype = "detailed", phi = 35, theta = -10
                )
                dev.off()
              }
              if (4 %in% subplot) {
                file <- paste("sel04_len_timevary_contour_flt", i, "sex", m, ".png", sep = "")
                caption <- paste("Countour plot of", main)
                plotinfo <- save_png(
                  plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                  pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                  caption = caption
                )
                contour(x, y, z,
                  nlevels = 5, xlab = labels[1], ylab = labels[3],
                  main = ifelse(mainTitle, main, ""), cex.main = cex.main,
                  col = ians_blues, lwd = lwd
                )
                dev.off()
              }
            }
          }
          # test for time-varying length retention
          time2 <- any(apply(plotret[-nrow(plotret), -(1:5)], 2, function(x) {
            any(x != x[1])
          }))
          if (time2) {
            x <- lbinspop
            subset <- intret[["Yr"]] >= minyr & intret[["Yr"]] <= maxyr
            y <- intret[["Yr"]][subset & intret[["Fleet"]] == i]
            z <- intret[subset & intret[["Fleet"]] == i, -(1:5)]
            z <- matrix(as.numeric(as.matrix(z)), ncol = ncol(z))
            z <- t(z)
            main <- paste(sextitle1, "varying retention for ", fleetnames[i], sep = "")
            if (plot) {
              if (5 %in% subplot) {
                persp(x, y, z,
                  col = "white", xlab = labels[1],
                  ylab = labels[3], zlab = labels[5], expand = 0.5,
                  box = TRUE, main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, ticktype = "detailed",
                  phi = 35, theta = -10
                )
              }
              if (6 %in% subplot) {
                contour(x, y, z,
                  nlevels = 5, xlab = labels[1],
                  ylab = labels[3], main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, col = ians_blues, lwd = lwd
                )
              }
            }
            if (print) {
              if (5 %in% subplot) {
                file <- paste("sel05_timevary_ret_surf_flt", i, "sex", m, ".png", sep = "")
                caption <- paste("Surface plot of", main)
                plotinfo <- save_png(
                  plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                  pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                  caption = caption
                )
                persp(x, y, z,
                  col = "white", xlab = labels[1],
                  ylab = labels[3], zlab = labels[5], expand = 0.5,
                  box = TRUE, main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, ticktype = "detailed",
                  phi = 35, theta = -10
                )
                dev.off()
              }
              if (6 %in% subplot) {
                file <- paste("sel06_timevary_ret_contour_flt", i, "sex", m, ".png", sep = "")
                caption <- paste("Countour plot of", main)
                plotinfo <- save_png(
                  plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                  pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                  caption = caption
                )
                contour(x, y, z,
                  nlevels = 5, xlab = labels[1],
                  ylab = labels[3], main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, col = ians_blues, lwd = lwd
                )
                dev.off()
              }
            }
          }
          # test for time-varying discard mortality rates
          time3 <- any(apply(plotmort[-nrow(plotmort), -(1:5)], 2, function(x) {
            any(x != x[1])
          }))
          if (time3) {
            x <- lbinspop
            subset <- intmort[["Yr"]] >= minyr & intmort[["Yr"]] <= maxyr
            y <- intmort[["Yr"]][subset & intmort[["Fleet"]] == i]
            z <- intmort[subset & intmort[["Fleet"]] == i, -(1:5)]
            z <- matrix(as.numeric(as.matrix(z)), ncol = ncol(z))
            z <- t(z)
            main <- paste(sextitle1, "varying discard mortality for ", fleetnames[i], sep = "")
            if (plot) {
              if (7 %in% subplot) {
                persp(x, y, z,
                  col = "white", xlab = labels[1], ylab = labels[3], zlab = labels[6],
                  expand = 0.5, box = TRUE,
                  main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, ticktype = "detailed",
                  phi = 35, theta = -10, zlim = c(0, max(z))
                )
              }
              if (8 %in% subplot) {
                contour(x, y, z,
                  nlevels = 5, xlab = labels[1], ylab = labels[3],
                  main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, col = ians_blues, lwd = lwd
                )
              }
            }
            if (print) {
              if (7 %in% subplot) {
                file <- paste("sel07_timevary_mort_surf_flt", i, "sex", m, ".png", sep = "")
                caption <- paste("Surface plot of", main)
                plotinfo <- save_png(
                  plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                  pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                  caption = caption
                )
                persp(x, y, z,
                  col = "white", xlab = labels[1],
                  ylab = labels[3], zlab = labels[6], expand = 0.5,
                  box = TRUE, main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, ticktype = "detailed", phi = 35, theta = -10
                )
                dev.off()
              }
              if (8 %in% subplot) {
                file <- paste("sel08_timevary_mort_contour_flt", i, "sex", m, ".png", sep = "")
                caption <- paste("Surface plot of", main)
                plotinfo <- save_png(
                  plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                  pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                  caption = caption
                )
                contour(x, y, z,
                  nlevels = 5, xlab = labels[1],
                  ylab = labels[3], main = ifelse(mainTitle, main, ""),
                  cex.main = cex.main, col = ians_blues, lwd = lwd
                )
                dev.off()
              }
            }
          }

          # make plot of end year selectivity (with retention and discard mortality if used)
          endselex <- plotselex[plotselex[["Yr"]] == endyr, -(1:5)]

          plotret <- plotret[nrow(plotret), -(1:5)] # final year only
          ylab <- labels[4]
          bins <- as.numeric(names(endselex))
          vals <- as.numeric(paste(endselex))
          retvals <- as.numeric(plotret)
          main <- paste(sextitle2, " year selectivity for ", fleetnames[i], sep = "")
          selfunc <- function() {
            # determine whether retention was used
            intret2 <- intret[intret[["Fleet"]] == i, ]
            retchecktemp <- as.vector(unlist(intret2[1, ]))
            retcheck <- as.numeric(retchecktemp[6:length(retchecktemp)])
            if (is.na(sum(retcheck))) retcheckuse <- 0
            # if minimum retention is less than 1, show additional stuff in plot
            if (!is.na(sum(retcheck))) retcheckuse <- 1 - min(retcheck)

            # make plot
            if (!add) {
              par(mar = mar)
              plot(bins, vals,
                xlab = labels[1], ylim = c(0, 1),
                main = ifelse(mainTitle, main, ""),
                cex.main = cex.main, ylab = "", type = "n"
              )
            }
            abline(h = 0, col = "grey")
            abline(h = 1, col = "grey")
            if (1 %in% selexlines) lines(bins, vals, type = "o", col = col2, cex = 1.1)
            if (retcheckuse > 0) {
              # if retention, then add additional lines & legend
              useret <- intret[intret[["Fleet"]] == i, ]
              usekeep <- intkeep[intkeep[["Fleet"]] == i, ]
              usemort <- intmort[intmort[["Fleet"]] == i, ]
              usedead <- intdead[intdead[["Fleet"]] == i, ]
              if (endyr %in% as.numeric(useret[["Yr"]])) {
                useyr <- endyr
              } else {
                useyr <- max(as.numeric(useret[["Yr"]]))
              }
              plotret <- useret[useret[["Yr"]] == useyr, ]
              plotkeep <- usekeep[usekeep[["Yr"]] == useyr, ]
              plotmort <- usemort[usemort[["Yr"]] == useyr, ]
              plotdead <- usedead[usedead[["Yr"]] == useyr, ]
              # compute discard as function of size: selectivity*(1 - retention)
              plotdisc <- plotret
              plotdisc[-(1:5)] <- vals * (1 - plotret[, -(1:5)])
              # add additional lines if requested
              if (2 %in% selexlines) {
                lines(as.numeric(as.vector(names(plotret)[-(1:5)])),
                  as.numeric(as.character(plotret[1, -(1:5)])),
                  col = "red", type = "o", pch = 3, cex = .9
                )
                ylab <- paste(ylab, ", Retention", sep = "")
              }
              if (3 %in% selexlines) {
                lines(as.numeric(as.vector(names(plotmort)[-(1:5)])),
                  as.numeric(as.character(plotmort[1, -(1:5)])),
                  col = "orange", type = "o", pch = 4, cex = .9
                )
                ylab <- paste(ylab, ", Mortality", sep = "")
              }
              if (4 %in% selexlines) {
                lines(as.numeric(as.vector(names(plotkeep)[-(1:5)])),
                  as.numeric(as.character(plotkeep[1, -(1:5)])),
                  col = "purple", type = "o", pch = 2, cex = .9
                )
              }
              if (5 %in% selexlines) {
                lines(as.numeric(as.vector(names(plotdead)[-(1:5)])),
                  as.numeric(as.character(plotdead[1, -(1:5)])),
                  col = "green3", type = "o", pch = 5, cex = .9
                )
              }
              if (6 %in% selexlines) {
                lines(as.numeric(as.vector(names(plotdead)[-(1:5)])),
                  as.numeric(as.character(plotdisc[1, -(1:5)])),
                  col = "grey50", type = "o", pch = 6, cex = .9
                )
              }
              # add legend
              legend(legendloc,
                inset = c(0, 0.05), bty = "n",
                c(
                  labels[4], labels[5], labels[6], "Keep = Sel*Ret",
                  "Dead = Sel*(Ret+(1-Ret)*Mort)", "Discard = Sel*(1-Ret)"
                )[selexlines],
                lty = 1, col = c("blue", "red", "orange", "purple", "green3", "grey50")[selexlines],
                pch = c(1, 3, 4, 2, 5, 6)[selexlines], pt.cex = c(1.1, .9, .9, .9, .9, .9)[selexlines]
              )
            }
            mtext(ylab, side = 2, line = 3)
          }
          # make plot if selectivity is not constant at 0 or 1 for all bins
          if ((min(vals) < 1 & max(vals) > 0) |
            (!is.na(diff(range(retvals))) && diff(range(retvals)) != 0)) {
            if (9 %in% subplot) {
              if (plot) selfunc()
              if (print) {
                file <- paste("sel09_len_flt", i, "sex", m, ".png", sep = "")
                caption <- main
                plotinfo <- save_png(
                  plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                  pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                  caption = caption
                )
                selfunc()
                dev.off()
              }
            }
          }
        } # sexes
      } # fleets
    } # check for any of the plots in this section requested

    ################################################################################
    ### loop over fleets and sexes to make individual plot of age-based patterns
    if (any(11:14 %in% subplot) & !is.null(ageselex)) {

      # Age based selex
      ylab <- labels[4]
      for (facnum in 1) {
        factor <- c("Asel", "Asel2")[facnum]
        for (i in fleets) {
          for (m in sexes) {
            if (m == 1 & nsexes == 1) sextitle1 <- "Time-"
            if (m == 1 & nsexes == 2) sextitle1 <- "Female time-"
            if (m == 2) sextitle1 <- "Male time-"
            if (m == 1 & nsexes == 1) sextitle2 <- "Ending"
            if (m == 1 & nsexes == 2) sextitle2 <- "Female ending"
            if (m == 2) sextitle2 <- "Male ending"
            ageselexcols <- (1:ncol(ageselex))[names(ageselex) %in% as.character(0:accuage)]
            plotageselex <- ageselex[ageselex[["Factor"]] == factor & ageselex[["Fleet"]] == i &
              ageselex[["Yr"]] != startyr - 3 & ageselex[["Sex"]] == m, ]

            # test for time-varying age selectivity
            time <- any(apply(
              plotageselex[-c(1, nrow(plotageselex)), ageselexcols],
              2, function(x) {
                any(x != x[1])
              }
            ))
            if (time) {
              if ((min(as.numeric(as.vector(t(plotageselex[, -(1:7)])))) < 1)) {
                subset <- as.numeric(plotageselex[["Yr"]]) >= minyr &
                  as.numeric(plotageselex[["Yr"]]) <= maxyr
                x <- seq(0, accuage, by = 1)
                y <- as.numeric(plotageselex[["Yr"]])[subset]
                z <- plotageselex[subset, -(1:7)]
                z <- matrix(as.numeric(as.matrix(z)), ncol = ncol(z))
                z <- t(z)
                main <- paste(sextitle1, "varying selectivity for ", fleetnames[i], sep = "")
                if (plot) {
                  if (11 %in% subplot) {
                    persp(x, y, z,
                      col = "white", xlab = labels[2],
                      ylab = labels[3], zlab = ylab, expand = 0.5,
                      box = TRUE, main = ifelse(mainTitle, main, ""), cex.main = cex.main,
                      ticktype = "detailed", phi = 35, theta = -10
                    )
                  }
                  if (12 %in% subplot) {
                    contour(x, y, z,
                      nlevels = 5, xlab = labels[2],
                      main = ifelse(mainTitle, main, ""),
                      cex.main = cex.main, col = ians_blues, lwd = lwd
                    )
                  }
                }
                if (print) {
                  if (11 %in% subplot) {
                    file <- paste("sel11_timevary_surf_flt", i, "sex", m, ".png", sep = "")
                    caption <- main
                    plotinfo <- save_png(
                      plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                      pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                      caption = caption
                    )
                    persp(x, y, z,
                      col = "white", xlab = labels[2],
                      ylab = labels[3], zlab = ylab, expand = 0.5,
                      box = TRUE, main = ifelse(mainTitle, main, ""),
                      cex.main = cex.main, ticktype = "detailed",
                      phi = 35, theta = -10
                    )
                    dev.off()
                  }
                  if (12 %in% subplot) {
                    file <- paste("sel12_timevary_contour_flt", i, "sex", m, ".png", sep = "")
                    caption <- main
                    plotinfo <- save_png(
                      plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                      pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                      caption = caption
                    )
                    contour(x, y, z,
                      nlevels = 5, xlab = labels[2],
                      main = ifelse(mainTitle, main, ""),
                      cex.main = cex.main, col = ians_blues, lwd = lwd
                    )
                    dev.off()
                  }
                }
                plotageselex2 <- plotageselex[plotageselex[["Yr"]] %in% c(max(as.numeric(plotageselex[["Yr"]]))), ]
                plotageselex2 <- plotageselex2[, -(1:7)]
                main <- paste(sextitle2, " year selectivity for ", fleetnames[i], sep = "")
                endselfunc <- function() {
                  if (!add) {
                    par(mar = mar)
                    plot(as.numeric(names(plotageselex2)),
                      as.numeric(paste(c(plotageselex2))),
                      xlab = labels[2], ylim = c(0, 1), main = ifelse(mainTitle, main, ""),
                      cex.main = cex.main, ylab = ylab,
                      type = "n", col = col2, cex = 1.1
                    )
                  }
                  lines(as.numeric(names(plotageselex2)),
                    as.numeric(paste(c(plotageselex2))),
                    type = "o", col = col2, cex = 1.1
                  )
                  abline(h = 0, col = "grey")
                }
                if (13 %in% subplot) {
                  if (plot) {
                    endselfunc()
                  }
                  if (print) {
                    file <- paste("sel13_age_flt", i, "sex", m, ".png", sep = "")
                    caption <- main
                    plotinfo <- save_png(
                      plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                      pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                      caption = caption
                    )
                    endselfunc()
                    dev.off()
                  }
                }
              }
            }
            if (!time) {
              plotageselex <- plotageselex[plotageselex[["Yr"]] == max(plotageselex[["Yr"]]), ]
              plotageselex <- plotageselex[, -(1:7)]
              vals <- as.numeric(paste(c(plotageselex)))
              # test for constant across all values
              doplot <- nrow(plotageselex) > 0 && diff(range(vals)) != 0
              if (doplot & skipAgeSelex10) {
                # skip if selectivity type 10 is used (0.0 for age-0, 1.0 for other ages)
                doplot <- !(vals[1] == 0 & all(vals[-1] == 1))
              }
              if (doplot) {
                main <- paste(sextitle2, " year selectivity for ", fleetnames[i], sep = "")

                endselfunc2 <- function() {
                  if (!add) {
                    par(mar = mar)
                    plot((as.numeric(names(plotageselex))), vals,
                      xlab = labels[2], ylim = c(0, 1),
                      main = ifelse(mainTitle, main, ""), cex.main = cex.main,
                      ylab = ylab, type = "n"
                    )
                  }
                  lines((as.numeric(names(plotageselex))), vals, type = "o", col = col2, cex = 1.1)
                  abline(h = 0, col = "grey")
                }

                if (14 %in% subplot) {
                  if (plot) endselfunc2()
                  if (print) {
                    file <- paste("sel14_age_flt", i, "sex", m, ".png", sep = "")
                    caption <- main
                    plotinfo <- save_png(
                      plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                      pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                      caption = caption
                    )
                    endselfunc2()
                    dev.off()
                  }
                }
              } # end if
            } # end if not time varying
          } # sexes
        } # fleets
        flush.console()
      } # factor (Asel vs. Asel2)
    } # check for any of the plots in this section requested

    ################################################################################
    ### Matrix of selectivity deviations for semi-parametric (2D-AR1) selectivity

    if (15 %in% subplot & !is.null(replist[["seldev_matrix"]])) {
      seldev_pars <- replist[["seldev_pars"]]
      seldev_matrix <- replist[["seldev_matrix"]]
      # define color palette
      devcol.fn <- colorRampPalette(colors = c("red", "white", "blue"))

      # define function to make an image plot
      seldev_func <- function(m, mar = c(4.1, 4.1, 1, 1)) {
        bins <- as.numeric(colnames(m))
        years <- as.numeric(rownames(m))
        par(mar = mar)
        image(
          x = bins, y = years, z = t(m), col = devcol.fn(10),
          xlab = names(dimnames(m))[2],
          ylab = names(dimnames(m))[1],
          axes = FALSE,
          ylim = rev(range(years) + c(-0.5, 0.5))
        )
        axis(1, at = bins)
        axis(2, at = years, las = 1)
        box()
      }

      for (imatrix in 1:length(seldev_matrix)) {
        label <- names(seldev_matrix)[imatrix]
        main <- gsub(pattern = "_", replacement = " ", x = label)
        main <- gsub(pattern = "seldevs", replacement = "selectivity deviations", x = main)

        if (plot) {
          seldev_func(m = seldev_matrix[[imatrix]], mar = c(5, 4, 4, 1) + 0.1)
          title(main = ifelse(mainTitle, main, ""))
        }
        if (print) {
          file <- paste("sel15_", label, ".png", sep = "")
          caption <- gsub(pattern = "selectivity ", replacement = "", x = main)
          caption <- paste0(
            caption,
            " for semi-parametric (2D-AR1) selectivity. ",
            "Blue value are positive deviations and red values negative. ",
            "The matrix of values is available in the list created by ",
            "<code>SS_output()</code> as <code>$seldev_matrix</code> which ",
            "is a list with an element for each combination of fleet and length or ",
            "age which uses the semi-parametric selectivity."
          )
          plotinfo <- save_png(
            plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
            pheight = pheight, punits = punits, res = res, ptsize = ptsize,
            caption = caption
          )
          seldev_func(m = seldev_matrix[[imatrix]])
          dev.off()
        }
      } # end loop over matrices
    } # end subplot


    ################################################################################
    ### Selectivity contours over age and length shown with growth curve
    if (21 %in% subplot &
      !is.null(ngpatterns) &&
      ngpatterns == 1 &
      !is.null(growdat) &
      !is.null(sizeselex) &
      !is.null(ageselex) &
      all(!is.na(lbinspop))
    ) { # need to connect growth patterns to fleets in future
      # subsetting for one season only. This could be replaced
      #   by info on the growth within the season when each fleet operates.
      growdat <- growdat[growdat[["Seas"]] == season, ]
      if (nseasons > 1) {
        message(
          "Warning: plots showing growth curve with selectivity are using season ",
          season, " growth, which may not match the timing of the fishery."
        )
      }

      # Mid year mean length at age with 95% range of lengths (by sex if applicable)
      growdatF <- growdat[growdat[["Sex"]] == 1 & growdat[["Morph"]] == mainmorphs[1], ]
      growdatF[["Sd_Size"]] <- growdatF[["SD_Mid"]]
      if (growthCVtype == "logSD=f(A)") { # lognormal distribution of length at age
        growdatF[["high"]] <- qlnorm(0.975, meanlog = log(growdatF[["Len_Mid"]]), sdlog = growdatF[["Sd_Size"]])
        growdatF[["low"]] <- qlnorm(0.025, meanlog = log(growdatF[["Len_Mid"]]), sdlog = growdatF[["Sd_Size"]])
      } else { # normal distribution of length at age
        growdatF[["high"]] <- qnorm(0.975, mean = growdatF[["Len_Mid"]], sd = growdatF[["Sd_Size"]])
        growdatF[["low"]] <- qnorm(0.025, mean = growdatF[["Len_Mid"]], sd = growdatF[["Sd_Size"]])
      }
      if (nsexes > 1) {
        growdatM <- growdat[growdat[["Sex"]] == 2 & growdat[["Morph"]] == mainmorphs[2], ]
        growdatM[["Sd_Size"]] <- growdatM[["SD_Mid"]]
        if (growthCVtype == "logSD=f(A)") { # lognormal distribution of length at age
          growdatM[["high"]] <- qlnorm(0.975, meanlog = log(growdatM[["Len_Mid"]]), sdlog = growdatM[["Sd_Size"]])
          growdatM[["low"]] <- qlnorm(0.025, meanlog = log(growdatM[["Len_Mid"]]), sdlog = growdatM[["Sd_Size"]])
        } else { # normal distribution of length at age
          growdatM[["high"]] <- qnorm(0.975, mean = growdatM[["Len_Mid"]], sd = growdatM[["Sd_Size"]])
          growdatM[["low"]] <- qnorm(0.025, mean = growdatM[["Len_Mid"]], sd = growdatM[["Sd_Size"]])
        }
      }

      xlab <- labels[2]
      ylab <- labels[1]
      zlab <- labels[4]
      for (i in fleets)
      {
        for (m in sexes)
        {
          if (m == 1 & nsexes == 1) sextitle2 <- "Ending"
          if (m == 1 & nsexes == 2) sextitle2 <- "Female ending"
          if (m == 2) sextitle2 <- "Male ending"
          plotlenselex <- as.numeric(sizeselex[sizeselex[["Factor"]] == "Lsel" &
            sizeselex[["Yr"]] == endyr &
            sizeselex[["Fleet"]] == i &
            sizeselex[["Sex"]] == m, -(1:5)])
          # test if there is any length-based selectivity (otherwise plot is uninformative)
          if (any(plotlenselex != 1)) {
            plotageselex <- as.numeric(ageselex[ageselex[["Factor"]] == "Asel" &
              ageselex[["Yr"]] == endyr &
              ageselex[["Fleet"]] == i &
              ageselex[["Sex"]] == m, -(1:7)])
            # x here should probably be replaced by $Age_Mid or some more informative value
            x <- seq(0, accuage, by = 1)
            y <- lbinspop
            z <- plotageselex %o% plotlenselex # outer product of age- and length-selectivity

            main <- paste(sextitle2, " year selectivity and growth for ", fleetnames[i], sep = "")

            agelenselcontour <- function() {
              contour(x, y, z,
                nlevels = 5, xlab = xlab, ylab = ylab,
                main = ifelse(mainTitle, main, ""), cex.main = cex.main, col = ians_blues, lwd = lwd
              )
              if (m == 1) {
                lines(x, growdatF[["Len_Mid"]], col = "white", lwd = 5)
                lines(x, growdatF[["Len_Mid"]], col = col1, lwd = 3)
                lines(x, growdatF[["high"]], col = "white", lwd = 1, lty = 1)
                lines(x, growdatF[["high"]], col = col1, lwd = 1, lty = "dashed")
                lines(x, growdatF[["low"]], col = "white", lwd = 1, lty = 1)
                lines(x, growdatF[["low"]], col = col1, lwd = 1, lty = "dashed")
              }
              if (m == 2) {
                lines(x, growdatM[["Len_Mid"]], col = "white", lwd = 5)
                lines(x, growdatM[["Len_Mid"]], col = col2, lwd = 3)
                lines(x, growdatM[["high"]], col = "white", lwd = 1, lty = 1)
                lines(x, growdatM[["high"]], col = col2, lwd = 1, lty = "dashed")
                lines(x, growdatM[["low"]], col = "white", lwd = 1, lty = 1)
                lines(x, growdatM[["low"]], col = col2, lwd = 1, lty = "dashed")
              }
            }
            if (plot) {
              if (21 %in% subplot) agelenselcontour()
            }
            if (print) {
              if (21 %in% subplot) {
                file <- paste("sel21_agelen_contour_flt", i, "sex", m, ".png", sep = "")
                caption <- main
                plotinfo <- save_png(
                  plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
                  pheight = pheight, punits = punits, res = res, ptsize = ptsize,
                  caption = caption
                )
                agelenselcontour()
                dev.off()
              }
            }
          } # if there is any length-based selectivity
        } # sexes
      } # fleets
    } # end subplot

    ################################################################################
    ### Plot selectivity with uncertainty if "Extra SD reporting" requested in control file

    if (22 %in% subplot) {
      # get values from Extra SD reporting if created by request at bottom of control file
      rows <- grep("Selex_std", derived_quants[["Label"]])
      if (length(rows) > 0) {
        sel <- derived_quants[rows, ]
        names <- sel[["Label"]]
        splitnames <- strsplit(names, "_")
        namesDF <- as.data.frame(matrix(unlist(strsplit(names, "_")), ncol = 6, byrow = T))
        sel[["Fleet"]] <- as.numeric(as.character(namesDF[["V3"]]))
        sel[["Sex"]] <- as.character(namesDF[["V4"]])
        sel[["agelen"]] <- as.character(namesDF[["V5"]])
        sel[["bin"]] <- as.numeric(as.character(namesDF[["V6"]]))
        sel[["lower"]] <- pmax(qnorm(0.025, mean = sel[["Value"]], sd = sel[["StdDev"]]), 0) # trim at 0
        sel[["upper"]] <- pmin(qnorm(0.975, mean = sel[["Value"]], sd = sel[["StdDev"]]), 1) # trim at 1
        i <- sel[["Fleet"]][1]
        agelen <- sel[["agelen"]][1]
        xlab <- labels[1:2][1 + (sel[["agelen"]][1] == "A")] # decide label between length and age

        for (m in intersect(unique(sel[["Sex"]]), c("Fem", "Mal")[sexes])) {
          seltemp <- sel[sel[["Sex"]] == m, ]
          if (m == "Fem" & nsexes == 1) sextitle3 <- ""
          if (m == "Fem" & nsexes == 2) sextitle3 <- "females"
          if (m == "Mal") sextitle3 <- "males"
          main <- paste("Uncertainty in selectivity for", fleetnames[i], sextitle3)
          no0 <- seltemp[["StdDev"]] > 0.001

          if (FALSE) {
            # Ian T.: this is the beginning of code to add the full selectivity line,
            #        including bins for which no uncertainty was requested
            if (agelen == "L") {
              plotselex <- sizeselex[sizeselex[["Factor"]] == "Lsel" &
                ageselex[["Fleet"]] == i &
                sizeselex[["Sex"]] == m, ]
            }
            if (agelen == "A") {
              plotselex <- ageselex[ageselex[["Factor"]] == "Asel" &
                ageselex[["Fleet"]] == i &
                ageselex[["Sex"]] == m, ]
            }
          }

          plot_extra_selex_SD <- function() {
            if (!add) {
              par(mar = mar)
              plot(seltemp[["bin"]], seltemp[["Value"]],
                xlab = xlab, ylim = c(0, 1),
                main = ifelse(mainTitle, main, ""), cex.main = cex.main,
                ylab = labels[4], type = "n", col = col2, cex = 1.1,
                xlim = c(0, max(seltemp[["bin"]])),
              )
            }
            lines(seltemp[["bin"]], seltemp[["Value"]],
              xlab = xlab, ylim = c(0, 1),
              main = ifelse(mainTitle, main, ""), cex.main = cex.main,
              ylab = labels[4], type = "o", col = col2, cex = 1.1,
              xlim = c(0, max(seltemp[["bin"]]))
            )
            arrows(
              x0 = seltemp[["bin"]][no0], y0 = seltemp[["lower"]][no0],
              x1 = seltemp[["bin"]][no0], y1 = seltemp[["upper"]][no0],
              length = 0.01, angle = 90, code = 3, col = col2
            )
            abline(h = 0, col = "grey")
          }
          if (plot) {
            plot_extra_selex_SD()
          }
          if (print) {
            file <- paste("sel22_uncertainty", "sex", m, ".png", sep = "")
            caption <- main
            plotinfo <- save_png(
              plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
              pheight = pheight, punits = punits, res = res, ptsize = ptsize,
              caption = caption
            )
            plot_extra_selex_SD()
            dev.off()
          }
        }
      } # end test for presence of selectivity uncertainty output
    } # check check for subplot in list

    # return info on any PNG files created
    if (!is.null(plotinfo)) plotinfo[["category"]] <- "Sel"
    return(invisible(list(infotable = infotable2, plotinfo = plotinfo)))
  }