### The function for estimating parameters of GLM or GSSM
### For using GSSM, excecutable file of gamma_ss4.tpl is needed.
### The function for simulation using prediction models estimated in est.coef
### 

est.coef <-  function(ndata1,models=list(formula=c(expression(cbind(N1-1,max.N1-N1)~ factor(fy)+factor(fm)+pK),
                                   expression(mcatch.wt > 0 ~ log(pcatch.wt+1)+log(biomass.c)+log(biomass.p)+factor(fm)),
                                   expression(mcatch.wt ~ log(N1)+log(pcatch.wt+1)+log(biomass.c)+log(biomass.p)+factor(fm))),
                       func=list(glm,glm,glm,glm,glm),
                       is.stepAIC=c(TRUE,TRUE,TRUE,TRUE,TRUE)),
           is.Z=FALSE,N1.order=2,x4.init=NULL,
           admb.mode=list(is.admb=rep(FALSE,5)),data.N1=NULL){

  if(is.null(models$func)) models$func <- list(glm,glm,glm,glm,glm)    
  library(MASS)

  # estimate parameters
  model.res <- list()
  # ship model
  if(is.null(data.N1)) data.N1 <- subset(ndata1, mcatch.wt > 0 & Z==0)
  model.res[[1]] <- models$func[[1]](eval(models$formula[1]),
                                     family=binomial,
                                     data=data.N1)
  if(admb.mode$is.admb[1]==TRUE){
    ytime <- data.N1$xtime
    xtime <- ndata1$xtime
    if(N1.order==2){
      admb.init1 <- c(model.res[[1]]$coef,0,
                      2.2,1.1,-1,rep(0,length(xtime))) # gamma & sigma1                      
    }
    else{
      if(N1.order==3){
        admb.init1 <- c(model.res[[1]]$coef,0,
                        2,1,0,-1.5,rep(0,length(xtime))) # gamma & sigma1        
      }
      else{
        if(N1.order==4){
          if(is.null(x4.init)) x4.init <- rep(0,length(xtime))
          admb.init1 <- c(model.res[[1]]$coef,0,
                          1,1,0,0,-1.5,x4.init) # gamma & sigma1                  
        }
      }
    }
    write(admb.init1,file="init.pin")
    model.res[[1]] <- glm.k(eval(models$formula[1]),
                            data=data.N1,phase=admb.mode$phase[[1]],
                            offset=NULL, 
                            dist=admb.mode$dist[[1]], ytime=ytime,  xtime=xtime,logit=TRUE,
                            pinfile="init.pin")
    model.res[[1]]$formula <- models$formula[1]
    model.res[[1]]$data <- data.N1
    model.res[[1]]$admb.mode <- admb.mode    
  }

  # 0/1 model
  model.res[[2]] <- models$func[[2]](eval(models$formula[2]),
                                     family=binomial,
                                     data=data01 <- subset(ndata1, Z==0 & K==0))
  if(admb.mode$is.admb[2]==TRUE){
    ytime <- data01$xtime01
    xtime <- ndata1$xtime01
    admb.init2 <- c(model.res[[2]]$coef,0,
                    1.8,0.5,-1,rep(0,length(unique(xtime)))) # gamma & sigma1
    write(admb.init2,file="init.pin")    
    model.res[[2]] <- glm.k(eval(models$formula[2]),
                            data=data01,phase=admb.mode$phase[[2]],
                            offset=NULL,nx=length(unique(xtime)),
                            dist=admb.mode$dist[[2]], ytime=ytime,  xtime=sort(unique(xtime)),logit=TRUE,
                            pinfile="init.pin")
    model.res[[2]]$formula <- models$formula[2]
    model.res[[2]]$data <- data01
    model.res[[2]]$admb.mode <- admb.mode    
  }
  
  # catch model
  model.res[[3]] <- models$func[[3]](eval(models$formula[3]),
                                     family=Gamma(link="log"),
                                     data=tmpdat <- subset(ndata1, Z==0 & K==0 & mcatch.wt>0))                                 
  if(admb.mode$is.admb[3]==TRUE){
    ytime <- tmpdat$xtime
    xtime <- ndata1$xtime
    if(!is.Z){
      admb.init3 <- c(model.res[[3]]$coef,log(summary(model.res[[3]])$dispersion),
                      2,-2,-1.5,rep(0,length(xtime)))
    }
    else{
      tmpz <- log(summary(model.res[[3]])$dispersion)
      mm <- length(model.res[[3]]$coef)
      admb.init3 <- c(model.res[[3]]$coef[1:(mm-5)],rep(tmpz,6),
                      2,-2,-1.5,rep(0,length(xtime)))      
    }
    write(admb.init3,file="init.pin")    
    model.res[[3]] <- glm.k(eval(models$formula[3]),
                            data=tmpdat,phase=admb.mode$phase[[3]],
                            offset=NULL,is.Z=is.Z,        
                            dist=admb.mode$dist[[3]], ytime=ytime,  xtime=xtime,
                            pinfile="init.pin")
    model.res[[3]]$formula <- models$formula[3]
    model.res[[3]]$data <- tmpdat
    model.res[[3]]$admb.mode <- admb.mode
  }

  # closure A model
  model.res[[4]] <- models$func[[4]](eval(models$formula[4]),
                                     family=binomial,
                                     data=subset(ndata1,V==1 & Z==0))

  # closure B model
  model.res[[5]] <- models$func[[5]](eval(models$formula[5]), 
                                     family=binomial,
                                     data=subset(ndata1,V==0 & cup<Inf & Z==0))

  for(i in 1:5){
    if(!models$is.stepAIC[i]){
#      if(i==3 & !is.null(admb.mode)){
      model.res[[i]] <- model.res[[i]]
     }
    else{
      model.res[[i]] <- stepAIC(model.res[[i]])
    }
  }

  return(list(res.eff.bi=model.res[[1]],
              res.all.pa=model.res[[2]],
              res.all.po=model.res[[3]],
              under.k=model.res[[4]],
              over.k=model.res[[5]],
              ndata1=ndata1))
}


