#Exercises Yield and bio-economic model II - capture fishery
#as is Per Recruit model
#Equilibrium stock for single age of first capture. Projection matrix method.
#tc and F based computations; fishing effort specification in output

#Author: Jorge Santos, Tromsų, v2013-2014

#Suitable to finely adjust the most appropriate fishing strategies
#Install necessary libraries (packages)
#Adjust 1ary and 2ary input according to species biology and fishery
#To run script highlight the whole script and Run (or ctrl-enter)
#Aggregated output in Results.mat file and 4 figures in terms of F 
                # and the same 4 figures in terms of f (effort)

#needs
library(Hmisc)

#age structured growth - INPUT
R=1                           #Number of recruits (constant) - Y/R => R=1
M=?                           #Inst rate natural mortality (year-1)
Linf=?                        #Linf von Bertalanffy (cm)
K=?                           #K von Bertalanffy
t0=?                          #t0 von Bertalanffy (y) (set to 0)
a=?                           #Constant weight length (g/cm-3)
b=?                           #Power weight length relation (unitless)
tr=?                          #Age recruitment to fishery (y)
ts=?                          #Age maturity (knife-edge) (y)
#fishery
tc=?                          #Age 1st capture (knife-edge) (y)
q=?                           #catchability (=F/(f^elast)); stock and fishery specific
elast=?                       #elasticity (exponent); denotes aggregation of stock (elast<1) or search time (>1); specific
#economics
price1=?                      #Sales price constant ($/kg), price=price1*price2*ln(Weight,kg)
price2=?                      #Sales price exponent (if =0, then price constant)
cost=?                        #Production (capture) cost ($/f unit); variable cost only


 
#time steps (2ary input)
ages=30                       #Maximum age of species (y), be generous (large ages)
years=30                      #Projection horizon (y)
step=1                        #Time step of analysis [0,1] (y), = number of cohorts/year
                              #=age groups/year; R=1 each step

###############
#Simulation control 

F.vec=c(seq(0,2,by=0.05))
Results.mat=matrix(NA,nrow=19,ncol=length(F.vec))
rownames(Results.mat)=c("F", "Survivors", "TotBio","SpawnBio",
                        "CatchN","Yield", "Ymax", "Fmax", "Yopt", "F_0.1",
                        "Costs","Revenue","Profit","MEY","Fmey",
                        "f","fYmax","fF0.1","fMEY")
