# Script to generate plots from the state-space assessment model
#  
# Anders Nielsen <anders@nielsensweb.org> Aug. 2012
oldwd<-getwd()
#system("rm ../*.zip") 
source('src/common.R')
source('src/datavalidator.R')

col.pref<-cbind(c(0,0,0),c(255,204,0),c(102,0,153),c(102,204,0),
                c(255,0,153),c(255,0,0),c(153,0,102),c(51,204,255),
                c(153,153,153))

palette(apply(rgb2hsv(col.pref),2,function(x)hsv(x[1],x[2],x[3])))

############################## read data ##############################


############################## read data ##############################

fit.current<-read.fit('run/sam')

#mybase<-"../../user3/nscod-2014-base/run/"
mybase<-"baserun/"

BASE.OK<-file.exists(paste(mybase,'sam.par',sep=""))
if(BASE.OK){
  fit.base<-read.fit(paste(mybase,'sam',sep=""))
}else{
  fit.base<-fit.current
  file.copy(paste('run',dir('run'), sep='/'),'baserun')
}



# Figure out a few basics 
minAge<-min(fit.current$res[,3])  
maxAge<-max(fit.current$res[,3])  
noN<-maxAge-minAge+1
noFleet<-max(fit.current$res[,2])  
conffile<-'run/model.cfg'
range.fbar<-fit.current$keys$fbarRange 
fbarlab=substitute(bar(F)[X-Y],list(X=range.fbar[1],Y=range.fbar[2]))
fleetnames<-c("Total catches",names(read.surveys('data/survey.dat')))


#######

# Read in data 
setwd('data')
catch.no<-read.ices('cn.dat')
land.no<-read.ices('lf.dat')
stock.mean.weight<-read.ices('sw.dat')
natural.mortality<-read.ices('nm.dat')
setwd('..')
# Modify to 7+ data
cutage<-maxAge
low<-minAge
GE<-which(as.numeric(colnames(catch.no))>=cutage)
E<-which(as.numeric(colnames(catch.no))==cutage)
w<-catch.no[,GE]/rowSums(catch.no[,GE])
wex<-rbind(w,w[nrow(w),])
catch.no[,E]<-rowSums(catch.no[,GE])
catch.no<-catch.no[,low:E]
land.no[,E]<-rowSums(land.no[,GE])
land.no<-land.no[,low:E]
lf<-ifelse(catch.no>0,land.no/catch.no,1)
stock.mean.weight[,E]<-rowSums(stock.mean.weight[,GE]*wex)
stock.mean.weight<-stock.mean.weight[,low:E]
natural.mortality[,E]<-rowSums(natural.mortality[,GE]*wex)
MM<-natural.mortality[,low:E]



#######


doextra<-FALSE
if(file.exists('conf/viewextra.cfg')){  
  extra<-read.table('conf/viewextra.cfg')  
  if(nrow(extra)>1){
    doextra<-TRUE
    extra<-extra[-1,]
    exNames<-paste("../../user3/",as.character(extra$V1),"/run/sam",sep="")
    exCol<-as.character(extra$V2)
    fit.extra<-lapply(exNames,function(x){if(file.exists(x))ret<-read.fit(x);ret}) 
  }
}  
  
  
LO.OK<-file.exists('run/LO.RData')
if(LO.OK){
  load('run/LO.RData')
}

RETRO.OK<-file.exists('run/RETRO.RData')
if(RETRO.OK){
  load('run/RETRO.RData')
}

SIM.OK<-file.exists('run/SIM.RData')
if(SIM.OK){
  load('run/SIM.RData')
}


# Figure out a few basics 
minAge<-min(fit.current$res[,3])  
maxAge<-max(fit.current$res[,3])  
noN<-maxAge-minAge+1
noFleet<-max(fit.current$res[,2])  
conffile<-'run/model.cfg'
range.fbar<-fit.current$keys$fbarRange 
fbarlab=substitute(bar(F)[X-Y],list(X=range.fbar[1],Y=range.fbar[2]))
fleetnames<-c("Total catches",names(read.surveys('data/survey.dat')))


plotcounter<-1 
tit.list<-list()

setcap<-function(title="", caption=""){   
 tit.list[length(tit.list)+1]<<-paste("# Title",plotcounter,"\n")
 tit.list[length(tit.list)+1]<<-paste(title,"\n")
 tit.list[length(tit.list)+1]<<-paste("# Caption",plotcounter,"\n")
 tit.list[length(tit.list)+1]<<-paste(caption,"\n")
 plotcounter<<-plotcounter+1 
}


