##--------------------------------------------------------
## CMSY analysis with estimation of total biomass, including Bayesian Schaefer 
## written by Rainer Froese with support from Gianpaolo Coro and Henning Winker
## Version 44 of 05.02.2015
##--------------------------------------------------------
library(R2jags)  # Interface with JAGS
library(coda) 

#-----------------------------------------
# Some general settings
#-----------------------------------------
# set.seed(999) # use for comparing results between runs
rm(list=ls(all=TRUE)) # clear previous variables etc
options(digits=3) # displays all numbers with three significant digits as default
graphics.off() # close graphics windows from previous sessions

#-----------------------------------------
# General settings for the analysis
#-----------------------------------------
sigR         <- 0.05 # overall process error for CMSY; 0.05 works reasonable for simulations, 0.02 for real data; 0 if deterministic model
n            <- 10000 # initial number of r-k pairs
ni           <- 3 # iterations for r-k-startbiomass combinations, to test different variability patterns; no improvement seen above 3
mgraphs      <- T # set to TRUE to produce additional graphs for management
write.output <- T # set to true if table of output is wanted
FutureCrash  <- "No" # initialize variable
FullSchaefer <- F    # initialize variable
qyr          <- c(NA,NA) # c(1985,2000) #  set year range with stable catch/biomass ratio, min 5 years; else set NA


#-----------------------------------------
# Start output to screen
#-----------------------------------------
cat("-------------------------------------------\n")
cat("Catch-MSY Analysis,", date(),"\n")
cat("-------------------------------------------\n")
#------------------------------------------
# Read data and assign to vectors
#------------------------------------------
 filename_1  <- "SmallTuna_Catch_2015.csv"
 filename_2  <- "SmallTuna_ID_2015.csv"
# filename_1  <- "ATUN_Catch.csv"
# filename_2  <- "ATUN_ID.csv"


# Define file names for output
outfile      <- paste("Out_",format(Sys.Date(),format="%B%d%Y_"),filename_2,sep="")
# create headers for data table file
if(write.output==T){
  outheaders = data.frame("SciName","Stock","StartYear","EndYear","start_k_lo","start_k_hi",
                          "MeanCatch","MaxCatch","LastCatch","MSY_S","lcl","ucl","r_S","lcl","ucl",
                          "k_S","lcl","ucl","q.1","q.2","q_S","lcl","ucl","r_CMSY","lcl","ucl",
                          "k_CMSY","lcl","ucl","MSY_CMSY","lcl","ucl",
                          "LastBio","LastBio_CMSY","2.5th","25th","97.5th","NextBio_CMSY","2.5th","25th","97.5th")
  write.table(outheaders,file=outfile, append = T, sep=",",row.names=F,col.names=F)
}
outfile.txt <- paste(outfile,".txt", sep="")

# Read data
cdat         <- read.csv(filename_1, header=T, dec=".", stringsAsFactors = FALSE)
cinfo        <- read.csv(filename_2, header=T, dec=".", stringsAsFactors = FALSE)
cat("Files", filename_1, ",", filename_2, "read successfully","\n")

#---------------------------------
# Select stocks to be analyzed
#---------------------------------
# stocks       <- sort(as.character(cinfo$stock)) # All stocks 
 stocks       <-  "WAH" # "BRS" "KGM"  "SMA" "DOL" "FRI" "POR" "BON_ATL" "BON_MED" "LTA_ATL"