for (j in 1:length(F.vec)) {
F=F.vec[j]
Results.mat["F",j]=F
  
#CALCULATIONS

#time steps
time=years/step
Age=c(seq(tr,(tr+ages),by=step))

#Survival and mortality (instantaneous and fractional)
Fage=c(rep(0,times=((tc-tr)/step)),rep(F*step,times=(tr+ages-tc)/step+1))
Mage=c(rep(M*step,times=(ages/step)+1))
Zage=numeric((ages/step)+1)
for(i in 1:(ages/step+1)) Zage[i]= Mage[i]+Fage[i]
SurvAge=numeric((ages/step)+1)
for(i in 1:(ages/step+1)) SurvAge[i]=exp(-Zage[i])

#Growth and weight (von Bertalanffy)
Wage=numeric((ages/step)+1)
for(i in 1:(ages/step+1)) Wage[i]=a*(Linf*(1-exp(-K*(Age[i]-t0))))^b


#(Lesley) A-matrix, N-matrix
A.mat=matrix(0,nrow=((ages/step)+1),ncol=((ages/step)+1))
A.mat[1,1]=1
for (i in 1:(ages/step)) A.mat[i+1,i]=SurvAge[i]
#Initial abundance vector
N0=c(R,rep(0,times=ages/step))
#Projection matrices
N.proj=matrix(0,nrow=nrow(A.mat),ncol=time+1)
N.proj[,1]=t(N0)
for (i in 1:time)N.proj[,i+1]=A.mat%*%N.proj[,i]
#matplot(0:time+1,t(N.proj),type="l",lty=1:3,
#        col=1,ylab="Age abundance N", xlab="Year")
N.totals=apply(N.proj,2,sum,na.rm=TRUE)
N=N.totals[time+1]
Results.mat["Survivors",j]=N

#Biomass matrix
B.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
B.proj=N.proj*Wage
B.totals=apply(B.proj,2,sum,na.rm=TRUE)
B=B.totals[time+1]
Results.mat["TotBio",j]=B

#Spawning biomass matrix
SB.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
for (i in ((ts-tr)/step+1):((ages/step)+1)) SB.proj[i,]=N.proj[i,]*Wage[i]
SB.totals=apply(SB.proj,2,sum,na.rm=TRUE)
SB=SB.totals[time+1]
Results.mat["SpawnBio",j]=SB

#Catch matrix
C.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
for (i in 1:(ages/step+1))C.proj[i,]=N.proj[i,]*(1-exp(-Fage[i]))
C.totals=apply(C.proj,2,sum,na.rm=TRUE)
Catch=C.totals[time+1]
Results.mat["CatchN",j]=Catch
  
#Yield matrix
Y.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
for (i in 1:(ages/step+1))Y.proj[i,]=C.proj[i,]*(Wage[i])/1000    #NEW (/1000) yield in kgs - not useful in per Recruit analysis
Y.totals=apply(Y.proj,2,sum,na.rm=TRUE)
Yield=Y.totals[time+1]
Results.mat["Yield",j]=Yield

#Economics (Price at size ($/kg), Profit=total cost-revenue)
Pricsize=numeric((ages/step)+1)

for(i in 1:(ages/step+1)) Pricsize[i]=price1+price2*log(Wage[i]/1000)
#Results.mat["Costs",]=Results.mat["F",]*cost                          #OLD costs function of F
Results.mat["Costs",]=(Results.mat["F",]/q)^(1/elast)*cost            #NEW costs function of f (=effort)


Revenue.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
for (i in 1:(ages/step+1))Revenue.proj[i,]=Y.proj[i,]*(Pricsize[i])*1000 #NEW *1000 for per recruit version
R.totals=apply(Revenue.proj,2,sum,na.rm=TRUE)
Revenue=R.totals[time+1]
Results.mat["Revenue",j]=Revenue

Results.mat["Profit",j]=
  (Results.mat["Revenue",j]-Results.mat["Costs",j])
Profit=max(Results.mat["Profit",j])


}

#transpose results matrix
Results.mat=t(Results.mat)

#Ymax, Fmax
iYmax=which.max(Results.mat[,"Yield"])
Ymax=round(max(Results.mat[,"Yield"]),0)
Results.mat[iYmax,"Ymax"]=Ymax
Fmax=max(Results.mat[iYmax,"F"])
Results.mat[iYmax,"Fmax"]=Fmax

#Yopt,F_0.1 (the F must be equally spaced; estimated from differences, not slopes)
Ydiffs=c(rep(NA, times=length(F.vec)))
for (i in 1:(length(F.vec)-1))
  Results.mat[i+1,"Yield"]-Results.mat[i,"Yield"]->Ydiffs[i+1]
Yratio=c(rep(NA, times=length(F.vec)))
for (i in 1:(length(F.vec)-1))
  Ydiffs[i+1]/max(Ydiffs[!is.na(Ydiffs)])->Yratio[i+1]
iYratio=which.max(Yratio<=0.1)
Yopt=max(round(Results.mat[iYratio,"Yield"],0))
Results.mat[iYratio,"Yopt"]=Yopt
F_0.1=max(Results.mat[iYratio,"F"])
Results.mat[iYratio,"F_0.1"]=F_0.1

