LOPlot<-function(){
    if(!(file.exists("../LO/LOrun2.RData") || file.exists("../LO/LOrun3.RData") || file.exists("../LO/LOrun4.RData") || file.exists("../LO/LOrun5.RData") )) return(NULL);
    for(i in 2:max(obs$fleet)) try(load(paste0("../LO/LOrun",i,".RData")))
    
    
    cols=c("red","green","cyan","orange","deeppink","darkorchid1","green")

    LOOK <- c(exists("LOrun2") && LOrun2$opt$convergence==0,
              exists("LOrun3") && LOrun3$opt$convergence==0,
              exists("LOrun4") && LOrun4$opt$convergence==0,
              exists("LOrun5") && LOrun5$opt$convergence==0,
              exists("LOrun6") && LOrun6$opt$convergence==0,
              exists("LOrun7") && LOrun7$opt$convergence==0)
    par(mfrow=c(3,1))

    time=data$times
    uncertaintyPlot(time, SSB/1000, SSB.lo/1000, SSB.hi/1000, xlab="Time", ylab="SSB", type="b", cex=.5, las=1)
    if(LOOK[1]) lines(time[1:length(LOrun2[[2]])], LOrun2[[2]]/1000,col=cols[2],lwd=1.5)
    if(LOOK[2]) lines(time[1:length(LOrun3[[2]])], LOrun3[[2]]/1000,col=cols[3],lwd=1.5)
    if(LOOK[3]) lines(time[1:length(LOrun4[[2]])], LOrun4[[2]]/1000,col=cols[4],lwd=1.5)
    if(LOOK[4]) lines(time[1:length(LOrun5[[2]])], LOrun5[[2]]/1000,col=cols[5],lwd=1.5)
    if(LOOK[5]) lines(time[1:length(LOrun6[[2]])], LOrun6[[2]]/1000,col=cols[6],lwd=1.5)
    if(LOOK[6]) lines(time[1:length(LOrun7[[2]])], LOrun7[[2]]/1000,col=cols[7],lwd=1.5)

    legend("topright",legend=paste("w/o",2:max(obs$fleet))[LOOK],col=cols[2:max(obs$fleet)][LOOK],lwd=1.5)
    
    uncertaintyPlot(time, FBAR, FBAR.lo, FBAR.hi, xlab="Time", ylab="FBAR", type="b", cex=.5, las=1)
    if(LOOK[1]) lines(time[1:length(LOrun2[[2]])], LOrun2[[3]],col=cols[2],lwd=1.5)
    if(LOOK[2]) lines(time[1:length(LOrun3[[2]])], LOrun3[[3]],col=cols[3],lwd=1.5)
    if(LOOK[3]) lines(time[1:length(LOrun4[[2]])], LOrun4[[3]],col=cols[4],lwd=1.5)
    if(LOOK[4]) lines(time[1:length(LOrun5[[2]])], LOrun5[[3]],col=cols[5],lwd=1.5)
    if(LOOK[5]) lines(time[1:length(LOrun6[[2]])], LOrun6[[3]],col=cols[6],lwd=1.5)
    if(LOOK[6]) lines(time[1:length(LOrun7[[2]])], LOrun7[[3]],col=cols[7],lwd=1.5)

    legend("topright",legend=paste("w/o",2:max(obs$fleet))[LOOK],col=cols[2:max(obs$fleet)][LOOK],lwd=1.5)
    
    sel = seq(3,ncol(pl$logN),by=4)
    time = data$time[sel]-0.5
    logrecr = pl$logN[1,sel]
    logrecr.sd = plsd$logN[1,sel]
    logrecr.lo = logrecr - 2*logrecr.sd
    logrecr.hi = logrecr + 2*logrecr.sd
    uncertaintyPlot(time, exp(logrecr), exp(logrecr.lo), exp(logrecr.hi), xlab="Time", ylab="Recruitment", type="b", cex=.5, las=1)
    if(LOOK[1]) lines(time[1:length(LOrun2[[2]])],exp(LOrun2[[4]]$logR[sel][1:length(LOrun2[[2]])]), type="l",col=cols[2],lwd=1.5)
    if(LOOK[2]) lines(time[1:length(LOrun3[[2]])],exp(LOrun3[[4]]$logR[sel][1:length(LOrun3[[2]])]), type="l",col=cols[3],lwd=1.5)
    if(LOOK[3]) lines(time[1:length(LOrun4[[2]])],exp(LOrun4[[4]]$logR[sel][1:length(LOrun4[[2]])]), type="l",col=cols[4],lwd=1.5)
    if(LOOK[4]) lines(time[1:length(LOrun5[[2]])],exp(LOrun5[[4]]$logR[sel][1:length(LOrun5[[2]])]), type="l",col=cols[5],lwd=1.5)
    if(LOOK[5]) lines(time[1:length(LOrun6[[2]])],exp(LOrun6[[4]]$logR[sel][1:length(LOrun6[[2]])]), type="l",col=cols[6],lwd=1.5)
    if(LOOK[6]) lines(time[1:length(LOrun7[[2]])],exp(LOrun7[[4]]$logR[sel][1:length(LOrun7[[2]])]), type="l",col=cols[7],lwd=1.5)

    legend("topright",legend=paste("w/o",2:max(obs$fleet))[LOOK],col=cols[2:max(obs$fleet)][LOOK],lwd=1.5)
    
}

setupMohnDF<-function(base, peel1, peel2, peel3, peel4, peel5, quarter){
    sel <- seq(quarter,ncol(pl$logN),by=4) 
    df <- data.frame( base = base[sel],
                     peel1 = peel1[sel], 
                     peel2 = peel2[sel],
                     peel3 = peel3[sel],
                     peel4 = peel4[sel],
                     peel5 = peel5[sel]
                     ,row.names=data$times[sel]
                     )
    df
}

