#Exercises Yield and economics pr recruit I - capture fishery
#as is Per Recruit model
#Yield: F and tc contour for a stock in equilibrium. Projection matrix method.

#Author: Jorge Santos, Tromsų, 2013, 2014

#Suitable for a general analysis of fishing strategy (intensity and pattern)
#Adjust 1ary and 2ary input according to species biology and fishery
#Choose range of ages of 1st capture. These will be combined with F range [0,2].
#Install required R-libraries when necessary
#To run model highlight the whole script and Run (or ctrl-enter)
#Aggregated output in YFtc.mat and PYtc.mat files and 2 figures (contours)

#needs
library(graphics)

#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-range, F~qf      #Age 1st capture (knife-edge) (y); relationship between fishing mortality and fishing effort
tcup=?                        #max age considered (y)
tcdown=?                      #min age considered (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
                                  #When R=1 use approximate value of cost (fit by eye); cost per unit of f doesnt make much sense.

#cost=200                       #Production (capture) cost ($/F unit)

#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 

tc.vec=c(seq(tcdown,tcup,length=20))
F.vec=c(seq(0,2,length=20))
YFtc.mat=matrix(NA,nrow=length(tc.vec),ncol=length(F.vec))
PFtc.mat=matrix(NA,nrow=length(tc.vec),ncol=length(F.vec))
Results.mat=matrix(NA,nrow=7,ncol=length(F.vec))
rownames(Results.mat)=c("F", "Survivors", "CatchN","Yield",
                        "Costs","Revenue","Profit")

for (t in 1:length(tc.vec)){
  tc=tc.vec[t]
  
  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

#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])          
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])    
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])

#pass values between matrices
YFtc.mat[t,j]=Yield
PFtc.mat[t,j]=Profit
  }
}



#PLOTS
#needs
library(graphics)

#Transpose matrices (required for contour plots in R)
YFtc.mat <- t(YFtc.mat)[,nrow(YFtc.mat):1]
PFtc.mat <- t(PFtc.mat)[,nrow(PFtc.mat):1]


# Filled contour plot of Yield
filled.contour(x = F.vec, y = tc.vec, z = YFtc.mat,color = topo.colors,
               plot.title = title(main = "Yield contour for combinations of F and tc",
               xlab = "F - fishing mortality (y-1)", ylab = "tc- age of 1st capture (y)",cex.main=0.6),
               key.title = title(main = "(g/R)",cex.main=0.6), 
               plot.axes = {axis(1); axis(2); points(10, 10)})

# Filled contour plot of Profit
filled.contour(x = F.vec, y = tc.vec, z = PFtc.mat,color = topo.colors,
               plot.title = title(main = "Profit contour for combinations of F and tc",
               xlab = "F - fishing mortality (y-1)", ylab = "tc- age of 1st capture (y)",cex.main=0.6),
               key.title = title(main = "($/R)",cex.main=0.6),
               plot.axes = { axis(1); axis(2); points(10, 10) })


#END 1

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

#to pass results files (e.g. PYtc.mat) to another software as text (.txt) 
#         or Excel (.xls) file do:
#set working directory (Files in bottom right pane of RStudio)
#
#write.table(PYtc.mat,file = "PYtc.txt")
#write.table(PYtc.mat,file = "PYtc.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.