set.seed(999)  ## for same random sequence
#require(hacks)

#setwd("C:/Users/Ye/Documents/Data poor fisheries/Martell Froese Method/")

## Read Data for stock, year=yr, catch=ct, and resilience=res. Expects space delimited file with header yr  ct and years in integer and catch in real with decimal point
## For example
## stock  res	 yr     ct       
## cap-icel	Medium 1984  1234.32 

## filename <- "RAM_MSY.csv"
##filename <- "ICESct2.csv"

cat("Step 1","\n")
TestRUN <- F  # if it is true, just run on the test samples, false will go for a formal run!

filename <- "D20.csv"
outfile  <- "CatchMSY_Output1.csv"
outfile2 <- paste("NonProcessedSpecies.csv",sep="")

#cdat <- read.csv2(filename, header=T, dec=".")
cdat1 <- read.csv(filename)
cat("\n", "File", filename, "read successfully","\n")

cat("Step 2","\n")
if(file.exists("cdat.RData"))
{load("cdat.RData")} else 
{
  
  dim(cdat1)
  yrs=1950:2012
  
  # to set NA as 0
  cdat1[is.na(cdat1)] <- 0
  nrow <- length(cdat1[,1])
  ndatColn <- length(cdat1[1,c(-1:-12)])
  rownames(cdat1) <- NULL
  
  cdat <- NULL
  for(i in 1:nrow)
  {#i=1
    #a <- ctotal3[i,-1]
    tmp=data.frame(stock=rep(as.character(cdat1[i,"Stock_ID"]),ndatColn),
                   species=rep(as.character(cdat1[i,"Scientific_name"]),ndatColn),
                   yr=yrs,ct=unlist(c(cdat1[i,c(-1:-12)])),
                   res=rep(cdat1[i,"ResilienceIndex"],ndatColn))
    
    cdat <- rbind(cdat,tmp)
    #edit(cdat)
  }
} 

StockList=unique(as.character(cdat$stock))

colnames(cdat)


#stock_id <- unique(as.character(cdat$stock)) 
#??
# stock_id <- "cod-2224" ## for selecting individual stocks
# stock=stock_id
#??

cat("Step 3","\n")

## FUNCTIONS are going to be used subsequently
.schaefer  <- function(theta)
{
  with(as.list(theta), {  ## for all combinations of ri & ki
    bt=vector()
    ell = 0  ## initialize ell
    J=0 #Ye
    for (j in startbt)
    {
      if(ell == 0) 
      {
        bt[1]=j*k*exp(rnorm(1,0, sigR))  ## set biomass in first year
        for(i in 1:nyr) ## for all years in the time series
        {
          xt=rnorm(1,0, sigR)
          bt[i+1]=(bt[i]+r*bt[i]*(1-bt[i]/k)-ct[i])*exp(xt) 
          ## calculate biomass as function of previous year's biomass plus net production minus catch
        }
        
        #Bernoulli likelihood, assign 0 or 1 to each combination of r and k
        ell = 0
        if(bt[nyr+1]/k>=lam1 && bt[nyr+1]/k <=lam2 && min(bt) > 0 && max(bt) <=k && bt[which(yr==interyr)]/k>=interbio[1] && bt[which(yr==interyr)]/k<=interbio[2]) 
          ell = 1
        J=j # Ye
      }	
    }
    return(list(ell=ell,J=J)) # Ye adding J=J
    
    
  })
}

sraMSY	<-function(theta, N)
{
  #This function conducts the stock reduction
  #analysis for N trials
  #args:
  #	theta - a list object containing:
  #		r (lower and upper bounds for r)
  #		k (lower and upper bounds for k)
  #		lambda (limits for current depletion)
  
  
  with(as.list(theta), 
{
  ri = exp(runif(N, log(r[1]), log(r[2])))  ## get N values between r[1] and r[2], assign to ri
  ki = exp(runif(N, log(k[1]), log(k[2])))  ## get N values between k[1] and k[2], assing to ki
  itheta=cbind(r=ri,k=ki, lam1=lambda[1],lam2=lambda[2], sigR=sigR) 
  ## assign ri, ki, and final biomass range to itheta
  M = apply(itheta,1,.schaefer) ## call Schaefer function with parameters in itheta
  i=1:N
  ## prototype objective function
  get.ell=function(i) M[[i]]$ell
  ell = sapply(i, get.ell) 
  get.J=function(i) M[[i]]$J # Ye
  J=sapply(i,get.J) # Ye
  return(list(r=ri,k=ki, ell=ell, J=J)) # Ye adding J=J	
})
}


getBiomass  <- function(r, k, j)
{
  BT <- NULL
  bt=vector()
  for (v in 1:length(r))
  {
    bt[1]=j[v]*k[v]*exp(rnorm(1,0, sigR))  ## set biomass in first year
    for(i in 1:nyr) ## for all years in the time series
    {
      xt=rnorm(1,0, sigR)
      bt[i+1]=(bt[i]+r[v]*bt[i]*(1-bt[i]/k[v])-ct[i])*exp(xt) 
      ## calculate biomass as function of previous year's biomass plus net production minus catch
    }
    BT=rbind(BT, t(t(bt)))        
  }
  return(BT)
}

## The End of Functions section

cat("Step 4","\n")
stockLoop <- StockList
# randomly select stocks from randomly selected 5 area codes first
if(TestRUN) 
{
  set.seed(999)
  AreaCodeList <- unique(cdat1$AREA_Code)
  sampledAC <- sample(AreaCodeList,size=5,replace=F)
  stockLoop <- cdat1[cdat1$AREA_Code %in% sampledAC,c("Stock_ID")]
}

#setup counters
counter1 <- 0
counter2 <- 0

cat("Step 4","\n")
## Loop through stocks
write.table("x",file=outfile,append = FALSE, row.names = FALSE,col.names=FALSE,sep=",")
write.table("x",file=outfile2,append = FALSE, row.names = FALSE,col.names=FALSE,sep=",")

for(stock in stockLoop) 
{
  t0<-Sys.time()
  xr <- runif(1, 1.0, 10000)
  x1<-c(paste("processed",xr,sep=","))
  xr <- runif(1, 1.0, 10000)
  x2<-c(paste("non processed",xr,sep=","))
  write.table(x1,file=outfile,append = T, row.names = FALSE,col.names=FALSE,sep=",")
  write.table(x2,file=outfile2,append = T, row.names = FALSE,col.names=FALSE,sep=",")
    
  cat("Elapsed: ",Sys.time()-t0," \n")
}