sprat-xsa-experimental

1
# ************************************* #
2
# Author: Piatinskii M., Azov-black sea bench of VNIRO
3
# &copy 2019
4
# **************************************#
5
# make cleanup
6
rm(list = ls(all.names=TRUE))
7
# load libraries
8
library("FLCore")
9
library("FLXSA")
10
library("FLAssess")
11
library("ggplotFL")
12
# library("FLBRP") # no support for R 4.0, require compilation from *nix
13
library("FLash")
14
#install.packages("devtools")
15
#install_github("ices-tools-prod/msy")
16
library("msy")
17
library("rmarkdown")
18
library("dplyr")
19
library("magrittr")
20
library("tibble")
21
library("icesAdvice")
22

23
# ****** GLOBAL CONFIGURATIONS SECTION *******#
24

25
# stock object configurations
26
config.stock.name <- "Sprattus sprattus" # fish title
27
config.stock.desc <- "Sprat stock information in FLCore format from 2000years" # fish stock description
28
config.stock.fbar.min <- 1 # minimum fish age become fishing (Fbar min)
29
config.stock.fbar.max <- 3 # maximum fish age at fishing (Fbar max)
30
config.stock.no_harvestspawn <- TRUE # Set TRUE if there is no information about fishing before spawning
31
config.stock.no_mspawn <- TRUE # Set TRUE if no information about natural mortality before spawning
32
config.stock.no_discards <- TRUE # Set TRUE if no separated information about landing/discards exist (only catch)
33
config.stock.window.start <- 1994 # model start
34
config.stock.window.end <- 2022 # model end
35

36
# configure surveys indices 
37
config.index1.name <- "Sprat crimea cpue index" 
38
config.index1.desc <- "Fishery index information about sprat abundance at Crimea shelf from real fishing surveys"
39
config.index1.startf <- 0.3 # start month of survey in format month/12, example: 1 feb = 2/12 = 0.166
40
config.index1.endf <- 0.7 # end month of survey in format month/12
41
config.index1.type <- "number" # type of survey, allowed: number, biomass
42

43
# same with previous
44
config.index2.name <- "Sprat caucasian cpue index" 
45
config.index2.desc <- "Fishery index information about sprat abundance at Caucasian shelf from real fishing surveys"
46
config.index2.startf <- 0.3
47
config.index2.endf <- 0.7
48
config.index2.type <- "number"
49

50
# pre-analysis diagnostic procedures. Useful to choice the model parametrization, convergence, etc
51
config.iter.fse = TRUE # Fse choice procedure
52
config.iter.retro = TRUE # retrospective procedure
53
config.iter.fseage = TRUE # Fse at age impact procedure
54

55
# configure retrospective depth
56
config.retro.horizon <- 5
57

58
# Shrinkage standard error max value. Default value Fse = 0.5
59
# Higher shrinkage se leads to high unrestraint of estimating procedure. 
60
config.tuning.xsa <- FLXSA.control(
61
  maxit = 100, # maximum covergence iterations
62
  fse = 1.5, # shrinkage s.e value
63
  shk.n = T,
64
  shk.f = T,
65
  shk.yrs = 4,
66
  tsrange = 15, # weighting time series length (from latest year)
67
  tspower = 3, # weighting procedure sequence (1 = linear, 2 = square, 3 = cubic)
68
  qage = 2, # age after catchability not depends on age
69
)
70

71
# sr model configuration
72
config.srmodel.eqsim.use = TRUE # should we use eqsim package insted of FLSR default?
73
config.srmodel.use_geomean = FALSE # set FALSE to use EQSIM SR model results to forecasting
74

75
# eqsim model configurations (if config.srmodel.eqsim.use = TRUE)
76
config.srmodel.eqsim.models = c("Segreg") # types of model to fit and compare
77
config.srmodel.eqsim.simulations = 10 # number of simulations
78
config.srmodel.eqsim.shift_years <- 1994:2005 # remove years from sr-model fitting
79

80
# flsr model configuration (if config.srmodel.eqsim.use = FALSE)
81
config.srmodel.type <- segregAR1() # SR model type: ricker(), bevholt(), geomean(), segreg(), shepherd(), cushing(), AR in model name = AutoCorrelation fix
82
config.srmodel.optimization = FALSE # Use optimization procedure for SR model fitting?
83

84
# configure MSY reference point calculation: Blim, Bpa, Fmsy
85
config.brp.blim_multiplier <- 0.3 # multiplier of SSBmax to get Blim ref.point, allowed: 0.2 - 0.4
86
config.brp.trim_uncertainity <- TRUE # trim uncerteinity extreme deviations? Optimization procedure for Bpa estimate
87
config.brp.current_range <- 2005:2022 # shift previous population states from BRP calculations
88

89
# Custom forecast scenarious. By default performed: F_SQ, F_01, F_msy scenarious.
90
# There you can define f(y+1), f(y+2), f(y+3) targets for 2 scenarious
91
# precautionary scenario: Fpa = (B-Blim/Bmsy-Blim) * Fmsy = 0.51
92
config.forecast.scenario1.f_year1 <- 0.51
93
config.forecast.scenario1.f_year2 <- 0.51
94
config.forecast.scenario1.f_year3 <- 0.51
95
config.forecast.scenario1.f_year4 <- 0.51
96
config.forecast.scenario1.f_year5 <- 0.51
97