### The function for simulation using prediction models estimated in est.coef
###
make.simu2 <- function(ndata1,glm.objs,vpa.res,type="reg",# or "nreg" 
                       N=1000,pK.effect=FALSE,admb.mode=rep(FALSE,3),
                       teiji.effect=FALSE,seed=1,oldver=NULL,adhoc.catch=rep(1,100),
                       float.dyn=TRUE,n=nrow(ndata1),obs.effort=FALSE,
                       detail.res=TRUE,is.Z=FALSE,fix01=FALSE,
                       start.y=1,K.error=TRUE,random.K=NULL){
  set.seed(seed)
  if(type!="reg") ndata1$cup  <- Inf      
  y <- start.y
  uyear <- sort(unique(ndata1$fy))
  uyear <- c(uyear,max(uyear)+1)
  yn <- as.character(uyear)
  nage <- nrow(vpa.res$naa)
  f.day <- c(which(ndata1$fm==1 & ndata1$date==1),n+1)
  fff <- function(x) tapply(x,ndata1$fy,sum)

 
  # statistical model results
  res.eff.bi <- glm.objs$res.eff.bi
  res.all.pa <- glm.objs$res.all.pa
  res.all.po <- glm.objs$res.all.po
  under.k <- glm.objs$under.k
  over.k  <- glm.objs$over.k

  # vpa res
  vnaa <- vpa.res$naa[as.character(uyear)]
  vmaa <- vpa.res$maa[as.character(uyear)]
  vwaa <- vpa.res$waa[as.character(uyear)]
  vsaa <- vpa.res$saa[as.character(uyear)]
  vM <- vpa.res$M[as.character(uyear)]
  vbiom <- vpa.res$biom[as.character(uyear)]
  vrps <- as.numeric(vpa.res$rps[as.character(uyear)])
  ocatch <- tapply(ndata1$ocatch.wt,ndata1$fy,sum)

  pF <- numeric()
  for(i in 1:length(ocatch)){
    pF <- cbind(pF,exp(caa.est(vnaa[,i],vsaa[,i],vwaa[,i],vM[,i],ocatch[i],)$x) * vsaa[,i]) # partial F for other fishery
    catch.ratio <- rep(1,length(uyear))
  }
  if(!is.null(oldver)){
    pF[] <- 0
    catch.ratio <- oldver
  }
  
  # N at age
  naa <- array(0,dim=c(nage,length(uyear),N)) # by age, day, replication, nonreg or reg
  
  vb <- biom <- matrix(0,n+1,N) 
  naa[,1,] <- unlist(vnaa[,1])
  vb[1,] <- sum(vbiom[,1] * vsaa[,1])
  biom[1,] <- sum(vbiom[,1])

  # explanatery variables
  K <- catch <- N1 <-  mgcatch <- matrix(0,n,N)
  biomass.p <- biomass.c <- pK <- pW <- V <- D <- vpacatch <- ssb <- matrix(0,n+1,N)
  X2 <- matrix(0,n+2,N)
  V[1,] <- ndata1$V[1] # Does the catch previous day exceed threshold?
  X2[1,] <- ndata1$X2[1] # Does the catch two days before exceed threshold & no day closure previous day?
  X2[2,] <- ndata1$X2[2] # 
  pK[1,] <- ndata1$pK[1] # Is the previous day closed?
  K[1,] <- ndata1$K[1] # Is today closed?
  D[1,] <- ndata1$D[1] # Differences between actual catch and threshold
  pW[1,] <- ndata1$pW[1] # Is the catch previous day zero?

  # biomass
  biomass.c[1,] <- ndata1$biomass.c[1] 
  biomass.p[1,] <- ndata1$biomass.p[1]
  ssb[1,] <- colSums(sweep(naa[,1,],1,vwaa[,y]*vmaa[,y],FUN="*"))/1000


  # define dispersion parameter
  if(!admb.mode[3]){
    disp <- summary(res.all.po)$dispersion
  }
  else{
    disp <- NULL
  }

  start.n <- 1

  # calculate day by day
  for(i in start.n:n){
    cat(i," ")
    if(float.dyn==FALSE){
      biomass.c[i,] <- ndata1$biomass.c[i]
      biomass.p[i,] <- ndata1$biomass.p[i]
    }

    # create data for prediction
    edata1 <- data.frame(pW=pW[i,],
                         biomass.c=biomass.c[i,],
                         biomass.p=biomass.p[i,],
                         fy=rep(ndata1$fy[i],N),                         
                         fm=rep(ndata1$fm[i],N),
                         pK=pK[i,])
    edata1$pW2 <- edata1$pW-edata1$pK
   
    if(admb.mode[1]==FALSE){
      e.p.obs <- predict(res.eff.bi, newdata=edata1,se=FALSE,
                         type="response")
    }
    else{
      tmpfunc <- function(dat){
        # dummy data
        dat.dammy <- dat[1:12,]
        dat.dammy$fm[1:12] <- 1:12
        dat.dammy$fy[1:6] <- 2003:2008
        dat <- rbind(dat.dammy,dat)

        dat$N1 <- 1
        dat$max.N1 <- 10        
        dat$year <- 1                
        write(res.eff.bi$allpar,file="tmp.pin")
        ares <- glm.k(eval(res.eff.bi$formula),dat=dat,est=FALSE,pinfile="tmp.pin",
                      dist=res.eff.bi$admb.mode$dist[1],logit=TRUE,
                      ytime=rep(ndata1$xtime[i],nrow(dat)),xtime=as.numeric(names(res.eff.bi$x)),
                      offset=eval(res.eff.bi$admb.mode$offset[[1]]))
        return(ares$report[-1:-12,2]) # return except for dummy data(-1:-12)
      }
      e.p.obs <- tmpfunc(edata1)
    }

    if(obs.effort==FALSE){
      N1[i,] <- rbinom(N, ndata1$max.N1[i]-1, prob=e.p.obs)+1
      edata1$N1 <- N1[i,]
    }
    else{
      edata1$N1 <- N1[i,] <- ifelse(is.na(ndata1$N1[i]),1,ndata1$N1[i])
    }

    if(is.null(random.K)){
      if(ndata1$cup[i]<Inf){
        catch.V1 <- rbinom(N, V[i,], predict(under.k, data.frame(D=D[i,]), type="response"))
        catch.V2 <- rbinom(N, 1-V[i,], predict(over.k, data.frame(D=D[i,],pK=pK[i,],X2=X2[i,]), type="response"))
        if(K.error==TRUE){
          K[i,] <- ifelse(V[i,]==1,catch.V1,catch.V2)
        }
        else{
          K[i,] <- V[i,]
        }
      }
      else{
        K[i,] <- rep(0,N) 
      }
    }
    else{
        K[i,] <- rbinom(N,1,random.K)
      }

    # keep the predicted catch 
    ifcatch.obs <- rand.pred.catch(res.all.pa, res.all.po, edata1,
                                   fix01=fix01,catch01=ndata1$mcatch.wt[i],
                                   disp=disp,admb.mode=admb.mode,
                                   day=ndata1$xtime[i],is.Z=is.Z,day01=ndata1$xtime01[i])

    if(ndata1$Z[i]==0) { 
      catch[i,] <- ifelse(K[i,]==1, 0, ifcatch.obs ) 
    }
    else{
      catch[i,] <- 0
    }

    #-------------- Update for next calculation ------------------
    # If Cup < chub mackerel + spotted mackerel catch
    mgcatch[i,] <- catch[i,] + catch[i,] * ndata1$goma.ratio[i]  
    V[i+1,] <- ifelse(mgcatch[i,] > ndata1$cup[i], 1, 0)
    
    if(i>1) X2[i+1,] <- ifelse(mgcatch[i-1,] > ndata1$cup[i-1] & K[i,]==0, 1, 0)
    D[i+1,] <- abs(mgcatch[i,] - ndata1$cup[i])/ndata1$cup[i]
    
    D[i+1,is.na(D[i+1,])] <- 0
    if(pK.effect==TRUE){
      pK[i+1,] <- 0
    }
    else{
      pK[i+1,] <- K[i,]      
    }
    pW[i+1,] <- as.numeric(catch[i,]==0)

    #----------------- naa update -------------
    if((ndata1$date[i+1]==1 && ndata1$fm[i+1]==1 && i > 1)|i==n){
      acatch <- apply(catch, 2, fff)[y,] 

      tmp <- caa.est.mat(naa.mat=naa[,y,],saa=vsaa[,y],waa=vwaa[,y],M=vM[,y],
                         vcatch=rep(0,N),pF=pF[,y])       # For the first time, the other fishery catch the fish
      tmp <- caa.est.mat(naa.mat=naa[,y,],saa=vsaa[,y],waa=vwaa[,y],M=vM[,y],
                         vcatch=(acatch*catch.ratio[y]*adhoc.catch[y] + tmp$vpacatch),pF=0)   # Then, purse seine fishery catch the fish
      naa[,y+1,] <- tmp$naa
      vpacatch[i+1,] <- tmp$vpacatch
      naa[1,y+1,] <- colSums(sweep(naa[,y+1,],1,vwaa[,y+1]*vmaa[,y+1],FUN="*")) * vrps[y+1]
      
      biomass.c[i+1,]  <- biom[i+1,] <-  colSums(sweep(naa[,y+1,],1,vwaa[,y+1],FUN="*"))/1000
      ssb[i+1,] <- colSums(sweep(naa[,y+1,],1,vwaa[,y+1]*vmaa[,y+1],FUN="*"))/1000
      if(y > 1){
        biomass.p[i+1,] <- biomass.c[i-365,]
      }
      else{
        biomass.p[i+1,] <- sum(vpa.res$biom[as.character(min(uyear)-1)])/1000
      }
    }
    else{
      biomass.c[i+1,]  <- biomass.c[i,]
      if(y > 1){
        biomass.p[i+1,] <- biomass.c[i-365,]
      }
      else{
        biomass.p[i+1,] <- sum(vpa.res$biom[as.character(min(uyear)-1)])/1000
      }
    }
    
    if((ndata1$date[i+1]==1 && ndata1$fm[i+1]==1 && i> 1)|i==n) y <- y+1
  }

  acatch <- apply(catch, 2, fff)

  if(isTRUE(detail.res)){ # simple output
    invisible(list(acatch=acatch,abiom=biomass.c[f.day,],
                   assb=ssb[f.day,],biomass.c=biomass.c,glm.objs=glm.objs,
                   biomass.p=biomass.p,pK=pK,
                   catch=catch,N1=N1,K=K,f.day=f.day,vpacatch=vpacatch[f.day,],
                   vpa.bio=vpa.res$vbiom[yn]/1000,vpa.ssb=vpa.res$vssb[yn]/1000,
                   vpa.catch.true=vpa.res$acatch[yn]))
  }
  else{
    if(detail.res==2){ # output for model averaging
      invisible(list(biomass.c=biomass.c,pW=pW,biomass.p=biomass.p,pK=pK,naa=naa,ssb=ssb,
                     catch=catch,K=K,f.day=f.day,acatch=acatch))
    }
    else{ # the least output
      invisible(list(acatch=acatch,abiom=biomass.c[f.day,],glm.objs=glm.objs,
                     assb=ssb[f.day,],vpacatch=vpacatch[f.day,],K=K,
                     vpa.bio=vpa.res$vbiom[yn]/1000,vpa.ssb=vpa.res$vssb[yn]/1000,
                     vpa.catch.true=vpa.res$acatch[yn]))
    }
  }}