############################## plots ##############################
plots<-function(){
  par(cex.lab=1.5, cex.axis=1.5, mar=c(5,5,1,1))


  # SSB plot
  baseplot(fit.current, 'ssb', trans=function(x)x/1000, ylab='SSB', drop=0)
  addlines(fit.base, 'ssb', col=gray(.5), trans=function(x)x/1000, drop=0, lwd=2)
  if(doextra){
    for(i in 1:length(fit.extra)){
      addlines(fit.extra[[i]], 'ssb', col=exCol[i], trans=function(x)x/1000, lwd=3, drop=0, ci=FALSE)
    }
  }
  NN<-exp(fit.current$stateEst[,1:noN])
  FF<-exp(fit.current$stateEst[,-c(1:noN)])[,fit.current$keys$keyLogFsta[1,]]
  EBM<-rowSums((FF/rowSums(FF))*NN*stock.mean.weight)
#  lines(fit.current$years, EBM/1000, col="grey", lwd=4)
#  lines(fit.current$years, EBM/1000, col="grey", lwd=4)
  
  
  stampit()
  setcap("Spawning stock biomass", "Spawning stock biomass. 
         Estimates from the current run and point wise 95% confidence 
         intervals are shown by black line and shaded area. 
         The perviously saved 'baserun' is the overlying gray lines.")

  # Fbar plot
  baseplot(fit.current, 'fbar', ylab=fbarlab)
  addlines(fit.base, 'fbar', col=gray(.5), lwd=2)
  if(doextra){
    for(i in 1:length(fit.extra)){
      addlines(fit.extra[[i]], 'fbar', col=exCol[i], lwd=3, ci=FALSE)
    }
  }
  stampit()
  setcap("Average F", "Average fishing mortalities. Estimates from the 
          current run and point wise 95% confidence intervals are shown 
          by black line and shaded area. The perviously saved 'baserun' 
          is the overlying gray lines.")


  # Recruit plot
  baseplot(fit.current, 'R', ylab='Recruits', trans=function(x)x/1000)
  addlines(fit.base, 'R', trans=function(x)x/1000, lwd=2, col=gray(.5))
  if(doextra){
    for(i in 1:length(fit.extra)){
      addlines(fit.extra[[i]], 'R', col=exCol[i], trans=function(x)x/1000, ci=FALSE, lwd=3)
    }
  }
  stampit()
  setcap("Recruits", "Number of recruits entering the population. Estimates 
          from the current run and point wise 95% confidence intervals are shown 
          by black line and shaded area. The perviously saved 'baserun' is the 
          overlying gray lines.")
  # scaleplot
  
  
  baseplot(fit.current, "logscale", xx=fit.current$keys$keyScaledYears, trans=exp, drop=0)
  
  #yy<-fit.current$keys$keyScaledYears
  #addlines(fit.base, "logscale", xx=yy[-length(yy)], trans=exp, lwd=2, col=gray(.5), drop=0)
  #lines(fit.current$keys$keyScaledYears,exp(fit.current$logscale[,3]), lty="dashed", lwd=2)
  #lines(fit.current$keys$keyScaledYears,exp(fit.current$logscale[,4]), lty="dashed", lwd=2)
  stampit()
  setcap("Catch scaling", "")
  
  
  # Residual plot 
  div<-list(c(1,1),c(2,1),c(3,1),c(2,2),c(3,2),c(3,2),c(3,3),c(3,3),c(3,3),c(4,3),c(4,3),c(4,3))
  op<-par(mfrow=div[[noFleet]], mar=c(4,4,1,2), mgp=c(2,1,0))
  #fleet.names<-paste(c('Total catches', paste('S',1:(noFleet-1),sep='')))
  for(f in 1:noFleet){
    cur<-fit.current$res[fit.current$res[,2]==f,]
    scale<-3
    bp(cur[,1],cur[,3],cur[,6], ylim=c(minAge-.5,maxAge+.5), xlim=range(fit.current$years), 
       xlab='year', ylab='Age', main=fleetnames[f], scale=scale)
  }
  par(op)
  setcap("Normalized residuals", "Normalized residuals for the current run. Blue circles 
          indicate positive residuals (obs larger than predicted) and filled green circles 
          indicate negative residuals.")
          
          
  
  for(f in 1:noFleet){
    curf<-fit.current$res[fit.current$res[,2]==f,]
    loca<-sort(unique(curf[,3]))
    op<-par(mfrow=div[[length(loca)]], mar=c(4,4,1,2), mgp=c(2,1,0))
    for(a in loca){
      sub<-curf[curf[,3]==a,]
      plot(sub[,1],sub[,4], xlab="year", ylab="logObs")
      lines(sub[,1],sub[,5])
      legend("topright", legend=paste("Fleet no ",f," , Age ",a,sep=""),bty="n")
    }
    par(op)
    setcap("Observed vs. fitted", "")
  }


  F<-exp(fit.current$stateEst[-nrow(fit.current$stateEst),-c(1:noN)])
  colnames(F)<-paste(1:ncol(F)+minAge-1,c(rep('',ncol(F)-1),'+'),sep='')
  rownames(F)<-fit.current$years[-length(fit.current$years)]
  matplot(rownames(F),F,lty=1:ncol(F),col=1:ncol(F), type='l', xlab='Year', lwd=5)
  legend('topleft', col=1:ncol(F), lty=1:ncol(F), legend=colnames(F), bty='n', lwd=5)
  stampit()
  setcap()
  
  ###
  F<-exp(fit.current$stateEst[,-c(1:noN)])[,fit.current$keys$keyLogFsta[1,]]
  F<-F[-nrow(F),]
  years<-fit.current$years
  fbar.land<-rowMeans((F*lf)[,2:4])
  fbar.dis<-rowMeans((F*(1-lf))[,2:4])
  plot(years[-length(years)],fbar.land, lwd=5, col='blue', type='l', ylim=c(0,max(fbar.land)), 
      las=1, xlab='Year', ylab=expression(bar(F)[2-4]))
  lines(years[-length(years)],fbar.dis, lwd=5, col='red')
  legend('topleft', lty='solid', col=c('blue','red'), legend=c('Landed', 'Discarded'), bty='n', lwd=5)
  stampit()
  setcap()
  
  
   ###
   
  NN<-exp(fit.current$stateEst[,1:noN])
  FF<-exp(fit.current$stateEst[,-c(1:noN)])[,fit.current$keys$keyLogFsta[1,]]
   
  predNN<-matrix(NA,nrow=nrow(FF), ncol=ncol(FF))
  sdNmat<-matrix(exp(fit.current$est[fit.current$names=="logSdLogN"][fit.current$keys$keyVarLogN]),nrow=nrow(FF), ncol=ncol(FF), byrow=T)
  
  NNsd<-sqrt((exp(sdNmat^2)-1)*exp(2*fit.current$stateEst[,1:noN]+sdNmat^2))
  for(y in 2:nrow(FF)){
    predNN[y,1]<-NN[y-1,1]   
    for(a in 2:ncol(FF)){
      predNN[y,a]<-NN[y-1,a-1]*exp(-FF[y-1,a-1]-MM[y-1,a-1])
    }
    predNN[y,a]<-predNN[y,a]+NN[y-1,ncol(FF)]*exp(-FF[y-1,ncol(FF)]-MM[y-1,ncol(FF)])
  }
  res<-(NN-predNN)/NNsd
  bp(as.vector(row(res)), as.vector(col(res)),(as.vector(res)), scale=2, main="(o-p)/sd")
  stampit()
  setcap()
  
  
  if(LO.OK){
    labels<-c('All data', paste('w.o. ', fleetnames[-1], sep=''))
    cols<-c('red',2+1:length(LO))

    # SSB plot 
    baseplot(fit.current, 'ssb', trans=function(x)x/1000, ylab='SSB')
    for(i in 1:length(LO)){
      addlines(LO[[i]], 'ssb', col=i+2, trans=function(x)x/1000, ci=FALSE)
    }
    legend('bottomleft', legend=labels, lty='solid', col=cols, bty='n')
    stampit()
    setcap()

    # Fbar plot
    baseplot(fit.current, 'fbar', ylab=fbarlab)
    for(i in 1:length(LO)){
      addlines(LO[[i]], 'fbar', col=i+2, ci=FALSE)
    }
    legend('bottomleft', legend=labels, lty='solid', col=cols, bty='n')
    stampit()
    setcap()

    # Recruit plot
    baseplot(fit.current, 'R', ylab='Recruits', trans=function(x)x/1000)
    for(i in 1:length(LO)){
      addlines(LO[[i]], 'R', col=i+2,trans=function(x)x/1000, ci=FALSE)
    }
    legend('bottomleft', legend=labels, lty='solid', col=cols, bty='n')
    stampit()
    setcap()
    
    #Ellipse plot
    #civec<-civec(fit.current)
    #plot(civec[5:6],xlim=c(civec[5]+c(-10,10)*sqrt(civec[1])),
    #                ylim=c(civec[6]+c(-10,10)*sqrt(civec[4])), 
    #     type='n', xlab=expression(bar(F)[2-4]), ylab='SSB')
    #.CI.reg(civec, col='red')
    #for(i in 1:length(LO)){
    #  .CI.reg(LO[[i]]$civec, col=i+2)
    #}
    #legend('bottomleft', legend=labels, lty='solid', col=cols, bty='n')

  }

  if(RETRO.OK){
    # SSB plot 
    baseplot(fit.current, 'ssb', trans=function(x)x/1000, ylab='SSB', drop=0)
    for(i in 1:length(RETRO)){
      addlines(RETRO[[i]], 'ssb', col=i+2, trans=function(x)x/1000, with.dot=TRUE, ci=FALSE, drop=0)
    }
    ssbrho<-mean(unlist(lapply(RETRO, function(ret){i<-length(ret$ssb[,1]);(1-ret$ssb[i,1]/fit.current$ssb[i,1])})))
    ssbrholab=substitute(rho == X,list(X=round(ssbrho,3)))
    legend("topright", legend=ssbrholab, cex=1.5, bty="n")

    stampit()
    setcap()

    # Fbar plot
    baseplot(fit.current, 'fbar', ylab=fbarlab)
    for(i in 1:length(RETRO)){
      addlines(RETRO[[i]], 'fbar', col=i+2, with.dot=TRUE, ci=FALSE)
    }
    Frho<-mean(unlist(lapply(RETRO, function(ret){i<-length(ret$fbar[,1])-1;(1-ret$fbar[i,1]/fit.current$fbar[i,1])})))
    Frholab=substitute(rho == X,list(X=round(Frho,3)))
    legend("topright", legend=Frholab, cex=1.5, bty="n")
    stampit()
    setcap()

    # Recruit plot
    baseplot(fit.current, 'R', ylab='Recruits', trans=function(x)x/1000)
    for(i in 1:length(RETRO)){
      addlines(RETRO[[i]], 'R', col=i+2, with.dot=TRUE,trans=function(x)x/1000, ci=FALSE)
    }
    Rrho<-mean(unlist(lapply(RETRO, function(ret){i<-length(ret$R[,1])-1;(1-ret$R[i,1]/fit.current$R[i,1])})))
    Rrholab=substitute(rho == X,list(X=round(Rrho,3)))
    legend("topright", legend=Rrholab, cex=1.5, bty="n")

    stampit()
    setcap()
  }

  if(SIM.OK){
    # SSB plot 
    baseplot(fit.current, 'ssb', trans=function(x)x/1000, ylab='SSB')
    for(i in 1:length(SIM)){
      addlines(SIM[[i]], 'ssb', col=i+2, trans=function(x)x/1000, with.dot=TRUE, ci=FALSE)
    }
    stampit()
    setcap()

    # Fbar plot
    baseplot(fit.current, 'fbar', ylab=fbarlab)
    for(i in 1:length(SIM)){
      addlines(SIM[[i]], 'fbar', col=i+2, with.dot=TRUE, ci=FALSE)
    }
    stampit()
    setcap()

    # Recruit plot
    baseplot(fit.current, 'R', ylab='Recruits', trans=function(x)x/1000)
    for(i in 1:length(SIM)){
      addlines(SIM[[i]], 'R', col=i+2, with.dot=TRUE,trans=function(x)x/1000, ci=FALSE)
    }
    stampit()
    setcap()
  }

}


setwd('res')
file.remove(dir(pattern='png$'))
stamp<-gsub('-[[:digit:]]{4}$','',gsub(':','.',gsub(' ','-',gsub('^[[:alpha:]]{3} ','',date()))))
png(filename = paste(stamp,"_%03d.png", sep=''), width = 480, height = 480,
    units = "px", pointsize = 10, bg = "white")
  plots()    
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()    
dev.off()

#pdf(onefile=FALSE, width = 8, height = 8)
#  plots()    
#dev.off()


setwd(oldwd)

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

# SSB and Fbar table 

setwd('res')
file.remove(dir(pattern='html$'))

R<-fit.current$R[,c(1,3,4)]
tsb<-exp(fit.current$logtsb[,c(1,3,4)])
ssb<-exp(fit.current$logssb[,c(1,3,4)])
fbar<-exp(fit.current$logfbar[,c(1,3,4)])

# Summary table 
tab1<-cbind(R,tsb,ssb,fbar)
tab1[nrow(tab1),-c(7:9)]<-NA
colnames(tab1)<-c('Recruits', 'Low', 'High', 'TSB', 'Low', 'High', 'SSB', 'Low', 'High', 
                  paste('F',range.fbar[1],range.fbar[2],sep=''), 'Low', 'High')
rownames(tab1)<-fit.current$years
xtab(tab1, caption=paste('Table 1. Estimated recruitment, total stock biomass (TBS), spawning stock biomass (SSB), 
                         and average fishing mortality for ages ',range.fbar[1], ' to ',range.fbar[2],
                         ' (F',range.fbar[1],range.fbar[2],')','.',sep=''), cornername='Year', 
     file=paste(stamp,'_tab1.html',sep=''), dec=c(0,0,0,0,0,0,0,0,0,3,3,3))



# N table
tab2<-exp(fit.current$stateEst[,1:noN])
rownames(tab2)<-fit.current$years
colnames(tab2)<-paste(1:noN+minAge-1,c(rep('',noN-1),'+'),sep='')
xtab(tab2, caption='Table 2. Estimated stock numbers.', cornername='Year\\Age', 
     file=paste(stamp,'_tab2.html',sep=''), dec=rep(0,ncol(tab2)))


# F table
tab3<-exp(fit.current$stateEst[,-c(1:noN)])
tab3<-tab3[-nrow(tab3),]
rownames(tab3)<-fit.current$years[1:nrow(tab3)]
colnames(tab3)<-paste(1:ncol(tab3)+minAge-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)))

# Scale table 
tab4<-exp(fit.current$logscale[,c(1,3,4)])
if(nrow(tab4)>0){
  localconffile<-paste('..',conffile,sep='/')
  idx1<-grep('Years in which catch data are to be scaled',readLines(localconffile))
  idx2<-grep('Define Fbar range',readLines(localconffile))
  num<-scan(localconffile, skip=idx1, comment.char='#', quiet=TRUE, nlines=idx2-idx1)
  n<-num[1]
  y<-num[2:(n+1)]
  key<-matrix(num[-c(1:(n+1))], nrow=length(y), byrow=TRUE)

  if(!all(apply(key,1,function(x)length(unique(x))==1))){
    tab4<-matrix(tab4[,1][key], nrow=length(y))
    colnames(tab4)<-paste(1:ncol(tab3)+minAge-1,c(rep('',ncol(tab3)-1),'+'),sep='')
    rownames(tab4)<-y
    xtab(tab4, caption='Table 4. Estimated catch scaling factors', cornername='Year\\Age', 
         file=paste(stamp,'_tab4.html',sep=''), dec=rep(2,ncol(tab4)))
  }else{ # importent special case
    colnames(tab4)<-c('Catch multiplier', 'Low', 'High') 
    rownames(tab4)<-y[key[,1]]
    xtab(tab4, caption='Table 4. Estimated catch scaling factors (same for all ages within each year)', 
         cornername='Year', file=paste(stamp,'_tab4.html',sep=''), dec=rep(2,ncol(tab4)))
  }
}

conflines<-readLines(paste('../',conffile,sep=''))
idx1<-grep('Coupling of catchability PARAMETERS',conflines)
idx2<-grep('Coupling of power law model EXPONENTS',conflines)
num<-scan(paste('../',conffile,sep=''), skip=idx1, comment.char='#', quiet=TRUE, nlines=idx2-idx1)
a<-scan(paste('../',conffile,sep=''), comment.char='#', quiet=TRUE, n=2)
key<-matrix(num, ncol=diff(a)+1, byrow=TRUE)
o<-order(key[key>0])
fleet<-row(key)[key>0][o]
age<-col(key)[key>0][o]
tab5<-cbind(fleet,age,exp(fit.current$logFpar[key[key>0][o],c(1,3,4)])*1000)
colnames(tab5)<-c('Fleet number', 'Age', 'Catchability', 'Low', 'High')
rownames(tab5)<-1:nrow(tab5)
xtab(tab5, caption='Table 5. Estimated catchabilities multiplied by 1000.', 
     cornername='Index',file=paste(stamp,'_tab5.html',sep=''), dec=c(0,0,5,5,5))
##  
##  
tab6<-rbind(c(fit.base$nlogl, fit.base$nopar), 
            c(fit.current$nlogl, fit.current$nopar))
rownames(tab6)<-c('Base', 'Current')
colnames(tab6)<-c('Negative log likelihood','Number of parameters')
if(tab6[1,2]!=tab6[2,2]){
  if(tab6[1,2]<tab6[2,2]){tab6<-tab6[2:1,]}
  tab6<-cbind(tab6,'Degrees of freedom'=c(NA,tab6[1,2]-tab6[2,2]))
  tab6<-cbind(tab6,'P value'=c(NA,1-pchisq(2*(tab6[2,1]-tab6[1,1]),tab6[2,3])))
  xtab(tab6, 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,'_tab6.html',sep=''), 
       dec=c(2,0,0,4))
}else{
  xtab(tab6, caption='Table 6. Likelihood values.', 
       cornername='Model',file=paste(stamp,'_tab6.html',sep=''), 
       dec=c(2,0))
}

yy<-fit.current$years[-length(fit.current$years)]
scale<-rep(1,length(yy))
# Get catch scaling
phi<-1/rowMeans(exp(matrix(fit.current$logscale[fit.current$keys$keyParScaledYA],nrow=fit.current$keys$noScaledYears)))
scaleYears<-fit.current$keys$keyScaledYears
scale[yy%in%scaleYears]<-phi

#scale<-c(rep(1,length(fit.current$logLand[,1])-nrow(tab4)-1),1/tab4[,1],1)
tab7<-cbind(exp(fit.current$logLand[,1])*scale, 
            exp(fit.current$logDis[,1])*scale, 
            exp(fit.current$logCatch[,1])*scale,exp(fit.current$logCatch[,1]))


rownames(tab7)<-yy
colnames(tab7)<-c('Landings', 'Discards', 'Catch', 'Total Removal')

xtab(tab7, caption='Table 7. Catch table.', cornername='Year', 
     file=paste(stamp,'_tab7.html',sep=''), dec=rep(0,ncol(tab7)))


sub<-fit.current$res[fit.current$res[,2]==1,c(1,3,5)]
names(sub)<-c('year','age','logPred')
mat<-reshape(sub, idvar='year', timevar="age", direction='wide')
rownames(mat)<-mat[,1]
tab8<-exp(mat[,-1])
colnames(tab8)<-c(1:(ncol(tab8)-1),'+')
xtab(tab8, caption='Table 8. Predicted catches including catch multiplier', cornername='Year\\Age', 
     file=paste(stamp,'_tab8.html',sep=''), dec=rep(0,ncol(tab2)))


# Summary table 
sel=1:fit.current$nopar
tab9<-cbind(fit.current$est[sel], fit.current$std[sel], fit.current$est[sel]-2*fit.current$std[sel], fit.current$est[sel]+2*fit.current$std[sel])

colnames(tab9)<-c('Estimate', 'Standard deviation', '2.5% quantile', '97.5% quantile')
rownames(tab9)<-fit.current$names[sel]
xtab(tab9, caption=paste('Table 9. Estimates and standard deviations. The negative log likelihood is',fit.current$nlogl, 
                         'and the number of parameters is', fit.current$nopar), cornername='Name', 
     file=paste(stamp,'_tab9.html',sep=''), dec=c(5,3,3,3))


## variance parameters
conflines<-readLines(paste('../',conffile,sep=''))
idx1<-grep('Coupling of OBSERVATION VARIANCES',conflines)
idx2<-grep('Stock recruitment model code',conflines)
num<-scan(paste('../',conffile,sep=''), skip=idx1, comment.char='#', quiet=TRUE, nlines=idx2-idx1)
a<-scan(paste('../',conffile,sep=''), comment.char='#', quiet=TRUE, n=2)
key<-matrix(num, ncol=diff(a)+1, byrow=TRUE)
o<-order(key[key>0])
fleet<-row(key)[key>0][o]
age<-col(key)[key>0][o]
tab9a<-cbind(fleet,age,exp(fit.current$est[fit.current$names=='logSdLogObs'][key[key>0][o]]))


colnames(tab9a)<-c('Fleet number', 'Age', 'Estimate')
rownames(tab9a)<-1:nrow(tab9a)
xtab(tab9a, caption='Table 9a. Estimated observation SDs.', 
     cornername='Index',file=paste(stamp,'_tab9a.html',sep=''), dec=c(0,0,5))
     

source('../src/sisam.R')

cat('',file='footprint')

setwd(oldwd)


############################## FORECAST#########################################

FORECAST.OK<-file.exists('run/FORECAST.RData')
if(FORECAST.OK){
  load('run/FORECAST.RData')
 
  plotfore<-function(fit.current, forecast){
    par(mfrow=c(2,2), cex.lab=1.5, cex.axis=1.5, mar=c(5,5,5,1))
    py<-forecast$last.year.used+0:4
    summ<-function(x){
      x<-as.vector(x)
      c(mean(x),sd(x),quantile(x,c(2.5,5,25,50,75,95,97.5)/100))
    }
  
    baseplot(fit.current, 'logfbar', trans=exp, ylab=fbarlab, 
             ci=FALSE, line=FALSE, xlim=c(min(fit.current$years), forecast$last.year.used+4))
    abline(v=py[1], col=grey(.8), lwd=3)
    addlines(fit.current, 'logfbar', trans=exp, col='red', lwd=2)
    fbar<-rbind(summ(forecast$last.fbar.sim), summ(forecast$s1.fbar.sim), summ(forecast$s2.fbar.sim), summ(forecast$s3.fbar.sim), summ(forecast$s4.fbar.sim))
    lines(py,fbar[,1], col='blue', lwd=2)
    lines(py,fbar[,3], col='blue', lty='dashed', lwd=2)
    lines(py,fbar[,9], col='blue', lty='dashed', lwd=2)
  
    baseplot(fit.current, 'ssb', trans=function(x)x/1000, ylab='SSB', 
             ci=FALSE, line=FALSE, xlim=c(min(fit.current$years), forecast$last.year.used+5))
    abline(v=py[1], col=grey(.8), lwd=3)
    addlines(fit.current, 'ssb', col='red', trans=function(x)x/1000, lwd=2)
    ssb<-rbind(summ(forecast$last.ssb.sim), summ(forecast$s1.ssb.sim), summ(forecast$s2.ssb.sim), summ(forecast$s3.ssb.sim), summ(forecast$s4.ssb.sim), summ

(forecast$s5.ssb.sim))
    pyl<-forecast$last.year.used+0:5
    lines(pyl,ssb[,1]/1000, col='blue', lwd=2)
    lines(pyl,ssb[,3]/1000, col='blue', lty='dashed', lwd=2)
    lines(pyl,ssb[,9]/1000, col='blue', lty='dashed', lwd=2)
  
    baseplot(fit.current, 'R', trans=function(x)x/1000, ylab='Recruitment', 
             ci=FALSE, line=FALSE, xlim=c(min(fit.current$years), forecast$last.year.used+4))
    abline(v=py[1], col=grey(.8), lwd=3)
    addlines(fit.current, 'R', col='red', trans=function(x)x/1000, lwd=2)
    R<-rbind(summ(forecast$last.state.sim[,1]), summ(forecast$s1.state.sim[,1]), summ(forecast$s2.state.sim[,1]), summ(forecast$s3.state.sim[,1]), summ(forecast

$s4.state.sim[,1]))
    lines(py,exp(R[,1])/1000, col='blue', lwd=2)
    lines(py,exp(R[,3])/1000, col='blue', lty='dashed', lwd=2)
    lines(py,exp(R[,9])/1000, col='blue', lty='dashed', lwd=2)
    
    baseplot(fit.current, 'logCatch', trans=function(x)exp(x)/1000, ylab='Catch', 
             ci=FALSE, line=FALSE, xlim=c(min(fit.current$years), forecast$last.year.used+4), drop=0)
    abline(v=py[1], col=grey(.8), lwd=3)
    addlines(fit.current, 'logCatch', col='red', trans=function(x)exp(x)/1000, lwd=2, drop=0)
    catch<-rbind(summ(forecast$last.catch.sim), summ(forecast$s1.catch.sim), summ(forecast$s2.catch.sim), summ(forecast$s3.catch.sim), summ(forecast

$s3.catch.sim))
    lines(py,catch[,1]/1000, col='blue', lwd=2)
    lines(py,catch[,3]/1000, col='blue', lty='dashed', lwd=2)
    lines(py,catch[,9]/1000, col='blue', lty='dashed', lwd=2)
    
    tac<-round(summ(forecast$s2.land.sim))
    title(paste('TAC = ',tac[6],' (',tac[3],';',tac[9],')', sep=''), outer=TRUE, line=-3, cex.main=.9)
  }
   
  
  setwd('res')
  
  unlink(dir(pattern='XXX'))
  
  png(filename = paste(stamp,"_XXX%03d.png", sep=''), width = 480, height = 480, 
      units = "px", pointsize = 10, bg = "white")
    lapply(1:length(forecast),function(i){plotfore(fit.current, forecast[[i]]); title(foretab[i,ncol(foretab)], outer=TRUE, line=-2)})
  dev.off()
  
  png(filename = paste("big_",stamp,"_XXX%03d.png", sep=''), width = 1200, height = 1200, 
      units = "px", pointsize = 20, bg = "white")
    lapply(1:length(forecast),function(i){plotfore(fit.current, forecast[[i]]); title(foretab[i,ncol(foretab)], outer=TRUE, line=-2)})
  
  dev.off()
  
  
  foreTab<-function(fore){
    Fm<-c(1,fore$fadj)
    summ<-function(x)c(quantile(x,c(5,25,50,75,95)/100), mean(x))
    # more robust version
    summr<-function(x){xx<-c(quantile(x,c(5,25,50+c(-2,2),75,95)/100), mean(x)); c(xx[1:2],mean(xx[3:4]),xx[5:7])}
    mean(exp(colMeans(forecast[[1]]$last.state.sim)[c(11,12,13,13)]))
    #fbar<-fore$fbar.approx[-1]
    fbar<-cbind(summ(fore$last.fbar.sim),summ(fore$s1.fbar.sim),summ(fore$s2.fbar.sim),summ(fore$s3.fbar.sim),summ(fore$s4.fbar.sim))[3,-1]
    ssb<-cbind(summ(fore$last.ssb.sim),summ(fore$s1.ssb.sim),summ(fore$s2.ssb.sim),summ(fore$s3.ssb.sim),summ(fore$s4.ssb.sim))[3,-1]
    catch<-cbind(summ(fore$last.catch.sim),summ(fore$s1.catch.sim),summ(fore$s2.catch.sim),summ(fore$s3.catch.sim),summ(fore$s4.catch.sim))[3,-1]
    land<-cbind(summ(fore$last.land.sim),summ(fore$s1.land.sim),summ(fore$s2.land.sim),summ(fore$s3.land.sim),summ(fore$s4.land.sim))[3,-1]
    fbarl<-cbind(summ(fore$last.fbarl.sim),summ(fore$s1.fbarl.sim),summ(fore$s2.fbarl.sim),summ(fore$s3.fbarl.sim),summ(fore$s4.fbarl.sim))[3,-1]
    R<-cbind(summr(exp(fore$last.state.sim[,1])),
                 summr(exp(fore$s1.state.sim[,1])),
                 summr(exp(fore$s2.state.sim[,1])),
                 summr(exp(fore$s3.state.sim[,1])),
                 summr(exp(fore$s4.state.sim[,1])))[3,-1]
    X<-c(fbar,ssb,catch,land,fbarl,R,summ(fore$s5.ssb.sim)[3])
    return(X)
  }
  
  tabX<-do.call(rbind,lapply(forecast, foreTab))
  
  colnames(tabX)<-c('Fbar15','Fbar16','Fbar17','Fbar18','SSB15','SSB16','SSB17','SSB18','C15','C16','C17','C18',  
                    'L15','L16','L17','L18','Fbar15(L)','Fbar16(L)','Fbar17(L)','Fbar18(L)','R15','R16','R17','R18','SSB19')
  rownames(tabX)<-foretab[,ncol(foretab)]
  
  tabXb<-cbind(tabX[,13:16],tabX[,9:12]-tabX[,13:16],tabX[,17:20], (tabX[,14]-oldtac)/oldtac*100)
  colnames(tabXb)<-c('L15','L16','L17','L18','D15','D16','D17','D18','Fbar15(L)','Fbar16(L)','Fbar17(L)','Fbar18(L)','TAC % change')
  rownames(tabXb)<-foretab[,ncol(foretab)]
  
  tabXa<-tabX[,c(1:8,25,9:12)]
  tabXa<-cbind(tabXa,(tabXa[,7]-tabXa[,6])/tabXa[,6]*100)
  colnames(tabXa)<-c('Fbar15','Fbar16','Fbar17','Fbar18','SSB15','SSB16','SSB17','SSB18','SSB19','C15','C16','C17','C18','SSB % change')
  
  xtab(tabXa[,c(1:9,14,10:13)], caption='Table X. Forecasts (Fs and Catches are including estimated unallocated removals)', cornername='Description\\Variable', 
       file=paste(stamp,'_tabX.html',sep=''), dec=c(2,2,2,2,0,0,0,0,0,0,0,0,0,0))
  
  tabXb<-tabXb[,c(1:12)]
  xtab(tabXb, caption='Table Xb. Forecasts (split in landings and discards)', cornername='Description\\Variable', 
       file=paste(stamp,'_tabXb.html',sep=''), dec=c(0,0,0,0,0,0,0,0,2,2,2,2))
  
  mult<-1.0#mean(exp(fit.current$logscale[nrow(fit.current$logscale)-0:2,1]))
  tabXc<-tabXb[,c(1:8)]/mult
  tabXc<-cbind(tabXc,tabX[,21:24], (tabXc[,2]-oldtac)/oldtac*100)
  colnames(tabXc)<-c('L15','L16','L17','L18','D15','D16','D17','D18','R15','R16','R17','R18','TAC % change')
  xtab(tabXc[,c(1:4,13,5:12)], caption='Table Xc. Forecasts (L and D without catch multiplier, and R)', cornername='Description\\Variable', 
       file=paste(stamp,'_tabXc.html',sep=''), dec=c(0,0,0,0,0,0,0,0,0,0,0,0,0))

  




##
foreTab<-function(fore){
    Fm<-c(1,fore$fadj)
    summ<-function(x)c(quantile(x,c(5,25,50,75,95)/100), mean(x))
    # more robust version
    summr<-function(x){xx<-c(quantile(x,c(5,25,50+c(-2,2),75,95)/100), mean(x)); c(xx[1:2],mean(xx[3:4]),xx[5:7])}
    mean(exp(colMeans(forecast[[1]]$last.state.sim)[c(11,12,13,13)]))
    #fbar<-fore$fbar.approx[-1]
    fbar<-cbind(summ(fore$last.fbar.sim),summ(fore$s1.fbar.sim),summ(fore$s2.fbar.sim),summ(fore$s3.fbar.sim),summ(fore$s4.fbar.sim))[1,-1]
    ssb<-cbind(summ(fore$last.ssb.sim),summ(fore$s1.ssb.sim),summ(fore$s2.ssb.sim),summ(fore$s3.ssb.sim),summ(fore$s4.ssb.sim))[1,-1]
    catch<-cbind(summ(fore$last.catch.sim),summ(fore$s1.catch.sim),summ(fore$s2.catch.sim),summ(fore$s3.catch.sim),summ(fore$s4.catch.sim))[1,-1]
    land<-cbind(summ(fore$last.land.sim),summ(fore$s1.land.sim),summ(fore$s2.land.sim),summ(fore$s3.land.sim),summ(fore$s4.land.sim))[3,-1]
    fbarl<-cbind(summ(fore$last.fbarl.sim),summ(fore$s1.fbarl.sim),summ(fore$s2.fbarl.sim),summ(fore$s3.fbarl.sim),summ(fore$s4.fbarl.sim))[1,-1]
    R<-cbind(summr(exp(fore$last.state.sim[,1])),
                 summr(exp(fore$s1.state.sim[,1])),
                 summr(exp(fore$s2.state.sim[,1])),
                 summr(exp(fore$s3.state.sim[,1])),
                 summr(exp(fore$s4.state.sim[,1])))[1,-1]
    X<-c(fbar,ssb,catch,land,fbarl,R,summ(fore$s5.ssb.sim)[1])
    return(X)
  }
  
  tabX<-do.call(rbind,lapply(forecast, foreTab))
  
  colnames(tabX)<-c('Fbar15','Fbar16','Fbar17','Fbar18','SSB15','SSB16','SSB17','SSB18','C15','C16','C17','C18',  
                    'L15','L16','L17','L18','Fbar15(L)','Fbar16(L)','Fbar17(L)','Fbar18(L)','R15','R16','R17','R18','SSB19')
  rownames(tabX)<-foretab[,ncol(foretab)]
  
  tabXb<-cbind(tabX[,13:16],tabX[,9:12]-tabX[,13:16],tabX[,17:20], (tabX[,14]-oldtac)/oldtac*100)
  colnames(tabXb)<-c('L15','L16','L17','L18','D15','D16','D17','D18','Fbar15(L)','Fbar16(L)','Fbar17(L)','Fbar18(L)','TAC % change')
  rownames(tabXb)<-foretab[,ncol(foretab)]
  
  tabXa<-tabX[,c(1:8,25,9:12)]
  tabXXa<-cbind(tabXa,(tabXa[,7]-tabXa[,6])/tabXa[,6]*100)
  colnames(tabXXa)<-c('Fbar15','Fbar16','Fbar17','Fbar18','SSB15','SSB16','SSB17','SSB18','SSB19','C15','C16','C17','C18','SSB % change')
  
  xtab(tabXXa[,c(1:9,10:13)], caption='Table XX. lower 5% quantiles where applies', cornername='Description\\Variable', 
       file=paste(stamp,'_tabXX.html',sep=''), dec=c(2,2,2,2,0,0,0,0,0,0,0,0,0))
  


}