retroPlot<-function(){
    if(!(file.exists("../RETRO/retrorun1.RData") || file.exists("../RETRO/retrorun2.RData") || file.exists("../RETRO/retrorun3.RData") || file.exists("../RETRO/retrorun4.RData") || file.exists("../RETRO/retrorun5.RData") )) return(NULL);

    for(i in 1:5) try(load(paste0("../RETRO/retrorun",i,".RData")))

     RETROOK <- c(exists("retrorun1") && retrorun1$opt$convergence==0,
              exists("retrorun2") && retrorun2$opt$convergence==0,
              exists("retrorun3") && retrorun3$opt$convergence==0,
              exists("retrorun4") && retrorun4$opt$convergence==0,
              exists("retrorun5") && retrorun5$opt$convergence==0)
    
    
    
    cols=c("red","magenta","brown","cyan","orange","deeppink","darkorchid1")
    par(mfrow=c(3,1))
    
    time=data$times
    uncertaintyPlot(time, SSB/1000, SSB.lo/1000, SSB.hi/1000, xlab="Time", ylab="SSB", type="b", cex=.5, las=1)
    if(RETROOK[1]) lines(time[1:length(retrorun1[[2]])], retrorun1[[2]]/1000,col=cols[1])
    if(RETROOK[2]) lines(time[1:length(retrorun2[[2]])], retrorun2[[2]]/1000,col=cols[2])
    if(RETROOK[3]) lines(time[1:length(retrorun3[[2]])], retrorun3[[2]]/1000,col=cols[3])
    if(RETROOK[4]) lines(time[1:length(retrorun4[[2]])], retrorun4[[2]]/1000,col=cols[4])
    if(RETROOK[5]) lines(time[1:length(retrorun5[[2]])], retrorun5[[2]]/1000,col=cols[5])

    ssbdf <- setupMohnDF(SSB/1000, retrorun1[[2]]/1000, retrorun2[[2]]/1000, retrorun3[[2]]/1000,retrorun4[[2]]/1000,retrorun5[[2]]/1000,quarter=3)
    mrho.ssb  <- icesAdvice::mohn(ssbdf)
    legend("topright", legend=paste0("rho = ",round(mrho.ssb*100),"%"))

    uncertaintyPlot(time, FBAR, FBAR.lo, FBAR.hi, xlab="Time", ylab="FBAR", type="b", cex=.5, las=1)
    if(RETROOK[1]) lines(time[1:length(retrorun1[[2]])], retrorun1[[3]],col=cols[1])
    if(RETROOK[2]) lines(time[1:length(retrorun2[[2]])], retrorun2[[3]],col=cols[2])
    if(RETROOK[3]) lines(time[1:length(retrorun3[[2]])], retrorun3[[3]],col=cols[3])
    if(RETROOK[4]) lines(time[1:length(retrorun4[[2]])], retrorun4[[3]],col=cols[4])
    if(RETROOK[5]) lines(time[1:length(retrorun5[[2]])], retrorun5[[3]],col=cols[5])

    
    fbardf <- setupMohnDF(FBAR, retrorun1[[3]], retrorun2[[3]], retrorun3[[3]],retrorun4[[3]],retrorun5[[3]],quarter=3)
    mrho.f  <- icesAdvice::mohn(fbardf)
    legend("topright", legend=paste0("rho = ",round(mrho.f*100),"%"))

    sel = seq(3,ncol(pl$logN),by=4)
    time = data$time[sel]-0.5
    logrecr = pl$logN[1,sel]
    logrecr.sd = plsd$logN[1,sel]
    logrecr.lo = logrecr - 2*logrecr.sd
    logrecr.hi = logrecr + 2*logrecr.sd
    uncertaintyPlot(time, exp(logrecr), exp(logrecr.lo), exp(logrecr.hi), xlab="Time", ylab="Recruitment", type="b", cex=.5, las=1)
    if(RETROOK[1]) lines(time[1:length(retrorun1[[2]])],exp(retrorun1[[4]]$logR[sel][1:length(retrorun1[[2]])]), type="l",col=cols[1])
    if(RETROOK[2]) lines(time[1:length(retrorun2[[2]])],exp(retrorun2[[4]]$logR[sel][1:length(retrorun2[[2]])]), type="l",col=cols[2])
    if(RETROOK[3]) lines(time[1:length(retrorun3[[2]])],exp(retrorun3[[4]]$logR[sel][1:length(retrorun3[[2]])]), type="l",col=cols[3])
    if(RETROOK[4]) lines(time[1:length(retrorun4[[2]])],exp(retrorun4[[4]]$logR[sel][1:length(retrorun4[[2]])]), type="l",col=cols[4])
    if(RETROOK[5]) lines(time[1:length(retrorun5[[2]])],exp(retrorun5[[4]]$logR[sel][1:length(retrorun5[[2]])]), type="l",col=cols[5])

    recdf<- setupMohnDF(exp(pl$logN[1,]), exp(retrorun1[[4]]$logR), exp(retrorun2[[4]]$logR), exp(retrorun3[[4]]$logR),exp(retrorun4[[4]]$logR),exp(retrorun5[[4]]$logR),quarter=3)
    mrho.r  <- icesAdvice::mohn(recdf)
    legend("topright", legend=paste0("rho = ",round(mrho.r*100),"%"))
    return(list(mrho.ssb=mrho.ssb,mrho.f=mrho.f,mrho.r=mrho.r))
}