98
# Scenario MSY at M = 0,95
99
# by old STECF (2010) results
100
config.forecast.scenario2.f_year1 <- 0.91
101
config.forecast.scenario2.f_year2 <- 0.91
102
config.forecast.scenario2.f_year3 <- 0.91
103
config.forecast.scenario2.f_year4 <- 0.91
104
config.forecast.scenario2.f_year5 <- 0.91
105

106
# additional scenario for manual Fmsy and F_01 estimates
107
config.forecast.manual_F01 <- 0.605
108
config.forecast.manual_FMSY <- 0.64
109

110
# custom calculations. Calculate partial recruitment for Y/R analysis
111
config.custom.partial.years <- 2004:2020
112

113
# custom calculations. Calculate F01 reference point from Y/R analysis
114
config.custom.f01.years <- 2004:2020 # years used to calculate mean(wt), mean(Fi/Fmax), mean(M)
115
config.custom.f01.plusgroup = TRUE # is last group is plus group?
116
config.custom.f01.oldest = 5 # oldest fish found in waters (not in stock, that means "even somewhen exist in sea")
117
# ******** GLOBAL CONFIGURATION DONE!!! **** #
118

119
# ********************************************************* #
120
# Do not change enything below if you have no idea what is going on
121
# ********************************************************* #
122

123
# simple write functions
124
savePic <- function(path, imgObj, width=1900, height=1600, res=200) {
125
  png(path, width=width, height=height, res=res, type="cairo")
126
  print(plot(imgObj) + theme_bw())
127
  dev.off()
128
}
129
saveText <- function(path, object) {
130
  con <- file(path, encoding = "UTF-8")
131
  sink(con)
132
  print(object)
133
  sink()
134
  close(con)
135
}
136

137
# get mbar from stock object
138
mbar <- function(stock) {
139
  stock.min <- as.numeric(stock@range["min"])
140
  stock.max <- as.numeric(stock@range["max"])
141
  
142
  bar.min <- as.numeric(stock@range["minfbar"])
143
  bar.max <- as.numeric(stock@range["maxfbar"])
144
  
145
  diff.min <- 1 + bar.min - stock.min
146
  diff.max <- 1 + bar.max - bar.min
147
  if (diff.max > stock.max+1)
148
    diff.max <- stock.max+1
149
  
150
  vec <- (m(stock)[diff.min:diff.max,])
151
  bar <- apply(vec, 2, mean)
152
  return(bar)
153
}
154

155

156

157
# read input data - stock
158
tmp.stock <- read.csv("input/stock.csv", sep=",")
159
data.stock <- as.FLStock(tmp.stock)
160
# read input data - index
161
tmp.idx1.idx <- read.csv("input/survey1/index.csv", sep=",")
162
tmp.idx1.eff <- read.csv("input/survey1/effort.csv", sep=",")
163
tmp.idx2.idx <- read.csv("input/survey2/index.csv", sep=",")
164
tmp.idx2.eff <- read.csv("input/survey2/effort.csv", sep=",")
165

166
data.index1 <- FLIndex(index=as.FLQuant(tmp.idx1.idx), effort=as.FLQuant(tmp.idx1.eff), catch.n=as.FLQuant(tmp.idx1.idx))
167
data.index2 <- FLIndex(index=as.FLQuant(tmp.idx2.idx), effort=as.FLQuant(tmp.idx2.eff), catch.n=as.FLQuant(tmp.idx2.idx))
168

169
remove(tmp.stock, tmp.idx1.idx, tmp.idx1.eff, tmp.idx2.idx, tmp.idx2.eff)
170
print("[XSA toolkit]: Load data: OK")
171
# Set name,desc from cfgs
172
name(data.stock) <- config.stock.name
173
desc(data.stock) <- config.stock.desc
174
# set fbar ranges
175
range(data.stock)[c("minfbar", "maxfbar")] <- c(config.stock.fbar.min, config.stock.fbar.max)
176

177
# shift data starting from 2000
178
data.stock <- window(data.stock, start = config.stock.window.start, end = config.stock.window.end)
179

180
# set harvest.spwn and m.spwn = 0 if no data exist (or XSA throw error)
181
if (config.stock.no_harvestspawn) {
182
  harvest.spwn(data.stock) <- 0
183
}
184
if (config.stock.no_mspawn) {
185
  m.spwn(data.stock) <- 0
186
}
187
# if no info about discards/landings separatly than set discards = 0, landings = catch
188
if (config.stock.no_discards) {
189
  discards(data.stock) <- 0
190
  discards.n(data.stock) <- 0
191
  discards.wt(data.stock) <- catch.wt(data.stock)
192
  landings(data.stock) <- catch(data.stock)
193
  landings.n(data.stock) <- catch.n(data.stock)
194
  landings.wt(data.stock) <- catch.wt(data.stock)
195
}
196
stock.wt(data.stock) <- catch.wt(data.stock)
197

198
# set cfgs for index1 and index2
199
type(data.index1) <- config.index1.type
200
name(data.index1) <- config.index1.name
201
desc(data.index1) <- config.index1.desc
202
range(data.index1)[c("startf", "endf")] <- c(config.index1.startf, config.index1.endf)
203