# analyze one stock after the other
for(stock in stocks) {
  cat("Processing",stock,",", as.character(cinfo$ScientificName[cinfo$stock==stock]),"\n")
  # assign data from cinfo to vectors
  res          <- as.character(cinfo$Resilience[cinfo$stock==stock])
  StartYear    <- as.numeric(cinfo$StartYear[cinfo$stock==stock])
  EndYear      <- as.numeric(cinfo$EndYear[cinfo$stock==stock])
  r_low        <- as.numeric(cinfo$r_low[cinfo$stock==stock])
  r_hi         <- as.numeric(cinfo$r_hi[cinfo$stock==stock])
  stb_low      <- as.numeric(cinfo$stb_low[cinfo$stock==stock])
  stb_hi       <- as.numeric(cinfo$stb_hi[cinfo$stock==stock])
  intyr        <- as.numeric(cinfo$intyr[cinfo$stock==stock])
  intbio_low   <- as.numeric(cinfo$intbio_low[cinfo$stock==stock])
  intbio_hi    <- as.numeric(cinfo$intbio_hi[cinfo$stock==stock])
  endbio_low   <- as.numeric(cinfo$endbio_low[cinfo$stock==stock])
  endbio_hi    <- as.numeric(cinfo$endbio_hi[cinfo$stock==stock])
  Btype        <- as.character(cinfo$Btype[cinfo$stock==stock])
  FutureCrash  <- as.character(cinfo$FutureCrash[cinfo$stock==stock])
  comment      <- as.character(cinfo$comment[cinfo$stock==stock])
  
  
  # extract data on stock
  yr           <- as.numeric(cdat$yr[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])
  ct.raw           <- as.numeric(cdat$ct[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])/1000  ## assumes that catch is given in tonnes, transforms to '000 tonnes
  if(Btype=="observed" | Btype=="CPUE" | Btype=="simulated") {
    bt <- as.numeric(cdat$TB[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])/1000  ## assumes that biomass is in tonnes, transforms to '000 tonnes
  } else {bt <- NA}
  if(filename_1=="SimCatchCPUE.csv") {
    bt_true <- as.numeric(cdat$TB_true[cdat$stock==stock])/1000  ## true biomass used in simulated CPUE stocks
  }
  nyr          <- length(yr) # number of years in the time series

# change catch to 3 years moving average 
# for first year use reported catch and for second year 2 year average
ma                <- function(x){filter(x,rep(1/3,3),sides=1)}
  ct              <- ma(ct.raw)
  not.na          <- which(is.na(ct.raw)==F)
  ct[not.na[1]]   <- ct.raw[not.na[1]]
  ct[not.na[2]]   <- (ct.raw[not.na[1]]+ct.raw[not.na[2]])/2
  
  # initialize vectors for viable r, k, bt, and all in a matrix
  mdat.all    <- matrix(data=vector(),ncol=2+nyr+1)
   
  #----------------------------------------------------
  # Determine initial ranges for parameters and biomass
  #----------------------------------------------------
  # initial range of r from input file
  if(is.na(r_low)==F & is.na(r_hi)==F) {
    start_r <- c(r_low,r_hi)
  } else {
    # initial range of r based on resilience
    if(res == "High") {
      start_r <- c(0.6,1.5)} else if(res == "Medium") {
        start_r <- c(0.2,0.8)}    else if(res == "Low") {
          start_r <- c(0.05,0.5)}  else { # i.e. res== "Very low"
            start_r <- c(0.015,0.1)} 
  }
  
# initial biomass range from input file
if(is.na(stb_low)==F & is.na(stb_hi)==F) {
  startbio <- c(stb_low,stb_hi)
} else {
  # us high biomass at start as default
  startbio <- c(0.4,0.8)
}

  MinYear     <- yr[which.min(ct.raw)]
  MaxYear     <- yr[which.max(ct.raw)]
  # use year and biomass range for intermediate biomass from input file
  if(is.na(intbio_low)==F & is.na(intbio_hi)==F) {
    intyr    <- intyr
    intbio   <- c(intbio_low,intbio_hi)
    # else if year of minimum catch is at least 3 years away from StartYear and EndYear of series, use min catch
  } else if((MinYear - StartYear) > 3 & (EndYear - MinYear) > 3 ) {
    # assume that biomass range in year before minimum catch was 0.01 - 0.4
    intyr    <- MinYear-1
    intbio <- c(0.01,0.4) 
    # else if year of max catch is at least 3 years away from StartYear and EndYear of series, use max catch  
  } else if((MaxYear - StartYear) > 3 & (EndYear - MaxYear) > 3 ) {
    # assume that biomass range in year before maximum catch was 0.3 - 0.9
    intyr    <- MaxYear-1
    intbio <- c(0.3,0.9)  
  } else {
    # assume uninformative range 0-1 in mid-year 
    intyr     <- as.integer(mean(c(StartYear, EndYear)))
    intbio    <- c(0,1) }
  # end of intbio setting
  
  # final biomass range from input file
  if(is.na(endbio_low)==F & is.na(endbio_hi)==F) {
    endbio   <- c(endbio_low,endbio_hi)
  } else {
    # else use smoothed Catch/maxCatch to estimate final biomass
    endbio  <- if(ct[nyr]/max(ct) > 0.5) {c(0.4,0.8)} else {c(0.01,0.4)} 
  } # end of final biomass setting

# initial range of k values, assuming min k will be larger than max catch / prior for r 
# and max catch will be at least MSY if prior mean biomass b/k < = 0.5
# and max catch will at least half of MSY otherwise; all divided by prior of r

  gm.start_r   <- exp(mean(log(start_r))) # get geometric mean of prior r range
  m.max.ct     <- sort(ct.raw,decreasing=T)[2] # median of 3 highest catches = 2nd highest catch
  if(mean(endbio) <= 0.5) {
    start_k <- c(m.max.ct/start_r[2],6*m.max.ct/gm.start_r)} else {
      start_k <- c(2*m.max.ct/start_r[2],8*m.max.ct/gm.start_r)} 

  #----------------------------------------------
  # MC with Schaefer Function filtering
  #----------------------------------------------
  Schaefer <- function(ri, ki, startbio, intyr, intbio, endbio, sigR, pt) {

    # if stock is not expected to crash within 3 years if last catch continues
    if(FutureCrash == "No") {
      yr.s     <- c(yr,EndYear+1,EndYear+2,EndYear+3)
      ct.s     <- c(ct,ct[yr==EndYear],ct[yr==EndYear],ct[yr==EndYear])
      nyr.s    <- length(yr.s)
    } else{
      yr.s     <- yr
      ct.s     <- ct
      nyr.s    <- nyr
    }
    
    # create vector for initial biomasses
    startbt     <- seq(from =startbio[1], to=startbio[2], by = (startbio[2]-startbio[1])/10)
    nstartbt    <- length(startbt)
    npoints     <- length(ri)

    # create vectors for viable r, k and bt
    mdat        <- matrix(data=NA, nrow = npoints*nstartbt, ncol = 2+nyr+1) 
    
    # get index of intermediate year
    intyr.i     <- which(yr.s==intyr) 
          
    #loop through r-k pairs
    for(i in 1 : npoints) {
      if (i%%1000==0)
        cat(".")
      
      # create empty vector for annual biomasses
      bt <- vector()
      # set flag for viable point to FALSE
      vp <- F
      
      # loop through range of relative start biomasses
      for(j in startbt) {    
        # set initial biomass, including process error
        bt[1]=j*ki[i]*exp(rnorm(1,0, sigR))  ## set biomass in first year
        
        # repeat test of r-k-startbt combination to allow for different random error
        for(re in 1:ni) {
        
        #loop through years in catch time series
        for(t in 1:nyr.s) {  # for all years in the time series
          xt=rnorm(1,0, sigR) # set new process error for every year  
          zt=rnorm(1,0, sigR/2) # set new observation error for every year, as half of process error
          
          # calculate biomass as function of previous year's biomass plus surplus production minus catch
          bt[t+1]=bt[t]+ri[i]*bt[t]*(1-bt[t]/ki[i])*exp(xt)-ct.s[t]*exp(zt) 
          
          # if biomass < 0.01 k or > 1.1 k, discard r-k-startbt combination
          if(bt[t+1] < 0.01*ki[i] || bt[t+1] > 1.1*ki[i]) { break } # stop looping through years, go to next upper level
          
          if ((t+1)==intyr.i && (bt[t+1]>(intbio[2]*ki[i]) || bt[t+1]<(intbio[1]*ki[i]))) { break }  #intermediate year check
          
        } # end of loop of years
            
        # if loop was broken or last biomass falls outside of expected ranges do not store results, go directly to next startbt
        if(t < nyr.s || bt[yr.s==EndYear] > (endbio[2]*ki[i]) || bt[yr.s==EndYear] < (endbio[1]*ki[i])) {
          next } else { 
            # store r, k, and bt, plot point, then go to next startbt
             mdat[((i-1)*nstartbt)+which(startbt==j),1]           <- ri[i]
             mdat[((i-1)*nstartbt)+which(startbt==j),2]           <- ki[i]
             mdat[((i-1)*nstartbt)+which(startbt==j),3:(2+nyr+1)] <- bt[1:(nyr+1)]/ki[i]
             vp <- T }
      
        } # end of repetition for random error
        
      } # end of j-loop of initial biomasses 
       # plot point to show progress
       if(pt==T & vp==T) points(x=ri[i], y=ki[i], pch=".", cex=4, col="gray")

    } # end of i-loop of r-k pairs 

    # return results
    mdat <- na.omit(mdat) # ignore rows filled with NA from initial matrix   
    cat("\n")
    return(list(mdat))
} # end of Schaefer function
    
  #------------------------------------------------------------------
  # Uniform sampling of the r-k space
  #------------------------------------------------------------------
  # get random set of r and k from log space distribution 
  ri1 = exp(runif(n, log(start_r[1]), log(start_r[2])))  
  ki1 = exp(runif(n, log(start_k[1]), log(start_k[2])))  
  
  #-----------------------------------------------------------------
  # Plot data and progress
  #-----------------------------------------------------------------
  windows(14,9)
  par(mfcol=c(2,3))
  # plot catch
  plot(x=yr, y=ct.raw, ylim=c(0,1.2*max(ct.raw)), type ="l", bty="l", main=paste(stock,"catch"), xlab="Year", 
       ylab="Catch", lwd=2)
  points(x=yr[which.max(ct.raw)], y=max(ct.raw), col="red", lwd=2)
  points(x=yr[which.min(ct.raw)], y=min(ct.raw), col="red", lwd=2)
  
  # plot r-k graph
  plot(x=ri1, y=ki1, xlim = start_r, ylim = start_k, log="xy", xlab="r", ylab="k", 
       main="Finding viable r-k", pch=".", cex=3, bty="l", col="gray95")
  
  #---------------------------------------------------------------------
  # 1 - Call CMSY-Schaefer function to preliminary explore the r-k space
  #---------------------------------------------------------------------
  cat("First Monte Carlo filtering of r-k space with ",n," points\n")
  MCA <-  Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T)
   mdat.all <- rbind(mdat.all,MCA[[1]])
   rv.all   <- mdat.all[,1]
   kv.all   <- mdat.all[,2]
   btv.all  <- mdat.all[,3:(2+nyr+1)]
   # count viable trajectories and r-k pairs 
   n.viable.b   <- length(mdat.all[,1])
   n.viable.pt <- length(unique(mdat.all[,1]))
   cat("Found ",n.viable.b," viable trajectories for", n.viable.pt," r-k pairs\n")
 
  #----------------------------------------------------------------------- 
  # 2 - if the lower bound of k is too high, reduce it by half and rerun
  #-----------------------------------------------------------------------
  if(length(kv.all[kv.all < 1.1*start_k[1] & rv.all < mean(start_r)]) > 10) {
   cat("Reducing lower bound of k, resampling area with",n,"additional points\n")
   start_k <- c(0.5*start_k[1],start_k[2])
   ri1 = exp(runif(n, log(start_r[1]), log(start_r[2])))  
   ki1 = exp(runif(n, log(start_k[1]), log(start_k[2])))  
   MCA <-  Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T)
   mdat.all <- rbind(mdat.all,MCA[[1]])
   rv.all   <- mdat.all[,1]
   kv.all   <- mdat.all[,2]
   btv.all  <- mdat.all[,3:(2+nyr+1)]
   n.viable.b   <- length(mdat.all[,1])
   n.viable.pt <- length(unique(mdat.all[,1]))
   cat("Found altogether",n.viable.b," viable trajectories for", n.viable.pt," r-k pairs\n")
 }

  #-------------------------------------------------------------------
  # 3 - if few points were found then resample and shrink the log k space
  #-------------------------------------------------------------------
  if (n.viable.b <= 1000){
    log.start_k.new  <- log(start_k) 
    max_attempts = 3
    current_attempts = 1
    while (n.viable.b <= 1000 && current_attempts <= max_attempts){
      if(n.viable.pt > 0) {
      log.start_k.new[1] <- mean(c(log.start_k.new[1], min(log(kv.all))))
      log.start_k.new[2] <- mean(c(log.start_k.new[2], max(log(kv.all)))) }
      n.new=n*current_attempts #add more points
      ri1 = exp(runif(n.new, log(start_r[1]), log(start_r[2])))  
      ki1 = exp(runif(n.new, log.start_k.new[1], log.start_k.new[2]))
      cat("Shrinking k space: repeating Monte Carlo in the interval [",exp(log.start_k.new[1]),",",exp(log.start_k.new[2]),"]\n")
      cat("Attempt ",current_attempts," of ",max_attempts," with ",n.new," additional points","\n")
      MCA <-  Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T)
      mdat.all <- rbind(mdat.all,MCA[[1]])
      rv.all   <- mdat.all[,1]
      kv.all   <- mdat.all[,2]
      btv.all  <- mdat.all[,3:(2+nyr+1)]
      n.viable.b   <- length(mdat.all[,1])
      n.viable.pt <- length(unique(mdat.all[,1]))
      cat("Found altogether",n.viable.b," viable trajectories for", n.viable.pt," r-k pairs\n")
      current_attempts=current_attempts+1 #increment the number of attempts
    }
  }

  #------------------------------------------------------------------
  # 4 - if tip of viable r-k pairs is 'thin', do extra sampling there
  #------------------------------------------------------------------
  gm.rv <- exp(mean(log(rv.all)))
  if(length(rv.all[rv.all > 0.9*start_r[2]]) < 5) { 
  l.sample.r        <- (gm.rv + max(rv.all))/2
  cat("Final sampling in the tip area above r =",l.sample.r,"with",3*n,"additional points\n")
  log.start_k.new <- c(log(0.8*min(kv.all)),log(max(kv.all[rv.all > l.sample.r])))
  if(is.na(gm.rv) | length(rv.all) < 3){  # if no viable r-k were found, try with 100000 points 
    n <- 33333
    l.sample.r <- start_r[1]
    log.start_k.new <- log(start_k)
    cat("No viable r-k found, resampling with 100,000 new points \n")
  }
  ri1 = exp(runif(3*n, log(l.sample.r), log(start_r[2])))  
  ki1 = exp(runif(3*n, log.start_k.new[1], log.start_k.new[2]))
  MCA <-  Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T)
  mdat.all <- rbind(mdat.all,MCA[[1]])
  rv.all   <- mdat.all[,1]
  kv.all   <- mdat.all[,2]
  btv.all  <- mdat.all[,3:(2+nyr+1)]
  n.viable.b   <- length(mdat.all[,1])
  n.viable.pt <- length(unique(mdat.all[,1]))
  cat("Found altogether",n.viable.b," viable trajectories for", n.viable.pt," r-k pairs\n")
  }


