library(FLash) library(FLAssess) data(ple4) mn.rec<-exp(mean(log(stock.n(ple4)[1,ac(1990:2000)]))) ## set up 5 years for projection ple4<-stf(ple4,nyrs=5) ## absolute options ############################################################ ## F target #################################################################### ## set up fwdTarget object to determine projection options Target<-fwdTarget(year=2002:2006,value=.3,quantity="f") proj.f <-fwd(ple4, Target, sr.model="mean", sr.param=mn.rec) ## check fbar(proj.f)[,ac(2000:2006)] ## Catch target ################################################################ ## set up fwdTarget object to determine projection options Target<-fwdTarget(year=2002:2006,value=50000,quantity="catch") proj.catch <-fwd(ple4, Target, sr.model="mean", sr.param=250000) ## check fbar(proj.catch)[,ac(2000:2006)] ## SSB target ################################################################## # SSB targets are at end of year if harvest.spwn =0.0, otherwise in year Target<-fwdTarget(year=2001:2005,value=2:6+200000,quantity="ssb") proj.ssb <-fwd(ple4, Target, sr.model="mean", sr.param=250000) ## check ssb(proj.ssb)[,ac(2000:2006)] plot(FLStocks(proj.f,proj.catch,proj.ssb)) ## SSB in year ple4@harvest.spwn[]<-0.5 proj.ssb.rel)[,ac(2000:2006)] ## relative options ############################################################ # set catch to 1, 0.5, 0.2, 0.1, .1 times that in 1995 Target<-fwdTarget(year=2002:2006,rel=2001:2005,value=c(.85,.85,.85,.85,.85),quantity="catch") proj.ssb.rel <-fwd(ple4, Target, sr.model="mean", sr.param=250000) ## check computeCatch(proj.ssb.rel)[,ac(2001:2006)] ################################################################################ HCRs ################################################################################ ## Recovery plans ############################################################## # 30% reduction in F until F=F0.1 f0.1<-refpts(brp(FLBRP(ple4)))["f0.1","harvest"] Target<-fwdTarget(year=rep(2002:2006,each=2),rel=rep(c(0,NA),5)+rep(2001:2005,each=2),value=c(.70,NA),min=c(NA,f0.1),quantity="f") proj.recovery.1 <-fwd(ple4, Target, sr.model="mean", sr.param=250000) plot(FLStocks(proj.recovery.1)) # 10% increase in SSB until F=F0.1 f0.1<-refpts(brp(FLBRP(ple4)))["f0.1","harvest"] Target<-fwdTarget(year=rep(2002:2006,each=2),rel=rep(c(0,NA),5)+rep(2001:2005,each=2),value=c(1.10,NA),min=c(NA,f0.1),quantity=c("ssb","f")) ## bug in if can't achieve SSB F ->0 proj.recovery.2 <-fwd(ple4, Target, sr.model="mean", sr.param=250000) plot(FLStocks(proj.recovery.1,proj.recovery.2)) ################################################################################ ## rebuild to BMSY ############################################################# ## Rebuilding target, a proxy for BMSY ## i.e. to 0.5*BMSY by 2006 ssbTarget<-refpts(brp(FLBRP(ple4)))["f0.1","ssb"]*mn.rec*.5 ## function to minimise f <- function(x,stk,ssb.target,target,sr.model,sr.param) { target[,"value"]<-x ssb.<-c(ssb(fwd(stk,target,sr.model=sr.model,sr.param=sr.param))[,"2006"]) return((ssb.-ssb.target)^2) } ## Recover stock to BMY in 2006 with a constant F strategy Target<-fwdTarget(year=2002:2006,value=.5,rel=2001,quantity="f") xmin <-optimize(f, c(0.1, 1.0), tol = 0.0000001, stk=ple4, ssb.target=ssbTarget, target=Target, sr.model="mean", sr.param=mn.rec) Target[,"value"]<-xmin$minimum ple4.rebuild.f <-fwd(ple4, Target,sr.model="mean",sr.param=mn.rec) plot(FLStocks(ple4.rebuild.f)) ## Recover stock to BMY in 2020 with a constant catch strategy Target<-fwdTarget(year=2002:2006,value=.5,rel=2001,quantity="catch") xmin <-optimize(f, c(0.1, 1.0), tol = 0.0000001, stk=ple4, ssb.target=ssbTarget, target=Target, sr.model="mean", sr.param=mn.rec) Target[,"value"]<-xmin$minimum ple4.rebuild.catch <-fwd(ple4, Target,sr.model="mean",sr.param=mn.rec) plot(FLStocks(ple4.rebuild.f,ple4.rebuild.catch),yrs=1996:2006) ## Recover stock to 0.5*BMY in 2006 with an F reduction strategy Target<-fwdTarget(year=2002:2006,value=.9,rel=2001:2005,quantity="f") xmin <-optimize(f, c(0.1, 1.0), tol = 0.0000001, stk=ple4, ssb.target=ssbTarget, target=Target, sr.model="mean", sr.param=mn.rec) Target[,"value"]<-xmin$minimum ple4.rebuild.fred <-fwd(ple4, Target,sr.model="mean",sr.param=mn.rec) plot(FLStocks(ple4.rebuild.fred,ple4.rebuild.f,ple4.rebuild.catch),yrs=1996:2006) ################################################################################ ## Assessment upto and including 2001 data(ple4) black.bird <-stf(ple4,nyrs=2) # set courtship and egg laying in Autumn black.bird@m.spwn[] <-0.66 black.bird@harvest.spwn[]<-0.66 # assessment is in year 2002, set catch constraint in 2002 and a first guess for F in 2003 target <-fwdTarget(year=2002:2003,value=c(85000,.5),quantity=c("catch","f")) black.bird <-fwd(black.bird, target, sr.model="mean", sr.param=25000) # HCR specifies F=0.1 if ssb<100000, F=0.5 if ssb>300000 # otherwise linear increase as SSB increases min.ssb<-100000 max.ssb<-300000 min.f <-0.1 max.f <-0.5 # slope of HCR a. <-(max.f-min.f)/(max.ssb-min.ssb) b. <-min.f-a.*min.ssb # plot of HCR plot(c(0.0,min.ssb,max.ssb,max.ssb*2),c(min.f,min.f,max.f,max.f),type="l",ylim=c(0,max.f*1.25),xlim=c(0,max.ssb*2)) ## find F through iteration t. <-999 i <-0 while (abs(target[2,"value"]-t.)>10e-6 & i<50) { t.<-target[2,"value"] ## save last value of F # calculate new F based on SSB last iter target[2,"value"]<-a.*c(ssb(black.bird)[,"2003"])+b. black.bird<-fwd(black.bird, target, sr.model="mean", sr.param=25000) # 'av a gander points(c(ssb(black.bird)[,"2003"]),c(target[2,"value"]),cex=1.25,pch=19,col=i) print(c(ssb(black.bird)[,"2003"])) print(c(target[2,"value"])) i<-i+1 } # F bounds black.bird <-fwd(black.bird, target, sr.model="mean", sr.param=25000) plot(FLStocks(black.bird)) ################################################################################ ## TAC constraint ############################################################## target<-fwdTarget(year =c(iyr,rep(iyr+1,2),iyr+2), quantity=c(rep( "f",2),"catch","f"), value =c(mean(Fsq),mean(c(pmin(pmax(Fsq,F0.1),Fpa))),NA,mean(Fsq)), min =c(NA,NA,mean(TAC*(1.0-TAC.bound)),NA), max =c(NA,NA,mean(TAC*(1.0+TAC.bound)),NA)) ################################################################################ ## PA limit check ############################################################## ## If SSB+2 < Bpa ## ensure SSBy+2 >= SSBy+1 and Fy+2=<Fy+1 target<-fwdTarget(year =c(iyr+1,iyr+1), quantity=c("ssb", "f"), rel =c(iyr, iyr), min =c( 1, NA), max =c( NA, 1)) ################################################################################ ## Monte Carlo versions ######################################################## ## F status quo Fsq <-c(fbar(stk)[,yr.now]) ## target options target<-fwdTarget(year =c(iyr,rep(iyr+1,2),iyr+2), quantity=c(rep( "f",2),"catch","f"), value =c(mean(Fsq),mean(c(pmin(pmax(Fsq,F0.1),Fpa))),NA,mean(Fsq)), min =c(NA,NA,mean(TAC*(1.0-TAC.bound)),NA), max =c(NA,NA,mean(TAC*(1.0+TAC.bound)),NA)) cat("\n===================","Target"," ====================\n") print(target[,1:5]) cat("==================================================\n") ## targets by iter value <-array(NA,c(dim(target)[1], dims(stk)$iter)) min <-value max <-value value[1:2,]<-c(Fsq) value[2,]<-c(pmin(pmax(Fsq,F0.1),Fpa)) min[3,] <-c(TAC*(1.0-TAC.bound)) max[3,] <-c(TAC*(1.0+TAC.bound)) ## run mn.rec<-exp(mean(log(stk@stock.n[1,yr.recruits]),na.rm=TRUE))/exp(var(log(stk@stock.n[1,yr.recruits]),na.rm=TRUE)/2) t.<-system.time(stk<-fwd(stk,target[-3,],value=value[-3,], min=min[-3,], max=max[-3,], sr.model="mean",sr.param=mn.rec))