## Note, this plotscript assumes that sesam.R has been run already in the same R session.
source("src/extraPlots.R")
####################
## OSA residuals
####################
epsfleet = data$eps[data$fleet]
if(doOSA && !((exists("is.RETRO") && is.RETRO)) || ((exists("is.LO") && is.LO)) ){
    ## condition on these first data points to make OSA work
    mycond=1:12
    set.seed(123)
    system.time(OSA<-oneStepPredict(obj,"logobs","keep",discrete=FALSE,parallel=FALSE,subset=(max(mycond)+1):length(data$obs),conditional=mycond)) 
    
    ## obs less than eps need special residual calc and randomization
    OSA$mnll = diff(c(0,OSA$nll))
    OSA$Px = exp( - OSA$mnll)
    
    OSA$isLTeps = !is.na(data$obs[-c(mycond)]) & data$obs[-c(mycond)]<epsfleet[-c(mycond)] 
    
    OSA$Px[OSA$isLTeps][ OSA$Px[OSA$isLTeps] > 1 ] = 1  
    nr = sum(OSA$isLTeps)
    rand = runif(nr,max=OSA$Px[OSA$isLTeps])
    OSA$residual[ OSA$isLTeps ] = qnorm(rand)
}

AIC.sesam <- function(x) 2*x$objective + 2*length(x$par)

##save(data,param,pl,plsd,rep,jointrep,sdrep,opt,OSArep, file="run/samtmb.RData")
#################################
## Survival process residuals
#################################
rmvnorm <- function(n = 1, mu, Sigma) {
        p <- length(mu)
        if (!all(dim(Sigma) == c(p, p))) {
            stop("incompatible arguments")
        }
        idx <- diag(Sigma) > .Machine$double.xmin
        L <- matrix(0, p, p)
        if (any(idx)) {
            L[idx, idx] <- chol(Sigma[idx, idx])
        }
        X <- matrix(rnorm(p * n), n)
        X <- drop(mu) + t(X %*% L)
        if (n == 1) {
            drop(X)
        }
        else {
            t(X)
        }
}

idx <- which(names(sdrep$value) == "residN")
    residN <- rmvnorm(1, mu = sdrep$value[idx], Sigma = sdrep$cov[idx, 
        idx])
residNmat <- matrix(residN,ncol=4)
sel0 <- which(residNmat[,1]==0)
residNmat[sel0,1] = NA
residNmat <- residNmat[-1,]
rownames(residNmat)<-data$times[-1]
colnames(residNmat)<-c(0:3)

residQ = (as.numeric(rownames(residNmat)) - floor(as.numeric(rownames(residNmat))))*4 +1

process.bias.p <- rbind(apply(residNmat[residQ==1,],2,function(x) ifelse(any(!is.na(x)),t.test(na.omit(x))$p.value,NA)),
apply(residNmat[residQ==2,],2,function(x) ifelse(any(!is.na(x)),t.test(na.omit(x))$p.value,NA)),
apply(residNmat[residQ==3,],2,function(x) ifelse(any(!is.na(x)),t.test(na.omit(x))$p.value,NA)),
apply(residNmat[residQ==4,],2,function(x) ifelse(any(!is.na(x)),t.test(na.omit(x))$p.value,NA)))

rownames(process.bias.p)<-1:4

## Stock-recruitment
ssb = sum[rownames(sum)=="ssb",]
ssb.lo = ssb[,1] - 2*ssb[,2]
ssb.hi = ssb[,1] + 2*ssb[,2]
selq1 = seq(1,ncol(pl$logN),by=4)
selq1 = selq1 [ selq1 < max(sel) ] 
selq3 = seq(3,ncol(pl$logN),by=4)

#####################
## Plotting
#####################