204
type(data.index2) <- config.index2.type
205
name(data.index2) <- config.index2.name
206
desc(data.index2) <- config.index2.desc
207
range(data.index2)[c("startf", "endf")] <- c(config.index2.startf, config.index2.endf)
208
print("[XSA toolkit]: FLStock and FLIndex object configured successful")
209

210
data.indices <- FLIndices(data.index1, data.index2)
211

212
# make Fse comparision test
213
if (config.iter.fse) {
214
  print("[XSA toolkit]: Fse impact factor estimate")
215
  steps <- seq(0.5, 2.5, by = 0.5)
216
  fseIters <- c()
217
  for (i in 1:length(steps)) {
218
    ctl <- config.tuning.xsa
219
    ctl@fse <- steps[i]
220
    res <- FLXSA(data.stock, data.indices, ctl) + data.stock
221
    fseIters <- c(fseIters, res)
222
  }
223
  
224
  result.fse.summary <- FLStocks(fseIters)
225
  names(result.fse.summary) <- steps
226
  # save output to pic
227
  savePic("output/1_FseChoose/summary_all.png", result.fse.summary)
228
  # get latest 5 year window
229
  lastYear <- result.fse.summary[[1]]@range["maxyear"]
230
  tmpObj <- window(result.fse.summary, lastYear-config.retro.horizon, lastYear)
231
  savePic("output/1_FseChoose/summary_short.png", tmpObj)
232
  remove(lastYear, tmpObj)
233
  # save text output - ssb
234
  con <- file("output/1_FseChoose/ssb.txt", encoding = "UTF-8")
235
  sink(con)
236
  
237
  for (i in 1:length(steps)) {
238
    print(sprintf("============= SSB at Fse = %s", steps[i]))
239
    print(ssb(result.fse.summary[[i]]))
240
    cat("\r\n")
241
  }
242
  
243
  sink()
244
  close(con)
245
  # save text output - rec
246
  con <- file("output/1_FseChoose/rec.txt", encoding = "UTF-8")
247
  sink(con)
248
  
249
  for (i in 1:length(steps)) {
250
    print(sprintf("============= Recruitment at Fse = %s", steps[i]))
251
    print(rec(result.fse.summary[[i]]))
252
    cat("\r\n")
253
  }
254
  
255
  sink()
256
  close(con)
257
  # save text output - fbar
258
  con <- file("output/1_FseChoose/fbar.txt", encoding = "UTF-8")
259
  sink(con)
260
  
261
  for (i in 1:length(steps)) {
262
    print(sprintf("============= Fbar at Fse = %s", steps[i]))
263
    print(fbar(result.fse.summary[[i]]))
264
    cat("\r\n")
265
  }
266
  
267
  sink()
268
  close(con)
269
  
270
  remove(steps, fseIters, ctl, res, con)
271
}
272

273
if (config.iter.retro) {
274
  print("[XSA toolkit]: Retrospective analysis")
275
  # make retrospective analysis
276
  result.retro.summary <- list()
277
  endYear <- data.stock@range["maxyear"]
278
  startYear <- endYear-config.retro.horizon
279
  years <- startYear:endYear
280
  frange <- seq(0.1, 2.5, by=0.1)
281
  for (i in frange) {
282
    res <- tapply(years, 1:length(years), function(x) {
283
      ctl <- config.tuning.xsa
284
      ctl@fse <- i
285
      window (data.stock, end=x) + FLXSA(window(data.stock, end=x), data.indices, ctl)
286
    })
287
    res <- FLStocks(res)
288
    res@names <- ac(c(years))
289
    result.retro.summary[[ac(i)]] <- res
290
    
291
    # save output as img
292
    savePic(sprintf("output/2_RetroFse/%s.png", i), res)
293
    remove(res)
294
  }
295
}
296

297
print("[XSA toolkit]: Retrospective Mohn's rho tests")
298
# compare XSA estimates in terminal year by -5 year retrospective
299
lastYear <- as.numeric(range(data.stock)["maxyear"])
300
ssb.base <- as.vector(ssb(window(data.stock+FLXSA(data.stock, data.indices, config.tuning.xsa), start=lastYear-5)))
301
f.base <- as.vector(fbar(window(data.stock+FLXSA(data.stock, data.indices, config.tuning.xsa), start=lastYear-5)))
302
ssb.retro <- list()
303
f.retro <- list()
304
for (i in 1:config.retro.horizon) { # prepare retrospective matrix for mohn's rho test
305
  # make retro
306
  tmp.xsa <- FLXSA(window(data.stock, end=(lastYear-i)), window(data.indices, end=(lastYear-i)), config.tuning.xsa)
307
  tmp.xsa.st <- data.stock + tmp.xsa
308
  
309
  ssb.t <- as.vector(ssb(window(tmp.xsa.st, start=lastYear-config.retro.horizon)))
310
  f.t <- as.vector(fbar(window(tmp.xsa.st, start=lastYear-config.retro.horizon)))
311
  
312
  
313
  ssb.retro[[i]] <- ssb.t
314
  f.retro[[i]] <- f.t
315
}
316

