R/kobe-FLSAM.R

In flr/FLSAM: An Implementation of the State-Space Assessment Model for FLR

kobeFLSamFn=function(object,what=c("trks","smry","pts","ellipse")[1],prob=c(0.75,0.5,.25),nits=1000,bmsy=1,fmsy=1){
  
  require(ellipse)
  require(reshape2)
  require(plyr)
  
  trmnl=c(max(seq(dim(rescov(object))[1])[dimnames(rescov(object))[[1]]=="logssb"]),
          max(seq(dim(rescov(object))[1])[dimnames(rescov(object))[[1]]=="logfbar"]))
  
  currentCov=rescov(object)[trmnl,trmnl]
  
  ssb =params(object)[seq(dim(params(object))[1])[params(object)[,"name"]=="logssb"],2:3]
  fbar=params(object)[seq(dim(params(object))[1])[params(object)[,"name"]=="logfbar"],2:3]
  
  trks=cbind(melt(data.frame(year=range(object)["minyear"]:range(object)["maxyear"],ssb),id="year"),
             melt(fbar)[,2])
  names(trks)[2:4]=c("quantity","stock","harvest")
  trks[,c("stock","harvest")] <- exp(trks[,c("stock","harvest")])
  
  ellipse=as.data.frame(ellipse(currentCov,level=prob[1]))
  
  names(ellipse)=c("stock","harvest")
  ellipse=transform(ellipse, stock=stock+rev(ssb[,1])[1],harvest=harvest+rev(fbar[,1])[1])
  ellipse <- exp(ellipse)
  
  rnd=c(t(mvrnorm(nits,c(rev(ssb[,1])[1],rev(fbar[,1])[1]),currentCov)))
  rnd=as.data.frame(t(array(c(rnd),c(2,nits))))
  names(rnd)=c("stock","harvest")
  rnd <- exp(rnd)
  
  kb=object=llply(list(trks=trks,ellipse=ellipse,pts=rnd), function(x,bmsy,fmsy) transform(x, stock=stock/bmsy, harvest=harvest/fmsy), bmsy=bmsy, fmsy=fmsy)
  
  return(object)}

#### exported function
kobeFLSAM=function(object,what=c("trks","smry","pts","ellipse")[1],prob=c(0.75,0.5,.25),nits=1000,bmsy=1,fmsy=1){
  require(plyr)
  
  if (length(object)==1){
    res=kobeFLSamFn(object,what=what,prob=prob,nits=nits,bmsy=bmsy,fmsy=fmsy)
  } else {  
    res=mlply(object,   function(x,what=what,prob=prob,nits=nits,bmsy=bmsy,fmsy=fmsy)
      kobeFLSamFn(x,what=what,prob=prob,nits=nits,bmsy=bmsy,fmsy=fmsy),
              what=what,prob=prob,nits=nits,bmsy=bmsy,fmsy=fmsy)
    
    res=list(trks   =ldply(res, function(x) x$trks),
             pts    =ldply(res, function(x) x$pts),
             smry   =ldply(res, function(x) x$smry),
             ellipse=ldply(res, function(x) x$ellipse))
  }
  
  if (length(what)==1) return(res[[what]]) else return(res[what])}

#library(ggplot2)
#kb=kobeFLSamFn(sam)
#ggplot(kb$ellipse)+geom_path(aes(stock,harvest))+geom_path(aes(stock,harvest),data=subset(kb$trks,quantity=="value"))