library(unittest)
library(gadget3)
prey_a <- g3_stock('prey_a', seq(1, 10)) |> g3s_age(1,3)
prey_b <- g3_stock('prey_b', seq(1, 10)) |> g3s_age(1,3)
pred_a <- g3_stock('pred_a', seq(50, 80, by = 10)) |> g3s_age(0, 10)
otherfood <- g3_stock('otherfood', 0)
pred_a_catch_obs <- expand.grid(
year = 2000:2005,
length = c(0,5,10),
stock = c('prey_a', 'prey_b', 'otherfood'),
predator = c('pred_a'), # NB: Not essential if there's only one predator, only for testing
predator_length = c(50,70),
predator_age = c("[0,5)", "[6,10)"), # ((
number = 0 )
pred_a_preypref_obs <- expand.grid(
year = 2000:2005,
predator_length = c(50,70),
stock = c('prey_a', 'prey_b', 'otherfood'),
number = 0 )
pred_a_sizepref_obs <- expand.grid(
year = 2000:2005,
predator_age = c("[0,5)", "[6,10)"), # ((
length = seq(1, 10), # i.e. the union of prey_a / prey_b lengths
number = 0 )
actions <- list(
g3a_time(2000, 2005, step_lengths = c(6,6), project_years = 0),
g3a_age(prey_a),
g3a_age(prey_b),
g3a_age(pred_a),
g3a_initialconditions(prey_a, ~1e10 + 0 * prey_a__midlen, ~100),
g3a_initialconditions(prey_b, ~2e10 + 0 * prey_b__midlen, ~200),
g3a_initialconditions(pred_a, ~1e5 + 0 * pred_a__midlen, ~1000),
g3a_otherfood(otherfood, ~1e10 + 0 * otherfood__midlen, 100),
g3a_predate(
pred_a,
list(prey_a, prey_b, otherfood),
suitabilities = list(prey_a = 0.1, prey_b = 0.1, otherfood = 0.01),
catchability_f = g3a_predate_catchability_predator(
temperature = g3_parameterized('temp', value = 0, by_year = TRUE, optimise = FALSE)) ),
g3l_understocking(list(prey_a, prey_b)),
g3l_catchdistribution(
'pred_a_catch',
pred_a_catch_obs,
predators = list(pred_a),
stocks = list(prey_a, prey_b, otherfood),
function_f = g3l_distribution_sumofsquares(),
nll_breakdown = TRUE,
report = TRUE ),
g3l_catchdistribution(
'pred_a_preypref',
pred_a_preypref_obs,
predators = list(pred_a),
stocks = list(prey_a, prey_b, otherfood),
function_f = g3l_distribution_sumofsquares(),
nll_breakdown = TRUE,
report = TRUE ),
g3l_catchdistribution(
'pred_a_sizepref',
pred_a_sizepref_obs,
predators = list(pred_a),
stocks = list(prey_a, prey_b), # NB: otherfood missing
function_f = g3l_distribution_sumofsquares(),
nll_breakdown = TRUE,
report = TRUE ),
# NB: Dummy parameter so model will compile in TMB
~{nll <- nll + g3_param("x", value = 0)} )
actions <- c(actions, list(
g3a_report_history(actions, var_re = "__feedinglevel$|__totalsuit$"),
g3a_report_detail(actions) ))
model_fn <- g3_to_r(actions)
model_cpp <- g3_to_tmb(actions)
ok_group("Default params") ########
params <- attr(model_fn, 'parameter_template')
result <- model_fn(params)
def_r <- attributes(result)
gadget3:::ut_tmb_r_compare2(model_fn, model_cpp, params)
######## Default params
ok_group("Energy content") ########
params <- attr(model_fn, 'parameter_template')
params$prey_a.energycontent <- 2
params$prey_b.energycontent <- 1
result <- model_fn(params)
r <- attributes(result)
ok(all.equal(
def_r$detail_prey_a_pred_a__suit[] * 2,
r$detail_prey_a_pred_a__suit[],
tolerance = 1e-6), "detail_prey_a_pred_a__suit: Doubled thanks to energycontent")
gadget3:::ut_tmb_r_compare2(model_fn, model_cpp, params)
######## Energy content
ok_group("consumption.m3") ########
ok(ut_cmp_equal(
diff(def_r$detail_prey_a_pred_a__cons[1,1,,1,1]),
c("60:70" = 0, "70:80" = 0, "80:Inf" = 0),
tolerance = 1e-8 ), "def_r$detail_prey_a_pred_a__cons: No difference in pred.length consumption")
params <- attr(model_fn, 'parameter_template')
params$pred_a.consumption.m3 <- 2
params$pred_a.halffeeding <- 1e15 # Needs to be proportion higher than total_predsuit to see any result
result <- model_fn(params) ; r <- attributes(result)
ok(all.equal(
def_r$detail_prey_a_pred_a__suit,
r$detail_prey_a_pred_a__suit,
tolerance = 1e-3 ), "r$detail_prey_a_pred_a__suit: consumption.m3 results in approximately equal __suit")
ok(ut_cmp_equal(
r$detail_prey_a_pred_a__cons[1,1,2,1,1],
r$detail_prey_a_pred_a__cons[1,1,1,1,1] / 55^2 * 65^2,
tolerance = 1e-8 ), "r$detail_prey_a_pred_a__cons: Jump in consumption 55 -> 65")
ok(ut_cmp_equal(
r$detail_prey_a_pred_a__cons[1,1,3,1,1],
r$detail_prey_a_pred_a__cons[1,1,2,1,1] / 65^2 * 75^2,
tolerance = 1e-8 ), "r$detail_prey_a_pred_a__cons: Jump in consumption 65 -> 75")
ok(ut_cmp_equal(
r$detail_prey_a_pred_a__cons[1,1,4,1,1],
r$detail_prey_a_pred_a__cons[1,1,3,1,1] / 75^2 * 85^2,
tolerance = 1e-8 ), "r$detail_prey_a_pred_a__cons: Jump in consumption 75 -> 85")
ok(ut_cmp_identical(
dimnames(r$cdist_sumofsquares_pred_a_catch_model__num)$predator,
c("pred_a") ), "cdist_sumofsquares_pred_a_catch_model__num: Single predator dimension for our one predator")
ok(gadget3:::ut_cmp_array(drop(r$cdist_sumofsquares_pred_a_catch_model__num), '
length stock predator_length predator_age time Freq
1 0:5 prey_a 50:70 0:4 2000 84464.2913
2 5:10 prey_a 50:70 0:4 2000 105580.3641
3 10:Inf prey_a 50:70 0:4 2000 21116.0728
4 0:5 prey_b 50:70 0:4 2000 168928.5826
5 5:10 prey_b 50:70 0:4 2000 211160.7283
6 10:Inf prey_b 50:70 0:4 2000 42232.1457
7 0:5 otherfood 50:70 0:4 2000 703.8698
# NB: All otherfood are length 0
8 5:10 otherfood 50:70 0:4 2000 0.0000
9 10:Inf otherfood 50:70 0:4 2000 0.0000
10 0:5 prey_a 70:Inf 0:4 2000 149705.6750
11 5:10 prey_a 70:Inf 0:4 2000 187132.0937
12 10:Inf prey_a 70:Inf 0:4 2000 37426.4187
13 0:5 prey_b 70:Inf 0:4 2000 299411.3499
14 5:10 prey_b 70:Inf 0:4 2000 374264.1874
15 10:Inf prey_b 70:Inf 0:4 2000 74852.8375
16 0:5 otherfood 70:Inf 0:4 2000 1247.5486
17 5:10 otherfood 70:Inf 0:4 2000 0.0000
18 10:Inf otherfood 70:Inf 0:4 2000 0.0000
19 0:5 prey_a 50:70 6:9 2000 67571.4331
20 5:10 prey_a 50:70 6:9 2000 84464.2913
21 10:Inf prey_a 50:70 6:9 2000 16892.8583
22 0:5 prey_b 50:70 6:9 2000 135142.8661
23 5:10 prey_b 50:70 6:9 2000 168928.5826
24 10:Inf prey_b 50:70 6:9 2000 33785.7165
25 0:5 otherfood 50:70 6:9 2000 563.0959
26 5:10 otherfood 50:70 6:9 2000 0.0000
27 10:Inf otherfood 50:70 6:9 2000 0.0000
28 0:5 prey_a 70:Inf 6:9 2000 119764.5400
29 5:10 prey_a 70:Inf 6:9 2000 149705.6750
30 10:Inf prey_a 70:Inf 6:9 2000 29941.1350
31 0:5 prey_b 70:Inf 6:9 2000 239529.0799
32 5:10 prey_b 70:Inf 6:9 2000 299411.3499
33 10:Inf prey_b 70:Inf 6:9 2000 59882.2700
34 0:5 otherfood 70:Inf 6:9 2000 998.0389
35 5:10 otherfood 70:Inf 6:9 2000 0.0000
36 10:Inf otherfood 70:Inf 6:9 2000 0.0000
37 0:5 prey_a 50:70 0:4 2001 67571.1514
38 5:10 prey_a 50:70 0:4 2001 84463.9393
39 10:Inf prey_a 50:70 0:4 2001 16892.7879
40 0:5 prey_b 50:70 0:4 2001 135142.3028
41 5:10 prey_b 50:70 0:4 2001 168927.8785
42 10:Inf prey_b 50:70 0:4 2001 33785.5757
43 0:5 otherfood 50:70 0:4 2001 563.0960
44 5:10 otherfood 50:70 0:4 2001 0.0000
45 10:Inf otherfood 50:70 0:4 2001 0.0000
46 0:5 prey_a 70:Inf 0:4 2001 119764.0408
47 5:10 prey_a 70:Inf 0:4 2001 149705.0510
48 10:Inf prey_a 70:Inf 0:4 2001 29941.0102
49 0:5 prey_b 70:Inf 0:4 2001 239528.0816
50 5:10 prey_b 70:Inf 0:4 2001 299410.1019
51 10:Inf prey_b 70:Inf 0:4 2001 59882.0204
52 0:5 otherfood 70:Inf 0:4 2001 998.0390
53 5:10 otherfood 70:Inf 0:4 2001 0.0000
54 10:Inf otherfood 70:Inf 0:4 2001 0.0000
55 0:5 prey_a 50:70 6:9 2001 67571.1514
56 5:10 prey_a 50:70 6:9 2001 84463.9393
57 10:Inf prey_a 50:70 6:9 2001 16892.7879
58 0:5 prey_b 50:70 6:9 2001 135142.3028
59 5:10 prey_b 50:70 6:9 2001 168927.8785
60 10:Inf prey_b 50:70 6:9 2001 33785.5757
61 0:5 otherfood 50:70 6:9 2001 563.0960
62 5:10 otherfood 50:70 6:9 2001 0.0000
63 10:Inf otherfood 50:70 6:9 2001 0.0000
64 0:5 prey_a 70:Inf 6:9 2001 119764.0408
65 5:10 prey_a 70:Inf 6:9 2001 149705.0510
66 10:Inf prey_a 70:Inf 6:9 2001 29941.0102
67 0:5 prey_b 70:Inf 6:9 2001 239528.0816
68 5:10 prey_b 70:Inf 6:9 2001 299410.1019
69 10:Inf prey_b 70:Inf 6:9 2001 59882.0204
70 0:5 otherfood 70:Inf 6:9 2001 998.0390
71 5:10 otherfood 70:Inf 6:9 2001 0.0000
72 10:Inf otherfood 70:Inf 6:9 2001 0.0000
73 0:5 prey_a 50:70 0:4 2002 50678.1523
74 5:10 prey_a 50:70 0:4 2002 63347.6904
75 10:Inf prey_a 50:70 0:4 2002 12669.5381
76 0:5 prey_b 50:70 0:4 2002 101356.3047
77 5:10 prey_b 50:70 0:4 2002 126695.3808
78 10:Inf prey_b 50:70 0:4 2002 25339.0762
79 0:5 otherfood 50:70 0:4 2002 422.3220
80 5:10 otherfood 50:70 0:4 2002 0.0000
81 10:Inf otherfood 50:70 0:4 2002 0.0000
82 0:5 prey_a 70:Inf 0:4 2002 89822.6562
83 5:10 prey_a 70:Inf 0:4 2002 112278.3202
84 10:Inf prey_a 70:Inf 0:4 2002 22455.6640
85 0:5 prey_b 70:Inf 0:4 2002 179645.3124
86 5:10 prey_b 70:Inf 0:4 2002 224556.6405
87 10:Inf prey_b 70:Inf 0:4 2002 44911.3281
88 0:5 otherfood 70:Inf 0:4 2002 748.5294
89 5:10 otherfood 70:Inf 0:4 2002 0.0000
90 10:Inf otherfood 70:Inf 0:4 2002 0.0000
91 0:5 prey_a 50:70 6:9 2002 67570.8698
92 5:10 prey_a 50:70 6:9 2002 84463.5872
93 10:Inf prey_a 50:70 6:9 2002 16892.7174
94 0:5 prey_b 50:70 6:9 2002 135141.7395
95 5:10 prey_b 50:70 6:9 2002 168927.1744
96 10:Inf prey_b 50:70 6:9 2002 33785.4349
97 0:5 otherfood 50:70 6:9 2002 563.0960
98 5:10 otherfood 50:70 6:9 2002 0.0000
99 10:Inf otherfood 50:70 6:9 2002 0.0000
100 0:5 prey_a 70:Inf 6:9 2002 119763.5416
101 5:10 prey_a 70:Inf 6:9 2002 149704.4270
102 10:Inf prey_a 70:Inf 6:9 2002 29940.8854
103 0:5 prey_b 70:Inf 6:9 2002 239527.0832
104 5:10 prey_b 70:Inf 6:9 2002 299408.8540
105 10:Inf prey_b 70:Inf 6:9 2002 59881.7708
106 0:5 otherfood 70:Inf 6:9 2002 998.0392
107 5:10 otherfood 70:Inf 6:9 2002 0.0000
108 10:Inf otherfood 70:Inf 6:9 2002 0.0000
109 0:5 prey_a 50:70 0:4 2003 33785.2941
110 5:10 prey_a 50:70 0:4 2003 42231.6176
111 10:Inf prey_a 50:70 0:4 2003 8446.3235
112 0:5 prey_b 50:70 0:4 2003 67570.5881
113 5:10 prey_b 50:70 0:4 2003 84463.2352
114 10:Inf prey_b 50:70 0:4 2003 16892.6470
115 0:5 otherfood 50:70 0:4 2003 281.5480
116 5:10 otherfood 50:70 0:4 2003 0.0000
117 10:Inf otherfood 50:70 0:4 2003 0.0000
118 0:5 prey_a 70:Inf 0:4 2003 59881.5212
119 5:10 prey_a 70:Inf 0:4 2003 74851.9015
120 10:Inf prey_a 70:Inf 0:4 2003 14970.3803
121 0:5 prey_b 70:Inf 0:4 2003 119763.0424
122 5:10 prey_b 70:Inf 0:4 2003 149703.8030
123 10:Inf prey_b 70:Inf 0:4 2003 29940.7606
124 0:5 otherfood 70:Inf 0:4 2003 499.0196
125 5:10 otherfood 70:Inf 0:4 2003 0.0000
126 10:Inf otherfood 70:Inf 0:4 2003 0.0000
127 0:5 prey_a 50:70 6:9 2003 67570.5881
128 5:10 prey_a 50:70 6:9 2003 84463.2352
129 10:Inf prey_a 50:70 6:9 2003 16892.6470
130 0:5 prey_b 50:70 6:9 2003 135141.1763
131 5:10 prey_b 50:70 6:9 2003 168926.4703
132 10:Inf prey_b 50:70 6:9 2003 33785.2941
133 0:5 otherfood 50:70 6:9 2003 563.0961
134 5:10 otherfood 50:70 6:9 2003 0.0000
135 10:Inf otherfood 50:70 6:9 2003 0.0000
136 0:5 prey_a 70:Inf 6:9 2003 119763.0424
137 5:10 prey_a 70:Inf 6:9 2003 149703.8030
138 10:Inf prey_a 70:Inf 6:9 2003 29940.7606
139 0:5 prey_b 70:Inf 6:9 2003 239526.0848
140 5:10 prey_b 70:Inf 6:9 2003 299407.6060
141 10:Inf prey_b 70:Inf 6:9 2003 59881.5212
142 0:5 otherfood 70:Inf 6:9 2003 998.0393
143 5:10 otherfood 70:Inf 6:9 2003 0.0000
144 10:Inf otherfood 70:Inf 6:9 2003 0.0000
145 0:5 prey_a 50:70 0:4 2004 16892.5766
146 5:10 prey_a 50:70 0:4 2004 21115.7208
147 10:Inf prey_a 50:70 0:4 2004 4223.1442
148 0:5 prey_b 50:70 0:4 2004 33785.1532
149 5:10 prey_b 50:70 0:4 2004 42231.4416
150 10:Inf prey_b 50:70 0:4 2004 8446.2883
151 0:5 otherfood 50:70 0:4 2004 140.7740
152 5:10 otherfood 50:70 0:4 2004 0.0000
153 10:Inf otherfood 50:70 0:4 2004 0.0000
154 0:5 prey_a 70:Inf 0:4 2004 29940.6358
155 5:10 prey_a 70:Inf 0:4 2004 37425.7948
156 10:Inf prey_a 70:Inf 0:4 2004 7485.1590
157 0:5 prey_b 70:Inf 0:4 2004 59881.2716
158 5:10 prey_b 70:Inf 0:4 2004 74851.5895
159 10:Inf prey_b 70:Inf 0:4 2004 14970.3179
160 0:5 otherfood 70:Inf 0:4 2004 249.5099
161 5:10 otherfood 70:Inf 0:4 2004 0.0000
162 10:Inf otherfood 70:Inf 0:4 2004 0.0000
163 0:5 prey_a 50:70 6:9 2004 67570.3065
164 5:10 prey_a 50:70 6:9 2004 84462.8831
165 10:Inf prey_a 50:70 6:9 2004 16892.5766
166 0:5 prey_b 50:70 6:9 2004 135140.6130
167 5:10 prey_b 50:70 6:9 2004 168925.7662
168 10:Inf prey_b 50:70 6:9 2004 33785.1532
169 0:5 otherfood 50:70 6:9 2004 563.0962
170 5:10 otherfood 50:70 6:9 2004 0.0000
171 10:Inf otherfood 50:70 6:9 2004 0.0000
172 0:5 prey_a 70:Inf 6:9 2004 119762.5432
173 5:10 prey_a 70:Inf 6:9 2004 149703.1790
174 10:Inf prey_a 70:Inf 6:9 2004 29940.6358
175 0:5 prey_b 70:Inf 6:9 2004 239525.0865
176 5:10 prey_b 70:Inf 6:9 2004 299406.3581
177 10:Inf prey_b 70:Inf 6:9 2004 59881.2716
178 0:5 otherfood 70:Inf 6:9 2004 998.0394
179 5:10 otherfood 70:Inf 6:9 2004 0.0000
180 10:Inf otherfood 70:Inf 6:9 2004 0.0000
# NB: Ages 0..4 have all grown up by this point, none left
181 0:5 prey_a 50:70 0:4 2005 0.0000
182 5:10 prey_a 50:70 0:4 2005 0.0000
183 10:Inf prey_a 50:70 0:4 2005 0.0000
184 0:5 prey_b 50:70 0:4 2005 0.0000
185 5:10 prey_b 50:70 0:4 2005 0.0000
186 10:Inf prey_b 50:70 0:4 2005 0.0000
187 0:5 otherfood 50:70 0:4 2005 0.0000
188 5:10 otherfood 50:70 0:4 2005 0.0000
189 10:Inf otherfood 50:70 0:4 2005 0.0000
190 0:5 prey_a 70:Inf 0:4 2005 0.0000
191 5:10 prey_a 70:Inf 0:4 2005 0.0000
192 10:Inf prey_a 70:Inf 0:4 2005 0.0000
193 0:5 prey_b 70:Inf 0:4 2005 0.0000
194 5:10 prey_b 70:Inf 0:4 2005 0.0000
195 10:Inf prey_b 70:Inf 0:4 2005 0.0000
196 0:5 otherfood 70:Inf 0:4 2005 0.0000
197 5:10 otherfood 70:Inf 0:4 2005 0.0000
198 10:Inf otherfood 70:Inf 0:4 2005 0.0000
199 0:5 prey_a 50:70 6:9 2005 67570.0249
200 5:10 prey_a 50:70 6:9 2005 84462.5311
201 10:Inf prey_a 50:70 6:9 2005 16892.5062
202 0:5 prey_b 50:70 6:9 2005 135140.0497
203 5:10 prey_b 50:70 6:9 2005 168925.0621
204 10:Inf prey_b 50:70 6:9 2005 33785.0124
205 0:5 otherfood 50:70 6:9 2005 563.0962
206 5:10 otherfood 50:70 6:9 2005 0.0000
207 10:Inf otherfood 50:70 6:9 2005 0.0000
208 0:5 prey_a 70:Inf 6:9 2005 119762.0440
209 5:10 prey_a 70:Inf 6:9 2005 149702.5551
210 10:Inf prey_a 70:Inf 6:9 2005 29940.5110
211 0:5 prey_b 70:Inf 6:9 2005 239524.0881
212 5:10 prey_b 70:Inf 6:9 2005 299405.1101
213 10:Inf prey_b 70:Inf 6:9 2005 59881.0220
214 0:5 otherfood 70:Inf 6:9 2005 998.0395
215 5:10 otherfood 70:Inf 6:9 2005 0.0000
216 10:Inf otherfood 70:Inf 6:9 2005 0.0000
'), "cdist_sumofsquares_pred_a_catch_model__num: Baseline output")
for (t in dimnames(r$hist_pred_a__feedinglevel)$time) {
ok(all(r$hist_pred_a__feedinglevel[,,time = t] == r$hist_pred_a__feedinglevel[1,1,time = t]), paste0(
"r$hist_pred_a__feedinglevel: Not dependent on length/age", t))
}
ok(ut_cmp_equal(r$hist_pred_a__feedinglevel[1,1,], c(
"2000-01" = 0.0291450651443661,
"2000-02" = 0.0291450044465519,
"2001-01" = 0.0291449437488566,
"2001-02" = 0.0291448830512803,
"2002-01" = 0.0291448223538228,
"2002-02" = 0.0291447616564842,
"2003-01" = 0.0291447009592646,
"2003-02" = 0.0291446402621638,
"2004-01" = 0.029144579565182,
"2004-02" = 0.029144518868319,
"2005-01" = 0.029144458171575,
"2005-02" = 0.0291443974749498,
NULL ), tolerance = 1e-8), "r$hist_pred_a__feedinglevel[1,1,]: Gentrly dropping as predators age out")
ok(ut_cmp_identical(dimnames(r$cdist_sumofsquares_pred_a_sizepref_model__num), list(
length = c("1:2", "2:3", "3:4", "4:5", "5:6", "6:7", "7:8", "8:9", "9:10", "10:Inf"),
predator_age = c("0:4", "6:9"),
time = c("2000", "2001", "2002", "2003", "2004", "2005") )), "r$cdist_sumofsquares_pred_a_sizepref_model__num: expected dimensions")
ok(ut_cmp_identical(dimnames(r$cdist_sumofsquares_pred_a_preypref_model__num), list(
length = "0:Inf",
stock = c("prey_a", "prey_b", "otherfood"),
predator_length = c("50:70", "70:Inf"),
time = c("2000", "2001", "2002", "2003", "2004", "2005") )), "r$cdist_sumofsquares_pred_a_preypref_model__num: expected dimensions")
gadget3:::ut_tmb_r_compare2(model_fn, model_cpp, params)
######## consumption.m3
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.