317
result.retro.diagnostic.ssb <- data.frame(ssb.base, ssb.retro[[1]], ssb.retro[[2]], ssb.retro[[3]], ssb.retro[[4]], ssb.retro[[5]])
318
rownames(result.retro.diagnostic.ssb) <- (lastYear-config.retro.horizon):lastYear
319
result.retro.diagnostic.f <- data.frame(f.base, f.retro[[1]], f.retro[[2]], f.retro[[3]], f.retro[[4]], f.retro[[5]])
320
rownames(result.retro.diagnostic.f) <- (lastYear-config.retro.horizon):lastYear
321
names(result.retro.diagnostic.ssb) <- c("ssb.base", "-1", "-2", "-3", "-4", "-5")
322
names(result.retro.diagnostic.f) <- c("f.base", "-1", "-2", "-3", "-4", "-5")
323

324
# replace NaN's (x/0) to NA (not available)
325
is.nan.data.frame <- function(x) {
326
  return (do.call(cbind, lapply(x, is.nan)))
327
}
328

329
result.retro.diagnostic.f[is.nan(result.retro.diagnostic.f)] <- NA
330
result.retro.diagnostic.ssb[is.nan(result.retro.diagnostic.ssb)] <- NA
331

332
result.retro.mohnrho.ssb <- mohn(result.retro.diagnostic.ssb, peels = config.retro.horizon)
333
result.retro.mohnrho.fbar <- mohn(result.retro.diagnostic.f, peels = config.retro.horizon)
334
remove(ssb.base, f.base, ssb.retro, f.retro)
335

336
print(paste0("Rho ssb=", round(result.retro.mohnrho.ssb, 3)))
337
print(paste0("Rho F=", round(result.retro.mohnrho.fbar, 3)))
338

339
# calculate Fse impact to different age groups
340
if (config.iter.fseage) {
341
  print("[XSA toolkit]: Fse impact to F-at-age estimates")
342
  steps <- seq(0.5, 2.5, by = 0.5)
343
  f.res <- propagate(harvest(data.stock), length(steps))
344
  for (i in 1:length(steps)) {
345
    raw.control <- config.tuning.xsa
346
    raw.control@fse <- steps[i]
347
    iter(f.res, i) <- harvest(FLXSA(data.stock, data.indices, raw.control))
348
  }
349
  endYear <- data.stock@range["maxyear"]
350
  startYear <- endYear - 9 
351
  tmpObj <- xyplot(data ~ year | age, groups = iter, data = f.res, ylab="F", type = "l", xlim = c(startYear:endYear))
352
  savePic("output/3_FseImpactFAtAge/summary_all.png", tmpObj)
353
  remove(startYear, endYear, steps, f.res, raw.control, tmpObj)
354
}
355

356
# save survey regression age-vs-age test 
357
suppressWarnings(savePic("output/4_SurveysRegressions/summary_all.png", data.indices))
358

359
suppressWarnings(savePic("output/4_SurveysRegressions/survey_1.png", data.indices[[1]]))
360
suppressWarnings(savePic("output/4_SurveysRegressions/survey_2.png", data.indices[[2]]))
361

362
print("[XSA toolkit]: Performing XSA analysis...")
363
# fit xsa
364
result.xsa.fit <- FLXSA(data.stock, data.indices, config.tuning.xsa)
365
# build final stock + xsa estimates
366
result.xsa.stock <- data.stock + result.xsa.fit
367

368
# estimate uncertainty right now!
369
# get f vals from xsa results
370
tmp.fvals <- as.vector(fbar(result.xsa.stock))
371
# get whiskers statistics for distribution
372
tmp.wstats <- boxplot.stats(tmp.fvals)$stats
373
# find outliers and remove from fvals
374
tmp.fvals <- tmp.fvals[tmp.fvals > tmp.wstats[[1]]]
375
tmp.fvals <- tmp.fvals[tmp.fvals < tmp.wstats[[5]]]
376
# get fixed uncerteinity estimation
377
result.xsa.sigma <- sd(tmp.fvals)
378
remove(tmp.fvals, tmp.wstats)
379

380
# save results 
381
endYear <- result.xsa.stock@range["maxyear"]
382
F_SQ = mean(fbar(window(result.xsa.stock, (endYear-2), endYear)))
383
savePic("output/5_XsaEstimates/summary_all.png", result.xsa.stock)
384
savePic("output/5_XsaEstimates/stock_numbers_at_age.png", stock.n(result.xsa.stock))
385
savePic("output/5_XsaEstimates/f_at_age.png", harvest(result.xsa.stock))
386

387
saveText("output/5_XsaEstimates/ssb.txt", ssb(result.xsa.stock))
388
saveText("output/5_XsaEstimates/rec.txt", rec(result.xsa.stock))
389
saveText("output/5_XsaEstimates/fbar.txt", fbar(result.xsa.stock))
390
saveText("output/5_XsaEstimates/f.txt", harvest(result.xsa.stock))
391
saveText("output/5_XsaEstimates/stock.n.txt", stock.n(result.xsa.stock))
392
saveText("output/5_XsaEstimates/diagnostics.txt", diagnostics(result.xsa.fit))
393

