# Script to generate plots from the state-space assessment model
# 
# Anders Nielsen <anders@nielsensweb.org> Aug. 2012
oldwd<-getwd()
 
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 ##############################

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


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

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

  
  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.base, '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')
    for(i in 1:length(RETRO)){
      addlines(RETRO[[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(RETRO)){
      addlines(RETRO[[i]], 'fbar', col=i+2, with.dot=TRUE, ci=FALSE)
    }
    stampit()
    setcap()

    # Recruit plot
    baseplot(fit.base, '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)
    }
    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.base, '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))
}

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:3
    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+3))
    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))
    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+3))
    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))
    lines(py,ssb[,1]/1000, col='blue', lwd=2)
    lines(py,ssb[,3]/1000, col='blue', lty='dashed', lwd=2)
    lines(py,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+3))
    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]))
    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+3))
    abline(v=py[1], col=grey(.8), lwd=3)
    addlines(fit.current, 'logCatch', col='red', trans=function(x)exp(x)/1000, lwd=2)
    catch<-rbind(summ(forecast$last.catch.sim), summ(forecast$s1.catch.sim), summ(forecast$s2.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))
    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))[3,-1]
    ssb<-cbind(summ(fore$last.ssb.sim),summ(fore$s1.ssb.sim),summ(fore$s2.ssb.sim),summ(fore$s3.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))[3,-1]
    land<-cbind(summ(fore$last.land.sim),summ(fore$s1.land.sim),summ(fore$s2.land.sim),summ(fore$s3.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))[3,-1]
    R<-cbind(summ(exp(fore$last.state.sim[,1])),
                 summ(exp(fore$s1.state.sim[,1])),
                 summ(exp(fore$s2.state.sim[,1])),
                 summ(exp(fore$s3.state.sim[,1])))[3,-1]
    X<-c(fbar,ssb,catch,land,fbarl,R)
    return(X)
  }
  
  tabX<-do.call(rbind,lapply(forecast, foreTab))
  
  colnames(tabX)<-c('Fbar11','Fbar12','Fbar13','SSB11','SSB12','SSB13','C11','C12','C13',  
                    'L11','L12','L13','Fbar11(L)','Fbar12(L)','Fbar13(L)','R11','R12','R13')
  rownames(tabX)<-foretab[,ncol(foretab)]
  
  tabXb<-cbind(tabX[,10:12],tabX[,7:9]-tabX[,10:12],tabX[,13:15], (tabX[,11]-oldtac)/oldtac*100)
  colnames(tabXb)<-c('L11','L12','L13','D11','D12','D13','Fbar11(L)','Fbar12(L)','Fbar13(L)','TAC % change')
  rownames(tabXb)<-foretab[,ncol(foretab)]
  
  tabXa<-tabX[,1:9]
  tabXa<-cbind(tabXa,(tabXa[,6]-tabXa[,5])/tabXa[,5]*100)
  colnames(tabXa)<-c('Fbar11','Fbar12','Fbar13','SSB11','SSB12','SSB13','C11','C12','C13','SSB % change')
  
  xtab(tabXa[,c(1:6,10,7:9)], 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,0,0,0,0,0,0,0))
  
  tabXb<-tabXb[,c(1:3,4:9)]
  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,2,2,2))
  
  mult<-1.0#mean(exp(fit.current$logscale[nrow(fit.current$logscale)-0:2,1]))
  tabXc<-tabXb[,c(1:3,4:6)]/mult
  tabXc<-cbind(tabXc,tabX[,16:18], (tabXc[,2]-oldtac)/oldtac*100)
  colnames(tabXc)<-c('L11','L12','L13','D11','D12','D13','R11','R12','R13','TAC % change')
  xtab(tabXc[,c(1:3,10,4:9)], 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))

}  
setwd(oldwd)