# accomodate shorter time series of CPUE
if(Btype=="CPUE") {
cpue     <- bt[is.na(bt)==F]
cpue.yr  <- yr[is.na(bt)==F]
cpue.nyr <- length(cpue)
}

 # ------------------------------------------------------------
 # Bayesian analysis of catch & biomass (or CPUE) with Schaefer model
 # ------------------------------------------------------------
FullSchaefer <- F
if(Btype != "None" & length(bt[is.na(bt)==F])>=10) {
FullSchaefer <- T
# accomodate shorter time series of biomass
bt.dat     <- bt[is.na(bt)==F]
bt.yr      <- yr[is.na(bt)==F]
bt.nyr     <- length(bt[is.na(bt)==F])
bt.ct      <- ct[is.na(bt)==F]

cat("Running Schaefer MCMC analysis....\n")
mcmc.burn <- as.integer(30000)
mcmc.chainLength <- as.integer(60000)  # burn-in plus post-burn
mcmc.thin = 10 # to reduce autocorrelation
mcmc.chains = 3 # needs to be at least 2 for DIC

if(Btype == "observed" | Btype=="simulated") {
    # Data to be passed on to JAGS
   jags.data        <- c('bt.ct','bt.dat','bt.nyr', 'start_r','startbio','start_k')
   # Parameters to be returned by JAGS
   jags.save.params <- c('r','k','sigma', 'alpha', 'sigma.log.r') # 
   
   # JAGS model
   Model = "model{
   # to avoid crash due to 0 values
   eps<- 0.01
    
     # define process (tau) and observation (sigma) variances as inversegamma prios

    itau2 ~ dgamma(4,0.01)
    tau2 <- 1/itau2
    tau <- pow(tau2,0.5)
    
    Pmean[1] <- log(alpha)
    P[1] ~ dlnorm(Pmean[1],itau2)
   
    for (t in 2:bt.nyr)
    {
    Pmean[t] <- log(max(P[t-1] + r*P[t-1]*(1-P[t-1]) - bt.ct[t-1]/k,eps))  # Process equation
    P[t] ~ dlnorm(Pmean[t],itau2) # Process error
    }
 
    # Now add observation error
    for (t in 1:bt.nyr){
    Bm[t] <- log(P[t]*k); #ignoring q=1
    bt.dat[t] ~ dlnorm(Bm[t],isigma2);
    }
  
   # priors
  # search in the alpha space from the center of the range. Allow high variability
  log.alpha               <- log((startbio[1]+startbio[2])/2)
  sd.log.alpha            <- (log.alpha-log(startbio[1]))/4
  tau.log.alpha           <- pow(sd.log.alpha,-2)
  alpha                   ~  dlnorm(log.alpha,tau.log.alpha)
  
  #inverse cubic root relationship between the range of viable r and the height of the density function 
  inverseRangeFactor <- 1/((start_r[2]-start_r[1])^1/3)
  #give sigma some variability in the inverse relationship 
  sigma.log.r         ~  dunif(0.001*inverseRangeFactor,0.02*inverseRangeFactor)
  tau.log.r           <- pow(sigma.log.r,-2)
  log.rm              <- log((start_r[1]+start_r[2])/2)  
  r                   ~  dlnorm(log.rm,tau.log.r)

  # search in the k space from the center of the range. Allow high variability
  log.km              <- log((start_k[1]+start_k[2])/2)
  sd.log.k            <- (log.km-log(start_k[1]))/4
  tau.log.k           <- pow(sd.log.k,-2)
  k                   ~  dlnorm(log.km,tau.log.k)
  
  isigma2 ~ dgamma(4,0.01)
  sigma2 <- 1/isigma2
  sigma <- pow(sigma2,0.5) 

  }" # end of model for Btype = observed or simulated


# ---------------------------------------------------------------------
# Schaefer model for Catch & CPUE
# ---------------------------------------------------------------------
} else {  

# get prior for q from stable catch/biomass period
if(is.na(qyr[1])==F) {
  mean.last.ct      <-mean(ct[yr >= qyr[1] & yr <= qyr[2]],na.rm=T) # get mean catch of last years
  mean.last.cpue    <-mean(bt[yr >= qyr[1] & yr <= qyr[2]],na.rm=T) # get mean of CPUE of last years
} else { 
# get prior range for q from mean catch and mean CPUE in recent years 
  lyr               <- ifelse(mean(start_r)>=0.5,5,10)  # determine number of last years to use, 5 for normal and 10 for slow growing fish    
  mean.last.ct      <-mean(ct[(length(ct)-lyr):length(ct)],na.rm=T) # get mean catch of last years
  mean.last.cpue    <-mean(bt[(length(ct)-lyr):length(ct)],na.rm=T) # get mean of CPUE of last years

}
if(mean(endbio) >= 0.5) {  # if biomass is high  
q.1           <- mean.last.cpue*0.25*gm.start_r/mean.last.ct
q.2           <- mean.last.cpue*0.5*start_r[2]/mean.last.ct
} else {
  q.1         <- mean.last.cpue*0.5*gm.start_r/mean.last.ct
  q.2         <- mean.last.cpue*start_r[2]/mean.last.ct
} 
q.prior       <- c(q.1,q.2)



  # Data to be passed on to JAGS
  jags.data        <- c('ct','cpue','cpue.nyr', 'start_r', 'start_k', 'startbio', 'q.prior')
  # Parameters to be returned by JAGS
  jags.save.params <- c('r','k','sigma', 'alpha','q') # 
    
# JAGS model
  Model = "model{
   # to avoid crash due to 0 values
   eps<-0.01
  
    Pmean[1] <- log(alpha)
    P[1] ~ dlnorm(Pmean[1],itau2)
    for (t in 2:cpue.nyr)
    {
    Pmean[t] <- log(max(P[t-1] + r*P[t-1]*(1-P[t-1]) - ct[t-1]/k,eps))  # Process equation
    P[t] ~ dlnorm(Pmean[t],itau2) # Process error
    }
    
    for (t in 1:cpue.nyr){
    cpuem[t] <- log(q*P[t]*k);
    cpue[t] ~ dlnorm(cpuem[t],isigma2);
    }
 
   # priors
#   alpha  ~  dunif(startbio[1],startbio[2]) # needed for fit of first biomass
   log.alpha               <- log((startbio[1]+startbio[2])/2)
   sd.log.alpha            <- (log.alpha-log(startbio[1]))/4
   tau.log.alpha           <- pow(sd.log.alpha,-2)
   alpha                   ~  dlnorm(log.alpha,tau.log.alpha)

  #inverse cubic root relationship between the range of viable r and the height of the density function
  inverseRangeFactor <- 1/((start_r[2]-start_r[1])^1/3)
  # give sigma some variability in the inverse relationship
  sigma.log.r         ~  dunif(0.001*inverseRangeFactor,0.02*inverseRangeFactor)
  tau.log.r           <- pow(sigma.log.r,-2)
  log.rm              <- log((start_r[1]+start_r[2])/2) 
  r                   ~  dlnorm(log.rm,tau.log.r)

  # search in the k space from the center of the range. Allow high variability
  log.km              <- log((start_k[1]+start_k[2])/2)
  sd.log.k            <- (log.km-log(start_k[1]))/4
  tau.log.k           <- pow(sd.log.k,-2)
  k                   ~  dlnorm(log.km,tau.log.k)
 
  # set realistic prior for q
  log.qm              <- mean(log(q.prior))
  sd.log.q            <- (log.qm-log(q.prior[1]))/4
  tau.log.q           <- pow(sd.log.q,-2)
  q                   ~  dlnorm(log.qm,tau.log.q)

  # define process (tau) and observation (sigma) variances as inversegamma prios
  
  itau2 ~ dgamma(4,0.01)
  tau2 <- 1/itau2
  tau <- pow(tau2,0.5)
  
  isigma2 ~ dgamma(4,0.01)
  sigma2 <- 1/isigma2
  sigma <- pow(sigma2,0.5)
  
}" # end model

} 


   # Write JAGS model to file
   cat(Model, file="r2jags.bug")  
   
   ### random seed
   set.seed(runif(1,1,500)) # needed in JAGS 
   
   ### run model
   jags_outputs <- jags(data=jags.data, 
                        working.directory=NULL, inits=NULL, 
                        parameters.to.save= jags.save.params, 
                        model.file="r2jags.bug", n.chains = mcmc.chains, 
                        n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength,
                        refresh=mcmc.burn/20, )
   
   # ------------------------------------------------------
   # Results from JAGS Schaefer
   # ------------------------------------------------------
   r_out            <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$r))
   k_out            <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$k))
   alpha_out        <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$alpha))
   
   median.r.jags        <- median(r_out) 
   lcl.r.jags       <- as.numeric(quantile(r_out,0.025))
   ucl.r.jags       <- as.numeric(quantile(r_out,0.975))
   median.k.jags        <- median(k_out)
   lcl.k.jags       <- as.numeric(quantile(k_out,0.025))
   ucl.k.jags       <- as.numeric(quantile(k_out,0.975))
   MSY.posterior 	   <- r_out*k_out/4 # simpler
   median.MSY.jags      <- median(MSY.posterior)
   lcl.MSY.jags     <- as.numeric(quantile(MSY.posterior,0.025))
   ucl.MSY.jags     <- as.numeric(quantile(MSY.posterior,0.975))

  if(Btype=="CPUE") {
    q_out           <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$q))
    median.q            <- median(q_out)
    lcl.q           <- as.numeric(quantile(q_out,0.025))
    ucl.q           <- as.numeric(quantile(q_out,0.975))


  }
} # end of MCMC Schaefer loop 