########################################################
### FUNCTIONS used inside the above functions
########################################################

expit <- function(x) 1/(1+exp(-x))

rand.pred.catch <- function(res.pa,res.po,dat,disp=NULL,admb.mode=FALSE,day=1,fix01=FALSE,catch01=NULL,
                            resample=NULL,is.Z=FALSE,day01=NULL){
  n <- nrow(dat)

  # dummy data for ADMB
  dat.dammy <- dat[1:12,]
  dat.dammy$fm <- 1:12
  dat.dammy$fy[1:6] <- 2003:2008
  dat <- rbind(dat.dammy,dat)
  dat$catch.wt <- 1
  dat$mcatch.wt <- 1
  dat$year <- 1  

  # 0/1 catch
  if(admb.mode[2]==FALSE){
    if(!fix01){
      pred.c <- rbinom(n, 1, expit(predict(res.pa, newdata=dat)))
    }
    else{
      pred.c <- rep(as.numeric(catch01>0),n)
    }
  }
  else{
    write(res.pa$allpar,file="tmp.pin")
    pred.c <- glm.k(eval(res.pa$formula),dat=dat,est=FALSE,pinfile="tmp.pin",
                  dist=res.pa$admb.mode$dist[2],logit=TRUE,
                  ytime=rep(day01,nrow(dat)),xtime=as.numeric(names(res.pa$x)),
                  offset=eval(res.po$admb.mode$offset[[2]]),is.Z=FALSE)
    pred.c <- rbinom(n, 1, pred.c$report[-1:-12,2])
  }

  # positive catch
  if(admb.mode[3]==FALSE){
    m <- exp(predict(res.po, newdata=dat))
    if(is.null(disp)){
      disp <- summary(res.po)$dispersion
    }  
  }
  else{
    write(res.po$allpar,file="tmp.pin")
    ares <- glm.k(eval(res.po$formula),dat=dat,est=FALSE,pinfile="tmp.pin",
                  dist=res.po$admb.mode$dist[3],
                  ytime=rep(day,nrow(dat)),xtime=as.numeric(names(res.po$x)),
                  offset=eval(res.po$admb.mode$offset[[3]]),is.Z=is.Z)
    if(!is.Z){
      disp <- exp(ares$pars["dispersion"])
    }
    else{
      disp <- ares$pars[seq(x <- which(names(ares$pars)=="dispersion"),x+5)]
      disp <- exp(ifelse(dat$fy==2003,disp[1],disp[dat$fy-2002] + disp[1]))
    }
    m <- ares$report[-1:-12,2] # remove dummy data
  }

  v <- disp*m^2
  scale.p <- ifelse(m==0, 0.000001, v/m)
  shape.p <- m/scale.p
  pred.n <- ifelse(m > 0 & !is.na(m), rgamma(n, shape=shape.p, scale=scale.p), 0)
  shape.p <- m/scale.p    
  pred.n <- rgamma(n, shape=shape.p, scale=scale.p)  

  pred <- pred.c*pred.n
  pred
}

