diff --git a/.gitignore b/.gitignore
index 7036b5eb..a602a79c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,10 +7,11 @@
*.htp
*.obj
*.log
-ss.tpl
-ss3.tpl
-ss_opt.tpl
-ss_trans.tpl
+Compile/ss.tpl
+Compile/ss_opt.tpl
+Compile/ss_trans.tpl
+Compile/ss3.tpl
+Compile/Make_SS_warn.bat
~$*.*
Compile/ss.log
Compile/ss3.log
diff --git a/Compile/Make_SS_warn.bat b/Compile/Make_SS_warn.bat
index f53510ee..af9b1430 100644
--- a/Compile/Make_SS_warn.bat
+++ b/Compile/Make_SS_warn.bat
@@ -28,6 +28,6 @@ tpl2cpp ss3
g++ -c -std=c++17 -O2 -D_FILE_OFFSET_BITS=64 -DUSE_ADMB_CONTRIBS -D_USE_MATH_DEFINES -I. -I"C:\ADMB-13.2\include" -I"C:\ADMB-13.2\include\contrib" -Wall -Wextra -o ss3.obj ss3.cpp
-g++ -static -o ss3.exe ss3.obj "C:\ADMB-13.2\lib\libadmb-contrib-mingw64-g++12.a"
+g++ -static -o ss3.exe ss3.obj "C:\ADMB-13.2\lib\libadmb-contrib-mingw64-g++13.a"
dir *.exe
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index d2e5025e..74e2f7e2 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -4,6 +4,15 @@
// SS_Label_file # * get_forecast() // calculates forecast quantities, includes all popdy characteristics of the time series, writes forecast-report.sso
// SS_Label_file #
+// Terminology
+// SSB refers to spawning stock biomass, calculated from reproductive output at age (fec()) and numbers-at-age at spawn_month in spawn_seas
+// SSBpR refers to SSB per recruit calculated with equilibrium age composition in equil_calc
+// SPR refers to spawner potential ratio which is the ratio of SSBpR at some level of F to SSBpR with F = 0
+
+// SSBpR_virgin calculated and reported, but never used
+// SSBpR_virgin_adj used only in the alpha-beta spawner-recruitment. _adj means could be updated if SR_update_SSBpR0_timeseries == 1. Also used to get alpha in equil_spawn_recr B-H
+// but note that also uses SSB_virgin_use. So need to align _use with _adj
+
FUNCTION void setup_Benchmark() // and forecast
{
// SS_Label_Info_7.5 #Get averages from selected years to use in forecasts
@@ -113,12 +122,7 @@ FUNCTION void setup_Benchmark() // and forecast
}
}
t = styr + (endyr + 1 - styr) * nseas + spawn_seas - 1;
- fec = Wt_Age_t(t, -2);
- // for (g=1;g<=gmorph;g++)
- // if(use_morph(g)>0 && sx(g)==1)
- // {
- // fec(g)=save_sel_num(t,0,g);
- // }
+// fec = Wt_Age_t(t, -2); this will always be overwritten, so deleting
if (Fcast_Loop_Control(3) == 3) // using mean recr_dist from range of years
{
@@ -536,7 +540,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
dvariable last_F1;
dvariable Closer;
dvariable Vbio1_unfished;
- dvariable SPR_unfished;
+ dvariable SSBpR_unf;
dvariable Vbio_MSY;
dvariable Vbio1_MSY;
dvariable junk;
@@ -569,7 +573,6 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
eq_yr = y;
t_base = y + (y - styr) * nseas - 1;
bio_t_base = styr + (bio_yr - styr) * nseas - 1;
-
// set the Hrate for bycatch fleets so not scaled with other fleets
// bycatch_F(f,s) is created here for use in forecast
for (f = 1; f <= Nfleet; f++)
@@ -643,6 +646,8 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Make_AgeLength_Key(s, subseas);
// SPAWN-RECR: call make_fecundity for benchmarks
+ // this means that any calculation of SSB in benchmark will use the updated fec
+ // but SSBpR0 will only use that updated fec if SR_update_SSBpR0_bmark == 1
if (s == spawn_seas)
{
{
@@ -744,21 +749,41 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SR_parm_work(j) = temp / (Bmark_Yr(10) - Bmark_Yr(9) + 1.);
}
}
- Fishon = 0;
- Recr_unf = mfexp(SR_parm_work(1));
- Do_Equil_Calc(Recr_unf);
- SSB_unf = SSB_equil;
- SR_parm_work(N_SRparm2 + 1) = SSB_unf;
- if (show_MSY == 1)
- report5 << "SR_parm for benchmark: " << SR_parm_work << endl
- << "for years: " << Bmark_Yr(9) << " " << Bmark_Yr(10) << " SSB_virgin was: " << SSB_virgin << endl;
+ if (SR_update_SSBpR0_bmark == 0) // use virgin biology for the spawner-recruitment R0,h calculations in bmark
+ {
+ Recr_unf = Recr_virgin;
+ SSB_unf = SSB_virgin;
+ SSBpR_unf = SSB_unf / Recr_unf;
+ }
+ else // use benchmark biology in the spawner-recruitment R0,h calculations
+ {
+ Fishon = 0;
+ Recr_unf = mfexp(SR_parm_work(1));
+ Do_Equil_Calc(Recr_unf);
+ SSB_unf = SSB_equil; // equilibrium unfished SSB using the benchmark averaged Recr_unf and benchmark averaged biology
+ SSBpR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
+ SSBpR_virgin_adj = SSB_unf / Recr_unf; // update
+ }
if (show_MSY == 1)
- report5 << "Repro_output_by_age_for_morph_1: " << fec(1) << endl;
+ {
+ report5 << "SR_parms for benchmark: " << SR_parm_work << endl
+ << "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl <<
+ "input.SR_update_SSBpR0_rd: " << SR_update_SSBpR0_rd << "flag for updating SSBpR0_Bmark: " << SR_update_SSBpR0_bmark << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Virgin SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SSBpR_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
+ if ( SR_update_SSBpR0_bmark == 1)
+ {
+ report5 << "SPR0 for equilibrium spawner-recruit based on benchmark biology, not virgin biology" << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Benchmark SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " " << SSBpR_unf << " equil: " << Equ_SpawnRecr_Result << endl;
+ }
+
+ }
+ SR_parm_work(N_SRparm2 + 1) = SSB_unf;
Mgmt_quant(1) = SSB_unf;
Mgmt_quant(2) = totbio;
Mgmt_quant(3) = smrybio;
Mgmt_quant(4) = Recr_unf;
-
// find Fspr SS_Label_710
{
if (show_MSY == 1)
@@ -787,7 +812,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SPR_target100 = SPR_target * 100.;
Do_Equil_Calc(equ_Recr);
- SPR_unfished = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
+ SSBpR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
Vbio1_unfished = smrybio; // gets value from equil_calc
if (show_MSY == 1)
{
@@ -833,7 +858,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Fishon = 1;
Do_Equil_Calc(equ_Recr);
- yld1(ii) = 100. * SSB_equil / SPR_unfished; // spawning potential ratio
+ yld1(ii) = 100. * SSB_equil / SSBpR_unf; // spawning potential ratio
}
SPR_actual = yld1(1); // spawning potential ratio
@@ -902,7 +927,8 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
// SPAWN-RECR: calc equil spawn-recr in YPR; need to make this area-specific
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SSB_equil); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSBpR_temp = SSB_equil; // based on most recent call to Do_Equil_Calc
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Bspr = Equ_SpawnRecr_Result(1);
Bspr_rec = Equ_SpawnRecr_Result(2);
@@ -985,7 +1011,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
last_F1 = F1(1);
if (show_MSY == 1)
{
- report5 << j << " " << F1(1) << " " << equ_F_std << " " << SSB_equil / SPR_unfished << " " << YPR_opt << " " << F01_actual << " " << F01_second << " last F1 " << last_F1 << " Closer " << Closer << " delta " << (F01_origin * 0.1 - F01_actual) / (F01_second) << endl;
+ report5 << j << " " << F1(1) << " " << equ_F_std << " " << SSB_equil / SSBpR_unf << " " << YPR_opt << " " << F01_actual << " " << F01_second << " last F1 " << last_F1 << " Closer " << Closer << " delta " << (F01_origin * 0.1 - F01_actual) / (F01_second) << endl;
}
F1(1) += (F01_origin * 0.1 - F01_actual) / (F01_second);
F1(1) = (1. - Closer) * F1(1) + Closer * last_F1;
@@ -1016,8 +1042,8 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Btgt_Fmult = F1(1);
if (rundetail > 0 && mceval_counter == 0 && show_MSY == 1)
echoinput << "Calculated F0.1: " << Btgt_Fmult << endl;
- SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Btgt = Equ_SpawnRecr_Result(1);
Btgt_Rec = Equ_SpawnRecr_Result(2);
YPR_Btgt_enc = YPR_enc; // total encountered yield per recruit
@@ -1028,7 +1054,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
YPR_Btgt_revenue = (PricePerF * YPR_val_vec) * Btgt_Rec; // vector*vector*scalar
// YPR_Btgt_revenue = Price*YPR_ret*Btgt_Rec;
YPR_Btgt_profit = YPR_Btgt_revenue - Cost;
- SPR_Btgt = SSB_equil / SPR_unfished;
+ SPR_Btgt = SSB_equil / SSBpR_unf;
Vbio_Btgt = totbio;
Vbio1_Btgt = smrybio;
Mgmt_quant(7) = equ_F_std;
@@ -1043,9 +1069,12 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
// ******************************************************
if (show_MSY == 1)
{
- report5 << "#" << endl
- << "Find_target_SSB/Bzero; where Bzero is for Bmark years, not Virgin" << endl
- << "Iter Fmult ann_F SPR Catch SSB Recruits SSB/Bzero Tot_catch";
+ report5 << "#" << endl;
+ if (SR_update_SSBpR0_bmark == 0) // use virgin biology for the spawner-recruitment R0,h calculations
+ {report5 << "Find_target_SSB/Bzero; where Bzero is Virgin SSB:" << SSB_unf << " where SSBpR_unf = " << SSBpR_unf << endl;}
+ else
+ {report5 << "Find_target_SSB/Bzero; where Bzero is for Bmark biology and updated SPR0: " << SSB_unf << " where SSBpR_unf = " << SSBpR_unf << endl;}
+ report5 << "Iter Fmult ann_F SPR Catch SSB Recruits SSB/Bzero Tot_catch";
for (p = 1; p <= pop; p++)
for (gp = 1; gp <= N_GP; gp++)
{
@@ -1071,8 +1100,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Nloops = 28;
}
- // Btgttgt=BTGT_target*SSB_virgin; // this is relative to virgin, not to the average biology from benchmark years
- Btgttgt = BTGT_target * SSB_unf; // now relative to Bmark
+ Btgttgt = BTGT_target * SSB_unf;
for (j = 0; j <= Nloops; j++) // loop find Btarget
{
@@ -1107,11 +1135,11 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
// else Hrate for bycatch fleets already set
}
- Do_Equil_Calc(equ_Recr); // where equ_Recr=1.0, so returned SSB_equil is a SSB/R,
- SPR_Btgt = SSB_equil / SPR_unfished;
+ Do_Equil_Calc(equ_Recr); // where equ_Recr=1.0, so returned SSB_equil is in units of SSB/R,
+ SSBpR_temp = SSB_equil;
+ SPR_Btgt = SSBpR_temp / SSBpR_unf; // where SSBpR_unf = SSB_unf / Recr_unf so units of SSB/R; so result is SPR_Btgt = (fished SSB/R) / (unfished SSB/R)
// SPAWN-RECR: calc equil spawn-recr for Btarget calcs; need to make area-specific
- SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
yld1(ii) = Equ_SpawnRecr_Result(1);
}
@@ -1151,6 +1179,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
report5 << " " << SSB_equil_pop_gp(p, gp) * Equ_SpawnRecr_Result(2);
}
report5 << endl;
+// " SSB_unf " << SSB_unf << " recr " < 0)
@@ -2197,6 +2231,8 @@ FUNCTION void Get_Forecast()
for (int Fcast_Loop1 = 1; Fcast_Loop1 <= jloop; Fcast_Loop1++) // for different forecast conditions
{
+// report5 << Fcast_Loop1 << " y: " << 0 << " Repro_output_by_age_for_morph_1 top_forecast: " << fec(1) << endl;
+
switch (Fcast_Loop1) // select which ABC_loops to use
{
case 1: // do OFL only
@@ -2229,7 +2265,7 @@ FUNCTION void Get_Forecast()
if (HarvestPolicy == 0)
report5 << "pop year ABC_Loop season No_buffer bio-all bio-Smry SpawnBio Depletion recruit-0 ";
if (HarvestPolicy <= 2)
- report5 << "pop year ABC_Loop season Ramp&Buffer bio-all bio-Smry SpawnBio Depletion recruit-0 ";
+ report5 << "pop year ABC_Loop season Ramp&Buffer Buffer2 bio-all bio-Smry SpawnBio Depletion recruit-0 ";
if (HarvestPolicy >= 3)
report5 << "pop year ABC_Loop season Ramp bio-all bio-Smry SpawnBio Depletion recruit-0 ";
for (int ff = 1; ff <= N_catchfleets(0); ff++)
@@ -2267,7 +2303,7 @@ FUNCTION void Get_Forecast()
get_growth3(y, t, s, subseas); // in case needed for Lorenzen M
Make_AgeLength_Key(s, subseas); // which also updates Wt_Age_beg, etc.
}
- if (s == spawn_seas)
+// if (s == spawn_seas) //
{
if (WTage_rd == 1)
{
@@ -2278,12 +2314,12 @@ FUNCTION void Get_Forecast()
}
else
{
- get_mat_fec();
+ get_mat_fec(); // does spawnseas and stores in wt_Age_t(t, -2)
}
}
}
}
-
+// report5 << Fcast_Loop1 << " y: " << y << " updated_Repro_output_by_age_for_morph_1 endyr: " << fec(1) << endl;
for (y = endyr + 1; y <= YrMax; y++)
{
t_base = styr + (y - styr) * nseas - 1;
@@ -2473,6 +2509,7 @@ FUNCTION void Get_Forecast()
Wt_Age_mid(s, g) = ALK(ALK_idx, g) * wt_len(s, GP(g)); // use for fisheries with no size selectivity
}
}
+// report5 << Fcast_Loop1 << " y: " << y << " updated_Repro_output_by_age_for_morph_1 annual: " << fec(1) << endl;
Wt_Age_t(t, 0) = Wt_Age_beg(s);
for (g = 1; g <= gmorph; g++)
if (use_morph(g) > 0)
@@ -2554,23 +2591,24 @@ FUNCTION void Get_Forecast()
if (Hermaphro_Option != 0) // get male biomass
{
- MaleSPB(y).initialize();
+ MaleSSB(y).initialize();
for (p = 1; p <= pop; p++)
{
for (g = 1; g <= gmorph; g++)
if (sx(g) == 2 && use_morph(g) > 0) // male; all assumed to be mature
{
- MaleSPB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * natage(t, p, g); // accumulates SSB by area and by growthpattern
+ MaleSSB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * natage(t, p, g); // accumulates SSB by area and by growthpattern
}
}
- if (Hermaphro_maleSPB > 0.0) // add MaleSPB to female SSB
+ if (Hermaphro_maleSSB > 0.0) // add MaleSSB to female SSB
{
- SSB_current += Hermaphro_maleSPB * sum(MaleSPB(y));
+ SSB_current += Hermaphro_maleSSB * sum(MaleSSB(y));
SSB_yr(y) = SSB_current;
}
}
// SPAWN-RECR: get recruitment in forecast; needs to be area-specific
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ // SR_fxn
+ if (SR_update_SSBpR0_timeseries == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
@@ -2593,7 +2631,6 @@ FUNCTION void Get_Forecast()
}
SSB_use = SSB_equil;
}
-
Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr
if (SR_fxn != 7) apply_recdev(Recruits, R0_use); // apply recruitment deviation
if (Fcast_Loop1 < Fcast_Loop_Control(2)) // use expected recruitment this should include environ effect - CHECK THIS
@@ -2628,9 +2665,9 @@ FUNCTION void Get_Forecast()
}
else if (ABC_Loop == 2 && s == 1) // Calc the buffer in season 1, will use last year's spawnbio if multiseas and spawnseas !=1
{
- temp = SSB_unf;
- join1 = 1. / (1. + mfexp(10. * (SSB_current - H4010_bot * temp)));
- join2 = 1. / (1. + mfexp(10. * (SSB_current - H4010_top * temp)));
+
+ join1 = 1. / (1. + mfexp(10. * (SSB_current - H4010_bot * HCR_anchor)));
+ join2 = 1. / (1. + mfexp(10. * (SSB_current - H4010_top * HCR_anchor)));
switch (HarvestPolicy)
{
@@ -2643,8 +2680,8 @@ FUNCTION void Get_Forecast()
// ramp scales catch as f(B) and buffer (H4010_scale) applied to F
{
ABC_buffer(y) = H4010_scale_vec(y) *
- ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
- + (0.0001 + (1.0 - 0.0001) * (H4010_top * temp / SSB_current) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+ ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+ + (0.0001 + (1.0 - 0.0001) * (H4010_top * HCR_anchor / SSB_current) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
) *
(join2) // scale combo
+
@@ -2655,8 +2692,8 @@ FUNCTION void Get_Forecast()
// ramp scales F as f(B) and buffer (H4010_scale) applied to F
{
ABC_buffer(y) = H4010_scale_vec(y) *
- ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
- + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+ ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+ + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
) *
(join2) // scale combo
+
@@ -2667,8 +2704,8 @@ FUNCTION void Get_Forecast()
// ramp scales catch as f(B) and buffer (H4010_scale) applied to catch
{
ABC_buffer(y) = 1.0 *
- ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
- + (0.0001 + (1.0 - 0.0001) * (H4010_top * temp / SSB_current) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+ ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+ + (0.0001 + (1.0 - 0.0001) * (H4010_top * HCR_anchor / SSB_current) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
) *
(join2) // scale combo
+
@@ -2679,8 +2716,8 @@ FUNCTION void Get_Forecast()
// ramp scales F as f(B) and buffer (H4010_scale) applied to catch
{
ABC_buffer(y) = 1.0 *
- ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
- + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+ ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+ + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
) *
(join2) // scale combo
+
@@ -3100,7 +3137,7 @@ FUNCTION void Get_Forecast()
}
if (show_MSY == 1)
{
- report5 << p << " " << y << " " << ABC_Loop << " " << s << " " << ABC_buffer(y) << " " << totbio << " " << smrybio << " ";
+ report5 << p << " " << y << " " << ABC_Loop << " " << s << " " << ABC_buffer(y) << " " << H4010_scale_vec(y) << " " << totbio << " " << smrybio << " ";
if (s == spawn_seas)
{
report5 << SSB_current << " ";
@@ -3192,24 +3229,25 @@ FUNCTION void Get_Forecast()
if (Hermaphro_Option != 0) // get male biomass
{
- MaleSPB(y).initialize();
+ MaleSSB(y).initialize();
for (p = 1; p <= pop; p++)
{
for (g = 1; g <= gmorph; g++)
if (sx(g) == 2 && use_morph(g) > 0) // male; all assumed to be mature
{
- MaleSPB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * elem_prod(natage(t, p, g), mfexp(-Z_rate(t, p, g) * spawn_time_seas)); // accumulates SSB by area and by growthpattern
+ MaleSSB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * elem_prod(natage(t, p, g), mfexp(-Z_rate(t, p, g) * spawn_time_seas)); // accumulates SSB by area and by growthpattern
}
}
- if (Hermaphro_maleSPB > 0.0) // add MaleSPB to female SSB
+ if (Hermaphro_maleSSB > 0.0) // add MaleSSB to female SSB
{
- SSB_current += Hermaphro_maleSPB * sum(MaleSPB(y));
+ SSB_current += Hermaphro_maleSSB * sum(MaleSSB(y));
SSB_yr(y) = SSB_current;
}
}
// SS_Label_Info_24.3.4.1 #Get recruitment from this spawning biomass
- // SPAWN-RECR: calc recruitment in time series; need to make this area-specififc
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ // SPAWN-RECR: calc recruitment in time series; need to make this area-specific
+ // SR_fxn
+ if (SR_update_SSBpR0_timeseries == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
@@ -3484,7 +3522,7 @@ FUNCTION void Get_Forecast()
f = fish_fleet_area(0, ff);
if (fleet_type(f) == 1)
{
- if (ABC_Loop == 2 && HarvestPolicy >= 3)
+ if (ABC_Loop == 2 && HarvestPolicy >= 3) // alternative ABC_buffer approach
{
catch_fleet(t, f) *= H4010_scale_vec(y);
}
@@ -3641,14 +3679,23 @@ FUNCTION void Get_Forecast()
Smry_Table(y, 4) = Mgmt_quant(Fcast_catch_start + y - endyr);
eq_yr = y;
equ_Recr = Recr_unf;
- bio_yr = endyr;
+ bio_yr = y;
Fishon = 0;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
Smry_Table(y, 11) = SSB_equil;
Smry_Table(y, 13) = GenTime;
+ if( SR_fxn == 10 )
+ {
+ temp = SSB_equil / equ_Recr; // current year's SSB/R with current biology at age
+ alpha = mfexp(SR_parm_work(3));
+ beta = mfexp(SR_parm_work(4));
+ SR_parm_byyr(y, 2) = alpha * temp / (4. + alpha * temp); // implied steepness
+ SR_parm_byyr(y, 1) = log( 1. / beta * (alpha - (1. / temp))); // implied ln_R0
+ }
Fishon = 1;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
+// warning < 0)
SPR_std(STD_Yr_Reverse_Ofish(y)) = SSB_equil / Smry_Table(y, 11);
Smry_Table(y, 9) = totbio;
diff --git a/SS_global.tpl b/SS_global.tpl
index 4fa5dcfb..c18dc5af 100644
--- a/SS_global.tpl
+++ b/SS_global.tpl
@@ -1270,10 +1270,10 @@ REPORT_SECTION
if (Do_TG > 0)
report << " TG-fleetcomp " << TG_like1 << endl
<< " TG-negbin " << TG_like2 << endl;
- report << " -log(L): " << obj_fun << " Spbio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr)) << endl;
+ report << " -log(L): " << obj_fun << " SSBio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr)) << endl;
report << endl
- << "Year Spbio Recruitment" << endl;
+ << "Year SSBio Recruitment" << endl;
report << "Virg " << SSB_yr(styr - 2) << " " << exp_rec(styr - 2, 4) << endl;
report << "Init " << SSB_yr(styr - 1) << " " << exp_rec(styr - 1, 4) << endl;
for (y = styr; y <= endyr; y++)
diff --git a/SS_param.tpl b/SS_param.tpl
index 9ddf9a86..b6f9341b 100644
--- a/SS_param.tpl
+++ b/SS_param.tpl
@@ -159,13 +159,14 @@ PARAMETER_SECTION
3darray recr_dist_endyr(1,N_GP*gender,1,N_settle_timings,1,pop);
!!// SS_Label_Info_5.1.2 #Create SR_parm vector, recruitment vectors
init_bounded_number_vector SR_parm(1,N_SRparm3,SR_parm_LO,SR_parm_HI,SR_parm_PH)
- matrix SR_parm_byyr(styr-3,YrMax,1,N_SRparm2+1) // R0, steepness, parm3, sigmar, rec_dev_offset, R1, rho, SPB Time_vary implementation of spawner-recruitment
+ matrix SR_parm_byyr(styr-3,YrMax,1,N_SRparm2+1) // R0, steepness, parm3, sigmar, rec_dev_offset, R1, rho, SSB Time_vary implementation of spawner-recruitment
vector SR_parm_virg(1,N_SRparm2+1)
vector SR_parm_work(1,N_SRparm2+1)
number two_sigmaRsq;
number half_sigmaRsq;
number sigmaR;
- number SPR_virgin;
+ number SSBpR_virgin;
+ number SSBpR_virgin_adj;
number regime_change;
number rho;
number dirichlet_Parm;
@@ -234,21 +235,25 @@ PARAMETER_SECTION
init_bounded_vector Fcast_recruitments(recdev_end+1,s,recdev_LO,recdev_HI,Fcast_recr_PH2)
init_bounded_vector Fcast_impl_error(endyr+1,j,-1,1,k)
vector ABC_buffer(endyr+1,YrMax);
+ number HCR_anchor // basis (denominator) for inflection in control rule. Select virgin SSB or benchmark SSB
// SPAWN-RECR: define some spawning biomass and recruitment entities
number SSB_virgin
number Recr_virgin
number SSB_vir_LH
- number SSB_unf
+ number SSB_unf // SSB unfished, based on benchmark biology
number Recr_unf
+ number SSB_use
+ number R0_use; // annually updated value if SR_update_SSBpR0_timeseries == 1
+ number SSB_deplete // SSB that will be used as denominator for depletion calculations and as basis for control rule inflection
number SSB_current; // Spawning biomass
number SSB_equil;
number SPR_trial
number SPR_actual;
- number SPR_temp; // used to pass quantity into Equil_SpawnRecr
+ number SSBpR_temp; // SSB per Recruit; used to pass quantity into Equil_SpawnRecr
number Recruits; // Age0 Recruits
number equ_mat_bio
number equ_mat_num
@@ -321,7 +326,7 @@ PARAMETER_SECTION
!!k=0;
!!if(Hermaphro_Option!=0) k=1;
- 3darray MaleSPB(styr-3,YrMax*k,1,pop,1,N_GP) //Male Spawning biomass
+ 3darray MaleSSB(styr-3,YrMax*k,1,pop,1,N_GP) //Male Spawning biomass
matrix SSB_equil_pop_gp(1,pop,1,N_GP);
matrix MaleSSB_equil_pop_gp(1,pop,1,N_GP);
@@ -416,6 +421,7 @@ PARAMETER_SECTION
// note that bycatch_F(1,Nfleet,1,nseas) has similar role
number alpha;
number beta;
+ number steepness;
number GenTime;
number Yield;
number Adj4010;
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index e6c64a53..6b51ce91 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -337,17 +337,42 @@ FUNCTION void get_initial_conditions()
Fishon = 0;
virg_fec = fec;
Recr.initialize(); // will store recruitment by area
+
+ // SPAWN-RECR: get expected recruitment globally or by area
+ if (recr_dist_area == 1 || pop == 1) // do global spawn_recruitment calculations
+ {
+ equ_Recr = 1.0;
+ Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation. Returns SPR because R = 1.0
+ SSBpR_virgin = SSB_equil; // spawners per recruit. Needed for Sr_fxn = 10
+ SSBpR_virgin_adj = SSB_equil; // also needed for Sr_fxn 10. Will get revised in benchmark to use averaged biology if requested.
+ if(SR_fxn == 10) // B-H with a,b
+ {
+ // WHAM based on R = A*S/(1+B*S)
+ // 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)));
+ // h = a * SPR0 / (4. + a * SPR0)
+ // R0 = 1/b * (a-1/SPR0)
+
+ alpha = mfexp(SR_parm(3));
+ beta = mfexp(SR_parm(4));
+ steepness = alpha * SSBpR_virgin_adj / (4. + alpha * SSBpR_virgin_adj);
+ Recr_virgin = 1. / beta * (alpha - (1. / SSBpR_virgin_adj));
+// warning << " before AB_calcs " << "parm " << SR_parm(1) << " calc " << log(Recr_virgin) << endl;
+ SR_parm(1) = log(Recr_virgin);
+ SR_parm(2) = steepness;
+ }
+ else {
+ Recr_virgin = mfexp(SR_parm(1));
+ }
+
for (int i = 1; i <= N_SRparm2; i++)
{
SR_parm_byyr(eq_yr, i) = SR_parm(i);
SR_parm_virg(i) = SR_parm(i);
SR_parm_work(i) = SR_parm(i);
}
-
- // SPAWN-RECR: get expected recruitment globally or by area
- if (recr_dist_area == 1 || pop == 1) // do global spawn_recruitment calculations
- {
- Recr_virgin = mfexp(SR_parm(1));
+// if (SR_fxn == 3) warning << "tester_A: " << SR_parm_work(1) << " base: " << SR_parm(1) << endl;
+// if (SR_fxn == 10) warning << "tester_A: " << SR_parm_work(4) << " base: " << SR_parm(4) << endl;
equ_Recr = Recr_virgin;
exp_rec(eq_yr, 1) = Recr_virgin; // expected Recr from s-r parms
exp_rec(eq_yr, 2) = Recr_virgin;
@@ -355,22 +380,22 @@ FUNCTION void get_initial_conditions()
exp_rec(eq_yr, 4) = Recr_virgin;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
SSB_virgin = SSB_equil;
- SPR_virgin = SSB_equil / Recr_virgin; // spawners per recruit
- if(Do_Benchmark==0)
+// SSBpR_virgin = SSB_equil / Recr_virgin; // spawners per recruit already calculated
+ if(Do_Benchmark==0) // assign values that would be created in benchmark section
{
- Mgmt_quant(1)=SSB_virgin;
- SSB_unf=SSB_virgin;
- Recr_unf=Recr_virgin;
- Mgmt_quant(2)=totbio; // from equil calcs
- Mgmt_quant(3)=smrybio; // from equil calcs
- Mgmt_quant(4)=Recr_virgin;
+ Mgmt_quant(1) = SSB_virgin; // can be overwritten in benchmark by updated SSB_unf
+ SSB_unf = SSB_virgin;
+ Recr_unf = Recr_virgin;
+ Mgmt_quant(2) = totbio; // from equil calcs
+ Mgmt_quant(3) = smrybio; // from equil calcs
+ Mgmt_quant(4) = Recr_virgin;
}
Smry_Table(styr - 2, 1) = totbio; // from equil calcs
Smry_Table(styr - 2, 2) = smrybio; // from equil calcs
Smry_Table(styr - 2, 3) = smrynum; // from equil calcs
SSB_pop_gp(eq_yr) = SSB_equil_pop_gp; // dimensions of pop x N_GP
if (Hermaphro_Option != 0)
- MaleSPB(eq_yr) = MaleSSB_equil_pop_gp;
+ MaleSSB(eq_yr) = MaleSSB_equil_pop_gp;
SSB_yr(eq_yr) = SSB_equil;
SR_parm_byyr(eq_yr, N_SRparm2 + 1) = SSB_equil;
SR_parm_virg(N_SRparm2 + 1) = SSB_equil;
@@ -454,6 +479,9 @@ FUNCTION void get_initial_conditions()
else
{
SR_parm_work(f) = parm_timevary(SR_parm_timevary(f), eq_yr);
+// warning << "tester_B: " << SR_parm_work(f) << " timevary " << " base " << SR_parm(f) < 0)
Morphcomp_exp.initialize();
// SS_Label_Info_24.0 #Retrieve spawning biomass and recruitment from the initial equilibrium
- // SPAWN-RECR: begin of time series, retrieve last spbio and recruitment
+ // SPAWN-RECR: begin of time series, retrieve last SSBio and recruitment
SSB_current = SSB_yr(styr); // need these initial assignments in case recruitment distribution occurs before spawnbio&recruits
if (recdev_doit(styr - 1) > 0)
{
@@ -704,6 +730,7 @@ FUNCTION void get_time_series()
else
{
SR_parm_work(f) = parm_timevary(SR_parm_timevary(f), y);
+// warning << "tester_C: " << SR_parm_work(f) << " timevary_year " << endl;
}
SR_parm_byyr(y, f) = SR_parm_work(f);
}
@@ -954,7 +981,7 @@ FUNCTION void get_time_series()
}
}
// SS_Label_Info_24.2.2 #Compute spawning biomass if this is spawning season so recruits could occur later this season
- // SPAWN-RECR: calc SPB in time series if spawning is at beginning of the season
+ // SPAWN-RECR: calc SSB in time series if spawning is at beginning of the season
if (s == spawn_seas && spawn_time_seas < 0.0001) // compute spawning biomass if spawning at beginning of season so recruits could occur later this season
{
SSB_pop_gp(y).initialize();
@@ -980,25 +1007,26 @@ FUNCTION void get_time_series()
if (Hermaphro_Option != 0) // get male biomass
{
- MaleSPB(y).initialize();
+ MaleSSB(y).initialize();
for (p = 1; p <= pop; p++)
{
for (g = 1; g <= gmorph; g++)
if (sx(g) == 2 && use_morph(g) > 0) // male; all assumed to be mature
{
- MaleSPB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * natage(t, p, g); // accumulates SSB by area and by growthpattern
+ MaleSSB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * natage(t, p, g); // accumulates SSB by area and by growthpattern
}
}
- if (Hermaphro_maleSPB > 0.0) // add MaleSPB to female SSB
+ if (Hermaphro_maleSSB > 0.0) // add MaleSSB to female SSB
{
- SSB_current += Hermaphro_maleSPB * sum(MaleSPB(y));
+ SSB_current += Hermaphro_maleSSB * sum(MaleSSB(y));
SSB_yr(y) = SSB_current;
}
}
- // SS_Label_Info_24.2.3 #Get the total recruitment produced by this spawning biomass
+ // SS_Label_Info_24.2.3 #Get the total recruitment produced by this spawning biomass at the beginning of the season
// SPAWN-RECR: calc recruitment in time series; need to make this area-specific
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ // SR_Fxn relevant keyword
+ if (SR_update_SSBpR0_timeseries == 0) // SRparm are not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
@@ -1010,7 +1038,9 @@ FUNCTION void get_time_series()
Fishon = 0;
eq_yr = y;
bio_yr = y;
- Do_Equil_Calc(R0_use); // call function to do equilibrium calculation
+ Do_Equil_Calc(R0_use); // call function to do equilibrium calculation with current year's biology and adjusted R0
+ SSB_use = SSB_equil;
+ SSBpR_virgin_adj = SSB_use / R0_use; // update
if (fishery_on_off == 1)
{
Fishon = 1;
@@ -1019,9 +1049,8 @@ FUNCTION void get_time_series()
{
Fishon = 0;
}
- SSB_use = SSB_equil;
}
- Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr
+ Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr using either virgin or adjusted R0 and SSB0
if (SR_fxn != 7) apply_recdev(Recruits, R0_use); // apply recruitment deviation
// distribute Recruitment of age 0 fish among the current and future settlements; and among areas and morphs
// use t offset for each birth event: Settlement_offset(settle)
@@ -1440,36 +1469,38 @@ FUNCTION void get_time_series()
SSB_yr(y) = SSB_current;
if (Hermaphro_Option != 0) // get male biomass
{
- MaleSPB(y).initialize();
+ MaleSSB(y).initialize();
for (p = 1; p <= pop; p++)
{
for (g = 1; g <= gmorph; g++)
if (sx(g) == 2 && use_morph(g) > 0) // male; all assumed to be mature
{
- MaleSPB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * elem_prod(natage(t, p, g), mfexp(-Z_rate(t, p, g) * spawn_time_seas)); // accumulates SSB by area and by growthpattern
+ MaleSSB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * elem_prod(natage(t, p, g), mfexp(-Z_rate(t, p, g) * spawn_time_seas)); // accumulates SSB by area and by growthpattern
}
}
- if (Hermaphro_maleSPB > 0.0) // add MaleSPB to female SSB
+ if (Hermaphro_maleSSB > 0.0) // add MaleSSB to female SSB
{
- SSB_current += Hermaphro_maleSPB * sum(MaleSPB(y));
+ SSB_current += Hermaphro_maleSSB * sum(MaleSSB(y));
SSB_yr(y) = SSB_current;
}
}
- // SS_Label_Info_24.3.4.1 #Get recruitment from this spawning biomass
+ // SS_Label_Info_24.3.4.1 #Get recruitment from this spawning biomass at some time during the season
// SPAWN-RECR: calc recruitment in time series; need to make this area-specific
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ // SR_fxn
+ if (SR_update_SSBpR0_timeseries == 0) // SR parms are not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
}
- else
+ else // update SSB_use and R0_use where will update SSBpR0 inside the Spawn_recr fxn
{
- R0_use = mfexp(SR_parm_work(1));
- equ_Recr = R0_use;
+ R0_use = mfexp(SR_parm_work(1)); // check to be sure this works when R0 is derived from B-H with alpha, beta parameters
Fishon = 0;
eq_yr = y;
bio_yr = y;
Do_Equil_Calc(R0_use); // call function to do equilibrium calculation
+ SSB_use = SSB_equil;
+ SSBpR_virgin_adj = SSB_use / R0_use; // update
if (fishery_on_off == 1)
{
Fishon = 1;
@@ -1478,7 +1509,6 @@ FUNCTION void get_time_series()
{
Fishon = 0;
}
- SSB_use = SSB_equil;
}
Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr
@@ -1764,6 +1794,14 @@ FUNCTION void get_time_series()
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation with current year's biology
Smry_Table(y, 11) = SSB_equil;
Smry_Table(y, 13) = GenTime;
+ if( SR_fxn == 10 )
+ {
+ temp = SSB_equil / Recr_virgin; // current year's SSB/R with current biology at age
+ alpha = mfexp(SR_parm_work(3));
+ beta = mfexp(SR_parm_work(4));
+ SR_parm_byyr(y, 2) = alpha * temp / (4. + alpha * temp); // implied steepness
+ SR_parm_byyr(y, 1) = log( 1. / beta * (alpha - (1. / temp))); // implied ln_R0
+ }
Fishon = 1;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation with current year's biology and F
if (STD_Yr_Reverse_Ofish(y) > 0)
@@ -2244,9 +2282,9 @@ FUNCTION void Do_Equil_Calc(const prevariable& equ_Recr)
SSB_equil = sum(SSB_equil_pop_gp);
GenTime /= SSB_equil;
smryage /= smrynum;
- if (Hermaphro_maleSPB > 0.0) // add MaleSPB to female SSB
+ if (Hermaphro_maleSSB > 0.0) // add MaleSSB to female SSB
{
- SSB_equil += Hermaphro_maleSPB * sum(MaleSSB_equil_pop_gp);
+ SSB_equil += Hermaphro_maleSSB * sum(MaleSSB_equil_pop_gp);
}
} // end equil calcs
diff --git a/SS_proced.tpl b/SS_proced.tpl
index 534adf53..3e545ba7 100644
--- a/SS_proced.tpl
+++ b/SS_proced.tpl
@@ -206,15 +206,15 @@ PROCEDURE_SECTION
temp = sqrt((temp + 0.0000001) / (double(recdev_end - recdev_start + 1)));
if (mcmc_counter == 0 && mceval_counter == 0)
{
- cout << current_phase() << " " << niter << " -log(L): " << obj_fun << " Spbio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
+ cout << current_phase() << " " << niter << " -log(L): " << obj_fun << " SSBio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
}
else if (mcmc_counter > 0)
{
- cout << " MCMC: " << mcmc_counter << " -log(L): " << obj_fun << " Spbio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
+ cout << " MCMC: " << mcmc_counter << " -log(L): " << obj_fun << " SSBio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
}
else if (mceval_counter > 0)
{
- cout << " MCeval: " << mceval_counter << " -log(L): " << obj_fun << " Spbio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
+ cout << " MCeval: " << mceval_counter << " -log(L): " << obj_fun << " SSBio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
}
if (F_Method > 1 && sum(catch_like) > 0.01)
{
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index e33d5173..054f7c74 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -1148,11 +1148,11 @@
int Hermaphro_seas;
int Hermaphro_firstage;
number Hermaphro_seas_rd;
- number Hermaphro_maleSPB;
+ number Hermaphro_maleSSB;
LOCAL_CALCS
// clang-format on
Hermaphro_seas = 0;
- Hermaphro_maleSPB = 0.0;
+ Hermaphro_maleSSB = 0.0;
Hermaphro_firstage = 0;
MGparm_Hermaphro = 0;
@@ -1175,8 +1175,8 @@
// so 2.3 will do switch in season 2 beginning with age 3.
echoinput << Hermaphro_seas << " Hermaphro_season (-1 means all seasons)" << endl;
echoinput << Hermaphro_firstage << " Hermaphro_firstage (from decimal part of seas input; note that firstage can only be a single digit, so 9 is max" << endl;
- *(ad_comm::global_datafile) >> Hermaphro_maleSPB; // read as a fraction (0.0 to 1.0) of the male SSB added into the total SSB
- echoinput << Hermaphro_maleSPB << " Hermaphro_maleSPB " << endl;
+ *(ad_comm::global_datafile) >> Hermaphro_maleSSB; // read as a fraction (0.0 to 1.0) of the male SSB added into the total SSB
+ echoinput << Hermaphro_maleSSB << " Hermaphro_maleSSB " << endl;
}
// clang-format off
END_CALCS
@@ -1656,6 +1656,7 @@
!!// SS_Label_Info_4.5.4 #Set up time-varying parameters for MG parms
int timevary_used;
+ int timevary_MG_firstyr;
int timevary_parm_cnt_MG;
int timevary_parm_start_MG;
@@ -1704,6 +1705,7 @@
timevary_parm_start_MG = 0;
timevary_parm_cnt_MG = 0;
timevary_used = 0;
+ timevary_MG_firstyr = YrMax;
MGparm_timevary.initialize();
ivector block_design_null(1, 1);
block_design_null.initialize();
@@ -1804,6 +1806,7 @@
{
MG_active(f) = 1;
timevary_MG(y, 0) = 1; // tracks active status for all MG types
+ if(timevary_MG_firstyr == YrMax) timevary_MG_firstyr = y; // save for reporting in MSY and spawn_recruit output
}
}
// timevary growth or maturity and Maunder M refers to that maturity
@@ -1937,13 +1940,26 @@
!!// SS_Label_Info_4.6 #Read setup for Spawner-Recruitment parameters
// SPAWN-RECR: read setup for SR parameters: LO, HI, INIT, PRIOR, PRtype, CV, PHASE
init_int SR_fxn
-!!echoinput< 0 || SR_parm_timevary(4) > 0) ) {SR_update_parm = 1;} // alpha or beta is time-varying
+ else if ((SR_parm_timevary(1) > 0 || SR_parm_timevary(2) > 0) ) {SR_update_parm = 1;} // R0 or steepness is time-varying
}
}
N_SRparm3 = N_SRparm2;
@@ -2122,6 +2155,35 @@
echoinput << " SR timevary_parm_cnt start and end " << timevary_parm_start_SR << " " << timevary_parm_cnt_SR << endl;
echoinput << "link to timevary parms: " << SR_parm_timevary << endl;
}
+
+// SR_update_SSBpR0_rd values
+// 1 best: update SSBpR0 for benchmark and for time series only if SRparm R0 or h (not regime) is set to have time-varying property
+// 2 incorrect, relic: always update SSBpR0 for benchmark's use of spawner-recruitment (old, incorrect SS3 approach), but only for the time series if there is a timevary SR parm
+// 3 option: do not update SSBpR0 (do keep start year SPR0), even if R0 or h is set to have time-varying property
+
+ switch (SR_update_SSBpR0_rd)
+ {
+ case 0:
+ if(timevary_MG_firstyr < YrMax) // timevary biology exists and SR_update not set
+ {
+ warnstream << "user must select 1, 2, or 3 for updating SPR0 flag (formerly labelled future feature in SR input) because there is time-varying biology";
+ write_message (FATAL, 0); // EXIT!
+ }
+ break;
+ case 1:
+ SR_update_SSBpR0_bmark = 1 * SR_update_parm; // but conditional on SRparm_timevary, so update value there
+ SR_update_SSBpR0_timeseries = 1 * SR_update_parm;
+ break;
+ case 2:
+ SR_update_SSBpR0_bmark = 1;
+ SR_update_SSBpR0_timeseries = 0;
+ break;
+ case 3:
+ SR_update_SSBpR0_bmark = 0;
+ SR_update_SSBpR0_timeseries = 0;
+ break;
+ }
+ echoinput << "SRupdate " << SR_update_parm << " use " << SR_update_SSBpR0_bmark << " rd " << SR_update_SSBpR0_rd<< endl;
echoinput << "SR_Npar and N_SRparm2 and N_SRparm3: " << N_SRparm(SR_fxn) << " " << N_SRparm2 << " " << N_SRparm3 << endl;
// clang-format off
END_CALCS
@@ -5754,7 +5816,7 @@
Extra_Std_N += YrMax - (styr - 2) + 1;
if (More_Std_Input(12) == 2) Extra_Std_N += YrMax - (styr - 2) + 1; // for recruitment
}
- // add 3 values for ln(Spbio)
+ // add 3 values for ln(SSBio)
// (years are automatically generated as startyr, mid-point, and endyr)
Do_se_LnSSB = Extra_Std_N + 1;
Extra_Std_N += 3;
@@ -6808,23 +6870,23 @@
}
}
- // output ln(SPB) std for selected years
- echoinput << " do ln(SPB) std labels for 3 years" << endl;
+ // output ln(SSB) std for selected years
+ echoinput << " do ln(SSB) std labels for 3 years" << endl;
CoVar_Count++;
j++;
active_parm(CoVar_Count) = j;
sprintf(onenum, "%d", styr);
- ParmLabel += "ln(SPB)_" + onenum + CRLF(1);
+ ParmLabel += "ln(SSB)_" + onenum + CRLF(1);
CoVar_Count++;
j++;
active_parm(CoVar_Count) = j;
sprintf(onenum, "%d", int((endyr + styr) / 2));
- ParmLabel += "ln(SPB)_" + onenum + CRLF(1);
+ ParmLabel += "ln(SSB)_" + onenum + CRLF(1);
CoVar_Count++;
j++;
active_parm(CoVar_Count) = j;
sprintf(onenum, "%d", endyr);
- ParmLabel += "ln(SPB)_" + onenum + CRLF(1);
+ ParmLabel += "ln(SSB)_" + onenum + CRLF(1);
if (Do_se_smrybio > 0)
{
diff --git a/SS_readdata_330.tpl b/SS_readdata_330.tpl
index d212bcde..96658403 100644
--- a/SS_readdata_330.tpl
+++ b/SS_readdata_330.tpl
@@ -4345,8 +4345,8 @@
"even when the base is set to the mean of earlier recruitments" << endl;
}
- echoinput << Fcast_Loop_Control(5) << " #echo: loop control 5 not used" << endl;
-
+ echoinput << Fcast_Loop_Control(5) << " #control rule anchor: 1=unfished_benchmark_SSB(old_approach), 2=virgin_SSB " << endl;
+
echoinput << "#next enter year in which Fcast loop 3 caps and allocations begin to be applied" << endl;
*(ad_comm::global_datafile) >> Fcast_Cap_FirstYear;
echoinput << Fcast_Cap_FirstYear << " # echoed value" << endl;
@@ -5035,10 +5035,10 @@
if (STD_Yr_Reverse(y) > 0)
{
j++;
- STD_Yr_Reverse(y) = j; // use for SPB and recruitment
+ STD_Yr_Reverse(y) = j; // use for SSB and recruitment
if (y >= styr)
{
- // depletion must start in year AFTER first catch. It could vary earlier if recdevs happened enough earlier to change spbio, but this is not included
+ // depletion must start in year AFTER first catch. It could vary earlier if recdevs happened enough earlier to change SSBio, but this is not included
if ((depletion_basis > 0 && y > first_catch_yr) || y == endyr)
{
N_STD_Yr_Dep++;
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 93d3fe73..645dbcf3 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -6,7 +6,7 @@
//********************************************************************
/* SS_Label_FUNCTION 43 Spawner-recruitment function */
// SPAWN-RECR: function: to calc R from S
-FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariable& Recr_virgin_adj, const prevariable& SSB_current)
+FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_use, const prevariable& Recr_virgin_use, const prevariable& SSB_current)
{
RETURN_ARRAYS_INCREMENT();
dvariable NewRecruits;
@@ -22,6 +22,7 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
dvariable SRZ_0;
dvariable srz_min;
dvariable SRZ_surv;
+// warning << y << " Tester_R0 " << Recr_virgin_use << " SSB0 " << SSB_virgin_use << " SSB_curr: " << SSB_current << endl;
// SS_Label_43.1 add 0.1 to input spawning biomass value to make calculation more rebust
SSB_curr_adj = SSB_current + 0.100; // robust
@@ -29,7 +30,7 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
regime_change = SR_parm_work(N_SRparm2 - 1); // this is a persistent deviation off the S/R curve
// SS_Label_43.3 calculate expected recruitment from the input spawning biomass and the SR curve
- // functions below use Recr_virgin_adj,SSB_virgin_adj which could have been adjusted adjusted above from R0,SSB_virgin
+ // functions below use Recr_virgin_use,SSB_virgin_use which could have been adjusted adjusted above from R0,SSB_virgin
switch (SR_fxn)
{
case 1: // previous placement for B-H constrained
@@ -42,36 +43,30 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 2: // ricker
{
steepness = SR_parm_work(2);
- NewRecruits = Recr_virgin_adj * SSB_curr_adj / SSB_virgin_adj * mfexp(steepness * (1. - SSB_curr_adj / SSB_virgin_adj));
+ NewRecruits = Recr_virgin_use * SSB_curr_adj / SSB_virgin_use * mfexp(steepness * (1. - SSB_curr_adj / SSB_virgin_use));
break;
}
// SS_Label_43.3.3 Beverton-Holt
case 3: // Beverton-Holt
{
steepness = SR_parm_work(2);
- alpha = 4.0 * steepness * Recr_virgin / (5. * steepness - 1.);
- beta = (SSB_virgin_adj * (1. - steepness)) / (5. * steepness - 1.);
- NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_curr_adj) /
- (SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);
+ NewRecruits = (4. * steepness * Recr_virgin_use * SSB_curr_adj) /
+ (SSB_virgin_use * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);
break;
}
- // Beverton-Holt with alpha beta
- /*
- case 3: // Beverton-Holt
+ case 10: // Beverton-Holt with alpha beta per WHAM: R = A*S/(1+B*S)
{
- steepness=SR_parm_work(2);
- alpha = 4.0 * steepness*Recr_virgin / (5.*steepness-1.);
- beta = (SSB_virgin_adj*(1.-steepness)) / (5.*steepness-1.);
- NewRecruits = (alpha*SSB_curr_adj) / (beta+SSB_curr_adj);
+ dvariable alpha = mfexp(SR_parm_work(3));
+ dvariable beta = mfexp(SR_parm_work(4));
+ NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
break;
}
- */
// SS_Label_43.3.4 constant expected recruitment
case 4: // none
{
- NewRecruits = Recr_virgin_adj;
+ NewRecruits = Recr_virgin_use;
break;
}
// SS_Label_43.3.5 Hockey stick
@@ -79,8 +74,8 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
// the 3rd parameter allows for a minimum recruitment level
{
steepness = SR_parm_work(2);
- temp = SR_parm_work(3) * Recr_virgin_adj + SSB_curr_adj / (steepness * SSB_virgin_adj) * (Recr_virgin_adj - SR_parm_work(3) * Recr_virgin_adj); // linear decrease below steepness*SSB_virgin_adj
- NewRecruits = Join_Fxn(0.0 * SSB_virgin_adj, SSB_virgin_adj, steepness * SSB_virgin_adj, SSB_curr_adj, temp, Recr_virgin_adj);
+ temp = SR_parm_work(3) * Recr_virgin_use + SSB_curr_adj / (steepness * SSB_virgin_use) * (Recr_virgin_use - SR_parm_work(3) * Recr_virgin_use); // linear decrease below steepness*SSB_virgin_use
+ NewRecruits = Join_Fxn(0.0 * SSB_virgin_use, SSB_virgin_use, steepness * SSB_virgin_use, SSB_curr_adj, temp, Recr_virgin_use);
break;
}
@@ -88,31 +83,32 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 6: //Beverton-Holt constrained
{
steepness = SR_parm_work(2);
- alpha = 4.0 * steepness * Recr_virgin / (5. * steepness - 1.);
- beta = (SSB_virgin_adj * (1. - steepness)) / (5. * steepness - 1.);
- if (SSB_curr_adj > SSB_virgin_adj)
+// dvariable SPR = SSB_virgin_use / Recr_virgin;
+// alpha = ((4.0 * steepness) / (1. - steepness)) / SPR ;
+// beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
+ if (SSB_curr_adj > SSB_virgin_use)
{
- SSB_BH1 = SSB_virgin_adj;
+ SSB_BH1 = SSB_virgin_use;
}
else
{
SSB_BH1 = SSB_curr_adj;
}
- NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_BH1) / (SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_BH1);
+ NewRecruits = (4. * steepness * Recr_virgin_use * SSB_BH1) / (SSB_virgin_use * (1. - steepness) + (5. * steepness - 1.) * SSB_BH1);
break;
}
// SS_Label_43.3.7 survival based
case 7: // survival based, so constrained such that recruits cannot exceed fecundity
{
- // PPR_0=SSB_virgin_adj/Recr_virgin_adj; // pups per recruit at virgin
+ // PPR_0=SSB_virgin_use/Recr_virgin_use; // pups per recruit at virgin
// Surv_0=1./PPR_0; // recruits per pup at virgin
- // Pups_0=SSB_virgin_adj; // total population fecundity is the number of pups produced
+ // Pups_0=SSB_virgin_use; // total population fecundity is the number of pups produced
// Sfrac=SR_parm(2);
- SRZ_0 = log(1.0 / (SSB_virgin_adj / Recr_virgin_adj));
+ SRZ_0 = log(1.0 / (SSB_virgin_use / Recr_virgin_use));
steepness = SR_parm_work(2);
srz_min = SRZ_0 * (1.0 - steepness);
- SRZ_surv = mfexp((1. - pow((SSB_curr_adj / SSB_virgin_adj), SR_parm_work(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
+ SRZ_surv = mfexp((1. - pow((SSB_curr_adj / SSB_virgin_use), SR_parm_work(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
NewRecruits = SSB_curr_adj * SRZ_surv;
exp_rec(y, 1) = NewRecruits; // expected arithmetic mean recruitment
// SS_Label_43.3.7.1 Do variation in recruitment by adjusting survival
@@ -149,8 +145,8 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
Shepherd_c2 = pow(0.2, SR_parm_work(3));
Hupper = 1.0 / (5.0 * Shepherd_c2);
steepness = 0.2 + (SR_parm_work(2) - 0.2) / (0.8) * (Hupper - 0.2);
- temp = (SSB_curr_adj) / (SSB_virgin_adj);
- NewRecruits = (5. * steepness * Recr_virgin_adj * (1. - Shepherd_c2) * temp) /
+ temp = (SSB_curr_adj) / (SSB_virgin_use);
+ NewRecruits = (5. * steepness * Recr_virgin_use * (1. - Shepherd_c2) * temp) /
(1.0 - 5.0 * steepness * Shepherd_c2 + (5. * steepness - 1.) * pow(temp, Shepherd_c));
break;
}
@@ -160,19 +156,20 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
{
steepness = SR_parm_work(2);
dvariable RkrPower = SR_parm_work(3);
- temp = SSB_curr_adj / SSB_virgin_adj;
+ temp = SSB_curr_adj / SSB_virgin_use;
temp2 = posfun(1.0 - temp, 0.0000001, temp3);
temp = 1.0 - temp2; // Rick's new line to stabilize recruitment at R0 if B>B0
dvariable RkrTop = log(5.0 * steepness) * pow(temp2, RkrPower) / pow(0.8, RkrPower);
- NewRecruits = Recr_virgin_adj * temp * mfexp(RkrTop);
+ NewRecruits = Recr_virgin_use * temp * mfexp(RkrTop);
break;
}
+
}
RETURN_ARRAYS_DECREMENT();
return NewRecruits;
} // end spawner_recruitment
-FUNCTION void apply_recdev(prevariable& NewRecruits, const prevariable& Recr_virgin_adj)
+FUNCTION void apply_recdev(prevariable& NewRecruits, const prevariable& Recr_virgin_use)
{
RETURN_ARRAYS_INCREMENT();
// SS_Label_43.4 For non-survival based SRR, get recruitment deviations by adjusting recruitment itself
@@ -199,7 +196,7 @@ FUNCTION void apply_recdev(prevariable& NewRecruits, const prevariable& Recr_vir
{
if (do_recdev >= 3)
{
- NewRecruits = Recr_virgin_adj * mfexp(recdev(y)); // recruitment deviation
+ NewRecruits = Recr_virgin_use * mfexp(recdev(y)); // recruitment deviation
}
else if (SR_fxn != 7)
{
@@ -231,7 +228,7 @@ FUNCTION void apply_recdev(prevariable& NewRecruits, const prevariable& Recr_vir
}
case 2: // use multiplier of R0
{
- exp_rec(y, 2) = Recr_virgin_adj * Fcast_Loop_Control(4); // apply fcast multiplier to the virgin recruitment
+ exp_rec(y, 2) = Recr_virgin_use * Fcast_Loop_Control(4); // apply fcast multiplier to the virgin recruitment
NewRecruits = exp_rec(y, 2);
if (SR_fxn != 4)
NewRecruits *= mfexp(-biasadj(y) * half_sigmaRsq); // bias adjustment
@@ -270,8 +267,8 @@ FUNCTION void apply_recdev(prevariable& NewRecruits, const prevariable& Recr_vir
//********************************************************************
/* SS_Label_FUNCTION 44 Equil_Spawn_Recr_Fxn */
// SPAWN-RECR: function Equil_Spawn_Recr_Fxn
-FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prevariable& SRparm3,
- const prevariable& SSB_virgin, const prevariable& Recr_virgin, const prevariable& SPR_temp)
+FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
+ const prevariable& SSB_virgin_use, const prevariable& Recr_virgin_use, const prevariable& SSBpR_current)
{
RETURN_ARRAYS_INCREMENT();
dvar_vector Equil_Spawn_Recr_Calc(1, 2); // values to return 1 is B_equil, 2 is R_equil
@@ -286,8 +283,8 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
dvariable srz_min;
dvariable SRZ_surv;
- steepness = SRparm2; // common usage but some different
- // SS_Label_44.1 calc equilibrium SpawnBio and Recruitment from input SPR_temp, which is spawning biomass per recruit at some given F level
+ steepness = SRparm(2); // common usage but some different
+ // SS_Label_44.1 calc equilibrium SpawnBio and Recruitment from input SSBpR_current, which is spawning biomass per recruit at some given F level
switch (SR_fxn)
{
case 1: // previous placement for B-H constrained
@@ -296,60 +293,79 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
write_message (FATAL, 0); // EXIT!
break;
}
- // SS_Label_44.1.1 Beverton-Holt with flattop beyond Bzero
- case 6: //Beverton-Holt
- {
- alpha = 4.0 * steepness * Recr_virgin / (5. * steepness - 1.);
- beta = (SSB_virgin * (1. - steepness)) / (5. * steepness - 1.);
- B_equil = alpha * SPR_temp - beta;
- B_equil = posfun(B_equil, 0.0001, temp);
- R_equil = (4. * steepness * Recr_virgin * B_equil) / (SSB_virgin * (1. - steepness) + (5. * steepness - 1.) * B_equil);
- break;
- }
+
// SS_Label_44.1.2 Ricker
case 2: // Ricker
{
- B_equil = SSB_virgin * (1. + (log(Recr_virgin / SSB_virgin) + log(SPR_temp)) / steepness);
- R_equil = Recr_virgin * B_equil / SSB_virgin * mfexp(steepness * (1. - B_equil / SSB_virgin));
+ B_equil = SSB_virgin_use * (1. + (log(Recr_virgin_use / SSB_virgin_use) + log(SSBpR_current)) / steepness);
+ R_equil = Recr_virgin_use * B_equil / SSB_virgin_use * mfexp(steepness * (1. - B_equil / SSB_virgin_use));
break;
}
- // SS_Label_44.1.3 Beverton-Holt
+ // SS_Label_44.1.1 Beverton-Holt
+ case 6: //Beverton-Holt with flattop beyond Bzero, but no flattop in equil calcs
+ {
+ }
+ // SS_Label_44.1.3 Beverton-Holt
case 3: // same as case 6
{
- alpha = 4.0 * steepness * Recr_virgin / (5. * steepness - 1.);
- beta = (SSB_virgin * (1. - steepness)) / (5. * steepness - 1.);
- B_equil = alpha * SPR_temp - beta;
+ // from WHAM per Tim Miller:
+ // WHAM based on R = A*S/(1+B*S)
+ // 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)));
+
+ // SS3 previously used alternative formulation: R = A*S/(B+S)
+ // converting SS3 to align with WHAM
+ alpha = 4.0 * steepness / (SSBpR_virgin_adj * (1. - steepness));
+ beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin_use);
+ // " h " << steepness << " derive " << alpha * SSBpR_virgin / (4. + alpha * SSBpR_virgin) << " " << endl;
+ // " R0 " << Recr_virgin_use << " derive " << 1. / beta * (alpha - 1./SSBpR_virgin) << endl;
+// report5 <<" SSB_unf "< 0)
- temp += Hermaphro_maleSPB * sum(MaleSPB(y, p));
+ if (Hermaphro_maleSSB > 0)
+ temp += Hermaphro_maleSSB * sum(MaleSSB(y, p));
SS2out << temp;
}
else
@@ -1529,7 +1529,7 @@ FUNCTION void write_bigoutput()
{
SS2out << SSB_pop_gp(y, p);
if (Hermaphro_Option != 0)
- SS2out << MaleSPB(y, p);
+ SS2out << MaleSSB(y, p);
}
else
{
@@ -1710,7 +1710,7 @@ FUNCTION void write_bigoutput()
if (k < 4)
k = 4;
// quantities to store summary statistics
- dvector rmse(1, k); // used in the SpBio, Index, Lencomp and Agecomp reports
+ dvector rmse(1, k); // used in the SSBio, Index, Lencomp and Agecomp reports
dvector Hrmse(1, k);
dvector Rrmse(1, k);
dvector n_rmse(1, k);
@@ -1779,6 +1779,7 @@ FUNCTION void write_bigoutput()
var /= (n_rmse(1) + 1.0e-09);
dvariable steepness = SR_parm(2);
+
SS2out << endl
<< pick_report_name(19);
SS2out << " Function: " << SR_fxn << " RecDev_method: " << do_recdev << " sum_recdev: " << sum_recdev << endl
@@ -1835,7 +1836,117 @@ FUNCTION void write_bigoutput()
}
SS2out << endl;
- SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev" << endl;
+ SS2out << endl << "#New_Expanded_Spawn_Recr_report" << endl << pick_report_name(19) << endl;
+ SS2out << "SR_Function: " << SR_fxn << endl;
+ SS2out << "N_SRparms: " << N_SRparm2 << endl;
+ SS2out << "#" << endl << "parm parm_label value phase" << endl;
+ for (int j = 1; j <=N_SRparm2; j++)
+ {
+ SS2out << j << " " << ParmLabel(firstSRparm + j) << " " << SR_parm(j) << " " << SR_parm_PH(j);
+ if (SR_parm_timevary(j) > 0 && j <= 4 ) // timevary SRparm exists
+ {SS2out << " #_is_time_vary,_so_SRR_updates_base_SPR_annually";}
+ if (j == (N_SRparm2 - 1) && SR_parm_timevary(j) > 0) // timevary regime exists
+ {SS2out << " #_Regime_parameter_used_to_offset_from_SRR";}
+ SS2out << endl;
+ }
+
+ SS2out << "# " < Bmark_Yr(k))
+ {SS2out << "#_range_of_years_is_averaged,_so_reduces_standard_error_of_result;_do_this_only_when_timevarying_makes_necessary: " << k << " "<< k+1 << endl;}
+ }
+ SS2out << "SSBpR_unfished_benchmark: " << Mgmt_quant(1) / Mgmt_quant(4) << " #_based_on_averaging_biology_over_benchmark_year_range " << endl;
+ SS2out << "Bmsy/Bzero: "<< Bmsy / SSB_virgin << " # using styr bio for Bzero" << endl;
+ SS2out << "Bmsy/Bunf: "<< Bmsy / Mgmt_quant(1) << " # using MSY's averaged bio for Bunf" << endl;
+
+ SS2out << "#" << endl << "RecDev_method: " << do_recdev << endl << "sum_recdev: " << sum_recdev << endl << "recr_logL: " << recr_like << endl;
+ SS2out << "main_recdev:start_end: " << recdev_start << " " << recdev_end << endl
+ << "breakpoints_for_bias_adjustment_ramp: " < (0.01 + 2.0 * square(rmse(1)) / temp))
+ {
+ warnstream << "Main recdev biasadj is >2 times ratio of rmse to sigmaR";
+ SS2out << " # " << warnstream.str() ;
+ write_message (WARN, 0);
+ }
+ }
+ SS2out << endl;
+
+ SS2out << "early " << n_rmse(3) << " " << rmse(3) << " " << square(rmse(3)) / temp << " " << rmse(4);
+ if (wrote_bigreport == 0) // first time writing bigreport
+ {
+ if (rmse(3) < 0.5 * sigmaR && rmse(4) > (0.01 + 2.0 * square(rmse(3)) / temp))
+ {
+ warnstream << "Early recdev biasadj is >2 times ratio of rmse to sigmaR";
+ SS2out << " # " << warnstream.str();
+ write_message (WARN, 0);
+ }
+ }
+ SS2out << endl << "#" << endl << "Initial_equilibrium: " << init_equ_steepness << " # 0/1_to_use_spawner-recruitment_in_initial_equ_recruitment_calculation" << endl << "#" << endl;
+ if (SR_fxn == 10) SS2out << "#_Note:_h_curr_and_R0_curr_are_for_info_only;_calculated_from_alpha_beta_and_current_SPR0" << endl;
+ SS2out << "#_columns_with_P_will_show_time_vary_SR_parameters" << endl << "#" << endl;
+ SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SPR0_curr ";
+ if(SR_fxn == 10)
+ {SS2out << "h_curr R0_curr ";}
+ else
+ {SS2out << "NA1 NA2 ";}
+;
+ for (j = 1; j <= N_SRparm2; j++)
+ {SS2out << "P" << j << " ";}
+ SS2out << endl;
SS2out << "S/Rcurve " << SSB_virgin << " " << Recr_virgin << endl;
y = styr - 2;
SS2out << "Virg " << SSB_yr(y) << " " << exp_rec(y) << " - " << 0.0 << " Virg " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " 0.0 " << endl;
@@ -1886,7 +1997,16 @@ FUNCTION void write_bigoutput()
{
SS2out << " _ _ Fixed";
}
- SS2out << endl;
+ dvariable SPR_curr = Smry_Table(y, 11) / Recr_virgin;
+ SS2out << " " << SPR_curr << " ";
+ if (SR_fxn == 10)
+ {
+ SS2out << alpha * SPR_curr / (4. + alpha * SPR_curr) << " "; // steepness with current SPR
+ SS2out << 1. / beta * (alpha - (1. / SPR_curr)) << " "; // R0 with current SPR
+ }
+ else
+ {SS2out << " - - ";}
+ SS2out << SR_parm_byyr(y)(1,N_SRparm2) << endl;
}
// REPORT_KEYWORD SPAWN_RECR_CURVE
@@ -4613,33 +4733,11 @@ FUNCTION void SPR_profile()
for (s = 1; s <= nseas; s++)
{
t = styr - 3 * nseas + s - 1;
-
- if (MG_active(2) > 0 || MG_active(3) > 0 || save_for_report > 0 || WTage_rd > 0)
- {
- subseas = 1;
- ALK_idx = (s - 1) * N_subseas + subseas; // for midseason
- Make_AgeLength_Key(s, subseas); // for begin season
- subseas = mid_subseas;
- ALK_idx = (s - 1) * N_subseas + subseas; // for midseason
- Make_AgeLength_Key(s, subseas); // for midseason
- if (s == spawn_seas)
- {
- subseas = spawn_subseas;
- if (spawn_subseas != 1 && spawn_subseas != mid_subseas)
- {
- //don't call get_growth3(subseas) because using an average ave_size
- Make_AgeLength_Key(s, subseas); // spawn subseas
- }
- get_mat_fec();
- }
- }
-
- for (g = 1; g <= gmorph; g++)
- if (use_morph(g) > 0)
- {
- ALK_idx = (s - 1) * N_subseas + mid_subseas; // for midseason
- Make_FishSelex();
- }
+ Wt_Age_beg(s) = Wt_Age_t(t, 0); // used for smrybio
+ Wt_Age_mid(s) = Wt_Age_t(t, -1);
+ if (s == spawn_seas)
+ fec = Wt_Age_t(t, -2);
+// report5 << 0 << " y: " << y << " updated_Repro_output spr/ypr: " << fec(1) << endl;
}
SS2out << "SPRloop Iter Bycatch Fmult F_std SPR YPR_dead YPR_dead*Recr YPR_ret*Recr Revenue Cost Profit SSB Recruits SSB/Bzero Tot_Catch ";
@@ -4792,8 +4890,8 @@ FUNCTION void SPR_profile()
Do_Equil_Calc(equ_Recr);
// SPAWN-RECR: calc equil spawn-recr in the SPR loop
- SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Btgt_prof = Equ_SpawnRecr_Result(1);
Btgt_prof_rec = Equ_SpawnRecr_Result(2);
if (Btgt_prof < 0.001 || Btgt_prof_rec < 0.001)
@@ -4893,7 +4991,7 @@ FUNCTION void Global_MSY()
// GLOBAL_MSY with knife-edge age selection, then slot-age selection
SS2out << endl
<< pick_report_name(55) << endl;
- y = styr - 3; // stores the averaged
+ y = styr - 3; // stores the averaged biology and selectivity, etc. from benchmark
yz = y;
bio_yr = y;
eq_yr = y;
@@ -4972,7 +5070,9 @@ FUNCTION void Global_MSY()
SS2out << "Actual ";
show_MSY = 2; // invokes just brief output in benchmark
did_MSY = 0;
+// report5 << 0 << " y: " << y << " updated_Repro_output global_1: " << fec(1) << endl;
Get_Benchmarks(show_MSY);
+// report5 << 0 << " y: " << y << " updated_Repro_output global_2: " << fec(1) << endl;
did_MSY = 0;
}
}
diff --git a/SS_write_ssnew.tpl b/SS_write_ssnew.tpl
index 38a9a5ec..82e966aa 100644
--- a/SS_write_ssnew.tpl
+++ b/SS_write_ssnew.tpl
@@ -1548,7 +1548,7 @@ FUNCTION void write_nucontrol()
NuStart << final_conv << " # final convergence criteria (e.g. 1.0e-04) " << endl;
NuStart << retro_yr - endyr << " # retrospective year relative to end year (e.g. -4)" << endl;
NuStart << Smry_Age << " # min age for calc of summary biomass" << endl;
- NuStart << depletion_basis_rd << " # Depletion basis: denom is: 0=skip; 1=X*SPBvirgin; 2=X*SPBmsy; 3=X*SPB_styr; 4=X*SPB_endyr; 5=X*dyn_Bzero; values>=11 invoke N multiyr with 10s & 100s digit; append .1 to invoke log(ratio); e.g. 122.1 produces log(12 year trailing average of B/Bmsy)" << endl;
+ NuStart << depletion_basis_rd << " # Depletion basis: denom is: 0=skip; 1=X*SSBvirgin; 2=X*SSBmsy; 3=X*SSB_styr; 4=X*SSB_endyr; 5=X*dyn_Bzero; values>=11 invoke N multiyr with 10s & 100s digit; append .1 to invoke log(ratio); e.g. 122.1 produces log(12 yr trailing average of B/Bmsy)" << endl;
NuStart << depletion_level << " # Fraction (X) for Depletion denominator (e.g. 0.4)" << endl;
NuStart << SPR_reporting << " # SPR_report_basis: 0=skip; 1=(1-SPR)/(1-SPR_tgt); 2=(1-SPR)/(1-SPR_MSY); 3=(1-SPR)/(1-SPR_Btarget); 4=rawSPR" << endl;
NuStart << F_reporting << " # F_std_reporting_units: 0=skip; 1=exploitation(Bio); 2=exploitation(Num); 3=sum(Apical_F's); 4=mean F for range of ages (numbers weighted); 5=unweighted mean F for range of ages" << endl;
@@ -1644,6 +1644,7 @@ FUNCTION void write_nucontrol()
NuFore << H4010_top_rd << " # Control rule inflection for constant F (as frac of Bzero, e.g. 0.40); must be > control rule cutoff, or set to -1 to use Bmsy/SSB_unf " << endl;
NuFore << H4010_bot << " # Control rule cutoff for no F (as frac of Bzero, e.g. 0.10) " << endl;
NuFore << H4010_scale_rd << " # Buffer: enter Control rule target as fraction of Flimit (e.g. 0.75), negative value invokes list of [year, scalar] with filling from year to YrMax " << endl;
+ NuFore << "# Also see HCR_anchor below" << endl;
if (H4010_scale_rd < 0)
{
j = H4010_scale_vec_rd.size() - 1;
@@ -1665,7 +1666,7 @@ FUNCTION void write_nucontrol()
{
NuFore << Fcast_Loop_Control(4) << " # multiplier on base recruitment " << endl;
}
- NuFore << Fcast_Loop_Control(5) << " # not used" << endl << "#" << endl;
+ NuFore << Fcast_Loop_Control(5) << " # HCR_anchor: 0 or 1 uses unfished benchmark SSB (old hardwired approach), 2 = virgin SSB" << endl << "#" << endl;
NuFore << Fcast_Cap_FirstYear << " # FirstYear for caps and allocations (should be after years with fixed inputs) " << endl;
@@ -1941,7 +1942,7 @@ FUNCTION void write_nucontrol()
if (Hermaphro_Option != 0)
{
report4 << Hermaphro_seas_rd << " # Hermaphro_season.first_age (seas=-1 means all seasons; first_age must be 0 to 9)" << endl
- << Hermaphro_maleSPB << " # fraction_of_maleSSB_added_to_total_SSB " << endl;
+ << Hermaphro_maleSSB << " # fraction_of_maleSSB_added_to_total_SSB " << endl;
}
report4 << MGparm_def << " #_parameter_offset_approach for M, G, CV_G: 1- direct, no offset**; 2- male=fem_parm*exp(male_parm); 3: male=female*exp(parm) then old=young*exp(parm)" << endl;
@@ -2136,7 +2137,11 @@ FUNCTION void write_nucontrol()
report4 << "#" << endl;
report4 << SR_fxn << " #_Spawner-Recruitment; Options: 1=NA; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=survival_3Parm; 8=Shepherd_3Parm; 9=RickerPower_3parm" << endl;
report4 << init_equ_steepness << " # 0/1 to use steepness in initial equ recruitment calculation" << endl;
- report4 << sigmaR_dendep << " # future feature: 0/1 to make realized sigmaR a function of SR curvature" << endl;
+ report4 << SR_update_SSBpR0_rd << " # SR_update_SSBpR0" << endl <<
+ "# 0 - OK, but only if no timevary biology or SR parm" << endl <<
+ "# 1 - best: update SSBpR0 for benchmark and for time series only if SRparm R0 or h (not regime) is set to have time-varying property" << endl <<
+ "# 2 - incorrect (old, incorrect SS3 approach): always update SSBpR0 for benchmark's use of spawner-recruitment, but only for the time series if there is a timevary SR parm" << endl <<
+ "# 3 - option: do not update SSBpR0 (do keep start year SSBpR0), even if R0 or h is set to have time-varying property" << endl << "#" << endl;
report4 << "#_ LO HI INIT PRIOR PR_SD PR_type PHASE env-var use_dev dev_mnyr dev_mxyr dev_PH Block Blk_Fxn # parm_name" << endl;
report4.unsetf(std::ios_base::fixed);
report4.unsetf(std::ios_base::floatfield);
diff --git a/StockSynthesis.code-workspace b/StockSynthesis.code-workspace
index 2bfc56e8..14159ac5 100644
--- a/StockSynthesis.code-workspace
+++ b/StockSynthesis.code-workspace
@@ -10,7 +10,9 @@
"\"*.extension\":": "\"tpl\"",
"*.htp": "c",
"ostream": "c",
- "iosfwd": "c"
+ "xlocale": "c",
+ "iosfwd": "c",
+ "cmath": "c"
},
"explorer.excludeGitIgnore": true
}