R/biodyn-FLBRP2biodyn.R

In laurieKell/mpb: Management Procedures

FLBRP2biodyn=function(from,quantity=c("ssb","biomass","exploitable")[2]){

  warn=options()$warn
  options(warn=-1)

  #r  =lambda(leslie(from,f=c(FLBRP:::refpts(from)["crash","harvest"]))[drop=T])-1
  msy=from@refpts["msy","yield"]

  fbar(from)[,1]=0
  k=switch(tolower(substr(quantity[1],1,1)),
      s=from@refpts["virgin","ssb"],
      b=from@refpts["virgin","biomass"],
      e=apply(stock.n(from)[,1]%*%stock.wt(from)[,1]%*%catch.sel(from)%/%
        apply(catch.sel(from),c(2,6),max),6,sum))

  b0=switch(tolower(substr(quantity[1],1,1)),
        s=from@ssb.obs[,1]%/%k,
        b=from@stock.obs[,1]%/%k,
        e=from@stock.obs[,1]%/%from@refpts["virgin","biomass"])

  from@fbar[,1]=from@refpts["msy","harvest"]
  bmsy=switch(tolower(substr(quantity[1],1,1)),
           s=from@refpts["msy","ssb"],
           b=from@refpts["msy","biomass"],
           e=apply(stock.n(from)[,1]%*%stock.wt(from)[,1]%*%catch.sel(from)%/%
             apply(catch.sel(from),c(2,6),max),6,sum))

  # if (tolower(substr(fix,1,1))=="m")
  #   p=optimise(function(p,msy,r,k) {
  #             res=(msy/(r*k)-(1/(1+p))^(1/p+1))^2
  #             res},
  #              c(1e-12,1e12),
  #              msy=msy,
  #              r  =r,
  #              k  =k)$minimum
  # else
  p=getP(bmsy,k,p=c(0.001,5))
  r=msy/(k*(1/(1+p))^(1/p+1))
  
  bd=biodyn(catch=catch.obs(from))
  #par=pellat(from)

  bd@params=FLPar(c(r=r,k=k,p=p,b0=b0))

  bd@control["r", c("min","val","max")][]=c(params(bd)["r"])
  bd@control["k", c("min","val","max")]=params(bd)["k"]
  bd@control["p", c("min","val","max")]=params(bd)["p"]
  bd@control["b0",c("min","val","max")]=params(bd)["b0"]

  bd@control[,"min"]=bd@control[,"min"]*.1
  bd@control[,"max"]=bd@control[,"max"]*10

  bd@priors["r",   "a"]=params(bd)["r"]
  bd@priors["k",   "a"]=params(bd)["k"]
  bd@priors["p",   "a"]=params(bd)["p"]
  bd@priors["b0",  "a"]=params(bd)["b0"]
  bd@priors[ "msy","a"]=mpb::refpts(bd)["msy"]
  bd@priors["bmsy","a"]=mpb::refpts(bd)["bmsy"]
  bd@priors["fmsy","a"]=mpb::refpts(bd)["fmsy"]

  #bd=mpb::fwd(bd,catch=catch(bd))

  range(bd)["minyear"]=dims(bd@catch)$minyear
  range(bd)["maxyear"]=dims(bd@catch)$maxyear

  options(warn=warn)

  bd}

setAs('FLBRP','biodyn',
      function(from) FLBRP2biodyn(from,to))