R/MMSE_plots.R

In tcarruth/MSEtool: Management Strategy Evaluation Toolkit

#' A fairly tidy time-series quantile plot
#'
#' @param x Matrix. A time series quantity [simulation, year]
#' @param p Numeric vector. The percentiles that are plotted (LB2,LB1,UB1,UB2). LB2 and UB2 are the outer percentiles, LB1 and UB1 are the inner percentiles.
#' @param yrs Numeric vector. The years corresponding to the indexing of x
#' @param qcol Character, color. The color of the inner percentile range
#' @param lcol Character, color. The color of the outer percentile range.
#' @param addline Logical. Should two individual simulations be added to the percentile plots?
#' @param ablines Numeric vector. Horizontal lines to be added to the plot.
#' @author T.Carruthers
#' @export
plotquant<-function(x,p=c(0.05,0.25,0.75,0.95),yrs,qcol,lcol,addline=T,ablines=NA){
  ny<-length(yrs)
  qs<-apply(x,2,quantile,p=p[c(1,4)],na.rm=T)
  qsi<-apply(x,2,quantile,p=p[2:3],na.rm=T)
  polygon(c(yrs,yrs[ny:1]),c(qs[1,],qs[2,ny:1]),border=NA,col='#b3ecff')

  polygon(c(yrs,yrs[ny:1]),c(qsi[1,],qsi[2,ny:1]),border=NA,col=qcol)
  if(!is.na(ablines[1]))abline(h=ablines,col='#99999980')

  if(addline)for(i in 1:2)lines(yrs,x[i,],col=lcol,lty=i)
  lines(yrs,apply(x,2,quantile,p=0.5,na.rm=T),lwd=2,col="white")
}



#' @name plot.MMSE
#' @title Standard plot for an object of class MMSE (multi MSE)
#'
#' @description Plot the projected biomass, fishing, mortality rate and yield for all stocks and MPs
#'
#' @param x Object of class \linkS4class{MMSE}. A Multi-OM object created by \code{multiMSE(MOM, ...)}
#' @param maxcol Integer. The maximum number of columns (MPs) to be plotted in each plot
#' @param qcol Character, color. The color of the inner percentile range
#' @param lcol Character, color. The color of the outer percentile range.
#' @param quants Numeric vector. The percentiles that are plotted (LB2,LB1,UB1,UB2). LB2 and UB2 are the outer percentiles, LB1 and UB1 are the inner percentiles.
#' @param curyr Integer. The current year from which projections start.
#' @param addline Logical. Should two individual simulations be added to the percentile plots?
#' @author T.Carruthers
#' @aliases plot,MMSE,missing-method
#' @importFrom grDevices rgb
#' @export
setMethod("plot", signature(x = "MMSE", y = "missing"),
          function(x, maxcol = 6, qcol = rgb(0.4, 0.8, 0.95), lcol = "dodgerblue4",
                   quants = c(0.05, 0.25, 0.75, 0.95), curyr = 2018, addline = FALSE) {
            MMSE <- x
            if(is.na(maxcol))maxcol=ceiling(length(MMSE@MPs)/0.5) # defaults to portrait 1:2
            MPs<-MMSE@MPs
            MPrefs<-MMSE@MPrefs
            nMPs<-length(MPrefs[,1,1])
            yrs<-curyr+(1:MMSE@proyears)
            ns<-MMSE@nstocks
            nf<-MMSE@nfleets

            plots<-split(1:nMPs, ceiling(seq_along(1:nMPs)/maxcol))



            # --- Biomass projection ---------------------------------------------------

            B_BMSY<-MMSE@B_BMSY
            Blims <- c(0,quantile(B_BMSY,0.95))

            for(pp in 1:length(plots)){

              toplot<-plots[[pp]]
              nt<-length(toplot)
              par(mfcol=c(ns,nt),mai=c(0.3,0.3,0.3,0.01),omi=c(0.4,0.5,0.05,0.05))

              for(MP in toplot){

                for(ss in 1:ns){
                  plot(range(yrs),Blims,col="white",yaxs="i")
                  plotquant(B_BMSY[,ss,MP,],p=quants,yrs,qcol,lcol,ablines=c(0.5,1),addline=addline)
                  mtext(paste(paste0("F",1:nf),MPrefs[MP,,ss],collapse=", "),3,line=0.2,font=2,cex=0.7)

                  if(MP==toplot[1])mtext(MMSE@Snames[ss],2,line=2.5)
                }
              }
            }

            mtext("Projection Year",1,font=2,outer=T,line=1.2)
            mtext("Biomass / BMSY",2,font=2,outer=T,line=1.9)


            # --- F projection -----------------------------------------------------------

            F_FMSY<-MMSE@F_FMSY
            F_FMSYsum<-apply(F_FMSY,c(1,2,4,5),sum,na.rm=T)
            Flims<- c(0,quantile(F_FMSYsum,0.95,na.rm=T))

            for(pp in 1:length(plots)){

              toplot<-plots[[pp]]
              nt<-length(toplot)
              par(mfcol=c(ns,nt),mai=c(0.3,0.3,0.3,0.01),omi=c(0.4,0.5,0.05,0.05))

              for(MP in toplot){

                for(ss in 1:ns){
                  plot(range(yrs),Flims,col="white",yaxs="i")
                  plotquant(F_FMSYsum[,ss,MP,],p=quants,yrs,qcol,lcol,ablines=c(0.5,1),addline=addline)
                  mtext(paste(paste0("F",1:nf),MPrefs[MP,,ss],collapse=", "),3,line=0.2,font=2,cex=0.7)

                  if(MP==toplot[1])mtext(MMSE@Snames[ss],2,line=2.5)
                }
              }
            }

            mtext("Projection Year",1,font=2,outer=T,line=1.2)
            mtext("Total (all fleets) F / FMSY",2,font=2,outer=T,line=1.9)


            # --- Total yield projection -----------------------------------------------------

            Yd<-MMSE@C#MMSE@OM$RefY
            Ydsum<-apply(Yd, c(1,2,4,5),sum,na.rm=T)
            Ydsum<-Ydsum/array(rep(Ydsum[,,,1],MMSE@proyears),dim(Ydsum))
            Ylims<- c(0,quantile(Ydsum,0.95,na.rm=T))

            for(pp in 1:length(plots)){

              toplot<-plots[[pp]]
              nt<-length(toplot)
              par(mfcol=c(ns,nt),mai=c(0.3,0.3,0.3,0.01),omi=c(0.4,0.5,0.05,0.05))

              for(MP in toplot){

                for(ss in 1:ns){

                  plot(range(yrs),Ylims,col="white",yaxs="i")
                  plotquant(Ydsum[,ss,MP,],p=quants,yrs,qcol,lcol,ablines=c(0.5,1),addline=addline)
                  mtext(paste(paste0("F",1:nf),MPrefs[MP,,ss],collapse=", "),3,line=0.2,font=2,cex=0.7)

                  if(MP==toplot[1])mtext(MMSE@Snames[ss],2,line=2.5)
                }
              }
            }

            mtext("Projection Year",1,font=2,outer=T,line=1.2)
            mtext(paste("Total yield (all fleets) relative to",curyr),2,font=2,outer=T,line=1.9)


            # --- Yield by fleet projection --------------------------------------------------

            YdbyF<-apply(Yd, c(2,3,4,5),mean,na.rm=T)
            cols<-rep(c('black','red','green','blue','orange','grey','purple'),3)
            ltys<-rep(1:3,each=7)


            for(pp in 1:length(plots)){

              toplot<-plots[[pp]]
              nt<-length(toplot)
              par(mfcol=c(ns,nt),mai=c(0.3,0.3,0.3,0.01),omi=c(0.4,0.5,0.05,0.05))

              for(MP in toplot){
                for(ss in 1:ns){
                  matplot(yrs,t(matrix(YdbyF[ss,,MP,],nrow=nf)),type='l',col=cols,lty=ltys,yaxs="i",ylim=c(0,max(YdbyF[ss,,MP,])))
                  mtext(paste(paste0("F",1:nf),MPrefs[MP,,ss],collapse=", "),3,line=0.2,font=2,cex=0.7)
                  if(MP==toplot[1])mtext(MMSE@Snames[ss],2,line=2.5)
                  if(MP==toplot[1])legend("bottomleft",legend=paste(paste0("F",1:nf),MMSE@Fnames[,ss]),text.col=cols,cex=0.8,text.font=2,bty='n')
                }
              }
            }

            mtext("Projection Year",1,font=2,outer=T,line=1.2)
            mtext("Expected yield by fleet",2,font=2,outer=T,line=1.9)


          })