solv.Feq <- function(cvec,nvec,mvec){

  Fres <- rep(0,length(cvec))
 # cat(nvec," ")
  for(i in 1:length(cvec)){
    F0 <- cvec[i]/nvec[i]
    F1 <- cvec[i]*(F0+mvec[i])/nvec[i]/(1-exp(-F0-mvec[i]))
    if(round(cvec[i],6)<round(nvec[i],6)){
      while(abs(F0-F1)>0.0001 ){
        F0 <- F1
        F1 <- cvec[i]*(F0+mvec[i])/nvec[i]/(1-exp(-F0-mvec[i]))
        if(F0-F1==-Inf) cat("\n",cvec[i]," ",nvec[i]," \n")
      }
      Fres[i] <- F1
    }
    else{
      Fres[i] <- 10
      cat("Warning: catch exceeded tot_num at: ",i," ",
          round(cvec[i],6)," ",round(nvec[i],6),"\n")
      }
  }
  Fres
}

calc.rel.abund <- function(sel,Fr,na,M,waa,maa,max.age=Inf,masaba=FALSE){
  rel.abund <- rep(NA, na)
  rel.abund[1] <- 1
  for (i in 2:(na-1)) {
    rel.abund[i] <- rel.abund[i-1]*exp(-M[i-1]-sel[i-1]*Fr)
  }
  rel.abund[na] <- rel.abund[na-1]*exp(-M[na-1]-sel[na-1]*Fr)*(1-exp(-(max.age-(na-2))*(M[na]+sel[na]*Fr)))/(1-exp(-M[na]-sel[na]*Fr))

  ypr1 <- rel.abund*waa*(1-exp(-sel*Fr))*exp(-M/2)

  # tentative option to set masaba ------
  ypr2 <- numeric()
  for(i in 1:(na-1)){
    ypr2[i] <- (rel.abund[i]-rel.abund[i+1])*Fr*sel[i]/(Fr*sel[i]+M[i])
  }
  ypr2[na] <- (rel.abund[na])*Fr*sel[i]/(Fr*sel[i]+M[i])
  ypr2 <- ypr2*waa
  #--------------------------------------
  
  if(masaba==TRUE) ypr <- ypr2
  else ypr <- ypr1
#  browser()
  spr <- rel.abund*waa*maa
  return(list(rel.abund=rel.abund,ypr=ypr,spr=spr))
}


