library(FLCore)
library(FLBRP)
library(FLasher)
library(FLife)
library(ggplotFL)
library(FLCandy)
library(rfishbase)
library(SPMpriors)
library(FishLife)
library(plyr)
library(dplyr)
library(reshape)
library(ggplot2)
theme_set(theme_bw())
library(ggpubr)
library(ggcorrplot)
library(GGally)
h <-function(x) 1-exp(-x)
f <-function(x) -log(1-x)
M1Fn=function(x,y=0.025) FLQuant(y,dimnames=dimnames(m(x)))
M2Fn=function(x) m(x) - M1Fn(x)
pars=flmvn_traits(Genus="Limnothrissa",Species="miodon",Plot=FALSE)
[,1] [,2]
[1,] "K" "M"
[2,] "Winfinity" "Loo"
[3,] "tmax" "tm"
[4,] "Lm" "Temperature"
[5,] "ln_margsd" "rho"
[6,] "logitbound_h" "ln_r"
pars=FLPar(c(cast(pars[[1]][,c("trait","mu.sp")],.~trait)[,c("h","K","tm","Loo","sigR")]))
dimnames(pars)[[1]]=c("s","k","t50","linf","sigmaR")
L.midon=pars
pars=flmvn_traits(Genus="Stolothrissa",Species="tanganicae",Plot=FALSE)
[,1] [,2]
[1,] "K" "M"
[2,] "Winfinity" "Loo"
[3,] "tmax" "tm"
[4,] "Lm" "Temperature"
[5,] "ln_margsd" "rho"
[6,] "logitbound_h" "ln_r"
pars=FLPar(c(cast(pars[[1]][,c("trait","mu.sp")],.~trait)[,c("h","K","tm","Loo","sigR")]))
dimnames(pars)[[1]]=c("s","k","t50","linf","sigmaR")
S.tanganicae=pars
pars=flmvn_traits(Genus="Rastrineobola",Species="argentea",Plot=FALSE)
[,1] [,2]
[1,] "K" "M"
[2,] "Winfinity" "Loo"
[3,] "tmax" "tm"
[4,] "Lm" "Temperature"
[5,] "ln_margsd" "rho"
[6,] "logitbound_h" "ln_r"
pars=FLPar(c(cast(pars[[1]][,c("trait","mu.sp")],.~trait)[,c("h","K","tm","Loo","sigR")]))
dimnames(pars)[[1]]=c("s","k","t50","linf","sigmaR")
R.argentea=pars
save(R.argentea,S.tanganicae,L.midon,file="../data/pars.RData")
pars=lhPar(L.midon)
eq=lhEql(pars)
landings.sel(eq)["0"]=landings.sel(eq)["1"]*0.1
eq=brp(eq)
fbar(eq)=FLQuant(c(seq(0,0.5,length.out=100),seq(0.505,5.7,length.out=200)))
f=FLQuant(c(rep(0.1,60),seq(0.1,2.5,length.out=40)[-40],
seq(2.5,0.8,length.out=11),rep(0.8,51)))
f=f%*%refpts(eq)["msy","harvest"]
fbar(eq)=f
eq=brp(eq)
f=FLQuant(c(rep(0.1,60),seq(0.1,2.5,length.out=40)[-40],
seq(2.5,0.8,length.out=11),rep(0.8,51)))
fbar(eq)=f%*%refpts(eq)["msy","harvest"]
om=as(eq,"FLStock")
dimnames(om) <- list(year=an(dimnames(om)$year)+1950)
f=propagate(f,3)
dimnames(f)$year=dimnames(fbar(om))$year
om=fwd(om,sr=eq,fbar=f[,-1])
Figure 1 Vectors-at-age
fbar(eq)=FLQuant(seq(0,1,length.out=101))*max(refpts(eq)["crash","harvest"])
plot(eq,refpts="msy")
Figure 2 Equilibrium dynamics
plot(window(om,start=2000), metrics=list(Rec=rec,SSB=ssb,Catch=landings,F=fbar))+
theme_bw()+
geom_flpar(data=FLPars(Rec =FLPar("Rmsy"=refpts(eq)[c("msy"),"rec",1,drop=T]),
SSB =FLPar("Bmsy"=refpts(eq)[c("msy"),"ssb",1,drop=T]),
F =FLPar("Fmsy"=refpts(eq)[c("msy"),"harvest",1,drop=T]),
Catch=FLPar("MSY" =refpts(eq)[c("msy"),"yield",1,drop=T])),x=rep(c(2001),4))+
theme_bw()+xlab("Year")+theme(legend.position="bottom")
[1] "maxfbar has been changed to accomodate new plusgroup"
Figure 3 Operating model for 3 growth curves
save(eq,om,file="~/Desktop/tmp/t.RData")
load("~/Desktop/tmp/t.RData")
## Calculate reference points with seasonality
source("~/Desktop/flr/FLCandy/R/seasonalise.R")
#source("~/Desktop/flr/FLCandy/R/wtInterp.R")
## take South parameters
eq=iter(eq,3)
om=iter(om,3)
om4=seasonalise(window(iter(om,3), end=dims(om)$maxyear-1))
[1] "maxfbar has been changed to accomodate new plusgroup"
range(om4)[c("minfbar","maxfbar")]=39
mat(om4)[,,,c(1,3:4)]=NA
## make years+seasons into psuedo+years #######################################
## i.e. "annualize" ############################################################
harvest(om4)[1] =0
harvest(om4)[,,,1] =harvest(om4)[,,,1]*0.2e0
harvest(om4)[,,,2:3]=harvest(om4)[,,,2:3]*1e-6
harvest(om4)[,,,4] =harvest(om4)[,,,4]*0.8e0
stk=qapply(om4, function(x) {
if (dim(x)[4]==1) return(x)
if (dim(x)[1]==1) {
x=x[,,,1]
return(x)}
dnms=dimnames(x)
dnms[[4]]=1
dnms[[1]]=(as.numeric(rep(dnms[[1]],each=dim(x)[4]))+seq(0,1,length.out=dim(x)[4]+1)[-(dim(x)[4]+1)])*dim(x)[4]
FLQuant(c(aperm(x,c(4,1,2,3,5:6))),dimnames=dnms,units=units(x))})
## tidy up
range(stk)["min"]=min(as.numeric(dimnames(m(stk))[[1]]))
range(stk)["max"]=max(as.numeric(dimnames(m(stk))[[1]]))
if (!is.na(range(stk)["plusgroup"]))
range(stk)["plusgroup"]=range(stk)["max"]
range(stk)[c("minfbar","maxfbar")]=range(stk)[c("minfbar","maxfbar")]*dim(m(om4))[4]
range(stk)["minfbar"]=range(stk)["maxfbar"]-3
## Calculate Reference points ##################################################
eq4=FLBRP(stk[-1],sr=list(model=model(eq),params=params(eq)))
eq4=brp(eq4
)
fbar(eq4)=FLQuant(seq(0,refpts(eq4)["crash","harvest",1]/refpts(eq4)["msy","harvest",1],length.out=101))%*%refpts(eq4)["msy","harvest",1]
stk=qapply(om4, function(x) {
if (dim(x)[4]==1) return(x)
if (dim(x)[1]==1) {
x=x[,,,1]
return(x)}
dnms=dimnames(x)
dnms[[4]]=1
dnms[[1]]=(as.numeric(rep(dnms[[1]],each=dim(x)[4]))+seq(0,1,length.out=dim(x)[4]+1)[-(dim(x)[4]+1)])*dim(x)[4]
FLQuant(c(aperm(x,c(4,1,2,3,5:6))),dimnames=dnms,units=units(x))})
## tidy up
range(stk)["min"]=min(as.numeric(dimnames(m(stk))[[1]]))
range(stk)["max"]=max(as.numeric(dimnames(m(stk))[[1]]))
if (!is.na(range(stk)["plusgroup"]))
range(stk)["plusgroup"]=range(stk)["max"]
range(stk)[c("minfbar","maxfbar")]=range(stk)[c("minfbar","maxfbar")]*dim(m(om4))[4]
range(stk)["minfbar"]=range(stk)["maxfbar"]-3
m(stk)=m(stk) #*1.5
## Calculate Reference points ##################################################
eq4=FLBRP(stk[-1],sr=list(model=model(eq),params=params(eq)))
eq4=brp(eq4)
fbar(eq4)=FLQuant(seq(0,refpts(eq4)["crash","harvest",1]/refpts(eq4)["msy","harvest",1],length.out=101))%*%refpts(eq4)["msy","harvest",1]
dat=as.data.frame(FLQuants(eq4,"Mass"=stock.wt,"M1"=M1Fn,"M2"=M2Fn,"Mat"=mat,"Selection"=catch.sel))
dat=transform(dat,Season=ac(age-((age-1)%/%4)*4))
#dat=subset(dat,!(Season%in%c(1,3:4)&qname=="Mat"))
## vectors-at-age
ggplot(dat)+
geom_line(aes(age,data,col=Season))+
facet_wrap(~factor(qname,levels=c("Mass","M1","Selection","Mat","M2")),scale="free")+
scale_x_continuous(limits=c(0,40))+
theme(legend.position="bottom")
Figure 4 Vectors from seasonal operating model
fbar(eq4)=FLQuant(c(seq(0,0.1,length.out=100),seq(0.1,1,length.out=101)[-1]))*max(refpts(eq4)["crash","harvest"])
plot(eq4,refpts="msy")
Figure 5 Equilibrium from seasonal operating model
sr.par=FLPar(NA,
dimnames=list(params=dimnames(params(eq))$params, season=1:4, iter=1))
sr.par[,2]=params(eq)[,1]
sr.par[1,2]=sr.par[1,2]
sr.par[2,2]=sr.par[2,2]
sr4=predictModel(model=model(eq), params=sr.par)
#lnds=1+0*FLQuant(sample(lnd$landings,100,T),dimnames=list(year=2:100))/mean(lnd$landings)
#ctc=FLQuant(0,dimnames=dimnames(catch(om4)[,2:100]))
#ctc[,,,1:2]=0
#ctc[,,,3]=ctc[,,,3]%*%lnds*50
#ctc[,,,4]=ctc[,,,4]%*%lnds*200
#control=as(FLQuants("catch"=ctc),"fwdControl")
f=c(seq(0,1,length.out=51),
seq(1,c(refpts(eq4)["crash","harvest",1]%/%refpts(eq4)["msy","harvest",1]),length.out=51)[-1])*
refpts(eq4)["msy","harvest",1,drop=T]*4
f=FLQuant(rep(f,each=prod(dim(om4)[c(2,4)])),dimnames=list(year=dimnames(om4)$year,season=1:4,iter=seq(101)))
harvest(om4)[1]=0
f[,,,2:3]=f[,,,2:3]*1e-6
f[,,,1] =f[,,,1]*0.2
f[,,,4] =f[,,,4]*0.8
control=as(FLQuants("f"=f[,-1]),"fwdControl")
mat(om4)[is.na(mat(om4))]=0
om4=fwd(propagate(om4,101),control=control,sr=sr4)
print(sr4)
An object of class "FLQuants": EMPTY
model:
rec ~ a * ssb/(b + ssb)
<environment: R_EmptyEnv>
params:
An object of class "FLPar"
season
params 1 2 3 4
a NA 12951.2 NA NA
b NA 57.2 NA NA
units: NA
mod=model.frame(FLQuants(om4[,40],
Rec =function(x) rec(x)[,,,2],
SSB =function(x) ssb(x)[,,,2],
Catch =function(x) apply(catch(x),c(2,6), sum),
F =function(x) apply(fbar(x), c(2,6),mean),
Forage=function(x) apply(stock.wt(x)%*%(m(x)-0.1)/(z(x))%*%(1-exp(-z(x))),c(2,6),sum)))
a=ggplot(mod)+geom_line(aes(SSB,Rec))
b=ggplot(mod)+geom_line(aes(SSB,Catch))
c=ggplot(mod)+geom_line(aes(F,Catch))
d=ggplot(mod)+geom_line(aes(F,SSB))
e=ggplot(mod)+geom_line(aes(F,Forage))
ggarrange(d,a,c,b,e,
widths = c(5, 5), heights = c(5, 5),
nrow = 3, ncol = 2,
labels = c("A","B","C","D"),
common.legend = TRUE)+theme_bw()
Figure 6 Equilibrium curves from running seasonal operating model to equilibrium
p=plot(window(om4,start=1960,end=2110),metrics=list(
Rec =function(x) rec(x)[,,,3],
SSB =function(x) ssb(x)[,,,2],
Catch=function(x) apply(catch(x),c(2,6), sum),
F =function(x) apply(fbar(x), c(2,6),mean)))+
geom_flpar(data=FLPars(Rec =FLPar("Rmsy"=refpts(eq4)[c("msy"),"rec",1,drop=T])*0.6718, #M at age 0 season 2
SSB =FLPar("Bmsy"=refpts(eq4)[c("msy"),"ssb",1,drop=T]),
F =FLPar("Fmsy"=refpts(eq4)[c("msy"),"harvest",1,drop=T]),
Catch=FLPar("MSY" =refpts(eq4)[c("msy"),"yield",1,drop=T])),x=rep(c(1960),4))+
theme_bw()+xlab("Year")+theme(legend.position="bottom")
p
msyIter=(1:101)[c(apply(catch(om4[,100]),6,sum)==max(apply(catch(om4[,100]),6,sum)))]
fbar(eq4)=fbar(eq4)[,1];fbar(eq4)[]=refpts(eq4)["msy","harvest",1]
dt2=qp(om4[,dim(om4)[2],,,,msyIter],eq4)
ggplot(dt2)+
geom_line(aes(FLBRP,FLStock,col=season))+
facet_wrap(~What,scale="free")
library(FLCore)
library(FLBRP)
library(FLife)
eq=lhEql(lhPar(FLPar(linf=100)))
om=as(eq,"FLStock")
om=fwd(om,f=fbar(om)[,-1],sr=eq)
om=FLCore:::adjust( om )
catch.n( om )
landings.n( om )
discards.n( om )
n =100
object=m(ns)
devs =mDev
rBoot<-function(n,object,devs){
object=propagate(object,n)
rtn=mdply(data.frame(year=dimnames(object)$year), function(year){
object[,year]%*%devs[,,,,,sample(dimnames(devs)$iter,n,TRUE)]})
rtn}
dev=rBoot(100,m(ns),mDev)
library(RTseries)
fn<-function(x) c(acf(x,plot=T)$acf)
dat=FLCohort((stock.wt(ple4)%-%apply(stock.wt(ple4),1,mean))%/%apply(stock.wt(ple4),1,mean))
dt2=apply(dat,2,function(x) c(acf(x,plot=F,na.action=na.pass)$acf))
dt3=apply(as(dt2,"FLQuant"),1,mean,na.rm=T)[,drop=T]
data.table(as.data.frame(FLCohort(dat)))[!is.na(data), function(x) c(acf(x,plot=T)$acf), by=.(cohort)]
om4=iter(om4,-1)
## Cohort
rlcohort<-function(x, sd=1, b=0){
as(apply(FLCohort(x), c(2:3,5:6), function(y)
FLife:::noiseFn(prod(dim(y)),sd=sd,b=b)),"FLQuant")
}
rlcohort<-function(x, sd=1, b=0){
res=adply(FLCohort(x), c(2:3,5:6), function(y)
data.frame(expand.grid(dimnames(y)[2:1]),data=FLife:::noiseFn(prod(dim(y)),sd=sd,b=b)))
as.FLQuant(mutate(res,year=an(ac(cohort))+an(ac(age)))[,-1])[,dimnames(x)$year]
}
devM1 =rlnoise(dim(om4)[6],m(om4)[1,,,1,,1]%=%0,sd=0.3,b=0.0)
devM2 =rlcohort(propagate(m( om4),dim(om4)[6]), sd=0.3,b=0.6)
devShock=FLQuant(sample(c(rep(0,39),-log(1-0.75)),prod(dim(rec(om4)))*100,T),dimnames=dimnames(m(om4)[1]))
ggplot(iter(devM2,1))+
geom_point(aes(year,age,size=abs(data),col=(data>0)))+
facet_grid(season~.)
devM2=exp(devM2)
devRecs =rlnorm(dim(om4)[6],iter(rec(om4),1)%=%0,0.3)
burnin =propagate(window(iter(om4,51),start=1980,end=2040),100)
burnin=FLStocks("Base" =burnin,
"M Cohort"=burnin,
"M Annual"=burnin,
"M Shock" =burnin)
devM1 =devM2[ ,dimnames(burnin[["Base"]])$year]
devM2 =devM2[ ,dimnames(burnin[["Base"]])$year]
devShock=devShock[,dimnames(burnin[["Base"]])$year]
Ms=list("Base" =FLQuants(M1=M1Fn(burnin[["Base"]]), M2=M2Fn(burnin[["Base"]])),
"M Cohort"=FLQuants(M1=M1Fn(burnin[["Base"]]), M2=M2Fn(burnin[["Base"]])%*%devM2),
"M Annual"=FLQuants(M1=M1Fn(burnin[["Base"]])%*%devM1, M2=M2Fn(burnin[["Base"]])%*%devM2),
"M Shock" =FLQuants(M1=M1Fn(burnin[["Base"]])%*%devM1%+%devShock,M2=M2Fn(burnin[["Base"]])%*%devM2))
m(burnin[["M Cohort"]])=Ms[["M Cohort"]][["M1"]]%+%Ms[["M Cohort"]][["M2"]]
m(burnin[["M Annual"]])=Ms[["M Annual"]][["M1"]]%+%Ms[["M Annual"]][["M2"]]
m(burnin[["M Shock"]]) =Ms[["M Shock" ]][["M1"]]%+%Ms[["M Shock" ]][["M2"]]
targetF=FLQuant(c(refpts(eq4)["msy","harvest",1])/4,dimnames=dimnames(fbar(burnin[["Base"]])[,-1]))
control=as(FLQuants("f"=targetF),"fwdControl")
burnin =FLStocks(llply(burnin, function(x) window(fwd(x,control=control,sr=sr4,residuals=devRecs),start=1981)))
control=as(FLQuants("f"=targetF[,ac(2021:2040)]*1.5),"fwdControl")
kick =FLStocks(llply(burnin, function(x) window(fwd(x,control=control,sr=sr4,residuals=devRecs),start=2001)))
fnl=c(seq(0,1,length.out=51)[-51],seq(1,4,length.out=50))
fnl=t(maply(fnl, function(x) seq(1,x,length.out=20)))
fnl=FLQuant(c(fnl),dimnames=dimnames(targetF[,ac(2021:2040),,1]))
targetF[,ac(2021:2040)]=targetF[,ac(2021:2040)]%*%fnl
control =as(FLQuants("f"=targetF[,ac(2021:2040)]),"fwdControl")
ftrend=FLStocks(llply(burnin, function(x) window(fwd(x,control=control,sr=sr4,residuals=devRecs),start=2001)))
library(ggplotFL)
p=plot(window(burnin,start=2010),metrics=list(
Rec =function(x) rec(x)[,,,3],
SSB =function(x) ssb(x)[,,,2],
Catch =function(x) apply(catch(x),c(2,6), sum),
F =function(x) apply(fbar(x), c(2,6),mean),
Forage=function(x) apply(stock.wt(x)%*%stock.n(x)%*%(m(x)-0.1)/(z(x))%*%(1-exp(-z(x))),c(2,6),sum)),iter=1)+
geom_flpar(data=FLPars(Rec =FLPar("Rmsy"=refpts(eq4)[c("msy"),"rec",1,drop=T])*0.6718, #M at age 0 season 2,
SSB =FLPar("Bmsy"=refpts(eq4)[c("msy"),"ssb",1,drop=T]),
F =FLPar("Fmsy"=refpts(eq4)[c("msy"),"harvest",1,drop=T])/4,
Catch=FLPar("MSY" =refpts(eq4)[c("msy"),"yield",1,drop=T])*4),
x=rep(ISOdate(2012,1,1),4))+
theme_bw()+xlab("Year")+theme(legend.position="bottom")
p
p1=plot(window(FLStocks("Status Quo"=burnin[[1]],"Kicked"=kick[[1]],
"F Trend"=ftrend[[1]]),start=2010),metrics=list(
Rec =function(x) rec(x)[,,,3],
SSB =function(x) ssb(x)[,,,2],
Catch=function(x) apply(catch(x),c(2,6), sum),
F =function(x) apply(fbar(x), c(2,6),mean),
Forage=function(x) apply(stock.wt(x)%*%stock.n(x)%*%(m(x)-0.1)/(z(x))%*%(1-exp(-z(x))),c(2,6),sum)),iter=1)+
geom_flpar(data=FLPars(Rec =FLPar("Rmsy"=refpts(eq4)[c("msy"),"rec",1,drop=T])*0.6718, #M at age 0 season 2,
SSB =FLPar("Bmsy"=refpts(eq4)[c("msy"),"ssb",1,drop=T]),
F =FLPar("Fmsy"=refpts(eq4)[c("msy"),"harvest",1,drop=T])/4,
Catch=FLPar("MSY" =refpts(eq4)[c("msy"),"yield",1,drop=T])*4),
x=rep(ISOdate(2012,1,1),4))+
theme_bw()+xlab("Year")+theme(legend.position="bottom")
p1
save(burnin,kick,ftrend,
eq4,sr4,pars,
Ms,
devM1,devM2,devShock,
devRecs,
file=file.path("../data","omScenarios.RData"))
nseason=4
nunit =2
## Set up OM
par=lhPar(FLPar(linf=40))
## Adjust parameters
dmns=dimnames(par)
dmns[["season"]]=seq(nseason)
dmns[["unit"]] =seq(nunit)
dmns=dmns[c(1,4,3,2)]
par42=FLPar(array(c(par),dim =laply(dmns,length),
dimnames=dmns))
## Units
par42[,2]=lhPar(FLPar(linf=30))
## Seasons
par42["t0"]=par42["t0"]-seq(0,dim(par42)[3]-1)/dim(par42)[3]
eq=lhEql(par)
om=as(eq,"FLStock")
om=window(fwd(om,fbar=m(om)[4,-1],sr=eq),start=19)
[1] "maxfbar has been changed to accomodate new plusgroup"
dimnames(om)$year =an(dimnames(om)$year)+1980
range(om)["plusgroup"]=range(om)["maxage"]
om=fwd(om,fbar=window(fbar(om),start=2021)%*%
FLQuant(seq(1,10,length.out=61)),sr=eq)
## Make it seasonal with units
source("~/Desktop/active/FLCandy/R/seasonalise.R")
om42=seasonalise(om)
om42=FLCore:::expand(om, unit=seq(nunit),season=seq(nseason))
source("~/Desktop/flr/FLCandy/R/saveTmp.R")
## ALKs
res =list()
ak42=list()
for (i in seq(dim(par42)[2])){
for (j in seq(dim(par42)[3]))
res[[j]]=invALK(par42[,i,j],cv=0.1,age=an(dimnames(om)$age),bin=0.1)
ak42[[i]]=FLPars(res)}
An object of class "FLPar"
, , season = 1
unit
params 1 2
linf 40 30
, , season = 2
unit
params 1 2
linf 40 30
, , season = 3
unit
params 1 2
linf 40 30
, , season = 4
unit
params 1 2
linf 40 30
units: NA
R version 4.2.1 (2022-06-23)
knitr: 1.43
FLCore: 2.6.19.9105
FLBRP: 2.5.9.9022
FLasher: 0.7.1.9221
FLife: 3.4.0
ggplotFL: 2.7.0.9132
FLCandy: 0.1.0
rfishbase: 3.1.9.99
SPMpriors: 1.0.4
FishLife: 3.0.0
plyr: 1.8.9
dplyr: 1.1.4
reshape: 0.8.9
ggplot2: 3.4.4
ggpubr: 0.6.0
ggcorrplot: 0.1.3
GGally: 2.1.2
Compiled: Sun Jan 21 15:28:48 2024