#NEW expressed in terms of fishing effo                            #NEW
Results.mat[,"f"]=round((Results.mat[,"F"]/q)^(1/elast))
Results.mat[,"fYmax"]=round((Results.mat[,"Fmax"]/q)^(1/elast))
fYmax=max(Results.mat[,"fYmax"],na.rm = TRUE)                      
Results.mat[,"fF0.1"]=round((Results.mat[,"F_0.1"]/q)^(1/elast))
fF0.1=max(Results.mat[,"fF0.1"],na.rm = TRUE)

#Economics
#Results.mat[,"Costs"]=Results.mat[,"F"]*cost
#Results.mat[,"Revenue"]=round(Results.mat[,"Yield"]*price,0)
#Results.mat[,"Profit"]=round(Results.mat[,"Revenue"]-Results.mat[,"Costs"],0)
iMEY=which.max(Results.mat[,"Profit"])
MEY=max(Results.mat[,"Profit"])
FMEY=max(Results.mat[iMEY,"F"])
Results.mat[iMEY,"Fmey"]=FMEY
Results.mat[iMEY,"MEY"]=MEY
Results.mat[,"fMEY"]=round((Results.mat[,"Fmey"]/q)^(1/elast))   #NEW (f-based)
fMEY=max(Results.mat[,"fMEY"],na.rm = TRUE)                      #NEW

#subset economics matrix                                         
EconResults.mat=subset(Results.mat,select=Costs:Profit)


#PLOTS
#needs
library(Hmisc)

#Plot, Yield full recruit (F-based)
plot(Results.mat[,"F"], Results.mat[,"Yield"],type="l",lty=1:3,        
        col=1,ylab="Yield/R (kg)", xlab="F (fishing mortality, Year-1)")    
minor.tick(nx=5, ny=5, tick.ratio=0.5)
legend("bottomright",legend=c(paste("yield at tc",tc,"y"),
      paste("Fmax=",Fmax),paste("F0.1=",F_0.1)), lty=1:3,col=1,bty="n")
title(main=paste("Yield per Recruit in equilibrium, at tc", tc,"y"),cex.main=0.8)
abline(v=Fmax,lty=2)
abline(v=F_0.1,lty=3)


#Plot, Spawning biomass full recruit (F-based)
plot(Results.mat[,"F"], Results.mat[,"SpawnBio"],type="l",lty=1:3,     
     col=1,ylab="SB/R (g)", xlab="F (fishing mortality, Year-1)")
minor.tick(nx=5, ny=5, tick.ratio=0.5)
legend("topright",legend=c(paste("Spawners at tc",tc,"y"),
     "MBAL40%","MBAL20%"), lty=1:3,col=1,bty="n")
title(main="Spawning biomass/R in equilibrium",cex.main=0.8)
abline(h=0.4*Results.mat[1,"SpawnBio"],lty=2)                         
abline(h=0.2*Results.mat[1,"SpawnBio"],lty=3)  

#Plot, Total biomass recruited stock full recruit (F-based)
plot(Results.mat[,"F"], Results.mat[,"TotBio"],type="l",lty=1:3,      
     col=1,ylab="Biomass / R (g)", xlab="F (fishing mortality, Year-1)")
minor.tick(nx=5, ny=5, tick.ratio=0.5)
legend("topright",legend=c(paste("Stock at tc",tc,"y"),
                           "Stock 50%","Stock 40%"), lty=1:3,col=1,bty="n")
title(main="Biomass of recruited stock in equilibrium",cex.main=0.8)
abline(h=0.5*Results.mat[1,"TotBio"],lty=2)                           
abline(h=0.4*Results.mat[1,"TotBio"],lty=3)                          

#Plot, Economics full recruit (F based)
matplot(Results.mat[,"F"], EconResults.mat,type="l",lty=1:3,          
     col=1,ylab="Flow / R ($)", xlab="F (fishing mortality, Year-1)")
minor.tick(nx=5, ny=5, tick.ratio=0.5)
legend("bottomright",legend=c("Costs", "Revenue", paste("Profit, Fmey=",FMEY)), 
       lty=1:3,col=1,bty="n")
