## Note, this plotscript assumes that sesam.R has been run already in the same R session.

####################
## OSA residuals
####################
epsfleet = data$eps[data$fleet]
if(doOSA && !(exists("is.RETRO") && is.RETRO)){
    ## condition on these first data points to make OSA work
    mycond=1:12
    
    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)
}

##save(data,param,pl,plsd,rep,jointrep,sdrep,opt,OSArep, file="run/samtmb.RData")

#####################
## 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=c(2,3))
    if(doOSA){
        for(f in unique(obs$fleet)){
            sel = which(data$fleet==f)
            sel=setdiff(sel,mycond)
            bp(data$t1[sel],data$ageFrom[sel],OSA$residual[sel]*(!is.na(data$obs[sel]) & (data$obs[sel]>epsfleet[sel])),scale=3,main=paste0("OSA fleet ",f),ylim=c(-1,4))
            v = OSA$residual[sel]*(!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=c(2,3))
    }
    
    ## full conditional residuals
    for(f in 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))
    }


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


    ## 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)
    
    ##plot(data$time[sel]-0.5,exp(logrecr),type="b",ylim=range(exp(c(logrecr.lo,logrecr.hi))),ylab="Recruitment",xlab="time")
    ##lines(data$time[sel]-0.5,exp(logrecr.lo),lty=2)
    ##lines(data$time[sel]-0.5,exp(logrecr.hi),lty=2)
    #lines(xsa$Year,xsa$recr,col="green",lwd=2)

    
    ## 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) )

    plot(ssb[selq1],exp(logrecr),ylim=c(0,max(exp(logrecr))),xlim=c(0,max(ssb[selq1])),xlab="SSB Q1",ylab="recruitment",pch=16,col=cols)
    disply = 0.02*max(exp(logrecr))
    displx = 0.02*max(ssb[selq1])
    text(ssb[selq1]+displx,exp(logrecr)+disply,labels=sapply(as.character(data$times[selq1]),substr,start=3,stop=4 ))

    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 ))

    if(doOSA){
        sel = which(obs$fleet==1)
        tab=xtabs(OSA$residual[sel] ~ data$ageFrom[sel]+ data$t1[sel])

        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")
            sel = which(obs$fleet==f)
            tab=xtabs(OSA$residual[sel] ~ data$ageFrom[sel]+ data$t1[sel])
            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))
            
        }

        ## residual corr

        par(mfrow=c(2,3))
        ccolors <- c("#A50F15","#DE2D26","#FB6A4A","#FCAE91","#FEE5D9","white",
                     "#EFF3FF","#BDD7E7","#6BAED6","#3182BD","#08519C")
        for(f in 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])
        }  

    }

    ## observed vs predicted
    par(mfrow=c(2,3))
    for(f in 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)
    }



    ## 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")

}## end plots    


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

do.forecast<-function(mylags.forward,myseed,mynosim,PM,SW,NM,ssbmin,Fscaling=NULL,oldstyle=TRUE,prob=0.05){
    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)
    ftimes = attr(sim.tac,"times")
    qftimes=y2q(ftimes)[1:length(ftimes)]

    ##browser()
    
    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)

    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
ssbminvec=numeric(4)
ssbminvectimes=numeric(4)
for(qq in 1:4){
    endq = qq
    refssbsel = which( y2q(data$times) == endq ) 
    ssbminvec[qq] = min(ssb[refssbsel])
    ssbminvectimes[qq] = data$times[ which(ssb == ssbminvec[qq]) ]
}
ssbmin = ssbminvec[2]
ssbmin = 42573 ## Q4 (2005) reference point from 2016 benchmark assessment run with revised IBTS data in 2020
last.yearly.fbar = mean( tail(exp(logfbar[,1]),4))

source("src/forecast.R")
fc = do.forecast(4,123,4000,PM,SW,NM,ssbmin,oldstyle=FALSE) ## tac solving forecast
fc.0 = do.forecast(4,123,4000,PM,SW,NM,ssbmin=NULL,Fscaling=0,oldstyle=FALSE) ## F zero forecast
fc.sq = do.forecast(4,123,4000,PM,SW,NM,ssbmin=NULL,Fscaling=1,oldstyle=FALSE) ## F status quo forecast
fc.p50 = do.forecast(4,123,4000,PM,SW,NM,ssbmin,oldstyle=FALSE,prob=0.5) ## solve for P(SSB<blim)=0.5
bpa=70000 ##Based on Bpa=Blim*(exp(0.3*1.645))
fc.bpa = do.forecast(4,123,4000,PM,SW,NM,bpa,oldstyle=FALSE,prob=0.5) ## 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)
fcapfcs = list()
for(fcr in 1:length(fcaps)){
  fcapfcs[[fcr]] = do.forecast(4,123,4000,PM,SW,NM,ssbmin=NULL,Fscaling=fcaps[fcr]/last.yearly.fbar,oldstyle=FALSE)
}	


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()    
plot.forecast(fc,ssbmin)
plot.forecast(fc.0,ssbmin)
plot.forecast(fc.sq,ssbmin)
plot.forecast(fc.p50,ssbmin)
plot.forecast(fc.bpa,bpa)
for(fcr in 1:length(fcaps)){
	plot.forecast(fcapfcs[[fcr]],ssbmin)
}
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()     
plot.forecast(fc,ssbmin)
plot.forecast(fc.0,ssbmin)
plot.forecast(fc.sq,ssbmin)
plot.forecast(fc.p50,ssbmin)
plot.forecast(fc.bpa,bpa)
for(fcr in 1:length(fcaps)){
	plot.forecast(fcapfcs[[fcr]],ssbmin)
}
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)), tapply(SSB,qq,mean),tapply(TSB,qq,mean),
                       mean(fbaryall,na.rm=TRUE)))
rownames(tab1a)<-c("Avg. recruitment","Geom. mean 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)), 
            c(opt$objective, length(opt$par)))
rownames(tab5)<-c('Base', 'Current')
colnames(tab5)<-c('Negative log likelihood','Number of parameters')
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))
}

## forecast tables
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)))
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)))
}
## Reference point table
##xtab( data.frame(SSB=ssbminvec,Time=ssbminvectimes),caption="Table 9. Minimum of estimated historic SSB by quarter",file=paste(stamp,'_tab9.html',sep=''),dec=rep(2,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$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)))