caa.est.mat <- function(naa.mat,vcatch,pF,...){
  naa.mat2 <- naa.mat
  vpacatch <- numeric()
  for(i in 1:length(vcatch)){
    tmp <- caa.est(naa=naa.mat[,i],catch.obs=vcatch[i],pF=pF,...)
    naa.mat2[,i] <- tmp$naa
    vpacatch[i] <- tmp$vpacatch    
  }
  return(list(naa=naa.mat2,vpacatch=vpacatch))
}

caa.est <- function(naa,saa,waa,M,catch.obs,plus.age=TRUE,pF=0){
    tmpfunc <- function(lx,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,out=FALSE){
    x <- exp(lx)
    caa <- naa*(1-exp(-saa*x))*exp(-M/2)     # Pope
    wcaa <- caa*waa

    if(out==FALSE){
      return((sum(wcaa)/1000-catch.obs/1000)^2)
    }
    else{
      return(caa)
    }
  }

  tmp <- optimize(tmpfunc,c(-20,10),catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,
                  M=M,out=FALSE,tol=0.00000001)
  tmp2 <- tmpfunc(lx=tmp$minimum,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,out=TRUE)

  faa <- exp(tmp$minimum) * saa
  naa2 <- update.number(naa,faa+pF,M,plus.age=plus.age)
  vpacatch <- sum(naa*(1-exp(-faa-pF))*exp(-M/2)*waa)
  return(list(x=tmp$minimum,caa=tmp2,naa=naa2,vpacatch=vpacatch))
}

update.number <- function(naa,faa,M,plus.age=TRUE){
  nage <- length(naa)
  naa2 <- rep(0,nage)  
  if(plus.age==TRUE){
    naa2[2:nage] <- naa[1:(nage-1)] * exp(-M[1:(nage-1)]-faa[1:(nage-1)])
    naa2[nage] <- naa2[nage]+ naa[nage]*exp(-M[nage]-faa[nage])
  }
  else{
    naa2 <- naa * exp(-M-faa)
  }
  naa2
}