#------------------------------------
# get results from CMSY
#------------------------------------
# get estimate of most probable r as median of mid log.r-classes above cut-off
# get unique combinations of r-k
unique.rk         <- unique(mdat.all[,1:2])
# get remaining viable log.r and log.k 
rem            <- which(unique.rk[,1] > gm.rv)
rem.log.r      <- log(unique.rk[,1][rem])
rem.log.k      <- log(unique.rk[,2][rem])
# get vectors with numbers of r and mid values in about 25 classes
hist.log.r      <- hist(x=rem.log.r, breaks=25, plot=F)
log.r.counts    <- hist.log.r$counts
log.r.mids      <- hist.log.r$mids
# get most probable log.r as mean of mids with counts > 0
log.r.est       <- median(log.r.mids[which(log.r.counts > 0)])
lq.log.r        <- as.numeric(quantile(x=log.r.mids[which(log.r.counts > 0)], 0.25))
uq.log.r        <- as.numeric(quantile(x=log.r.mids[which(log.r.counts > 0)], 0.75))
lcl.log.r       <- as.numeric(quantile(x=log.r.mids[which(log.r.counts > 0)], 0.025))
ucl.log.r       <- as.numeric(quantile(x=log.r.mids[which(log.r.counts > 0)], 0.975))
r.est           <- exp(log.r.est)
lcl.r.est       <- exp(lcl.log.r)
ucl.r.est       <- exp(ucl.log.r)

# do linear regression of log k ~ log r with slope fixed to -1 (from Schaefer)
reg         <- lm(rem.log.k ~ 1 + offset(-1*rem.log.r))
int.reg     <- as.numeric(reg[1])
sd.reg      <- sd(resid(reg))
# get estimate of log(k) from y where x = log.r.est
log.k.est   <- int.reg + (-1) * log.r.est
# get estimates of CL of log.k.est from y +/- SD where x = lcl.log.r or ucl.log.r 
lcl.log.k   <- int.reg + (-1) * ucl.log.r - sd.reg
ucl.log.k   <- int.reg + (-1) * lcl.log.r + sd.reg
k.est       <- exp(log.k.est)
lcl.k.est   <- exp(lcl.log.k)
ucl.k.est   <- exp(ucl.log.k)

