tests/testthat/test-FLStock_projection.R
In iagomosqueira/FLasher: Projection and Forecasting of Fish Populations, Stocks and Fleets
# Maintainer: Finlay Scott, JRC
# Distributed under the terms of the European Union Public Licence (EUPL) V.1.1.
context("Projections with FLStock")
# If these don't work we are in trouble!
source("expect_funs.R")
test_that("Catch target, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_val <- rlnorm(n=nyears, mean=log(min(catch(ple4)/10)), sd=0.1)
control=fwdControl(data.frame(year=years, quant="catch", value=catch_val))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
catch_out <- c(catch(res)[,ac(years)])
# Test to within tolerance of 1.5e-4
expect_equal(catch_out, catch_val, tolerance=1.5e-4)
})
test_that("Catch target, multiple iters",{
data(ple4)
niters <- round(runif(1, min=10, max=50))
ple4p <- propagate(ple4, niters)
# Spoof up stock.n, landings.wt, landings.n and harvest
harvest(ple4p)[] <- rlnorm(n=prod(dim(harvest(ple4p))), mean=log(c(harvest(ple4p))), sd=0.001)
stock.n(ple4p)[] <- rlnorm(n=prod(dim(stock.n(ple4p))), mean=log(c(stock.n(ple4p))), sd=0.001)
landings.wt(ple4p)[] <- rlnorm(n=prod(dim(landings.wt(ple4p))), mean=log(c(landings.wt(ple4p))), sd=0.001)
landings.n(ple4p)[] <- rlnorm(n=prod(dim(landings.n(ple4p))), mean=log(c(landings.n(ple4p))), sd=0.001)
year_range <- range(ple4p)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_val <- rlnorm(n=nyears*niters, mean=log(min(catch(ple4p)/10)), sd=0.1)
control=fwdControl(data.frame(year=years, quant="catch", value=0), iters=niters)
control@iters[,"value",] <- catch_val
res <- fwd(ple4p, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
catch_out <- c(catch(res)[,ac(years)])
# Test to within tolerance of 1.5e-8
expect_equal(catch_out, catch_val, tol=1e-6)
})
test_that("Fbar target, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
f_val <- rlnorm(n=nyears, mean=log(min(fbar(ple4)/10)), sd=0.1)
control=fwdControl(data.frame(year=years, quant="fbar", value=f_val))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
f_out <- c(fbar(res)[,ac(years)])
# Test to within tolerance of 1.5e-8
expect_equal(f_out, f_val)
})
test_that("Fbar target, multiple iters",{
data(ple4)
niters <- round(runif(1, min=10, max=50))
ple4p <- propagate(ple4, niters)
# Spoof up stock.n, landings.wt, landings.n and harvest
harvest(ple4p)[] <- rlnorm(n=prod(dim(harvest(ple4p))), mean=log(c(harvest(ple4p))), sd=0.001)
stock.n(ple4p)[] <- rlnorm(n=prod(dim(stock.n(ple4p))), mean=log(c(stock.n(ple4p))), sd=0.001)
landings.wt(ple4p)[] <- rlnorm(n=prod(dim(landings.wt(ple4p))), mean=log(c(landings.wt(ple4p))), sd=0.001)
landings.n(ple4p)[] <- rlnorm(n=prod(dim(landings.n(ple4p))), mean=log(c(landings.n(ple4p))), sd=0.001)
year_range <- range(ple4p)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
f_val <- rlnorm(n=nyears*niters, mean=log(min(fbar(ple4p)/10)), sd=0.1)
control=fwdControl(data.frame(year=years, quant="fbar", value=0), iters=niters)
control@iters[,"value",] <- f_val
res <- fwd(ple4p, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
f_out <- c(fbar(res)[,ac(years)])
# Test to within tolerance of 1.5e-8
expect_equal(f_out, f_val)
})
test_that("Catch target with min limit, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_val <- rlnorm(n=nyears, mean=log(min(catch(ple4)/10)), sd=0.1)
catch_lim <- mean(catch_val)
# First project without limit
control=fwdControl(list(year=years, quant="catch", value=catch_val))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean", params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# Get resulting fbar
fout <- fbar(res)[,ac(years)]
limit_year <- sample(years,1)
flim <- c(fout[,ac(limit_year)] * 1.5)
# New projection
control=fwdControl(list(year=years, quant="catch", value=catch_val),
list(year=limit_year, quant="fbar", min=flim))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean", params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# All years apart from lim_year met catch OK
catch_out <- catch(res)[,ac(years)]
non_lim_years <- !(years %in% limit_year)
expect_equal(c(catch_out[,ac(years[non_lim_years])]), catch_val[non_lim_years],
tolerance=0.0005)
# lim year has f = flim
expect_equal(c(fbar(res)[,ac(limit_year)]), flim)
})
test_that("Catch target with max limit, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_val <- rlnorm(n=nyears, mean=log(min(catch(ple4)/10)), sd=0.1)
catch_lim <- mean(catch_val)
# First project without limit
control=fwdControl(list(year=years, quant="catch", value=catch_val))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean", params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# Get resulting fbar
fout <- fbar(res)[,ac(years)]
limit_year <- sample(years,1)
flim <- c(fout[,ac(limit_year)] * 0.5)
# New projection
control=fwdControl(list(year=years, quant="catch", value=catch_val),
list(year=limit_year, quant="fbar", max=flim))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean", params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# All years apart from lim_year met catch OK
catch_out <- catch(res)[,ac(years)]
non_lim_years <- !(years %in% limit_year)
expect_true(all(c( catch_out[,ac(years[non_lim_years])]) -
catch_val[non_lim_years] < 1e-8))
# lim year has f = flim
expect_equal(c(fbar(res)[,ac(limit_year)]), flim)
})
test_that("Catch relative target, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1]+1, max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_rel_val <- runif(nyears, min=0.5, max=0.9)
control=fwdControl(data.frame(year=years, quant="catch", value=catch_rel_val, relYear=years-1))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
catch_out <- c(catch(res)[,ac(years)])
catch_trg <- c(catch(ple4)[,ac(years[1]-1)]) * cumprod(catch_rel_val)
# Test to within tolerance of 1.5e-8
expect_equal(catch_out, catch_trg, tol=1.5e-6)
})
# DEBUG
# test_that("Fbar relative target, single iter",{
# data(ple4)
# year_range <- range(ple4)[c("minyear","maxyear")]
# min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
# min_year[1] <- min_year[1] + 1
# nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
# years <- min_year:(min_year+nyears-1)
# f_rel_val <- runif(nyears, min=0.5, max=0.9)
# control <- fwdControl(data.frame(year=years, quant="fbar",
# value=f_rel_val, relYear=years-1))
# res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
# params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# f_out <- c(fbar(res)[,ac(years)])
# f_trg <- c(fbar(ple4)[,ac(years[1]-1)]) * cumprod(f_rel_val)
# # Test to within tolerance of 1.5e-8
# expect_equal(f_out, f_trg)
# })
#
#---------------------
test_that("Tests from Running Medium Term Forecasts with FLasher tutorial",{
data(ple4)
ple4_mtf <- stf(ple4, nyears = 10)
ple4_sr <- fmle(as.FLSR(ple4, model="bevholt"), control=list(trace=0))
# Simple F
f_status_quo <- mean(fbar(ple4)[,as.character(2005:2008)])
ctrl_f <- fwdControl(list(year=2009:2018, quant="f", value=f_status_quo))
ple4_f_sq <- fwd(ple4_mtf, control = ctrl_f, sr = ple4_sr)
expect_equal(c(fbar(ple4_f_sq)[,ac(2009:2018)]), rep(f_status_quo,10))
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_f_sq)[,ac(2008:2017)])), c(rec(ple4_f_sq)[,ac(2009:2018)]))
# Decreasing catch
future_catch <- c(catch(ple4)[,"2008"]) * 0.9^(1:10)
ctrl_catch <- fwdControl(list(year=2009:2018, quant = "catch", value=future_catch))
ple4_catch <- fwd(ple4_mtf, control = ctrl_catch, sr = ple4_sr)
expect_equal(c(catch(ple4_catch)[,ac(2009:2018)]), future_catch)
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_catch)[,ac(2008:2017)])), c(rec(ple4_catch)[,ac(2009:2018)]))
# SSB end
final_ssb <- 100000
ctrl_ssb <- fwdControl(list(year=2009, quant = "ssb_end", value=final_ssb))
ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
survivors <- stock.n(ple4_ssb) * exp(-harvest(ple4_ssb) - m(ple4_ssb))
expect_equal(c(quantSums((survivors * stock.wt(ple4_ssb) * mat(ple4_ssb))[,ac(2009)])), final_ssb)
# SSB at spawn
spawn_ssb <- 200000
ctrl_ssb <- fwdControl(list(year=2009, quant = "ssb_spawn", value=spawn_ssb))
expect_warning(fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr))
m.spwn(ple4_mtf)[,ac(2009)] <- 0.5
harvest.spwn(ple4_mtf)[,ac(2009)] <- 0.5
ctrl_ssb <- fwdControl(data.frame(year=2009, quant = "ssb_spawn", value=spawn_ssb))
ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
expect_equal(c(ssb(ple4_ssb)[,ac(2009)]), spawn_ssb)
# SRP
srp <- 200000
ctrl_ssb <- fwdControl(data.frame(year=2009, quant = "srp", value=srp))
ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
expect_equal(c(ssb(ple4_ssb)[,ac(2009)]), srp)
# FLash-like
m.spwn(ple4_mtf)[,ac(2009)] <- 0.5
harvest.spwn(ple4_mtf)[,ac(2009)] <- 0.5
flash_ssb <- 150000
ctrl_ssb <- fwdControl(data.frame(year=2009, quant = "ssb_flash", value=flash_ssb))
ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
expect_equal(flash_ssb, c(ssb(ple4_ssb)[,ac(2009)]))
# m.spwn(ple4_mtf)[,ac(2009)] <- 0.0
# harvest.spwn(ple4_mtf)[,ac(2009)] <- 0.0
# flash_ssb <- 150000
# ctrl_ssb <- fwdControl(data.frame(year=2009, quant = "ssb_flash",
# value=flash_ssb))
# ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
# expect_equal(c(ssb(ple4_ssb)[,ac(2010)]), flash_ssb)
#
# Longer SSB
m.spwn(ple4_mtf)[,ac(dims(ple4)$maxyear + 1)] <- 0.0
harvest.spwn(ple4_mtf)[,ac(dims(ple4)$maxyear + 1)] <- 0.0
future_ssb <- c(ssb(ple4)[,'2017']) * 0.8
ctrl_ssb <- fwdControl(data.frame(year=seq(2018, 2027), quant = "ssb_flash", value=future_ssb))
expect_warning(ple4_ssb <- fwd(ple4_mtf, control = ctrl_ssb, sr = ple4_sr))
expect_equal(c(ssb(ple4_ssb)[,ac(2019:2027)]), rep(future_ssb, 9))
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_ssb)[,ac(2019:2026)])), c(rec(ple4_ssb)[,ac(2020:2027)]))
# Relative catch
ctrl_rel_catch <- fwdControl(data.frame(
year = 2009:2018, quant = "catch", value = 0.9, relYear = 2008:2017))
ple4_rel_catch <- fwd(ple4_mtf, control = ctrl_rel_catch, sr = ple4_sr)
expect_equal(c(catch(ple4_rel_catch)[,ac(2009:2018)] / catch(ple4_rel_catch)[,ac(2008:2017)]),
rep(0.9, 10))
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_rel_catch)[,ac(2008:2017)])),
c(rec(ple4_rel_catch)[,ac(2009:2018)]))
# Min Max
f01 <- 0.1
min_catch <- mean(catch(ple4_mtf)[,as.character(2006:2008)])
ctrl_min_catch <- fwdControl(
list(year=2009:2018, quant="f", value=f01), list(year=2009:2018, quant="catch", min=min_catch))
ple4_min_catch <- fwd(ple4_mtf, control = ctrl_min_catch, sr = ple4_sr)
ple4c <- catch(ple4_min_catch)[,ac(2009:2018)]
# Expect a solving tolerance of better than 1.5e-8
expect_true(all((ple4c - min_catch) < 1.5e-8))
# Relative targets and bounds
current_fbar <- c(fbar(ple4)[,"2008"])
f_target <- c(seq(from = current_fbar, to = f01, length = 8)[-1], rep(f01, 3))
rel_catch_bound <- 0.10
ctrl_rel_min_max_catch <- fwdControl(
list(year=2009:2018, quant="f", value=f_target),
list(year=2009:2018, quant="catch", relYear=2008:2017, max=1+rel_catch_bound, min=1-rel_catch_bound))
recovery<-fwd(ple4_mtf, control=ctrl_rel_min_max_catch, sr=ple4_sr)
ple4c <- catch(recovery)[,ac(2009:2018)]
relcatch <- catch(recovery)[,ac(2009:2018)] / catch(recovery)[,ac(2008:2017)]
# All relative catches within bounds within tolerance
expect_true(all(relcatch > (1-rel_catch_bound - 1e6) & relcatch < (1+rel_catch_bound + 1e6)))
# Prepare iterations
niters <- 20
ple4_mtf <- stf(ple4, nyears = 10)
ple4_mtf <- propagate(ple4_mtf, niters)
# deviances
rec_deviances <- FLQuant(NA, dimnames = list(year=2009:2018, iter=1:niters))
sample_years <- sample(dimnames(residuals(ple4_sr))$year, niters * 10, replace = TRUE)
rec_deviances[] <- exp(residuals(ple4_sr)[,sample_years])
ple4_stoch_rec <- fwd(ple4_mtf, control = ctrl_catch, sr = ple4_sr, deviances = rec_deviances)
expect_equal(c(catch(ple4_stoch_rec)[,ac(2009:2018)]),rep(future_catch,niters))
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_stoch_rec)[,ac(2008:2017)]) %*% rec_deviances), c(rec(ple4_stoch_rec)[,ac(2009:2018)]))
# Stochastic target
stoch_catch <- rlnorm(10*niters, meanlog=log(future_catch), sdlog=0.3)
ctrl_catch_iters <- fwdControl(list(year=2009:2018, quant="catch", value=stoch_catch))
ple4_catch_iters <- fwd(ple4_mtf, control=ctrl_catch_iters, sr = ple4_sr)
expect_equal(c(catch(ple4_catch_iters)[,ac(2009:2018)]), stoch_catch)
# Iterations in recruitment parameters
sr_iters <- FLPar(NA, dimnames=list(params=c("a","b"), iter=1:niters))
aiters <- rlnorm(niters, meanlog=log(params(ple4_sr)["a"]), sdlog=0.5)
biters <- rlnorm(niters, meanlog=log(params(ple4_sr)["b"]), sdlog=0.01)
sr_iters["a"] <- aiters
sr_iters["b"] <- biters
ple4_sr_iters <- fwd(ple4_mtf, control=ctrl_catch, sr = list(model="bevholt", params=sr_iters))
expect_equal(c(catch(ple4_sr_iters)[,ac(2009:2018)]), rep(future_catch, niters))
srmod <- predictModel(model=bevholt()$model, params=sr_iters)
expect_equal(c(predict(srmod, ssb=ssb(ple4_sr_iters)[,ac(2008:2017)])), c(rec(ple4_sr_iters)[,ac(2009:2018)]))
# Stochastic targets and recruitment
ple4_iters <- fwd(ple4_mtf, control=ctrl_catch_iters, sr = list(model="bevholt", params=sr_iters), deviances = rec_deviances)
expect_equal(c(catch(ple4_iters)[,ac(2009:2018)]), stoch_catch)
expect_equal(c(predict(srmod, ssb=ssb(ple4_iters)[,ac(2008:2017)]) %*% rec_deviances), c(rec(ple4_iters)[,ac(2009:2018)]))
})
test_that("Catch target, maxf limiting",{
data(ple4)
control=fwdControl(year=1990, quant="catch", value=150000)
res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))), maxF=0.1)
f_out <- c(fbar(res)[,"1990"])
expect_equal(f_out, 0.1)
})
test_that("control as FLQuants works", {
data(ple4)
# catch
res <- fwd(ple4, catch=catch(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(catch(ple4)[,ac(2000:2008)], catch(res)[,ac(2000:2008)])
# landings
res <- fwd(ple4, landings=landings(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(landings(ple4)[,ac(2000:2008)], landings(res)[,ac(2000:2008)])
# discards
res <- fwd(ple4, discards=discards(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(discards(ple4)[,ac(2000:2008)], discards(res)[,ac(2000:2008)])
# fbar
res <- fwd(ple4, fbar=fbar(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(fbar(ple4)[,ac(2000:2008)], fbar(res)[,ac(2000:2008)])
# f
res <- fwd(ple4, f=fbar(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(fbar(ple4)[,ac(2000:2008)], fbar(res)[,ac(2000:2008)])
})
iagomosqueira/FLasher documentation built on Jan. 11, 2024, 12:58 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Maintainer: Finlay Scott, JRC
# Distributed under the terms of the European Union Public Licence (EUPL) V.1.1.
context("Projections with FLStock")
# If these don't work we are in trouble!
source("expect_funs.R")
test_that("Catch target, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_val <- rlnorm(n=nyears, mean=log(min(catch(ple4)/10)), sd=0.1)
control=fwdControl(data.frame(year=years, quant="catch", value=catch_val))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
catch_out <- c(catch(res)[,ac(years)])
# Test to within tolerance of 1.5e-4
expect_equal(catch_out, catch_val, tolerance=1.5e-4)
})
test_that("Catch target, multiple iters",{
data(ple4)
niters <- round(runif(1, min=10, max=50))
ple4p <- propagate(ple4, niters)
# Spoof up stock.n, landings.wt, landings.n and harvest
harvest(ple4p)[] <- rlnorm(n=prod(dim(harvest(ple4p))), mean=log(c(harvest(ple4p))), sd=0.001)
stock.n(ple4p)[] <- rlnorm(n=prod(dim(stock.n(ple4p))), mean=log(c(stock.n(ple4p))), sd=0.001)
landings.wt(ple4p)[] <- rlnorm(n=prod(dim(landings.wt(ple4p))), mean=log(c(landings.wt(ple4p))), sd=0.001)
landings.n(ple4p)[] <- rlnorm(n=prod(dim(landings.n(ple4p))), mean=log(c(landings.n(ple4p))), sd=0.001)
year_range <- range(ple4p)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_val <- rlnorm(n=nyears*niters, mean=log(min(catch(ple4p)/10)), sd=0.1)
control=fwdControl(data.frame(year=years, quant="catch", value=0), iters=niters)
control@iters[,"value",] <- catch_val
res <- fwd(ple4p, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
catch_out <- c(catch(res)[,ac(years)])
# Test to within tolerance of 1.5e-8
expect_equal(catch_out, catch_val, tol=1e-6)
})
test_that("Fbar target, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
f_val <- rlnorm(n=nyears, mean=log(min(fbar(ple4)/10)), sd=0.1)
control=fwdControl(data.frame(year=years, quant="fbar", value=f_val))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
f_out <- c(fbar(res)[,ac(years)])
# Test to within tolerance of 1.5e-8
expect_equal(f_out, f_val)
})
test_that("Fbar target, multiple iters",{
data(ple4)
niters <- round(runif(1, min=10, max=50))
ple4p <- propagate(ple4, niters)
# Spoof up stock.n, landings.wt, landings.n and harvest
harvest(ple4p)[] <- rlnorm(n=prod(dim(harvest(ple4p))), mean=log(c(harvest(ple4p))), sd=0.001)
stock.n(ple4p)[] <- rlnorm(n=prod(dim(stock.n(ple4p))), mean=log(c(stock.n(ple4p))), sd=0.001)
landings.wt(ple4p)[] <- rlnorm(n=prod(dim(landings.wt(ple4p))), mean=log(c(landings.wt(ple4p))), sd=0.001)
landings.n(ple4p)[] <- rlnorm(n=prod(dim(landings.n(ple4p))), mean=log(c(landings.n(ple4p))), sd=0.001)
year_range <- range(ple4p)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
f_val <- rlnorm(n=nyears*niters, mean=log(min(fbar(ple4p)/10)), sd=0.1)
control=fwdControl(data.frame(year=years, quant="fbar", value=0), iters=niters)
control@iters[,"value",] <- f_val
res <- fwd(ple4p, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
f_out <- c(fbar(res)[,ac(years)])
# Test to within tolerance of 1.5e-8
expect_equal(f_out, f_val)
})
test_that("Catch target with min limit, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_val <- rlnorm(n=nyears, mean=log(min(catch(ple4)/10)), sd=0.1)
catch_lim <- mean(catch_val)
# First project without limit
control=fwdControl(list(year=years, quant="catch", value=catch_val))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean", params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# Get resulting fbar
fout <- fbar(res)[,ac(years)]
limit_year <- sample(years,1)
flim <- c(fout[,ac(limit_year)] * 1.5)
# New projection
control=fwdControl(list(year=years, quant="catch", value=catch_val),
list(year=limit_year, quant="fbar", min=flim))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean", params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# All years apart from lim_year met catch OK
catch_out <- catch(res)[,ac(years)]
non_lim_years <- !(years %in% limit_year)
expect_equal(c(catch_out[,ac(years[non_lim_years])]), catch_val[non_lim_years],
tolerance=0.0005)
# lim year has f = flim
expect_equal(c(fbar(res)[,ac(limit_year)]), flim)
})
test_that("Catch target with max limit, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_val <- rlnorm(n=nyears, mean=log(min(catch(ple4)/10)), sd=0.1)
catch_lim <- mean(catch_val)
# First project without limit
control=fwdControl(list(year=years, quant="catch", value=catch_val))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean", params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# Get resulting fbar
fout <- fbar(res)[,ac(years)]
limit_year <- sample(years,1)
flim <- c(fout[,ac(limit_year)] * 0.5)
# New projection
control=fwdControl(list(year=years, quant="catch", value=catch_val),
list(year=limit_year, quant="fbar", max=flim))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean", params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# All years apart from lim_year met catch OK
catch_out <- catch(res)[,ac(years)]
non_lim_years <- !(years %in% limit_year)
expect_true(all(c( catch_out[,ac(years[non_lim_years])]) -
catch_val[non_lim_years] < 1e-8))
# lim year has f = flim
expect_equal(c(fbar(res)[,ac(limit_year)]), flim)
})
test_that("Catch relative target, single iter",{
data(ple4)
year_range <- range(ple4)[c("minyear","maxyear")]
min_year <- round(runif(1, min=year_range[1]+1, max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
min_year[1] <- min_year[1] + 1
nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
years <- min_year:(min_year+nyears-1)
catch_rel_val <- runif(nyears, min=0.5, max=0.9)
control=fwdControl(data.frame(year=years, quant="catch", value=catch_rel_val, relYear=years-1))
res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
catch_out <- c(catch(res)[,ac(years)])
catch_trg <- c(catch(ple4)[,ac(years[1]-1)]) * cumprod(catch_rel_val)
# Test to within tolerance of 1.5e-8
expect_equal(catch_out, catch_trg, tol=1.5e-6)
})
# DEBUG
# test_that("Fbar relative target, single iter",{
# data(ple4)
# year_range <- range(ple4)[c("minyear","maxyear")]
# min_year <- round(runif(1, min=year_range[1], max=(year_range[1] + round((year_range[2] - year_range[1])/2))))
# min_year[1] <- min_year[1] + 1
# nyears <- round(runif(n=1,min=3,max=(round((year_range[2] - year_range[1])/2)-5)))
# years <- min_year:(min_year+nyears-1)
# f_rel_val <- runif(nyears, min=0.5, max=0.9)
# control <- fwdControl(data.frame(year=years, quant="fbar",
# value=f_rel_val, relYear=years-1))
# res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
# params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))))
# f_out <- c(fbar(res)[,ac(years)])
# f_trg <- c(fbar(ple4)[,ac(years[1]-1)]) * cumprod(f_rel_val)
# # Test to within tolerance of 1.5e-8
# expect_equal(f_out, f_trg)
# })
#
#---------------------
test_that("Tests from Running Medium Term Forecasts with FLasher tutorial",{
data(ple4)
ple4_mtf <- stf(ple4, nyears = 10)
ple4_sr <- fmle(as.FLSR(ple4, model="bevholt"), control=list(trace=0))
# Simple F
f_status_quo <- mean(fbar(ple4)[,as.character(2005:2008)])
ctrl_f <- fwdControl(list(year=2009:2018, quant="f", value=f_status_quo))
ple4_f_sq <- fwd(ple4_mtf, control = ctrl_f, sr = ple4_sr)
expect_equal(c(fbar(ple4_f_sq)[,ac(2009:2018)]), rep(f_status_quo,10))
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_f_sq)[,ac(2008:2017)])), c(rec(ple4_f_sq)[,ac(2009:2018)]))
# Decreasing catch
future_catch <- c(catch(ple4)[,"2008"]) * 0.9^(1:10)
ctrl_catch <- fwdControl(list(year=2009:2018, quant = "catch", value=future_catch))
ple4_catch <- fwd(ple4_mtf, control = ctrl_catch, sr = ple4_sr)
expect_equal(c(catch(ple4_catch)[,ac(2009:2018)]), future_catch)
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_catch)[,ac(2008:2017)])), c(rec(ple4_catch)[,ac(2009:2018)]))
# SSB end
final_ssb <- 100000
ctrl_ssb <- fwdControl(list(year=2009, quant = "ssb_end", value=final_ssb))
ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
survivors <- stock.n(ple4_ssb) * exp(-harvest(ple4_ssb) - m(ple4_ssb))
expect_equal(c(quantSums((survivors * stock.wt(ple4_ssb) * mat(ple4_ssb))[,ac(2009)])), final_ssb)
# SSB at spawn
spawn_ssb <- 200000
ctrl_ssb <- fwdControl(list(year=2009, quant = "ssb_spawn", value=spawn_ssb))
expect_warning(fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr))
m.spwn(ple4_mtf)[,ac(2009)] <- 0.5
harvest.spwn(ple4_mtf)[,ac(2009)] <- 0.5
ctrl_ssb <- fwdControl(data.frame(year=2009, quant = "ssb_spawn", value=spawn_ssb))
ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
expect_equal(c(ssb(ple4_ssb)[,ac(2009)]), spawn_ssb)
# SRP
srp <- 200000
ctrl_ssb <- fwdControl(data.frame(year=2009, quant = "srp", value=srp))
ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
expect_equal(c(ssb(ple4_ssb)[,ac(2009)]), srp)
# FLash-like
m.spwn(ple4_mtf)[,ac(2009)] <- 0.5
harvest.spwn(ple4_mtf)[,ac(2009)] <- 0.5
flash_ssb <- 150000
ctrl_ssb <- fwdControl(data.frame(year=2009, quant = "ssb_flash", value=flash_ssb))
ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
expect_equal(flash_ssb, c(ssb(ple4_ssb)[,ac(2009)]))
# m.spwn(ple4_mtf)[,ac(2009)] <- 0.0
# harvest.spwn(ple4_mtf)[,ac(2009)] <- 0.0
# flash_ssb <- 150000
# ctrl_ssb <- fwdControl(data.frame(year=2009, quant = "ssb_flash",
# value=flash_ssb))
# ple4_ssb <- fwd(ple4_mtf, control=ctrl_ssb, sr = ple4_sr)
# expect_equal(c(ssb(ple4_ssb)[,ac(2010)]), flash_ssb)
#
# Longer SSB
m.spwn(ple4_mtf)[,ac(dims(ple4)$maxyear + 1)] <- 0.0
harvest.spwn(ple4_mtf)[,ac(dims(ple4)$maxyear + 1)] <- 0.0
future_ssb <- c(ssb(ple4)[,'2017']) * 0.8
ctrl_ssb <- fwdControl(data.frame(year=seq(2018, 2027), quant = "ssb_flash", value=future_ssb))
expect_warning(ple4_ssb <- fwd(ple4_mtf, control = ctrl_ssb, sr = ple4_sr))
expect_equal(c(ssb(ple4_ssb)[,ac(2019:2027)]), rep(future_ssb, 9))
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_ssb)[,ac(2019:2026)])), c(rec(ple4_ssb)[,ac(2020:2027)]))
# Relative catch
ctrl_rel_catch <- fwdControl(data.frame(
year = 2009:2018, quant = "catch", value = 0.9, relYear = 2008:2017))
ple4_rel_catch <- fwd(ple4_mtf, control = ctrl_rel_catch, sr = ple4_sr)
expect_equal(c(catch(ple4_rel_catch)[,ac(2009:2018)] / catch(ple4_rel_catch)[,ac(2008:2017)]),
rep(0.9, 10))
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_rel_catch)[,ac(2008:2017)])),
c(rec(ple4_rel_catch)[,ac(2009:2018)]))
# Min Max
f01 <- 0.1
min_catch <- mean(catch(ple4_mtf)[,as.character(2006:2008)])
ctrl_min_catch <- fwdControl(
list(year=2009:2018, quant="f", value=f01), list(year=2009:2018, quant="catch", min=min_catch))
ple4_min_catch <- fwd(ple4_mtf, control = ctrl_min_catch, sr = ple4_sr)
ple4c <- catch(ple4_min_catch)[,ac(2009:2018)]
# Expect a solving tolerance of better than 1.5e-8
expect_true(all((ple4c - min_catch) < 1.5e-8))
# Relative targets and bounds
current_fbar <- c(fbar(ple4)[,"2008"])
f_target <- c(seq(from = current_fbar, to = f01, length = 8)[-1], rep(f01, 3))
rel_catch_bound <- 0.10
ctrl_rel_min_max_catch <- fwdControl(
list(year=2009:2018, quant="f", value=f_target),
list(year=2009:2018, quant="catch", relYear=2008:2017, max=1+rel_catch_bound, min=1-rel_catch_bound))
recovery<-fwd(ple4_mtf, control=ctrl_rel_min_max_catch, sr=ple4_sr)
ple4c <- catch(recovery)[,ac(2009:2018)]
relcatch <- catch(recovery)[,ac(2009:2018)] / catch(recovery)[,ac(2008:2017)]
# All relative catches within bounds within tolerance
expect_true(all(relcatch > (1-rel_catch_bound - 1e6) & relcatch < (1+rel_catch_bound + 1e6)))
# Prepare iterations
niters <- 20
ple4_mtf <- stf(ple4, nyears = 10)
ple4_mtf <- propagate(ple4_mtf, niters)
# deviances
rec_deviances <- FLQuant(NA, dimnames = list(year=2009:2018, iter=1:niters))
sample_years <- sample(dimnames(residuals(ple4_sr))$year, niters * 10, replace = TRUE)
rec_deviances[] <- exp(residuals(ple4_sr)[,sample_years])
ple4_stoch_rec <- fwd(ple4_mtf, control = ctrl_catch, sr = ple4_sr, deviances = rec_deviances)
expect_equal(c(catch(ple4_stoch_rec)[,ac(2009:2018)]),rep(future_catch,niters))
expect_equal(c(predict(ple4_sr,ssb=ssb(ple4_stoch_rec)[,ac(2008:2017)]) %*% rec_deviances), c(rec(ple4_stoch_rec)[,ac(2009:2018)]))
# Stochastic target
stoch_catch <- rlnorm(10*niters, meanlog=log(future_catch), sdlog=0.3)
ctrl_catch_iters <- fwdControl(list(year=2009:2018, quant="catch", value=stoch_catch))
ple4_catch_iters <- fwd(ple4_mtf, control=ctrl_catch_iters, sr = ple4_sr)
expect_equal(c(catch(ple4_catch_iters)[,ac(2009:2018)]), stoch_catch)
# Iterations in recruitment parameters
sr_iters <- FLPar(NA, dimnames=list(params=c("a","b"), iter=1:niters))
aiters <- rlnorm(niters, meanlog=log(params(ple4_sr)["a"]), sdlog=0.5)
biters <- rlnorm(niters, meanlog=log(params(ple4_sr)["b"]), sdlog=0.01)
sr_iters["a"] <- aiters
sr_iters["b"] <- biters
ple4_sr_iters <- fwd(ple4_mtf, control=ctrl_catch, sr = list(model="bevholt", params=sr_iters))
expect_equal(c(catch(ple4_sr_iters)[,ac(2009:2018)]), rep(future_catch, niters))
srmod <- predictModel(model=bevholt()$model, params=sr_iters)
expect_equal(c(predict(srmod, ssb=ssb(ple4_sr_iters)[,ac(2008:2017)])), c(rec(ple4_sr_iters)[,ac(2009:2018)]))
# Stochastic targets and recruitment
ple4_iters <- fwd(ple4_mtf, control=ctrl_catch_iters, sr = list(model="bevholt", params=sr_iters), deviances = rec_deviances)
expect_equal(c(catch(ple4_iters)[,ac(2009:2018)]), stoch_catch)
expect_equal(c(predict(srmod, ssb=ssb(ple4_iters)[,ac(2008:2017)]) %*% rec_deviances), c(rec(ple4_iters)[,ac(2009:2018)]))
})
test_that("Catch target, maxf limiting",{
data(ple4)
control=fwdControl(year=1990, quant="catch", value=150000)
res <- fwd(ple4, control=control, sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2006:2008)]))), maxF=0.1)
f_out <- c(fbar(res)[,"1990"])
expect_equal(f_out, 0.1)
})
test_that("control as FLQuants works", {
data(ple4)
# catch
res <- fwd(ple4, catch=catch(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(catch(ple4)[,ac(2000:2008)], catch(res)[,ac(2000:2008)])
# landings
res <- fwd(ple4, landings=landings(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(landings(ple4)[,ac(2000:2008)], landings(res)[,ac(2000:2008)])
# discards
res <- fwd(ple4, discards=discards(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(discards(ple4)[,ac(2000:2008)], discards(res)[,ac(2000:2008)])
# fbar
res <- fwd(ple4, fbar=fbar(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(fbar(ple4)[,ac(2000:2008)], fbar(res)[,ac(2000:2008)])
# f
res <- fwd(ple4, f=fbar(ple4)[,ac(2000:2008)],
sr=predictModel(model="geomean",
params=FLPar(a=yearMeans(rec(ple4)[, ac(2000:2008)]))))
expect_equal(fbar(ple4)[,ac(2000:2008)], fbar(res)[,ac(2000:2008)])
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.