394
# investigate uncertainty factor (build conf.intervals for plots)
395
pics <- FLQuants(
396
  Rec = rnorm(1000, rec(result.xsa.stock), result.xsa.sigma * rec(result.xsa.stock)),
397
  SSB = rnorm(1000, ssb(result.xsa.stock), result.xsa.sigma * ssb(result.xsa.stock)),
398
  F = rnorm(1000, fbar(result.xsa.stock), result.xsa.sigma * fbar(result.xsa.stock)),
399
  Catch = catch(result.xsa.stock)
400
)
401

402
savePic("output/5_XsaEstimates/summary_uncertainty.png", pics)
403
remove(endYear, pics)
404

405
# perform residual diagnostics
406
pic <- bubbles(age~year, data=result.xsa.fit@index.res, warning=FALSE)
407
suppressWarnings(savePic("output/6_Residuals/summary_all.png", pic))
408
remove(pic)
409

410
pic <- bubbles(age~year, data=result.xsa.fit@index.res[1], warning=FALSE)
411
suppressWarnings(savePic("output/6_Residuals/summary_survey1.png", pic))
412
remove(pic)
413

414
pic <- bubbles(age~year, data=result.xsa.fit@index.res[2], warning=FALSE)
415
suppressWarnings(savePic("output/6_Residuals/summary_survey2.png", pic))
416
remove(pic)
417

418
pic <- xyplot(
419
  data ~ year | ac(age) + qname, 
420
  data=index.res(result.xsa.fit), 
421
  panel = function(x,y, ...){
422
    panel.xyplot(x, y, ...)
423
    panel.loess(x,y, ...)
424
    panel.abline(h=0, col="gray", lty=2)
425
  }, 
426
  main="Surveys log catchability residuals at age")
427
suppressWarnings(savePic("output/6_Residuals/catch_at_age_summary.png", pic))
428
remove(pic)
429

430
con <- file("output/6_Residuals/residuals.txt", encoding = "UTF-8")
431
sink(con)
432
print(result.xsa.fit@index.res)
433
cat("\r\n")
434
print("===== Summary statistics")
435
print(summary(result.xsa.fit@index.res))
436
sink()
437
close(con)
438

439
# analyse recruits-per-ssb unit 
440
result.sr.recssb <- rec(result.xsa.stock) / ssb(result.xsa.stock)
441
tmp <- as.data.frame(log(result.sr.recssb))
442
png("output/7_SR_BRP/rec_per_ssb.png", width=1900, height=1600, res=200, type="cairo")
443
print(plot(tmp$year, tmp$data, type="b", xlab="Year", ylab = "Ln(Recruit/SSB)"))
444
dev.off()
445
remove(tmp)
446

447

448
# calculate stock-recruitment relation
449
# be careful to use bevholt or ricker approximation: there is different approach between eqsim and flsr packages!!!
450
# Bevholt (eqsim): rec ~ log(a * ssb / (1 + b*ssb))
451
# Bevholt (flsr): rec ~ (a * ssb / (b+ssb))
452
# ricker (eqsim): log(a) + log(b) - b*ssb
453
# ricker (flsr): a * ssb * exp(-b * ssb)
454
# segreg(flsr,eqsim): ifelse(ssb < b, ssb*b, a*b)
455
if (config.srmodel.eqsim.use) {
456
  # build SR model based on EQSIM package and MSY approach
457
  print("[XSA toolkit]: Built SR model by EQSIM")
458
  
459
  result.eqsim.srmodel <- eqsr_fit(result.xsa.stock, 
460
                                   nsamp = config.srmodel.eqsim.simulations, 
461
                                   models = config.srmodel.eqsim.models,
462
                                   remove.years = config.srmodel.eqsim.shift_years
463
  )
464
  
465
  savePic("output/7_SR_BRP/eqsim_model.png", eqsr_plot(result.eqsim.srmodel, ggPlot = TRUE))
466
  saveText("output/7_SR_BRP/eqsim_model.txt", result.eqsim.srmodel$sr.det)
467
} else {
468
  # built Stock-Recruitment model based on FLSR and fmle()
469
  print("[XSA toolkit]: Built SR model by fmle()")
470
  result.sr.stock <- as.FLSR(window(result.xsa.stock, 
471
                                    start = config.stock.window.start, 
472
                                    end = config.stock.window.end))
473
  model(result.sr.stock) <- config.srmodel.type
474
  if (config.srmodel.optimization) {
475
    print("[XSA toolkit]: Performing SANN (Nelder-Mead, quasi-Newton, etc) optimization SR model. This can take a lot of time ...")
476
    result.sr.model <- fmle(result.sr.stock, method="L-BFGS-B", control=list(trace=0))
477
  } else {
478
    result.sr.model <- fmle(result.sr.stock, control=list(trace=0))
479
  }
480
  # save output
481
  png("output/7_SR_BRP/flsr_model.png", width=1900, height=1600, res=200, type="cairo")
482
  suppressWarnings(suppressMessages(plot(result.sr.model)))
483
  dev.off()
484
  
485
  png("output/7_SR_BRP/flsr_profile.png", width=1900, height=1600, res=200, type="cairo")
486
  suppressWarnings(profile(result.sr.model, plot=TRUE, trace = FALSE))
487
  dev.off()
488
  
489
  saveText("output/7_SR_BRP/flsr_model.txt", summary(result.sr.model)) 
490
}
491

492
print("[XSA toolkit]: Calculate BRP")
493