# get MSY from remaining log r-k pairs
log.MSY.est     <- mean(rem.log.r + rem.log.k - log(4))
sd.log.MSY.est  <- sd(rem.log.r + rem.log.k - log(4))
lcl.log.MSY.est <- log.MSY.est - 1.96*sd.log.MSY.est 
ucl.log.MSY.est <- log.MSY.est + 1.96*sd.log.MSY.est
MSY.est         <- exp(log.MSY.est)
lcl.MSY.est     <- exp(lcl.log.MSY.est)
ucl.MSY.est     <- exp(ucl.log.MSY.est)

# get predicted biomass vectors as median and quantiles 
# only use biomass trajectories from r-k pairs within the confidence limits
rem.btv.all <- mdat.all[which(mdat.all[,1] > lcl.r.est & mdat.all[,1] < ucl.r.est 
                & mdat.all[,2] > lcl.k.est & mdat.all[,2] < ucl.k.est),3:(2+nyr+1)]
median.btv <- apply(rem.btv.all,2, median)
lastyr.bt <- median.btv[length(median.btv)-1]
nextyr.bt <- median.btv[length(median.btv)]
lcl.btv    <- apply(rem.btv.all,2, quantile, probs=0.025)
q.btv      <- apply(rem.btv.all,2, quantile, probs=0.25)
ucl.btv    <- apply(rem.btv.all,2, quantile, probs=0.975)
lcl.lastyr.bt <- lcl.btv[length(lcl.btv)-1]
ucl.lastyr.bt <- ucl.btv[length(lcl.btv)-1]  
lcl.nextyr.bt <- lcl.btv[length(lcl.btv)]
ucl.nextyr.bt <- ucl.btv[length(lcl.btv)]

# -----------------------------------------
# Plot results 
# -----------------------------------------
# Analysis of viable r-k plot
# ----------------------------
plot(x=rv.all, y=kv.all, xlim=start_r, 
     ylim=c(0.9*min(kv.all, ifelse(Btype == "observed",k_out,NA), na.rm=T), max(kv.all)), 
     pch=16, col="gray",log="xy", bty="l",
     xlab="r", ylab="k", main="Analysis of viable r-k")
abline(v=gm.rv, lty="dashed")

# plot points and best estimate from full Schaefer analysis
if(FullSchaefer==T) {
  # plot r-k pairs from MCMC
  points(x=r_out, y=k_out, pch=16,cex=0.5)
  # plot best r-k pair from MCMC
  points(x=median.r.jags, y=median.k.jags, pch=19, col="green")  
  lines(x=c(lcl.r.jags, ucl.r.jags),y=c(median.k.jags,median.k.jags), col="green")
  lines(x=c(median.r.jags,median.r.jags),y=c(lcl.k.jags, ucl.k.jags), col="green")
}

# if data are from simulation, plot true r and k
if(Btype=="simulated" | substr(filename_1,1,3)=="Sim") {
  l.stock <- nchar(stock) # get length of sim stock name
  r.char  <- substr(stock,l.stock-1,l.stock)  # get last character of sim stock name
  r.sim   <- NA # initialize vector for r used in simulation
  if(r.char=="_H") {r.sim=1; lcl.r.sim=0.8; ucl.r.sim=1.25} else 
    if(r.char=="_M") {r.sim=0.5;lcl.r.sim=0.4;ucl.r.sim=0.62} else
      if(r.char=="_L") {r.sim=0.25;lcl.r.sim=0.2;ucl.r.sim=0.31} else {r.sim=0.05;lcl.r.sim=0.04;ucl.r.sim=0.062}
  # plot true r-k point with error bars
  points(x=r.sim, y=1000, pch=19, col="red")
  # add +/- 20% error bars  
  lines(x=c(lcl.r.sim,ucl.r.sim), y=c(1000,1000), col="red")
  lines(x=c(r.sim,r.sim), y=c(800,1250), col="red")
}

# plot blue dot for proposed r-k, with 95% CL lines 
points(x=r.est, y=k.est, pch=19, col="blue")
lines(x=c(lcl.r.est, ucl.r.est),y=c(k.est,k.est), col="blue")
lines(x=c(r.est,r.est),y=c(lcl.k.est, ucl.k.est), col="blue")

# Pred. biomass plot
#--------------------
# determine k to use for red line in b/k plot 
if(substr(filename_1,1,3)=="Sim") {k2use <- 1000} else 
      if(FullSchaefer==T)  {k2use <- median.k.jags} else {k2use <- k.est}
# determine hight of y-axis in plot
max.y  <- max(c(ifelse(Btype!="None",max(bt/k2use,na.rm=T),NA), 
              max(ucl.btv),0.6,startbio[2], endbio[2]),
              ifelse(FullSchaefer==T & Btype=="observed",max(ucl.k.jags/median.k.jags*bt/k2use,na.rm=T),NA),
              ifelse(FullSchaefer==T & Btype=="CPUE",1.1*max(cpue/(median.q*lcl.k.jags),na.rm=T),NA), na.rm=T)

# Main plot of relative CMSY biomass
plot(x=yr,y=median.btv[1:nyr], lwd=1.5, xlab="Year", ylab="Relative biomass b/k", type="l",
     ylim=c(0,max.y), bty="l", main=paste("Pred. biomass vs ", Btype,sep=""))
lines(x=yr, y=lcl.btv[1:nyr],type="l",lty="dotted")
lines(x=yr, y=ucl.btv[1:nyr],type="l",lty="dotted")
# plot lines for 0.5 and 0.25 biomass
abline(h=0.5, lty="dashed")
abline(h=0.25, lty="dotted")
# plot biomass windows
lines(x=c(yr[1],yr[1]), y=startbio, col="blue")
lines(x=c(intyr,intyr), y=intbio, col="blue")  
lines(x=c(max(yr),max(yr)), y=endbio, col="blue")  

# plot 25th percentile in final year if < 0.5 k and data are not simulated
if(Btype !="simulated" & substr(filename_1,1,3)!="Sim" & q.btv[yr==EndYear]<0.5) {
  points(x=EndYear,y=q.btv[yr==EndYear], col="purple", cex=1.5, lwd=2)}

# if observed biomass is available, plot red biomass line (use non-smoothed bt)
if(Btype=="observed" & FullSchaefer==T) {
  lines(x=bt.yr, y=bt[is.na(bt)==F]/median.k.jags,type="l", col="red", lwd=1) 
  lines(x=bt.yr, y=bt[is.na(bt)==F]/ucl.k.jags,type="l",col="red",lty="dotted") 
  lines(x=bt.yr, y=bt[is.na(bt)==F]/lcl.k.jags,type="l",col="red",lty="dotted")
}

# if observed CPUE is available, plot red biomass line
if(Btype=="CPUE" & FullSchaefer==T & substr(filename_1,1,3)!="Sim") {
  lines(x=cpue.yr, y=cpue/(median.q*median.k.jags),type="l", col="red", lwd=1) 
  lines(x=cpue.yr, y=cpue/(median.q*ucl.k.jags),type="l",col="red",lty="dotted") 
  lines(x=cpue.yr, y=cpue/(median.q*lcl.k.jags),type="l",col="red",lty="dotted")
}

if(Btype=="simulated") {
  # plot true biomass in red
  lines(x=yr, y=bt/k2use,type="l", col="red", lwd=1) 
  # plot BSM estimate of biomass in green 
  lines(x=yr, y=bt/median.k.jags,type="l", col="green", lwd=1) 
  lines(x=yr, y=bt/lcl.k.jags,type="l",col="green",lty="dotted") 
  lines(x=yr, y=bt/ucl.k.jags,type="l",col="green",lty="dotted")
}

if(filename_1=="SimCatchCPUE.csv") {
  # plot true biomass in red
  lines(x=cpue.yr, y=bt_true/1000,type="l", col="red", lwd=1) 
  # plot BSM estimate of biomass in green 
  lines(x=cpue.yr, y=cpue/(median.q*median.k.jags),type="l", col="green", lwd=1) 
  lines(x=cpue.yr, y=cpue/(median.q*lcl.k.jags),type="l",col="green",lty="dotted") 
  lines(x=cpue.yr, y=cpue/(median.q*ucl.k.jags),type="l",col="green",lty="dotted")
}