### The function for graphical output
show.res <- function(x1,x2=NULL,ndata1=NULL,pick=TRUE,yscale=c(1,1,1,1)){
  x <- t(x1$vpacatch[-1,pick])
  boxplot(x,ylim=c(0,max(x1$vpacatch[,pick])*yscale[1]),ylab="VPA catch",boxwex=0.3,at=1:ncol(x))
  if(!is.null(x2))boxplot(t(x2$vpacatch[-1,pick]),add=T,boxwex=0.3,col=2,at=1:ncol(x)+0.3)
  points(x1$vpa.catch.true[1:ncol(x)],col=3,type="b",lwd=3)

  x <- t(x1$acatch[,pick])
  boxplot(x,ylim=c(0,max(x1$acatch)*yscale[2]),ylab="m catch",boxwex=0.3,at=1:ncol(x))
  if(!is.null(x2))boxplot(t(x2$acatch[,pick]),add=T,boxwex=0.3,col=2,at=1:ncol(x)+0.3)
  if(!is.null(ndata1)){
    tmp <- tapply(ndata1$mcatch.wt,ndata1$fy,sum)
    points(1:ncol(x),tmp[names(tmp)%in%colnames(x)],col=3,type="b",lwd=2)
  }

  x <- t(x1$abiom[,pick])
  boxplot(x,ylim=c(0,max(x1$abiom)*yscale[3]),ylab="Biomass",boxwex=0.3,at=1:ncol(x))
  if(!is.null(x2))boxplot(t(x2$abiom[,pick]),add=T,boxwex=0.3,col=2,at=1:ncol(x)+0.3)
  points(x1$vpa.bio,col=3,type="b",lwd=3)
  x <- x1$abiom/x2$abiom
  title(paste("Percentage ratio=",round(median(x[nrow(x),]),2)))

  x <- t(x1$assb[,pick])
  boxplot(x,ylim=c(0,max(x1$assb)*yscale[4]),ylab="ssb",boxwex=0.3,at=1:ncol(x))
  if(!is.null(x2))boxplot(t(x2$assb[,pick]),add=T,boxwex=0.3,col=2,at=1:ncol(x)+0.3)
  x <- x1$assb/x2$assb
  points(x1$vpa.ssb,col=3,type="b",lwd=3)
  title(paste("Percentage ratio=",round(median(x[nrow(x),]),2)))  
}