494
# calculate Bvirgin in assumtion of SSBmax + Ci
495
#startYear <- as.numeric(data.stock@range["minyear"])
496
#terminalYear <- as.numeric(data.stock@range["maxyear"])
497
startYear <- min(config.brp.current_range)
498
terminalYear <- max(config.brp.current_range)
499
df <- data.frame(year = startYear:terminalYear, ssb = as.vector(ssb(window(result.xsa.stock, start = startYear, end = terminalYear))), catch = as.vector(catch(window(result.xsa.stock, start = startYear, end = terminalYear))))
500
tmp.maxssb <- df[df$ssb == max(df$ssb),]
501
tmp.biomass <- tmp.maxssb$ssb + tmp.maxssb$catch
502

503
# process manual BRP points
504
result.brp.blim <- tmp.biomass * config.brp.blim_multiplier
505
# get Bpa based on Blim and uncerteinity sd
506
result.brp.bpa <- result.brp.blim * exp(qnorm(0.95) * result.xsa.sigma)
507

508
# estimate BRP from EQSIM
509
lastYear <- as.numeric(result.xsa.stock@range["maxyear"])
510
if (config.srmodel.eqsim.use) {
511
  result.eqsim.brp <- eqsim_run(result.eqsim.srmodel, 
512
                               bio.years = c(lastYear-4, lastYear),
513
                               bio.const = TRUE,
514
                               extreme.trim = c(0.1, 0.9), 
515
                               Bpa = result.brp.bpa,
516
                               Blim = result.brp.blim) 
517
  result.brpsr.blim <- as.numeric(result.eqsim.srmodel$sr.det["b"])
518
  result.brpsr.bpa <- as.numeric(result.eqsim.srmodel$sr.det["b"]*exp(qnorm(0.95)*result.xsa.sigma))
519
}
520

521
if (config.srmodel.use_geomean == TRUE) {
522
  # create geomean SR model by latest 5 years geomean recruitment numbers
523
  tmp.rec.mean <- exp(mean(log(window(rec(result.xsa.stock), lastYear-5, lastYear-1))))
524
  # build geomean SR model
525
  result.sr.model <- FLSR(model = "geomean", params=FLPar(tmp.rec.mean))
526
  remove(tmp.rec.mean, tmp.maxssb, df, tmp.biomass)
527
} else {
528
  # todo: process model
529
  result.sr.model <- FLSR(model = "segreg")
530
  result.sr.model@params["a"] <- as.numeric(result.eqsim.srmodel$sr.det["a"])
531
  result.sr.model@params["b"] <- as.numeric(result.eqsim.srmodel$sr.det["b"])
532
}
533

534
# plot SSB graphic with reference points Blim, Bpa
535
minYear <- as.numeric(result.xsa.stock@range["minyear"])
536
png("output/7_SR_BRP/ssb_blim_bpa_manual.png", width=1900, height=1600, res=200, type="cairo")
537
print(
538
plot(ssb(result.xsa.stock)) + 
539
  geom_hline(aes(yintercept=result.brp.blim), linetype=2, col="red") + 
540
  annotate(geom = "text", x = minYear+1, y = result.brp.blim+result.brp.blim*0.05, label = sprintf("Blim=%s", round(result.brp.blim, 2)), color = "red") +
541
  geom_hline(aes(yintercept=result.brp.bpa), linetype=2, col="blue") +
542
  annotate(geom = "text", x = minYear+1, y = result.brp.bpa + result.brp.bpa*0.05, label=sprintf("Bpa=%s", round(result.brp.bpa, 2)), color="blue") + 
543
  theme_bw()
544
)
545
      
546
dev.off()
547
remove(minYear)
548

549
# plot SSB with reference points by SR-model
550
# @todo
551
minYear <- as.numeric(result.xsa.stock@range["minyear"])
552
png("output/7_SR_BRP/ssb_blim_bpa_srmodel.png", width=800, height=600, res=150, type="cairo")
553
print(
554
  plot(ssb(result.xsa.stock)) + 
555
    geom_hline(aes(yintercept=result.brpsr.blim), linetype=2, col="red") + 
556
    annotate(geom = "text", x = minYear+2, y = result.brpsr.blim*1.2, label = sprintf("Blim=%s", round(result.brpsr.blim, 2)), color = "red") +
557
    geom_hline(aes(yintercept=result.brpsr.bpa), linetype=2, col="blue") +
558
    annotate(geom = "text", x = minYear+2, y =result.brpsr.bpa*1.2, label=sprintf("Bpa=%s", round(result.brpsr.bpa, 2)), color="blue") + 
559
    theme_bw()
560
)
561

562
dev.off()
563
remove(minYear)
564

565
# plot Fbar with F0.1 reference point
566
minYear <- as.numeric(result.xsa.stock@range["minyear"])
567
png("output/7_SR_BRP/fbar_f01_manual.png", width=1900, height=1600, res=200, type="cairo")
568
print(plot(fbar(result.xsa.stock)) + 
569
        geom_hline(aes(yintercept=config.forecast.manual_F01), linetype=2, col="red") + 
570
        annotate(geom="text", x = minYear, y = config.forecast.manual_F01+0.1, label = sprintf("F0.1=%s", round(config.forecast.manual_F01,3)), color="red") + 
571
        theme_bw()
572
      )