# if biomass or CPUE data are available but fewer than 10 years, plot on second axis
if(Btype != "None" & FullSchaefer==F) {
  par(new=T) # prepares for new plot on top of previous
  plot(x=yr, y=bt, type="l", col="red", lwd=1, 
    ann=F,axes=F,ylim=c(0,1.2*max(bt, na.rm=T))) # forces this plot on top of previous one
  axis(4, col="red", col.axis="red")
}

# Parabola plot
#-------------------------
max.y <- max(c(max(ct/MSY.est),ifelse(Btype=="observed"|Btype=="simulated",max(ct/median.MSY.jags),NA),1.2),na.rm=T)
# plot parabola
x=seq(from=0,to=2,by=0.001)
y=4*x-(2*x)^2
plot(x=x, y=y, xlim=c(1,0), ylim=c(0,max.y), type="l", bty="l",xlab="Relative biomass b/k", 
     ylab="Catch / MSY", main="Equilibrium curve")

# plot catch against CMSY estimates of relative biomass
points(x=median.btv[1:nyr], y=ct/MSY.est, pch=16, col="grey")

# plot catch scaled by BSM MSY against observed biomass scaled by BSM k
if(Btype == "observed") {
  points(x=bt/median.k.jags, y=ct/median.MSY.jags, pch=16, cex=0.5, col="black")
}

# for observed or simulated CPUE, plot catch scaled by BSM MSY against observed biomass derived as q * CPUE scaled by BSM k
if(FullSchaefer==T & Btype=="CPUE") {
  points(x=cpue/(median.q*median.k.jags), y=ct[is.na(bt)==F]/median.MSY.jags, pch=16, cex=0.5, col="black")
}

# for simulated biomass, scale catch by true MSY and scale simulated biomass by true k, show BSM estimate
if(Btype=="simulated") {
  points(x=bt/1000, y=ct/(1000*r.sim/4), pch=16, cex=0.5, col="red")
  points(x=bt/median.k.jags, y=ct/median.MSY.jags, pch=16, cex=0.5, col="black")
}

# for simulated CPUE, scale catch by true MSY and use true biomass scaled by true k
if(filename_1=="SimCatchCPUE.csv") {
  points(x=bt_true/1000, y=ct/(1000*r.sim/4), pch=16, cex=0.5, col="red")
}

# Exploitation rate plot
# -------------------------
# get u derived from predicted CMSY biomass
u.CMSY <- ct.raw/(median.btv[1:nyr]*k.est)
u.msy.CMSY  <- r.est/2 # Fmsy from CMSY expressed as exploitation rate

# get u from observed biomass
if(FullSchaefer==T & Btype != "CPUE") {
  u.bt       <- ct.raw/bt
  u.msy.bt   <- median.r.jags/2
}

# get u from simulated biomass
if(Btype == "simulated") {
  u.bt.sim   <- ct.raw/bt
  u.msy.sim   <- r.sim/2
}

# get u from CPUE
if(FullSchaefer==T & Btype == "CPUE") {
  u.bt.cpue   <- median.q*ct.raw/bt
  u.msy.cpue   <- median.r.jags/2
}

# get u from simulated CPUE
if(filename_1 == "SimCatchCPUE.csv") {
  u.bt.cpue.sim   <- ct.raw/bt_true
  u.msy.cpue.sim   <- r.sim/2
}

# if CPUE data are available but fewer than 10 years, plot on second axis
if(Btype == "CPUE") {
  q=1/(max(median.btv[1:nyr][is.na(bt)==F],na.rm=T)*k.est/max(bt,na.rm=T))
  u.cpue      <- q*ct/bt
}

# determine upper bound of Y-axis
max.y <- max(c(1.5, 1.2*u.CMSY/u.msy.CMSY,ct[yr==EndYear]/(q.btv[yr==EndYear]*k.est)/u.msy.CMSY,
               ifelse(Btype=="observed" & FullSchaefer==T,max(u.bt[is.na(u.bt)==F]/u.msy.bt),0),
               ifelse(Btype=="simulated",max(u.bt.sim[is.na(u.bt.sim)==F]/u.msy.sim),0),
               ifelse(FullSchaefer==T & Btype=="CPUE",max(u.bt.cpue/u.msy.cpue,na.rm=T),0),
               ifelse(filename_1 == "SimCatchCPUE.csv",max(u.bt.cpue.sim/u.msy.cpue.sim,na.rm=T),0),
               na.rm=T))
# plot u from CMSY
plot(x=yr,y=u.CMSY/u.msy.CMSY, type="l", bty="l", ylim=c(0,max.y), xlab="Year", 
     ylab="u / u_msy", main="Exploitation rate")
abline(h=1, lty="dashed")

# plot u from observed biomass
if(Btype == "observed" & FullSchaefer==T) lines(x=yr, y=u.bt/u.msy.bt, col="red")

# plot u from observed CPUE
if(FullSchaefer==T & Btype == "CPUE" & substr(filename_1,1,3)!="Sim") lines(x=yr, y=u.bt.cpue/u.msy.cpue, col="red")

# plot BSM u and true u from simulated biomass
if(Btype=="simulated")  {
  lines(x=yr, y=u.bt.sim/u.msy.sim, col="red")
  lines(x=yr, y=u.bt/u.msy.bt, col="green")
}
# plot BSM u from simulated CPUE
if(filename_1 == "SimCatchCPUE.csv") {
  lines(x=yr, y=u.bt.cpue.sim/u.msy.cpue.sim, col="red")
  lines(x=yr, y=u.bt.cpue/u.msy.cpue, col="green")
}
# plot u from CPUE on second axis if less than 10 years
if(FullSchaefer==F & Btype == "CPUE") {
  par(new=T) # prepares for new plot on top of previous
  plot(x=yr, y=u.cpue, type="l", col="red", ylim=c(0, 1.2*max(u.cpue,na.rm=T)),ann=F,axes=F)
  axis(4, col="red", col.axis="red")
}