glm.k <- function(formula, data, offset, hessian = FALSE, phase = NULL, inits.x = NULL,
                  inits.z=NULL,inits.w=NULL,inits.phi=NULL,inits.psi=NULL,invis = TRUE,
                  contrasts=NULL,  dist = "gamma",is.Z=FALSE, est=TRUE,use.parfile=FALSE,
                  logit=FALSE,nx=length(xtime),weight=1,
                  pinfile=NULL,ytime=NULL,xtime=NULL)
{

  if (is.null(phase)){
    np <- 6   # 1st for x, 2nd for z, 3-5th for random effects
    phase <- rep(1,np)
  }

# initial settings

  N <- nrow(data)
  mf <- match.call(expand.dots = FALSE)
  m <- match(c("formula",  "data", "offset"), names(mf), 0)
  mf <- mf[c(1L, m)]
  mf$drop.unused.levels <- TRUE
  ff <- formula
  if(is.Z==TRUE){  
    formula[[3]][1] <- call("+")
    mf$formula <- formula
    ffc <- . ~ .
    ffz <- ~.
    ffc[[2]] <- ff[[2]]
    ffc[[3]] <- ff[[3]][[2]]
    ffz[[3]] <- ff[[3]][[3]]
    ffz[[2]] <- NULL
    mf[[1L]] <- as.name("model.frame")
    mf <- eval(mf, parent.frame())
  
    mtX <- terms(ffc, data = data)
    X <- model.matrix(mtX, mf, contrasts)
    mtZ <- terms(ffz, data = data)
    Z <- model.matrix(mtZ, mf, contrasts)
    Y <- model.response(mf, "any")
    npx <- ncol(X)
    npz <- ncol(Z)
  }
  else{
    ffc <- ff
    mf[[1L]] <- as.name("model.frame")
    mf <- eval(mf, parent.frame())
    mtX <- terms(ffc, data = data)
    X <- model.matrix(mtX, mf, contrasts)
    Y <- model.response(mf, "any")
    npx <- ncol(X)
    npz <- 1
  }
  
  offset <- as.vector(model.offset(mf))
   
  if (is.null(offset)) offset <- rep(0,N)
   
  # data generation for ADMB
  
  file.name <- dist
  system(paste2("mv ",file.name,".par"," back.par"))
  system(paste2("mv ",file.name,".rep"," back.rep"))
  
  npar <- c(npx, npz)
  cat("# model type", file = paste(file.name,".dat",sep=""), sep ="\n")
  cat(ifelse(logit==TRUE,2,1),"\n",file = paste(file.name,".dat",sep=""), sep =" ", append=T)  
  
  cat("# number of observations and parameters", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)
  
  cat(c(N,npar), file = paste(file.name,".dat",sep=""), sep =" ", append=T)
  
  cat(file = paste(file.name,".dat",sep=""), sep="\n",append=T)  
  
  cat("# phase setting", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)
  cat("# X Z beta1 beta2 sigma_x x", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)  
  cat(phase, file = paste(file.name,".dat",sep=""), sep =" ", append=T)
  
  cat(file = paste(file.name,".dat",sep=""), sep="\n",append=T)  
    
  cat("# response variable", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)
  
  write.table(Y, file = paste(file.name,".dat",sep=""),append=T,sep=" ",col.names=F, row.names=F)
  
  cat(file = paste(file.name,".dat",sep=""),sep="\n",append=T)  
  
  cat("# eplanatory variables for mu", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)
  
  write.table(X,file = paste(file.name,".dat",sep=""),append=T, row.names=F, col.names=F)

  if(is.Z==TRUE){
    write.table(Z,file = paste(file.name,".dat",sep=""),append=T, row.names=F, col.names=F)
  }
  else{
    write.table(rep(1,N),file = paste(file.name,".dat",sep=""),append=T, row.names=F, col.names=F)    
  }
 
  cat("# offset variable", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)
  
  cat(offset, file = paste(file.name,".dat",sep=""), append=T, sep=" ")

  if(is.null(ytime)){
    ytime <- 1:length(Y)
    xtime <- 1:length(Y)
  }
  cat("\n# obs time", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)    
  write.table(ytime,file = paste(file.name,".dat",sep=""),append=T, row.names=F, col.names=F)

  cat("\n# x time length & data", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)    
  write.table(length(xtime),file = paste(file.name,".dat",sep=""),append=T, row.names=F, col.names=F)
  write.table(xtime,file = paste(file.name,".dat",sep=""),append=T, row.names=F, col.names=F)
  cat("\n# length of random effects", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)
  write.table(nx,file = paste(file.name,".dat",sep=""),append=T, row.names=F, col.names=F)
  cat("# initial value of x1 and x2", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)
  cat(c(0,0), file = paste(file.name,".dat",sep=""), sep ="\n", append=T)    
  cat("# for debuging", file = paste(file.name,".dat",sep=""), sep ="\n", append=T)
  cat(999, file = paste(file.name,".dat",sep=""), sep ="\n", append=T)  

# initial values
  pinfile0 <- paste(file.name,".pin",sep="")
  if(is.null(pinfile)){
    if(use.parfile==FALSE){
      if(is.null(inits.x)) inits.x <- rep(0,npx)
      cat("# initial parameters for PX", file = pinfile0, sep ="\n")
      cat(inits.x, file = pinfile0,append=T,sep=" ")
    cat(file = pinfile0,sep="\n",append=T)
      if(is.null(inits.z)) inits.z <- rep(0,npz)
      cat("# initial parameters for Pz", file = pinfile0, sep ="\n", append=T)
      cat(inits.z, file = pinfile0,append=T,sep=" ")
      cat(file = pinfile0,sep="\n",append=T)
    }
    else{
      file.copy(from=paste(file.name,".par",sep=""),to=pinfile0,overwrite=TRUE)
    }
  }
  else{
      file.copy(from=pinfile,to=pinfile0,overwrite=TRUE)    
  }
 
# model run
  if(hessian == FALSE) hes <- " -nohess" else hes <- ""
  if(est==FALSE) hes <- paste(hes," -maxfn 0")
  
  if(.Platform$OS.type=="unix"){
    system(paste("./",file.name, " -mno 15000 ",hes,sep=""))
  }
  else{
    shell(paste(file.name, ".exe -mno 15000",hes,sep=""),invisible=invis)
  }
 
# outputs
  res <- list()
  if(file.exists(paste(file.name,".par",sep=""))){
    p1 <- c(na.omit(as.numeric(scan(paste(file.name,".par",sep=""),what="",quiet=TRUE)[-(1:16)])))
    p2 <- scan(paste(file.name,".par",sep=""),what="character",quiet=TRUE)
    np <- length(p1)
    og <- scan(paste(file.name,".par",sep=""),what="",n=16,quiet=TRUE)
    obj.f <- as.numeric(og[[11]])
    grad.f <- as.numeric(og[[16]])
    res$allpar <- p1
    names(res$allpar) <- c(colnames(X),rep("dispersion",npz),"beta1","beta2","sigma1",xtime)
    res$x <- as.numeric(p2[(which(p2=="x:")+1):length(p2)])
    names(res$x) <- 1:length(res$x)
    res$pars <- res$allpar[1:(which(as.numeric(p2[(which(p2=="x:")+1)])==p1)-1)]
    res$objective <- obj.f
    res$max.gradient <- grad.f
    res$aic <- 2*obj.f+2*sum(npar)
    res$report <- read.table(paste(file.name,".rep",sep=""))
    if (hessian == TRUE) {
      cor0 <- NULL
      for (i in 1:np){
        cor1 <- scan(paste(file.name,".cor",sep=""),what="",nline=1,skip=i+1,quiet=TRUE)
        cor1 <- c(as.numeric(cor1[-c(1:2)]),rep(0,np-i))
        cor0 <- rbind(cor0,cor1)
      }
      std.error <- cor0[,2]
      res$std.error <- as.numeric(std.error)
      correl.mat <- as.matrix(cor0[,-(1:2)],rownames=1:np)
      row.names(correl.mat) <- NULL
      res$correl.mat <- as.dist(correl.mat)
    }
    res$convergence <- "TRUE"
  }
  else{
    res$convergence <- "FALSE"
  }
  return(res)
}

paste2 <- function (x, ...) 
{
    paste(x, sep = "", ...)
}

make.gmat <- function(ndata){
  ndata$ym <- paste(ndata$fy,ndata$fm,sep="-")
  x <- tapply(ndata$gcpue,ndata$ym,function(x)x[!is.na(x)])
  n <- lapply(x,length)
  gmat <- matrix(NA,length(ndata$gcpue),max(unlist(n)))
  for(i in 1:nrow(ndata)){
    tmp <- x[[which(names(x)==ndata$ym[i])]]        
    if(length(tmp)>0){
      gmat[i,1:length(tmp)] <- tmp
    }
    else{
      gmat[i,1] <- 0
    }
  }
  return(gmat)
}


