In laurieKell/biodyn: Biomass Dynamic Model

library(knitr)

opts_chunk$set(comment=NA, 
               fig.width =6, 
               fig.height=6,
               fig.path  ="../tex/",
               warning=FALSE, 
               message=FALSE, 
               error  =FALSE, 
               echo   =FALSE, 
               eval   =TRUE,
               cache  =TRUE)

#install.packages("C:/Users/laurie.kell/Documents/MEGAsync/temp/aspic_2.0.1.tar.gz",repos=NULL,type="source",INSTALL_opts=c("--no-multiarch"))

library(FLCore)
library(ggplotFL)
library(scales)
library(plyr)
library(reshape)
library(FLBRP)
library(biodyn)
library(aspic)
library(plotrix)
library(corrplot)
library(diags)
library(ucminf)
library(numDeriv)
library(FLFishery)
library(FLasher)
library(testthat)

source('~/Desktop/flr/git/biodyn/R/biodyn-F.R')

theme_set(theme_bw(10))

dirMy ="c:/temp/albn"
dirDat=paste(dirMy,"/data",sep="")
dirInp="http://rscloud.iccat.int/kobe/Inputs/albn/2013/aspic/"
dirInp="/home/laurie/Desktop/rfmos/iccat/kobe/Inputs/albn/2013/aspic"

# aspic objects
asp=aspics(paste(dirInp,"/run",1:7,"/aspic.inp",sep=""))
control(asp[[4]])["k","val"]=par[4,"K"]
control(asp[[4]])["msy","val"]=par[4,"MSY"]
control(asp[[4]])["k","fit"]=0

asp=aspics(llply(asp,fit))

plot(asp)+theme_bw()

Figure 1 ASPIC runs

i=c(1,2,4,5)[2] #Troll, CT, JLL, All
# 1&5 OK, 2&4 bad

bd  =as(asp[[i]],"biodyn")
cpue=index(asp[[i]],FALSE)
params(bd)[c("r")]=1.0*params(bd)[c("r")]
control(bd)[c("r"),"val"]=params(bd)[c("r")]

bd1=fwd(bd,catch=catch(bd))
bd1=fit(bd1,cpue)
plot(biodyns("ASPIC"=bd,"biodyn"=bd1))+
  theme_bw()

control(bd1)[c("p"),"val"]=1
bd1=fit(bd1,cpue)

plot(biodyns("ASPIC"=bd,"biodyn"=bd1))+
  theme_bw()

Figure 2 ASPIC runs

lag<-function(x) {
  res=FLQuant(NA,dimnames=dimnames(x))
  res[,-dim(res)[2]][]=c(x[,-1])

  res}
i=2
tm=model.frame(mcf(FLQuants(
         lag    =lag(stock(asp[[i]])),
         biomass=stock(asp[[i]]),
         catch  =catch(asp[[i]]))),drop=T)
tm=transform(tm,prd=lag-biomass+catch)
ggplot(tm)+
  geom_path(aes(biomass,prd,col=year))+
  geom_point(aes(biomass,prd,col=year))+
  theme_bw()

Figure 3 Productivity for run 2

Newton Rhapson

Estimate stock biomass next year.

#PELLA , J. J., 1967 A study of methods to estimate the Schaefer model parameters with special reference to the yellowfin tuna fishery in the eastern tropical Pacific ocean. University of Washington, Seattle.

# PELLA , J. J., and P. K. T OMLINSON , 1969 A generalized stock production model. Bulletin of the Inter-American Tropical Tuna Commission 13: 419-496.

# PRAGER , M. H., 1994 A suite of extensions to a nonequilibrium surplus-production model. U. S. Fishery Bulletin 92: 374-389.


C=catch(asp[["2"]])
H=C/stock(asp[["2"]])[,dimnames(C)$year]

dat2=dget("/home/laurie/Desktop/rfmos/iccat/kobe/Inputs/albn/2013/aspic/run2/aspic.rdat")

F=window(as.FLQuant(transform(dat2$t.series[c("year","F.total")],data=F.total)[,-2]),
         end=2011)
B=window(as.FLQuant(transform(dat2$t.series[c("year","b")],      data=b      )[,-2]),
         end=2011)
F.=F
B.=B

r=dat2$estimates["r"]
K=dat2$estimates["K"]
b0=dat2$estimates["B1.K"]

mm2=nr(C,H,B,r,K,tolVal=1e-6,niter=20)
plot(mm2[["F"]])+
  geom_line(aes(year,data),col="red",data=as.data.frame(F))
plot(mm2[["B"]])+
  geom_line(aes(year,data),col="red",data=as.data.frame(B))
plot(yield(mm2[["F"]],mm[["B"]],r,K))+
  geom_line(aes(year,data),col="red",data=as.data.frame(C))

mm1=slv(C,H,B,r,K)
plot(mm1[["F"]])+
  geom_line(aes(year,data),col="red",data=as.data.frame(F))
plot(mm1[["B"]])+
  geom_line(aes(year,data),col="red",data=as.data.frame(B))
plot(yield(mm1[["F"]],mm1[["B"]],r,K))+
  geom_line(aes(year,data),col="red",data=as.data.frame(C))

hvt=nr(catch(bd1),harvest(bd1),stock(bd),params(bd1)["r"],params(bd)["k"],
tolVal=1e-6,niter=20)

ggplot(plot(FLQuants(
             #"Differential"=fCPP(harvest(bd), catch(bd), 
              #                             bd@stock, as.FLQuant(params(bd)[c("r","k")])),
             "NR"=hvt[["F"]],
              #"ASPIC"=harvest(asp[[2]]),
              "Difference"  =harvest(bd1))
              )$data)+
  geom_line(aes(year,`50%`,col=qname))

ggplot(plot(FLQuants("Differential"=biodyn:::stockCPP(harvest(bd), catch(bd), 
                                           bd@stock, as.FLQuant(params(bd)[c("r","k")])),
              "ASPIC"=stock(asp[[2]]),
              "Difference"  =stock(bd)))$data)+
  geom_line(aes(year,`50%`,col=qname))