-
Notifications
You must be signed in to change notification settings - Fork 16
/
wham_v0.cpp
1403 lines (1334 loc) · 65.2 KB
/
wham_v0.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#define TMB_LIB_INIT R_init_wham
#include <TMB.hpp>
#include <iostream>
#include "helper_functions.hpp"
#include "age_comp_osa.hpp"
#include "age_comp_sim.hpp"
template<class Type>
Type objective_function<Type>::operator() ()
{
using namespace density; // necessary to use AR1, SCALE, SEPARABLE
DATA_INTEGER(n_years_catch); //same as n_years_model
DATA_INTEGER(n_years_indices); //same as n_years_model
DATA_INTEGER(n_years_model);
DATA_INTEGER(n_ages);
DATA_INTEGER(n_fleets);
DATA_INTEGER(n_indices);
DATA_INTEGER(n_selblocks);
DATA_IVECTOR(selblock_models); // for each block: 1 = age-specific, 2 = logistic, 3 = double-logistic, 4 = logistic (declining)
DATA_IVECTOR(selblock_models_re); // for each block: 1 = none, 2 = IID, 3 = ar1, 4 = ar1_y, 5 = 2dar1
DATA_IVECTOR(n_selpars);
DATA_IMATRIX(selpars_est); // n_blocks x (n_pars + n_ages), is the selpar estimated in this block?
DATA_IVECTOR(n_selpars_est); // of the selpars, how many are actually estimated (not fixed at 0 or 1)
DATA_IVECTOR(n_years_selblocks); // for each block, number of years the block covers
DATA_IMATRIX(selblock_years); // n_years_model x n_selblocks, = 1 if block covers year, = 0 if not
DATA_IMATRIX(selblock_pointer_fleets);
DATA_IMATRIX(selblock_pointer_indices);
DATA_IVECTOR(age_comp_model_fleets);
DATA_IVECTOR(age_comp_model_indices);
DATA_VECTOR(fracyr_SSB);
DATA_MATRIX(mature);
DATA_IVECTOR(waa_pointer_fleets);
DATA_INTEGER(waa_pointer_totcatch);
DATA_IVECTOR(waa_pointer_indices);
DATA_INTEGER(waa_pointer_ssb);
DATA_INTEGER(waa_pointer_jan1);
DATA_ARRAY(waa);
DATA_MATRIX(agg_catch);
DATA_IMATRIX(use_agg_catch);
DATA_MATRIX(agg_catch_sigma);
DATA_ARRAY(catch_paa); //n_fleets x n_years x n_ages
DATA_IMATRIX(use_catch_paa);
DATA_MATRIX(catch_Neff);
DATA_IVECTOR(units_indices);
DATA_MATRIX(fracyr_indices);
DATA_MATRIX(agg_indices);
DATA_IMATRIX(use_indices);
DATA_MATRIX(agg_index_sigma);
DATA_IVECTOR(units_index_paa);
DATA_ARRAY(index_paa); //n_indices x n_years x n_ages
DATA_IMATRIX(use_index_paa);
DATA_MATRIX(index_Neff);
DATA_VECTOR(q_lower);
DATA_VECTOR(q_upper);
DATA_IVECTOR(use_q_prior);
DATA_VECTOR(logit_q_prior_sigma);
DATA_IVECTOR(use_q_re); //n_indices, 0= no re, >0 = use re
DATA_MATRIX(selpars_lower);
DATA_MATRIX(selpars_upper);
DATA_INTEGER(n_NAA_sigma); // 0 = SCAA, 1 = logR only, 2 = full state-space with shared sig_a for a > 1
DATA_IVECTOR(NAA_sigma_pointers);
DATA_INTEGER(recruit_model);
DATA_INTEGER(n_M_a);
DATA_INTEGER(M_model); // 1: "constant", 2: "age-specific", 3: "weight-at-age"
DATA_INTEGER(N1_model); //0: just age-specific numbers at age, 1: 2 pars: log_N_{1,1}, log_F0, age-structure defined by equilibrium NAA calculations
DATA_INTEGER(M_re_model); // 1 = none, 2 = IID, 3 = ar1_a, 4 = ar1_y, 5 = 2dar1
// M_est and n_M_est are not use
//DATA_IVECTOR(M_est); // Is mean M estimated for each age? If M-at-age, dim = length(n_M_a). If constant or weight-at-age M, dim = 1.
//DATA_INTEGER(n_M_est); // How many mean M pars are estimated?
DATA_INTEGER(use_b_prior);
DATA_IVECTOR(which_F_age); //which age of F to use for full total F for msy/ypr calculations and projections (n_years_model + n_years_proj)
DATA_INTEGER(use_steepness); // which parameterization to use for BH/Ricker S-R, if needed.
DATA_INTEGER(bias_correct_pe); //bias correct lognormal process error?
DATA_INTEGER(bias_correct_oe); //bias correct lognormal observation error?
DATA_IVECTOR(Fbar_ages);
DATA_IVECTOR(simulate_state); //vector (0/1) if 1 then state parameters (NAA, MAA, sel, Ecov, q) in that order) will be simulated.
DATA_IVECTOR(simulate_data); //vector (0/1) if 1 then data type (catch, indices, Ecov obs) will be simulated.
DATA_IVECTOR(simulate_period); //vector (0/1) if 1 then period (model years, projection years) will be simulated.
DATA_SCALAR(percentSPR); // percentage to use for SPR-based reference points. Default = 40.
DATA_SCALAR(percentFXSPR); // percent of F_XSPR to use for calculating catch in projections. For example, GOM cod uses F = 75% F_40%SPR, so percentFXSPR = 75 and percentSPR = 40. Default = 100.
DATA_SCALAR(percentFMSY); // percent of FMSY to use for calculating catch in projections.
DATA_INTEGER(XSPR_R_opt); //1(3): use annual R estimates(predictions) for annual SSB_XSPR, 2(4): use average R estimates(predictions). See next line for years to average over.
DATA_IVECTOR(XSPR_R_avg_yrs); // model year indices (TMB, starts @ 0) to use for averaging recruitment when defining SSB_XSPR (if XSPR_R_opt = 2,4)
DATA_VECTOR(FXSPR_init); // annual initial values to use for newton steps to find FXSPR (n_years_model+n_proj_years)
DATA_VECTOR(FMSY_init); // annual initial values to use for newton steps to find FMSY (n_years_model+n_proj_years)
// data for one-step-ahead (OSA) residuals
DATA_INTEGER(do_osa); //whether to do osa residuals. For efficiency reasons with age comp likelihoods.
DATA_VECTOR(obsvec); // vector of all observations for OSA residuals
DATA_IVECTOR(agesvec); // vector of ages associated with paa observations for OSA residuals. same length as obsvec
DATA_VECTOR_INDICATOR(keep, obsvec); // for OSA residuals
DATA_IMATRIX(keep_C); // indices for catch obs, can loop years/fleets with keep(keep_C(y,f))
DATA_IMATRIX(keep_I);
DATA_IMATRIX(keep_E); // Ecov
DATA_IARRAY(keep_Cpaa);
DATA_IARRAY(keep_Ipaa);
DATA_IVECTOR(do_post_samp); //length = 5, whether to ADREPORT posterior residuals for NAA, M, selectivity, Ecov, q.
// data for environmental covariate(s), Ecov
DATA_INTEGER(n_Ecov); // also = 1 if no Ecov
DATA_INTEGER(n_years_Ecov); // num years in Ecov process model
DATA_IMATRIX(Ecov_use_obs); // all 0 if no Ecov
DATA_MATRIX(Ecov_obs);
//Below is not used anymore
//DATA_IVECTOR(Ecov_lag);
DATA_IVECTOR(Ecov_how); // specific to recruitment effects. 0 = no effect, 1 = controlling, 2 = limiting, 3 = lethal, 4 = masking, 5 = directive
//Below is not used anymore
//DATA_IMATRIX(Ecov_poly); // n_Ecov x 2+n_indices. polynomial order for ecov effects (1 = linear, 2 = quadratic, 3 = cubic, ...)
DATA_IMATRIX(Ecov_where); // n_Ecov x 2+n_indices. 0/1 values with columns corresponding to recruit, mortality, indices in that order
DATA_IVECTOR(Ecov_model); // 0 = no Ecov, 1 = RW, 2 = AR1
//DATA_INTEGER(year1_Ecov); // first year Ecov
//DATA_INTEGER(year1_model); // first year model
DATA_IMATRIX(ind_Ecov_out_start); // n_Ecov x (2 + n_indices) index of Ecov_x to use for Ecov_out (operates on pop model, lagged effects specific the multiple types of effects each Ecov can have)
DATA_IMATRIX(ind_Ecov_out_end); // n_Ecov x (2 + n_indices) index of Ecov_x to use for Ecov_out (operates on pop model, lagged effects specific the multiple types of effects each Ecov can have)
DATA_IVECTOR(Ecov_obs_sigma_opt); // n_Ecov, 1 = given, 2 = estimate 1 value, shared among obs, 3 = estimate for each obs, 4 = estimate for each obs as random effects
DATA_IMATRIX(Ecov_use_re); // 0/1: use Ecov_re? If yes, add to nll. (n_years_Ecov + n_years_proj_Ecov) x n_Ecov
// data for projections
DATA_INTEGER(do_proj); // 1 = yes, 0 = no
DATA_INTEGER(n_years_proj); // number of years to project
DATA_INTEGER(n_years_proj_Ecov); // number of years to project Ecov
DATA_IVECTOR(avg_years_ind); // model year indices (TMB, starts @ 0) to use for averaging MAA, waa, maturity, and F (if use.avgF = TRUE)
DATA_IVECTOR(proj_F_opt); // for each projection year, 1 = last year F (default), 2 = average F, 3 = F at X% SPR, 4 = user-specified F, 5 = calculate F from user-specified catch
DATA_VECTOR(proj_Fcatch); // user-specified F or catch in projection years, only used if proj_F_opt = 4 or 5
DATA_INTEGER(proj_M_opt); // 1 = continue M_re (check for time-varying M_re on R side), 2 = average M (over avg_years_ind)
DATA_SCALAR(logR_mean); // empirical mean recruitment in model years, used for SCAA recruit projections
DATA_SCALAR(logR_sd); // empirical sd recruitment in model years, used for SCAA recruit projections
DATA_VECTOR(F_proj_init); // annual initial values to use for newton steps to find F for use in projections (n_years_proj)
//static brp info
DATA_INTEGER(which_F_age_static); //which age of F to use for full total F for static brps (max of average FAA_tot over avg_years_ind)
DATA_SCALAR(static_FXSPR_init); // initial value to use for newton steps to find FXSPR_static
//DATA_IVECTOR(static_brp_years_sel); // model year indices (TMB, starts @ 0) to use for averaging selectivity for static biological reference points
//DATA_IVECTOR(static_brp_years_mat); // model year indices (TMB, starts @ 0) to use for averaging maturity for static biological reference points
//DATA_IVECTOR(static_brp_years_waa_ssb); // model year indices (TMB, starts @ 0) to use for averaging SSB weight at age for static biological reference points
//DATA_IVECTOR(static_brp_years_waa_catch); // model year indices (TMB, starts @ 0) to use for averaging catch weight at age for static biological reference points
//DATA_IVECTOR(static_brp_years_M); // model year indices (TMB, starts @ 0) to use for averaging nat mort. for static biological reference points
//DATA_IVECTOR(static_brp_years_R); // ******just use XSPR_R_avg_yrs********model year indices (TMB, starts @ 0) to use for averaging recruitment for static biological reference points
// parameters - general
PARAMETER_VECTOR(mean_rec_pars);
PARAMETER_VECTOR(logit_q); //n_indices (mean/constant q pars)
PARAMETER_VECTOR(q_prior_re); //n_indices (if a prior is used for q, this is used instead of logit_q)
PARAMETER_ARRAY(q_re); //n_years x n_indices (time series of)
PARAMETER_MATRIX(q_repars) //n_indices x 2 (sigma, rho)
PARAMETER_VECTOR(log_F1);
PARAMETER_MATRIX(F_devs);
PARAMETER_VECTOR(log_N1_pars); //length = n_ages or 2
PARAMETER_VECTOR(log_NAA_sigma);
PARAMETER_VECTOR(trans_NAA_rho); // rho_a, rho_y (length = 2)
PARAMETER_ARRAY(log_NAA); // numbers-at-age random effects, (dim = n.yrs-1, n.ages)
PARAMETER_VECTOR(logR_proj); // recruitment (random effects) in proj years, only if SCAA
// PARAMETER_ARRAY(NAA_re); // random effects / deviations from pred NAA for year y, age a (dim = n.yrs-1, n.ages)
PARAMETER_MATRIX(logit_selpars); // mean selectivity, dim = n_selblocks x n_ages + 6 (n_ages for by age, 2 for logistic, 4 for double-logistic)
PARAMETER_VECTOR(selpars_re); // deviations in selectivity parameters (random effects), length = sum(n_selpars)*n_years per block
PARAMETER_MATRIX(sel_repars); // fixed effect parameters controlling selpars_re, dim = n_blocks, 3 (sigma, rho, rho_y)
//PARAMETER_VECTOR(catch_paa_pars);
PARAMETER_MATRIX(catch_paa_pars); //n_fleets x 3
//PARAMETER_VECTOR(index_paa_pars);
PARAMETER_MATRIX(index_paa_pars); //n_indices x 3
PARAMETER_VECTOR(M_a); // mean M-at-age, fixed effects, length = n_ages if M_model = 2 (age-specific), length = 1 if M_model = 1 (constant) or 3 (weight-at-age M)
PARAMETER_ARRAY(M_re); // random effects for year- and age-varying M deviations from mean M_a), dim = n_years x n_M_a
PARAMETER_VECTOR(M_repars); // parameters controlling M_re, length = 3 (sigma_M, rho_M_a, rho_M_y)
PARAMETER(log_b);
PARAMETER_VECTOR(log_catch_sig_scale) //n_fleets
PARAMETER_VECTOR(log_index_sig_scale) //n_indices
// parameters - environmental covariate ("Ecov")
PARAMETER_MATRIX(Ecov_re); // nrows = n_years_Ecov, ncol = N_Ecov
PARAMETER_ARRAY(Ecov_beta); // dim = (2 + n_indices) x n_poly x n_ecov x n_ages, beta_R in eqns 4-5, Miller et al. (2016)
PARAMETER_MATRIX(Ecov_process_pars); // nrows = RW: 2 par (Ecov1, sig), AR1: 3 par (mu, sig, phi); ncol = N_ecov
PARAMETER_MATRIX(Ecov_obs_logsigma); // N_Ecov_years x n_Ecov. options: just given (data), or fixed effect(s)
PARAMETER_MATRIX(Ecov_obs_logsigma_re); // N_Ecov_years x n_Ecov. columns of random effects used if Ecov_obs_sigma_opt = 4
PARAMETER_MATRIX(Ecov_obs_sigma_par); // ncol = N_Ecov, nrows = 2 (mean, sigma of random effects)
Type nll= 0.0; //negative log-likelihood
vector<int> any_index_age_comp(n_indices);
vector<int> any_fleet_age_comp(n_fleets);
vector<Type> SSB(n_years_model + n_years_proj);
matrix<Type> F(n_years_model,n_fleets);
matrix<Type> log_F(n_years_model,n_fleets);
array<Type> pred_CAA(n_years_model+n_years_proj,n_fleets,n_ages);
array<Type> pred_catch_paa(n_years_model+n_years_proj,n_fleets,n_ages);
matrix<Type> pred_catch(n_years_model+n_years_proj,n_fleets);
matrix<Type> pred_log_catch(n_years_model+n_years_proj,n_fleets);
array<Type> pred_IAA(n_years_model+n_years_proj,n_indices,n_ages);
array<Type> pred_index_paa(n_years_model+n_years_proj,n_indices,n_ages);
matrix<Type> pred_indices(n_years_model+n_years_proj,n_indices); // not bias corrected
matrix<Type> pred_log_indices(n_years_model+n_years_proj,n_indices); // bias corrected
array<Type> FAA(n_years_model+n_years_proj,n_fleets,n_ages);
array<Type> log_FAA(n_years_model+n_years_proj,n_fleets,n_ages);
matrix<Type> FAA_tot(n_years_model + n_years_proj,n_ages);
matrix<Type> ZAA(n_years_model + n_years_proj,n_ages);
array<Type> QAA(n_years_model+n_years_proj,n_indices,n_ages);
vector<matrix<Type> > selAA(n_selblocks); // selAA(b)(y,a) gives selectivity by block, year, age; selAA(b) is matrix with dim = n_years x n_ages;
matrix<Type> q(n_years_model+n_years_proj,n_indices);
vector<Type> t_paa(n_ages);
vector<Type> t_pred_paa(n_ages);
int n_toavg = avg_years_ind.size();
//Type SR_a, SR_b, SR_R0, SR_h;
for(int i = 0; i < n_indices; i++)
{
any_index_age_comp(i) = 0;
for(int y = 0; y < n_years_indices; y++) if(use_index_paa(y,i) == 1) any_index_age_comp(i) = 1;
}
for(int i = 0; i < n_fleets; i++)
{
any_fleet_age_comp(i) = 0;
for(int y = 0; y < n_years_catch; y++) if(use_catch_paa(y,i) == 1) any_fleet_age_comp(i) = 1;
}
// Selectivity --------------------------------------------------------------
vector<array<Type> > selpars_re_mats(n_selblocks); // gets selectivity deviations (RE vector, selpars_re) as vector of matrices (nyears x npars), one for each block
vector<matrix<Type> > selpars(n_selblocks); // selectivity parameter matrices for each block, nyears x npars
Type nll_sel = Type(0);
int istart = 0;
int ct = 0;
for(int b = 0; b < n_selblocks; b++){
array<Type> tmp2(n_years_model,n_selpars(b));
tmp2.setZero();
selpars_re_mats(b) = tmp2;
int jstart = 0; // offset for indexing selectivity pars, depends on selectivity model for block b: n_ages (age-specific) + 2 (logistic) + 4 (double-logistic)
if(selblock_models(b) == 2) jstart = n_ages;
if(selblock_models(b) == 3) jstart = n_ages + 2;
if(selblock_models_re(b) > 1){
// fill in sel devs from RE vector, selpars_re (fixed at 0 if RE off)
array<Type> tmp(n_years_selblocks(b), n_selpars_est(b));
for(int j=0; j<n_selpars_est(b); j++){
tmp.col(j) = selpars_re.segment(istart,n_years_selblocks(b));
istart += n_years_selblocks(b);
}
//question: is it faster here to just work on the selectivity parameters as re rather than the deviations?
// likelihood of RE sel devs (if turned on)
Type sigma; // sd selectivity deviations (fixed effect)
Type rho; // among-par correlation selectivity deviations (fixed effect)
Type rho_y; // among-year correlation selectivity deviations (fixed effect)
Type Sigma_sig_sel;
sigma = exp(sel_repars(b,0));
rho = rho_trans(sel_repars(b,1)); // rho_trans ensures correlation parameter is between -1 and 1, see helper_functions.hpp
rho_y = rho_trans(sel_repars(b,2)); // rho_trans ensures correlation parameter is between -1 and 1, see helper_functions.hpp
if((selblock_models_re(b) == 2) | (selblock_models_re(b) == 5)){
// 2D AR1 process on selectivity parameter deviations
Sigma_sig_sel = pow(pow(sigma,2) / ((1-pow(rho_y,2))*(1-pow(rho,2))),0.5);
nll_sel += SCALE(SEPARABLE(AR1(rho),AR1(rho_y)), Sigma_sig_sel)(tmp);
SIMULATE if(simulate_state(2) == 1) if(sum(simulate_period) > 0) SEPARABLE(AR1(rho),AR1(rho_y)).simulate(tmp);
} else {
// 1D AR1 process on selectivity parameter deviations
if(selblock_models_re(b) == 3){ // ar1 across parameters in selblock, useful for age-specific pars.
vector<Type> tmp0 = tmp.matrix().row(0); //random effects are constant across years
Sigma_sig_sel = pow(pow(sigma,2) / (1-pow(rho,2)),0.5);
nll_sel += SCALE(AR1(rho), Sigma_sig_sel)(tmp0);
SIMULATE if(simulate_state(2) == 1) if(sum(simulate_period) > 0)
{
AR1(rho).simulate(tmp0);
for(int y = 0; y < tmp.rows(); y++) for(int i = 0; i < tmp0.size(); i++) tmp(y,i) = tmp0(i);
}
} else { // selblock_models_re(b) = 4, ar1_y, not sure if this one really makes sense.
vector<Type> tmp0 = tmp.matrix().col(0); //random effects are constant within years
Sigma_sig_sel = pow(pow(sigma,2) / (1-pow(rho_y,2)),0.5);
//Sigma_sig_sel = sigma;
nll_sel += SCALE(AR1(rho_y), Sigma_sig_sel)(tmp0);
SIMULATE if(simulate_state(2) == 1) if(sum(simulate_period) > 0)
{
AR1(rho_y).simulate(tmp0);
for(int a = 0; a < tmp.cols(); a++) tmp.col(a) = tmp0;
}
}
}
SIMULATE if(simulate_state(2) == 1) if(sum(simulate_period) > 0) {
tmp = tmp * Sigma_sig_sel;
istart -= n_selpars_est(b) * n_years_selblocks(b); //bring it back to the beginning for this selblock
for(int j=0; j<n_selpars_est(b); j++){
for(int y = 0; y < n_years_selblocks(b); y++){
selpars_re(istart) = tmp(y,j);
istart++;
}
}
}
// construct deviations array with full dimensions (n_years_model instead of n_years_selblocks, n_selpars instead of n_selpars_est)
for(int j=0; j<n_selpars(b); j++){
for(int y=0; y<n_years_model; y++){
if((selblock_years(y,b) == 1) & (selpars_est(b,j+jstart) > 0)){
selpars_re_mats(b)(y,j) = selpars_re(ct);
ct++;
}
}
}
}
// get selpars = mean + deviations
matrix<Type> tmp1(n_years_model, n_selpars(b));
for(int j=jstart; j<(jstart+n_selpars(b)); j++){ // transform from logit-scale
for(int i=0; i<n_years_model; i++){
tmp1(i,j-jstart) = selpars_lower(b,j) + (selpars_upper(b,j) - selpars_lower(b,j)) / (1.0 + exp(-(logit_selpars(b,j) + selpars_re_mats(b).matrix()(i,j-jstart))));
}
}
selpars(b) = tmp1;
}
REPORT(selpars);
REPORT(sel_repars);
REPORT(selpars_re); //even if not simulated
if(do_post_samp(2) == 1) ADREPORT(selpars_re);
REPORT(logit_selpars);
REPORT(nll_sel);
selAA = get_selectivity(n_years_model, n_ages, n_selblocks, selpars, selblock_models); // Get selectivity by block, age, year
nll += nll_sel;
// Environmental covariate process model --------------------------------------
matrix<Type> Ecov_x(n_years_Ecov + n_years_proj_Ecov, n_Ecov); // 'true' estimated Ecov (x_t in Miller et al. 2016 CJFAS)
matrix<Type> nll_Ecov(n_years_Ecov + n_years_proj_Ecov, n_Ecov); // nll contribution each Ecov_re
nll_Ecov.setZero();
// Ecov_model == 0) no Ecov
for(int i = 0; i < n_Ecov; i++){ // loop over Ecovs
// Ecov model option 1: RW
if(Ecov_model(i) == 1){
Type Ecov_sig; // sd (sig_x in Eq1, pg 1262, Miller et al. 2016)
Ecov_sig = exp(Ecov_process_pars(1,i));
Type Ecov1; // Ecov_x in year 1 (fixed effect)
Ecov1 = Ecov_process_pars(0,i);
Ecov_x(0,i) = Ecov1;
nll_Ecov(1,i) -= dnorm(Ecov_re(1,i), Ecov1, Ecov_sig, 1); // Ecov_re(0,i) set to NA
SIMULATE if((simulate_state(3) == 1) & (Ecov_use_re(1,i) == 1)) {
if(simulate_period(0) == 1) {
Ecov_re(1,i) = rnorm(Ecov1, Ecov_sig);
}
}
Ecov_x(1,i) = Ecov_re(1,i);
for(int y = 2; y < n_years_Ecov + n_years_proj_Ecov; y++){
nll_Ecov(y,i) -= dnorm(Ecov_re(y,i), Ecov_re(y-1,i), Ecov_sig, 1);
SIMULATE if((simulate_state(3) == 1) & (Ecov_use_re(y,i) == 1)) {
if(((simulate_period(0) == 1) & (y < n_years_Ecov)) | ((simulate_period(1) == 1) & (y > n_years_Ecov-1))) {
Ecov_re(y,i) = rnorm(Ecov_re(y-1,i), Ecov_sig);
}
}
Ecov_x(y,i) = Ecov_re(y,i);
}
}
// Ecov model option 2: AR1
if(Ecov_model(i) == 2){
Type Ecov_mu; // mean
Type Ecov_phi; // autocorrelation
Type Ecov_sig; // conditional sd
Ecov_mu = Ecov_process_pars(0,i);
Ecov_phi = -Type(1) + Type(2)/(Type(1) + exp(-Ecov_process_pars(2,i)));
Ecov_sig = exp(Ecov_process_pars(1,i));
nll_Ecov(0,i) -= dnorm(Ecov_re(0,i), Type(0), Ecov_sig*exp(-Type(0.5) * log(Type(1) - pow(Ecov_phi,Type(2)))), 1);
SIMULATE if((simulate_state(3) == 1) & (Ecov_use_re(0,i) == 1)) {
if(simulate_period(0) == 1) {
Ecov_re(0,i) = rnorm(Type(0), Ecov_sig*exp(-Type(0.5) * log(Type(1) - pow(Ecov_phi,Type(2)))));
}
}
for(int y = 1; y < n_years_Ecov + n_years_proj_Ecov; y++)
{
nll_Ecov(y,i) -= dnorm(Ecov_re(y,i), Ecov_phi * Ecov_re(y-1,i), Ecov_sig, 1);
SIMULATE if((simulate_state(3) == 1) & (Ecov_use_re(y,i) == 1)) {
if(((simulate_period(0) == 1) & (y < n_years_Ecov)) | ((simulate_period(1) == 1) & (y > n_years_Ecov-1))) {
Ecov_re(y,i) = rnorm(Ecov_phi * Ecov_re(y-1,i), Ecov_sig);
}
}
}
for(int y = 0; y < n_years_Ecov + n_years_proj_Ecov; y++) Ecov_x(y,i) = Ecov_mu + Ecov_re(y,i);
}
// add to nll if estimated (option in projection years to fix Ecov at last or average value)
for(int y = 0; y < n_years_Ecov + n_years_proj_Ecov; y++){
if(Ecov_use_re(y,i) == 1){
nll += nll_Ecov(y,i);
}
}
} // end loop over Ecovs
SIMULATE if(simulate_state(3) == 1) if(sum(simulate_period) > 0) REPORT(Ecov_re);
if(Ecov_model.sum() > 0) if(do_post_samp(3) == 1) ADREPORT(Ecov_re);
// Environmental covariate observation model -------------------------------------
//TODO: Ecov obs are not yet simulated in projection years!!!!!!!!
Type nll_Ecov_obs = Type(0);
Type nll_Ecov_obs_sig = Type(0); // Ecov obs sigma random effects (opt = 4)
matrix<Type> Ecov_obs_sigma(n_years_Ecov, n_Ecov);
for(int i = 0; i < n_Ecov; i++){
for(int y = 0; y < n_years_Ecov; y++){
if(Ecov_obs_sigma_opt(i) == 4){
Type mu_logsigma = Ecov_obs_sigma_par(0,i);
Type sd_logsigma = exp(Ecov_obs_sigma_par(1,i));
nll_Ecov_obs_sig -= dnorm(Ecov_obs_logsigma_re(y,i), mu_logsigma, sd_logsigma, 1);
SIMULATE if(simulate_data(2) == 1) if(simulate_period(0) == 1) {
Ecov_obs_logsigma_re(y,i) = rnorm(mu_logsigma, sd_logsigma);
}
Ecov_obs_sigma(y,i) = exp(Ecov_obs_logsigma_re(y,i));
} else{
Ecov_obs_sigma(y,i) = exp(Ecov_obs_logsigma(y,i));
}
if(Ecov_use_obs(y,i) == 1){
nll_Ecov_obs -= keep(keep_E(y,i)) * dnorm(obsvec(keep_E(y,i)), Ecov_x(y,i), Ecov_obs_sigma(y,i), 1);
nll_Ecov_obs -= keep.cdf_lower(keep_E(y,i)) * log(squeeze(pnorm(obsvec(keep_E(y,i)), Ecov_x(y,i), Ecov_obs_sigma(y,i))));
nll_Ecov_obs -= keep.cdf_upper(keep_E(y,i)) * log(1.0 - squeeze(pnorm(obsvec(keep_E(y,i)), Ecov_x(y,i), Ecov_obs_sigma(y,i))));
SIMULATE if(simulate_data(2) ==1) if(simulate_period(0) == 1) {
Ecov_obs(y,i) = rnorm(Ecov_x(y,i), Ecov_obs_sigma(y,i));
obsvec(keep_E(y,i)) = Ecov_obs(y,i);
}
}
}
}
nll += nll_Ecov_obs_sig;
nll += nll_Ecov_obs;
REPORT(nll_Ecov_obs_sig);
SIMULATE if(simulate_data(2) ==1) if(simulate_period(0) == 1) {
REPORT(Ecov_obs);
REPORT(Ecov_obs_logsigma);
}
// Lag environmental covariates -------------------------------------
// Then use Ecov_out(t) for processes in year t, instead of Ecov_x
int n_effects = Ecov_beta.dim(0); // 2 + n_indices (recruitment, mortality and any catchabilities)
array<Type> Ecov_out(n_years_model + n_years_proj, n_effects, n_Ecov); // Pop model uses Ecov_out(t) for processes in year t (Ecov_x shifted by lag and padded)
Ecov_out.setZero(); // set Ecov_out = 0
for(int i = 0; i < n_Ecov; i++){
for(int t = 0; t < n_effects; t++){
int ct = 0;
for(int y = ind_Ecov_out_start(i,t); y < ind_Ecov_out_end(i,t) + 1 + n_years_proj; y++){
Ecov_out(ct,t,i) = Ecov_x(y,i);
ct++;
}
}
}
// Calculate ecov link model (b1*ecov + b2*ecov^2 + ...) --------------------
// ecov_beta is now 4D array, dim = (2 + n_indices) x n_poly x n_ecov x n_ages
int n_poly = Ecov_beta.dim(1); // now a 4D array dim: (n_effects,n_poly,n_Ecov,n_ages) is second dimension
//vector<matrix<Type>> Ecov_lm(n_Ecov)(n_effects); // ecov linear model for each Ecov, dim = n_years_model + n_years_proj, n_ages
// Ecov_lm.setZero();
// Ecov_lm stores the linear models for each Ecov and where it is used. dim = n_Ecov, n_effects, n_years_model + n_years_proj, n_ages
// n_effects dimension is: 0: recruitment, 1: M, 2-1+n_indices: which catchability it affects
array<Type> Ecov_lm(n_Ecov, n_effects,n_years_model + n_years_proj, n_ages);
//vector<matrix<Type> > Ecov_lm(n_Ecov); // ecov linear model for each Ecov, dim = n_years_model + n_years_proj, n_ages
for(int i = 0; i < n_Ecov; i++){
for(int t = 0; t < n_effects; t++){
vector<Type> thecol(n_years_model + n_years_proj);
for(int y = 0; y < n_years_model + n_years_proj; y++) thecol(y) = Ecov_out(y,t,i);
matrix<Type> X_poly(n_years_model + n_years_proj, n_poly);
X_poly.setZero();
if(n_poly == 1){ // n_poly = 1 if ecov effect is none or linear
X_poly = thecol.matrix();
} else { // n_poly > 1, get poly transformation for ith ecov
X_poly = poly_trans(thecol, n_poly, n_years_model, n_years_proj);
}
for(int y = 0; y < n_years_model + n_years_proj; y++){
for(int a = 0; a < n_ages; a++){
for(int j = 0; j < n_poly; j++){
Ecov_lm(i,t,y,a) += Ecov_beta(t,j,i,a) * X_poly(y,j); // poly transformation returns design matrix, don't need to take powers
}
}
}
}
}
// --------------------------------------------------------------------------
// Calculate mortality (M, F, then Z)
// Natural mortality process model
Type nll_M = Type(0);
if(M_re_model > 1) // random effects on M, M_re = 2D AR1 deviations on M(year,age), dim = n_years x n_M_a
{
Type sigma_M = exp(M_repars(0));
Type rho_M_a = rho_trans(M_repars(1));
Type rho_M_y = rho_trans(M_repars(2));
Type Sigma_M;
// likelihood of M deviations, M_re
if((M_re_model == 2) | (M_re_model == 5)){ //2D AR1: age, year
Sigma_M = pow(pow(sigma_M,2) / ((1-pow(rho_M_y,2))*(1-pow(rho_M_a,2))),0.5);
nll_M += SCALE(SEPARABLE(AR1(rho_M_a),AR1(rho_M_y)), Sigma_M)(M_re); // must be array, not matrix!
SIMULATE if(simulate_state(1) == 1) if(sum(simulate_period) > 0) {
array<Type> Mre_tmp = M_re;
SEPARABLE(AR1(rho_M_a),AR1(rho_M_y)).simulate(Mre_tmp);
Mre_tmp = Sigma_M * Mre_tmp;
for(int y = 0; y < n_years_model + n_years_proj; y++){
if(((simulate_period(0) == 1) & (y < n_years_model)) | ((simulate_period(1) == 1) & (y > n_years_model-1))){
for(int a = 0; a < n_M_a; a++) M_re(y,a) = Mre_tmp(y,a);
}
}
}
} else {
if(M_re_model == 3){ // 1D ar1_a
vector<Type> Mre0 = M_re.matrix().row(0);
Sigma_M = pow(pow(sigma_M,2) / (1-pow(rho_M_a,2)),0.5);
nll_M += SCALE(AR1(rho_M_a), Sigma_M)(Mre0);
SIMULATE if(simulate_state(1) == 1) if(sum(simulate_period) > 0) {
AR1(rho_M_a).simulate(Mre0);
for(int i = 0; i < Mre0.size(); i++) Mre0(i) = Sigma_M * Mre0(i);
for(int y = 0; y < n_years_model + n_years_proj; y++){
for(int i = 0; i < Mre0.size(); i++){
M_re(y,i) = Mre0(i);
}
}
}
} else { // M_re_model = 4, 1D ar1_y
vector<Type> Mre0 = M_re.matrix().col(0);
Sigma_M = pow(pow(sigma_M,2) / (1-pow(rho_M_y,2)),0.5);
nll_M += SCALE(AR1(rho_M_y), Sigma_M)(Mre0);
SIMULATE if(simulate_state(1) == 1) if(sum(simulate_period) > 0) {
AR1(rho_M_y).simulate(Mre0);
for(int i = 0; i < Mre0.size(); i++) Mre0(i) = Sigma_M * Mre0(i);
for(int y = 0; y < n_years_model + n_years_proj; y++){
if(((simulate_period(0) == 1) & (y < n_years_model)) | ((simulate_period(1) == 1) & (y > n_years_model-1))){
for(int a = 0; a < n_M_a; a++) M_re(y,a) = Mre0(y); //all ages mapped to the same annual RE
}
}
}
}
}
if(do_post_samp.sum()==0){
ADREPORT(sigma_M);
ADREPORT(rho_M_a);
ADREPORT(rho_M_y);
}
}
REPORT(nll_M);
nll += nll_M;
REPORT(M_re); //even if M_re not simulated.
if(do_post_samp(1) == 1) ADREPORT(M_re);
REPORT(M_a);
REPORT(M_repars);
// Construct mortality-at-age (MAA)
matrix<Type> MAA(n_years_model + n_years_proj,n_ages);
if(M_model == 2){ // age-specific M
for(int a = 0; a < n_ages; a++) for(int y = 0; y < n_years_model; y++) MAA(y,a) = exp(M_a(a) + M_re(y,a));
} else {
if(M_model == 1){ // constant M
for(int a = 0; a < n_ages; a++) for(int y = 0; y < n_years_model; y++) MAA(y,a) = exp(M_a(0) + M_re(y,a));
} else { // M_model = 3, M is allometric function of weight
for(int a = 0; a < n_ages; a++) for(int y = 0; y < n_years_model; y++) MAA(y,a) = exp(M_a(0) + M_re(y,a) - exp(log_b) * log(waa(waa_pointer_jan1-1,y,a)));
}
}
// add to MAA in projection years
if(do_proj == 1){
if(proj_M_opt == 2){ // use average MAA over avg.yrs
matrix<Type> MAA_toavg(n_toavg,n_ages);
for(int a = 0; a < n_ages; a++){
for(int i = 0; i < n_toavg; i++){
MAA_toavg(i,a) = MAA(avg_years_ind(i),a);
}
}
vector<Type> MAA_proj = MAA_toavg.colwise().mean();
for(int y = n_years_model; y < n_years_model + n_years_proj; y++){
MAA.row(y) = MAA_proj;
}
} else { // proj_M_opt == 1, use M_re and/or ecov_lm in projection years
if(M_model == 2){ // age-specific M
for(int a = 0; a < n_ages; a++) for(int y = n_years_model; y < n_years_model + n_years_proj; y++) MAA(y,a) = exp(M_a(a) + M_re(y,a));
} else {
if(M_model == 1){ // constant M
for(int a = 0; a < n_ages; a++) for(int y = n_years_model; y < n_years_model + n_years_proj; y++) MAA(y,a) = exp(M_a(0) + M_re(y,a));
} else { // M_model = 3, M is allometric function of weight
for(int a = 0; a < n_ages; a++) for(int y = n_years_model; y < n_years_model + n_years_proj; y++) MAA(y,a) = exp(M_a(0) + M_re(y,a) - exp(log_b) * log(waa(waa_pointer_jan1-1,y,a)));
}
}
}
}
// add ecov effect on M (by year, shared across ages)
for(int i=0; i < n_Ecov; i++){
if(Ecov_where(i,1) == 1) { //not sure why Ecov_how is needed for M. if(Ecov_how(i) == 1){ // if ecov i affects M
for(int a = 0; a < n_ages; a++){
for(int y = 0; y < n_years_model + n_years_proj; y++) MAA(y,a) *= exp(Ecov_lm(i,1,y,a));
}
}
}
// prior on M(WAA) coefficient
if(use_b_prior == 1)
{
Type mu = log(0.305);
if(bias_correct_pe == 1) mu -= 0.5*exp(2*log(0.08));
Type lprior_b = dnorm(log_b, mu, Type(0.08), 1);
SIMULATE
{
if(simulate_state(1) == 1) if(sum(simulate_period) > 0) log_b = rnorm(mu, Type(0.08));
REPORT(log_b);
}
REPORT(lprior_b);
nll -= lprior_b;
}
// Survey catchability models
matrix<Type> nll_q(n_years_model+n_years_proj,n_indices);
nll_q.setZero();
vector<Type> sigma_q(n_indices);
sigma_q.setZero();
vector<Type> rho_q(n_indices);
rho_q.setZero();
vector<Type> nll_q_prior(n_indices);
nll_q_prior.setZero();
matrix<Type> logit_q_mat(n_years_model+n_years_proj, n_indices);
logit_q_mat.setZero();
for(int i = 0; i < n_indices; i++) {
//use prior for q? q_prior_re are random effects with mean logit_q (fixed) and sd = logit_q_prior_sigma.
if(use_q_prior(i) == 1){
nll_q_prior(i) -= dnorm(q_prior_re(i), logit_q(i), logit_q_prior_sigma(i), 1);
SIMULATE if(simulate_state(4) == 1) if(sum(simulate_period) > 0){
q_prior_re(i) = rnorm(logit_q(i), logit_q_prior_sigma(i));
}
for(int y = 0; y < n_years_model + n_years_proj; y++) logit_q_mat(y,i) += q_prior_re(i);
}
else for(int y = 0; y < n_years_model + n_years_proj; y++) logit_q_mat(y,i) += logit_q(i);
if(use_q_re(i) > 0) // random effects on q, q_re = AR1 deviations on (year,age), dim = n_years x n_M_a
{
sigma_q(i) = exp(q_repars(i,0)); // conditional sd
rho_q(i) = rho_trans(q_repars(i,1)); // autocorrelation
nll_q(0,i) -= dnorm(q_re(0,i), Type(0), sigma_q(i)*exp(-0.5 * log(1 - pow(rho_q(i),Type(2)))), 1);
SIMULATE if((simulate_state(4) == 1) & (simulate_period(0) == 1)) {
q_re(0,i) = rnorm(Type(0), sigma_q(i)*exp(-0.5 * log(1 - pow(rho_q(i),Type(2)))));
}
logit_q_mat(0,i) += q_re(0,i); //add in q random effects.
for(int y = 1; y < n_years_model + n_years_proj; y++)
{
nll_q(y,i) -= dnorm(q_re(y,i), rho_q(i) * q_re(y-1,i), sigma_q(i), 1);
SIMULATE if(simulate_state(4) == 1) {
if(((simulate_period(0) == 1) & (y < n_years_model)) | ((simulate_period(1) == 1) & (y > n_years_model-1))) {
q_re(y,i) = rnorm(rho_q(i) * q_re(y-1,i), sigma_q(i));
}
}
logit_q_mat(y,i) += q_re(y,i); //add in q random effects.
}
}
}
nll += nll_q.sum() + nll_q_prior.sum();
int Ecov_effects_on_q = 0;
for(int y = 0; y < n_years_model + n_years_proj; y++) {
for(int ind = 0; ind < n_indices; ind++) {
for(int i=0; i < n_Ecov; i++){
if(Ecov_where(i,2+ind) == 1){ // if ecov i affects q and which index
logit_q_mat(y,ind) += Ecov_lm(i,2+ind,y,0);
Ecov_effects_on_q++;
}
}
q(y,ind) = q_lower(ind) + (q_upper(ind) - q_lower(ind))/(1 + exp(-logit_q_mat(y,ind)));
}
}
// Construct survey catchability-at-age (QAA)
for(int i = 0; i < n_indices; i++)
{
// add ecov effect on M (by year, shared across ages)
for(int y = 0; y < n_years_model; y++)
{
for(int a = 0; a < n_ages; a++) QAA(y,i,a) = q(y,i) * selAA(selblock_pointer_indices(y,i)-1)(y,a);
}
//just use last years selectivity for now
if(do_proj == 1) for(int y = n_years_model; y < n_years_model + n_years_proj; y++) for(int a = 0; a < n_ages; a++) {
QAA(y,i,a) = q(y,i) * selAA(selblock_pointer_indices(n_years_model-1,i)-1)(n_years_model-1,a);
}
}
REPORT(logit_q_mat);
if(use_q_re.sum()>0 || Ecov_effects_on_q>0) if(do_post_samp.sum()< 1) ADREPORT(logit_q_mat);
REPORT(sigma_q);
REPORT(rho_q);
REPORT(nll_q);
REPORT(nll_q_prior);
REPORT(q_prior_re); //even if q_prior_re not simulated
REPORT(q_re);
if(do_post_samp(4)==1) ADREPORT(q_re); //even if q_re not simulated.
REPORT(q);
REPORT(QAA);
// Construct fishing mortality-at-age (FAA)
FAA_tot.setZero();
for(int f = 0; f < n_fleets; f++)
{
log_F(0,f) = log_F1(f);
F(0,f) = exp(log_F(0,f));
for(int a = 0; a < n_ages; a++)
{
FAA(0,f,a) = F(0,f) * selAA(selblock_pointer_fleets(0,f)-1)(0,a);
log_FAA(0,f,a) = log(FAA(0,f,a));
FAA_tot(0,a) = FAA_tot(0,a) + FAA(0,f,a);
}
for(int y = 1; y < n_years_model; y++)
{
log_F(y,f) = log_F(y-1,f) + F_devs(y-1,f);
F(y,f) = exp(log_F(y,f));
for(int a = 0; a < n_ages; a++)
{
FAA(y,f,a) = F(y,f) * selAA(selblock_pointer_fleets(y,f)-1)(y,a);
log_FAA(y,f,a) = log(FAA(y,f,a));
FAA_tot(y,a) = FAA_tot(y,a) + FAA(y,f,a);
}
}
}
// Total mortality, Z = F + M (non-projection years only)
for(int y = 0; y < n_years_model; y++) ZAA.row(y) = FAA_tot.row(y) + MAA.row(y);
// ---------------------------------------------------------------------------------
// Set up population model
// Year 1 initialize population
SSB.setZero();
matrix<Type> NAA(n_years_model + n_years_proj,n_ages);
NAA.setZero();
matrix<Type> pred_NAA(n_years_model + n_years_proj,n_ages);
pred_NAA.setZero();
for(int a = 0; a < n_ages; a++)
{
if(N1_model == 0) NAA(0,a) = exp(log_N1_pars(a));
else
{
if(a==0) NAA(0,0) = exp(log_N1_pars(0));
else
{
if(a == n_ages-1) NAA(0,a) = NAA(0,a-1)/(1.0 + exp(-MAA(0,a) - exp(log_N1_pars(1)) * FAA_tot(0,a)/FAA_tot(0,which_F_age(0)-1)));
else NAA(0,a) = NAA(0,a-1)* exp(-MAA(0,a) - exp(log_N1_pars(1)) * FAA_tot(0,a)/FAA_tot(0,which_F_age(0)-1));
}
}
SSB(0) += NAA(0,a) * waa(waa_pointer_ssb-1,0,a) * mature(0,a) * exp(-ZAA(0,a)*fracyr_SSB(0));
pred_NAA(0,a) = NAA(0,a);
}
// get SPR0
vector<Type> M(n_ages), sel(n_ages), mat(n_ages), waacatch(n_ages), waassb(n_ages), log_SPR0(n_years_model + n_years_proj);
int na = n_years_model + n_years_proj;
vector<Type> log_SR_a(na), log_SR_b(na), SR_h(na), SR_R0(na);
for(int y = 0; y < n_years_model + n_years_proj; y++)
{
for(int a = 0; a < n_ages; a++)
{
M(a) = MAA(y,a);
waassb(a) = waa(waa_pointer_ssb-1,y,a);
mat(a) = mature(y,a);
}
log_SPR0(y) = log(get_SPR_0(M, mat, waassb, fracyr_SSB(y)));
}
REPORT(log_SPR0);
// calculate stock-recruit parameters (steepness, R0, a, b)
if(recruit_model > 2) //BH or Ricker SR
{
vector<Type> SR_h_tf(SR_h.size()); //different transformations for BH and Ricker
if(recruit_model == 3) //BH stock recruit
{
if(use_steepness == 1)
{
SR_h.fill(0.2 + 0.8/(1+exp(-mean_rec_pars(0)))); //SR_a * SPR0/(4.0 + SR_a*SPR0);
SR_R0.fill(exp(mean_rec_pars(1))); //(SR_a - 1/SPR0) / SR_b;
log_SR_a = log(4 * SR_h/(exp(log_SPR0)*(1 - SR_h)));
log_SR_b = log((5*SR_h - 1)/((1-SR_h)*SR_R0*exp(log_SPR0)));
}
else
{
log_SR_a.fill(mean_rec_pars(0));
log_SR_b.fill(mean_rec_pars(1));
}
for(int i=0; i < n_Ecov; i++){
if(Ecov_where(i,0) == 1){ // if ecov i affects recruitment
for(int y = 0; y < n_years_model + n_years_proj; y++)
{
// (1) "controlling" = dens-indep mortality or (4) "masking" = metabolic/growth (decreases dR/dS)
if((Ecov_how(i) == 1) | (Ecov_how(i) == 4))
{
log_SR_a(y) += Ecov_lm(i,0,y,0);
}
// (2) "limiting" = carrying capacity or (4) "masking" = metabolic/growth (decreases dR/dS)
if((Ecov_how(i) == 2) | (Ecov_how(i) == 4))
{
log_SR_b(y) += Ecov_lm(i,0,y,0);
}
}
}
}
if(use_steepness != 1)
{
SR_h = exp(log_SR_a) * exp(log_SPR0)/(4.0 + exp(log_SR_a + log_SPR0));
SR_R0 = (exp(log_SR_a) - 1/exp(log_SPR0)) / exp(log_SR_b);
}
SR_h_tf = log(SR_h - 0.2) - log(1 - SR_h);
}
if(recruit_model>3) //Ricker stock recruit
{
if(use_steepness == 1)
{
SR_h.fill(0.2 + exp(mean_rec_pars(0)));
SR_R0.fill(exp(mean_rec_pars(1)));
log_SR_a = 1.25*log(5*SR_h) - log_SPR0;
log_SR_b = log(1.25*log(5*SR_h)/(SR_R0*exp(log_SPR0)));
}
else
{
log_SR_a.fill(mean_rec_pars(0));
log_SR_b.fill(mean_rec_pars(1));
}
for(int i=0; i < n_Ecov; i++){
if(Ecov_where(i,0) == 1){ // if ecov i affects recruitment
for(int y = 0; y < n_years_model + n_years_proj; y++)
{
if(Ecov_how(i) == 1) // "controlling" = dens-indep mortality
{
log_SR_a(y) += Ecov_lm(i,0,y,0);
}
if(Ecov_how(i) == 4) // "masking" = metabolic/growth (decreases dR/dS)
{ //NB: this is not identical to Iles and Beverton (1998), but their definition can give negative values of "b"
log_SR_b(y) += 1.0 + Ecov_lm(i,0,y,0);
}
}
}
}
if(use_steepness != 1)
{
SR_h = 0.2 * exp(0.8*log(exp(log_SR_a) * exp(log_SPR0)));
SR_R0 = log(exp(log_SR_a + log_SPR0))/(exp(log_SR_b + log_SPR0));
}
SR_h_tf = log(SR_h - 0.2);
}
vector<Type> log_SR_R0 = log(SR_R0);
if(do_post_samp.sum()==0){
ADREPORT(log_SR_a);
ADREPORT(log_SR_b);
ADREPORT(SR_h_tf);
ADREPORT(log_SR_R0);
}
REPORT(log_SR_a);
REPORT(log_SR_b);
REPORT(SR_h_tf);
REPORT(log_SR_R0);
}
// ---------------------------------------------------------------------------------
// Population model (get NAA, numbers-at-age, for all years)
array<Type> NAA_devs(n_years_model+n_years_proj-1, n_ages);
NAA_devs.setZero();
for(int y = 1; y < n_years_model + n_years_proj; y++)
{
pred_NAA.row(y) = get_pred_NAA_y(y, recruit_model, mean_rec_pars, SSB, NAA, log_SR_a,
log_SR_b, Ecov_where, Ecov_how, Ecov_lm, ZAA);
// calculate NAA
if(n_NAA_sigma > 1){
// all NAA are estimated (random effects)
for(int a = 0; a < n_ages; a++) NAA(y,a) = exp(log_NAA(y-1,a));
// calculate mean-0 deviations of log NAA (possibly bias-corrected)
for(int a = 0; a < n_ages; a++) NAA_devs(y-1,a) = log_NAA(y-1,a) - log(pred_NAA(y,a));
} else { // only recruitment estimated (either fixed or random effects)
for(int a = 1; a < n_ages; a++) NAA(y,a) = pred_NAA(y,a); // for ages > 1 survival is deterministic
if((n_NAA_sigma == 0) && (y > n_years_model-1)){ //recruit FE, but recruit RE in projection years
NAA(y,0) = exp(logR_proj(y-n_years_model)); // SCAA recruit in projections use diff object (random effect)
for(int a = 1; a < n_ages; a++) NAA_devs(y-1,a) = log(NAA(y,a)) - log(pred_NAA(y,a));
NAA_devs(y-1,0) = logR_proj(y-n_years_model) - log(pred_NAA(y,0));
} else { //recruit RE, or (recruit FE and not projection year)
NAA(y,0) = exp(log_NAA(y-1,0));
for(int a = 0; a < n_ages; a++) NAA_devs(y-1,a) = log(NAA(y,a)) - log(pred_NAA(y,a));
}
}
// calculate F and Z in projection years, here bc need NAA(y) if using F from catch
if(do_proj == 1){ // now need FAA by fleet for projections, use total of average FAA by fleet over avg.yrs
// get selectivity using average over avg.yrs
if(y > n_years_model-1){
waacatch = get_waa_y(waa, y, n_ages, waa_pointer_totcatch);
waassb = get_waa_y(waa, y, n_ages, waa_pointer_ssb);
//n_fleets x n_ages: projected full F is sum of (means across years at age) across fleets
matrix<Type> FAA_proj = get_F_proj(y, n_fleets, proj_F_opt, FAA, NAA, MAA, mature, waacatch, waassb, fracyr_SSB,
log_SPR0, avg_years_ind, n_years_model, which_F_age, percentSPR, proj_Fcatch, percentFXSPR, F_proj_init(y-n_years_model),
log_SR_a, log_SR_b, recruit_model, percentFMSY);
FAA_tot.row(y) = FAA_proj.colwise().sum();
for(int f = 0; f < n_fleets; f++) for(int a = 0; a < n_ages; a++) FAA(y,f,a) = FAA_proj(f,a);
//FAA_tot.row(y) = get_F_proj(y, proj_F_opt, FAA_tot, NAA, MAA, mature, waacatch, waassb, fracyr_SSB, log_SPR0, avg_years_ind, n_years_model,
// which_F_age, percentSPR, proj_Fcatch);
ZAA.row(y) = FAA_tot.row(y) + MAA.row(y);
}
} // end proj F
SSB(y) = get_SSB(NAA,ZAA,waa, mature,y, waa_pointer_ssb, fracyr_SSB);
} // end pop model loop
// --------------------------------------------------------------------------
// NAA random effects
Type nll_NAA = Type(0);
Type NAA_rho_a = rho_trans(trans_NAA_rho(0));
Type NAA_rho_y = rho_trans(trans_NAA_rho(1));
vector<Type> NAA_sigma = exp(log_NAA_sigma);
vector<Type> sigma_a_sig(n_ages);
sigma_a_sig.setZero();
if(n_NAA_sigma == 1){
sigma_a_sig(0) = NAA_sigma(0) / pow((1-pow(NAA_rho_y,2)),0.5);
}
if(n_NAA_sigma > 1){
for(int a=0; a<n_ages; a++) sigma_a_sig(a) = NAA_sigma(NAA_sigma_pointers(a)-1) / pow((1-pow(NAA_rho_y,2))*(1-pow(NAA_rho_a,2)),0.5);
}
if(n_NAA_sigma > 0){
if(do_post_samp.sum()==0){
ADREPORT(NAA_sigma);
ADREPORT(NAA_rho_a);
ADREPORT(NAA_rho_y);
}
if(do_post_samp(0) == 1){
ADREPORT(log_NAA);
}
}
// likelihood of NAA deviations
if((n_NAA_sigma == 0) && (do_proj == 1)){ // SCAA treats recruitment in proj years as random effects with fixed mean, SD
for(int y = 0; y < n_years_proj; y++){
nll_NAA -= dnorm(logR_proj(y), logR_mean, logR_sd, 1);
}
SIMULATE if((simulate_state(0) == 1) & (simulate_period(1) == 1)){ // proj years only
for(int y = 0; y < n_years_proj; y++){
logR_proj(y) = rnorm(logR_mean, logR_sd);
NAA_devs(y+n_years_model-1,0) = logR_proj(y) - log(pred_NAA(y+n_years_model,0));
}
}
REPORT(logR_proj);
}
if(n_NAA_sigma == 1){
if(bias_correct_pe == 1) NAA_devs.col(0) += 0.5*pow(sigma_a_sig(0),2); //make sure this is ok when just recruitment is random.
nll_NAA += SCALE(AR1(NAA_rho_y),sigma_a_sig(0))(NAA_devs.col(0));
SIMULATE if(simulate_state(0) == 1) {
vector<Type> NAAdevs0 = NAA_devs.col(0);
AR1(NAA_rho_y).simulate(NAAdevs0); // sigma = 1, scale below
NAAdevs0 = sigma_a_sig(0) * NAAdevs0;
if(bias_correct_pe == 1) NAAdevs0 -= 0.5*pow(sigma_a_sig(0),2);
for(int y = 0; y < n_years_model + n_years_proj - 1; y++){
if(((simulate_period(0) == 1) & (y < n_years_model - 1)) | ((simulate_period(1) == 1) & (y > n_years_model - 2))){
NAA_devs(y,0) = NAAdevs0(y);
}
}
}
}
if(n_NAA_sigma > 1){
if(bias_correct_pe == 1) for(int a = 0; a < n_ages; a++) NAA_devs.col(a) += 0.5*pow(sigma_a_sig(a),2);
nll_NAA += SEPARABLE(VECSCALE(AR1(NAA_rho_a), sigma_a_sig),AR1(NAA_rho_y))(NAA_devs);
SIMULATE if(simulate_state(0) == 1) {
array<Type> NAAdevs = NAA_devs;
SEPARABLE(VECSCALE(AR1(NAA_rho_a), sigma_a_sig),AR1(NAA_rho_y)).simulate(NAAdevs); // scaled here
if(bias_correct_pe == 1) for(int a = 0; a < n_ages; a++) NAAdevs.col(a) -= 0.5*pow(sigma_a_sig(a),2);
for(int y = 0; y < n_years_model + n_years_proj - 1; y++){
if(((simulate_period(0) == 1) & (y < n_years_model - 1)) | ((simulate_period(1) == 1) & (y > n_years_model - 2))){
for(int a = 0; a < n_ages; a++) NAA_devs(y,a) = NAAdevs(y,a);
}
}
}
}
if((n_NAA_sigma > 0) | (do_proj == 1)) SIMULATE if(simulate_state(0) == 1){ // if n_NAA_sigma = 0 (SCAA), recruitment now random effects in projections
matrix<Type> sims = sim_pop(NAA_devs, recruit_model, mean_rec_pars, SSB, NAA, log_SR_a, log_SR_b, Ecov_where, Ecov_how, Ecov_lm,
n_NAA_sigma, do_proj, proj_F_opt, FAA, FAA_tot, MAA, mature, waa, waa_pointer_totcatch, waa_pointer_ssb, fracyr_SSB, log_SPR0,
avg_years_ind, n_years_model, n_fleets, which_F_age, percentSPR, proj_Fcatch, percentFXSPR, F_proj_init, percentFMSY);
SSB = sims.col(sims.cols()-1);
for(int a = 0; a < n_ages; a++)
{
NAA.col(a) = sims.col(a);
pred_NAA.col(a) = sims.col(a+n_ages);
for(int y = 1;y < n_years_model + n_years_proj; y++) {
log_NAA(y-1,a) = log(NAA(y,a));
for(int f = 0; f < n_fleets; f++){
FAA(y,f,a) = sims(y,2*n_ages + f*n_ages + a);
}
}
}
for(int y = 1;y < n_years_model + n_years_proj; y++) for(int a = 0; a < n_ages; a++){
Type tot = 0;
for(int f = 0; f < n_fleets; f++) tot += FAA(y,f,a);
FAA_tot(y,a) = tot;
}
ZAA = MAA + FAA_tot;
REPORT(sims);
REPORT(log_NAA);
REPORT(log_NAA_sigma);
REPORT(trans_NAA_rho);
}
REPORT(NAA_devs);
REPORT(nll_NAA);
nll += nll_NAA;
// Catch data likelihood
matrix<Type> nll_agg_catch(n_years_model,n_fleets), nll_catch_acomp(n_years_model,n_fleets), agg_catch_proj(n_years_proj,n_fleets);
array<Type> catch_paa_proj(n_fleets, n_years_proj, n_ages);
nll_agg_catch.setZero();
nll_catch_acomp.setZero();
for(int y = 0; y < n_years_model+n_years_proj; y++)
{
//for now just use uncertainty from last year of catch
int usey = y;
if(y > n_years_model-1) usey = n_years_model-1;
//int acomp_par_count = 0;
for(int f = 0; f < n_fleets; f++)
{
pred_catch(y,f) = 0.0;
Type tsum = 0.0;
for(int a = 0; a < n_ages; a++)