573

574
dev.off()
575
remove(minYear)
576

577
# plot E = 0.4 BRP reference point (ICES advice)
578
# E = F/(F+M)
579
# E = exploitation rate (*100 = percentage of usage)
580
result.er.table <- fbar(result.xsa.stock) / (fbar(result.xsa.stock) + mbar(result.xsa.stock))
581
png("output/7_SR_BRP/e_rate04.png", width=1900, height=1600, res=200, type="cairo")
582
print(plot(result.er.table) + ylab(label = "E = F/Z") + 
583
  geom_hline(yintercept = 0.4, color = "red", linetype="dashed") + 
584
  annotate(geom = "text", x = as.numeric(result.xsa.stock@range["minyear"]), y = 0.41, label="E = 0.4", color = "red") + 
585
  theme_bw())
586
dev.off()
587

588
# Prepare forecast procedure
589
print("[XSA toolkit]: Forecast procedure with different scenarious")
590
data.stf <- stf(result.xsa.stock, nyears=5)
591
lastYear <- as.numeric(result.xsa.stock@range["maxyear"])
592

593
# perform 3-year short-term forecasts by F_SQ
594
print("[XSA toolkit]: Short-term forecast by F_SQ")
595
target.tmp.F_SQ <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=F_SQ)
596
target.forecast.F_SQ <- fwdControl(target.tmp.F_SQ)
597
result.forecast.F_SQ <- fwd(data.stf, sr=result.sr.model, control=target.forecast.F_SQ)
598
savePic("output/8_Forecasts/F_SQ/summary_all.png", result.forecast.F_SQ)
599
savePic("output/8_Forecasts/F_SQ/summary_short.png", window(result.forecast.F_SQ, lastYear, lastYear+5))
600

601
saveText("output/8_Forecasts/F_SQ/ssb.txt", ssb(result.forecast.F_SQ))
602
saveText("output/8_Forecasts/F_SQ/rec.txt", rec(result.forecast.F_SQ))
603
saveText("output/8_Forecasts/F_SQ/fbar.txt", fbar(result.forecast.F_SQ))
604
saveText("output/8_Forecasts/F_SQ/f.txt", harvest(result.forecast.F_SQ))
605
saveText("output/8_Forecasts/F_SQ/stock.n.txt", stock.n(result.forecast.F_SQ))
606
saveText("output/8_Forecasts/F_SQ/ssb.txt", ssb(result.forecast.F_SQ))
607
saveText("output/8_Forecasts/F_SQ/catch.txt", catch(result.forecast.F_SQ))
608

609
remove(target.tmp.F_SQ)
610

611
# perform 3-year short-term forecasts by manual scenarius F_msy
612
print("[XSA toolkit]: Manual short-term forecast at F_MSY")
613
target.tmp.F_01manual <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=config.forecast.manual_F01)
614
target.forecast.F_01manual <- fwdControl(target.tmp.F_01manual)
615
result.forecast.F_01manual <- fwd(data.stf, sr=result.sr.model, control=target.forecast.F_01manual)
616

617
savePic("output/8_Forecasts/F_01manual/summary_all.png", result.forecast.F_01manual)
618
savePic("output/8_Forecasts/F_01manual/summary_short.png", window(result.forecast.F_01manual, lastYear, lastYear+5))
619

620
saveText("output/8_Forecasts/F_01manual/ssb.txt", ssb(result.forecast.F_01manual))
621
saveText("output/8_Forecasts/F_01manual/rec.txt", rec(result.forecast.F_01manual))
622
saveText("output/8_Forecasts/F_01manual/fbar.txt", fbar(result.forecast.F_01manual))
623
saveText("output/8_Forecasts/F_01manual/f.txt", harvest(result.forecast.F_01manual))
624
saveText("output/8_Forecasts/F_01manual/stock.n.txt", stock.n(result.forecast.F_01manual))
625
saveText("output/8_Forecasts/F_01manual/ssb.txt", ssb(result.forecast.F_01manual))
626
saveText("output/8_Forecasts/F_01manual/catch.txt", catch(result.forecast.F_01manual))
627

628
# perform 3-year short-term forecasts by manual scenarius F_MSY=config
629
print("[XSA toolkit]: Manual short-term forecast at F_MSY")
630
target.tmp.F_MSY <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=config.forecast.manual_FMSY)
631
target.forecast.F_MSY <- fwdControl(target.tmp.F_MSY)
632
result.forecast.F_MSY <- fwd(data.stf, sr=result.sr.model, control=target.forecast.F_MSY)
633

634
savePic("output/8_Forecasts/F_MSY/summary_all.png", result.forecast.F_MSY)
635
savePic("output/8_Forecasts/F_MSY/summary_short.png", window(result.forecast.F_MSY, lastYear, lastYear+5))
636

637
saveText("output/8_Forecasts/F_MSY/ssb.txt", ssb(result.forecast.F_MSY))
638
saveText("output/8_Forecasts/F_MSY/rec.txt", rec(result.forecast.F_MSY))
639
saveText("output/8_Forecasts/F_MSY/fbar.txt", fbar(result.forecast.F_MSY))
640
saveText("output/8_Forecasts/F_MSY/f.txt", harvest(result.forecast.F_MSY))
641
saveText("output/8_Forecasts/F_MSY/stock.n.txt", stock.n(result.forecast.F_MSY))
642
saveText("output/8_Forecasts/F_MSY/ssb.txt", ssb(result.forecast.F_MSY))
643
saveText("output/8_Forecasts/F_MSY/catch.txt", catch(result.forecast.F_MSY))
644