#---------------------------------------------
# Plot Management-Graphs if desired
#---------------------------------------------
if(mgraphs==T) {
  # open window for plot of four panels
  windows(14,12) 
  par(mfrow=c(2,2))  
  # make margins narrower
  par(mar=c(3.1,4.1,2.1,2.1))
  
  # plot catch with MSY
  #---------------------
  max.y <- max(c(1.1*max(ct.raw),1.1*MSY.est))
  plot(x=yr, y=ct.raw, ylim=c(0,max.y), type ="l", bty="l", main=paste("Catch",stock), 
       ylab="Catch in 1000 t", lwd=2)
  abline(h=MSY.est,lty="dotted", lwd=1.5)
  abline(h=lcl.MSY.est, lty="dotted")

  # plot biomass with Bmsy
  #-----------------------
  # determine k to use for red line in b/k plot 
  k2use <- k.est
  # determine hight of y-axis in plot
  max.y.b  <- 2*max(c( max(ucl.btv),0.6,startbio[2], endbio[2]),
              ifelse(FullSchaefer==T,max(cpue/(median.q*median.k.jags),na.rm=T),1),na.rm=T)
  
  # Main plot of relative CMSY biomass
  plot(x=yr,y=2*median.btv[1:nyr], lwd=1.5, ylab="B / Bmsy", type="l",
       ylim=c(0,max.y.b), bty="l", main="Relative biomass")
  lines(x=yr, y=2*lcl.btv[1:nyr],type="l",lty="dotted")
  lines(x=yr, y=2*ucl.btv[1:nyr],type="l",lty="dotted")
  # plot lines for Bmsy and 0.5 Bmsy
  abline(h=1, lty="dotted", lwd=1.5)
  abline(h=0.5, lty="dotted")
  
  # plot 25th percentile in final year if B < Bmsy 
  if(q.btv[yr==EndYear]<0.5) {
    points(x=EndYear,y=2*q.btv[yr==EndYear], cex=1.5, lwd=2)}
  
  # if observed biomass is available, plot red biomass line (use non-smoothed bt)
  if(Btype=="observed" & FullSchaefer==T) {
    lines(x=bt.yr, y=2*bt[is.na(bt)==F]/median.k.jags,type="l", lty="dashed", lwd=1.5, col="red") 
  }
  
  # if observed CPUE is available, plot dashed biomass curve
  if(Btype=="CPUE" & FullSchaefer==T) {
    lines(x=cpue.yr, y=2*cpue/(median.q*median.k.jags),type="l", lty="dashed", lwd=1.5, col="red") 
  }
  
  # if biomass or CPUE data are available but fewer than 10 years, plot on second axis
  if(Btype != "None" & FullSchaefer==F) {
    par(new=T) # prepares for new plot on top of previous
    plot(x=yr, y=bt, type="l", lty="dashed", lwd=1.5, col="red",
         ann=F,axes=F,ylim=c(0,1.2*max(bt, na.rm=T))) # forces this plot on top of previous one
    axis(4)
  }
  
  # plot exploitation rate
  # -------------------------
  # get u derived from predicted CMSY biomass
  u.CMSY <- ct.raw/(median.btv[1:nyr]*k.est)
  u.msy.CMSY  <- r.est/2 # Fmsy from CMSY expressed as exploitation rate
  
  # get u from observed biomass
  if(FullSchaefer==T & Btype != "CPUE") {
    u.bt       <- ct.raw/bt
    u.msy.bt   <- median.r.jags/2
  }
  
  # get u from CPUE
  if(FullSchaefer==T & Btype == "CPUE") {
    u.bt.cpue   <- median.q*ct.raw/bt
    u.msy.cpue  <- median.r.jags/2
  }
  
  # if CPUE data are available but fewer than 10 years, plot on second axis
  if(Btype == "CPUE") {
    q=1/(max(median.btv[1:nyr][is.na(bt)==F],na.rm=T)*k.est/max(bt,na.rm=T))
    u.cpue      <- q*ct/bt
  }
  
  # determine upper bound of Y-axis
  max.y.u <- max(c(1.5, max(1.2*u.CMSY/u.msy.CMSY,na.rm=T),ct[yr==EndYear]/(q.btv[yr==EndYear]*k.est)/u.msy.CMSY,
                 ifelse(Btype=="observed" & FullSchaefer==T,max(u.bt[is.na(u.bt)==F]/u.msy.bt),0),
                 ifelse(FullSchaefer==T & Btype=="CPUE",max(u.bt.cpue/u.msy.cpue,na.rm=T),0),na.rm=T))
  # plot u from CMSY
  plot(x=yr,y=u.CMSY/u.msy.CMSY, type="l", bty="l", ylim=c(0,max.y.u), xlab="Year", 
       ylab="u / u.msy", main="Exploitation rate", lwd=1.5)
  abline(h=1, lty="dotted", lwd=1.5)
  
  # plot u from observed biomass
  if(Btype == "observed" & FullSchaefer==T) lines(x=yr, y=u.bt/u.msy.bt, lty="dashed", lwd=1.5, col="red")
  
  # plot u from observed CPUE
  if(FullSchaefer==T & Btype == "CPUE") lines(x=yr, y=u.bt.cpue/u.msy.cpue, lty="dashed", lwd=1.5, col="red")

  # plot u from CPUE on second axis if less than 10 years
  if(FullSchaefer==F & Btype == "CPUE") {
    par(new=T) # prepares for new plot on top of previous
    plot(x=yr, y=u.cpue, type="l",lty="dashed", lwd=1.5, ylim=c(0, 1.2*max(u.cpue,na.rm=T)),ann=F,axes=F, col="red")
    axis(4)
  }

  # plot stock-status graph
  #---------------------------
  max.x <- max(2, max(1.1*u.CMSY/u.msy.CMSY,na.rm=T))
  plot(x=u.CMSY/u.msy.CMSY, y=2*median.btv[1:nyr], xlim=c(0,max.x),ylim=c(0,2),ylab="",
       xlab="", main="Stock status", type="l", lwd=1.5, bty="l")
  mtext("u / u.msy",side=1, line=2)
  mtext("B / Bmsy",side=2, line=2)
  abline(h=1,lty="dotted",lwd=1.5)
  abline(h=0.5,lty="dotted")
  abline(v=1,lty="dotted",lwd=1.5)
  points(x=u.CMSY[yr==EndYear]/u.msy.CMSY,y=2*median.btv[yr==EndYear],cex=1.5,lwd=2, pch=5)
  points(x=0.75*max.x,y=2.0,cex=1.5,lwd=2, pch=5)
  text(EndYear,x=0.9*max.x, y=2.0)
  
}


# ------------------------------------------
# print input and results to screen
cat("---------------------------------------\n")

cat("Species:", cinfo$ScientificName[cinfo$stock==stock], ", stock:",stock,"\n")
cat("Name and region:", cinfo$EnglishName[cinfo$stock==stock], ",", cinfo$Region[cinfo$stock==stock], "\n")
cat("Catch data used from years", min(yr),"-", max(yr),", biomass =", Btype, "\n")
cat("Prior initial relative biomass =", startbio[1], "-", startbio[2], "\n")
cat("Prior intermediate rel. biomass=", intbio[1], "-", intbio[2], "in year", intyr, "\n")
cat("Prior final relative biomass   =", endbio[1], "-", endbio[2], "\n")
cat("If current catches continue, is the stock likely to crash within 3 years?",FutureCrash,"\n")
cat("Prior range for r =", format(start_r[1],digits=2), "-", format(start_r[2],digits=2),  
    ", prior range for k =", start_k[1], "-", start_k[2],"\n")
# if Schaefer on CPUE, print prior range of q
if(FullSchaefer==T & Btype=="CPUE") {
  cat("Prior range of q =",q.prior[1],"-",q.prior[2],"\n")
}

# if data are simulated, print true r-k
if(substr(filename_1,1,3)=="Sim") {
cat("True r =", r.sim, ", true k = 1000 (true values known because data were simulated)\n")
cat("True MSY =", 1000*r.sim/4,", true mean catch / MSY ratio =", mean(ct)/(1000*r.sim/4),"\n")
if(Btype=="simulated") cat("True biomass in last year =",bt[length(bt)],"or",bt[length(bt)]/1000,"k \n")
if(filename_1=="SimCatchCPUE.csv") cat("True biomass in last year =",bt_true[length(bt)],"or",bt_true[length(bt)]/1000,"k \n")
}

# print results from full Schaefer if available
if(FullSchaefer==T) {
cat("Results from Bayesian Schaefer model using catch &",Btype,"biomass\n")
cat("r =", median.r.jags,", 95% CL =", lcl.r.jags, "-", ucl.r.jags,", k =", median.k.jags,", 95% CL =", lcl.k.jags, "-", ucl.k.jags,"\n")
cat("MSY =", median.MSY.jags,", 95% CL =", lcl.MSY.jags, "-", ucl.MSY.jags,"\n")
if(Btype=="observed") {cat("Biomass in last year =",bt[length(bt)],"or",bt[length(bt)]/median.k.jags,"k \n")}

if(Btype == "CPUE") {
  cat("q =", median.q,", lcl =", lcl.q, ", ucl =", ucl.q,"\n")
  cat("Biomass in last year from CPUE/q =",cpue[length(cpue)]/median.q,"or",cpue[length(cpue)]/(median.k.jags*median.q),"k \n")
  }

}
# indicate if less than 10 years of biomass or CPUE are available
if(Btype !="None" & length(bt[is.na(bt)==F])<10) {
  cat("Less than 10 years with abundance data available, shown on second axis\n")
}

