#############################
# SYNPOP2, version 3.1 beta #
#############################
# WITH A SIMPLE TABLE OF INPUTS
#####
##
#
# With error terms (% C.V. of Linf and K, e.g. 10% and 3%)
# 4 cohorts
# seasonal recruitment (e.g. 90 days per year)
######## ######## ########
######## SYNPOP2 ########
######## ######## ########
# (Copyright: R. Schwamborn, 2010, 2017)
##### Simple Steps:
### A. INPUTS: input manually all parameters one by one (or keep defaults)
### B. Define Functions for growth and mortality
### C. Calculate and Plot raw growth trajectories
### D. Sample monthly (LFD histograms)
### E. Apply recruitment and gear selection
### F. Apply mortality
### G. Plot the Catch Curve`
### H. Outputs to *.csv Files (monthly histograms and Catch Curve)
##############
## A. INPUTS ##
##############
# INPUT Cohort stength
NCoh_t0 <- 2000 # N individuals in each Cohort, i.e. growth trajectories per year
# INPUT Growth Parameters
Linf_input = 10
cvLinf_input = 3 # INPUT Coef. of Variation,in % (e.g. 3 = 3%)
K_input = 0.8
cvK_input = 20 # INPUT Coef. of Variation,in % (e.g. 30 = 30%)
t0_input = -0.1
# INPUT Recruitment season
cvt0_input = 50 # INPUT Coef. of Variation,in % (e.g. 300 = 300% = +-110 days = 220 days DURATION) -> used for Recruitment variab.
sdt0 <- abs(t0_input * (cvt0_input/100)) # standard deviation of the recruitment curve, in years
(Days.season = round((2* (sdt0 <- abs(t0_input * (cvt0_input/100)))*365), 0)) # Recruitment season (days)
## [1] 36
# Days.season: Recruitment season (days) = 2 * standard deviation of the recruitment curve, in days.
# INPUT Recruitment season
# Mortality
Z_input = 1.8 # TOTAL Mortality Z
###### ###### ###
# End of INPUTS #
###### ###### ###
################################################################################
### Calculation of vectors
sdLinf <- Linf_input* (cvLinf_input/100)
Linfe <- rnorm(1, mean = Linf_input, sd = sdLinf)
sdK <- K_input * (cvK_input/100)
Ke <- rnorm(1, mean = K_input, sd = sdK)
sdt0 <- abs(t0_input * (cvt0_input/100))
t0e <- rnorm(1, mean = t0_input, sd = sdt0)
########################
# B. Define Functions #
########################
########################
## B.a Define the Mortality-at-length-function
## "L", "Z", "Nzero",Linf, K, t0 # output: N(L)
Mort.at.L.fun3 <- function(L.t, Z_inp, N_0inp, Linf_in, K_in, tzero) {
N.0 <- N_0inp
M <- Z_inp
inv.VBGF.fun.simple <- function(L) {
Linf <- Linf_in
K <- K_in
t0 <- -tzero
t = ( log(Linf/(Linf-L))*(1/K) )+t0
;t
}
age.2 <- inv.VBGF.fun.simple(L.t)
N.m <- N.0 * exp(-(age.2 * M))
; N.m
} ### EOF ### end of Mort.at.L.fun2
Mort.at.L.fun3(7, Z_input, 1000, Linf_input, K_input, 0)
## [1] 66.60745
## "L", total Mort. "Z", Nzero, Linf, K, tzero # output: N(L)
########################
## B.b Define growth function "VBG.fun.ecv"
########################
VBG.fun.ecv <- function(t, Linf,cvLinf, K, cvK, t0, cvt0) {
sdLinf <- Linf* (cvLinf/100)
Linfe <- rnorm(1, mean = Linf, sd = sdLinf)
sdK <- K * (cvK/100)
Ke <- rnorm(1, mean = K, sd = sdK)
sdt0 <- abs(t0 * (cvt0/100))
t0e <- rnorm(1, mean = t0, sd = sdt0)
L_t = (Linfe * (1 - exp(-Ke*(t-t0e))))
; L_t
} ### EOF ### end of function VBG.fun.ecv
VBG.fun.ecv(t = 8, Linf = 10, 3, 0.5, 3, 0,3)# input: t, Linf,cvLinf (%), K, cvK(%), t0, output: L(t)
## [1] 10.42419
##########################################################################
### C. Calculate and Plot raw growth trajectories
##################
## Calculations ##
##################
### C.a. Generates Input tables
#
age = (seq(0,1, length.out = NCoh_t0))
# ages of cohort 1
inp.table.1 <- data.frame(t = seq(0,1, length.out = NCoh_t0),Linf = rep(Linf_input, NCoh_t0), cvLinf = rep(cvLinf_input, NCoh_t0),K = rep(K_input, NCoh_t0), cvK = rep(cvK_input, NCoh_t0),to = rep(t0_input, NCoh_t0), cvto = rep(cvt0_input, NCoh_t0))
# View(inp.table.1)
# ages of cohort 2
inp.table.2 <- data.frame(cbind (t = seq(1,2, length.out = NCoh_t0),Linf = rep(Linf_input, NCoh_t0), cvLinf = rep(cvLinf_input, NCoh_t0),K = rep(K_input, NCoh_t0), cvK = rep(cvK_input, NCoh_t0),to = rep(t0_input, NCoh_t0), cvto = rep(cvt0_input, NCoh_t0)))
# View(inp.table.2)
# ages of cohort 3
inp.table.3 <- data.frame(cbind (t = seq(2,3, length.out = NCoh_t0),Linf = rep(Linf_input, NCoh_t0), cvLinf = rep(cvLinf_input, NCoh_t0),K = rep(K_input, NCoh_t0), cvK = rep(cvK_input, NCoh_t0),to = rep(t0_input, NCoh_t0), cvto = rep(cvt0_input, NCoh_t0)))
# View(inp.table.3)
# ages of cohort 4
inp.table.4 <- data.frame(cbind (t = seq(3,4, length.out = NCoh_t0),Linf = rep(Linf_input, NCoh_t0), cvLinf = rep(cvLinf_input, NCoh_t0),K = rep(K_input, NCoh_t0), cvK = rep(cvK_input, NCoh_t0),to = rep(t0_input, NCoh_t0), cvto = rep(cvt0_input, NCoh_t0)))
# View(inp.table.3)
#### Simulating n = 4,000 Growth trajectories with unique individual growth paramters
### 4 cohorts
# Generates a complete table (Outputs + Inputs, puts everything into a file: "outputXX2.csv")
# RUN
w <-1
outp.obj.1 <-matrix(data = 0, ncol= 1, nrow = NCoh_t0)
outp.obj.2 <-matrix(data = 0,ncol= 1, nrow = NCoh_t0)
outp.obj.3 <-matrix(data = 0,ncol= 1, nrow = NCoh_t0)
outp.obj.4 <-matrix(data = 0,ncol= 1, nrow = NCoh_t0)
###########
# C.b. LOOP # Generates Growth trajectories
###########
for(w in 1:NCoh_t0) {
#coh1
outp.obj.1[w,1] <- VBG.fun.ecv(inp.table.1[w,1], inp.table.1[w,2], inp.table.1[w,(3)], inp.table.1[w,(4)], inp.table.1[w,(5)], inp.table.1[w,(6)], inp.table.1[w,(7)])
#View(outp.obj.4)
#View(inp.table.1)
#coh2
outp.obj.2[w,1] <- VBG.fun.ecv(inp.table.2[w,1], inp.table.2[w,2], inp.table.2[w,(3)], inp.table.2[w,(4)], inp.table.2[w,(5)], inp.table.2[w,(6)], inp.table.2[w,(7)])
#coh3
outp.obj.3[w,1] <- VBG.fun.ecv(inp.table.3[w,1], inp.table.3[w,2], inp.table.3[w,(3)], inp.table.3[w,(4)], inp.table.3[w,(5)], inp.table.3[w,(6)], inp.table.3[w,(6)])
#coh4
outp.obj.4[w,1] <- VBG.fun.ecv(inp.table.4[w,1], inp.table.4[w,2], inp.table.4[w,(3)], inp.table.4[w,(4)], inp.table.4[w,(5)], inp.table.4[w,(6)], inp.table.4[w,(7)])
# write.table(outp.obj.all, file = "output_x.csv", append = TRUE, sep = ",", row.names = FALSE, col.names = FALSE)
}
###############
# end of loop #
###############
# Read full data (10,001 lines) from file into results.1:
# C.c # POOL all 4 cohorts together
inp.obj.all <- rbind(inp.table.1, inp.table.2, inp.table.3, inp.table.4)
outp.obj.all_temp <- rbind(outp.obj.1, outp.obj.2, outp.obj.3, outp.obj.4)
outp.obj.all_temp2 <- cbind(outp.obj.all_temp,inp.obj.all )
outp.obj.all <- as.data.frame(outp.obj.all_temp2)
outp.obj.all$L <- outp.obj.all[,1]
results.1 <- outp.obj.all
# View(outp.obj.all)
# results.1 <- read.csv("output_x.csv", col.names = c( "age", "Linf","cvLinf", "K","cvK","t0", "L1", "L2", "L3", "L4"))
results.1 <- subset(results.1, select = -c(outp.obj.all_temp) )
attach(results.1)
# View(results.1)
# remove temporary csv file from disk (optional):
# file.remove("output_x.csv")
#Output for later analysis
# setwd("~/to/path")
# filename =paste("SYNPOP2_Run_2") # NAME of RUN
# write.table(results.1, file = "xxx.csv")
#write.table(results.1, file = paste("filename","_Z",Z_input,"_K",K_input,
# "_CVK", cvK_input,"_CVLinf", cvLinf_input,
# "CV_t0", cvt0_input, "raw1.csv"),
# append = FALSE, sep = ",", row.names = TRUE, col.names = TRUE)
# View(results.1) # view results
##########################################
##########################################
#### C.d PLOT ALL 4 cohorts in 1 graph
## PLots the 4 cohorts without mortality N= NCoh_t0 individuals
# View(Res.all)
### Plots 4 cohorts
t.new <- rep(age,4)
L.all <- outp.obj.all$L
Res.all <- data.frame( age = t.new, L = L.all)
###
plot(Res.all, xlab = "t (years)", ylab = "L (cm)", col = "white",
main = paste("4 * ", NCoh_t0, " Growth Trajectories"))
points( Res.all, col=rgb(0, 0.1, 1, 0.17), pch=1, cex = 0.7)
legend("bottomright", legend =
paste("Linf: ",Linf_input,"cm +-",cvLinf_input,"%, K: ",
K_input," y-1 +-",cvK_input,"%, to: ",t0_input,
"y+-",cvt0_input, "% (", Days.season, "days y-1)", sep= ""),
cex = 0.54, bty = "n")
###
# INPUT Growth Parameters
Linf_input = 10
cvLinf_input = 3 # INPUT Coef. of Variation,in % (e.g. 3 = 3%)
K_input = 0.8
cvK_input = 20 # INPUT Coef. of Variation,in % (e.g. 30 = 30%)
t0_input = -0.1
# INPUT Recruitment season
cvt0_input = 50 # INPUT Coef. of Variation,in % (e.g. 300 = 300% = +-110 days = 220 days DURATION) -> used for Recruitment variab.
(Days.season = 2* (sdt0 <- abs(t0_input * (cvt0_input/100)))*365) # Recruitment season (days)
## [1] 36.5
# INPUT Recruitment season
# Mortality
Z_input = 1.8 # TOTAL Mortality Z
##############################################################
##############################################################
##############################################################
#### HISTOGRAMS For ELEFAN AND CATCH CURVE ANALYSIS #
#### 4 cohorts
###############################
#### D. Sample monthlY
#### E. Apply a gear selection and recruitment curve function
#### F. Apply mortality
#### G. Plot the Catch Curve
#### H. Outputs to *.csv files
attach(Res.all)
## The following object is masked _by_ .GlobalEnv:
##
## age
## The following object is masked from results.1:
##
## L
# View(results.1)
# View(Res.all)
length(t)
## [1] 8000
length(L)
## [1] 8000
max(Res.all$L)
## [1] 10.63087
##########################################################################
### D. Sample monthly (LFD histograms)
########
# Length-Frequency-Distribs for FISAT and Histograms
## This selects only data with L < Linf
# Reason: the inverted VBGF function is not defined when L > Linf
# windows(5,5,); plot(Res.all)
attach(Res.all)
## The following object is masked _by_ .GlobalEnv:
##
## age
## The following objects are masked from Res.all (pos = 3):
##
## age, L
## The following object is masked from results.1:
##
## L
Res.all <- subset (Res.all, L < Linf_input)
Res.all <- subset (Res.all, L > 0)
###########
##D.a Subsets by age (5 days sampled per month)
#
#jan
a0 <- subset(Res.all, age < (1/12))
a0 <- subset(a0, age > ((1/12)-(5/365)))
a0b <- subset(Res.all, age > (1/12))
length(a0$L)
## [1] 108
length(a0[[1]])
## [1] 108
#plot(a0)
#feb
a1 <- subset(a0b, age < (2/12))
a1 <- subset(a1, age > ((2/12)-(5/365)))
a1b <- subset(Res.all, age > (2/12))
length(a1[[1]])
## [1] 112
#plot(a1)
#mar
a2 <- subset(a1b, age < (3/12))
a2 <- subset(a2, age > ((3/12)-(5/365)))
a2b <- subset(a1b, age > (3/12))
#plot(a2)
#apr
a3 <- subset(a2b, age < (4/12))
a3 <- subset(a3, age > ((4/12)-(5/365)))
a3 <- subset(a3, age > ((4/12)-(5/365)))
a3b <- subset(a2b, age > (4/12))
#may
a4 <- subset(a3b, age < (5/12))
a4 <- subset(a4, age > ((5/12)-(5/365)))
a4b <- subset(a3b, age > (5/12))
#jun
a5 <- subset(a4b, age < (6/12))
a5 <- subset(a5, age > ((6/12)-(5/365)))
a5b <- subset(a4b, age > (6/12))
#jul
a6 <- subset(a5b, age < (7/12))
a6 <- subset(a6, age > ((7/12)-(5/365)))
a6b <- subset(a5b, age > (7/12))
#aug
a7 <- subset(a6b, age < (8/12))
a7 <- subset(a7, age > ((8/12)-(5/365)))
a7b <- subset(a6b, age > (8/12))
#sep
a8 <- subset(a7b, age < (9/12))
a8 <- subset(a8, age > ((9/12)-(5/365)))
a8b <- subset(a7b, age > (9/12))
#oct
a9 <- subset(a8b, age < (10/12))
a9 <- subset(a9, age > ((10/12)-(5/365)))
a9b <- subset(a8b, age > (10/12))
#plot(a9)
#nov
a10 <- subset(a9b, age < (11/12))
a10 <- subset(a10, age > ((11/12)-(5/365)))
a10b <- subset(a9b, age > (11/12))
# plot(a10)
#length(a10[[1]])
#dec
a11 <- subset(a10b, age < (12/12))
a11 <- subset(a11, age > ((12/12)-(5/365)))
length(a11[[1]])
## [1] 102
# plot(a11)
# summary(a11)
#length(a7[[1]])
br.1 <- seq(0,(Linf_input), by = 0.5) ## This selects only data with L < Linf (the inverted VBGF function is not defined when L > Linf)
#histograms by month
out11 <- hist(a11[[2]], breaks = br.1)
out10 <- hist(a10[[2]], breaks = br.1)
out9 <- hist(a9[[2]], breaks = br.1)
out8 <- hist(a8[[2]], breaks = br.1)
out7 <- hist(a7[[2]], breaks = br.1)
out6 <- hist(a6[[2]], breaks = br.1)
out5 <- hist(a5[[2]], breaks = br.1)
out4 <- hist(a4[[2]], breaks = br.1)
out3 <- hist(a3[[2]], breaks = br.1)
out2 <- hist(a2[[2]], breaks = br.1)
out1 <- hist(a1[[2]], breaks = br.1)
out0 <- hist(a0[[2]], breaks = br.1)
## D.b. "pseudo catch curve" (with zero mortality)
# Raw "pseudo catch curve" (with zero mortality) of ALL Data
outall <- hist(Res.all[[2]], breaks = br.1)
log.outall <- hist(log(Res.all[[2]]))
# "pseudo catch curve" (with zero mortality)
plot(outall$counts ~ outall$mids)
plot(log(outall$counts) ~ outall$mids)# beware of pseudo-mortality after L > 9
length(outall$mids)
## [1] 20
max(Res.all$L)
## [1] 9.999741
max(outall$mids)
## [1] 9.75
output.for.FisatIIIMOrtZero <- data.frame(Lavg = out2$mids, jan = out0$counts,
feb =out1$counts, mar = out2$counts, apr = out3$counts, may = out4$counts, jun = out5$counts, jul = out6$counts,
aug =out7$counts, sep = out8$counts, oct = out9$counts, nov = out10$counts, dec = out11$counts)
# View(output.for.FisatIIIMOrtZero)
# write.table(output.for.FisatIIIMOrtZero, file = "outputFisatIIIMOrtZeroxx.csv", append = FALSE, sep = ",", row.names = FALSE, col.names = TRUE)
#############
# Now apply Mortality (Z_input) and Recruitment to L-Freq (Histograms)
#############
######
# E. Recruitment + gear selection: logistic curve with midpont at L = 2 cm
######
recruit <- SSfpl(out2$mids, 0, 1, 2,0.5) # recruitment + gear selection: logistic curve with midpont at L = 2 cm
plot(recruit ~ out2$mids)
plot(recruit ~ out2$mids, xlim = c(0,5), main= "Recruitment & Gear Selection", ylab= "Recruitment (0 to 1)", xlab = "L (cm)")
legend("bottomright", legend = "midpont=2 scale=0.5")
######
# F. Apply Mortality (using the inverted VGBF function: "Mort.at.L.fun3")
######
Z_input # Z_input
## [1] 1.8
###
## F.a Apply Mortality only
#
Mort.at.L.fun3(7, Z_input, 1000, Linf_input, K_input, 0)
## [1] 66.60745
## "L", total Mort."Z", Nzero, Linf, K, tzero # output: N(L)
Mort.at.L.fun3(out2$mids[1], Z_input, 1, Linf_input, K_input, 0)
## [1] 0.9446271
m.t.1_vec <- rep(0,39)
for (w in 1:39) {
m.t.1_vec[w] <- round(Mort.at.L.fun3(out2$mids[w], Z_input, 1, Linf_input, K_input, 0),2)
}
M.t.F <- subset(m.t.1_vec, m.t.1_vec!= "NA" | m.t.1_vec!= 0) # Vector (without "NA"s)giving N at each length class midpoint
M.t.F <- subset(M.t.F, M.t.F != 0) # Vector (no Zeros) giving N at each length class midpoint
length(out2$mids)
## [1] 20
length(m.t.1_vec)
## [1] 39
length(M.t.F)
## [1] 18
out2$mids <- out2$mids[1:length(M.t.F)] #crop (adjust length of vector)
length(out2$mids)
## [1] 18
### PLot
# M.t.F # Vector giving Mortality (as a factor) at each length class midpoint*100
# plot(M.t.F ~ out2$mids)
# plot(log(M.t.F) ~ out2$mids)
plot(M.t.F ~ out2$mids, xlim = c(0,12), main= "Mortality only", ylab= "Mortality(%)", xlab = "L (cm)")
###
#####
## F.b. Apply Mortality and Recruitment
### Integrate recruitment and mortality into one vector (M.t.2)
#
recruit <- recruit[1:length(M.t.F)] #crop (adjust length of vector)
M.t.2 <- M.t.F*recruit
#######################################################
### G. Plot the Catch Curve
# Final catch curve #
### Plots the Recruitment, Gear Selection (G.S.) and Mortality Factor
plot(M.t.2 ~ out2$mids, xlim = c(0,10), main= "Factor, Recruitment, G.S. & Mortality", ylab= "Mortality(%)", xlab = "L (cm)" )
legend("topright",
legend = paste("Z: ", Z_input,"K: ",K_input), cex = 0.8 )
### Plots the Catch Curve
# "pseudo catch curve" (with zero mortality)
# plot(outall$counts ~ outall$mids)
# length(outall$counts)
outall$mids <- outall$mids[1:length(M.t.F)] #crop (adjust length of vector)
outall$counts <- outall$counts[1:length(M.t.F)] #crop (adjust length of vector)
M.t.CC <- M.t.2 * outall$counts
plot(M.t.CC ~ out2$mids, xlim = c(0,10), main= "Catch Curve", ylab= "Catch", xlab = "L (cm)" )
legend("topright",
legend = paste("Z: ", Z_input,"K: ",K_input), cex = 0.8 )
# write csv file on disk
# output.for.FisatIII.catch_curveA <- data.frame(N = M.t.2, L = out2$mids)
# write.table(output.for.FisatIII.catch_curveA, file = "outputforFisatIII-catch-curveAiii.csv", append = FALSE, sep = ",", row.names = FALSE, col.names = TRUE)
##########################################
#### H. OUTPUTS
### Output to *.csv Files (monthly histograms and Catch Curve)
## Ia. Save Catch Curve to disk
# write csv file on disk
# output.for.FisatIII.catch_curveA <- data.frame(N = M.t.2, L = out2$mids)
# write.table(output.for.FisatIII.catch_curveA, file = "outputforFisatIII-catch-curveAiii.csv", append = FALSE, sep = ",", row.names = FALSE, col.names = TRUE)
##########################################
## Ib. Outputs for ELEFAN: monthly histograms
# Apply Gear Selection and Mortality M.t.2 to Monthly Subsets by Length, round to zero digits
out11M <- out11
out11M$counts <- out11M$counts[1:length(M.t.2)] * M.t.2
out11M$counts <- round(out11M$counts, 3)
out10M <- out10
out10M$counts <- out10M$counts[1:length(M.t.2)] * M.t.2
out10M$counts <- round(out10M$counts, 3)
out9M <- out9
out9M$counts <- out9M$counts[1:length(M.t.2)] * M.t.2
out9M$counts <- round(out9M$counts, 3)
out8M <- out8
out8M$counts <- out8M$counts[1:length(M.t.2)] * M.t.2
out8M$counts <- round(out8M$counts, 3)
out7M <- out7
out7M$counts <- out7M$counts[1:length(M.t.2)] * M.t.2
out7M$counts <- round(out7M$counts, 3)
out6M <- out6
out6M$counts <- out6M$counts[1:length(M.t.2)] * M.t.2
out6M$counts <- round(out6M$counts, 3)
out5M <- out5
out5M$counts <- out5M$counts[1:length(M.t.2)] * M.t.2
out5M$counts <- round(out5M$counts, 3)
out4M <- out4
out4M$counts <- out4M$counts[1:length(M.t.2)] * M.t.2
out4M$counts <- round(out4M$counts, 3)
out3M <- out3
out3M$counts <- out3M$counts[1:length(M.t.2)] * M.t.2
out3M$counts <- round(out3M$counts, 3)
out2M <- out2
out2M$counts <- out2M$counts[1:length(M.t.2)] * M.t.2
out2M$counts <- round(out2M$counts, 3)
out1M <- out1
out1M$counts <- out1M$counts[1:length(M.t.2)] * M.t.2
out1M$counts <- round(out1M$counts, 3)
out0M <- out0
out0M$counts <- out0M$counts[1:length(M.t.2)] * M.t.2
out0M$counts <- round(out0M$counts, 3)
output.for.FisatIIIMort18.K09.SD08 <- data.frame(Lavg = out2M$mids, jan =out0M$counts,
feb =out1M$counts, mar = out2M$counts, apr = out3M$counts, may = out4M$counts, jun = out5M$counts, jul = out6M$counts,
aug =out7M$counts, sep = out8M$counts, oct = out9M$counts, nov = out10M$counts, dec = out11M$counts)
# View(output.for.FisatIIIMort18.K09.SD08)
# write.table(output.for.FisatIIIMort18.K09.SD08, file = "output.for.FisatIIIMort18-K09-SD08.csv", append = F, sep = ",", row.names = F, col.names = T)
###############################################################################
###############################################################################
###############################################################################