## Retro run script for SESAM
for(i in 1:5){
    RETROnumber=i
    is.RETRO=TRUE
    source("src/sesam.R")
}
is.RETRO=FALSE
source("src/sesam.R")
setwd("..")
for(i in 1:5) load(paste0("RETRO/retrorun",i,".RData"))
cols=rainbow(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(){
    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)
    lines(time[1:length(retrorun1[[2]])], retrorun1[[2]]/1000,col=cols[1])
    lines(time[1:length(retrorun2[[2]])], retrorun2[[2]]/1000,col=cols[2])
    lines(time[1:length(retrorun3[[2]])], retrorun3[[2]]/1000,col=cols[3])
    lines(time[1:length(retrorun4[[2]])], retrorun4[[2]]/1000,col=cols[4])
    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  <- icesAdvice::mohn(ssbdf)
    legend("topright", legend=paste0("rho = ",round(mrho*100),"%"))

    uncertaintyPlot(time, FBAR, FBAR.lo, FBAR.hi, xlab="Time", ylab="FBAR", type="b", cex=.5, las=1)
    lines(time[1:length(retrorun1[[2]])], retrorun1[[3]],col=cols[1])
    lines(time[1:length(retrorun2[[2]])], retrorun2[[3]],col=cols[2])
    lines(time[1:length(retrorun3[[2]])], retrorun3[[3]],col=cols[3])
    lines(time[1:length(retrorun4[[2]])], retrorun4[[3]],col=cols[4])
    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  <- icesAdvice::mohn(fbardf)
    legend("topright", legend=paste0("rho = ",round(mrho*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)
    lines(time[1:length(retrorun1[[2]])],exp(retrorun1[[4]]$logR[sel][1:length(retrorun1[[2]])]), type="l",col=cols[1])
    lines(time[1:length(retrorun2[[2]])],exp(retrorun2[[4]]$logR[sel][1:length(retrorun2[[2]])]), type="l",col=cols[2])
    lines(time[1:length(retrorun3[[2]])],exp(retrorun3[[4]]$logR[sel][1:length(retrorun3[[2]])]), type="l",col=cols[3])
    lines(time[1:length(retrorun4[[2]])],exp(retrorun4[[4]]$logR[sel][1:length(retrorun4[[2]])]), type="l",col=cols[4])
    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  <- icesAdvice::mohn(recdf)
    legend("topright", legend=paste0("rho = ",round(mrho*100),"%"))
    
}




setwd("res")
stamp=paste0(stamp)
png(filename = paste(stamp,"_retro%03d.png", sep=''), width = 480, height = 480,
    units = "px", pointsize = 10, bg = "white")
retroPlot()
dev.off()


png(filename = paste("big_",stamp,"_retro%03d.png", sep=''), width = 1200, height = 1200, 
    units = "px", pointsize = 20, bg = "white")
retroPlot()     
dev.off()