# results of CMSY analysis
cat("Results of CMSY analysis \n")
cat("Altogether", n.viable.b, "viable trajectories for", n.viable.pt," r-k pairs were found \n")
cat(length(rem.log.r),"r-k pairs above r =", gm.rv, "and",length(rem.btv.all[,1]),"trajectories within r-k CLs were analyzed\n")
cat("r =", r.est,", 95% CL =", lcl.r.est, "-", ucl.r.est,", k =", k.est,", 95% CL =", lcl.k.est, "-", ucl.k.est,"\n")
cat("MSY =", MSY.est,", 95% CL =", lcl.MSY.est, "-", ucl.MSY.est,"\n")
cat("Predicted biomass in last year =", lastyr.bt, ", 2.5th perc =", lcl.lastyr.bt, 
    "25th perc =",q.btv[length(median.btv)-1]," 97.5th perc =", ucl.lastyr.bt,"\n")
cat("Predicted biomass in next year =", nextyr.bt, ", 2.5th perc =", lcl.nextyr.bt, 
    "25th perc =",q.btv[length(median.btv)], ", 97.5th perc =", ucl.nextyr.bt,"\n")
cat("Comment:", comment,"\n")
cat("----------------------------------------------------------\n")

## Write some results into csv outfile
if(write.output == TRUE) {

# write data into csv file
  output = data.frame(cinfo$ScientificName[cinfo$stock==stock], stock, StartYear, EndYear,
           start_k[1],start_k[2],mean(ct.raw),max(ct.raw),ct[length(ct)],
    ifelse(FullSchaefer==T,median.MSY.jags,NA), # full Schaefer
    ifelse(FullSchaefer==T,lcl.MSY.jags,NA),
    ifelse(FullSchaefer==T,ucl.MSY.jags,NA),
    ifelse(FullSchaefer==T,median.r.jags,NA), 
    ifelse(FullSchaefer==T,lcl.r.jags,NA),
    ifelse(FullSchaefer==T,ucl.r.jags,NA),
    ifelse(FullSchaefer==T,median.k.jags,NA),                                          
    ifelse(FullSchaefer==T,lcl.k.jags,NA),                    
    ifelse(FullSchaefer==T,ucl.k.jags,NA),
    ifelse(FullSchaefer==T & Btype=="CPUE",q.1,NA),
    ifelse(FullSchaefer==T & Btype=="CPUE",q.2,NA),
    ifelse(FullSchaefer==T & Btype=="CPUE",median.q,NA),
    ifelse(FullSchaefer==T & Btype=="CPUE",lcl.q,NA),
    ifelse(FullSchaefer==T & Btype=="CPUE",ucl.q,NA),
    r.est, lcl.r.est, ucl.r.est, # CMSY r
    k.est, lcl.k.est, ucl.k.est, # CMSY k     
    MSY.est, lcl.MSY.est, ucl.MSY.est, # CMSY MSY
    ifelse(Btype!="None" & filename_1!="SimCatchCPUE.csv",bt[length(bt)],ifelse(filename_1=="SimCatchCPUE.csv",bt_true[length(bt)],NA)), # last biomass on record 
    lastyr.bt, lcl.lastyr.bt,q.btv[length(median.btv)-1],ucl.lastyr.bt, # CMSY last year bio: median, 2.5, 25, 97.5 perc
    nextyr.bt, lcl.nextyr.bt, q.btv[length(median.btv)], ucl.nextyr.bt) # CMSY last year + 1 bio

    write.table(output, file=outfile, append = T, sep = ",", 
              dec = ".", row.names = FALSE, col.names = FALSE)

  
# write screen text into text outfile.txt

cat("Species:", cinfo$ScientificName[cinfo$stock==stock], ", stock:",stock,"\n",
  "Name and region:", cinfo$Name[cinfo$stock==stock], ",", cinfo$Region[cinfo$stock==stock], "\n",
  "Catch data used from years", min(yr),"-", max(yr),", biomass =", Btype, "\n",
  "Prior initial relative biomass =", startbio[1], "-", startbio[2], "\n",
  "Prior intermediate rel. biomass=", intbio[1], "-", intbio[2], "in year", intyr, "\n",
  "Prior final relative biomass   =", endbio[1], "-", endbio[2], "\n",
  "If current catches continue, is the stock likely to crash within 3 years?",FutureCrash,"\n",
  "Prior range for r =", format(start_r[1],digits=2), "-", format(start_r[2],digits=2),  
  ", prior range for k  =", start_k[1], "-", start_k[2],"\n",
  file=outfile.txt,append=T)

  if(FullSchaefer==T & Btype=="CPUE") {
    cat("Prior range of q =",q.prior[1],"-",q.prior[2],"\n",file=outfile.txt,append=T)
  }

  if(substr(filename_1,1,3)=="Sim") {
  cat("True r =", r.sim, ", true k = 1000 (true values known because data were simulated)\n",
  "True MSY =", 1000*r.sim/4,", true mean catch / MSY ratio =", mean(ct)/(1000*r.sim/4),"\n",
  file=outfile.txt,append=T)
  }

  if(Btype=="simulated") { 
  cat("True biomass in last year =",bt[length(bt)],"or",bt[length(bt)]/1000,"k \n",
  file=outfile.txt,append=T) }

  if(filename_1=="SimCatchCPUE.csv") { 
    cat("True biomass in last year =",bt_true[length(bt)],"or",bt_true[length(bt)]/1000,"k \n",
    file=outfile.txt,append=T) }


if(FullSchaefer==T) {
  cat("Results from Bayesian Schaefer model using catch &",Btype,"biomass\n",
  "r =", median.r.jags,", 95% CL =", lcl.r.jags, "-", ucl.r.jags,
  ", k =", median.k.jags,", 95% CL =", lcl.k.jags, "-", ucl.k.jags,"\n",
  "MSY =", median.MSY.jags,", 95% CL =", lcl.MSY.jags, "-", ucl.MSY.jags,"\n",
  file=outfile.txt,append=T)
  
    if(Btype=="observed") {cat("Biomass in last year =",bt[length(bt)],
      "or",bt[length(bt)]/median.k.jags,"k \n", file=outfile.txt,append=T)}
  
    if(Btype == "CPUE") {
      cat("q =", median.q,", lcl =", lcl.q, ", ucl =", ucl.q,"\n",
          "Biomass in last year from q*CPUE =",cpue[length(cpue)]/median.q,
          "or",(cpue[length(cpue)])/(median.k.jags*median.q),"k \n",  
          file=outfile.txt,append=T)
    }
  }

  if(Btype !="None" & length(bt[is.na(bt)==F])<10) {
    cat("Less than 10 years with abundance data available, shown on second axis\n",file=outfile.txt,append=T)
  }

  cat(" Results of CMSY analysis with altogether",n.viable.b, "viable trajectories for", n.viable.pt,"r-k pairs \n",
  length(rem.log.r),"r-k pairs above r =", gm.rv, "and",length(rem.btv.all[,1]),"trajectories within r-k CLs were analyzed\n",
  "r =", r.est,", 95% CL =", lcl.r.est, "-", ucl.r.est,
  ", k =", k.est,", 95% CL =", lcl.k.est, "-", ucl.k.est,"\n",
  "MSY =", MSY.est,", 95% CL =", lcl.MSY.est, "-", ucl.MSY.est,"\n",
  "Predicted biomass last year=", lastyr.bt, ", 2.5th =", lcl.lastyr.bt, 
  ", 25th =",q.btv[length(median.btv)-1],", 97.5th =", ucl.lastyr.bt,"\n",
  "Predicted biomass next year=", nextyr.bt, ", 2.5th =", lcl.nextyr.bt, 
  ", 25th =",q.btv[length(median.btv)], ", 97.5th =", ucl.nextyr.bt,"\n",
  "----------------------------------------------------------\n",
  file=outfile.txt,append=T)  

  }

} # end of stocks loop