###################################### 
### FUNCTIONS for outputs
######################################

make.table <- function(res0a,res0n,N=10000,func=median){
  res0a$E <- res0a$vpacatch[-1,]/res0a$abiom[1:5,]/1000
  res0n$E <- res0n$vpacatch[-1,]/res0n$abiom[1:5,]/1000

  set.seed(1)
  TBR <- sample(res0a$abiom[6,],10000,replace=TRUE)/sample(res0n$abiom[6,],10000,replace=TRUE)
  TER <- sample(apply(res0a$E,2,mean),10000,replace=TRUE)/sample(apply(res0n$E,2,mean),10000,replace=TRUE)
  TCR <- sample(apply(res0a$vpacatch[-1,],2,mean),10000,replace=TRUE)/sample(apply(res0n$vpacatch[-1,],2,mean),10000,replace=TRUE)    
  
  probs <- numeric()

  probs[c(1,2)] <- quantile(TCR,probs=c(0.05,0.95),na.rm=T)
  probs[c(3,4)] <- quantile(TER,probs=c(0.05,0.95),na.rm=T)
  probs[c(5,6)] <- quantile(TBR,probs=c(0.05,0.95),na.rm=T)

  med <- numeric()
  med[1] <- median(TCR,na.rm=T)
  med[2] <- median(TER,na.rm=T)
  med[3] <- median(TBR,na.rm=T)
  K <- apply(res0a$K,2,sum)
    
  return(c(med,probs,
           median(K),quantile(K,probs=c(0.05,0.95))))
}


show.res <- function(x1,x2=NULL,ndata1=NULL,pick=TRUE){
  x <- t(x1$acatch[,pick])/1000
  boxplot(x,ylim=c(0,max(x)),ylab="Annual catch of chub mackerels (1000 tons)",xlab="Year",boxwex=0.3,at=1:ncol(x),xaxt="n",pch=20,cex=0.5)
  if(!is.null(x2))boxplot(t(x2$acatch[,pick])/1000,add=T,boxwex=0.3,col="gray",at=1:ncol(x)+0.3,xaxt="n",pch=20,cex=0.5)
  axis(side=1,at=1:ncol(x)+0.15,labels=2004:2008)
  #points(x1$vpa.catch.true,col=2,type="b")
  if(!is.null(ndata1)){
    tmp <- tapply(ndata1$mcatch.wt,ndata1$fy,sum)
    points(1:ncol(x),tmp[names(tmp)%in%colnames(x)]/1000,col="gray",type="p",lwd=1,pch=4,cex=1.5)
  }

  x1$E <- x1$vpacatch[-1,]/x1$abiom[1:5,]/1000
  x2$E <- x2$vpacatch[-1,]/x2$abiom[1:5,]/1000  
  x <- t(x1$E[,pick])
  boxplot(x,ylim=c(0,max(x2$E)),ylab="Exploitation rates",boxwex=0.3,at=1:ncol(x),xaxt="n",xlab="Year",pch=20,cex=0.5)
  axis(side=1,at=1:ncol(x)+0.15,labels=2004:2008)    
  if(!is.null(x2))boxplot(t(x2$E[,pick]),add=T,boxwex=0.3,col="gray",at=1:ncol(x)+0.3,xaxt="n",pch=20,cex=0.5)

  tmp <- x1$vpa.catch.true
  tmp <- as.numeric(tmp[as.character(2004:2008)])/colSums(vpa.res$biom[as.character(2004:2008)])
  points(tmp,col="gray",type="p",lwd=1,pch=4,cex=1.5)

  x <- t(x1$abiom[,pick])
  boxplot(x,ylim=c(0,max(x1$abiom)),ylab="Biomass",boxwex=0.3,at=1:ncol(x),xaxt="n",xlab="Year",pch=20,cex=0.5)
  axis(side=1,at=1:ncol(x)+0.15,labels=2004:2009)  
  if(!is.null(x2))boxplot(t(x2$abiom[,pick]),add=T,boxwex=0.3,col="gray",at=1:ncol(x)+0.3,xaxt="n",pch=20,cex=0.5)
  points(as.numeric(x1$vpa.bio[as.character(2004:2009)]),col="gray",type="p",lwd=1.5,pch=4,cex=1.5)
  x <- x1$abiom/x2$abiom
  title(paste("kouka=",k1 <- round(median(x[nrow(x),]),2)),line=3)

  x <- t(x1$assb[,pick])
  boxplot(x,ylim=c(0,max(x1$assb)),ylab="Spawning biomass",boxwex=0.3,at=1:ncol(x),xaxt="n",xlab="Year",pch=20,cex=0.5)
  axis(side=1,at=1:ncol(x)+0.15,labels=2004:2009)    
  if(!is.null(x2))boxplot(t(x2$assb[,pick]),add=T,boxwex=0.3,col="gray",at=1:ncol(x)+0.3,xaxt="n",pch=20,cex=0.5)
  x <- x1$assb/x2$assb
  points(as.numeric(x1$vpa.ssb[as.character(2004:2009)]),col="gray",type="p",lwd=1,pch=4,cex=1.5)
  title(paste("kouka=",k2 <- round(median(x[nrow(x),]),2)),line=2)
  return(c(k1,k2))
}