sprat-xsa-experimental

1
f.retro[[i]] <- f.t
2
}
3
result.retro.diagnostic.ssb <- data.frame(ssb.base, ssb.retro[[1]], ssb.retro[[2]], ssb.retro[[3]])
4
config.retro.horizon <- 3
5
print("[XSA toolkit]: Retrospective Mohn's rho tests")
6
# compare XSA estimates in terminal year by -5 year retrospective
7
lastYear <- as.numeric(range(data.stock)["maxyear"])
8
ssb.base <- as.vector(ssb(window(data.stock+FLXSA(data.stock, data.indices, config.tuning.xsa), start=lastYear-3)))
9
f.base <- as.vector(fbar(window(data.stock+FLXSA(data.stock, data.indices, config.tuning.xsa), start=lastYear-3)))
10
ssb.retro <- list()
11
f.retro <- list()
12
for (i in 1:config.retro.horizon) { # prepare retrospective matrix for mohn's rho test
13
# make retro
14
tmp.xsa <- FLXSA(window(data.stock, end=(lastYear-i)), window(data.indices, end=(lastYear-i)), config.tuning.xsa)
15
tmp.xsa.st <- data.stock + tmp.xsa
16
ssb.t <- as.vector(ssb(window(tmp.xsa.st, start=lastYear-config.retro.horizon)))
17
f.t <- as.vector(fbar(window(tmp.xsa.st, start=lastYear-config.retro.horizon)))
18
ssb.retro[[i]] <- ssb.t
19
f.retro[[i]] <- f.t
20
}
21
result.retro.diagnostic.ssb <- data.frame(ssb.base, ssb.retro[[1]], ssb.retro[[2]], ssb.retro[[3]])
22
rownames(result.retro.diagnostic.ssb) <- (lastYear-config.retro.horizon):lastYear
23
result.retro.diagnostic.f <- data.frame(f.base, f.retro[[1]], f.retro[[2]], f.retro[[3]])
24
rownames(result.retro.diagnostic.f) <- (lastYear-config.retro.horizon):lastYear
25
names(result.retro.diagnostic.ssb) <- c("ssb.base", "-1", "-2", "-3", "-4", "-5")
26
result.retro.diagnostic.ssb <- data.frame(ssb.base, ssb.retro[[1]], ssb.retro[[2]], ssb.retro[[3]])
27
rownames(result.retro.diagnostic.ssb) <- (lastYear-config.retro.horizon):lastYear
28
result.retro.diagnostic.f <- data.frame(f.base, f.retro[[1]], f.retro[[2]], f.retro[[3]])
29
rownames(result.retro.diagnostic.f) <- (lastYear-config.retro.horizon):lastYear
30
names(result.retro.diagnostic.ssb) <- c("ssb.base", "-1", "-2", "-3", "-4", "-5")
31
names(result.retro.diagnostic.ssb) <- c("ssb.base", "-1", "-2", "-3")
32
names(result.retro.diagnostic.f) <- c("f.base", "-1", "-2", "-3")
33
# replace NaN's (x/0) to NA (not available)
34
is.nan.data.frame <- function(x) {
35
return (do.call(cbind, lapply(x, is.nan)))
36
}
37
result.retro.diagnostic.f[is.nan(result.retro.diagnostic.f)] <- NA
38
result.retro.diagnostic.ssb[is.nan(result.retro.diagnostic.ssb)] <- NA
39
result.retro.mohnrho.ssb <- mohn(result.retro.diagnostic.ssb, peels = config.retro.horizon)
40
result.retro.mohnrho.fbar <- mohn(result.retro.diagnostic.f, peels = config.retro.horizon)
41
remove(ssb.base, f.base, ssb.retro, f.retro)
42
result.retro.mohnrho.ssb
43
result.retro.mohnrho.fbar
44
source("~/Dev/R/BlackSea/sprat-xsa/run.R")
45
ls
46
# ************************************* #
47
# Author: Piatinskii M., Azov-black sea bench of VNIRO
48
# &copy 2019
49
# **************************************#
50
# make cleanup
51
rm(list = ls(all.names=TRUE))
52
# load libraries
53
library("FLCore")
54
library("FLXSA")
55
library("FLAssess")
56
library("ggplotFL")
57
# library("FLBRP") # no support for R 4.0, require compilation from *nix
58
library("FLash")
59
#install.packages("devtools")
60
#install_github("ices-tools-prod/msy")
61
library("msy")
62
library("rmarkdown")
63
library("dplyr")
64
library("magrittr")
65
library("tibble")
66
library("icesAdvice")
67
fit.xsa
68
source("~/Dev/R/bsea-models/sprat-xsa/run.R")
69
result.xsa.fit
70
AIC(result.xsa.fit)
71
BIC(result.xsa.fit)
72
source("http://flr-project.org/R/instFLR.R")
73
# ************************************* #
74
# Author: Piatinskii M., Azov-black sea bench of VNIRO
75
# &copy 2019
76
# **************************************#
77
# make cleanup
78
rm(list = ls(all.names=TRUE))
79
# load libraries
80
library("FLCore")
81
library("FLXSA")
82
library("FLAssess")
83
library("ggplotFL")
84
# library("FLBRP") # no support for R 4.0, require compilation from *nix
85
library("FLash")
86
#install.packages("devtools")
87
#install_github("ices-tools-prod/msy")
88
library("msy")
89
# stock object configurations
90
config.stock.name <- "Sprattus sprattus" # fish title
91
config.stock.desc <- "Sprat stock information in FLCore format from 2000years" # fish stock description
92
config.stock.fbar.min <- 1 # minimum fish age become fishing (Fbar min)
93
config.stock.fbar.max <- 3 # maximum fish age at fishing (Fbar max)
94
config.stock.no_harvestspawn <- TRUE # Set TRUE if there is no information about fishing before spawning
95
config.stock.no_mspawn <- TRUE # Set TRUE if no information about natural mortality before spawning
96
config.stock.no_discards <- TRUE # Set TRUE if no separated information about landing/discards exist (only catch)
97
config.stock.window.start <- 2001 # model start
98
config.stock.window.end <- 2022 # model end
99
# configure surveys indices
100
config.index1.name <- "Sprat crimea cpue index"
101
config.index1.desc <- "Fishery index information about sprat abundance at Crimea shelf from real fishing surveys"
102
config.index1.startf <- 0.3 # start month of survey in format month/12, example: 1 feb = 2/12 = 0.166
103
config.index1.endf <- 0.7 # end month of survey in format month/12
104
config.index1.type <- "number" # type of survey, allowed: number, biomass
105
# same with previous
106
config.index2.name <- "Sprat caucasian cpue index"
107
config.index2.desc <- "Fishery index information about sprat abundance at Caucasian shelf from real fishing surveys"
108
config.index2.startf <- 0.3
109
config.index2.endf <- 0.7
110
config.index2.type <- "number"
111
# pre-analysis diagnostic procedures. Useful to choice the model parametrization, convergence, etc
112
config.iter.fse = TRUE # Fse choice procedure
113
config.iter.retro = TRUE # retrospective procedure
114
config.iter.fseage = TRUE # Fse at age impact procedure
115
# configure retrospective depth
116
config.retro.horizon <- 5
117
# Shrinkage standard error max value. Default value Fse = 0.5
118
# Higher shrinkage se leads to high unrestraint of estimating procedure.
119
config.tuning.xsa <- FLXSA.control(
120
maxit = 100, # maximum covergence iterations
121
fse = 1.5, # shrinkage s.e value
122
shk.n = T,
123
shk.f = T,
124
shk.yrs = 4,
125
tsrange = 15, # weighting time series length (from latest year)
126
tspower = 3, # weighting procedure sequence (1 = linear, 2 = square, 3 = cubic)
127
qage = 2, # age after catchability not depends on age
128
)
129
# sr model configuration
130
config.srmodel.eqsim.use = TRUE # should we use eqsim package insted of FLSR default?
131
config.srmodel.use_geomean = FALSE # set FALSE to use EQSIM SR model results to forecasting
132
# eqsim model configurations (if config.srmodel.eqsim.use = TRUE)
133
config.srmodel.eqsim.models = c("Segreg") # types of model to fit and compare
134
config.srmodel.eqsim.simulations = 10 # number of simulations
135
config.srmodel.eqsim.shift_years <- c(1994:2003) # remove years from sr-model fitting
136
# flsr model configuration (if config.srmodel.eqsim.use = FALSE)
137
config.srmodel.type <- segregAR1() # SR model type: ricker(), bevholt(), geomean(), segreg(), shepherd(), cushing(), AR in model name = AutoCorrelation fix
138
config.srmodel.optimization = FALSE # Use optimization procedure for SR model fitting?
139
# configure MSY reference point calculation: Blim, Bpa, Fmsy
140
config.brp.blim_multiplier <- 0.2 # multiplier of SSBmax to get Blim ref.point, allowed: 0.2 - 0.4
141
config.brp.trim_uncertainity <- TRUE # trim uncerteinity extreme deviations? Optimization procedure for Bpa estimate
142
config.brp.current_range <- 2003:2021 # shift previous population states from BRP calculations
143
# Custom forecast scenarious. By default performed: F_SQ, F_01, F_msy scenarious.
144
# There you can define f(y+1), f(y+2), f(y+3) targets for 2 scenarious
145
config.forecast.scenario1.f_year1 <- 0.4
146
config.forecast.scenario1.f_year2 <- 0.4
147
config.forecast.scenario1.f_year3 <- 0.4
148
config.forecast.scenario2.f_year1 <- 0.5
149
config.forecast.scenario2.f_year2 <- 0.5
150
config.forecast.scenario2.f_year3 <- 0.5
151
# additional scenario for manual Fmsy and F_01 estimates
152
config.forecast.manual_F01 <- 0.605
153
config.forecast.manual_FMSY <- 0.64
154
# custom calculations. Calculate partial recruitment for Y/R analysis
155
config.custom.partial.years <- 2004:2020
156
# custom calculations. Calculate F01 reference point from Y/R analysis
157
config.custom.f01.years <- 2004:2020 # years used to calculate mean(wt), mean(Fi/Fmax), mean(M)
158
config.custom.f01.plusgroup = TRUE # is last group is plus group?
159
config.custom.f01.oldest = 5 # oldest fish found in waters (not in stock, that means "even somewhen exist in sea")
160
# ******** GLOBAL CONFIGURATION DONE!!! **** #
161
# ********************************************************* #
162
# Do not change enything below if you have no idea what is going on
163
# ********************************************************* #
164
# simple write functions
165
savePic <- function(path, imgObj, width=1900, height=1600, res=200) {
166
png(path, width=width, height=height, res=res, type="cairo")
167
print(plot(imgObj) + theme_bw())
168
dev.off()
169
}
170
saveText <- function(path, object) {
171
con <- file(path, encoding = "UTF-8")
172
sink(con)
173
print(object)
174
sink()
175
close(con)
176
}
177
# get mbar from stock object
178
mbar <- function(stock) {
179
stock.min <- as.numeric(stock@range["min"])
180
stock.max <- as.numeric(stock@range["max"])
181
bar.min <- as.numeric(stock@range["minfbar"])
182
bar.max <- as.numeric(stock@range["maxfbar"])
183
diff.min <- 1 + bar.min - stock.min
184
diff.max <- 1 + bar.max - bar.min
185
if (diff.max > stock.max+1)
186
diff.max <- stock.max+1
187
vec <- (m(stock)[diff.min:diff.max,])
188
bar <- apply(vec, 2, mean)
189
return(bar)
190
}
191
# read input data - stock
192
tmp.stock <- read.csv("input/stock.csv", sep=",")
193
data.stock <- as.FLStock(tmp.stock)
194
# read input data - index
195
tmp.idx1.idx <- read.csv("input/survey1/index.csv", sep=",")
196
tmp.idx1.eff <- read.csv("input/survey1/effort.csv", sep=",")
197
tmp.idx2.idx <- read.csv("input/survey2/index.csv", sep=",")
198
tmp.idx2.eff <- read.csv("input/survey2/effort.csv", sep=",")
199
data.index1 <- FLIndex(index=as.FLQuant(tmp.idx1.idx), effort=as.FLQuant(tmp.idx1.eff), catch.n=as.FLQuant(tmp.idx1.idx))
200
data.index2 <- FLIndex(index=as.FLQuant(tmp.idx2.idx), effort=as.FLQuant(tmp.idx2.eff), catch.n=as.FLQuant(tmp.idx2.idx))
201
remove(tmp.stock, tmp.idx1.idx, tmp.idx1.eff, tmp.idx2.idx, tmp.idx2.eff)
202
print("[XSA toolkit]: Load data: OK")
203
# Set name,desc from cfgs
204
name(data.stock) <- config.stock.name
205
desc(data.stock) <- config.stock.desc
206
# set fbar ranges
207
range(data.stock)[c("minfbar", "maxfbar")] <- c(config.stock.fbar.min, config.stock.fbar.max)
208
# shift data starting from 2000
209
data.stock <- window(data.stock, start = config.stock.window.start, end = config.stock.window.end)
210
# set harvest.spwn and m.spwn = 0 if no data exist (or XSA throw error)
211
if (config.stock.no_harvestspawn) {
212
harvest.spwn(data.stock) <- 0
213
}
214
if (config.stock.no_mspawn) {
215
m.spwn(data.stock) <- 0
216
}
217
# if no info about discards/landings separatly than set discards = 0, landings = catch
218
if (config.stock.no_discards) {
219
discards(data.stock) <- 0
220
discards.n(data.stock) <- 0
221
discards.wt(data.stock) <- catch.wt(data.stock)
222
landings(data.stock) <- catch(data.stock)
223
landings.n(data.stock) <- catch.n(data.stock)
224
landings.wt(data.stock) <- catch.wt(data.stock)
225
}
226
stock.wt(data.stock) <- catch.wt(data.stock)
227
# set cfgs for index1 and index2
228
type(data.index1) <- config.index1.type
229
name(data.index1) <- config.index1.name
230
desc(data.index1) <- config.index1.desc
231
range(data.index1)[c("startf", "endf")] <- c(config.index1.startf, config.index1.endf)
232
type(data.index2) <- config.index2.type
233
name(data.index2) <- config.index2.name
234
desc(data.index2) <- config.index2.desc
235
range(data.index2)[c("startf", "endf")] <- c(config.index2.startf, config.index2.endf)
236
print("[XSA toolkit]: FLStock and FLIndex object configured successful")
237
data.indices <- FLIndices(data.index1, data.index2)
238
data.stock
239
data.stock@catch
240
data.stock@catch.n
241
data.stock@catch.wt
242
format(data.stock@catch.wt, scientific = F)
243
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
244
library("devtools")
245
install_github("ices-tools-prod/msy")
246
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
247
data.stock
248
# make Fse comparision test
249
if (config.iter.fse) {
250
print("[XSA toolkit]: Fse impact factor estimate")
251
steps <- seq(0.5, 2.5, by = 0.5)
252
fseIters <- c()
253
for (i in 1:length(steps)) {
254
ctl <- config.tuning.xsa
255
ctl@fse <- steps[i]
256
res <- FLXSA(data.stock, data.indices, ctl) + data.stock
257
fseIters <- c(fseIters, res)
258
}
259
result.fse.summary <- FLStocks(fseIters)
260
names(result.fse.summary) <- steps
261
# save output to pic
262
savePic("output/1_FseChoose/summary_all.png", result.fse.summary)
263
# get latest 5 year window
264
lastYear <- result.fse.summary[[1]]@range["maxyear"]
265
tmpObj <- window(result.fse.summary, lastYear-config.retro.horizon, lastYear)
266
savePic("output/1_FseChoose/summary_short.png", tmpObj)
267
remove(lastYear, tmpObj)
268
# save text output - ssb
269
con <- file("output/1_FseChoose/ssb.txt", encoding = "UTF-8")
270
sink(con)
271
for (i in 1:length(steps)) {
272
print(sprintf("============= SSB at Fse = %s", steps[i]))
273
print(ssb(result.fse.summary[[i]]))
274
cat("\r\n")
275
}
276
sink()
277
close(con)
278
# save text output - rec
279
con <- file("output/1_FseChoose/rec.txt", encoding = "UTF-8")
280
sink(con)
281
for (i in 1:length(steps)) {
282
print(sprintf("============= Recruitment at Fse = %s", steps[i]))
283
print(rec(result.fse.summary[[i]]))
284
cat("\r\n")
285
}
286
sink()
287
close(con)
288
# save text output - fbar
289
con <- file("output/1_FseChoose/fbar.txt", encoding = "UTF-8")
290
sink(con)
291
for (i in 1:length(steps)) {
292
print(sprintf("============= Fbar at Fse = %s", steps[i]))
293
print(fbar(result.fse.summary[[i]]))
294
cat("\r\n")
295
}
296
sink()
297
close(con)
298
remove(steps, fseIters, ctl, res, con)
299
}
300
source("~/Dev/R/phdpm2/sprat-xsa-experimental/clean-output.R")
301
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
302
source("~/Dev/R/phdpm2/sprat-xsa-experimental/clean-output.R")
303
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
304
source("~/Dev/R/phdpm2/sprat-xsa-experimental/clean-output.R")
305
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
306
print(paste0("Rho ssb=", round(result.retro.mohnrho.ssb, 3)))
307
print(paste0("Rho F=", round(result.retro.mohnrho.fbar, 3)))
308
?FLXSA.control
309
data.stock
310
data.stock@catch.n
311
# ************************************* #
312
# Author: Piatinskii M., Azov-black sea bench of VNIRO
313
# &copy 2019
314
# **************************************#
315
# make cleanup
316
rm(list = ls(all.names=TRUE))
317
# load libraries
318
library("FLCore")
319
library("FLXSA")
320
library("FLAssess")
321
library("ggplotFL")
322
# library("FLBRP") # no support for R 4.0, require compilation from *nix
323
library("FLash")
324
#install.packages("devtools")
325
#install_github("ices-tools-prod/msy")
326
library("msy")
327
library("rmarkdown")
328
library("dplyr")
329
library("magrittr")
330
library("tibble")
331
library("icesAdvice")
332
# ****** GLOBAL CONFIGURATIONS SECTION *******#
333
# stock object configurations
334
config.stock.name <- "Sprattus sprattus" # fish title
335
config.stock.desc <- "Sprat stock information in FLCore format from 2000years" # fish stock description
336
config.stock.fbar.min <- 1 # minimum fish age become fishing (Fbar min)
337
config.stock.fbar.max <- 3 # maximum fish age at fishing (Fbar max)
338
config.stock.no_harvestspawn <- TRUE # Set TRUE if there is no information about fishing before spawning
339
config.stock.no_mspawn <- TRUE # Set TRUE if no information about natural mortality before spawning
340
config.stock.no_discards <- TRUE # Set TRUE if no separated information about landing/discards exist (only catch)
341
config.stock.window.start <- 1994 # model start
342
config.stock.window.end <- 2022 # model end
343
# configure surveys indices
344
config.index1.name <- "Sprat crimea cpue index"
345
config.index1.desc <- "Fishery index information about sprat abundance at Crimea shelf from real fishing surveys"
346
config.index1.startf <- 0.3 # start month of survey in format month/12, example: 1 feb = 2/12 = 0.166
347
config.index1.endf <- 0.7 # end month of survey in format month/12
348
config.index1.type <- "number" # type of survey, allowed: number, biomass
349
# same with previous
350
config.index2.name <- "Sprat caucasian cpue index"
351
config.index2.desc <- "Fishery index information about sprat abundance at Caucasian shelf from real fishing surveys"
352
config.index2.startf <- 0.3
353
config.index2.endf <- 0.7
354
config.index2.type <- "number"
355
# pre-analysis diagnostic procedures. Useful to choice the model parametrization, convergence, etc
356
config.iter.fse = TRUE # Fse choice procedure
357
config.iter.retro = TRUE # retrospective procedure
358
config.iter.fseage = TRUE # Fse at age impact procedure
359
# configure retrospective depth
360
config.retro.horizon <- 5
361
# Shrinkage standard error max value. Default value Fse = 0.5
362
# Higher shrinkage se leads to high unrestraint of estimating procedure.
363
config.tuning.xsa <- FLXSA.control(
364
maxit = 100, # maximum covergence iterations
365
fse = 1.5, # shrinkage s.e value
366
shk.n = T,
367
shk.f = T,
368
shk.yrs = 4,
369
tsrange = 15, # weighting time series length (from latest year)
370
tspower = 3, # weighting procedure sequence (1 = linear, 2 = square, 3 = cubic)
371
qage = 2, # age after catchability not depends on age
372
)
373
# sr model configuration
374
config.srmodel.eqsim.use = TRUE # should we use eqsim package insted of FLSR default?
375
config.srmodel.use_geomean = FALSE # set FALSE to use EQSIM SR model results to forecasting
376
# eqsim model configurations (if config.srmodel.eqsim.use = TRUE)
377
config.srmodel.eqsim.models = c("Segreg") # types of model to fit and compare
378
config.srmodel.eqsim.simulations = 10 # number of simulations
379
config.srmodel.eqsim.shift_years <- 1994:2005 # remove years from sr-model fitting
380
# flsr model configuration (if config.srmodel.eqsim.use = FALSE)
381
config.srmodel.type <- segregAR1() # SR model type: ricker(), bevholt(), geomean(), segreg(), shepherd(), cushing(), AR in model name = AutoCorrelation fix
382
config.srmodel.optimization = FALSE # Use optimization procedure for SR model fitting?
383
# configure MSY reference point calculation: Blim, Bpa, Fmsy
384
config.brp.blim_multiplier <- 0.2 # multiplier of SSBmax to get Blim ref.point, allowed: 0.2 - 0.4
385
config.brp.trim_uncertainity <- TRUE # trim uncerteinity extreme deviations? Optimization procedure for Bpa estimate
386
config.brp.current_range <- 2003:2021 # shift previous population states from BRP calculations
387
# Custom forecast scenarious. By default performed: F_SQ, F_01, F_msy scenarious.
388
# There you can define f(y+1), f(y+2), f(y+3) targets for 2 scenarious
389
# precautionary scenario: Fpa = FMSY * (B2022/Btr) = 0.64 * 66/87.5 = 0.482
390
config.forecast.scenario1.f_year1 <- 0.48
391
config.forecast.scenario1.f_year2 <- 0.48
392
config.forecast.scenario1.f_year3 <- 0.48
393
# Scenario MSY at M = 0,95
394
# by old STECF (2010) results
395
config.forecast.scenario2.f_year1 <- 0.95
396
config.forecast.scenario2.f_year2 <- 0.95
397
config.forecast.scenario2.f_year3 <- 0.95
398
# additional scenario for manual Fmsy and F_01 estimates
399
config.forecast.manual_F01 <- 0.605
400
config.forecast.manual_FMSY <- 0.64
401
# custom calculations. Calculate partial recruitment for Y/R analysis
402
config.custom.partial.years <- 2004:2020
403
# custom calculations. Calculate F01 reference point from Y/R analysis
404
config.custom.f01.years <- 2004:2020 # years used to calculate mean(wt), mean(Fi/Fmax), mean(M)
405
config.custom.f01.plusgroup = TRUE # is last group is plus group?
406
config.custom.f01.oldest = 5 # oldest fish found in waters (not in stock, that means "even somewhen exist in sea")
407
# ******** GLOBAL CONFIGURATION DONE!!! **** #
408
# ********************************************************* #
409
# Do not change enything below if you have no idea what is going on
410
# ********************************************************* #
411
# simple write functions
412
savePic <- function(path, imgObj, width=1900, height=1600, res=200) {
413
png(path, width=width, height=height, res=res, type="cairo")
414
print(plot(imgObj) + theme_bw())
415
dev.off()
416
}
417
saveText <- function(path, object) {
418
con <- file(path, encoding = "UTF-8")
419
sink(con)
420
print(object)
421
sink()
422
close(con)
423
}
424
# get mbar from stock object
425
mbar <- function(stock) {
426
stock.min <- as.numeric(stock@range["min"])
427
stock.max <- as.numeric(stock@range["max"])
428
bar.min <- as.numeric(stock@range["minfbar"])
429
bar.max <- as.numeric(stock@range["maxfbar"])
430
diff.min <- 1 + bar.min - stock.min
431
diff.max <- 1 + bar.max - bar.min
432
if (diff.max > stock.max+1)
433
diff.max <- stock.max+1
434
vec <- (m(stock)[diff.min:diff.max,])
435
bar <- apply(vec, 2, mean)
436
return(bar)
437
}
438
# read input data - stock
439
tmp.stock <- read.csv("input/stock.csv", sep=",")
440
data.stock <- as.FLStock(tmp.stock)
441
# read input data - index
442
tmp.idx1.idx <- read.csv("input/survey1/index.csv", sep=",")
443
tmp.idx1.eff <- read.csv("input/survey1/effort.csv", sep=",")
444
tmp.idx2.idx <- read.csv("input/survey2/index.csv", sep=",")
445
tmp.idx2.eff <- read.csv("input/survey2/effort.csv", sep=",")
446
data.index1 <- FLIndex(index=as.FLQuant(tmp.idx1.idx), effort=as.FLQuant(tmp.idx1.eff), catch.n=as.FLQuant(tmp.idx1.idx))
447
data.index2 <- FLIndex(index=as.FLQuant(tmp.idx2.idx), effort=as.FLQuant(tmp.idx2.eff), catch.n=as.FLQuant(tmp.idx2.idx))
448
remove(tmp.stock, tmp.idx1.idx, tmp.idx1.eff, tmp.idx2.idx, tmp.idx2.eff)
449
print("[XSA toolkit]: Load data: OK")
450
# Set name,desc from cfgs
451
name(data.stock) <- config.stock.name
452
desc(data.stock) <- config.stock.desc
453
# set fbar ranges
454
range(data.stock)[c("minfbar", "maxfbar")] <- c(config.stock.fbar.min, config.stock.fbar.max)
455
# shift data starting from 2000
456
data.stock <- window(data.stock, start = config.stock.window.start, end = config.stock.window.end)
457
# set harvest.spwn and m.spwn = 0 if no data exist (or XSA throw error)
458
if (config.stock.no_harvestspawn) {
459
harvest.spwn(data.stock) <- 0
460
}
461
if (config.stock.no_mspawn) {
462
m.spwn(data.stock) <- 0
463
}
464
# if no info about discards/landings separatly than set discards = 0, landings = catch
465
if (config.stock.no_discards) {
466
discards(data.stock) <- 0
467
discards.n(data.stock) <- 0
468
discards.wt(data.stock) <- catch.wt(data.stock)
469
landings(data.stock) <- catch(data.stock)
470
landings.n(data.stock) <- catch.n(data.stock)
471
landings.wt(data.stock) <- catch.wt(data.stock)
472
}
473
stock.wt(data.stock) <- catch.wt(data.stock)
474
# set cfgs for index1 and index2
475
type(data.index1) <- config.index1.type
476
name(data.index1) <- config.index1.name
477
desc(data.index1) <- config.index1.desc
478
range(data.index1)[c("startf", "endf")] <- c(config.index1.startf, config.index1.endf)
479
type(data.index2) <- config.index2.type
480
name(data.index2) <- config.index2.name
481
desc(data.index2) <- config.index2.desc
482
range(data.index2)[c("startf", "endf")] <- c(config.index2.startf, config.index2.endf)
483
print("[XSA toolkit]: FLStock and FLIndex object configured successful")
484
data.stock@catch.n
485
data.stock@catch.wt
486
data.stock@catch.wt*1000
487
data.stock@catch.wt*1000000
488
data.stock@catch
489
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
490
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
491
data.stock
492
data.stock@catch.n * data.stock@catch.wt
493
format(data.stock@catch.wt*1000000, scientific = F)
494
data.stock@catch.wt*1000000
495
data.stock@catch.wt*data.stock@catch.n
496
data.stock@catch.n
497
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
498
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
499
source("~/Dev/R/phdpm2/sprat-xsa-experimental/run.R")
500
result.forecast.F_scenario1
501
result.forecast.F_scenario1@stock.n
502
result.forecast.F_MSY@harvest
503
result.forecast.F_MSY@stock.n
504
result.forecast.F_MSY@harvest
505
19*0.75
506
27*0.75
507
26*0.25
508
26*0.75
509
15*0.75
510
20*0.75
511
46+15
512
61*0.75
513