#' A basic comparison of runMSE output (MSE) and multiMSE (MMSE)
#'
#' @param MSEsingle An object of class MSE arising from a run of runMSE(OM, ...)
#' @param MSEmulti An object of class MMSE arising from a run of multiMSE(MOM, ...)
#' @param p Integer. The stock number from the MSEmulti object (to be plotted)
#' @param f Integer. The fleet number from the MSEmulti object (to be plotted)
#' @param MPno Integer. The MP number from the MSEmulti and MSEsingle object (to be plotted)
#' @param maxsims Integer. The maximum number of simulations to plot.
#' @author T.Carruthers
#' @export
multidebug<-function(MSEsingle,MSEmulti,p=1,f=1,MPno=1,maxsims=4){

  MPno=1
  p=1
  f=1
  par(mfrow=c(2,2),mar=rep(2,4))

  SSBhs<-apply(slot(MSEsingle,"SSB_hist"),c(1,3),sum)
  SSBhm<-apply(slot(MSEmulti,"SSB_hist")[,p,,,],c(1,3),sum)
  matplot(t(SSBhs),type='l',main="SSBhist: runMSE")
  matplot(t(SSBhm),type='l',main="SSBhist: multiMSE")

  SSBs<-slot(MSEsingle,"SSB")
  SSBm<-slot(MSEmulti,"SSB")
  matplot(t(SSBs[,MPno,]),type='l',main="SSB: runMSE")
  matplot(t(SSBm[,p,MPno,]),type='l',main="SSB: multiMSE")

}


#' A basic SSB plot for debugging runMSE output
#'
#' @param MSEmulti An object of class MMSE arising from a run of multiMSE(MOM, ...)
#' @param maxsim Integer. The number of simulations to plot
#' @author T.Carruthers
#' @export
plotmulti<-function(MSEmulti,maxsim=8){

  np<-MSEmulti@nstocks
  nMPs<-MSEmulti@nMPs
  MPs<-MSEmulti@MPs
  Snames<-MSEmulti@Snames
  par(mfrow=c(np,nMPs),mai=c(0.7,0.7,0.1,0.1),omi=c(0.2,0.5,0.3,0.1))

  MSEmulti@MPrefs

  SSBm<-slot(MSEmulti,"SSB")
  nsim<-dim(SSBm)[1]
  maxsim<-min(nsim,maxsim)

  for(p in 1:np){
    for(m in 1:nMPs){
      matplot(t(SSBm[1:maxsim,p,m,]),type='l',main="",xlab="",ylab="")
      if(p==1)mtext(MPs[[p]][m],2,line=2.5,font=2)
      if(m==1)mtext(Snames[p],font=2,line=0.8)
    }
  }

  mtext("Projection Year",1,outer=T)

}