diff --git a/.github/workflows/run-ss3-no-est.yml b/.github/workflows/run-ss3-no-est.yml
index c77bf619..ccd0e537 100644
--- a/.github/workflows/run-ss3-no-est.yml
+++ b/.github/workflows/run-ss3-no-est.yml
@@ -111,14 +111,15 @@ jobs:
on.exit(system(paste0("cd ", wd)))
# delete old starter files, rename forecast.ss_new and starter.ss_new files
file.remove(file.path(dir, "starter.ss"))
- file.remove(file.path("forecast.ss"))
+ file.remove(file.path(dir, "forecast.ss"))
file.rename(file.path(dir, "starter.ss_new"), file.path(dir,"starter.ss"))
file.rename(file.path(dir, "forecast.ss_new"), file.path(dir,"forecast.ss"))
# rename control and data files to standardized names (from the starter files)
start <- readLines(file.path(dir, "starter.ss"))
- first_val_line <- grep("0=use init values in control file", start, fixed = TRUE)
- datname <- start[first_val_line-2]
- ctlname <- start[first_val_line-1]
+ datname <- start[grep("#_datfile", start, fixed = TRUE)]
+ datname <- gsub(" #_datfile", "", datname)
+ ctlname <- start[grep("#_ctlfile", start, fixed = TRUE)]
+ ctlname <- gsub(" #_ctlfile", "", ctlname)
print(datname)
print(ctlname)
file.remove(file.path(dir, datname))
diff --git a/.github/workflows/run-ss3-with-est.yml b/.github/workflows/run-ss3-with-est.yml
index cee49885..2225df50 100644
--- a/.github/workflows/run-ss3-with-est.yml
+++ b/.github/workflows/run-ss3-with-est.yml
@@ -1,5 +1,4 @@
# Build SS3 and run test models with estimation and hessian
-# Only runs manual run or if run-ss3-no-est.yml runs successfully
name: run-ss3-with-est
on:
workflow_dispatch:
@@ -8,21 +7,20 @@ on:
# - '**.tpl'
# branches:
# - main
- workflow_run:
- workflows: ["run-ss3-no-est"]
- types:
- - completed
- # pull_request:
- # types: ['opened', 'edited', 'reopened', 'synchronize', 'ready_for_review']
- # paths:
- # - '**.tpl'
- # branches:
- # - main
+ # workflow_run:
+ # workflows: ["run-ss3-no-est"]
+ # types:
+ # - completed
+ pull_request:
+ paths:
+ - '**.tpl'
+ branches:
+ - main
# Run fast running SS3 test models with estimation
jobs:
run-ss3-with-est:
- if: ${{ github.event_name == 'workflow_dispatch' || (github.event.pull_request.draft == 'false' && github.event.workflow_run.conclusion == 'success') }}
+ # if: ${{ github.event_name == 'workflow_dispatch' || (github.event.pull_request.draft == 'false'}}
runs-on: ubuntu-latest
env:
R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
diff --git a/.github/workflows/test-r4ss-with-ss3.yml b/.github/workflows/test-r4ss-with-ss3.yml
index e1d9e618..fcf49d09 100644
--- a/.github/workflows/test-r4ss-with-ss3.yml
+++ b/.github/workflows/test-r4ss-with-ss3.yml
@@ -52,7 +52,7 @@ jobs:
run: Rscript -e 'install.packages(c("remotes","parallely", "furrr", "future"))'
- name: Install r4ss
- run: Rscript -e 'remotes::install_github("r4ss/r4ss")'
+ run: Rscript -e 'remotes::install_github("r4ss/r4ss@expanded_SR")'
# - name: Get admb and put in path, linux
# run: |
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 cc9b6f97..87683d9d 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -4,6 +4,14 @@
// 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 is calculated in popdyn using the start year biology
+// SSBpR_virgin used to get alpha in equil_spawn_recr B-H
+
FUNCTION void setup_Benchmark() // and forecast
{
// SS_Label_Info_7.5 #Get averages from selected years to use in forecasts
@@ -113,12 +121,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
{
@@ -412,7 +415,7 @@ FUNCTION void setup_Benchmark() // and forecast
}
// recr_dist_unf is accumulated while doing the time_series
// then its mean is calculated in Get_Benchmarks and assigned to recr_dist
- // the SR_parm_bench is calculated from Bmark_yrs 9-10 in benchmark code using values stored in SR_parm_byyr
+ // the SRparm_bench is calculated from Bmark_yrs 9-10 in benchmark code using values stored in SRparm_byyr
}
// calc average selectivity to use in equil; store in styr-3
@@ -536,7 +539,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 +572,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 +645,7 @@ 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
if (s == spawn_seas)
{
{
@@ -730,34 +733,126 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
for (j = 1; j <= N_SRparm2; j++)
{
- if (SR_parm_timevary(j) == 0)
+ if (SRparm_timevary(j) == 0)
{
- SR_parm_work(j) = SR_parm(j);
+ SRparm_bench(j) = SRparm(j);
}
else
{
temp = 0.;
for (int y = Bmark_Yr(9); y <= Bmark_Yr(10); y++)
{
- temp += SR_parm_byyr(y, j);
+ temp += SRparm_byyr(y, j);
}
- SR_parm_work(j) = temp / (Bmark_Yr(10) - Bmark_Yr(9) + 1.);
+ SRparm_bench(j) = temp / (Bmark_Yr(10) - Bmark_Yr(9) + 1.);
}
}
+ SRparm_bench = SRparm_work;
+
+ /*
+ Flags:
+ timevary_MG_firstyr == YrMax // means that no biology is time-varying
+ timevary_SRparm_first > 0 // means that R0 or h (i.e. any except regime, sigmaR, autocorr) is time-varying, so SSBpR0 gets updated for time series and for bench
+ timevary_bio_4SRR is new user selected flag: 0 for legacy, vs 1 for improved use of timevary biology in SRR calcs
+
+ Legacy approach:
+ SSBpR0 set at start year using start year biology
+ SSBpR0 is not itself saved; instead R0_4_SRR and SSB0_4_SRR are saved and passed to the spawn_recruit functions
+ SSBpR0 updated during time series if there is time-varying R0, but does not call equil_spawn_recr_calc
+ SSBpR with benchmark biology used in benchmark calculations
+ SSBpR0 for benchmark always uses bench biology, which is incorrect
+ SSB_bench (aka SSB_unf) does not call equil_spawn_recr_calc, it is just R * SSBpR, so is incomplete accounting for effect of timevary bio
+ Btgttgt is a fraction of SSB_bench (no options)
+ Btgttgt2 can be fraction of SSB_MSY or fraction of SSB_virgin, but not SSB_bench
+ HCR inflection is a fraction of SSB_bench (no options)
+ depletion basis is user-selected as SSB_virgin or SSB_msy
+ SSB_msy uses equil_spawn_recr_calc in its creation, but SSB_bench does not. So they are inconsistent
+ none of the above is an issue if there is no timevarying biology
+ -------
+ Improved approach
+ SSBpR0 set at start year using start year biology
+ SSBpR0 updated during time series if there is time-varying R0; PLUS NEW: equil_spawn_recr_calc called to get new equilibrium R0, SSB0
+ SSBpR0 for benchmark stays at virgin unless timevary_bio_4SRR == 0, or if timevary_SRparm_first > 0
+ Btgttgt can now use either frac*SSB_bench or frac*SSB_virgin by using the existing flag for depletion basis
+ Btgttgt2 can be fraction of SSB_MSY, of SSB_virgin, or of SSB_bench
+ HCR inflection adds option to use SSB_virgin or SSB_bench
+ depletion adds option to use SSB_bench
+ */
+ Recr_unf = Recr_virgin; // default
+ SSB0_4_SRR = SSB_virgin; // default
+ R0_4_SRR = Recr_virgin;
Fishon = 0;
- Recr_unf = mfexp(SR_parm_work(1));
- Do_Equil_Calc(Recr_unf);
+ SSBpR_Calc(Recr_unf); // this returns SSB_equil using benchmark biology
+ // provides basis for values needed below
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;
+ SSBpR_bench = SSB_equil / Recr_unf;
+
+ if(timevary_SRparm_first == 0) // no timevary SRR parms
+ {
+ if( timevary_MG_firstyr == YrMax && WTage_rd == 0) // no time-varying biology
+ {
+ R0_4_SRR = Recr_virgin;
+ SSB0_4_SRR = SSB_virgin;
+ SSB_unf = SSB_virgin;
+ }
+ else // there is time-varying biology
+ {
+ R0_4_SRR = Recr_virgin; // same as Recr_virgin because no timevary SRparms
+ if( timevary_bio_4SRR == 0) // legacy approach; this switch is read from starter.ss
+ {
+ SSB_unf = SSB_equil;
+ SSB0_4_SRR = SSB_equil; // this is inaccurate legacy, as it moves equil off the SRR, rather than along the SRR
+ }
+ else
+ {
+ // get new equilibrium point using original SRR and SSBpR_bench
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SRparm_bench, SSB_virgin, Recr_virgin, SSBpR_bench); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSB_unf = Equ_SpawnRecr_Result(1);
+ Recr_unf = Equ_SpawnRecr_Result(2);
+ if (show_MSY == 1) report5 << " use virgin SSBpR0 in SRR - SSB: " << SSB_virgin << " Recr: " << Recr_virgin << " SPR: " << SSB_virgin / Recr_virgin << " bench SPR: " << SSBpR_bench << " new equil: " << Equ_SpawnRecr_Result << endl;
+ }
+ SSBpR_bench = SSB_unf / Recr_unf;
+ }
+ }
+ else // there are timevary SRR parms; use same code regardless of timevary biology. Legacy approach does not include new equilibrium
+ {
+ Recr_unf = mfexp(SRparm_bench(1)); // R0 to be used
+ // note that steepness will get updated when SRparm_bench is used in Equ_SpawnRecr_Result
+ SSBpR_Calc(Recr_unf); // this returns SSB_equil using benchmark biology
+ SSB_unf = SSB_equil;
+ SSBpR_bench = SSB_equil / Recr_unf;
+ if( timevary_bio_4SRR == 0) // legacy approach; this switch is read from starter.ss
+ {
+ R0_4_SRR = Recr_unf;
+ SSB0_4_SRR = SSB_equil; // this is legacy, but incorrect, as it moves equil off the SRR, rather than along the SRR
+ }
+ else // use improved approach with updated SRparms and benchmark biology
+ {
+ // get new equilibrium point for the benchmark SRR
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SRparm_bench, SSB_equil, Recr_unf, SSBpR_bench); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSB_unf = Equ_SpawnRecr_Result(1);
+ SSB0_4_SRR = Equ_SpawnRecr_Result(1);
+ Recr_unf = Equ_SpawnRecr_Result(2);
+ R0_4_SRR = Equ_SpawnRecr_Result(2);
+
+ if (show_MSY == 1) report5 << " use bench SSBpR0 in SRR - SSB: " << SSB_unf << " Recr: " << Recr_unf << " SPR: " << SSBpR_bench << " new equil: " << Equ_SpawnRecr_Result << endl;
+ }
+ }
+
if (show_MSY == 1)
- report5 << "Repro_output_by_age_for_morph_1: " << fec(1) << endl;
- Mgmt_quant(1) = SSB_unf;
- Mgmt_quant(2) = totbio;
- Mgmt_quant(3) = smrybio;
- Mgmt_quant(4) = Recr_unf;
+ {
+ SRparm_bench(N_SRparm2 + 1) = SSB_unf;
+ Mgmt_quant(1) = SSB_unf;
+ Mgmt_quant(2) = totbio; // this is calculated in Do_Equil_Calc
+ Mgmt_quant(3) = smrybio;
+ Mgmt_quant(4) = Recr_unf;
+ report5 << "SRparms for benchmark: " << SRparm_bench << endl
+ << "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl << endl;
+ Mgmt_quant(19) = SSB_unf; // placeholder for depletion denominator
+ Mgmt_quant(20) = SSB_unf; // placeholder to be replaced by SSB_HCR_infl
+ Mgmt_quant(21) = R0_4_SRR;
+ Mgmt_quant(22) = SSB0_4_SRR;
+ }
// find Fspr SS_Label_710
{
@@ -781,13 +876,13 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
last_calc = 0.;
Fchange = -4.0;
- equ_Recr = 1.0;
+ equ_Recr = 1.0; // so calls to Do_Equil_Calc will return values of SSBpR
Fishon = 0;
dvariable SPR_target100;
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_Calc(equ_Recr); // where equ_Recr has been set to 1.0
+ SSBpR_unf = SSB_equil / equ_Recr; // 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)
{
@@ -832,8 +927,8 @@ 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
+ SSBpR_Calc(equ_Recr);
+ yld1(ii) = 100. * SSB_equil / SSBpR_unf; // spawning potential ratio
}
SPR_actual = yld1(1); // spawning potential ratio
@@ -902,7 +997,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(SRparm_bench, SSB0_4_SRR, R0_4_SRR, 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);
@@ -943,7 +1039,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
// else Hrate for bycatch fleets already set
}
- Do_Equil_Calc(equ_Recr);
+ SSBpR_Calc(equ_Recr);
F01_origin = YPR_opt / Fmult;
BTGT_target = 0.1; // now relative to Bmark
@@ -975,7 +1071,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
} // else Hrate for bycatch fleets set above
}
- Do_Equil_Calc(equ_Recr);
+ SSBpR_Calc(equ_Recr);
yld1(ii) = YPR_opt;
}
@@ -985,7 +1081,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 +1112,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(SRparm_bench, SSB0_4_SRR, R0_4_SRR, 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 +1124,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;
@@ -1041,17 +1137,27 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
else // find F giving Btarget SS_Label_720
{
// ******************************************************
- if (show_MSY == 1)
+
+ if (depletion_basis == 1)
+ {Btgttgt = BTGT_target * SSB_virgin;}
+ else
+ {Btgttgt = BTGT_target * SSB_unf;} // current SS3 approach uses benchmark biology
+
+ 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";
- for (p = 1; p <= pop; p++)
- for (gp = 1; gp <= N_GP; gp++)
- {
- report5 << " SSB_Area:" << p << "_GP:" << gp;
- }
- report5 << endl;
+ report5 << "#" << endl;
+ if (depletion_basis == 1)
+ {report5 << "Find_target_SSB as fraction: " << BTGT_target << " of Virgin SSB:" << SSB_virgin << endl;}
+ else
+ {report5 << "Find_target_SSB as fraction: " << BTGT_target << " of SSB_unf: "<< SSB_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++)
+ {
+ report5 << " SSB_Area:" << p << "_GP:" << gp;
+ }
+ report5 << endl;
}
F1(1) = log(1.0e-3);
@@ -1071,9 +1177,6 @@ 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
-
for (j = 0; j <= Nloops; j++) // loop find Btarget
{
if (fabs(Fchange) <= Closer2)
@@ -1107,11 +1210,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;
+ SSBpR_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(SRparm_bench, SSB0_4_SRR, R0_4_SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
yld1(ii) = Equ_SpawnRecr_Result(1);
}
@@ -1287,7 +1390,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
if (F_Method == 1)
{
- Fmax = (Btgt_Fmult + SPR_Fmult) * 0.5 * SR_parm_work(2) / 0.05;
+ Fmax = (Btgt_Fmult + SPR_Fmult) * 0.5 * SRparm_bench(2) / 0.05;
} // previously /0.18
F1(1) = -log(Fmax / Btgt_Fmult - 1.);
F2(1) = -log(Fmax / Btgt_Fmult - 1.);
@@ -1327,12 +1430,12 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
// else Hrate for bycatch fleets already set
}
- Do_Equil_Calc(equ_Recr);
+ SSBpR_Calc(equ_Recr);
// SPAWN-RECR: calc spawn-recr for MSY calcs; need to make area-specific
- MSY_SPR = SSB_equil / SPR_unfished;
- 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
- Bmsy = Equ_SpawnRecr_Result(1);
+ MSY_SPR = SSB_equil / SSBpR_unf;
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SRparm_bench, SSB0_4_SRR, R0_4_SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ Bmsy = Equ_SpawnRecr_Result(1); // with MSY set to SPR, not directly estimated
Recr_msy = Equ_SpawnRecr_Result(2);
yld1(1) = YPR_opt * Recr_msy;
YPR_msy_enc = YPR_enc;
@@ -1358,12 +1461,13 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
report5 << endl;
}
- Mgmt_quant(15) = yld1(1);
Mgmt_quant(12) = Bmsy;
Mgmt_quant(13) = MSY_SPR;
Mgmt_quant(14) = equ_F_std;
+ Mgmt_quant(15) = yld1(1);
Mgmt_quant(16) = YPR_ret * Recr_msy;
Mgmt_quant(17) = Bmsy / SSB_unf;
+ Mgmt_quant(18) = Recr_msy;
Vbio1_MSY = smrybio;
Vbio_MSY = totbio;
}
@@ -1420,12 +1524,12 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
} // else Hrate for bycatch fleets set above
}
- Do_Equil_Calc(equ_Recr);
+ SSBpR_Calc(equ_Recr);
// SPAWN-RECR: calc spawn-recr for MSY calcs; need to make area-specific
- MSY_SPR = SSB_equil / SPR_unfished;
- 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
- Bmsy = Equ_SpawnRecr_Result(1);
+ MSY_SPR = SSB_equil / SSBpR_unf;
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SRparm_bench, SSB0_4_SRR, R0_4_SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ Bmsy = Equ_SpawnRecr_Result(1); // MSY is directly estimated
Recr_msy = Equ_SpawnRecr_Result(2);
Profit = (PricePerF * YPR_val_vec) * Recr_msy - Cost;
if (Do_MSY == 2) // dead catch without excluded bycatch fleets
@@ -1496,12 +1600,14 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
YPR_msy_profit = YPR_msy_revenue - Cost;
MSY = yld1(1);
MSY_Fmult = Fmult;
- Mgmt_quant(15) = yld1(1);
Mgmt_quant(12) = Bmsy;
Mgmt_quant(13) = MSY_SPR;
Mgmt_quant(14) = equ_F_std;
+ Mgmt_quant(15) = yld1(1);
Mgmt_quant(16) = YPR_ret * Recr_msy;
Mgmt_quant(17) = Bmsy / SSB_unf;
+ Mgmt_quant(18) = Recr_msy;
+
Vbio1_MSY = smrybio;
Vbio_MSY = totbio;
@@ -1569,7 +1675,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
Fishon = 0;
- Do_Equil_Calc(equ_Recr);
+ SSBpR_Calc(equ_Recr);
report5 << "Equil_N_at_age_M_only_Recr_MSY" << endl
<< "Seas Area GP Sex subM" << age_vector << endl;
for (s = 1; s <= nseas; s++)
@@ -1675,11 +1781,16 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Btgttgt2 = Blim_frac * Bmsy;
Blim_report = value(Bmsy);
}
- else
+ else if (depletion_basis == 1)
{
Btgttgt2 = -Blim_frac * SSB_virgin;
Blim_report = value(SSB_virgin);
}
+ else
+ {
+ Btgttgt2 = -Blim_frac * SSB_unf;
+ Blim_report = value(SSB_unf);
+ }
for (j = 0; j <= Nloops; j++) // loop find Btarget
{
@@ -1714,11 +1825,11 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
// else Hrate for bycatch fleets already set
}
- Do_Equil_Calc(equ_Recr);
- SPR_Btgt2 = SSB_equil / SPR_unfished;
+ SSBpR_Calc(equ_Recr);
+ SPR_Btgt2 = SSB_equil / SSBpR_unf;
// 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
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SRparm_bench, SSB0_4_SRR, R0_4_SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
yld1(ii) = Equ_SpawnRecr_Result(1);
}
@@ -1787,9 +1898,9 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Btgt_Fmult2 = Fmult;
if (rundetail > 0 && mceval_counter == 0 && show_MSY == 1)
echoinput << "Calculated F_Blimit " << Btgt_Fmult2 << " " << Btgt2 / Blim_report << endl;
- Mgmt_quant(18) = Btgt2;
- Mgmt_quant(19) = equ_F_std;
- Mgmt_quant(20) = sum(equ_catch_fleet(2)) * Equ_SpawnRecr_Result(2);
+ Mgmt_quant(N_STD_Mgmt_Quant - 2) = Btgt2;
+ Mgmt_quant(N_STD_Mgmt_Quant - 1) = equ_F_std;
+ Mgmt_quant(N_STD_Mgmt_Quant) = sum(equ_catch_fleet(2)) * Equ_SpawnRecr_Result(2);
} // end finding F for Blimit
if (rundetail > 0 && mceval_counter == 0 && show_MSY == 1)
@@ -1800,8 +1911,9 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
{
report5 << "#" << endl
<< "Management_report" << endl;
- report5 << "Steepness_Recr_SSB_virgin(R0) " << SR_parm(2) << " " << Recr_virgin << " " << SSB_virgin << endl;
- report5 << "Steepness_Recr_SSB_benchmark " << SR_parm_work(2) << " " << Recr_unf << " " << SSB_unf << endl;
+ report5 << "Virgin: Steepness_Recr_SSB " << SRparm(2) << " " << Recr_virgin << " " << SSB_virgin << endl;
+ report5 << "Bench: Steepness_Recr_SSB " << SRparm_bench(2) << " " << R0_4_SRR << " " << SSB0_4_SRR << endl;
+ report5 << "unf : Steepness_Recr_SSB " << SRparm_bench(2) << " " << Recr_unf << " " << SSB_unf << endl;
report5 << "#" << endl
<< "Summary_age: " << Smry_Age << endl;
report5 << "#_Bmark_years: beg_bio, end_bio, beg_selex, end_selex, beg_relF, end_relF, beg_recr_dist, end_recr_dist, beg_SRparm, end_SRparm" << endl
@@ -1826,7 +1938,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
report5 << "ann_F " << Mgmt_quant(10) << endl;
report5 << "Exploit(Catch_dead/B_smry) " << YPR_spr_dead / Vbio1_spr << endl;
report5 << "Recruits " << Bspr_rec << endl;
- report5 << "SPBio " << Bspr << " " << Bspr / Bspr_rec << endl;
+ report5 << "SSBio " << Bspr << " " << Bspr / Bspr_rec << endl;
report5 << "Catch_encountered " << YPR_spr_enc * Bspr_rec << " " << YPR_spr_enc << endl;
report5 << "Catch_dead " << YPR_spr_dead * Bspr_rec << " " << YPR_spr_dead << endl;
report5 << "Catch_retain " << YPR_spr_ret * Bspr_rec << " " << YPR_spr_ret << endl;
@@ -1846,7 +1958,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
report5 << "ann_F " << Mgmt_quant(7) << endl;
report5 << "Exploit(Catch_dead/B_smry) " << YPR_Btgt_dead / Vbio1_Btgt << endl;
report5 << "Recruits@F0.1 " << Btgt_Rec << endl;
- report5 << "SPBio " << Btgt << " " << Btgt / Btgt_Rec << endl;
+ report5 << "SSBio " << Btgt << " " << Btgt / Btgt_Rec << endl;
report5 << "Catch_encountered " << YPR_Btgt_enc * Btgt_Rec << " " << YPR_Btgt_enc << endl;
report5 << "Catch_dead " << YPR_Btgt_dead * Btgt_Rec << " " << YPR_Btgt_dead << endl;
report5 << "Catch_retain " << YPR_Btgt_ret * Btgt_Rec << " " << YPR_Btgt_ret << endl;
@@ -1866,7 +1978,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
report5 << "ann_F " << Mgmt_quant(7) << endl;
report5 << "Exploit(Catch_dead/B_smry) " << YPR_Btgt_dead / Vbio1_Btgt << endl;
report5 << "Recruits " << Btgt_Rec << endl;
- report5 << "SPBio " << Btgt << " " << Btgt / Btgt_Rec << endl;
+ report5 << "SSBio " << Btgt << " " << Btgt / Btgt_Rec << endl;
report5 << "Catch_encountered " << YPR_Btgt_enc * Btgt_Rec << " " << YPR_Btgt_enc << endl;
report5 << "Catch_dead " << YPR_Btgt_dead * Btgt_Rec << " " << YPR_Btgt_dead << endl;
report5 << "Catch_retain " << YPR_Btgt_ret * Btgt_Rec << " " << YPR_Btgt_ret << endl;
@@ -1883,9 +1995,9 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
report5 << "ann_F " << Mgmt_quant(14) << endl;
report5 << "Exploit(Catch/Bsmry) " << MSY / (Vbio1_MSY * Recr_msy) << endl;
report5 << "Recruits@MSY " << Recr_msy << endl;
- report5 << "SPBmsy " << Bmsy << " " << Bmsy / Recr_msy << endl;
- report5 << "SPBmsy/SPB_virgin " << Bmsy / SSB_virgin << endl;
- report5 << "SPBmsy/SPB_unfished " << Bmsy / SSB_unf << endl;
+ report5 << "SSBmsy " << Bmsy << " " << Bmsy / Recr_msy << endl;
+ report5 << "SSBmsy/SSB_virgin " << Bmsy / SSB_virgin << endl;
+ report5 << "SSBmsy/SSB_unfished " << Bmsy / SSB_unf << endl;
report5 << "MSY_for_optimize " << MSY << " " << MSY / Recr_msy << endl;
report5 << "MSY_encountered " << YPR_msy_enc * Recr_msy << " " << YPR_msy_enc << endl;
report5 << "MSY_dead " << YPR_msy_dead * Recr_msy << " " << YPR_msy_dead << endl;
@@ -1911,6 +2023,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
<< " " << Vbio1_MSY * Recr_msy << " # " << endl;
}
write_bodywt = write_bodywt_save;
+// report5 << "Repro_output_by_age_for_morph_1_after_benchmark: " << fec(1) << endl;
} // end benchmarks
FUNCTION void Get_Forecast()
@@ -1922,8 +2035,6 @@ FUNCTION void Get_Forecast()
dvariable OFL_catch;
dvariable Fcast_Crash;
dvariable totcatch;
- dvariable R0_use;
- dvariable SSB_use;
dvar_matrix catage_w(1, gmorph, 0, nages);
dvar_vector tempcatch(1, Nfleet);
imatrix Do_F_tune(t_base, TimeMax_Fcast_std, 1, Nfleet); // flag for doing F from catch
@@ -2135,23 +2246,41 @@ FUNCTION void Get_Forecast()
}
}
- if (H4010_top_rd < 0.0)
- {
- H4010_top = Bmsy / SSB_unf;
- if (H4010_bot > 0.25)
- {
- warnstream << "control rule cutoff is large (" << H4010_bot << "); so may not be < calculated Bmsy/SSB_unf (" << H4010_top << ")";
- write_message (WARN, 0);
- }
- }
- else
+ if (Fcast_Loop_Control(5) == 1)
+ {HCR_anchor = SSB_virgin;}
+ else if (Fcast_Loop_Control(5) == 2)
+ {HCR_anchor = SSB_unf;}
+ else if (Fcast_Loop_Control(5) == 3)
+ {HCR_anchor = Bmsy;} // so H4010_top_rd should be 1.0;
+
+ if (H4010_top_rd < 0.0) // legacy approach. This has already been converted to new approach in readdata
+ {
+ H4010_top = Bmsy / HCR_anchor; // convert to fraction of anchor
+ if (H4010_bot > 0.25)
{
- H4010_top = H4010_top_rd;
+ warnstream << "control rule cutoff is large (" << H4010_bot << "); so may not be < calculated Bmsy/SSB_unf (" << H4010_top << ")";
+ write_message (WARN, 0);
}
+ }
+ else
+ {
+ H4010_top = H4010_top_rd;
+ }
+ if (Fcast_Loop_Control(5) == 3 && H4010_top_rd != 1.0)
+ {
+ warnstream << "HCR_anchor is BMSY; so H4010_top normally is 1.0; are you sure you want: " << H4010_top;
+ write_message (WARN, 0);
+ }
+
+ Mgmt_quant(20) = H4010_top * HCR_anchor;
report5 << "#" << endl;
report5 << "N_forecast_yrs: " << N_Fcast_Yrs << endl;
- report5 << "OY_Control_Rule "
- << " Inflection: " << H4010_top << " Intercept: " << H4010_bot << " Scale: " << H4010_scale_vec(endyr + 1) << "; ";
+ report5 << "OY_Control_Rule: Inflection: " << H4010_top << " Intercept: " << H4010_bot << " Scale: " << H4010_scale_vec(endyr + 1) << endl
+ << "Control_rule_anchor_approach: " << Fcast_Loop_Control(5) << " HCR_anchor: " << HCR_anchor << endl
+ << "intercept(SSB): " << H4010_bot * HCR_anchor << endl
+ << "Inflection(SSB): " << H4010_top * HCR_anchor << endl
+ << "#" << endl;
+
switch (HarvestPolicy)
{
case 0: // none
@@ -2181,8 +2310,7 @@ FUNCTION void Get_Forecast()
break;
}
}
- report5 << "#" << endl;
- }
+ }
int jloop;
if (fishery_on_off == 1 || Do_Dyn_Bzero > 0)
@@ -2197,6 +2325,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 +2359,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 +2397,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,26 +2408,26 @@ 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;
for (f = 1; f <= N_SRparm2; f++)
{
- if (SR_parm_timevary(f) == 0)
+ if (SRparm_timevary(f) == 0)
{
- // no change to SR_parm_work
+ // no change to SRparm_work
}
else
{
- SR_parm_work(f) = parm_timevary(SR_parm_timevary(f), y);
+ SRparm_work(f) = parm_timevary(SRparm_timevary(f), y);
}
- SR_parm_byyr(y, f) = SR_parm_work(f);
+ SRparm_byyr(y, f) = SRparm_work(f);
}
env_data(y, -1) = log(SSB_current / SSB_yr(styr - 1)); // store most recent value for density-dependent effects, NOTE - off by a year if recalc'ed at beginning of season 1
env_data(y, -2) = recdev(y); // store for density-dependent effects
@@ -2473,6 +2603,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)
@@ -2555,48 +2686,63 @@ 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
{
natage(t, p, g, 0) = 0.0; // these fish do not yet exist
- 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
- {
- R0_use = Recr_virgin;
- SSB_use = SSB_virgin;
- }
- else
- {
- R0_use = mfexp(SR_parm_work(1));
- equ_Recr = R0_use;
- Fishon = 0;
- eq_yr = y;
- bio_yr = y;
- Do_Equil_Calc(R0_use); // call function to do equilibrium calculation
- if (fishery_on_off == 1)
- {
- Fishon = 1;
- }
- else
- {
- Fishon = 0;
- }
- SSB_use = SSB_equil;
- }
+ // SPAWN-RECR: get recruitment in at beginning of a season in forecast;
+ if (timevary_SRparm(y) == 0) // SRparm use virgin values (but regime still could be)
+ {
+ R0_use = Recr_virgin;
+ SSB_use = SSB_virgin;
+// warning << y << " virgin_SRR; SSB_use: "< 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
+ // SS_Label_Info_24.3.4.1 #Get recruitment from this spawning biomass after start of the season
+ // SPAWN-RECR
+ if (timevary_SRparm(y) == 0) // SRparm use virgin values (but regime still could be)
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
}
- else
+ else if (timevary_SRparm(y) == 1) // update R0_use and SSB_use in this year
+ // values will carry forward into subsequent years
{
- R0_use = mfexp(SR_parm_work(1));
+ R0_use = mfexp(SRparm_work(1));
+ // timevary steepness is in SRparm_work(2) and will be applied inside of Equil_Spawn_Recr_Fxn() and Spawn_Recr()
equ_Recr = R0_use;
Fishon = 0;
eq_yr = y;
bio_yr = y;
- Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
+ SSBpR_Calc(R0_use); // call function to do per recruit calculation with current year's biology and adjusted R0
+ SSB_use = SSB_equil;
+// warning << y << " update_SRR; SSB_use: "< 0)
@@ -3486,7 +3639,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);
}
@@ -3643,16 +3796,24 @@ 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
+ SSBpR_Calc(equ_Recr); // call function to do per recruit 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(SRparm_work(3));
+ beta = mfexp(SRparm_work(4));
+ SRparm_byyr(y, 2) = alpha * temp / (4. + alpha * temp); // implied steepness
+ SRparm_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
+ SSBpR_Calc(equ_Recr); // call function to do per recruit calculation
if (STD_Yr_Reverse_Ofish(y) > 0)
- SPR_std(STD_Yr_Reverse_Ofish(y)) = SSB_equil / Smry_Table(y, 11);
+ {SPR_std(STD_Yr_Reverse_Ofish(y)) = SSB_equil / Smry_Table(y, 11);}
Smry_Table(y, 9) = totbio;
Smry_Table(y, 10) = smrybio;
Smry_Table(y, 12) = SSB_equil;
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_objfunc.tpl b/SS_objfunc.tpl
index 82b91498..7d6315e3 100644
--- a/SS_objfunc.tpl
+++ b/SS_objfunc.tpl
@@ -745,8 +745,8 @@ FUNCTION void evaluate_the_objective_function()
//The recruitment prior is assumed to be a lognormal pdf with expected
// value equal to the deterministic stock-recruitment curve // SS_Label_260
// R1 deviation is weighted by ave_age because R1 represents a time series of recruitments
- // SR_parm(N_SRparm+1) is sigmaR
- // SR_parm(N_SRparm+4) is rho, the autocorrelation coefficient
+ // SRparm(N_SRparm+1) is sigmaR
+ // SRparm(N_SRparm+4) is rho, the autocorrelation coefficient
// POP code from Ianelli
// if (rho>0)
// for (i=styr_rec+1;i<=endyr;i++)
@@ -767,7 +767,7 @@ FUNCTION void evaluate_the_objective_function()
}
else
{
- rho = SR_parm(N_SRparm2);
+ rho = SRparm(N_SRparm2);
recr_like += square(recdev(recdev_first)) / two_sigmaRsq;
for (y = recdev_first + 1; y <= recdev_end; y++)
{
@@ -778,7 +778,7 @@ FUNCTION void evaluate_the_objective_function()
}
else
{
- rho = SR_parm(N_SRparm2);
+ rho = SRparm(N_SRparm2);
dvariable dev;
dvariable dev_last;
if (recdev_first >= styr)
@@ -884,10 +884,10 @@ FUNCTION void evaluate_the_objective_function()
}
for (i = 1; i <= N_SRparm3; i++)
- if (SR_parm_PRtype(i) > 0 && (active(SR_parm(i)) || Do_all_priors > 0))
+ if (SRparm_PRtype(i) > 0 && (active(SRparm(i)) || Do_all_priors > 0))
{
- SR_parm_Like(i) = Get_Prior(SR_parm_PRtype(i), SR_parm_LO(i), SR_parm_HI(i), SR_parm_PR(i), SR_parm_CV(i), SR_parm(i));
- parm_like += SR_parm_Like(i);
+ SRparm_Like(i) = Get_Prior(SRparm_PRtype(i), SRparm_LO(i), SRparm_HI(i), SRparm_PR(i), SRparm_CV(i), SRparm(i));
+ parm_like += SRparm_Like(i);
}
// SS_Label_Info_25.14 #logL for recdev_cycle
if (recdev_cycle > 0)
@@ -1150,7 +1150,7 @@ FUNCTION void Process_STDquant()
switch (SPR_reporting)
{
- case 0: // keep as raw value
+ case 0: // no scaling. this option skips SPR_reporting
{
break;
}
@@ -1170,11 +1170,15 @@ FUNCTION void Process_STDquant()
SPR_std = (1. - SPR_std) / (1. - SPR_Btgt);
break;
}
- case 4:
+ case 4: // 1-%SPR
{
SPR_std = 1. - SPR_std;
break;
}
+ case 5: // raw %SPR
+ {
+ break;
+ }
}
switch (depletion_basis)
@@ -1187,21 +1191,25 @@ FUNCTION void Process_STDquant()
case 1:
{
depletion /= (depletion_level * SSB_virgin);
+ Mgmt_quant(19) = SSB_virgin;
break;
}
case 2:
{
depletion /= (depletion_level * Bmsy);
+ Mgmt_quant(19) = Bmsy;
break;
}
case 3:
{
depletion /= (depletion_level * SSB_yr(styr));
+ Mgmt_quant(19) = SSB_yr(styr);
break;
}
case 4:
{
depletion /= (depletion_level * SSB_yr(endyr));
+ Mgmt_quant(19) = SSB_yr(endyr);
break;
}
case 5: // dynamic Bzero
@@ -1213,6 +1221,13 @@ FUNCTION void Process_STDquant()
depletion(STD_Yr_Reverse_Dep(y)) /= ( depletion_level * Extra_Std(Do_Dyn_Bzero + y - (styr - 2)));
// warning< 0)
diff --git a/SS_param.tpl b/SS_param.tpl
index 9ddf9a86..d8ff05f9 100644
--- a/SS_param.tpl
+++ b/SS_param.tpl
@@ -157,22 +157,26 @@ PARAMETER_SECTION
4darray recr_dist(styr-3,YrMax,1,N_GP*gender,1,N_settle_timings,1,pop);
3darray recr_dist_unf(1,N_GP*gender,1,N_settle_timings,1,pop);
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
- vector SR_parm_virg(1,N_SRparm2+1)
- vector SR_parm_work(1,N_SRparm2+1)
+!!// SS_Label_Info_5.1.2 #Create SRparm vector, recruitment vectors
+ init_bounded_number_vector SRparm(1,N_SRparm3,SRparm_LO,SRparm_HI,SRparm_PH)
+ matrix SRparm_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 SRparm_virg(1,N_SRparm2+1)
+ vector SRparm_work(1,N_SRparm2+1)
+ vector SRparm_bench(1,N_SRparm2+1)
number two_sigmaRsq;
number half_sigmaRsq;
number sigmaR;
- number SPR_virgin;
+ number SSBpR_virgin;
+ number SSBpR_bench;
+ number SSB0_4_SRR;
+ number R0_4_SRR;
number regime_change;
number rho;
number dirichlet_Parm;
LOCAL_CALCS
// clang-format on
Ave_Size.initialize();
- // if(SR_parm(N_SRparm2)!=0.0 || SR_parm_PH(N_SRparm2)>0) {SR_autocorr=1;} else {SR_autocorr=0;} // flag for recruitment autocorrelation
+ // if(SRparm(N_SRparm2)!=0.0 || SRparm_PH(N_SRparm2)>0) {SR_autocorr=1;} else {SR_autocorr=0;} // flag for recruitment autocorrelation
if (do_recdev == 1)
{
k = recdev_start;
@@ -234,21 +238,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
+ 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 +329,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 +424,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;
@@ -625,7 +634,7 @@ PARAMETER_SECTION
vector init_F_Like(1,N_init_F)
vector Q_parm_Like(1,Q_Npar2)
vector selparm_Like(1,N_selparm2)
- vector SR_parm_Like(1,N_SRparm3)
+ vector SRparm_Like(1,N_SRparm3)
vector recdev_cycle_Like(1,recdev_cycle)
!! k=Do_TG*(3*N_TG+2*Nfleet1);
vector TG_parm_Like(1,k);
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index e6c64a53..90e8cf2e 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -1,14 +1,14 @@
// SS_Label_file #12. **SS_popdyn.tpl**
// SS_Label_file # * setup_recdevs()
-// SS_Label_file # * get_initial_conditions() // does virgin and initial year by calling Do_Equil_Calc() with F=0, then F=init_F
+// SS_Label_file # * get_initial_conditions() // does virgin and initial year by calling SSBpR_Calc() with F=0, then F=init_F
// SS_Label_file # * get_time_series() // loops the years, calling biology, selectivity and spawn-recr functions as needed
-// SS_Label_file # * Do_Equil_Calc() // does per-recruit calculations and returns SSB/R and Y/R
+// SS_Label_file # * SSBpR_Calc() // does per-recruit calculations and returns SSB/R and Y/R
// SS_Label_file #
FUNCTION void setup_recdevs()
{
// SS_Label_Info_7.1 #Set up recruitment bias_adjustment vector
- sigmaR = SR_parm(N_SRparm(SR_fxn) + 1);
+ sigmaR = SRparm(N_SRparm(SR_fxn) + 1);
two_sigmaRsq = 2.0 * sigmaR * sigmaR;
half_sigmaRsq = 0.5 * sigmaR * sigmaR;
@@ -269,6 +269,7 @@ FUNCTION void get_initial_conditions()
}
Wt_Age_t(t, 0) = Wt_Age_beg(s);
+ Wt_Age_t(t, -1) = Wt_Age_mid(s);
for (g = 1; g <= gmorph; g++)
if (use_morph(g) > 0)
@@ -337,44 +338,67 @@ FUNCTION void get_initial_conditions()
Fishon = 0;
virg_fec = fec;
Recr.initialize(); // will store recruitment by area
- 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));
+ equ_Recr = 1.0;
+ SSBpR_Calc(equ_Recr); // call function to do per recruit calculation. Returns SPR because R = 1.0
+ SSBpR_virgin = SSB_equil; // spawners per recruit. Needed for Sr_fxn = 10
+ 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(SRparm(3));
+ beta = mfexp(SRparm(4));
+ steepness = alpha * SSBpR_virgin / (4. + alpha * SSBpR_virgin);
+ Recr_virgin = 1. / beta * (alpha - (1. / SSBpR_virgin));
+ SRparm(1) = log(Recr_virgin);
+ SRparm(2) = steepness;
+ }
+ else
+ {
+ Recr_virgin = mfexp(SRparm(1));
+ }
+
+ for (int i = 1; i <= N_SRparm2; i++)
+ {
+ SRparm_byyr(eq_yr, i) = SRparm(i);
+ SRparm_virg(i) = SRparm(i);
+ SRparm_work(i) = SRparm(i);
+ }
+// if (SR_fxn == 3) warning << "tester_A: " << SRparm_work(1) << " base: " << SRparm(1) << endl;
+// if (SR_fxn == 10) warning << "tester_A: " << SRparm_work(4) << " base: " << SRparm(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;
exp_rec(eq_yr, 3) = Recr_virgin;
exp_rec(eq_yr, 4) = Recr_virgin;
- Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
+ SSBpR_Calc(equ_Recr); // call function to do per recruit calculation
SSB_virgin = SSB_equil;
- SPR_virgin = SSB_equil / Recr_virgin; // spawners per recruit
- if(Do_Benchmark==0)
+ 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;
+ SSB_unf = SSB_virgin;
+ Mgmt_quant(1) = SSB_unf; // will be overwritten in benchmark
+ Recr_unf = Recr_virgin; // will be overwritten in benchmark
+ Mgmt_quant(2) = totbio; // from Do_Equil_Calc
+ Mgmt_quant(3) = smrybio; // from Do_Equil_Calc
+ 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;
- SR_parm_work(N_SRparm2 + 1) = SSB_equil;
+ SRparm_byyr(eq_yr, N_SRparm2 + 1) = SSB_equil;
+ SRparm_virg(N_SRparm2 + 1) = SSB_equil;
+ SRparm_work(N_SRparm2 + 1) = SSB_equil;
t = styr - 2 * nseas - 1;
for (s = 1; s <= nseas; s++)
for (p = 1; p <= pop; p++)
@@ -447,15 +471,18 @@ FUNCTION void get_initial_conditions()
for (f = 1; f <= N_SRparm2; f++)
{
- if (SR_parm_timevary(f) == 0)
+ if (SRparm_timevary(f) == 0)
{
- // no change to SR_parm_work
+ // no change to SRparm_work
}
else
{
- SR_parm_work(f) = parm_timevary(SR_parm_timevary(f), eq_yr);
+ SRparm_work(f) = parm_timevary(SRparm_timevary(f), eq_yr);
+// warning << "tester_B: " << SRparm_work(f) << " timevary " << " base " << SRparm(f) < 0) // timevary regime exists
+ if (SRparm_timevary(N_SRparm2 - 1) > 0) // timevary regime exists
{
- regime_change = mfexp(SR_parm_work(N_SRparm2 - 1));
+ regime_change = mfexp(SRparm_work(N_SRparm2 - 1));
}
if (init_equ_steepness == 0) // Adjustments do not include spawner-recruitment steepness
@@ -490,8 +517,8 @@ FUNCTION void get_initial_conditions()
exp_rec(eq_yr, 2) = R1;
exp_rec(eq_yr, 3) = R1;
exp_rec(eq_yr, 4) = R1;
- equ_Recr = R1; // equ_Recr is used inside of Do_Equil_Calc
- Do_Equil_Calc(equ_Recr);
+ equ_Recr = R1; // equ_Recr is used inside of SSBpR_Calc
+ SSBpR_Calc(equ_Recr);
CrashPen += Equ_penalty;
}
else
@@ -499,16 +526,16 @@ FUNCTION void get_initial_conditions()
// SS_Label_Info_23.5.1.2 #Adjustments include spawner-recruitment function
// do initial equilibrium with R1 based on offset from spawner-recruitment curve, using same approach as the benchmark calculations
// first get SPR for this init_F
- // SPAWN-RECR: calc initial equilibrium pop, SPB, Recruitment
+ // SPAWN-RECR: calc initial equilibrium pop, SSB, Recruitment
// equ_Recr=Recr_virgin;
- equ_Recr = R1_exp * regime_change;
-
- Do_Equil_Calc(equ_Recr);
+// equ_Recr = R1_exp * regime_change; // NOTE: seems wrong to apply regime here
+ equ_Recr = R1_exp;
+ SSBpR_Calc(equ_Recr);
CrashPen += Equ_penalty;
- SPR_temp = SSB_equil / equ_Recr; // spawners per recruit at initial F
- // get equilibrium SSB and recruitment from SPR_temp, Recr_virgin and virgin steepness
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm(2), SR_parm(3), SSB_virgin, Recr_virgin, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
-
+ SSBpR_temp = SSB_equil / equ_Recr; // spawners per recruit at initial F
+ // get equilibrium SSB and recruitment from SSBpR_temp, Recr_virgin and virgin steepness
+ // this is the initial year, so no time-vary effects available, so uses _virgin quantities for spawner-recruitment
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SRparm_work, SSB_virgin, Recr_virgin, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
R1_exp = Equ_SpawnRecr_Result(2); // set the expected recruitment equal to this equilibrium
exp_rec(eq_yr, 1) = R1_exp;
@@ -517,7 +544,7 @@ FUNCTION void get_initial_conditions()
exp_rec(eq_yr, 3) = equ_Recr;
exp_rec(eq_yr, 4) = equ_Recr;
R1 = equ_Recr;
- Do_Equil_Calc(equ_Recr); // calculated SSB_equil
+ SSBpR_Calc(equ_Recr); // calculated SSB_equil
CrashPen += Equ_penalty;
}
Smry_Table(styr - 1, 1) = totbio; // from equil calcs
@@ -526,10 +553,10 @@ FUNCTION void get_initial_conditions()
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_work(N_SRparm2 + 1) = SSB_equil;
+ SRparm_byyr(eq_yr, N_SRparm2 + 1) = SSB_equil;
+ SRparm_work(N_SRparm2 + 1) = SSB_equil;
SSB_yr(styr) = SSB_equil;
env_data(styr - 1, -1) = 0.0;
env_data(styr - 1, -2) = 0.0;
@@ -669,14 +696,12 @@ FUNCTION void get_time_series()
dvariable crashtemp1;
dvariable interim_tot_catch;
dvariable Z_adjuster;
- dvariable R0_use;
- dvariable SSB_use;
if (Do_Morphcomp > 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)
{
@@ -697,15 +722,16 @@ FUNCTION void get_time_series()
for (f = 1; f <= N_SRparm2; f++)
{
- if (SR_parm_timevary(f) == 0)
+ if (SRparm_timevary(f) == 0)
{
- // no change to SR_parm_work
+ // no change to SRparm_work
}
else
{
- SR_parm_work(f) = parm_timevary(SR_parm_timevary(f), y);
+ SRparm_work(f) = parm_timevary(SRparm_timevary(f), y);
+// warning << "tester_C: " << SRparm_work(f) << " timevary_year " << endl;
}
- SR_parm_byyr(y, f) = SR_parm_work(f);
+ SRparm_byyr(y, f) = SRparm_work(f);
}
// SS_Label_Info_24.1.1 #store begin of year quantities for use in density-dependent processes
@@ -900,6 +926,7 @@ FUNCTION void get_time_series()
}
Wt_Age_t(t, 0) = Wt_Age_beg(s);
+ Wt_Age_t(t, -1) = Wt_Age_mid(s);
if (y > styr) // because styr is done as part of initial conditions
{
@@ -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,37 +1007,43 @@ 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
- // SPAWN-RECR: calc recruitment in time series; need to make this area-specific
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ // 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
+ if (timevary_SRparm(y) == 0) // SRparm use virgin values (but regime still could be)
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
+// warning << y << " virgin_SRR; SSB_use: "< 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 (timevary_SRparm(y) == 0) // SRparm use virgin values (but regime still could be)
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
+// warning << y << " virgin_SRR; SSB_use: "< 0)
{
SPR_std(STD_Yr_Reverse_Ofish(y)) = SSB_equil / Smry_Table(y, 11);
@@ -1796,7 +1850,8 @@ FUNCTION void get_time_series()
//********************************************************************
/* SS_Label_FUNCTION 30 Do_Equil_Calc */
-FUNCTION void Do_Equil_Calc(const prevariable& equ_Recr)
+ // This function does per recruit calculations, so produces an age composition that is in equilibrium with M+F
+FUNCTION void SSBpR_Calc(const prevariable& equ_Recr)
{
int t_base;
int t;
@@ -2244,9 +2299,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_prelim.tpl b/SS_prelim.tpl
index b23dff21..c37b8fc0 100644
--- a/SS_prelim.tpl
+++ b/SS_prelim.tpl
@@ -689,9 +689,9 @@
for (i = 1; i <= N_SRparm3; i++)
{
- SR_parm(i) = SR_parm_RD(i);
+ SRparm(i) = SRparm_RD(i);
}
- echoinput << " SRR_parms read from ctl " << SR_parm << endl;
+ echoinput << " SRR_parms read from ctl " << SRparm << endl;
if (recdev_cycle > 0)
{
@@ -701,12 +701,9 @@
}
}
- if (recdev_do_early > 0)
- recdev_early.initialize();
- if (Do_Forecast > 0 && do_recdev != 0)
- Fcast_recruitments.initialize();
- if (Do_Impl_Error > 0)
- Fcast_impl_error.initialize();
+ if (recdev_do_early > 0) recdev_early.initialize();
+ if (Do_Forecast > 0 && do_recdev != 0) Fcast_recruitments.initialize();
+ if (Do_Impl_Error > 0) Fcast_impl_error.initialize();
if (do_recdev == 1)
{
@@ -863,13 +860,13 @@
MGparm_use = value(MGparm);
echoinput << endl
- << " now check SR_parm bounds and priors and do jitter if requested " << endl;
+ << " now check SRparm bounds and priors and do jitter if requested " << endl;
for (i = 1; i <= N_SRparm3; i++)
{
- SR_parm(i) = Check_Parm(i, SR_parm_PH(i), SR_parm_LO(i), SR_parm_HI(i), SR_parm_PRtype(i), SR_parm_PR(i), SR_parm_CV(i), jitter, SR_parm(i));
+ SRparm(i) = Check_Parm(i, SRparm_PH(i), SRparm_LO(i), SRparm_HI(i), SRparm_PRtype(i), SRparm_PR(i), SRparm_CV(i), jitter, SRparm(i));
}
- echoinput << " SRR_parms after check " << SR_parm << endl;
- SR_parm_use = value(SR_parm);
+ echoinput << " SRR_parms after check " << SRparm << endl;
+ SRparm_use = value(SRparm);
recdev_use.initialize();
if (recdev_cycle > 0)
diff --git a/SS_proced.tpl b/SS_proced.tpl
index 534adf53..8d4896aa 100644
--- a/SS_proced.tpl
+++ b/SS_proced.tpl
@@ -24,8 +24,8 @@ PROCEDURE_SECTION
{
if (mcmc_counter == 0)
{
- SR_parm(1) += MCMC_bump;
- cout << mcmc_counter << " adjusted SR_parm in first mcmc call " << SR_parm(1) << " by " << MCMC_bump << endl;
+ SRparm(1) += MCMC_bump;
+ cout << mcmc_counter << " adjusted SRparm in first mcmc call " << SRparm(1) << " by " << MCMC_bump << endl;
}
mcmc_counter++;
@@ -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)
{
@@ -255,11 +255,11 @@ PROCEDURE_SECTION
ParmTrace << " " << MGparm(j);
}
}
- for (j = 1; j <= SR_parm_PH.indexmax(); j++)
+ for (j = 1; j <= SRparm_PH.indexmax(); j++)
{
- if (SR_parm_PH(j) >= 0)
+ if (SRparm_PH(j) >= 0)
{
- ParmTrace << " " << SR_parm(j);
+ ParmTrace << " " << SRparm(j);
}
}
if (recdev_cycle > 0)
@@ -384,7 +384,7 @@ PROCEDURE_SECTION
ParmTrace << current_phase() << " " << niter << " " << obj_fun << " " << obj_fun - last_objfun
<< " " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr)) << " " << biasadj(styr) << " " << max(biasadj) << " " << biasadj(endyr);
ParmTrace << " " << MGparm << " ";
- ParmTrace << SR_parm << " ";
+ ParmTrace << SRparm << " ";
if (recdev_cycle > 0)
ParmTrace << recdev_cycle_parm;
if (recdev_do_early > 0)
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index bc42d9b5..88992625 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -52,6 +52,7 @@
init_int WTage_rd // 0 means do not read wtatage.ss; 1 means read and use wtatage.ss and also read and use growth parameters
// future option 2 will suppress reading and use of growth
!!echoinput< 0) timevary_MG_firstyr = styr;
init_int N_GP // number of growth patterns (morphs)
!!echoinput<> 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
@@ -1663,6 +1664,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;
@@ -1711,6 +1713,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();
@@ -1811,14 +1814,31 @@
{
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
if ((timevary_MG(y, 2) > 0 || timevary_MG(y, 3) > 0) && natM_type == 5 && natM_5_opt < 3) timevary_MG(y, 1) = 1;
echoinput << y << " timevary_MG: " << timevary_MG(y) << endl;
}
+ if( timevary_MG_firstyr < YrMax || WTage_rd == 1) // time-varying biology
+ {
+ if( timevary_bio_4SRR == 0) // legacy approach; this switch is read from starter.ss
+ {
+ warnstream << "There is timevary biology and the legacy approach to benchmark calculations is being used; user should be aware of possible impacts to benchmark results";
+ write_message(WARN, 0);
+ if( timevary_bio_4SRR_rd == -1) // older starter file did not contain necessary flag
+ {
+ warnstream << "There is timevary biology, so the flag for timevary_bio_4SRR must be set to 0 (old default) or 1 (new improved) in starter.ss";
+// make this a WARN while testing, then change to FATAL for operational code
+ write_message(WARN, 0);
+ }
+ }
+ }
+
for (y = endyr + 1; y <= YrMax; y++)
{
for (f = 1; f <= 7; f++)
@@ -1944,18 +1964,20 @@
!!// 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) varparm_estimated(2) = 1; // sigmaR is estimated so need sd_offset=1
+ if (SRparm_1(N_SRparm2 - 2, 7) > 0) varparm_estimated(2) = 1; // sigmaR is estimated so need sd_offset=1
- if (SR_parm_1(N_SRparm2, 3) != 0.0 || SR_parm_1(N_SRparm2, 7) > 0)
+ if (SRparm_1(N_SRparm2, 3) != 0.0 || SRparm_1(N_SRparm2, 7) > 0)
{
SR_autocorr = 1;
}
@@ -2051,23 +2087,25 @@
}
// flag for recruitment autocorrelation
echoinput << " Do recruitment_autocorr: " << SR_autocorr << endl;
+
+ // note that the regime parameter seems to bypass use of timevary_SRparm, but timevary_SRparm is used for R0, h beginning 3.30.24
timevary_used = 0;
- for (j = 1; j <= N_SRparm(SR_fxn) + 2; j++)
- if (j != N_SRparm(SR_fxn) + 1) // because sigmaR and autocorr cannot be time-varying
+ for (j = 1; j <= N_SRparm2 - 1; j++) // so omits autocorr
+ if (j != N_SRparm2 - 2) // because sigmaR cannot be time-varying
{
- if (SR_parm_1(j, 13) == 0 && SR_parm_1(j, 8) == 0 && SR_parm_1(j, 9) == 0)
+ if (SRparm_1(j, 13) == 0 && SRparm_1(j, 8) == 0 && SRparm_1(j, 9) == 0)
{
// no time-vary parameter effects
}
else // set up a timevary parameter definition
{
+ timevary_used = 1;
ivector timevary_setup(1, 14); // temporary vector for timevary specs
timevary_setup.initialize();
- if (timevary_parm_start_SR == 0) timevary_parm_start_SR = timevary_parm_cnt + 1;
+ if (timevary_SRparm_first == 0) timevary_SRparm_first = timevary_parm_cnt + 1; // cumulative index for first timevary SRparm
echoinput << " timevary for SR parm: " << j << endl;
- timevary_used = 1;
timevary_cnt++; // count parameters with time-vary effect
- SR_parm_timevary(j) = timevary_cnt; // base SR parameter will use this timevary specification
+ SRparm_timevary(j) = timevary_cnt; // base SR parameter will use this timevary specification
timevary_setup(1) = 2; // indicates a SR parm
if (autogen_timevary(2) == 0)
{
@@ -2080,96 +2118,109 @@
timevary_setup(2) = j; // index of base parm within that type of parameter
timevary_setup(13) = firstSRparm + j; // index of base parm relative to ParCount which is continuous across all types of parameters
timevary_setup(3) = timevary_parm_cnt + 1; // first parameter within total list of all timevary parms
- timevary_pass = 1; // placeholder; not used for SR parms
+ timevary_pass = 0; // placeholder; not used for SR parms
// set up env link info
- echoinput << " check for env " << SR_parm_1(j, 8) << endl;
- k = int(abs(SR_parm_1(j, 8)) / 100); // find the env link code
+ echoinput << " check for env " << SRparm_1(j, 8) << endl;
+ k = int(abs(SRparm_1(j, 8)) / 100); // find the env link code
timevary_setup(6) = k; // link code for env
- if (SR_parm_1(j, 8) > 0) // env variable used
+ if (SRparm_1(j, 8) > 0) // env variable used
{
- timevary_setup(7) = int(abs(SR_parm_1(j, 8))) - k * 100;
+ timevary_setup(7) = int(abs(SRparm_1(j, 8))) - k * 100;
k = timevary_setup(7);
// for(y=styr-1;y<=YrMax;y++) env_data_pass(y)=env_data_RD(y,k);
env_data_pass(1) = env_data_minyr(k);
env_data_pass(2) = env_data_maxyr(k);
}
- else if (abs(SR_parm_1(j, 8) > 0)) // density-dependence
+ else if (abs(SRparm_1(j, 8) > 0)) // density-dependence
{
- timevary_setup(7) = -int(abs(SR_parm_1(j, 8)) - k * 100);
+ timevary_setup(7) = -int(abs(SRparm_1(j, 8)) - k * 100);
do_densitydependent = 1;
k = 0;
env_data_pass.initialize();
}
- if (SR_parm_1(j, 13) > 0) // doing blocks
+ if (SRparm_1(j, 13) > 0) // doing blocks
{
- if (SR_parm_1(j, 13) > N_Block_Designs)
+ if (SRparm_1(j, 13) > N_Block_Designs)
{
warnstream << "SR block request exceeds N_block patterns";
write_message (FATAL, 0); // EXIT!
}
- create_timevary(SR_parm_1(j), timevary_setup, timevary_pass, autogen_timevary(timevary_setup(1)), f, Block_Design(SR_parm_1(j, 13)), env_data_pass, N_parm_dev, finish_starter);
+ create_timevary(SRparm_1(j), timevary_setup, timevary_pass, autogen_timevary(timevary_setup(1)), f, Block_Design(SRparm_1(j, 13)), env_data_pass, N_parm_dev, finish_starter);
}
else
{
- create_timevary(SR_parm_1(j), timevary_setup, timevary_pass, autogen_timevary(timevary_setup(1)), f, block_design_null, env_data_pass, N_parm_dev, finish_starter);
+ create_timevary(SRparm_1(j), timevary_setup, timevary_pass, autogen_timevary(timevary_setup(1)), f, block_design_null, env_data_pass, N_parm_dev, finish_starter);
}
timevary_def.push_back(timevary_setup(1, 14));
- for (y = styr - 3; y <= YrMax + 1; y++) {
- timevary_SRparm(y) = timevary_pass(y);
- } // year vector for this category og MGparm
+ int SRflag;
+ SRflag = 0;
+ for (y = styr - 3; y <= YrMax + 1; y++)
+ {
+ if (timevary_pass(y) > 0 && j != N_SRparm2 - 1)
+ {
+ timevary_SRparm(y) = timevary_pass(y); // set timevary flag, except for regime parameter
+ timevary_SRparm_first_yr = y;
+ SRflag = 1; // first change point
+ }
+ else if(SRflag == 1)
+ {
+ timevary_SRparm(y) = 2; // flag to carry forward current SRR info
+ }
+ }
}
}
+
N_SRparm3 = N_SRparm2;
- if (timevary_parm_start_SR > 0)
+ if (timevary_SRparm_first > 0)
{
- timevary_parm_cnt_SR = timevary_parm_cnt;
+ timevary_parm_SR_last = timevary_parm_cnt;
if (timevary_used == 1) autogen_timevary(2) = 1; // indicate that some parameter is time-varying
- N_SRparm3 += (timevary_parm_cnt_SR - timevary_parm_start_SR + 1);
- 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;
+ N_SRparm3 += (timevary_parm_SR_last - timevary_SRparm_first + 1);
+ echoinput << " SR timevary_parm_cnt start and end " << timevary_SRparm_first << " " << timevary_parm_SR_last << endl;
+ echoinput << "link to timevary parms: " << SRparm_timevary << endl;
}
echoinput << "SR_Npar and N_SRparm2 and N_SRparm3: " << N_SRparm(SR_fxn) << " " << N_SRparm2 << " " << N_SRparm3 << endl;
// clang-format off
END_CALCS
- vector SR_parm_LO(1,N_SRparm3)
- vector SR_parm_HI(1,N_SRparm3)
- vector SR_parm_RD(1,N_SRparm3)
- vector SR_parm_PR(1,N_SRparm3)
- ivector SR_parm_PRtype(1,N_SRparm3)
- vector SR_parm_CV(1,N_SRparm3)
- ivector SR_parm_PH(1,N_SRparm3)
+ vector SRparm_LO(1,N_SRparm3)
+ vector SRparm_HI(1,N_SRparm3)
+ vector SRparm_RD(1,N_SRparm3)
+ vector SRparm_PR(1,N_SRparm3)
+ ivector SRparm_PRtype(1,N_SRparm3)
+ vector SRparm_CV(1,N_SRparm3)
+ ivector SRparm_PH(1,N_SRparm3)
LOCAL_CALCS
// clang-format on
for (i = 1; i <= N_SRparm2; i++)
{
- SR_parm_LO(i) = SR_parm_1(i, 1);
- SR_parm_HI(i) = SR_parm_1(i, 2);
- SR_parm_RD(i) = SR_parm_1(i, 3);
- SR_parm_PR(i) = SR_parm_1(i, 4);
- SR_parm_CV(i) = SR_parm_1(i, 5);
- SR_parm_PRtype(i) = SR_parm_1(i, 6);
- SR_parm_PH(i) = SR_parm_1(i, 7);
+ SRparm_LO(i) = SRparm_1(i, 1);
+ SRparm_HI(i) = SRparm_1(i, 2);
+ SRparm_RD(i) = SRparm_1(i, 3);
+ SRparm_PR(i) = SRparm_1(i, 4);
+ SRparm_CV(i) = SRparm_1(i, 5);
+ SRparm_PRtype(i) = SRparm_1(i, 6);
+ SRparm_PH(i) = SRparm_1(i, 7);
}
- if (timevary_parm_start_SR > 0)
+ if (timevary_SRparm_first > 0)
{
j = N_SRparm2;
- for (f = timevary_parm_start_SR; f <= timevary_parm_cnt_SR; f++)
+ for (f = timevary_SRparm_first; f <= timevary_parm_SR_last; f++)
{
j++;
echoinput << f << " " << j << " " << timevary_parm_rd[f] << endl;
- SR_parm_LO(j) = timevary_parm_rd[f](1);
- SR_parm_HI(j) = timevary_parm_rd[f](2);
- SR_parm_RD(j) = timevary_parm_rd[f](3);
- SR_parm_PR(j) = timevary_parm_rd[f](4);
- SR_parm_PRtype(j) = timevary_parm_rd[f](6);
- SR_parm_CV(j) = timevary_parm_rd[f](5);
- SR_parm_PH(j) = timevary_parm_rd[f](7);
+ SRparm_LO(j) = timevary_parm_rd[f](1);
+ SRparm_HI(j) = timevary_parm_rd[f](2);
+ SRparm_RD(j) = timevary_parm_rd[f](3);
+ SRparm_PR(j) = timevary_parm_rd[f](4);
+ SRparm_PRtype(j) = timevary_parm_rd[f](6);
+ SRparm_CV(j) = timevary_parm_rd[f](5);
+ SRparm_PH(j) = timevary_parm_rd[f](7);
}
}
- echoinput << "SR_parm_RD: " << SR_parm_RD << endl;
+ echoinput << "SRparm_RD: " << SRparm_RD << endl;
// clang-format off
END_CALCS
@@ -5766,7 +5817,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;
@@ -5809,7 +5860,7 @@
// clang-format on
if (Do_Benchmark > 0)
{
- N_STD_Mgmt_Quant = 17;
+ N_STD_Mgmt_Quant = 22;
if (Do_Benchmark == 3) N_STD_Mgmt_Quant += 3; // for Blimit
}
else
@@ -5881,19 +5932,19 @@
}
}
- for (j = 1; j <= SR_parm_PH.indexmax(); j++)
+ for (j = 1; j <= SRparm_PH.indexmax(); j++)
{
ParCount++;
- if (SR_parm_PH(j) == -9999) {
- SR_parm_1(j, 3) = prof_var(prof_var_cnt);
- SR_parm_RD(j, 3) = SR_parm_1(j, 3);
+ if (SRparm_PH(j) == -9999) {
+ SRparm_1(j, 3) = prof_var(prof_var_cnt);
+ SRparm_RD(j, 3) = SRparm_1(j, 3);
prof_var_cnt += 1;
}
- if (depletion_fleet > 0 && depletion_type < 2 && SR_parm_PH(j) > 0) SR_parm_PH(j)++; // add 1 to phase if using depletion fleet
- if (depletion_fleet > 0 && depletion_type < 2 && j == 1) SR_parm_PH(1) = 1; // R0 active in phase 1, unless type==2
- if (SR_parm_PH(j) > Turn_off_phase2) SR_parm_PH(j) = -1;
- if (SR_parm_PH(j) > max_phase) max_phase = SR_parm_PH(j);
- if (SR_parm_PH(j) >= 0)
+ if (depletion_fleet > 0 && depletion_type < 2 && SRparm_PH(j) > 0) SRparm_PH(j)++; // add 1 to phase if using depletion fleet
+ if (depletion_fleet > 0 && depletion_type < 2 && j == 1) SRparm_PH(1) = 1; // R0 active in phase 1, unless type==2
+ if (SRparm_PH(j) > Turn_off_phase2) SRparm_PH(j) = -1;
+ if (SRparm_PH(j) > max_phase) max_phase = SRparm_PH(j);
+ if (SRparm_PH(j) >= 0)
{
active_count++;
active_parm(active_count) = ParCount;
@@ -6616,6 +6667,28 @@
CoVar_Count++;
j++;
active_parm(CoVar_Count) = j;
+// add quantities needed when time-vary life history is used; but report here for all cases; elements 18-21 of mgmt_quant
+ ParmLabel += "Recr_MSY_bmarkbio" + CRLF(1);
+ CoVar_Count++;
+ j++;
+ active_parm(CoVar_Count) = j;
+ ParmLabel += "Depletion_denom" + CRLF(1);
+ CoVar_Count++;
+ j++;
+ active_parm(CoVar_Count) = j;
+ ParmLabel += "HCR_inflect" + CRLF(1);
+ CoVar_Count++;
+ j++;
+ active_parm(CoVar_Count) = j;
+ ParmLabel += "R0_for_SRR_bench" + CRLF(1);
+ CoVar_Count++;
+ j++;
+ active_parm(CoVar_Count) = j;
+ ParmLabel += "SSB_for_SRR_bench" + CRLF(1);
+ CoVar_Count++;
+ j++;
+ active_parm(CoVar_Count) = j;
+
if (Do_Benchmark == 3)
{
ParmLabel += "SSB_Blim" + CRLF(1);
@@ -6820,23 +6893,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)
{
@@ -6915,6 +6988,11 @@
depletion_basis_label += " " + onenum + "%*Dyn_Bzero";
break;
}
+ case 6:
+ {
+ depletion_basis_label += " " + onenum + "%*Bmark_Biomass";
+ break;
+ }
}
if (depletion_log == 1) depletion_basis_label += ";log";
if (depletion_multi > 1)
@@ -6925,7 +7003,12 @@
switch (SPR_reporting)
{
- case 0: // keep as raw value
+ case 0: // skip SPR reporting
+ {
+ SPR_report_label += " raw_SPR";
+ break;
+ }
+ case 5: // keep as raw %SPR value
{
SPR_report_label += " raw_SPR";
break;
@@ -7034,7 +7117,7 @@
// containers for parameter values after jitter
vector MGparm_use(1,N_MGparm2)
- vector SR_parm_use(1,N_SRparm3);
+ vector SRparm_use(1,N_SRparm3);
vector recdev_cycle_use(1,recdev_cycle);
vector recdev_use(recdev_first,YrMax);
vector recdev_RD(recdev_first,YrMax);
diff --git a/SS_readdata_330.tpl b/SS_readdata_330.tpl
index d212bcde..0ba0045d 100644
--- a/SS_readdata_330.tpl
+++ b/SS_readdata_330.tpl
@@ -4281,8 +4281,8 @@
echoinput << "HarvestPolicy=0, so values for top, bottom, buffer will be ignored" << endl;
echoinput << HarvestPolicy << " # echoed HarvestPolicy " << endl;
- *(ad_comm::global_datafile) >> H4010_top_rd; // use -1 to set H4010_top to Bmsy/SSB_unf
- echoinput << H4010_top_rd << " # echoed control rule inflection (-1 to set to Bmsy/SSB_unf)" << endl;
+ *(ad_comm::global_datafile) >> H4010_top_rd; // as fraction of HCR_anchor; use -1 as legacy approach to set H4010_top to Bmsy/SSB_unf
+ echoinput << H4010_top_rd << " # echoed control rule inflection" << endl;
*(ad_comm::global_datafile) >> H4010_bot;
echoinput << H4010_bot << " # echoed control rule cutoff " << endl;
*(ad_comm::global_datafile) >> H4010_scale_rd;
@@ -4345,8 +4345,25 @@
"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;
-
+ if (Fcast_Loop_Control(5) <= 0) // default before 3.30.24
+ {
+ echoinput << "basis for HCR anchor was not set; setting to 2 to match default before 3.30.24" << endl;
+ warnstream << "basis for HCR anchor was not set; setting to 2 to match default before 3.30.24";
+ write_message(ADJUST, 0);
+ Fcast_Loop_Control(5) = 2;
+ }
+ if (H4010_top_rd < 0) // convert old legacy approach to new approach for using Bmsy
+ {
+ Fcast_Loop_Control(5) = 3;
+ H4010_top_rd = 1.0;
+ }
+ echoinput << Fcast_Loop_Control(5) << " #control rule anchor: 1=virgin_SSB; 2=unfished_benchmark_SSB(old_approach); 3=Bmsy" << endl;
+ if (depletion_basis == 1 && Fcast_Loop_Control(5) == 2)
+ {
+ warnstream << "depletion_basis is using virgin but HCR anchor is using SSB_unf from benchmark. Are you sure?";
+ write_message(WARN, 0);
+ }
+
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;
@@ -4814,16 +4831,16 @@
warnstream << "Set F_std_basis=0 because no benchmark or forecast";
write_message(WARN, 0);
}
- if (depletion_basis == 2)
+ if (depletion_basis == 2 || depletion_basis == 6 )
{
depletion_basis = 1;
- warnstream << "Change depletion basis to 1 because benchmarks are off";
+ warnstream << "Change depletion basis to 1 because benchmarks were not requested";
write_message(WARN, 0);
}
if (SPR_reporting >= 1 && SPR_reporting <= 3)
{
- SPR_reporting = 4;
- warnstream << "Change SPR_reporting to 4 because benchmarks are off";
+ SPR_reporting = 5;
+ warnstream << "Change SPR_reporting to 5 (raw %SPR) because benchmarks were not requested";
write_message(WARN, 0);
}
}
@@ -5035,10 +5052,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_readstarter.tpl b/SS_readstarter.tpl
index c9b7724e..348b4a70 100644
--- a/SS_readstarter.tpl
+++ b/SS_readstarter.tpl
@@ -643,7 +643,7 @@
int depletion_basis;
int depletion_multi;
int depletion_log;
- init_number depletion_basis_rd; // 0=skip; 1=B0; 2=Bmsy; 3=B_styr; 4=B_endyr; 5=dynamic_Bzero; values >=11 invoke multiyr with 10's digit; append .1 to invoke log(ratio) with hundreds digit
+ init_number depletion_basis_rd; // 0=skip; 1=B0; 2=Bmsy; 3=B_styr; 4=B_endyr; 5=dynamic_Bzero; 6=Bmark_SSB_unf; values >=11 invoke multiyr with 10's digit; append .1 to invoke log(ratio) with hundreds digit
LOCAL_CALCS
// clang-format on
echoinput << depletion_basis_rd << " depletion_basis as read; this is also known as Bratio and is a std quantity; has multi-yr and log(ratio) options" << endl;
@@ -669,7 +669,7 @@
init_number depletion_level;
!!echoinput << depletion_level << " depletion_level" << endl;
- init_int SPR_reporting; // 0=skip; 1=SPR; 2=SPR_MSY; 3=SPR_Btarget; 4=(1-SPR)
+ init_int SPR_reporting; // 0=skip; 1=SPR; 2=SPR_MSY; 3=SPR_Btarget; 4=(1-SPR); 5=SPR
!!echoinput << SPR_reporting << " SPR_reporting" << endl;
init_int F_reporting; // 0=skip; 1=exploit(Bio); 2=exploit(Num); 3=sum(frates); 4=true F for range of ages; 5=unweighted avg F for range of ages
LOCAL_CALCS
@@ -705,6 +705,8 @@
int ender;
int irand_seed;
int irand_seed_rd;
+ int timevary_bio_4SRR; // flag in 3.30.24 for impact of timevary biology on benchmark SRR calculations
+ int timevary_bio_4SRR_rd; // flag in 3.30.24 for impact of timevary biology on benchmark SRR calculations
int F_std_multi; // for multi-year averaging of F_std
int F_std_log; // for log(ratio) of F_std
int F_std_basis;
@@ -806,6 +808,24 @@
echoinput << "random number seed: " << irand_seed << endl;
tempin = 0;
}
+
+ echoinput << "now read flag for dealing with impact of time-varying biology on benchmark SRR calculations" << endl;
+ timevary_bio_4SRR = 0;
+ timevary_bio_4SRR_rd = -1;
+ *(ad_comm::global_datafile) >> tempin;
+ if (tempin == 3.30) // starter file does not contain the new line for timevary_bio_4SRR, so assign default
+ {
+ ender = 1;
+ timevary_bio_4SRR = 0;
+ }
+ else // new input line beginning 3.30.24
+ {
+ timevary_bio_4SRR = int(tempin);
+ timevary_bio_4SRR_rd = 0; // indicates that line was read
+ echoinput << "Compatibility flag for legacy (0) vs improved (1) impact of timevary biology on benchmark SRR calcs: " << timevary_bio_4SRR << endl;
+ tempin = 0;
+ }
+
if (ender == 0)
{
*(ad_comm::global_datafile) >> tempin;
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 93d3fe73..eaf4c980 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 dvar_vector& SRparm_work, const prevariable& SSB_virgin_use, const prevariable& Recr_virgin_use, const prevariable& SSB_current)
{
RETURN_ARRAYS_INCREMENT();
dvariable NewRecruits;
@@ -22,14 +22,15 @@ 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
- regime_change = SR_parm_work(N_SRparm2 - 1); // this is a persistent deviation off the S/R curve
+ regime_change = SRparm_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
@@ -41,82 +42,77 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
// SS_Label_43.3.2 Ricker
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));
+ steepness = SRparm_work(2);
+ 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);
+ steepness = SRparm_work(2);
+ 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(SRparm_work(3));
+ dvariable beta = mfexp(SRparm_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
case 5: // hockey stick where "steepness" is now the fraction of B0 below which recruitment declines linearly
// 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);
+ steepness = SRparm_work(2);
+ temp = SRparm_work(3) * Recr_virgin_use + SSB_curr_adj / (steepness * SSB_virgin_use) * (Recr_virgin_use - SRparm_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;
}
// SS_Label_43.3.6 Beverton-Holt, with constraint to have constant R about Bzero
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)
+ steepness = SRparm_work(2);
+// 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
- // Sfrac=SR_parm(2);
- SRZ_0 = log(1.0 / (SSB_virgin_adj / Recr_virgin_adj));
- steepness = SR_parm_work(2);
+ // Pups_0=SSB_virgin_use; // total population fecundity is the number of pups produced
+ // Sfrac=SRparm(2);
+ SRZ_0 = log(1.0 / (SSB_virgin_use / Recr_virgin_use));
+ steepness = SRparm_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), SRparm_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
- // if(SR_env_target==1) SRZ_surv*=mfexp(SR_parm(N_SRparm2-2)* env_data(y,SR_env_link)); // environ effect on survival
+ // if(SR_env_target==1) SRZ_surv*=mfexp(SRparm(N_SRparm2-2)* env_data(y,SR_env_link)); // environ effect on survival
if (recdev_cycle > 0)
{
gg = y - (styr + (int((y - styr) / recdev_cycle)) * recdev_cycle) + 1;
@@ -145,12 +141,12 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
// SS_Label_43.3.8 Shepherd
case 8: // Shepherd 3-parameter SRR. per Punt & Cope 2017
{
- Shepherd_c = SR_parm_work(3);
- Shepherd_c2 = pow(0.2, SR_parm_work(3));
+ Shepherd_c = SRparm_work(3);
+ Shepherd_c2 = pow(0.2, SRparm_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) /
+ steepness = 0.2 + (SRparm_work(2) - 0.2) / (0.8) * (Hupper - 0.2);
+ 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;
}
@@ -158,21 +154,22 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
// SS_Label_43.3.8 Ricker-power
case 9: // Ricker power 3-parameter SRR. per Punt & Cope 2017
{
- steepness = SR_parm_work(2);
- dvariable RkrPower = SR_parm_work(3);
- temp = SSB_curr_adj / SSB_virgin_adj;
+ steepness = SRparm_work(2);
+ dvariable RkrPower = SRparm_work(3);
+ 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
@@ -180,7 +177,7 @@ FUNCTION void apply_recdev(prevariable& NewRecruits, const prevariable& Recr_vir
// exp_rec(y,2) is with regime shift or other env effect;
// exp_rec(y,3) is with bias adjustment
// exp_rec(y,4) is with dev
- regime_change = SR_parm_work(N_SRparm2 - 1); // this is a persistent deviation off the S/R curve
+ regime_change = SRparm_work(N_SRparm2 - 1); // this is a persistent deviation off the S/R curve
if (recdev_cycle > 0)
{
@@ -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
@@ -285,9 +282,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
dvariable SRZ_0;
dvariable srz_min;
dvariable SRZ_surv;
+ dvariable SSBpR_virgin_use;
- 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
+ SSBpR_virgin_use = SSB_virgin_use / Recr_virgin_use;
+ 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 +295,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_use * (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)
report2 << recdev_cycle_parm << " ";
if (recdev_do_early > 0)
@@ -180,14 +180,14 @@ FUNCTION void write_summaryoutput()
}
}
- report2 << runnumber << " SR_parm ";
+ report2 << runnumber << " SRparm ";
for (i = 1; i <= N_SRparm3 + recdev_cycle; i++)
{
NP++;
report2 << " " << ParmLabel(NP);
}
report2 << endl
- << runnumber << " SR_parm " << SR_parm << " ";
+ << runnumber << " SRparm " << SRparm << " ";
if (recdev_cycle > 0)
report2 << recdev_cycle_parm;
report2 << endl;
@@ -453,8 +453,8 @@ FUNCTION void write_SS_summary()
for (j = 1; j <= N_SRparm3; j++)
{
NP++;
- SS_smry << ParmLabel(NP) << " " << SR_parm(j) << " ";
- if (active(SR_parm(j)))
+ SS_smry << ParmLabel(NP) << " " << SRparm(j) << " ";
+ if (active(SRparm(j)))
{
active_count++;
SS_smry << CoVar(active_count, 1) << " Act ";
@@ -463,7 +463,7 @@ FUNCTION void write_SS_summary()
{
SS_smry << 0.0 << " Fix ";
}
- SS_smry << (SR_parm(j) - SR_parm_LO(j)) / (SR_parm_HI(j) - SR_parm_LO(j) + 1.0e-6) << endl;
+ SS_smry << (SRparm(j) - SRparm_LO(j)) / (SRparm_HI(j) - SRparm_LO(j) + 1.0e-6) << endl;
}
if (recdev_cycle > 0)
@@ -989,7 +989,7 @@ FUNCTION void write_rebuilder_output()
rebuilder << SSB_yr(styr - 2) << " " << SSB_yr(styr, k) << " #SpawnBio" << endl;
//i. steepness; SigmaR; rho
- rebuilder << SR_parm(2) << " " << sigmaR << " " << SR_parm(N_SRparm2) << " # spawn-recr steepness, sigmaR, autocorr" << endl;
+ rebuilder << SRparm(2) << " " << sigmaR << " " << SRparm(N_SRparm2) << " # spawn-recr steepness, sigmaR, autocorr" << endl;
if (mceval_counter == 0)
{
@@ -1011,7 +1011,7 @@ FUNCTION void write_rebuilder_output()
rebuild_dat << exp_rec(y, 4) << " ";
}
rebuild_dat << " #recruits; first value is R0 (virgin)" << endl;
- rebuild_dat << SSB_yr(styr - 2) << " " << SSB_yr(styr, k) << " #spbio; first value is SSB_virgin (virgin)" << endl;
+ rebuild_dat << SSB_yr(styr - 2) << " " << SSB_yr(styr, k) << " #SSBio; first value is SSB_virgin (virgin)" << endl;
rebuild_dat << 1 << " ";
for (y = styr; y <= k; y++)
rebuild_dat << 0 << " ";
@@ -1042,7 +1042,7 @@ FUNCTION void write_rebuilder_output()
rebuild_dat << "# Which probability to product detailed results for (1=0.5; 2=0.6; etc.)" << endl
<< 3 << endl;
rebuild_dat << "# Steepness sigma-R Auto-correlation" << endl
- << SR_parm(2) << " " << sigmaR << " " << 0 << endl;
+ << SRparm(2) << " " << sigmaR << " " << 0 << endl;
rebuild_dat << "# Target SPR rate (FMSY Proxy); manually change to SPR_MSY if not using SPR_target" << endl
<< SPR_target << endl;
rebuild_dat << "# Discount rate (for cumulative catch)" << endl
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index f51aa464..5f6fa39b 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1,6 +1,6 @@
// SS_Label_file #19. **SS_write_report.tpl**
// SS_Label_file # * write_bigoutput() // produces *report.sso* and *compreport.sso*
-// SS_Label_file # * SPR_profile() // calls Do_Equil_Calc() and Equil_Spawn_Recr_Fxn() over a range of F to get SPR, YPR, and SSB and catch curves
+// SS_Label_file # * SPR_profile() // calls SSBpR_Calc() and Equil_Spawn_Recr_Fxn() over a range of F to get SPR, YPR, and SSB and catch curves
// SS_Label_file # * global_MSY() // similar to SPR_profile but first changes all selectivities to knife edge and profiles on age-at-entry
// SS_Label_file #
@@ -77,7 +77,9 @@ FUNCTION void write_bigoutput()
SS2out << "Data_File: " << datfilename << endl;
SS2out << "Control_File: " << ctlfilename << endl;
if (readparfile >= 1)
- SS2out << "Start_parm_values_from_SS.PAR" << endl;
+ {SS2out << "Start_parm_values_from_SS.PAR" << endl;}
+ else
+ {SS2out << "Start_parm_values_from_control_file" << endl;}
SS2out << endl
<< "Convergence_Level: " << objective_function_value::pobjfun->gmax << " is_final_gradient" << endl;
temp = get_ln_det_value();
@@ -647,12 +649,12 @@ FUNCTION void write_bigoutput()
{
NP++;
Activ = 0;
- if (active(SR_parm(j)))
+ if (active(SRparm(j)))
{
active_count++;
Activ = 1;
}
- Report_Parm(NP, active_count, Activ, SR_parm(j), SR_parm_LO(j), SR_parm_HI(j), SR_parm_RD(j), SR_parm_use(j), SR_parm_PR(j), SR_parm_CV(j), SR_parm_PRtype(j), SR_parm_PH(j), SR_parm_Like(j));
+ Report_Parm(NP, active_count, Activ, SRparm(j), SRparm_LO(j), SRparm_HI(j), SRparm_RD(j), SRparm_use(j), SRparm_PR(j), SRparm_CV(j), SRparm_PRtype(j), SRparm_PH(j), SRparm_Like(j));
}
if (recdev_cycle > 0)
@@ -1516,8 +1518,8 @@ FUNCTION void write_bigoutput()
if (s == spawn_seas)
{
temp = sum(SSB_pop_gp(y, p));
- if (Hermaphro_maleSPB > 0)
- temp += Hermaphro_maleSPB * sum(MaleSPB(y, p));
+ if (Hermaphro_maleSSB > 0)
+ temp += Hermaphro_maleSSB * sum(MaleSSB(y, p));
SS2out << temp;
}
else
@@ -1529,7 +1531,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 +1712,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);
@@ -1778,11 +1780,11 @@ FUNCTION void write_bigoutput()
Durbin /= (var + 1.0e-09);
var /= (n_rmse(1) + 1.0e-09);
- dvariable steepness = SR_parm(2);
+ dvariable steepness = SRparm(2);
SS2out << endl
<< pick_report_name(19);
SS2out << " Function: " << SR_fxn << " RecDev_method: " << do_recdev << " sum_recdev: " << sum_recdev << endl
- << SR_parm(1) << " Ln(R0) " << mfexp(SR_parm(1)) << endl
+ << SRparm(1) << " Ln(R0) " << mfexp(SRparm(1)) << endl
<< steepness << " steepness" << endl
<< Bmsy / SSB_virgin << " Bmsy/Bzero ";
if (SR_fxn == 8)
@@ -1790,22 +1792,22 @@ FUNCTION void write_bigoutput()
dvariable Shepherd_c;
dvariable Shepherd_c2;
dvariable Hupper;
- Shepherd_c = SR_parm(3);
+ Shepherd_c = SRparm(3);
Shepherd_c2 = pow(0.2, Shepherd_c);
Hupper = 1.0 / (5.0 * Shepherd_c2);
- temp = 0.2 + (SR_parm(2) - 0.2) / (0.8) * (Hupper - 0.2);
+ temp = 0.2 + (SRparm(2) - 0.2) / (0.8) * (Hupper - 0.2);
SS2out << " Shepherd_c: " << Shepherd_c << " steepness_limit: " << Hupper << " Adjusted_steepness: " << temp;
}
else if (SR_fxn == 9)
{
- SS2out << " Ricker_Power: " << SR_parm(3);
+ SS2out << " Ricker_Power: " << SRparm(3);
}
SS2out << endl;
SS2out << sigmaR << " sigmaR" << endl;
SS2out << init_equ_steepness << " # 0/1 to use steepness in initial equ recruitment calculation" << endl;
- SS2out << SR_parm(N_SRparm2 - 1) << " init_eq: see below" << endl
+ SS2out << SRparm(N_SRparm2 - 1) << " init_eq: see below" << endl
<< recdev_start << " " << recdev_end << " main_recdev:start_end" << endl
<< recdev_adj(1) << " " << recdev_adj(2, 5) << " breakpoints_for_bias_adjustment_ramp " << endl;
@@ -1889,21 +1891,189 @@ FUNCTION void write_bigoutput()
SS2out << endl;
}
+ SS2out << endl << "#Expanded_Spawn_Recr_report" << endl << pick_report_name(19) << endl;
+ SS2out << SR_fxn << " # SR_Function" << endl;
+ SS2out << N_SRparm2 << " # N_SRparms" << endl;
+ SS2out << "#" << endl << "#_SRparm parm_label value phase TV_year" << endl;
+ for (int j = 1; j <=N_SRparm2; j++)
+ {
+ SS2out << "# " << j << " " << ParmLabel(firstSRparm + j) << " " << SRparm(j) << " " << SRparm_PH(j);
+ if (SRparm_timevary(j) > 0 && j <= 4 ) // timevary SRparm exists
+ {SS2out << " " << timevary_SRparm_first_yr;} else {SS2out << " -";}
+ SS2out << endl;
+ }
+ SS2out << "#" << endl;
+ if (SRparm_timevary(N_SRparm2 - 1) > 0) // timevary regime exists
+ {SS2out << " #_Regime_used_to_offset_from_SRR";}
+ SS2out << timevary_bio_4SRR << " # timevary_bio_4SRR #_Compatibility_flag_for_legacy_(0)_vs_improved_(1)_impact_of_timevary_biology_on_benchmark_SRR_calcs" << endl;
+ if( timevary_MG_firstyr == YrMax)
+ { SS2out << " #_No_timevary_MGparm" << endl; }
+ else
+ { SS2out << timevary_MG_firstyr << " #_first year_timevary_MGparm_(or_any_year_EWAA) " << endl; }
+
+ SS2out << "#" << endl << "Quantities for MSY and other benchmark calculations " << endl
+ << "Benchmark_index: 1 2 3 4 5 6 7 8 9 10" << endl
+ << "Benchmark_label: beg_bio end_bio beg_selex end_selex beg_relF end_relF beg_recr_dist end_recr_dist beg_SRparm end_SRparm" << endl
+ << "Benchmark_years: " << Bmark_Yr << endl;
+ for ( int k = 1; k <=9; k+=2)
+ { if (Bmark_Yr(k+1) > Bmark_Yr(k))
+ {SS2out << "#_NOTE:_using_range_of_years_can_reduce_standard_error_of_result;_do_this_only_when_timevarying_makes_necessary: " << k << " " << Bmark_Yr(k) << " " << Bmark_Yr(k+1) << endl;}
+ }
+ SS2out << "SSBpR0_virgin: " << SSBpR_virgin << " #_uses_biology_from_start_year: " << styr < (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 << init_equ_steepness << " # Initial_equilibrium:_0/1_to_use_spawner-recruitment_in_initial_equ_recruitment_calculation" << endl << "#" << endl;
+ SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SSBpR(yr) " << endl;
+
+ y = styr - 2;
+ SS2out << "Virg " << SSB_yr(y) << " " << exp_rec(y) << " _ " << 0.0 << " Virg " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " _ " << SSBpR_virgin << endl;
+ y = styr - 1;
+ SS2out << "Init " << SSB_yr(y) << " " << exp_rec(y) << " _ " << 0.0 << " Init " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " _ " << SSBpR_virgin << endl;
+
+ if (recdev_first < styr)
+ {
+ for (y = recdev_first; y <= styr - 1; y++)
+ {
+ SS2out << y << " " << SSB_yr(styr - 1) << " " << exp_rec(styr - 1, 1) << " " << exp_rec(styr - 1, 2) << " " << exp_rec(styr - 1, 3) * mfexp(-biasadj(y) * half_sigmaRsq) << " " << exp_rec(styr - 1, 3) * mfexp(recdev(y) - biasadj(y) * half_sigmaRsq) << " "
+ << recdev(y) << " " << biasadj(y) << " Init_age " << SSB_B_yr(styr - 1) << " " << SSB_N_yr(styr - 1) << " " << recdev(y) << endl; // newdev approach uses devs for initial agecomp directly
+ }
+ }
+ for (y = styr; y <= YrMax; y++)
+ {
+ SS2out << y << " " << SSB_yr(y) << " " << exp_rec(y) << " ";
+ if (recdev_do_early > 0 && y >= recdev_early_start && y <= recdev_early_end)
+ {
+ SS2out << log(exp_rec(y, 4) / exp_rec(y, 3)) << " " << biasadj(y) << " Early " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " " << recdev(y);
+ }
+ else if (y >= recdev_start && y <= recdev_end)
+ {
+ SS2out << log(exp_rec(y, 4) / exp_rec(y, 3)) << " " << biasadj(y) << " Main " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " " << recdev(y);
+ }
+ else if (Do_Forecast > 0 && y > recdev_end)
+ {
+ SS2out << log(exp_rec(y, 4) / exp_rec(y, 3)) << " " << biasadj(y);
+ if (y <= endyr)
+ {
+ SS2out << " Late ";
+ }
+ else
+ {
+ SS2out << " Fore ";
+ }
+ SS2out << SSB_B_yr(y) << " " << SSB_N_yr(y) << " ";
+ if (do_recdev > 0)
+ {
+ SS2out << Fcast_recruitments(y);
+ }
+ else
+ {
+ SS2out << " 0.0";
+ }
+ }
+ else
+ {
+ SS2out << " _ _ Fixed " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " _ ";
+ }
+ SS2out << " " << Smry_Table(y, 11) / Recr_virgin << endl;
+ }
+
// REPORT_KEYWORD SPAWN_RECR_CURVE
if (pick_report_use(20) == "Y")
{
{
- SS2out << endl
- << pick_report_name(20) << endl;
+ SRparm_work = SRparm_byyr(styr);
+ y = styr;
+ SS2out << endl
+ << pick_report_name(20) << endl
+ << "# using_virgin_SR_parameters: " << SRparm_work << endl;
SS2out << "SSB/SSB_virgin SSB Recruitment" << endl;
- y = styr;
- SR_parm_work = SR_parm_byyr(styr);
for (f = 1; f <= 120; f++)
{
SSB_current = double(f) / 100. * SSB_virgin;
- temp = Spawn_Recr(SSB_virgin, Recr_virgin, SSB_current);
+ temp = Spawn_Recr(SRparm_work, SSB_virgin, Recr_virgin, SSB_current);
SS2out << SSB_current / SSB_virgin << " " << SSB_current << " " << temp << endl;
}
+ SS2out << endl
+ << "SPAWN_RECR_CURVE report:20 Benchmark" << endl // revise this name per r4ss needs
+ << "# using_benchmark_SR_parameters: " << SRparm_bench << endl;
+ SS2out << "SSB/SSB_benchmark SSB Recruitment" << endl;
+ for (f = 1; f <= 120; f++)
+ {
+ SSB_current = double(f) / 100. * SSB0_4_SRR;
+ temp = Spawn_Recr(SRparm_bench, SSB0_4_SRR, R0_4_SRR, SSB_current);
+ SS2out << SSB_current / SSB0_4_SRR << " " << SSB_current << " " << temp << endl;
+ }
+
}
}
@@ -4613,41 +4783,20 @@ FUNCTION void SPR_profile()
bio_t_base = styr + (bio_yr - styr) * nseas - 1;
// SPAWN-RECR: call make_fecundity for benchmark bio for SPR loop
- // this code section that creates fecundity and selectivity seems antiquated; why is it different from the averages used for benchmark?
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)
+ Wt_Age_beg(s) = Wt_Age_t(t, 0); // used for smrybio
+ Wt_Age_mid(s) = Wt_Age_t(t, -1); // used in global MSY
+ if (s == spawn_seas)
{
- 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();
- }
+ fec = Wt_Age_t(t, -2);
+ SS2out << " repro_output for 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 ";
+// do not recalculate here so force using values from benchmark
+ SS2out << "unfished values for SRR: SSB " << SSB0_4_SRR << " R " << R0_4_SRR << " SSBpR " << " SSBpR: " << SSB0_4_SRR / R0_4_SRR << endl;
+ SS2out << "SPRloop Iter Bycatch Fmult F_std SSBpR YpR_dead YpR_dead*Recr YpR_ret*Recr Revenue Cost Profit SSB Recruits SSB/Bzero Tot_Catch ";
for (f = 1; f <= Nfleet; f++)
{
if (fleet_type(f) <= 2)
@@ -4681,7 +4830,7 @@ FUNCTION void SPR_profile()
SPRloop1_end = 7;
}
int SPRloops;
- Do_Equil_Calc(equ_Recr);
+ SSBpR_Calc(equ_Recr);
if (N_bycatch == 0)
{
k = 0;
@@ -4691,20 +4840,28 @@ FUNCTION void SPR_profile()
k = 1;
}
for (int with_BYC = 0; with_BYC <= k; with_BYC++)
- for (int SPRloop1 = 0; SPRloop1 <= SPRloop1_end; SPRloop1++)
+ for (int SPRloop1 = -1; SPRloop1 <= SPRloop1_end; SPRloop1++)
{
Fmultchanger1 = value(pow(0.0001 / Fcrash, 0.025));
Fmultchanger2 = value(Fcrash / 39.);
SPRloops = 40;
switch (SPRloop1)
{
+ case -1:
+ {
+ SPRloops = 1;
+ Fmult2 = 0.0;
+ break;
+ }
case 0:
{
+ SPRloops = 40;
Fmult2 = maxpossF;
break;
}
case 1:
{
+ SPRloops = 40;
Fmult2 = Fcrash;
break;
}
@@ -4793,10 +4950,10 @@ FUNCTION void SPR_profile()
}
Fishon = 1;
- Do_Equil_Calc(equ_Recr);
+ SSBpR_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(SRparm_bench , SSB0_4_SRR, R0_4_SRR, 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)
@@ -4806,9 +4963,9 @@ FUNCTION void SPR_profile()
if (SPRloop1 == 0)
Fcrash = Fmult2;
}
- SS2out << SPRloop1 << " " << SPRloop << " " << with_BYC << " " << Fmult2 << " " << equ_F_std << " " << SSB_equil / (SSB_unf / Recr_unf) << " " << YPR_dead << " "
+ SS2out << SPRloop1 << " " << SPRloop << " " << with_BYC << " " << Fmult2 << " " << equ_F_std << " " << SSB_equil / (SSB0_4_SRR / R0_4_SRR) << " " << YPR_dead << " "
<< YPR_dead * Btgt_prof_rec << " " << YPR_ret * Btgt_prof_rec << " " << (PricePerF * YPR_val_vec) * Btgt_prof_rec
- << " " << Cost << " " << (PricePerF * YPR_val_vec) * Btgt_prof_rec - Cost << " " << Btgt_prof << " " << Btgt_prof_rec << " " << Btgt_prof / SSB_unf
+ << " " << Cost << " " << (PricePerF * YPR_val_vec) * Btgt_prof_rec - Cost << " " << Btgt_prof << " " << Btgt_prof_rec << " " << Btgt_prof / SSB0_4_SRR
<< " " << value(sum(equ_catch_fleet(2)) * Btgt_prof_rec);
for (f = 1; f <= Nfleet; f++)
if (fleet_type(f) <= 2)
@@ -4896,7 +5053,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;
@@ -4975,7 +5132,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 04c450d2..99ff2299 100644
--- a/SS_write_ssnew.tpl
+++ b/SS_write_ssnew.tpl
@@ -1509,11 +1509,11 @@ FUNCTION void write_nucontrol()
<< version_info2 << endl;
if (N_SC > 0)
NuStart << Starter_Comments << endl;
- NuStart << datfilename << endl
- << ctlfilename << endl;
- NuStart << readparfile << " # 0=use init values in control file; 1=use ss.par" << endl;
- NuStart << rundetail << " # run display detail (0 = minimal; 1=one line per iter; 2=each logL)" << endl;
- NuStart << reportdetail << " # detailed output (0=minimal for data-limited, 1=high (w/ wtatage.ss_new), 2=brief, 3=custom) " << endl;
+ NuStart << datfilename << " #_datfile" << endl
+ << ctlfilename << " #_ctlfile" << endl;
+ NuStart << readparfile << " #_init_values_src: 0 (use init values in control file); 1 (use ss3.par)" << endl;
+ NuStart << rundetail << " #_screen_display: 0 (minimal); 1 (one line per iter); 2 (each logL)" << endl;
+ NuStart << reportdetail << " #_report_table_selection: 0 (minimal; no wtatage.ss_new); 1 (all tables); 2 (brief), 3 (custom, read list) " << endl;
if (reportdetail == 3)
{
NuStart << "# custom report options: -100 to start with minimal; -101 to start with all; -number to remove, +number to add, -999 to end" << endl;
@@ -1524,34 +1524,36 @@ FUNCTION void write_nucontrol()
}
else
{
- NuStart << "#COND: custom report options: -100 to start with minimal; -101 to start with all; -number to remove, +number to add, -999 to end" << endl;
+ NuStart << "# COND: custom report options: -100 to start with minimal; -101 to start with all; -number to remove, +number to add, -999 to end" << endl;
}
- NuStart << docheckup << " # write 1st iteration details to echoinput.sso file (0,1) " << endl;
- NuStart << Do_ParmTrace << " # write parm values to ParmTrace.sso (0=no,1=good,active; 2=good,all; 3=every_iter,all_parms; 4=every,active)" << endl;
- NuStart << Do_CumReport << " # write to cumreport.sso (0=no,1=like×eries; 2=add survey fits)" << endl;
- NuStart << Do_all_priors << " # Include prior_like for non-estimated parameters (0,1) " << endl;
- NuStart << SoftBound << " # Use Soft Boundaries to aid convergence (0,1) (recommended)" << endl;
+ NuStart << docheckup << " #_checkup: write more 1st iteration age-specific details to echoinput.sso file (0,1) " << endl;
+ NuStart << Do_ParmTrace << " #_parmtrace: write parm values to ParmTrace.sso: 0 (no); 1 (good_iter,active_parms); 2 (good,all); 3 (every,all); 4 (every,active)" << endl;
+ NuStart << Do_CumReport << " #_cumreport: write to cumreport.sso: 0 (no); 1 (like×eries); 2 (add survey fits)" << endl;
+ NuStart << Do_all_priors << " #_prior_like: include prior_like for non-estimated parameters (0,1) " << endl;
+ NuStart << SoftBound << " #_soft_bounds: Use Soft Boundaries to aid convergence (0,1) (recommended)" << endl;
NuStart << "#" << endl
- << N_nudata_read << " # Number of datafiles to produce: 0 turns off all *.ss_new; 1st is data_echo.ss_new, 2nd is data_expval.ss, 3rd and higher are data_boot_**N.ss," << endl;
- NuStart << Turn_off_phase_rd << " # Turn off estimation for parameters entering after this phase" << endl;
+ << N_nudata_read << " #_N_bootstraps: Number of datafiles to produce: 0 turns off all *.ss_new; 1st is data_echo.ss_new, 2nd is data_expval.ss, 3rd and higher are data_boot_**N.ss," << endl;
+ NuStart << Turn_off_phase_rd << " #_last_estimation_phase: turn off estimation for parameters entering after this phase" << endl;
NuStart << "#" << endl
- << burn_intvl << " # MCeval burn interval" << endl;
- NuStart << thin_intvl << " # MCeval thin interval" << endl;
- NuStart << jitter << " # jitter initial parm value by this fraction" << endl;
- NuStart << STD_Yr_min_rd << " # min year for sdreport outputs (-1 for styr); #_" << STD_Yr_min << endl;
- NuStart << STD_Yr_max_rd << " # max year for sdreport outputs (-1 for endyr+1; -2 for endyr+Nforecastyrs); #_" << STD_Yr_max << endl;
- NuStart << N_STD_Yr_RD << " # N individual STD years " << endl;
+ << burn_intvl << " #_MCMCburn" << endl;
+ NuStart << thin_intvl << " #_MCMCthin" << endl;
+ NuStart << jitter << " # jitter_fraction: jitter within parameter bounds" << endl;
+
+ NuStart << STD_Yr_min_rd << " #_minyr_sdreport: min year for sdreport outputs (-1 for styr); #_" << STD_Yr_min << endl;
+ NuStart << STD_Yr_max_rd << " #_maxyr_sdreport: max year for sdreport outputs (-1 for endyr+1; -2 for endyr+Nforecastyrs); #_" << STD_Yr_max << endl;
+ NuStart << N_STD_Yr_RD << " #_N_STD_yrs: N individual STD years " << endl;
NuStart << "#COND: vector of year values if N>0" << endl
<< STD_Yr_RD << endl;
- 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_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;
+ NuStart << final_conv << " #_converge_criterion: (e.g. 1.0e-04) " << endl;
+ NuStart << retro_yr - endyr << " #_retro_yr: retrospective year relative to end year (e.g. -4)" << endl;
+ NuStart << Smry_Age << " #_min_age_summary_bio" << endl;
+ NuStart << depletion_basis_rd << " #_depl_basis: Depletion denom is: 0=skip; 1=X*SSB_virgin; 2=X*SSB_msy; 3=X*SSB_styr; 4=X*SSB_endyr; 5=X*dyn_Bzero; 6=X*SSB_unf_bmark; 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 << "# If value = 1, then Btarget in benchmark will be a fraction of SSB_virgin, else will be a fraction of SSB_benchmark" << endl;
+ NuStart << depletion_level << " #_depl_denom_frac: fraction (X) for Depletion denominator (e.g. 0.4)" << endl;
+ NuStart << SPR_reporting << " #_SPR_basis: 0 (skip); 1 (1-SPR)/(1-SPR_tgt); 2 (1-SPR)/(1-SPR_MSY); 3 (1-SPR)/(1-SPR_Btarget); 4 (1-SPR); 5 (SPR)" << endl;
+ NuStart << F_reporting << " # F_std_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;
if (F_reporting == 4 || F_reporting == 5)
{
NuStart << F_reporting_ages << " # min and max age over which mean F will be calculated, with F=Z-M" << endl;
@@ -1560,11 +1562,12 @@ FUNCTION void write_nucontrol()
{
NuStart << "#COND 10 15 #_min and max age over which mean F will be calculated with F_reporting=4 or 5" << endl;
}
- NuStart << F_std_basis_rd << " # F_std_scaling: 0=no scaling; 1=F/Fspr; 2=F/Fmsy; 3=F/Fbtgt; where F means annual F_std, Fmsy means F_std@msy; values >=11 invoke N multiyr using 10s and 100s digit; append .1 to invoke log(ratio)" << endl;
- NuStart << double(mcmc_output_detail) + MCMC_bump << " # MCMC output detail: integer part (0=default; 1=adds obj func components; 2= +write_report_for_each_mceval); and decimal part (added to SR_LN(R0) on first call to mcmc)" << endl;
- NuStart << ALK_tolerance << " # ALK tolerance ***disabled in code" << endl;
- NuStart << irand_seed_rd << " # random number seed for bootstrap data (-1 to use long(time) as seed): # " << irand_seed << endl;
- NuStart << "3.30 # check value for end of file and for version control" << endl;
+ NuStart << F_std_basis_rd << " # F_std_basis: 0=no scaling; 1 (F/Fspr); 2 (F/Fmsy); 3 (F/Fbtgt); where F means annual F_std, Fmsy means F_std@msy; values >=11 invoke N multiyr using 10s and 100s digit; append .1 to invoke log(ratio)" << endl;
+ NuStart << double(mcmc_output_detail) + MCMC_bump << " #_MCMC_output_detail: integer part (0=default; 1=adds obj func components; 2= +write_report_for_each_mceval); and decimal part (added to SR_LN(R0) on first call to mcmc)" << endl;
+ NuStart << ALK_tolerance << " #_deprecated: ALK tolerance ***disabled in code" << endl;
+ NuStart << irand_seed_rd << " #_seed: random number seed for bootstrap data (-1 to use long(time) as seed): # " << irand_seed << endl;
+ NuStart << timevary_bio_4SRR << " #_Compatibility: flag for legacy (0) vs improved (1) impact of timevary biology on benchmark SRR calcs >=3.30.24" << endl;
+ NuStart << "3.30 #_final: check value for end of file and for version control" << endl;
NuStart.close();
echoinput << "Write forecast.ss_new file " << endl;
@@ -1574,7 +1577,7 @@ FUNCTION void write_nucontrol()
if (N_FC > 0)
NuFore << Forecast_Comments << endl;
NuFore << "# for all year entries except rebuilder; enter either: actual year, -999 for styr, 0 for endyr, neg number for rel. endyr" << endl;
- NuFore << Do_Benchmark << " # Benchmarks: 0=skip; 1=calc F_spr,F_btgt,F_msy; 2=calc F_spr,F0.1,F_msy; 3=add F_Blimit; " << endl;
+ NuFore << Do_Benchmark << " # Benchmarks: 0=skip; 1=calc F_spr,F_Btgt,F_msy; 2=calc F_spr,F0.1,F_msy; 3=add F_Blimit; " << endl;
NuFore << Do_MSY << " # Do_MSY: 1= set to F(SPR); 2=calc F(MSY); 3=set to F(Btgt) or F0.1; 4=set to F(endyr); 5=calc F(MEY) with MSY_unit options" << endl;
NuFore << "# if Do_MSY=5, enter MSY_Units; then list fleet_ID, cost/F, price/mt, include_in_Fmey_scaling; # -fleet_ID to fill; -9999 to terminate" << endl;
if (Do_MSY == 5)
@@ -1591,14 +1594,14 @@ FUNCTION void write_nucontrol()
}
NuFore << SPR_target << " # SPR target (e.g. 0.40)" << endl;
- NuFore << BTGT_target << " # Biomass target (e.g. 0.40)" << endl;
- if (Do_Benchmark == 3)
- NuFore << Blim_frac << " # COND: Do_Benchmark==3; Blimit as fraction of Bmsy (neg value to use as frac of Bzero) (e.g. 0.50)" << endl;
- NuFore << "#_Bmark_years: beg_bio, end_bio, beg_selex, end_selex, beg_relF, end_relF, beg_recr_dist, end_recr_dist, beg_SRparm, end_SRparm (enter actual year, or values of 0 or -integer to be rel. endyr)" << endl
+ NuFore << BTGT_target << " # Biomass target (e.g. 0.40) as fraction of SSB_virgin if depletion basis = 1, else as fraction of SSB_unfished in benchmark" << endl;
+ if (Do_Benchmark != 3) NuFore << "#";
+ NuFore << Blim_frac << " # COND: Do_Benchmark==3; Blimit as fraction of Bmsy (neg value to use as frac of SSB_virgin or SSB_unfished) (e.g. 0.50)";
+ NuFore << "#" << endl << "# Bmark_years: beg_bio, end_bio, beg_selex, end_selex, beg_relF, end_relF, beg_recr_dist, end_recr_dist, beg_SRparm, end_SRparm (enter actual year, or values of 0 or -integer to be rel. endyr)" << endl
<< Bmark_Yr_rd << endl
<< "# " << Bmark_Yr << endl;
NuFore << "# value <0 convert to endyr-value; except -999 converts to start_yr; must be >=start_yr and <=endyr" << endl;
- NuFore << Bmark_RelF_Basis << " #Bmark_relF_Basis: 1 = use year range; 2 = set relF same as forecast below" << endl;
+ NuFore << Bmark_RelF_Basis << " # Bmark_relF_Basis: 1 = use year range; 2 = set relF same as forecast below" << endl;
NuFore << "#" << endl
<< Do_Forecast_rd << " # Forecast: -1=none; 0=simple_1yr; 1=F(SPR); 2=F(MSY) 3=F(Btgt) or F0.1; 4=Ave F (uses first-last relF yrs); 5=input annual F scalar" << endl;
NuFore << "# where none and simple require no input after this line; simple sets forecast F same as end year F" << endl;
@@ -1619,7 +1622,7 @@ FUNCTION void write_nucontrol()
}
// else
{ // new list based format for Fcast years
- NuFore << anystring << endl << anystring << "-12345 # code to invoke new format for expanded fcast year controls" << endl
+ NuFore << anystring << endl << anystring << " -12345 # code to invoke new format for expanded fcast year controls" << endl
<< "# biology and selectivity vectors are updated annually in the forecast according to timevary parameters, so check end year of blocks and dev vectors" << endl
<< "# input in this section directs creation of means over historical years to override any time_vary changes" << endl
<< "# Factors implemented so far: 1=M, 4=recr_dist, 5=migration, 10=selectivity, 11=rel_F, 12=recruitment" << endl
@@ -1640,10 +1643,11 @@ FUNCTION void write_nucontrol()
}
NuFore << HarvestPolicy << " # Control rule method (0: none; 1: ramp does catch=f(SSB), buffer on F; 2: ramp does F=f(SSB), buffer on F; 3: ramp does catch=f(SSB), buffer on catch; 4: ramp does F=f(SSB), buffer on catch) " << endl;
- NuFore << "# values for top, bottom and buffer exist, but not used when Policy=0" << endl;
- 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 << "# values for top, bottom and buffer required, but not used when Policy=0" << endl;
+ NuFore << H4010_top_rd << " # Control rule inflection for constant F (as frac of HCR_anchor, see below); must be > control rule cutoff" << endl;
+ NuFore << H4010_bot << " # Control rule cutoff for no F (as frac of HCR_anchor, 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]. -year fills from year to YrMax " << endl;
+ NuFore << "# Also see HCR_anchor below to use virgin vs benchmark SSB or Bmsy as basis for inflection and cutoff" << endl;
if (H4010_scale_rd < 0)
{
j = H4010_scale_vec_rd.size() - 1;
@@ -1665,7 +1669,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 2 uses unfished benchmark SSB (old hardwired approach); 1 = virgin SSB; 3 = BMSY" << endl << "#" << endl;
NuFore << Fcast_Cap_FirstYear << " # FirstYear for caps and allocations (should be after years with fixed inputs) " << endl;
@@ -1941,7 +1945,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;
@@ -2134,20 +2138,20 @@ 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 << 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; 10=B-H_ab" << 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 << " 0 # not_used" << 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);
for (f = 1; f <= N_SRparm2; f++)
{
NP++;
- SR_parm_1(f, 3) = value(SR_parm(f));
+ SRparm_1(f, 3) = value(SRparm(f));
for (j = 1; j <= 6; j++)
- report4 << setw(14) << SR_parm_1(f, j);
+ report4 << setw(14) << SRparm_1(f, j);
for (j = 7; j <= 14; j++)
- report4 << setw(11) << SR_parm_1(f, j);
+ report4 << setw(11) << SRparm_1(f, j);
report4 << " # " << ParmLabel(NP) << endl;
}
report4.unsetf(std::ios_base::fixed);
@@ -2155,7 +2159,7 @@ FUNCTION void write_nucontrol()
if (N_SRparm3 > N_SRparm2)
{
report4 << "# timevary SR parameters" << endl;
- for (f = timevary_parm_start_SR; f <= timevary_parm_cnt_SR; f++)
+ for (f = timevary_SRparm_first; f <= timevary_parm_SR_last; f++)
{
NP++;
timevary_parm_rd[f](3) = value(timevary_parm(f));
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
}