645

646
# short-term forecast at scenario1
647
print("[XSA toolkit]: Short-term forecast by Scenario1")
648
target.tmp.F_scen1 <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=c(
649
  config.forecast.scenario1.f_year1,
650
  config.forecast.scenario1.f_year2,
651
  config.forecast.scenario1.f_year3,
652
  config.forecast.scenario1.f_year4,
653
  config.forecast.scenario1.f_year5
654
))
655
target.forecast.Fsc1 <- fwdControl(target.tmp.F_scen1)
656
result.forecast.F_scenario1 <- fwd(data.stf, sr=result.sr.model, control=target.forecast.Fsc1)
657

658
savePic("output/8_Forecasts/F_scenario1/summary_all.png", result.forecast.F_scenario1)
659
savePic("output/8_Forecasts/F_scenario1/summary_short.png", window(result.forecast.F_scenario1, lastYear, lastYear+5))
660

661
saveText("output/8_Forecasts/F_scenario1/ssb.txt", ssb(result.forecast.F_scenario1))
662
saveText("output/8_Forecasts/F_scenario1/rec.txt", rec(result.forecast.F_scenario1))
663
saveText("output/8_Forecasts/F_scenario1/fbar.txt", fbar(result.forecast.F_scenario1))
664
saveText("output/8_Forecasts/F_scenario1/f.txt", harvest(result.forecast.F_scenario1))
665
saveText("output/8_Forecasts/F_scenario1/stock.n.txt", stock.n(result.forecast.F_scenario1))
666
saveText("output/8_Forecasts/F_scenario1/ssb.txt", ssb(result.forecast.F_scenario1))
667
saveText("output/8_Forecasts/F_scenario1/catch.txt", catch(result.forecast.F_scenario1))
668

669
remove(target.tmp.F_scen1)
670

671
# st forecast scenario2
672
print("[XSA toolkit]: Short-term forecast by Scenario2")
673
target.tmp.F_scen2 <- data.frame(year=(lastYear+1):(lastYear+5), quantity="f", val=c(
674
  config.forecast.scenario2.f_year1,
675
  config.forecast.scenario2.f_year2,
676
  config.forecast.scenario2.f_year3,
677
  config.forecast.scenario2.f_year4,
678
  config.forecast.scenario2.f_year5
679
))
680
target.forecast.Fsc2 <- fwdControl(target.tmp.F_scen2)
681
result.forecast.F_scenario2 <- fwd(data.stf, sr=result.sr.model, control=target.forecast.Fsc2)
682

683
savePic("output/8_Forecasts/F_scenario2/summary_all.png", result.forecast.F_scenario2)
684
savePic("output/8_Forecasts/F_scenario2/summary_short.png", window(result.forecast.F_scenario2, lastYear, lastYear+5))
685

686
saveText("output/8_Forecasts/F_scenario2/ssb.txt", ssb(result.forecast.F_scenario2))
687
saveText("output/8_Forecasts/F_scenario2/rec.txt", rec(result.forecast.F_scenario2))
688
saveText("output/8_Forecasts/F_scenario2/fbar.txt", fbar(result.forecast.F_scenario2))
689
saveText("output/8_Forecasts/F_scenario2/f.txt", harvest(result.forecast.F_scenario2))
690
saveText("output/8_Forecasts/F_scenario2/stock.n.txt", stock.n(result.forecast.F_scenario2))
691
saveText("output/8_Forecasts/F_scenario2/ssb.txt", ssb(result.forecast.F_scenario2))
692
saveText("output/8_Forecasts/F_scenario2/catch.txt", catch(result.forecast.F_scenario2))
693

694
remove(target.tmp.F_scen2)
695

696
# define stock names to summary plot
697
name(result.forecast.F_01manual) <- sprintf("F(0.1)=%s", config.forecast.manual_F01)
698
name(result.forecast.F_MSY) <- sprintf("F(MSY)=%s", config.forecast.manual_FMSY)
699
name(result.forecast.F_scenario1) <- sprintf("F=%s", config.forecast.scenario1.f_year1)
700
name(result.forecast.F_scenario2) <- sprintf("F=%s", config.forecast.scenario2.f_year1)
701
name(result.forecast.F_SQ) <- sprintf("F(SQ)=%s", round(F_SQ, 2))
702

703
# compile summary object with all forecasts
704
result.forecast.all <- FLStocks(result.forecast.F_scenario1, result.forecast.F_scenario2, result.forecast.F_SQ, result.forecast.F_01manual, result.forecast.F_MSY)
705

706
# save output
707
savePic("output/8_Forecasts/summary_all_in_one.png", result.forecast.all)
708
savePic("output/8_Forecasts/short_all_in_one.png", window(result.forecast.all, start=lastYear))
709

710

711
print("[XSA toolkit]: Procedure Done! Take a look on /output/ folder")
712

713
render("Result.Rmd")
714

715