plots<-function(){


    ## SSB
    par(mfrow=c(1,1), mar=c(4,4,.5,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(xsa$Year,xsa$SSB/1000,col=3,lwd=2)

    ## TSB
    uncertaintyPlot(time, TSB/1000, TSB.lo/1000, TSB.hi/1000, xlab="Time", ylab="TSB", type="b", cex=.5, las=1)


    ##Fbar
    uncertaintyPlot(time, FBAR, FBAR.lo, FBAR.hi, xlab="Time", ylab="FBAR", type="b", cex=.5, las=1)
    #lines(xsa$Year,xsa$Fbar,col=3,lwd=2)


    ## Fbar yearly
    minY = min(floor(time))
    time.y=seq(minY,minY+length(fbar.y)-1,by=1)
    lines(time.y,fbar.y,col="cyan",lwd=2)
    lines(time.y,fbar.y.lo,col="cyan",lwd=2,lty=2)
    lines(time.y,fbar.y.hi,col="cyan",lwd=2,lty=2)


    bpAdd<-function(x,y,v, scale=3, ...){
        points(x,y,cex=sqrt(abs(v))*scale, ...)
    }
    
    bplegendAdd<-function(scale=3){
      zscale <- c(1,2,3,4)
      xx <- 1:length(zscale)
      yy <- rep(1, length(zscale))
      plot(xx,yy,cex=sqrt(zscale)*scale, pch=1, col=1,axes=FALSE,xlim=c(0,max(xx)+1),xlab="",ylab="")
      text(xx,yy,labels=zscale)
    }

    ## OSA residuals
    par(mfrow=n2mfrow(max(obs$fleet)+1))
    if(doOSA){
        OSA$fleet = obs$fleet[-c(mycond)]
        OSA$year = obs$t1[-c(mycond)]
        OSA$age = obs$ageFrom[-c(mycond)]
        OSA$quarter = ((as.numeric(OSA$year) - floor(as.numeric(OSA$year)))*4)+1
        for(f in sort(unique(obs$fleet))){
            sel = which(data$fleet==f)
            sel=setdiff(sel,mycond)
            selOSA = which(OSA$fleet==f)
            ##browser()
            bp(data$t1[sel],data$ageFrom[sel],OSA$residual[selOSA]*(!is.na(data$obs[sel]) & (data$obs[sel]>epsfleet[sel])),scale=3,main=paste0("OSA fleet ",f),ylim=c(-1,4))
            v = OSA$residual[selOSA]*(!is.na(data$obs[sel]) & (data$obs[sel]<epsfleet[sel]))
            bpAdd(data$t1[sel],data$ageFrom[sel],v,scale=3,pch=3,col=ifelse(v<0,'tomato2','blue'))
        }
        bplegendAdd()
        par(mfrow=n2mfrow(max(obs$fleet)))
    }
    
    ## full conditional residuals
    for(f in sort(unique(obs$fleet))){
        sel = which(obs$fleet==f) 
        bp(data$t1[sel],data$ageFrom[sel],(!is.na(data$obs[sel]) & (data$obs[sel]>epsfleet[sel]))*(data$logobs[sel]-rep$predObs[sel])/rep$predSd[sel],scale=3,main=paste0("full cond resid f= ",f),ylim=c(-1,4))
    }
    
    
    if(doOSA){
        
        
        restab <- xtabs(residual ~ year + age + fleet,data=OSA)
        restab[ restab==0 ] <- NA

        restab.q <- xtabs(residual ~ year + age + quarter,data=subset(OSA,fleet==1))
        restab.q[ restab.q==0 ] <- NA
        
        
        var.interval<- function(data, conf.level = 0.95) {
            df = length(data) - 1
            chilower = qchisq((1 - conf.level)/2, df)
            chiupper = qchisq((1 - conf.level)/2, df, lower.tail = FALSE)
            v = var(data)
            c(df * v/chiupper, df * v/chilower)
        }

        mean.interval<-function(data, conf.level = 0.95) {
            n = length(data);
            mu = mean(data);
            s = sd(data);
            error = qt(0.975,df=n-1)*s/sqrt(n)
            c(mu-error, mu+error)
        }

        calc.limval <- function(p, acf) qnorm((2-p)/2)/sqrt(acf$n.used)

        for(f in unique(obs$fleet)){
            cat("======== fleet ",f," ==========\n")
            tab=xtabs(residual ~ age + year,data=subset(OSA,fleet==f)) 
            cat("mean intv:\n")
            mi = apply(tab,1,mean.interval)
            print(mi)
            print(mi[1,]<0 & mi[2,]>0)
            cat("sd intv:\n")
            si = apply(tab,1,var.interval)
            print(si)
            print(si[1,]<1 & si[2,]>1)
            cat("acf test:\n");
            acftab = apply(tab,1,acf,plot=FALSE,lag.max=5)
            acfok = lapply(acftab,function(x) !any(abs(x$acf[-1])>calc.limval(0.05,x)))
            print(unlist(acfok))
            
        }

        
        ages <- sort(unique(obs$ageFrom))
        fleets <- sort(unique(obs$fleet))
        quarters <- sort(unique(OSA$quarter))
        bias <- variance <- bias.p <- variance.p <- matrix(NA,nrow=dim(restab)[3],ncol=max(obs$ageFrom) - min(obs$ageFrom) + 1,dimnames=list(Fleet=fleets,Age=ages))

        bias.q <- variance.q <- bias.q.p <- variance.q.p <- matrix(NA,nrow=dim(restab.q)[3],ncol=max(obs$ageFrom) - min(obs$ageFrom) + 1,dimnames=list(Quarter=quarters,Age=ages))

        chisqtest<-function(x,testvar=1){
            n <- length(x)
            s2 <- var(x)
            chsq <- (n-1)*s2/testvar
            pval <- 1 - pchisq(chsq,n-1)
            pval
        }
    

        for(fl in fleets){
            for(a in ages){
                colu <- a+1
                if(!all(is.na(restab[,colu,fl]))){
                    
                    bias[fl,colu] <- mean(restab[,colu,fl],na.rm=TRUE)
                    variance[fl,colu] <- var(na.omit(restab[,colu,fl]))
                    bias.p[fl,colu] <- t.test(restab[,colu,fl])$p.value
                    variance.p[fl,colu] <- chisqtest(na.omit(restab[,colu,fl]))
                }
            }
        }
        for(qq in quarters){
            for(a in ages){
                colu <- a+1
                if(!all(is.na(restab.q[,colu,qq]))){
                    
                    bias.q[qq,colu] <- mean(restab.q[,colu,qq],na.rm=TRUE)
                    variance.q[qq,colu] <- var(na.omit(restab.q[,colu,qq]))
                    bias.q.p[qq,colu] <- t.test(na.omit(restab.q[,colu,qq]))$p.value
                    variance.q.p[qq,colu] <- chisqtest(na.omit(restab.q[,colu,qq]))
                }
            }
        }
        
        
        library(plot.matrix)
        palette=c("#273871", "#5087C1", "#ACCCE4", "#F4FAFE") ## hcl.colors(4,palette="Blues")
        p.breaks=c(0,0.001,0.05,0.1,1)
        op <- par(mfrow=n2mfrow(2),las=1,mar=c(5,8,5,5))
        plot(apply(t(bias.p),2,rev),breaks=p.breaks,col=palette,fmt.cell = "%.4f",xlab="Fleet",ylab="Age",main="Bias")
        plot(apply(t(variance.p),2,rev),breaks=p.breaks,col=palette,fmt.cell = "%.4f",xlab="Fleet",ylab="Age",main="Variance",key=NULL)

        ##plot(apply(t(bias.q.p),2,rev),breaks=p.breaks,col=palette,fmt.cell = "%.4f",xlab="Quarter",ylab="Age",main="Bias")
        ##plot(apply(t(variance.q.p),2,rev),breaks=p.breaks,col=palette,fmt.cell = "%.4f",xlab="Quarter",ylab="Age",main="Variance",key=NULL)


        plot(process.bias.p,breaks=p.breaks,col=palette,fmt.cell = "%.4f",xlab="Age",ylab="Quarter",main="Survival process bias")
        
        par(op)
        ## residual corr
        par(mfrow=n2mfrow(max(obs$fleet)))
        ccolors <- c("#A50F15","#DE2D26","#FB6A4A","#FCAE91","#FEE5D9","white",
                     "#EFF3FF","#BDD7E7","#6BAED6","#3182BD","#08519C")
        for(f in sort(unique(obs$fleet))){sel = which(obs$fleet==f)
            tab=xtabs(OSA$residual[sel] ~ data$t1[sel] + data$ageFrom[sel])
            tcor=cor(tab,use="pairwise.complete.obs")
            plotcorr(tcor,col=ccolors[5*tcor+6],main=paste("Fleet",f))
        }  

    }

    ## F by quarter
    par(mfrow=c(4,4))
    for(aa in 1:4){
        for(qq in 1:4){
            if(!(aa==1 && qq<3)){
            myf = pl$logF[aa,seq(qq,ncol(pl$logF),by=4)]
            sdmyf = plsd$logF[aa,seq(qq,ncol(pl$logF),by=4)]
            plot(data$times[seq(qq,ncol(pl$logF),by=4)],myf,main=paste0("age ",aa-1," q",qq),type="b",ylab="logF",xlab="time")
            lines(data$times[seq(qq,ncol(pl$logF),by=4)],myf-2*sdmyf,lty=2)
            lines(data$times[seq(qq,ncol(pl$logF),by=4)],myf+2*sdmyf,lty=2)
            } else {
                plot(1,1,type="n")
            }
        }}

    ## Time-steps with large jumps in F
    par(mfrow=c(1,1))
    plot(data$times,jump,col=bigjump+1,pch=16)
    abline(h=JUMPLIMIT)


    ## Recruitment 
    par(mfrow=c(1,1))
    sel = seq(3,ncol(pl$logN),by=4)
    logrecr = pl$logN[1,sel]
    logrecr.sd = plsd$logN[1,sel]
    logrecr.lo = logrecr - 2*logrecr.sd
    logrecr.hi = logrecr + 2*logrecr.sd
    uncertaintyPlot(data$time[sel]-0.5, exp(logrecr), exp(logrecr.lo), exp(logrecr.hi), xlab="Time", ylab="Recruitment", type="b", cex=.5, las=1)
    
    ## Stock-recruitment
    ssb = sum[rownames(sum)=="ssb",]
    ssb.lo = ssb[,1] - 2*ssb[,2]
    ssb.hi = ssb[,1] + 2*ssb[,2]
    selq1 = seq(1,ncol(pl$logN),by=4)
    selq1 = selq1 [ selq1 < max(sel) ] 
    selq3 = seq(3,ncol(pl$logN),by=4)
        
    cols = colorRampPalette(c("blue", "red"))( length(selq1) )

    q1ssb = ssb[selq1] 
    q1times = data$times[selq1]
    
    plot(q1ssb,exp(logrecr),ylim=c(0,max(exp(logrecr))),xlim=c(0,max(q1ssb)),xlab="SSB Q1",ylab="recruitment",pch=16,col=cols)
    disply = 0.02*max(exp(logrecr))
    displx = 0.02*max(q1ssb)
    text(q1ssb+displx,exp(logrecr)+disply,labels=sapply(as.character(data$times[selq1]),substr,start=3,stop=4 ))
    medrec = median(exp(logrecr))
    yabovei = exp(logrecr)>medrec & q1times<2025
    yabove = q1times[yabovei]
    ##blim1 = min(q1ssb[yabovei] )
    ##blim1y = data$times[selq1][which.min( 
    ##blim3 = mean( head( sort(q1ssb[yabovei], decreasing = FALSE ), 3 ))
    abline(h=medrec,lty=2)
    ##abline(v=blim1)
    ##abline(v=blim3)
    ##cat("blim1:", blim1, " blim3:",blim3,"\n")
    
    plot(ssb[selq3],exp(logrecr),ylim=c(0,max(exp(logrecr))),xlim=c(0,max(ssb[selq3])),xlab="SSB Q3",ylab="recruitment",pch=16,col=cols)
    disply = 0.02*max(exp(logrecr))
    displx = 0.02*max(ssb[selq3])
    text(ssb[selq3]+displx,exp(logrecr)+disply,labels=sapply(as.character(data$times[selq3]),substr,start=3,stop=4 ))


    ## observed vs predicted
    par(mfrow=c(3,3))
    for(f in sort(unique(obs$fleet))){
        sel = which(obs$fleet==f & data$obs>data$eps[f]) 
        plot( data$obs[sel], exp(rep$predObs[sel]),main=paste0("fleet ",f),xlab="observed",ylab="predicted")
        abline(0,1)
    }

    par(mfrow=c(max(obs$fleet),4))
    for(f in sort(unique(obs$fleet))){
        for(aa in sort(unique(obs$ageFrom))){
            sel = which(obs$fleet==f & data$obs>data$eps[f] & obs$ageFrom==aa)
            if(length(sel)>0){
                plot(obs$t1[sel],log(data$obs[sel]),main=paste0("fleet ",f," age ",aa),xlab="Time",ylab="log obs")
                lines(obs$t1[sel],rep$predObs[sel])
            } else {
                plot(1,1,type="n",axes=FALSE)
            }
        }
    }

    ## Total catch obs versus pred
    par(mfrow=c(2,1))
    predtotcatch<<-numeric(length(unique(data$t1)))
    totcatch<<-predtotcatch
    cc=0
    for(tt in unique(data$t1)){
        cc=cc+1
        sel = which(data$fleet==1 & data$t1==tt)
        auxsel = which(data$auxt1==tt)
        predtotcatch[cc] <<- sum(exp(rep$predObs[sel])*data$auxCW[auxsel]);
        totcatch[cc] <<- sum(data$obs[sel]*data$auxCW[auxsel]);
    }
    plot(unique(data$t1),totcatch,ylim=range(c(totcatch,predtotcatch)),xlab="time",ylab="total catch in weight")
    lines(unique(data$t1),predtotcatch,pch=3,type="b")
    legend("topright",pch=c(1,3),legend=c("observed","predicted"),lty=c(NA,1))


    ## Do by year
    totcatchy = aggregate(totcatch,by=list(floor(unique(data$t1))),FUN=sum)
    predtotcatchy = aggregate(predtotcatch,by=list(floor(unique(data$t1))),FUN=sum)
    plot(unique(floor(data$t1)),totcatchy$x,xlab="time",ylab="total catch in weight",ylim=range(c(totcatchy,predtotcatchy)))
    legend("topright",pch=c(1,3),legend=c("observed","predicted"),lty=c(NA,1))
    lines(unique(floor(data$t1)),predtotcatchy$x,pch=3,type="b")

    ## Q plot
    par(mfrow=c(1,1))
    qsel = grep("logQ",names(sdrep$par.fixed))
    qsds = sqrt(diag(sdrep$cov.fixed)[qsel])
    qs = sdrep$par.fixed[qsel]
    qlo = qs - 2*qsds
    qhi = qs + 2*qsds
    yl = range(qlo,qhi)
    plot(which(!is.na(data$keyLogQ[2,]))-1,qs[na.omit(data$keyLogQ[2,])+1],ylim=yl,xlim=c(0,3.5),xlab="age",ylab="logQ",cex=1.5,pch=16)
    arrows(which(!is.na(data$keyLogQ[2,]))-1, qlo[na.omit(data$keyLogQ[2,])+1], y1 = qhi[na.omit(data$keyLogQ[2,])+1], angle = 90, code = 3, length = 0.1,lwd=1.5)
    for(ff in 3:max(data$fleet)){
        points(which(!is.na(data$keyLogQ[ff,]))-1+0.05*(ff-2),qs[na.omit(data$keyLogQ[ff,])+1],ylim=yl,col=ff-1,cex=1.5,pch=16)
        arrows(which(!is.na(data$keyLogQ[ff,]))-1+0.05*(ff-2), qlo[na.omit(data$keyLogQ[ff,])+1], y1 = qhi[na.omit(data$keyLogQ[ff,])+1], angle = 90, code = 3, length = 0.1,col=ff-1,lwd=1.5)
    }
    legend("bottomright",legend=paste("fleet",2:6),col=1:6,pch=16)
    
}## end plots    

############################
    ## Forecasting
############################

do.forecast<-function(mylags.forward,myseed,mynosim,PM,SW,NM,ssbmin,Fscaling=NULL,oldstyle=TRUE,prob=0.05,Myearrange=NULL){
    tsel=which(data$times<=max(data$t1))

    sim.tac = forecast.sesam(obj,jointrep,data,stepsback=0,lags.forward=mylags.forward,noSim=100,seed=myseed,scale=rep(1,mylags.forward),logrecruitVec="historic",oldstyle=oldstyle,Myearrange=Myearrange)
    ftimes = attr(sim.tac,"times")
    qftimes=y2q(ftimes)[1:length(ftimes)]
    
    CWpred=CWpredict(aux)
    if(!is.null(ssbmin)){
        FtacScales = get.sim.scaled.to.blim(obj,jointrep,data,stepsback=0,lags.forward=mylags.forward,noSim=mynosim,seed=myseed,blimtarget=c(ssbmin),SW=SW[,qftimes],PM=PM[,qftimes],trace=TRUE,logrecruitVec="historic",prob=prob)
    } else {
        FtacScales = Fscaling
    }
    
    sim.tac = forecast.sesam(obj,jointrep,data,stepsback=0,lags.forward=mylags.forward,noSim=mynosim,seed=myseed,scale=rep(FtacScales,mylags.forward),logrecruitVec="historic",oldstyle=oldstyle,Myearrange=Myearrange)

    ssb.tac<-get.ssb.sim(sim.tac,SW[,qftimes],PM[,qftimes])
    fbar.tac = get.fbar.sim(sim.tac,fbarrange=2:3)
    median.fbar.tac = apply( do.call("rbind",fbar.tac),2,median)
    median.ssb.tac = apply( do.call("rbind",ssb.tac),2,median)
    ssb5 = apply( do.call("rbind",ssb.tac),2,quantile,probs=0.05)
    catch.tac = get.catch.sim(sim.tac,CWpred$CW[,qftimes],NM[,qftimes])
    median.catch.tac = apply( do.call("rbind",catch.tac),2,median)
    
    forecasttab <- data.frame(fbar=median.fbar.tac,ssb=median.ssb.tac,ssb5=ssb5,median.catch.tac)
    colnames(forecasttab)<-c(paste('F',data$fbarrange[1],data$fbarrange[2],sep=''),
                              "SSB","SSB 5th quantile","median catch")
    
    if(oldstyle) {
        forecasttab = rbind(forecasttab,c(NA,NA,NA,sum(forecasttab[,4])))
    }else{
        forecasttab[nrow(forecasttab),c(1,4)]=NA
        forecasttab = rbind(forecasttab,c(NA,NA,NA,sum(forecasttab[,4],na.rm=TRUE)))
    }                               
    rownames(forecasttab)<-c(ftimes,"Sum")
    list( sim.tac=sim.tac, ftimes=ftimes, forecasttab=forecasttab, CWpred=CWpred,SW=SW,NM=NM,PM=PM )

}

plot.forecast<-function(f,blim,fbarcols=2:3){
    par(mfrow=c(3,1))
    forecast.plotFbar(data,f$sim.tac,sdrep,fbarcols)
    forecast.plotSSB(data,f$sim.tac,sdrep,SW=f$SW,PM=f$PM,blim=blim)
    forecast.plotcatch(data,f$sim.tac,sdrep,aux,pl,CW=f$CWpred$CW,NM=f$NM)
}

aux$quarter = y2q(aux$auxt1)
PM = xtabs(auxPM ~ auxage + quarter,data=aux)/xtabs( ~ auxage + quarter,data=aux)
SW = xtabs(auxSW ~ auxage + quarter,data=aux)/xtabs( ~ auxage + quarter,data=aux)
NM = xtabs(auxM ~ auxage + quarter,data=aux)/xtabs( ~ auxage + quarter,data=aux)

## reference points
refyear = 2012 ## lowest biomass that gave above median recruitment 
blim = ssb[which( data$times==(refyear + 0.75))] ## Q4 biomass because forecast ends in Q4
bpa = blim*(exp(0.3*1.645))  
last.yearly.fbar = mean( tail(exp(logfbar[,1]),4))

nosims = 4000
fastforecast = FALSE
Myearrange=c(2011,2022.9) ## years to average M over (both inclusive)


source("src/forecast.R")

    fc.sq = do.forecast(4,123,nosims,PM,SW,NM,ssbmin=NULL,Fscaling=1,oldstyle=FALSE,Myearrange=Myearrange) ## F status quo forecast

if(!fastforecast){
    fc = do.forecast(4,123,nosims,PM,SW,NM,ssbmin=blim,oldstyle=FALSE,Myearrange=Myearrange) ## tac solving forecast

    fc.0 = do.forecast(4,123,nosims,PM,SW,NM,ssbmin=NULL,Fscaling=0,oldstyle=FALSE,Myearrange=Myearrange) ## F zero forecast
    fc.p50 = do.forecast(4,123,nosims,PM,SW,NM,ssbmin=blim,oldstyle=FALSE,prob=0.5,Myearrange=Myearrange) ## solve for P(SSB<blim)=0.5
    ##bpa=70000 ##Based on Bpa=Blim*(exp(0.3*1.645))
    
    fc.bpa = do.forecast(4,123,nosims,PM,SW,NM,bpa,oldstyle=FALSE,prob=0.5,Myearrange=Myearrange) ## solve for P(SSB<bpa)=0.5 (bpa=70000)
    
    
    ## F cap forecasts
    fcaps = c(0.2,0.3,0.4,0.5,0.6,0.7)
    ##fcaps = c(0.7)
    fcapfcs = list()
    for(fcr in 1:length(fcaps)){
        fcapfcs[[fcr]] = do.forecast(4,123,nosims,PM,SW,NM,ssbmin=NULL,Fscaling=fcaps[fcr]/last.yearly.fbar,oldstyle=FALSE,Myearrange=Myearrange)
    }	
}

stamp<-gsub('-[[:digit:]]{4}$','',gsub(':','.',gsub(' ','-',gsub('^[[:alpha:]]{3} ','',date()))))
setwd('res')
file.remove(dir(pattern='png$'))
file.remove(dir(pattern='html$'))

png(filename = paste(stamp,"_%03d.png", sep=''), width = 480, height = 480,
    units = "px", pointsize = 10, bg = "white")
plots()
par(mar=c(4,4,4,2))
if(file.exists("../run/mprof.RData")){
    load("../run/mprof.RData")
    plot(mults,lls,main="M likelihood profile",xlab="Multiplier",ylab="NLL")
}

if(!fastforecast){
    plot.forecast(fc,blim)
    plot.forecast(fc.0,blim)
  }

plot.forecast(fc.sq,blim)
    
if(!fastforecast){   
    plot.forecast(fc.p50,blim)
    plot.forecast(fc.bpa,bpa)
    for(fcr in 1:length(fcaps)){
	plot.forecast(fcapfcs[[fcr]],blim)
    }
}
dev.off()

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

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")
mrho = retroPlot()         
dev.off()    


##writeLines(unlist(tit.list),'titles.cfg') 

png(filename = paste("big_",stamp,"_%03d.png", sep=''), width = 1200, height = 1200, 
    units = "px", pointsize = 20, bg = "white")
plots()
par(mar=c(4,4,4,2))
if(exists("MPROF") && MPROF) plot(mults,lls,main="M likelihood profile",xlab="Multiplier",ylab="NLL")

if(!fastforecast){
    plot.forecast(fc,blim)
    plot.forecast(fc.0,blim)
    }

plot.forecast(fc.sq,blim)
    
if(!fastforecast){  
    plot.forecast(fc.p50,blim)
    plot.forecast(fc.bpa,bpa)
    for(fcr in 1:length(fcaps)){
	plot.forecast(fcapfcs[[fcr]],blim)
    }
}

dev.off()

##############
## Tables
##############

## Summary table
sel = seq(3,ncol(pl$logN),by=4)
logrecr = pl$logN[1,]
logrecr.sd = plsd$logN[1,]
logrecr.lo = logrecr - 2*logrecr.sd
logrecr.hi = logrecr + 2*logrecr.sd
logrecr[ -sel ] = NA
logrecr.lo[ -sel ] = NA
logrecr.hi[ -sel ] = NA

fbaryall<-fbaryall.hi<-fbaryall.lo<-rep(NA,length(FBAR))

sel = seq(1,ncol(pl$logN),by=4)
fbaryall[sel] = c(fbar.y,NA)
fbaryall.hi[sel] = c(fbar.y.hi,NA)
fbaryall.lo[sel] = c(fbar.y.lo,NA)


tab1<-cbind(exp(logrecr),exp(logrecr.lo),exp(logrecr.hi),SSB,SSB.lo,SSB.hi,TSB,TSB.lo,TSB.hi,fbaryall,fbaryall.lo,fbaryall.hi)

colnames(tab1)<-c('Recruits', 'Low', 'High', 'SSB', 'Low', 'High','TSB', 'Low', 'High', 
                  paste('F',data$fbarrange[1],data$fbarrange[2],sep=''), 'Low', 'High')
rownames(tab1)<-data$times
xtab(tab1, caption=paste('Table 1. Estimated recruitment, spawning stock biomass (SSB), 
                         and average fishing mortality for ages ',data$fbarrange[1], ' to ',data$fbarrange[2],
                         ' (F',data$fbarrange[1],data$fbarrange[1],')','.',sep=''), cornername='Time', 
     file=paste(stamp,'_tab1.html',sep=''), dec=c(0,0,0,0,0,0,0,0,0,3,3,3))


## Mean of summ. table
qq = factor(rep(1:4,length=length(SSB)))
tab1a = data.frame( value=c(mean( exp(logrecr),na.rm=TRUE),exp(mean(logrecr,na.rm=TRUE)),median(exp(logrecr),na.rm=TRUE),tapply(SSB,qq,mean),tapply(TSB,qq,mean),
                       mean(fbaryall,na.rm=TRUE)))
rownames(tab1a)<-c("Avg. recruitment","Geom. mean recruitment","Median recruitment",paste("Avg SSB Q",1:4),paste("Avg TSB Q",1:4),"Avg. FBAR");
xtab(tab1a, caption='Table 1a. Mean summary table.', cornername='', 
     file=paste(stamp,'_tab1a.html',sep=''), dec=rep(2,ncol(tab1a)))


## N table
tab2 = t(exp(pl$logN))
tab2[,1][ tab2[,1] < 1 ] = 0
colnames(tab2)<-data$ages
rownames(tab2)<-data$times
xtab(tab2, caption='Table 2. Estimated stock numbers.', cornername='Time\\Age', 
     file=paste(stamp,'_tab2.html',sep=''), dec=rep(0,ncol(tab2)))

## F table
tab3<-t(exp(pl$logF))
rownames(tab3)<-data$times
colnames(tab3)<-paste(1:ncol(tab3)+min(data$ages)-1,c(rep('',ncol(tab3)-1),'+'),sep='')
xtab(tab3, caption='Table 3. Estimated fishing mortalities.', cornername='Year\\Age', 
     file=paste(stamp,'_tab3.html',sep=''), dec=rep(3,ncol(tab3)))

## Q table
key=data$keyLogQ
o<-order(key[is.finite(key)])
fleet<-row(key)[is.finite(key)][o]
age<-data$ages[ col(key)[is.finite(key)][o] ]
logQ = sum[grep("logQ",rownames(sum)),1:2]
mylogQ = logQ[key[is.finite(key)][o]+1,]

tab4<-cbind(fleet,age,exp(mylogQ[,1]),exp(mylogQ[,1] -2*mylogQ[,2]), exp(mylogQ[,1] +2*mylogQ[,2])) 
colnames(tab4)<-c('Fleet number', 'Age', 'Catchability', 'Low', 'High')
rownames(tab4)<-1:nrow(tab4)
xtab(tab4, caption='Table 4. Estimated catchabilities.', 
     cornername='Index',file=paste(stamp,'_tab4.html',sep=''), dec=c(0,0,5,5,5))

tab5<-rbind(c(opt$objective, length(opt$par),AIC.sesam(opt)), 
            c(opt$objective, length(opt$par),AIC.sesam(opt)))
rownames(tab5)<-c('Base', 'Current')
colnames(tab5)<-c('Negative log likelihood','Number of parameters','AIC')
if(tab5[1,2]!=tab5[2,2]){
    if(tab5[1,2]<tab5[2,2]){tab5<-tab5[2:1,]}
    tab5<-cbind(tab5,'Degrees of freedom'=c(NA,tab5[1,2]-tab5[2,2]))
    tab5<-cbind(tab5,'P value'=c(NA,1-pchisq(2*(tab5[2,1]-tab5[1,1]),tab5[2,3])))
    xtab(tab5, caption='Table 6. Likelihood values and p-value for the likelihood 
                      ratio test for the model reduction. It is entirely the users 
                      responsibility to ensure that a likelihood ratio test is 
                      meaningful.', 
         cornername='Model',file=paste(stamp,'_tab5.html',sep=''), 
         dec=c(2,0,0,4))
}else{
    xtab(tab5, caption='Table 5. Likelihood values.', 
         cornername='Model',file=paste(stamp,'_tab5.html',sep=''), 
         dec=c(2,0,2))
}

## forecast tables

if(!fastforecast) xtab(fc$forecasttab,caption="Table 6. Forecast with F scaled such that the fifth percentile of the SSB distribution one year ahead equals Blim",file=paste(stamp,'_tab6.html',sep=''),dec=rep(4,ncol(x))) ## alternatively dec=c(4,rep(2,ncol(x)-1))

xtab(fc.0$forecasttab,caption="Table 7. Forecast with zero F",file=paste(stamp,'_tab7.html',sep=''),dec=rep(4,ncol(x)))

xtab(fc.sq$forecasttab,caption="Table 8. Forecast with status quo F",file=paste(stamp,'_tab8.html',sep=''),dec=rep(4,ncol(x)))

if(!fastforecast){
    xtab(fc.p50$forecasttab,caption="Table 9. Forecast with with F scaled such that the median of the SSB distribution one year ahead equals Blim",file=paste(stamp,'_tab9.html',sep=''),dec=rep(4,ncol(x)))
    xtab(fc.bpa$forecasttab,caption="Table 9a. Forecast with with F scaled such that the median of the SSB distribution one year ahead equals Bpa",file=paste(stamp,'_tab9a.html',sep=''),dec=rep(4,ncol(x)))
    for(fcr in 1:length(fcaps)){
        tabnam = paste0("9a",fcr)
        xtab(fcapfcs[[fcr]]$forecasttab,caption=paste0("Table ",tabnam,". Forecast with F=",fcaps[fcr],"."),file=paste(stamp,'_tab',tabnam,'.html',sep=''),dec=rep(4,ncol(x)))
    }
}


tc=data.frame(year=unique(data$t1),observed=totcatch,predicted=predtotcatch)
xtab(tc,caption="Table 9b. Observed and predicted total catches (mass) by quarter.",file=paste(stamp,'_tab9b.html',sep=''),dec=c(2,0,0))

## table with predicted catch weights
xtab(fc.sq$CWpred$CW,caption="Table 9c. Catch weights used in forecasts.",file=paste(stamp,'_tab9c.html',sep=''),cornername=c("Age/Quarter"),dec=rep(3,ncol(x)))

## Reference point table
xtab( x=t(data.frame(reference_year=refyear,reference_quarter=4,B_lim=blim,Bpa=bpa)),
    caption="Table 9d. Reference points",file=paste(stamp,'_tab9d.html',sep=''),dec=rep(2,ncol(x)))

avgM = aggregate(auxM ~ auxage+quarter,FUN=mean,data=subset(aux,auxt1 >= Myearrange[1] & auxt1 <= Myearrange[2]))
colnames(avgM) <- gsub("aux","",colnames(avgM))
## M used in forecast
xtab(avgM,caption="Table 9e. M used in forecasts.",file=paste(stamp,'_tab9e.html',sep=''),dec=c(0,0,2))

partab = data.frame(Estimate = sdrep$par.fixed, StdErr = sqrt(diag(sdrep$cov.fixed)))
rownames(partab) <- paste(1:nrow(partab),names(sdrep$par.fixed))
xtab(partab,caption="Table 9f. Fixed effects parameters.",file=paste(stamp,'_tab9f.html',sep=''),dec=rep(3,ncol(x)))




sel = seq(3,ncol(pl$logN),by=4)
logrecr = pl$logN[1,sel]
logrecr.sd = plsd$logN[1,sel]
logrecr.lo = logrecr - 2*logrecr.sd
logrecr.hi = logrecr + 2*logrecr.sd

sesamsum = list(data=data,param=param,pl=pl,plsd=plsd,rep=rep,opt=opt,logrecr=logrecr,logrecr.lo=logrecr.lo,logrecr.hi=logrecr.hi,SSB=SSB,SSB.lo=SSB.lo,SSB.hi=SSB.hi,FBAR=FBAR,FBAR.lo=FBAR.lo,FBAR.hi=FBAR.hi,mrho=mrho)
save(sesamsum,file="../run/sum.RData")