library(stockassessment)

oldwd<-setwd("data")

  cn<-read.ices("cn.dat")
  cw<-read.ices("cw.dat")
  dw<-read.ices("dw.dat")
  lw<-read.ices("lw.dat")
  mo<-read.ices("mo.dat")
  nm<-read.ices("nm.dat")
  pf<-read.ices("pf.dat")
  pm<-read.ices("pm.dat")
  sw<-read.ices("sw.dat")
  lf<-read.ices("lf.dat")
  surveys<-read.ices("survey.dat")
  surveys<-c(surveys, read.ices("Survey_DARDS.dat"))
  
 ### remove NaN #first check dw and lw in case there are NaN
  dw[which(!is.finite(dw))] <- 0
  lw[which(!is.finite(lw))] <- 0
  sw[which(!is.finite(sw))] <- 0
  cw[which(!is.finite(cw))] <- 0

setwd(oldwd)

dat<-setup.sam.data(surveys=surveys,
                    residual.fleet=cn, 
                    prop.mature=mo, 
                    stock.mean.weight=sw, 
                    catch.mean.weight=cw, 
                    dis.mean.weight=dw, 
                    land.mean.weight=lw,
                    prop.f=pf, 
                    prop.m=pm, 
                    natural.mortality=nm, 
                    land.frac=lf)
###################################################
### function for using different discard survival rates
set_survival <- function(survival_rate = 0, ### proportion of discards surviving
                         surveys, residual.fleet, prop.mature, 
                         stock.mean.weight, catch.mean.weight, dis.mean.weight,
                         land.mean.weight, prop.f, prop.m, natural.mortality,
                         land.frac, ...){
  
  ### separate landings and discards numbers
  landings.n <- residual.fleet * land.frac # catch.n *  landings.n/catch.n
  discards.n <- residual.fleet * (1 - land.frac)
  #all.equal((landings.n + discards.n), residual.fleet) ### worked
  
  ### apply survival rate to discards->only dead portion of discards is considered in the assessment
  discards.n <- discards.n * (1 - survival_rate)
  
  ### sum up to total catch numbers
  residual.fleet <- landings.n + discards.n
  
  ### update landings proportion, as discards only include now dead portion
  land.frac <- landings.n / (landings.n + discards.n)
  ### remove NaN
  land.frac[which(!is.finite(land.frac))] <- 1
  

  
  ### calculate catch weight as weighted mean.   
  catch.mean.weight <- land.mean.weight * (landings.n / residual.fleet) +
    dis.mean.weight * (discards.n / residual.fleet)
  catch.mean.weight[which(!is.finite(catch.mean.weight))] <- dis.mean.weight[which(!is.finite(catch.mean.weight))] 
  catch.mean.weight[which(!is.finite(catch.mean.weight))] <- land.mean.weight[which(!is.finite(catch.mean.weight))]
  
  ### set catch weight as stock weight
  stock.mean.weight <- catch.mean.weight
  
  ### set up SAM input object
  res <- setup.sam.data(surveys = surveys,
                        residual.fleet = residual.fleet,
                        prop.mature = prop.mature,
                        stock.mean.weight = stock.mean.weight,
                        catch.mean.weight = catch.mean.weight,
                        dis.mean.weight = dis.mean.weight,
                        land.mean.weight = land.mean.weight,
                        prop.f = prop.f,
                        prop.m = prop.m,
                        natural.mortality = natural.mortality,
                        land.frac = land.frac)
  
  ### return
  return(res)
  
}


dat_dead <- set_survival(survival_rate = 0.4,
                         surveys = surveys, residual.fleet = cn, prop.mature = mo,
                         stock.mean.weight = sw, catch.mean.weight = cw, 
                         dis.mean.weight = dw, land.mean.weight = lw, prop.f = pf,
                         prop.m = pm, natural.mortality = nm, land.frac = lf)
########################################################################################


#save(dat, file="run/data.RData")

save(dat_dead, file="run/dat_dead.RData")