R/Retrospective/Likelihood/02_summaryplot_likelihood.r

In mkapur/ssdiags:

# source this code to make plots
# wanted fishery names may need to change - these are the fisheries we want to compare accross retrospective
# e.g., currently no cpue

library(r4ss)
library(colorRamps)

setwd(dirname.Retrospective)
getwd()

csv.filename <- paste('likelihood_by_fleet',format(Sys.time(), "%Y%m%d_%H%M.csv"),sep='_')
csv.yn <- T
plot.dir <-'plots_2'
plot.SScompare.yn <- T 
#wantedlikelabel <- c('Surv_like','Length_like','SizeFreq_like:_1','SizeFreq_like:_2')
#wantedlambdalabel <- c('Surv_lambda','Length_lambda','SizeFreq_lambda:_1','SizeFreq_lambda:_2')
#wantedfishery.names <- c('F2','F3','F8','F11','S1','S2','S3','S4','S5','S6')
# DC Edits Begin SMA 2/18/2019 Recruitment?
wantedlikelabel <- c('Surv_like','Length_like')                                   
wantedlambdalabel <- c('Surv_lambda','Length_lambda')                             
wantedfishery.names <- c('F1_EU_LL',      # Estimated selectivity
                         'F2_JPN_LL',     # Estimated selectivity
                         'F3_CTP_LL',     # Estimated selectivity
                         'F4_USA_LL',     # Estimated selectivity
                         'F5_VEN_LL',     # Estimated selectivity
                         #'F6_CAN_LL',    # Mirrored selectivity
                         #'F7_MOR_LL',    # Mirrored selectivity
                         #'F8_USA_RR',    # Mirrored selectivity
                         #'F9_BEL_LL',    # Mirrored selectivity
                         #'F10_MOR_PS',   # Mirrored selectivity
                         #'F11_CPR_LL',   # Mirrored selectivity
                         #'F12_OTH',      # Mirrored selectivity
                         'S1_USA_LL_Log', # Mirrored selectivity
                         'S2_USA_LL_Obs', # Mirrored selectivity
                         'S3_JPN_LL',     # Mirrored selectivity
                         'S4_EU_POR_LL',  # Mirrored selectivity
                         'S5_EU_ESP_LL',  # Mirrored selectivity
                         'S6_CTP_LL'      # Mirrored selectivity
                         )     

#wantedindexfleets <- c(9,10,11,12,13,14)
wantedindexfleets <- c(13,14,15,16,17,18)
#legendlabels <- c('2015','2014','2013','2012','2011','2010')
legendlabels <- c("./retro-1", "./retro-2", "./retro-3", "./retro-4", "./retro-5", "./retro0") 
#workdir <- './'
#workdir <- dir.work 

nummodels <- length(start.retro:-end.retro)
numfisheries <- length(wantedfishery.names)
numlike <- length(wantedlikelabel)
SSreps <- retroModels
summaryoutput <- retroSummary


if (plot.SScompare.yn) {
  if (!(plot.dir %in% dir())) {
    dir.create(plot.dir)
  }
  SSplotComparisons(summaryoutput, plot=T, pdf=T, png=F, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=endyrvec)
  SSplotComparisons(summaryoutput, plot=F, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=endyrvec)
  for (ii in 1:length(wantedindexfleets)) {
    #ii <- 1
    SSplotComparisons(summaryoutput, subplot=11, plot=F, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, indexfleets=wantedindexfleets[ii], endyrvec=endyrvec)
    SSplotComparisons(summaryoutput, subplot=12, plot=F, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, indexfleets=wantedindexfleets[ii], endyrvec=endyrvec)
  }
}


dirvec <- NULL
#dirvec <- c('./00_2015')
#dirvec <- c( "./retro-1" "./retro-2" "./retro-3" "./retro-4" "./retro-5" "./retro0" )
#dirvec <- c()

if (is.null(dirvec)) {
  setwd(paste(dirname.Retrospective,"retrospectives",sep="\\"))
  dirvec <- list.dirs()
  dirvec <- dirvec[-1]
  setwd("../")
}

totlike <- summaryoutput$likelihoods[summaryoutput$likelihoods$Label == 'TOTAL',1:nummodels]
outdata <- as.data.frame(matrix(data=NA, nrow=nummodels, ncol=2+numfisheries*numlike*2))
colnames(outdata)[1] <- 'model'
colnames(outdata)[2] <- 'totlike' 


for (ii in 1:nummodels) {
  #ii <- 1 
  wantedrows.model <- summaryoutput$likelihoods_by_fleet$model == ii
  outdata$model[ii] <- dirvec[ii]
  outdata$totlike[ii] <- totlike[ii]
  outcol <- 3 # init colnum
  for (jj in 1:numfisheries) {
    #jj<-1
    wantedcols.fleet <- names(summaryoutput$likelihoods_by_fleet) == wantedfishery.names[jj]
    if (sum(wantedcols.fleet) == 0) {
      stop(paste(wantedfishery.names[jj],'does not match any fleetname'))
    }
    for (kk in 1:numlike) {
      #kk<-1
      wantedrows.like <- wantedrows.model & (summaryoutput$likelihoods_by_fleet$Label == wantedlikelabel[kk])
      wantedrows.lambda <- wantedrows.model & (summaryoutput$likelihoods_by_fleet$Label == wantedlambdalabel[kk])
      like <- summaryoutput$likelihoods_by_fleet[wantedrows.like,wantedcols.fleet]
      lambda <- summaryoutput$likelihoods_by_fleet[wantedrows.lambda,wantedcols.fleet]
      outdata[ii,outcol] <- as.numeric(lambda)
      outdata[ii,outcol+1] <- as.numeric(like)
      if (ii == 1) {
        colnames(outdata)[outcol] <- paste(wantedfishery.names[jj],wantedlambdalabel[kk])
        colnames(outdata)[outcol+1] <- paste(wantedfishery.names[jj],wantedlikelabel[kk])
      }
      outcol <- outcol + 2
    }
  }
}


out.lambdacols <- grep(pattern='lambda',x=names(outdata))
out.goodlambdacols <- out.lambdacols[which(colSums(outdata[,out.lambdacols], na.rm=T) > 0)]
out2 <- as.matrix(outdata[,sort(c(1,2,out.goodlambdacols,out.goodlambdacols+1))])

if (csv.yn) {
  setwd(plot.dir)
  write.csv(out2, file=csv.filename, row.names=F)
  setwd(dirname.Retrospective)
}
graphics.off()