title(main=paste("Economics of fishery in equilibrium, at tc", tc,"y"),cex.main=0.8)
abline(h=0,lty=4,col=2)
abline(v=FMEY,lty=5,col=4)

##### NEW f-based plots #############

#Plot, Total Yield (f-effort based)
plot(Results.mat[,"f"], Results.mat[,"Yield"],type="l",lty=1:3,
     col=1,ylab="Yield (kg)", xlab="f (fishing effort, unit)")    
minor.tick(nx=5, ny=5, tick.ratio=0.5)
legend("bottomright",legend=c(paste("Yield at tc",tc,"y"),
          paste("f(Fmax)=",fYmax),paste("f(F0.1)=",fF0.1)), lty=1:3,col=1,bty="n")
title(main="Yield/R in equilibrium",cex.main=0.8)
abline(v=fYmax,lty=2)
abline(v=fF0.1,lty=3)

#Plot, Spawning biomass full recruitment
plot(Results.mat[,"f"], Results.mat[,"SpawnBio"],type="l",lty=1:3,  
     col=1,ylab="SB/R (g)", xlab="f (fishing effort, unit)")
minor.tick(nx=5, ny=5, tick.ratio=0.5)
legend("topright",legend=c(paste("Spawners at tc",tc,"y"),
                           "MBAL40%","MBAL20%"), lty=1:3,col=1,bty="n")
title(main="Spawning biomass/R in equilibrium",cex.main=0.8)
abline(h=0.4*Results.mat[1,"SpawnBio"],lty=2)                       
abline(h=0.2*Results.mat[1,"SpawnBio"],lty=3)                       

#Plot, Total biomass full recruitment
plot(Results.mat[,"f"], Results.mat[,"TotBio"],type="l",lty=1:3,     
     col=1,ylab="Biomass/R (g)", xlab="f (fishing effort, unit)")
minor.tick(nx=5, ny=5, tick.ratio=0.5)
legend("topright",legend=c(paste("Stock at tc",tc,"y"),
                           "Stock 50%","Stock 40%"), lty=1:3,col=1,bty="n")
title(main="Biomass of recruited stock in equilibrium",cex.main=0.8)
abline(h=0.5*Results.mat[1,"TotBio"],lty=2)                       
abline(h=0.4*Results.mat[1,"TotBio"],lty=3)                       # (g->tons)


#Plot, Economics full recruitment (f-effort based)                         
matplot(Results.mat[,"f"], EconResults.mat,type="l",lty=1:3,     
        col=1,ylab="Flow / R (US$)", xlab="f (fishing effort, unit)")
minor.tick(nx=5, ny=5, tick.ratio=0.5)
legend("bottomright",legend=c("Costs", "Revenue", paste("Profit, fMEY=",fMEY)), 
       lty=1:3,col=1,bty="n")
title(main=paste("Economics of fishery in equilibrium, at tc", tc,"y"),cex.main=0.8)
abline(h=0,lty=1,col=2)
abline(v=fMEY,lty=5,col=4)





#END 1

# keep only necessary variables in workspace
#rm(list=(ls()[ls()!="Results.mat"]))
rm(list= ls()[!(ls() %in% c('Results.mat','EconResults.mat'))])

#to pass results files (e.g. Results.mat) to another software as text (.txt) 
#         or Excel (.xls) file do:
#set working directory (Files in bottom right pane of RStudio)
#write.table(Results.mat,file = "Results.txt")
#write.table(Results.mat,file = "Results.xls")
# obs: in Excel may have to import as space delimited file

#END 2 (not run)

#LICENSE & REFERENCE

#Santos, J. 2015. FI??H IT 1.0 - Student Manual: A Training System for Aquatic Resource Managers. Septentrio Educational 2015(3).
#DOI:http://dx.doi.org/10.7557/se.2015.3

#Chapter 7a R-script YR I PR  by projection contour F tc
#DOI: http://dx.doi.org/10.7557/8.3602


#This work is licensed under a Creative Commons Attribution 4.0 International License.