From 45e40fce830b89168f40b1f876baec10fd793f13 Mon Sep 17 00:00:00 2001
From: RickMethot <richard.methot@noaa.gov>
Date: Wed, 8 Feb 2023 16:24:40 -0800
Subject: [PATCH 01/29] first commit for B-H with ab
add S-R option 10 for a,b formulation
---
SS_benchfore.tpl | 12 ++---
SS_param.tpl | 1 +
SS_popdyn.tpl | 33 +++++++++++---
SS_readcontrol_330.tpl | 11 ++++-
SS_recruit.tpl | 85 ++++++++++++++++++-----------------
SS_write_report.tpl | 8 +++-
StockSynthesis.code-workspace | 1 +
7 files changed, 94 insertions(+), 57 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index ff3f5844..327810cc 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -902,7 +902,7 @@ 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
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSB_equil); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Bspr = Equ_SpawnRecr_Result(1);
Bspr_rec = Equ_SpawnRecr_Result(2);
@@ -1017,7 +1017,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
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
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_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
@@ -1111,7 +1111,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SPR_Btgt = SSB_equil / SPR_unfished;
// SPAWN-RECR: calc equil spawn-recr for Btarget calcs; need to make area-specific
SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
yld1(ii) = Equ_SpawnRecr_Result(1);
}
@@ -1331,7 +1331,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
// 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
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Bmsy = Equ_SpawnRecr_Result(1);
Recr_msy = Equ_SpawnRecr_Result(2);
yld1(1) = YPR_opt * Recr_msy;
@@ -1424,7 +1424,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
// 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
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Bmsy = Equ_SpawnRecr_Result(1);
Recr_msy = Equ_SpawnRecr_Result(2);
Profit = (PricePerF * YPR_val_vec) * Recr_msy - Cost;
@@ -1718,7 +1718,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SPR_Btgt2 = SSB_equil / SPR_unfished;
// SPAWN-RECR: calc equil spawn-recr for Btarget calcs; need to make area-specific
SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
yld1(ii) = Equ_SpawnRecr_Result(1);
}
diff --git a/SS_param.tpl b/SS_param.tpl
index d3c43b41..bd21a12b 100644
--- a/SS_param.tpl
+++ b/SS_param.tpl
@@ -416,6 +416,7 @@ PARAMETER_SECTION
// note that bycatch_F(1,Nfleet,1,nseas) has similar role
number alpha;
number beta;
+ number steepness;
number GenTime;
vector cumF(1,gmorph);
vector maxF(1,gmorph);
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index 664bdd5d..70bdc8f9 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -337,17 +337,38 @@ FUNCTION void get_initial_conditions()
Fishon = 0;
virg_fec = fec;
Recr.initialize(); // will store recruitment by area
+
+ // SPAWN-RECR: get expected recruitment globally or by area
+ if (recr_dist_area == 1 || pop == 1) // do global spawn_recruitment calculations
+ {
+ equ_Recr = 1.0;
+ Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation. Returns SPR because R = 1.0
+ SPR_virgin = SSB_equil; // spawners per recruit. Needed for Sr_fxn = 10
+
+ if(SR_fxn == 10) // B-H with a,b
+ {
+ alpha = mfexp(SR_parm(3));
+ beta = mfexp(SR_parm(4));
+ steepness = alpha * SPR_virgin / (4.00 + alpha * SPR_virgin);
+ Recr_virgin = 1.0 / beta * (alpha - ( 1.0 / SPR_virgin));
+ SR_parm(1) = log(Recr_virgin);
+ SR_parm(2) = steepness;
+ // warning<< SR_parm(1) << " Ln(R0) " << mfexp(SR_parm(1)) << endl
+ // << " alpha_beta: "<<alpha<<" "<<beta<<endl
+ // << " spr: "<<SPR_virgin<<endl
+ // << " der_steep: "<< alpha * SPR_virgin / (4.00 + alpha * SPR_virgin)<<endl
+ // << "der_R0: "<< 1.0 / beta * (alpha - ( 1.0 / SPR_virgin))<<endl;
+ }
+ else {
+ Recr_virgin = mfexp(SR_parm(1));
+ }
+
for (int i = 1; i <= N_SRparm2; i++)
{
SR_parm_byyr(eq_yr, i) = SR_parm(i);
SR_parm_virg(i) = SR_parm(i);
SR_parm_work(i) = SR_parm(i);
}
-
- // SPAWN-RECR: get expected recruitment globally or by area
- if (recr_dist_area == 1 || pop == 1) // do global spawn_recruitment calculations
- {
- Recr_virgin = mfexp(SR_parm(1));
equ_Recr = Recr_virgin;
exp_rec(eq_yr, 1) = Recr_virgin; // expected Recr from s-r parms
exp_rec(eq_yr, 2) = Recr_virgin;
@@ -507,7 +528,7 @@ FUNCTION void get_initial_conditions()
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
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm, SSB_virgin, Recr_virgin, SPR_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;
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index b1572680..62efbf6b 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -1937,13 +1937,13 @@
!!// 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<<SR_fxn<<" #_SR_function: 1=NA; 2=Ricker(2 parms); 3=BevHolt(2); 4=SCAA(2); 5=Hockey(3); 6=B-H_flattop(2); 7=Survival(3); 8=Shepherd(3); 9=Ricker_Power(3) "<<endl;
+!!echoinput<<SR_fxn<<" #_SR_function: 1=NA; 2=Ricker(2 parms); 3=BevHolt(2); 4=SCAA(2); 5=Hockey(3); 6=B-H_flattop(2); 7=Survival(3); 8=Shepherd(3); 9=Ricker_Power(3); 10=B-H_a,b(4)"<<endl;
init_int init_equ_steepness;
!!echoinput<<init_equ_steepness<<" # 0/1 to use steepness in initial equ recruitment calculation"<<endl;
init_int sigmaR_dendep;
!! echoinput<<sigmaR_dendep<<" # future feature: 0/1 to make realized sigmaR a function of SR curvature"<<endl;
ivector N_SRparm(1,10)
-!!N_SRparm.fill("{0,2,2,2,3,2,3,3,3,3}");
+!!N_SRparm.fill("{0,2,2,2,3,2,3,3,3,4}");
int N_SRparm2
int N_SRparm3 // with timevary links included
!!N_SRparm2=N_SRparm(SR_fxn)+3;
@@ -2026,6 +2026,13 @@
ParmLabel += "SR_RkrPower_gamma";
break;
}
+ case 10: // Bev-Holt a,b
+ {
+ ParmLabel += "SR_BH_steep_derived";
+ ParmLabel += "SR_BH_ln(alpha)";
+ ParmLabel += "SR_BH_ln(beta)";
+ break;
+ }
}
ParmLabel += "SR_sigmaR";
ParmLabel += "SR_regime";
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 93d3fe73..f6a58c32 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -49,25 +49,14 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
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.);
+ dvariable SPR = SSB_virgin_adj / Recr_virgin;
+ alpha = ((4.0 * steepness) / (1. - steepness)) / SPR ;
+ beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_curr_adj) /
(SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);
break;
}
- // Beverton-Holt with alpha beta
- /*
- 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 = (alpha*SSB_curr_adj) / (beta+SSB_curr_adj);
- break;
- }
- */
-
// SS_Label_43.3.4 constant expected recruitment
case 4: // none
{
@@ -88,8 +77,9 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 6: //Beverton-Holt constrained
{
steepness = SR_parm_work(2);
- alpha = 4.0 * steepness * Recr_virgin / (5. * steepness - 1.);
- beta = (SSB_virgin_adj * (1. - steepness)) / (5. * steepness - 1.);
+ dvariable SPR = SSB_virgin_adj / Recr_virgin;
+ alpha = ((4.0 * steepness) / (1. - steepness)) / SPR ;
+ beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
if (SSB_curr_adj > SSB_virgin_adj)
{
SSB_BH1 = SSB_virgin_adj;
@@ -167,6 +157,14 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
NewRecruits = Recr_virgin_adj * temp * mfexp(RkrTop);
break;
}
+
+ case 10: // Beverton-Holt with alpha beta
+ {
+ alpha = mfexp(SR_parm_work(3));
+ beta = mfexp(SR_parm_work(4));
+ NewRecruits = (alpha*SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
+ break;
+ }
}
RETURN_ARRAYS_DECREMENT();
return NewRecruits;
@@ -270,7 +268,7 @@ 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,
+FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
const prevariable& SSB_virgin, const prevariable& Recr_virgin, const prevariable& SPR_temp)
{
RETURN_ARRAYS_INCREMENT();
@@ -286,7 +284,7 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
dvariable srz_min;
dvariable SRZ_surv;
- steepness = SRparm2; // common usage but some different
+ steepness = SRparm(2); // 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
switch (SR_fxn)
{
@@ -296,16 +294,7 @@ 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
{
@@ -314,17 +303,32 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
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.);
+ dvariable SPR = SSB_virgin / Recr_virgin;
+ alpha = ((4.0 * steepness) / (1. - steepness)) / SPR ;
+ beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
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); //Beverton-Holt
break;
}
+ case 10: // Beverton-Holt with alpha and beta parameterization
+ {
+ alpha = mfexp(SRparm(3));
+ beta = mfexp(SRparm(4));
+ B_equil = alpha * SPR_temp - beta;
+ B_equil = posfun(B_equil, 0.0001, temp);
+ R_equil = (alpha * Recr_virgin * B_equil) / (1.0 + beta * Recr_virgin * B_equil);
+ break;
+ }
+
// SS_Label_44.1.4 constant recruitment
case 4: // constant; no bias correction
{
@@ -335,7 +339,7 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
// SS_Label_44.1.5 Hockey Stick
case 5: // hockey stick
{
- alpha = SRparm3 * Recr_virgin; // min recruitment level
+ alpha = SRparm(3) * Recr_virgin; // min recruitment level
// temp=SSB_virgin/R0*steepness; // spawners per recruit at inflection
beta = (Recr_virgin - alpha) / (steepness * SSB_virgin); // slope of recruitment on spawners below the inflection
B_equil = Join_Fxn(0.0 * SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin * steepness, SPR_temp, alpha / ((1. / SPR_temp) - beta), SPR_temp * Recr_virgin);
@@ -347,9 +351,8 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
{
SRZ_0 = log(1.0 / (SSB_virgin / Recr_virgin));
srz_min = SRZ_0 * (1.0 - steepness);
- B_equil = SSB_virgin * (1. - (log(1. / SPR_temp) - SRZ_0) / pow((srz_min - SRZ_0), (1. / SRparm3)));
- B_equil = posfun(B_equil, 0.0001, temp);
- SRZ_surv = mfexp((1. - pow((B_equil / SSB_virgin), SRparm3)) * (srz_min - SRZ_0) + SRZ_0); // survival
+ B_equil = SSB_virgin * (1. - (log(1. / SPR_temp) - SRZ_0) / pow((srz_min - SRZ_0), (1. / SRparm(3))));
+ SRZ_surv = mfexp((1. - pow((B_equil / SSB_virgin), SRparm(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
R_equil = B_equil * SRZ_surv;
break;
}
@@ -367,14 +370,14 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
// REC = (TOP/BOT)**(1.0/POWER)*SPRF0/SPR
// Power = exp(logC);
// Hupper = 1.0/(5.0 * pow(0.2,Power));
- Shepherd_c = SRparm3;
- Shepherd_c2 = pow(0.2, SRparm3);
+ Shepherd_c = SRparm(3);
+ Shepherd_c2 = pow(0.2, SRparm(3));
Hupper = 1.0 / (5.0 * Shepherd_c2);
- steepness = 0.2 + (SRparm2 - 0.2) / (0.8) * (Hupper - 0.2);
+ steepness = 0.2 + (SRparm(2) - 0.2) / (0.8) * (Hupper - 0.2);
Shep_top = 5.0 * steepness * (1.0 - Shepherd_c2) * (SPR_temp * Recr_virgin) / SSB_virgin - (1.0 - 5.0 * steepness * Shepherd_c2);
Shep_bot = 5.0 * steepness - 1.0;
Shep_top2 = posfun(Shep_top, 0.001, temp);
- R_equil = (SSB_virgin / SPR_temp) * pow((Shep_top2 / Shep_bot), (1.0 / SRparm3));
+ R_equil = (SSB_virgin / SPR_temp) * pow((Shep_top2 / Shep_bot), (1.0 / SRparm(3)));
B_equil = R_equil * SPR_temp;
break;
}
@@ -382,8 +385,8 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prev
// SS_Label_43.3.8 Ricker-power
case 9: // Ricker power 3-parameter SRR. per Punt & Cope 2017
{
- steepness = SRparm2;
- dvariable RkrPower = SRparm3;
+ steepness = SRparm(2);
+ dvariable RkrPower = SRparm(3);
temp = SSB_virgin / (SPR_temp * Recr_virgin);
dvariable RkrTop = pow(0.8, RkrPower) * log(temp) / log(5.0 * steepness);
RkrTop = posfun(RkrTop, 0.000001, CrashPen);
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index ca3cffec..a1567e8b 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1734,8 +1734,12 @@ FUNCTION void write_bigoutput()
<< 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
- << steepness << " steepness" << endl
+ << steepness << " steepness " << "#_derived_alpha_beta: "<<alpha<<" "<<beta<<endl
<< Bmsy / SSB_virgin << " Bmsy/Bzero ";
+ SS2out << " spr: "<<SPR_virgin<<endl;
+ SS2out << " der_steep: "<< alpha * SPR_virgin / (4.00 + alpha * SPR_virgin)<<endl;
+ SS2out << "der_R0: "<< 1.0 / beta * (alpha - ( 1.0 / SPR_virgin))<<endl;;
+
if (SR_fxn == 8)
{
dvariable Shepherd_c;
@@ -4743,7 +4747,7 @@ FUNCTION void SPR_profile()
Do_Equil_Calc(equ_Recr);
// SPAWN-RECR: calc equil spawn-recr in the SPR loop
SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_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)
diff --git a/StockSynthesis.code-workspace b/StockSynthesis.code-workspace
index 2bfc56e8..cdd30489 100644
--- a/StockSynthesis.code-workspace
+++ b/StockSynthesis.code-workspace
@@ -10,6 +10,7 @@
"\"*.extension\":": "\"tpl\"",
"*.htp": "c",
"ostream": "c",
+ "xlocale": "c",
"iosfwd": "c"
},
"explorer.excludeGitIgnore": true
From 94c7b88b9df1108a2caf512f04b3bae5ee61c5c8 Mon Sep 17 00:00:00 2001
From: RickMethot <richard.methot@noaa.gov>
Date: Mon, 13 Feb 2023 17:36:35 -0800
Subject: [PATCH 02/29] revise to use WHAM syntax
---
SS_benchfore.tpl | 5 ++++-
SS_popdyn.tpl | 5 ++---
SS_recruit.tpl | 43 ++++++++++++++++++++++++++++---------------
SS_write_report.tpl | 23 +++++++++++++++++------
4 files changed, 51 insertions(+), 25 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 327810cc..9c5c5e37 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -748,6 +748,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Recr_unf = mfexp(SR_parm_work(1));
Do_Equil_Calc(Recr_unf);
SSB_unf = SSB_equil;
+ report5<<" calc SSBunf "<<SSB_unf<<" recr "<<Recr_unf<<" work "<<SR_parm_work<<" parm "<<SR_parm<<endl;
SR_parm_work(N_SRparm2 + 1) = SSB_unf;
if (show_MSY == 1)
report5 << "SR_parm for benchmark: " << SR_parm_work << endl
@@ -902,7 +903,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, SSB_unf, Recr_unf, SSB_equil); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SPR_temp = SSB_equil; // based on most recent call to Do_Equil_Calc
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Bspr = Equ_SpawnRecr_Result(1);
Bspr_rec = Equ_SpawnRecr_Result(2);
@@ -1111,6 +1113,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SPR_Btgt = SSB_equil / SPR_unfished;
// 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, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
yld1(ii) = Equ_SpawnRecr_Result(1);
}
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index 70bdc8f9..c5f67be3 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -349,8 +349,8 @@ FUNCTION void get_initial_conditions()
{
alpha = mfexp(SR_parm(3));
beta = mfexp(SR_parm(4));
- steepness = alpha * SPR_virgin / (4.00 + alpha * SPR_virgin);
- Recr_virgin = 1.0 / beta * (alpha - ( 1.0 / SPR_virgin));
+ alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
+ beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
SR_parm(1) = log(Recr_virgin);
SR_parm(2) = steepness;
// warning<< SR_parm(1) << " Ln(R0) " << mfexp(SR_parm(1)) << endl
@@ -529,7 +529,6 @@ FUNCTION void get_initial_conditions()
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, SSB_virgin, Recr_virgin, SPR_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;
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index f6a58c32..3b2affb4 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -49,9 +49,10 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 3: // Beverton-Holt
{
steepness = SR_parm_work(2);
- dvariable SPR = SSB_virgin_adj / Recr_virgin;
- alpha = ((4.0 * steepness) / (1. - steepness)) / SPR ;
- beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
+// dvariable SPR = SSB_virgin_adj / Recr_virgin;
+// dvariable alpha = (4.0 * steepness) / ( SPR * (1. - steepness)) ;
+// dvariable beta = (5.0 * steepness - 1.0 )/((1.0 - steepness ) * Recr_virgin * SPR );
+// beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_curr_adj) /
(SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);
break;
@@ -77,9 +78,9 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 6: //Beverton-Holt constrained
{
steepness = SR_parm_work(2);
- dvariable SPR = SSB_virgin_adj / Recr_virgin;
- alpha = ((4.0 * steepness) / (1. - steepness)) / SPR ;
- beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
+// dvariable SPR = SSB_virgin_adj / Recr_virgin;
+// alpha = ((4.0 * steepness) / (1. - steepness)) / SPR ;
+// beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
if (SSB_curr_adj > SSB_virgin_adj)
{
SSB_BH1 = SSB_virgin_adj;
@@ -162,7 +163,7 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
{
alpha = mfexp(SR_parm_work(3));
beta = mfexp(SR_parm_work(4));
- NewRecruits = (alpha*SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
+ NewRecruits = (alpha * SSB_curr_adj) / (beta + SSB_curr_adj);
break;
}
}
@@ -310,22 +311,34 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS_Label_44.1.3 Beverton-Holt
case 3: // same as case 6
{
- dvariable SPR = SSB_virgin / Recr_virgin;
- alpha = ((4.0 * steepness) / (1. - steepness)) / SPR ;
- beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
- 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)));
+ // SR_h = 0.2 * exp(0.8*log(exp(log_SR_a) * exp(log_SPR0)));
+ // SR_R0 = log(exp(log_SR_a + log_SPR0))/(exp(log_SR_b + log_SPR0));
+
+ // SS3 previously used alternative formulation: R = A*S/(B+S)
+ // converting SS3 to align with WHAM
+ alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
+ beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
+ B_equil = (alpha * SPR_temp - 1.0) / 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); //Beverton-Holt
+ R_equil = alpha * B_equil / (1.0 + beta * B_equil);
+// report5<<" WHAM Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
+// R_equil = (4. * steepness * Recr_virgin * B_equil) / (SSB_virgin * (1. - steepness) + (5. * steepness - 1.) * B_equil); //Beverton-Holt
+// report5<<" SS3 Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
break;
}
- case 10: // Beverton-Holt with alpha and beta parameterization
+ case 10: // Beverton-Holt with alpha and beta parameterization using R = A*S/(B+S) approach
{
alpha = mfexp(SRparm(3));
beta = mfexp(SRparm(4));
- B_equil = alpha * SPR_temp - beta;
+ B_equil = (alpha * SPR_temp - 1.0) / beta;
B_equil = posfun(B_equil, 0.0001, temp);
- R_equil = (alpha * Recr_virgin * B_equil) / (1.0 + beta * Recr_virgin * B_equil);
+ R_equil = alpha * B_equil / (1.0 + beta * B_equil);
+// report5<<SPR_temp<<" Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
break;
}
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index a1567e8b..fcd322ed 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1734,13 +1734,24 @@ FUNCTION void write_bigoutput()
<< 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
- << steepness << " steepness " << "#_derived_alpha_beta: "<<alpha<<" "<<beta<<endl
+ << steepness << " steepness " << endl
<< Bmsy / SSB_virgin << " Bmsy/Bzero ";
- SS2out << " spr: "<<SPR_virgin<<endl;
- SS2out << " der_steep: "<< alpha * SPR_virgin / (4.00 + alpha * SPR_virgin)<<endl;
- SS2out << "der_R0: "<< 1.0 / beta * (alpha - ( 1.0 / SPR_virgin))<<endl;;
-
- if (SR_fxn == 8)
+ SS2out << " SPR_virgin "<<SPR_virgin<<endl;
+ if (SR_fxn == 10)
+ {
+ // SR_h = 0.2 * exp(0.8*log(exp(log_SR_a) * exp(log_SPR0)));
+ // SR_R0 = log(exp(log_SR_a + log_SPR0))/(exp(log_SR_b + log_SPR0));
+ SS2out << " derived_steepness "<< 0.2 * mfexp( 0.8 * log( alpha ) * SPR_virgin)<<" using_virgin_SPR; will be time-varying if biology is time-varying"<<endl;
+ SS2out << "derived_R0 "<< log( alpha * SPR_virgin ) / ( beta * SPR_virgin)<endl;
+ }
+ else if (SR_fxn == 3)
+ {
+ alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
+ beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
+ SS2out << " derived_alpha " << alpha << " using_virgin_SPR" << endl;
+ SS2out << "derived_beta " << beta <<endl;
+ }
+ else if (SR_fxn == 8)
{
dvariable Shepherd_c;
dvariable Shepherd_c2;
From 9e6f60fe09f41dc94fef24129a1dd676e231c4e9 Mon Sep 17 00:00:00 2001
From: RickMethot <richard.methot@noaa.gov>
Date: Mon, 13 Feb 2023 17:43:58 -0800
Subject: [PATCH 03/29] Update .gitignore
---
.gitignore | 7 +++----
1 file changed, 3 insertions(+), 4 deletions(-)
diff --git a/.gitignore b/.gitignore
index 2d39ecb0..a8becca9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,10 +7,9 @@
*.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/ss.log
Compile/ss3.log
From b66b15ff4d55a1cf5666c1a55aebe71b3acdffc6 Mon Sep 17 00:00:00 2001
From: RickMethot <richard.methot@noaa.gov>
Date: Fri, 19 May 2023 11:14:13 -0700
Subject: [PATCH 04/29] Srr=10 now matching SRR=3
---
SS_popdyn.tpl | 15 ++++++++-------
SS_readcontrol_330.tpl | 2 +-
SS_recruit.tpl | 5 ++---
SS_write_report.tpl | 10 ++++------
4 files changed, 15 insertions(+), 17 deletions(-)
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index c5f67be3..27f41507 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -347,17 +347,18 @@ FUNCTION void get_initial_conditions()
if(SR_fxn == 10) // B-H with a,b
{
+ // WHAM based on R = A*S/(1+B*S)
+ // log_SR_a = log(4 * SR_h/(exp(log_SPR0)*(1 - SR_h)));
+ // log_SR_b = log((5*SR_h - 1)/((1-SR_h)*SR_R0*exp(log_SPR0)));
+ // h = a * SPR0 / (4. + a * SPR0)
+ // R0 = 1/b * (a-1/SPR0)
+
alpha = mfexp(SR_parm(3));
beta = mfexp(SR_parm(4));
- alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
- beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
+ steepness = alpha * SPR_virgin / (4. + alpha * SPR_virgin);
+ Recr_virgin = 1. / beta * (alpha - (1. / SPR_virgin));
SR_parm(1) = log(Recr_virgin);
SR_parm(2) = steepness;
- // warning<< SR_parm(1) << " Ln(R0) " << mfexp(SR_parm(1)) << endl
- // << " alpha_beta: "<<alpha<<" "<<beta<<endl
- // << " spr: "<<SPR_virgin<<endl
- // << " der_steep: "<< alpha * SPR_virgin / (4.00 + alpha * SPR_virgin)<<endl
- // << "der_R0: "<< 1.0 / beta * (alpha - ( 1.0 / SPR_virgin))<<endl;
}
else {
Recr_virgin = mfexp(SR_parm(1));
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index 62efbf6b..b88046a0 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -1972,7 +1972,7 @@
SR_parm_timevary.initialize();
SR_env_link = 0;
SR_env_target = 0;
- //#_SR_function: 1=null; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=Survival_3Parm "<<endl;
+ //#_SR_function: 1=null; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=Survival_3Parm; 10=B-H with a,b "<<endl;
ParmLabel += "SR_LN(R0)";
switch (SR_fxn)
{
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 3b2affb4..86604f67 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -161,9 +161,10 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 10: // Beverton-Holt with alpha beta
{
+ // WHAM based on R = A*S/(1+B*S)
alpha = mfexp(SR_parm_work(3));
beta = mfexp(SR_parm_work(4));
- NewRecruits = (alpha * SSB_curr_adj) / (beta + SSB_curr_adj);
+ NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
break;
}
}
@@ -315,8 +316,6 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// 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)));
- // SR_h = 0.2 * exp(0.8*log(exp(log_SR_a) * exp(log_SPR0)));
- // SR_R0 = log(exp(log_SR_a + log_SPR0))/(exp(log_SR_b + log_SPR0));
// SS3 previously used alternative formulation: R = A*S/(B+S)
// converting SS3 to align with WHAM
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index fcd322ed..22e1df1b 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1739,17 +1739,15 @@ FUNCTION void write_bigoutput()
SS2out << " SPR_virgin "<<SPR_virgin<<endl;
if (SR_fxn == 10)
{
- // SR_h = 0.2 * exp(0.8*log(exp(log_SR_a) * exp(log_SPR0)));
- // SR_R0 = log(exp(log_SR_a + log_SPR0))/(exp(log_SR_b + log_SPR0));
- SS2out << " derived_steepness "<< 0.2 * mfexp( 0.8 * log( alpha ) * SPR_virgin)<<" using_virgin_SPR; will be time-varying if biology is time-varying"<<endl;
- SS2out << "derived_R0 "<< log( alpha * SPR_virgin ) / ( beta * SPR_virgin)<endl;
+ SS2out << " Ln_alpha_parameter: " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
+ SS2out << " Ln_beta_parameter: " << SR_parm(4) << " beta " << mfexp(SR_parm(4)) << endl;
}
else if (SR_fxn == 3)
{
alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
- SS2out << " derived_alpha " << alpha << " using_virgin_SPR" << endl;
- SS2out << "derived_beta " << beta <<endl;
+ SS2out << " Ln_alpha_derived: " << log(alpha) << " alpha " << alpha << endl;
+ SS2out << " Ln_beta_derived: " << log(beta) << " beta " << beta;
}
else if (SR_fxn == 8)
{
From 449e87b1f63a29e52091e7c97544cdc0f5b1be2e Mon Sep 17 00:00:00 2001
From: RickMethot <richard.methot@noaa.gov>
Date: Fri, 19 May 2023 17:06:15 -0700
Subject: [PATCH 05/29] adding reporting for timevary growth
---
SS_benchfore.tpl | 8 ++++++++
SS_popdyn.tpl | 8 ++++++++
SS_write_report.tpl | 14 ++++++++------
StockSynthesis.code-workspace | 3 ++-
4 files changed, 26 insertions(+), 7 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 9c5c5e37..42b0b332 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -3646,6 +3646,14 @@ FUNCTION void Get_Forecast()
Smry_Table(y, 11) = SSB_equil;
Smry_Table(y, 13) = GenTime;
+ if( SR_fxn == 10 )
+ {
+ temp = SSB_equil / equ_Recr; // current year's SPB/R with current biology at age
+ alpha = mfexp(SR_parm_work(3));
+ beta = mfexp(SR_parm_work(4));
+ SR_parm_byyr(y, 2) = alpha * temp / (4. + alpha * temp); // implied steepness
+ SR_parm_byyr(y, 1) = log( 1. / beta * (alpha - (1. / temp))); // implied ln_R0
+ }
Fishon = 1;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
if (STD_Yr_Reverse_Ofish(y) > 0)
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index 27f41507..b2dd8299 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -1782,6 +1782,14 @@ FUNCTION void get_time_series()
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation with current year's biology
Smry_Table(y, 11) = SSB_equil;
Smry_Table(y, 13) = GenTime;
+ if( SR_fxn == 10 )
+ {
+ temp = SSB_equil / Recr_virgin; // current year's SPB/R with current biology at age
+ alpha = mfexp(SR_parm_work(3));
+ beta = mfexp(SR_parm_work(4));
+ SR_parm_byyr(y, 2) = alpha * temp / (4. + alpha * temp); // implied steepness
+ SR_parm_byyr(y, 1) = log( 1. / beta * (alpha - (1. / temp))); // implied ln_R0
+ }
Fishon = 1;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation with current year's biology and F
if (STD_Yr_Reverse_Ofish(y) > 0)
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index 22e1df1b..bda7f9ac 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1740,7 +1740,7 @@ FUNCTION void write_bigoutput()
if (SR_fxn == 10)
{
SS2out << " Ln_alpha_parameter: " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
- SS2out << " Ln_beta_parameter: " << SR_parm(4) << " beta " << mfexp(SR_parm(4)) << endl;
+ SS2out << " Ln_beta_parameter: " << SR_parm(4) << " beta " << mfexp(SR_parm(4));
}
else if (SR_fxn == 3)
{
@@ -1764,9 +1764,7 @@ FUNCTION void write_bigoutput()
{
SS2out << " Ricker_Power: " << SR_parm(3);
}
-
- SS2out << endl;
- SS2out << sigmaR << " sigmaR" << endl;
+ SS2out << endl << 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
@@ -1799,7 +1797,7 @@ FUNCTION void write_bigoutput()
}
SS2out << endl;
- SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev" << endl;
+ SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SPR0 h R0" << endl;
SS2out << "S/Rcurve " << SSB_virgin << " " << Recr_virgin << endl;
y = styr - 2;
SS2out << "Virg " << SSB_yr(y) << " " << exp_rec(y) << " - " << 0.0 << " Virg " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " 0.0 " << endl;
@@ -1850,7 +1848,11 @@ FUNCTION void write_bigoutput()
{
SS2out << " _ _ Fixed";
}
- SS2out << endl;
+ temp = Smry_Table(y,11) / Recr_virgin;
+ alpha = mfexp(SR_parm_byyr(y,3));
+ beta = mfexp(SR_parm_byyr(y,4));
+ SS2out << " " << temp << " " << alpha * temp / (4. + alpha * temp) << " " << 1. / beta * (alpha - (1. / temp));
+ SS2out << SR_parm_byyr(y)(1,4) << endl;
}
// REPORT_KEYWORD SPAWN_RECR_CURVE
diff --git a/StockSynthesis.code-workspace b/StockSynthesis.code-workspace
index cdd30489..14159ac5 100644
--- a/StockSynthesis.code-workspace
+++ b/StockSynthesis.code-workspace
@@ -11,7 +11,8 @@
"*.htp": "c",
"ostream": "c",
"xlocale": "c",
- "iosfwd": "c"
+ "iosfwd": "c",
+ "cmath": "c"
},
"explorer.excludeGitIgnore": true
}
From 43afa5139b7ca3279d5a7193da4c08d3a74547cc Mon Sep 17 00:00:00 2001
From: RickMethot <richard.methot@noaa.gov>
Date: Wed, 31 May 2023 09:40:07 -0700
Subject: [PATCH 06/29] WIP
---
SS_benchfore.tpl | 19 ++++++++++++-------
SS_recruit.tpl | 5 ++++-
2 files changed, 16 insertions(+), 8 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 42b0b332..40f9bb89 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -748,18 +748,22 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Recr_unf = mfexp(SR_parm_work(1));
Do_Equil_Calc(Recr_unf);
SSB_unf = SSB_equil;
- report5<<" calc SSBunf "<<SSB_unf<<" recr "<<Recr_unf<<" work "<<SR_parm_work<<" parm "<<SR_parm<<endl;
- SR_parm_work(N_SRparm2 + 1) = SSB_unf;
+ SPR_unfished = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
if (show_MSY == 1)
+ {
report5 << "SR_parm for benchmark: " << SR_parm_work << endl
- << "for years: " << Bmark_Yr(9) << " " << Bmark_Yr(10) << " SSB_virgin was: " << SSB_virgin << endl;
- if (show_MSY == 1)
- report5 << "Repro_output_by_age_for_morph_1: " << fec(1) << endl;
+ << "mean from years: " << Bmark_Yr(9) << " " << Bmark_Yr(10) << endl;
+ SPR_virgin = SSB_virgin / Recr_virgin;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Virgin SPR0, SSB, R: " << SPR_virgin << " " << Equ_SpawnRecr_Result << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_unfished); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Benchmark SPR0, SSB, R: " << SPR_unfished << " " << Equ_SpawnRecr_Result << endl;
+ }
+ SR_parm_work(N_SRparm2 + 1) = SSB_unf;
Mgmt_quant(1) = SSB_unf;
Mgmt_quant(2) = totbio;
Mgmt_quant(3) = smrybio;
Mgmt_quant(4) = Recr_unf;
-
// find Fspr SS_Label_710
{
if (show_MSY == 1)
@@ -905,6 +909,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
// SPAWN-RECR: calc equil spawn-recr in YPR; need to make this area-specific
SPR_temp = SSB_equil; // based on most recent call to Do_Equil_Calc
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5<<SPR_temp<<" " << Equ_SpawnRecr_Result<<endl;
Bspr = Equ_SpawnRecr_Result(1);
Bspr_rec = Equ_SpawnRecr_Result(2);
@@ -1113,8 +1118,8 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SPR_Btgt = SSB_equil / SPR_unfished;
// 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, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5<<SPR_temp<<" " << Equ_SpawnRecr_Result<<endl;
yld1(ii) = Equ_SpawnRecr_Result(1);
}
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 86604f67..c24509c0 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -319,12 +319,15 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS3 previously used alternative formulation: R = A*S/(B+S)
// converting SS3 to align with WHAM
+ SPR_virgin = SSB_virgin / Recr_virgin;
alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
+ report5<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_virgin<<" SPR "<<SPR_virgin<<" steep "<<steepness<<endl;
B_equil = (alpha * SPR_temp - 1.0) / beta;
+// report5<<" spr "<<SPR_temp<<" Beq "<<B_equil<<endl;
B_equil = posfun(B_equil, 0.0001, temp);
R_equil = alpha * B_equil / (1.0 + beta * B_equil);
-// report5<<" WHAM Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
+// report5<<" Beq "<<B_equil<<" Req "<<R_equil<<endl;
// R_equil = (4. * steepness * Recr_virgin * B_equil) / (SSB_virgin * (1. - steepness) + (5. * steepness - 1.) * B_equil); //Beverton-Holt
// report5<<" SS3 Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
break;
From d6316c4b583678afb9be3c0fc521e89eb50611b3 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Wed, 22 May 2024 16:11:03 -0700
Subject: [PATCH 07/29] working version with a,b
---
SS_recruit.tpl | 41 +++++++++++++++++++----------------------
SS_write_report.tpl | 8 ++++----
2 files changed, 23 insertions(+), 26 deletions(-)
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index c24509c0..5553dd26 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -49,15 +49,19 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 3: // Beverton-Holt
{
steepness = SR_parm_work(2);
-// dvariable SPR = SSB_virgin_adj / Recr_virgin;
-// dvariable alpha = (4.0 * steepness) / ( SPR * (1. - steepness)) ;
-// dvariable beta = (5.0 * steepness - 1.0 )/((1.0 - steepness ) * Recr_virgin * SPR );
-// beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_curr_adj) /
(SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);
break;
}
+ case 10: // Beverton-Holt with alpha beta per WHAM: R = A*S/(1+B*S)
+ {
+ alpha = mfexp(SR_parm_work(3));
+ beta = mfexp(SR_parm_work(4));
+ NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
+ break;
+ }
+
// SS_Label_43.3.4 constant expected recruitment
case 4: // none
{
@@ -159,14 +163,6 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
break;
}
- case 10: // Beverton-Holt with alpha beta
- {
- // WHAM based on R = A*S/(1+B*S)
- alpha = mfexp(SR_parm_work(3));
- beta = mfexp(SR_parm_work(4));
- NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
- break;
- }
}
RETURN_ARRAYS_DECREMENT();
return NewRecruits;
@@ -322,18 +318,19 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
SPR_virgin = SSB_virgin / Recr_virgin;
alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
- report5<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_virgin<<" SPR "<<SPR_virgin<<" steep "<<steepness<<endl;
+
+// report5 <<" SSB_unf "<<SSB_virgin<<" SPR_unf "<<SPR_virgin<<" steep: "<<steepness<<" R0: "<<Recr_virgin << endl;
+// report5 <<" derive_alpha "<<alpha<<" derive_beta "<<beta << endl;
+// report5 << " deriv_h: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " derive_R0: " << 1. / beta * (alpha - (1. / SPR_virgin))<<endl;
B_equil = (alpha * SPR_temp - 1.0) / beta;
-// report5<<" spr "<<SPR_temp<<" Beq "<<B_equil<<endl;
B_equil = posfun(B_equil, 0.0001, temp);
R_equil = alpha * B_equil / (1.0 + beta * B_equil);
-// report5<<" Beq "<<B_equil<<" Req "<<R_equil<<endl;
-// R_equil = (4. * steepness * Recr_virgin * B_equil) / (SSB_virgin * (1. - steepness) + (5. * steepness - 1.) * B_equil); //Beverton-Holt
-// report5<<" SS3 Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
+// report5 << "SPR_input: " << SPR_temp << " B_equil: " << B_equil << " R_equil: "<<R_equil << endl<<endl;
+
break;
}
- case 10: // Beverton-Holt with alpha and beta parameterization using R = A*S/(B+S) approach
+ case 10: // Beverton-Holt with alpha and beta parameterization using R = A*S/(1+B*S) approach; same as WHAM
{
alpha = mfexp(SRparm(3));
beta = mfexp(SRparm(4));
@@ -354,11 +351,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS_Label_44.1.5 Hockey Stick
case 5: // hockey stick
{
- alpha = SRparm(3) * Recr_virgin; // min recruitment level
+ dvariable hockey_min = SRparm(3) * Recr_virgin; // min recruitment level
// temp=SSB_virgin/R0*steepness; // spawners per recruit at inflection
- beta = (Recr_virgin - alpha) / (steepness * SSB_virgin); // slope of recruitment on spawners below the inflection
- B_equil = Join_Fxn(0.0 * SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin * steepness, SPR_temp, alpha / ((1. / SPR_temp) - beta), SPR_temp * Recr_virgin);
- R_equil = Join_Fxn(0.0 * SSB_virgin, SSB_virgin, SSB_virgin * steepness, B_equil, alpha + beta * B_equil, Recr_virgin);
+ dvariable hockey_slope = (Recr_virgin - hockey_min) / (steepness * SSB_virgin); // slope of recruitment on spawners below the inflection
+ B_equil = Join_Fxn(0.0 * SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin * steepness, SPR_temp, hockey_min / ((1. / SPR_temp) - hockey_slope), SPR_temp * Recr_virgin);
+ R_equil = Join_Fxn(0.0 * SSB_virgin, SSB_virgin, SSB_virgin * steepness, B_equil, hockey_min + hockey_slope * B_equil, Recr_virgin);
break;
}
// SS_Label_44.1.7 3 parameter survival based
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index bda7f9ac..f1b33812 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1745,7 +1745,7 @@ FUNCTION void write_bigoutput()
else if (SR_fxn == 3)
{
alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
- beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
+ beta = (5.0 * steepness - 1.0) / ((1. - steepness) * SSB_virgin);
SS2out << " Ln_alpha_derived: " << log(alpha) << " alpha " << alpha << endl;
SS2out << " Ln_beta_derived: " << log(beta) << " beta " << beta;
}
@@ -1797,7 +1797,7 @@ FUNCTION void write_bigoutput()
}
SS2out << endl;
- SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SPR0 h R0" << endl;
+ SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SSBe/R0 h R0" << endl;
SS2out << "S/Rcurve " << SSB_virgin << " " << Recr_virgin << endl;
y = styr - 2;
SS2out << "Virg " << SSB_yr(y) << " " << exp_rec(y) << " - " << 0.0 << " Virg " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " 0.0 " << endl;
@@ -1848,10 +1848,10 @@ FUNCTION void write_bigoutput()
{
SS2out << " _ _ Fixed";
}
- temp = Smry_Table(y,11) / Recr_virgin;
+ dvariable SPR = Smry_Table(y, 11) / Recr_virgin;
alpha = mfexp(SR_parm_byyr(y,3));
beta = mfexp(SR_parm_byyr(y,4));
- SS2out << " " << temp << " " << alpha * temp / (4. + alpha * temp) << " " << 1. / beta * (alpha - (1. / temp));
+ SS2out << " " << SPR << " " << alpha * SPR / (4. + alpha * SPR) << " " << 1. / beta * (alpha - (1. / SPR));
SS2out << SR_parm_byyr(y)(1,4) << endl;
}
From 8558f5e207d35d1dc7c535ae404cc4cb3b981181 Mon Sep 17 00:00:00 2001
From: "Richard.Methot" <1365818351113302@mil>
Date: Tue, 28 May 2024 14:39:34 -0700
Subject: [PATCH 08/29] Update .gitignore
---
.gitignore | 2 ++
1 file changed, 2 insertions(+)
diff --git a/.gitignore b/.gitignore
index a8becca9..fac35ef1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,6 +10,8 @@
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
From 9886eb6edaaef28b47f0f66f3bf352cdfc2e3213 Mon Sep 17 00:00:00 2001
From: Rick-Methot-NOAA <richard.methot@noaa.gov>
Date: Wed, 12 Jun 2024 14:50:29 -0700
Subject: [PATCH 09/29] augment spawn_recr report and use to test time-vary
SR_parm approach
---
SS_benchfore.tpl | 12 ++--
SS_popdyn.tpl | 18 ++++--
SS_readcontrol_330.tpl | 9 +++
SS_recruit.tpl | 13 ++--
SS_write_report.tpl | 139 ++++++++++++++++++++++++++++++++++-------
5 files changed, 156 insertions(+), 35 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 40f9bb89..c76438e0 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -753,7 +753,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
{
report5 << "SR_parm for benchmark: " << SR_parm_work << endl
<< "mean from years: " << Bmark_Yr(9) << " " << Bmark_Yr(10) << endl;
- SPR_virgin = SSB_virgin / Recr_virgin;
+ // SPR_virgin = SSB_virgin / Recr_virgin; // already defined
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
report5 << " Virgin SPR0, SSB, R: " << SPR_virgin << " " << Equ_SpawnRecr_Result << endl;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_unfished); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
@@ -1906,15 +1906,15 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
else if (show_MSY == 2) // do brief output
{
- SS2out << SPR_actual / 100. << " " << SPR_Fmult << " " << Mgmt_quant(10) << " " << YPR_spr_dead / Vbio1_spr << " " << Bspr_rec << " "
+ report5 << SPR_actual / 100. << " " << SPR_Fmult << " " << Mgmt_quant(10) << " " << YPR_spr_dead / Vbio1_spr << " " << Bspr_rec << " "
<< Bspr << " " << YPR_spr_dead * Bspr_rec << " " << YPR_spr_ret * Bspr_rec
<< " " << Vbio1_spr * Bspr_rec << " # ";
- SS2out << SPR_Btgt << " " << Btgt / SSB_unf << " " << Btgt_Fmult << " " << Mgmt_quant(7) << " " << YPR_Btgt_dead / Vbio1_Btgt << " " << Btgt_Rec << " "
+ report5 << SPR_Btgt << " " << Btgt / SSB_unf << " " << Btgt_Fmult << " " << Mgmt_quant(7) << " " << YPR_Btgt_dead / Vbio1_Btgt << " " << Btgt_Rec << " "
<< Btgt << " " << YPR_Btgt_dead * Btgt_Rec << " " << YPR_Btgt_ret * Btgt_Rec
<< " " << Vbio1_Btgt * Btgt_Rec << " # ";
- SS2out << MSY_SPR << " " << Bmsy / SSB_unf << " " << MSY_Fmult << " " << Mgmt_quant(14) << " " << MSY / (Vbio1_MSY * Recr_msy) << " " << Recr_msy << " "
+ report5 << MSY_SPR << " " << Bmsy / SSB_unf << " " << MSY_Fmult << " " << Mgmt_quant(14) << " " << MSY / (Vbio1_MSY * Recr_msy) << " " << Recr_msy << " "
<< Bmsy << " " << MSY << " " << YPR_msy_dead * Recr_msy << " " << YPR_msy_ret * Recr_msy
<< " " << Vbio1_MSY * Recr_msy << " # " << endl;
}
@@ -2576,6 +2576,7 @@ FUNCTION void Get_Forecast()
}
}
// SPAWN-RECR: get recruitment in forecast; needs to be area-specific
+ // SR_fxn
if (SR_parm_timevary(1) == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
@@ -3212,7 +3213,8 @@ FUNCTION void Get_Forecast()
}
}
// 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
+ // SPAWN-RECR: calc recruitment in time series; need to make this area-specific
+ // SR_fxn
if (SR_parm_timevary(1) == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index b2dd8299..3891000a 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -357,6 +357,7 @@ FUNCTION void get_initial_conditions()
beta = mfexp(SR_parm(4));
steepness = alpha * SPR_virgin / (4. + alpha * SPR_virgin);
Recr_virgin = 1. / beta * (alpha - (1. / SPR_virgin));
+// warning << " before AB_calcs " << "parm " << SR_parm(1) << " calc " << log(Recr_virgin) << endl;
SR_parm(1) = log(Recr_virgin);
SR_parm(2) = steepness;
}
@@ -370,6 +371,8 @@ FUNCTION void get_initial_conditions()
SR_parm_virg(i) = SR_parm(i);
SR_parm_work(i) = SR_parm(i);
}
+// if (SR_fxn == 3) warning << "tester_A: " << SR_parm_work(1) << " base: " << SR_parm(1) << endl;
+// if (SR_fxn == 10) warning << "tester_A: " << SR_parm_work(4) << " base: " << SR_parm(4) << endl;
equ_Recr = Recr_virgin;
exp_rec(eq_yr, 1) = Recr_virgin; // expected Recr from s-r parms
exp_rec(eq_yr, 2) = Recr_virgin;
@@ -476,6 +479,9 @@ FUNCTION void get_initial_conditions()
else
{
SR_parm_work(f) = parm_timevary(SR_parm_timevary(f), eq_yr);
+// warning << "tester_B: " << SR_parm_work(f) << " timevary " << " base " << SR_parm(f) <<endl;
+// warning << parm_timevary(2) << endl;
+
}
SR_parm_byyr(eq_yr, f) = SR_parm_work(f);
}
@@ -529,7 +535,7 @@ FUNCTION void get_initial_conditions()
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, SSB_virgin, Recr_virgin, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_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;
@@ -725,6 +731,7 @@ FUNCTION void get_time_series()
else
{
SR_parm_work(f) = parm_timevary(SR_parm_timevary(f), y);
+// warning << "tester_C: " << SR_parm_work(f) << " timevary_year " << endl;
}
SR_parm_byyr(y, f) = SR_parm_work(f);
}
@@ -1017,6 +1024,7 @@ FUNCTION void get_time_series()
// 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
+ // SR_Fxn relevant keyword
if (SR_parm_timevary(1) == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
@@ -1025,6 +1033,7 @@ FUNCTION void get_time_series()
else
{
R0_use = mfexp(SR_parm_work(1));
+ warning << y << " set R0use to SRparm_work " << R0_use << " vir: " << Recr_virgin << endl;
equ_Recr = R0_use;
Fishon = 0;
eq_yr = y;
@@ -1475,14 +1484,15 @@ FUNCTION void get_time_series()
}
// 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-specific
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ // SR_fxn
+ if (SR_update_SPR0 == 0) // SR parms are not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
}
- else
+ else // update SSB_use and R0_use where will update SPR0 inside the Spawn_recr fxn
{
- R0_use = mfexp(SR_parm_work(1));
+ R0_use = mfexp(SR_parm_work(1)); // check to be sure this works when R0 is derived from B-H with alpha, beta parameters
equ_Recr = R0_use;
Fishon = 0;
eq_yr = y;
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index b88046a0..edd7a3a7 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -1656,6 +1656,7 @@
!!// SS_Label_Info_4.5.4 #Set up time-varying parameters for MG parms
int timevary_used;
+ int timevary_MG_firstyr;
int timevary_parm_cnt_MG;
int timevary_parm_start_MG;
@@ -1704,6 +1705,7 @@
timevary_parm_start_MG = 0;
timevary_parm_cnt_MG = 0;
timevary_used = 0;
+ timevary_MG_firstyr = YrMax;
MGparm_timevary.initialize();
ivector block_design_null(1, 1);
block_design_null.initialize();
@@ -1804,6 +1806,7 @@
{
MG_active(f) = 1;
timevary_MG(y, 0) = 1; // tracks active status for all MG types
+ if(timevary_MG_firstyr == YrMax) timevary_MG_firstyr = y; // save for reporting in MSY and spawn_recruit output
}
}
// timevary growth or maturity and Maunder M refers to that maturity
@@ -1961,6 +1964,7 @@
int timevary_parm_cnt_SR;
ivector timevary_SRparm(styr-3,YrMax+1);
ivector SR_parm_timevary(1,N_SRparm2);
+ int SR_update_SPR0 // 0/1 flag to control updating of SPR0 for timevary biology
LOCAL_CALCS
// clang-format on
@@ -1972,6 +1976,8 @@
SR_parm_timevary.initialize();
SR_env_link = 0;
SR_env_target = 0;
+ SR_update_SPR0 = 0;
+
//#_SR_function: 1=null; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=Survival_3Parm; 10=B-H with a,b "<<endl;
ParmLabel += "SR_LN(R0)";
switch (SR_fxn)
@@ -2118,6 +2124,9 @@
for (y = styr - 3; y <= YrMax + 1; y++) {
timevary_SRparm(y) = timevary_pass(y);
} // year vector for this category og MGparm
+ // special consideration for impact on spawner-recruitment SPR0 calculations
+ if (SR_fxn == 10 && (SR_parm_timevary(3) > 0 || SR_parm_timevary(4) > 0) ) {SR_update_SPR0 = 1;} // alpha or beta is time-varying
+ else if ((SR_parm_timevary(1) > 0 || SR_parm_timevary(2) > 0) ) {SR_update_SPR0 = 1;} // R0 or steepness is time-varying
}
}
N_SRparm3 = N_SRparm2;
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 5553dd26..2065b277 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -22,6 +22,7 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
dvariable SRZ_0;
dvariable srz_min;
dvariable SRZ_surv;
+// warning << y << " Tester_R0 " << Recr_virgin_adj << " SSB0 " << SSB_virgin_adj << " 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
@@ -56,8 +57,8 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 10: // Beverton-Holt with alpha beta per WHAM: R = A*S/(1+B*S)
{
- alpha = mfexp(SR_parm_work(3));
- beta = mfexp(SR_parm_work(4));
+ dvariable alpha = mfexp(SR_parm_work(3));
+ dvariable beta = mfexp(SR_parm_work(4));
NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
break;
}
@@ -315,10 +316,12 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS3 previously used alternative formulation: R = A*S/(B+S)
// converting SS3 to align with WHAM
- SPR_virgin = SSB_virgin / Recr_virgin;
+ // SPR_virgin = SSB_virgin / Recr_virgin; // this is already defined
alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
+ // " h " << steepness << " derive " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " " << endl;
+ // " R0 " << Recr_virgin << " derive " << 1. / beta * (alpha - 1./SPR_virgin) << endl;
// report5 <<" SSB_unf "<<SSB_virgin<<" SPR_unf "<<SPR_virgin<<" steep: "<<steepness<<" R0: "<<Recr_virgin << endl;
// report5 <<" derive_alpha "<<alpha<<" derive_beta "<<beta << endl;
// report5 << " deriv_h: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " derive_R0: " << 1. / beta * (alpha - (1. / SPR_virgin))<<endl;
@@ -332,8 +335,8 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
case 10: // Beverton-Holt with alpha and beta parameterization using R = A*S/(1+B*S) approach; same as WHAM
{
- alpha = mfexp(SRparm(3));
- beta = mfexp(SRparm(4));
+ dvariable alpha = mfexp(SRparm(3));
+ dvariable beta = mfexp(SRparm(4));
B_equil = (alpha * SPR_temp - 1.0) / beta;
B_equil = posfun(B_equil, 0.0001, temp);
R_equil = alpha * B_equil / (1.0 + beta * B_equil);
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index f1b33812..c5b13169 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1730,26 +1730,14 @@ FUNCTION void write_bigoutput()
var /= (n_rmse(1) + 1.0e-09);
dvariable steepness = SR_parm(2);
+
SS2out << endl
<< pick_report_name(19);
SS2out << " Function: " << SR_fxn << " RecDev_method: " << do_recdev << " sum_recdev: " << sum_recdev << endl
<< SR_parm(1) << " Ln(R0) " << mfexp(SR_parm(1)) << endl
- << steepness << " steepness " << endl
+ << steepness << " steepness" << endl
<< Bmsy / SSB_virgin << " Bmsy/Bzero ";
- SS2out << " SPR_virgin "<<SPR_virgin<<endl;
- if (SR_fxn == 10)
- {
- SS2out << " Ln_alpha_parameter: " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
- SS2out << " Ln_beta_parameter: " << SR_parm(4) << " beta " << mfexp(SR_parm(4));
- }
- else if (SR_fxn == 3)
- {
- alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
- beta = (5.0 * steepness - 1.0) / ((1. - steepness) * SSB_virgin);
- SS2out << " Ln_alpha_derived: " << log(alpha) << " alpha " << alpha << endl;
- SS2out << " Ln_beta_derived: " << log(beta) << " beta " << beta;
- }
- else if (SR_fxn == 8)
+ if (SR_fxn == 8)
{
dvariable Shepherd_c;
dvariable Shepherd_c2;
@@ -1764,7 +1752,9 @@ FUNCTION void write_bigoutput()
{
SS2out << " Ricker_Power: " << SR_parm(3);
}
- SS2out << endl << sigmaR << " sigmaR" << endl;
+
+ 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
@@ -1797,7 +1787,114 @@ FUNCTION void write_bigoutput()
}
SS2out << endl;
- SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SSBe/R0 h R0" << endl;
+ SS2out << endl << "#New_Expanded_Spawn_Recr_report" << endl << pick_report_name(19) << endl;
+ SS2out << "SR_Function: " << SR_fxn << endl;
+ SS2out << "#" << endl << "parm parm_label value phase" << endl;
+ for (int j = 1; j <=N_SRparm2; j++)
+ {
+ SS2out << j << " " << ParmLabel(firstSRparm + j) << " " << SR_parm(j) << " " << SR_parm_PH(j);
+ if (SR_parm_timevary(j) > 0 && j <= 4 ) // timevary SRparm exists
+ {SS2out << " #_is_time_vary,_so_SRR_updates_base_SPR_annually";}
+ if (j == (N_SRparm2 - 1) && SR_parm_timevary(j) > 0) // timevary regime exists
+ {SS2out << " #_Persistent_deviations_from_SRR_(e.g.regimes)_exist";}
+ SS2out << endl;
+ }
+
+ SS2out << "# " <<endl;
+ SS2out << "sigmaR: "<< sigmaR << endl;
+ SS2out << "SPR_virgin: " << SPR_virgin << " #_uses_biology_from_start_year: " << styr <<endl;
+
+ switch (SR_fxn)
+ {
+ case 3: // Beverton-Holt with steepness
+ {
+ alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
+ beta = (5.0 * steepness - 1.0) / ((1. - steepness) * SSB_virgin);
+ SS2out << "Ln(R0): " << SR_parm(1) << endl << "R0: " << mfexp(SR_parm(1)) << endl;
+ SS2out << "steepness: " << steepness << endl;
+ SS2out << "Ln(alpha)_derived: " << log(alpha) << " alpha " << alpha << endl;
+ SS2out << "Ln(beta)_derived: " << log(beta) << " beta " << beta;
+ break;
+ }
+ case 10: // Beverton-Holt with alpha, beta
+ {
+ SS2out << "Ln(alpha): " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
+ SS2out << "Ln(beta): " << SR_parm(4) << " beta " << mfexp(SR_parm(4)) << endl;
+ SS2out << "steepness_derived: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " " << SR_parm(3) << endl; // steepness virgin
+ SS2out << "R0_derived: " << 1. / beta * (alpha - (1. / SPR_virgin)) << " " << SR_parm(1) << endl; // virgin R0
+ break;
+ }
+ case 8:
+ {
+ dvariable Shepherd_c;
+ dvariable Shepherd_c2;
+ dvariable Hupper;
+ Shepherd_c = SR_parm(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);
+ SS2out << " Shepherd_c: " << Shepherd_c << " steepness_limit: " << Hupper << " Adjusted_steepness: " << temp << endl;
+ break;
+ }
+ case 9:
+ {
+ SS2out << " Ricker_Parm1: " << SR_parm(1) << endl;
+ SS2out << " Ricker_Parm2: " << SR_parm(2) << endl;
+ SS2out << " Ricker_Power: " << SR_parm(3) << endl;
+ break;
+ }
+ default:
+ {
+ SS2out << "default output needed " << endl;
+ break;
+ }
+ }
+
+ SS2out << endl << "#" << endl << "Quantities for MSY and other benchmark calculations " << endl << "Benchmark_years: 1_beg_bio 2_end_bio 3_beg_selex 4_end_selex 5_beg_relF 6_end_relF 7_beg_recr_dist 8_end_recr_dist 9_beg_SRparm 10_end_SRparm" << endl;
+ SS2out << "Benchmark_years: " << Bmark_Yr << endl;
+ if( timevary_MG_firstyr == YrMax)
+ { SS2out << "No_timevary_biology " << endl; }
+ else
+ { SS2out << "First_year_with_timevary_MG: " << timevary_MG_firstyr << endl; }
+ for ( int k = 1; k <=9; k+=2)
+ { if (Bmark_Yr(k+1) > Bmark_Yr(k))
+ {SS2out << "#_range_of_years_is_averaged,_so_reduces_standard_error_of_result;_do_this_only_when_timevarying_makes_necessary: " << k << " "<< k+1 << endl;}
+ }
+ SS2out << "SPR_unfished_benchmark: " << Mgmt_quant(1) / Mgmt_quant(4) << " #_based_on_averaging_biology_over_benchmark_year_range " << endl;
+ SS2out << "Bmsy/Bzero: "<< Bmsy / SSB_virgin << " # using styr bio for Bzero" << endl;
+ SS2out << "Bmsy/Bunf: "<< Bmsy / Mgmt_quant(1) << " # using MSY's averaged bio for Bunf" << endl;
+
+ SS2out << "#" << endl << "RecDev_method: " << do_recdev << endl << "sum_recdev: " << sum_recdev << endl << "recr_logL: " << recr_like << endl;
+ SS2out << recdev_start << " " << recdev_end << " main_recdev:start_end" << endl
+ << recdev_adj(1) << " " << recdev_adj(2, 5) << " breakpoints_for_bias_adjustment_ramp " << endl;
+
+ temp = sigmaR * sigmaR; // sigmaR^2
+ SS2out << "ERA N RMSE RMSE^2/sigmaR^2 mean_BiasAdj est_rho Durbin-Watson" << endl;
+ SS2out << "main " << n_rmse(1) << " " << rmse(1) << " " << square(rmse(1)) / temp << " " << rmse(2) << " " << cross / var << " " << Durbin;
+ if (wrote_bigreport == 0) // first time writing bigreport
+ {
+ if (rmse(1) < 0.5 * sigmaR && rmse(2) > (0.01 + 2.0 * square(rmse(1)) / temp))
+ {
+ warnstream << "Main recdev biasadj is >2 times ratio of rmse to sigmaR";
+ SS2out << " # " << warnstream.str() ;
+ write_message (WARN, 0);
+ }
+ }
+ SS2out << endl;
+
+ SS2out << "early " << n_rmse(3) << " " << rmse(3) << " " << square(rmse(3)) / temp << " " << rmse(4);
+ if (wrote_bigreport == 0) // first time writing bigreport
+ {
+ if (rmse(3) < 0.5 * sigmaR && rmse(4) > (0.01 + 2.0 * square(rmse(3)) / temp))
+ {
+ warnstream << "Early recdev biasadj is >2 times ratio of rmse to sigmaR";
+ SS2out << " # " << warnstream.str();
+ write_message (WARN, 0);
+ }
+ }
+ SS2out << endl << "#" << endl << "Initial_equilibrium: " << init_equ_steepness << " # 0/1_to_use_spawner-recruitment_in_initial_equ_recruitment_calculation" << endl << "#" << endl;
+
+ SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SPR0_curr h_curr R0_curr P1 P2 P3 P4" << endl;
SS2out << "S/Rcurve " << SSB_virgin << " " << Recr_virgin << endl;
y = styr - 2;
SS2out << "Virg " << SSB_yr(y) << " " << exp_rec(y) << " - " << 0.0 << " Virg " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " 0.0 " << endl;
@@ -1848,10 +1945,10 @@ FUNCTION void write_bigoutput()
{
SS2out << " _ _ Fixed";
}
- dvariable SPR = Smry_Table(y, 11) / Recr_virgin;
- alpha = mfexp(SR_parm_byyr(y,3));
- beta = mfexp(SR_parm_byyr(y,4));
- SS2out << " " << SPR << " " << alpha * SPR / (4. + alpha * SPR) << " " << 1. / beta * (alpha - (1. / SPR));
+ dvariable SPR_curr = Smry_Table(y, 11) / Recr_virgin;
+ SS2out << " " << SPR_curr << " ";
+ SS2out << alpha * SPR_curr / (4. + alpha * SPR_curr) << " "; // steepness with current SPR
+ SS2out << 1. / beta * (alpha - (1. / SPR_curr)) << " "; // R0 with current SPR
SS2out << SR_parm_byyr(y)(1,4) << endl;
}
From a9b353d42a74af9011c1b4011a9ae9f0b646a96d Mon Sep 17 00:00:00 2001
From: e-perl-NOAA <elizabeth.gugliotti@noaa.gov>
Date: Thu, 13 Jun 2024 08:46:20 -0400
Subject: [PATCH 10/29] try downgrading qemu
---
.github/workflows/build-ss3.yml | 2 ++
1 file changed, 2 insertions(+)
diff --git a/.github/workflows/build-ss3.yml b/.github/workflows/build-ss3.yml
index 108c745e..1d944ce7 100644
--- a/.github/workflows/build-ss3.yml
+++ b/.github/workflows/build-ss3.yml
@@ -217,6 +217,8 @@ jobs:
- name: Build stock synthesis for mac with admb docker image
if: matrix.config.os == 'macos-12'
run: |
+ curl -OSL https://raw.githubusercontent.com/Homebrew/homebrew-core/dc0669eca9479e9eeb495397ba3a7480aaa45c2e/Formula/qemu.rb
+ brew install ./qemu.rb
brew install docker
colima start
docker pull johnoel/admb-13.2:linux
From d3d0fb9f26245a5001cfbd079bc1429e97af2307 Mon Sep 17 00:00:00 2001
From: e-perl-NOAA <elizabeth.gugliotti@noaa.gov>
Date: Thu, 13 Jun 2024 09:01:54 -0400
Subject: [PATCH 11/29] try previous version of colima
---
.github/workflows/build-ss3.yml | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/.github/workflows/build-ss3.yml b/.github/workflows/build-ss3.yml
index 1d944ce7..83cf7228 100644
--- a/.github/workflows/build-ss3.yml
+++ b/.github/workflows/build-ss3.yml
@@ -217,8 +217,9 @@ jobs:
- name: Build stock synthesis for mac with admb docker image
if: matrix.config.os == 'macos-12'
run: |
- curl -OSL https://raw.githubusercontent.com/Homebrew/homebrew-core/dc0669eca9479e9eeb495397ba3a7480aaa45c2e/Formula/qemu.rb
- brew install ./qemu.rb
+ brew uninstall colima
+ sudo mkdir -p /usr/local/bin
+ sudo curl -L -o /usr/local/bin/colima https://github.com/abiosoft/colima/releases/download/v0.6.8/colima-Darwin-x86_64 && sudo chmod +x /usr/local/bin/colima
brew install docker
colima start
docker pull johnoel/admb-13.2:linux
From 123f5aa2d4e36ce104e1f0273f1219e37dbaaf25 Mon Sep 17 00:00:00 2001
From: e-perl-NOAA <elizabeth.gugliotti@noaa.gov>
Date: Thu, 13 Jun 2024 09:09:43 -0400
Subject: [PATCH 12/29] where colima
---
.github/workflows/build-ss3.yml | 5 ++---
1 file changed, 2 insertions(+), 3 deletions(-)
diff --git a/.github/workflows/build-ss3.yml b/.github/workflows/build-ss3.yml
index 83cf7228..8bb5448b 100644
--- a/.github/workflows/build-ss3.yml
+++ b/.github/workflows/build-ss3.yml
@@ -217,9 +217,8 @@ jobs:
- name: Build stock synthesis for mac with admb docker image
if: matrix.config.os == 'macos-12'
run: |
- brew uninstall colima
- sudo mkdir -p /usr/local/bin
- sudo curl -L -o /usr/local/bin/colima https://github.com/abiosoft/colima/releases/download/v0.6.8/colima-Darwin-x86_64 && sudo chmod +x /usr/local/bin/colima
+ brew update
+ brew where colima
brew install docker
colima start
docker pull johnoel/admb-13.2:linux
From 9d23824fe05321e58c362ea3926457bf2016e047 Mon Sep 17 00:00:00 2001
From: e-perl-NOAA <elizabeth.gugliotti@noaa.gov>
Date: Thu, 13 Jun 2024 09:12:38 -0400
Subject: [PATCH 13/29] fix where
---
.github/workflows/build-ss3.yml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/.github/workflows/build-ss3.yml b/.github/workflows/build-ss3.yml
index 8bb5448b..f8d3cebc 100644
--- a/.github/workflows/build-ss3.yml
+++ b/.github/workflows/build-ss3.yml
@@ -218,7 +218,7 @@ jobs:
if: matrix.config.os == 'macos-12'
run: |
brew update
- brew where colima
+ where colima
brew install docker
colima start
docker pull johnoel/admb-13.2:linux
From 8deb43345de04efe6a8c73c15c67991c7d5d4932 Mon Sep 17 00:00:00 2001
From: e-perl-NOAA <elizabeth.gugliotti@noaa.gov>
Date: Thu, 13 Jun 2024 09:22:32 -0400
Subject: [PATCH 14/29] try brew upgrade and --cpu-type max
---
.github/workflows/build-ss3.yml | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/.github/workflows/build-ss3.yml b/.github/workflows/build-ss3.yml
index f8d3cebc..da6b93e8 100644
--- a/.github/workflows/build-ss3.yml
+++ b/.github/workflows/build-ss3.yml
@@ -218,9 +218,9 @@ jobs:
if: matrix.config.os == 'macos-12'
run: |
brew update
- where colima
+ brew upgrade
brew install docker
- colima start
+ colima start --arch x86_64 --cpu-type max
docker pull johnoel/admb-13.2:linux
rm -rf SS330
From f74c36f1f976757bee9857095000f05e1932ff61 Mon Sep 17 00:00:00 2001
From: e-perl-NOAA <elizabeth.gugliotti@noaa.gov>
Date: Thu, 13 Jun 2024 09:33:33 -0400
Subject: [PATCH 15/29] remove brew upgrade
---
.github/workflows/build-ss3.yml | 2 --
1 file changed, 2 deletions(-)
diff --git a/.github/workflows/build-ss3.yml b/.github/workflows/build-ss3.yml
index da6b93e8..e3d4a1f3 100644
--- a/.github/workflows/build-ss3.yml
+++ b/.github/workflows/build-ss3.yml
@@ -217,8 +217,6 @@ jobs:
- name: Build stock synthesis for mac with admb docker image
if: matrix.config.os == 'macos-12'
run: |
- brew update
- brew upgrade
brew install docker
colima start --arch x86_64 --cpu-type max
docker pull johnoel/admb-13.2:linux
From 0468430ef93245adc3c39558e3a6954c8f76884a Mon Sep 17 00:00:00 2001
From: e-perl-NOAA <elizabeth.gugliotti@noaa.gov>
Date: Thu, 13 Jun 2024 09:42:25 -0400
Subject: [PATCH 16/29] try colima delete and reinstall
---
.github/workflows/build-ss3.yml | 5 ++++-
1 file changed, 4 insertions(+), 1 deletion(-)
diff --git a/.github/workflows/build-ss3.yml b/.github/workflows/build-ss3.yml
index e3d4a1f3..8784588d 100644
--- a/.github/workflows/build-ss3.yml
+++ b/.github/workflows/build-ss3.yml
@@ -217,8 +217,11 @@ jobs:
- name: Build stock synthesis for mac with admb docker image
if: matrix.config.os == 'macos-12'
run: |
+ brew update
+ colima delete
+ brew reinstall colima
brew install docker
- colima start --arch x86_64 --cpu-type max
+ colima start --arch x86_64
docker pull johnoel/admb-13.2:linux
rm -rf SS330
From f422412af05216471de3751067b5f2150cd9b4e3 Mon Sep 17 00:00:00 2001
From: e-perl-NOAA <elizabeth.gugliotti@noaa.gov>
Date: Thu, 13 Jun 2024 10:49:52 -0400
Subject: [PATCH 17/29] update comments
---
.github/workflows/build-ss3.yml | 7 +++++--
1 file changed, 5 insertions(+), 2 deletions(-)
diff --git a/.github/workflows/build-ss3.yml b/.github/workflows/build-ss3.yml
index 8784588d..0c1108c4 100644
--- a/.github/workflows/build-ss3.yml
+++ b/.github/workflows/build-ss3.yml
@@ -213,7 +213,6 @@ jobs:
mv SS330/ss3.exe SS330/ss3_win.exe
mv SS330/ss3_opt.exe SS330/ss3_opt_win.exe
-
- name: Build stock synthesis for mac with admb docker image
if: matrix.config.os == 'macos-12'
run: |
@@ -233,8 +232,12 @@ jobs:
colima stop
+ # brew update and colima uninstall and re-install needed to get colima
+ # to work with macOS-12 runner updates
+
# Once macos-12 no longer is kept up by GitHub actions, we will have to
- # build from source the macos-13 and the macos-14 x86_64 images to make versions
+ # build from source the macos-13 (or figure out how to get docker image
+ # to work with macos-13 and the macos-14 x86_64 images to make versions
# compatible with x86_64 machines
- name: Build stock synthesis for mac m2 with admb docker image
From 24c1260450370c9fe264645d6f4919e110769fb9 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Wed, 26 Jun 2024 15:05:50 -0700
Subject: [PATCH 18/29] adding benchmark change
---
SS_benchfore.tpl | 4 ++--
SS_write_report.tpl | 6 +++---
2 files changed, 5 insertions(+), 5 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index c76438e0..13db227f 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -1340,7 +1340,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
MSY_SPR = SSB_equil / SPR_unfished;
SPR_temp = SSB_equil;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- Bmsy = Equ_SpawnRecr_Result(1);
+ 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;
@@ -1433,7 +1433,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
MSY_SPR = SSB_equil / SPR_unfished;
SPR_temp = SSB_equil;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- Bmsy = Equ_SpawnRecr_Result(1);
+ 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
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index c5b13169..66963485 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1802,7 +1802,7 @@ FUNCTION void write_bigoutput()
SS2out << "# " <<endl;
SS2out << "sigmaR: "<< sigmaR << endl;
- SS2out << "SPR_virgin: " << SPR_virgin << " #_uses_biology_from_start_year: " << styr <<endl;
+ SS2out << "SPR0_(virgin): " << SPR_virgin << " #_uses_biology_from_start_year: " << styr <<endl;
switch (SR_fxn)
{
@@ -1820,8 +1820,8 @@ FUNCTION void write_bigoutput()
{
SS2out << "Ln(alpha): " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
SS2out << "Ln(beta): " << SR_parm(4) << " beta " << mfexp(SR_parm(4)) << endl;
- SS2out << "steepness_derived: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " " << SR_parm(3) << endl; // steepness virgin
- SS2out << "R0_derived: " << 1. / beta * (alpha - (1. / SPR_virgin)) << " " << SR_parm(1) << endl; // virgin R0
+ SS2out << "steepness_derived: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << endl; // steepness virgin
+ SS2out << "ln(R0)_derived: " << log( 1. / beta * (alpha - (1. / SPR_virgin))) << endl; // virgin R0
break;
}
case 8:
From 37a4ac97f69ebae928dba76ee2cfff600bfcf943 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Wed, 3 Jul 2024 15:30:57 -0700
Subject: [PATCH 19/29] format changes for spawn_recr report
---
SS_readcontrol_330.tpl | 6 +++++-
SS_write_report.tpl | 40 ++++++++++++++++++++++++----------------
2 files changed, 29 insertions(+), 17 deletions(-)
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index edd7a3a7..fbf18704 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -1979,7 +1979,11 @@
SR_update_SPR0 = 0;
//#_SR_function: 1=null; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=Survival_3Parm; 10=B-H with a,b "<<endl;
- ParmLabel += "SR_LN(R0)";
+ if (SR_fxn == 10)
+ {ParmLabel += "SR_LN(R0)_derived";}
+ else
+ {ParmLabel += "SR_LN(R0)";}
+
switch (SR_fxn)
{
case 1: // previous placement for B-H constrained
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index 66963485..aec89f6d 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1789,6 +1789,7 @@ FUNCTION void write_bigoutput()
SS2out << endl << "#New_Expanded_Spawn_Recr_report" << endl << pick_report_name(19) << endl;
SS2out << "SR_Function: " << SR_fxn << endl;
+ SS2out << "N_SRparms: " << N_SRparm2 << endl;
SS2out << "#" << endl << "parm parm_label value phase" << endl;
for (int j = 1; j <=N_SRparm2; j++)
{
@@ -1796,7 +1797,7 @@ FUNCTION void write_bigoutput()
if (SR_parm_timevary(j) > 0 && j <= 4 ) // timevary SRparm exists
{SS2out << " #_is_time_vary,_so_SRR_updates_base_SPR_annually";}
if (j == (N_SRparm2 - 1) && SR_parm_timevary(j) > 0) // timevary regime exists
- {SS2out << " #_Persistent_deviations_from_SRR_(e.g.regimes)_exist";}
+ {SS2out << " #_Regime_parameter_used_to_offset_from_SRR";}
SS2out << endl;
}
@@ -1820,8 +1821,8 @@ FUNCTION void write_bigoutput()
{
SS2out << "Ln(alpha): " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
SS2out << "Ln(beta): " << SR_parm(4) << " beta " << mfexp(SR_parm(4)) << endl;
- SS2out << "steepness_derived: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << endl; // steepness virgin
SS2out << "ln(R0)_derived: " << log( 1. / beta * (alpha - (1. / SPR_virgin))) << endl; // virgin R0
+ SS2out << "steepness_derived: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << endl; // steepness virgin
break;
}
case 8:
@@ -1833,19 +1834,12 @@ FUNCTION void write_bigoutput()
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);
- SS2out << " Shepherd_c: " << Shepherd_c << " steepness_limit: " << Hupper << " Adjusted_steepness: " << temp << endl;
- break;
- }
- case 9:
- {
- SS2out << " Ricker_Parm1: " << SR_parm(1) << endl;
- SS2out << " Ricker_Parm2: " << SR_parm(2) << endl;
- SS2out << " Ricker_Power: " << SR_parm(3) << endl;
+ SS2out << "Shepherd_c: " << Shepherd_c << endl << "Shepard_steepness_limit: " << Hupper << endl << "Shepard_adjusted_steepness: " << temp << endl;
break;
}
default:
{
- SS2out << "default output needed " << endl;
+ SS2out << "other_SRR " << endl;
break;
}
}
@@ -1865,8 +1859,8 @@ FUNCTION void write_bigoutput()
SS2out << "Bmsy/Bunf: "<< Bmsy / Mgmt_quant(1) << " # using MSY's averaged bio for Bunf" << endl;
SS2out << "#" << endl << "RecDev_method: " << do_recdev << endl << "sum_recdev: " << sum_recdev << endl << "recr_logL: " << recr_like << endl;
- SS2out << recdev_start << " " << recdev_end << " main_recdev:start_end" << endl
- << recdev_adj(1) << " " << recdev_adj(2, 5) << " breakpoints_for_bias_adjustment_ramp " << endl;
+ SS2out << "main_recdev:start_end: " << recdev_start << " " << recdev_end << endl
+ << "breakpoints_for_bias_adjustment_ramp: " <<recdev_adj(1,4) << endl << "max_bias_adj: " << recdev_adj(5) << endl;
temp = sigmaR * sigmaR; // sigmaR^2
SS2out << "ERA N RMSE RMSE^2/sigmaR^2 mean_BiasAdj est_rho Durbin-Watson" << endl;
@@ -1893,8 +1887,17 @@ FUNCTION void write_bigoutput()
}
}
SS2out << endl << "#" << endl << "Initial_equilibrium: " << init_equ_steepness << " # 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 SPR0_curr h_curr R0_curr P1 P2 P3 P4" << endl;
+ if (SR_fxn == 10) SS2out << "#_Note:_h_curr_and_R0_curr_are_for_info_only;_calculated_from_alpha_beta_and_current_SPR0" << endl;
+ SS2out << "#_columns_with_P_will_show_time_vary_SR_parameters" << endl << "#" << endl;
+ SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SPR0_curr ";
+ if(SR_fxn == 10)
+ {SS2out << "h_curr R0_curr ";}
+ else
+ {SS2out << "NA1 NA2 ";}
+;
+ for (j = 1; j <= N_SRparm2; j++)
+ {SS2out << "P" << j << " ";}
+ SS2out << endl;
SS2out << "S/Rcurve " << SSB_virgin << " " << Recr_virgin << endl;
y = styr - 2;
SS2out << "Virg " << SSB_yr(y) << " " << exp_rec(y) << " - " << 0.0 << " Virg " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " 0.0 " << endl;
@@ -1947,9 +1950,14 @@ FUNCTION void write_bigoutput()
}
dvariable SPR_curr = Smry_Table(y, 11) / Recr_virgin;
SS2out << " " << SPR_curr << " ";
+ if (SR_fxn == 10)
+ {
SS2out << alpha * SPR_curr / (4. + alpha * SPR_curr) << " "; // steepness with current SPR
SS2out << 1. / beta * (alpha - (1. / SPR_curr)) << " "; // R0 with current SPR
- SS2out << SR_parm_byyr(y)(1,4) << endl;
+ }
+ else
+ {SS2out << " - - ";}
+ SS2out << SR_parm_byyr(y)(1,N_SRparm2) << endl;
}
// REPORT_KEYWORD SPAWN_RECR_CURVE
From b67a60e8950f5a55355952d38fea664c314ec800 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Tue, 9 Jul 2024 15:13:53 -0700
Subject: [PATCH 20/29] fixes to ss_warn
---
Compile/Make_SS_warn.bat | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/Compile/Make_SS_warn.bat b/Compile/Make_SS_warn.bat
index 3fdb0d20..af9b1430 100644
--- a/Compile/Make_SS_warn.bat
+++ b/Compile/Make_SS_warn.bat
@@ -19,15 +19,15 @@ cd Compile
if exist ss3.exe (
if exist ss3_old.exe (
- del ss3old.exe
+ del ss3_old.exe
)
ren ss3.exe ss3_old.exe
)
tpl2cpp ss3
-g++ -c -std=c++14 -O2 -D_FILE_OFFSET_BITS=64 -DUSE_ADMB_CONTRIBS -D_USE_MATH_DEFINES -I. -I"C:\ADMB-13.1\include" -I"C:\ADMB-13.1\include\contrib" -Wall -Wextra -o ss3.obj ss3.cpp
+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.1\lib\libadmb-contrib-mingw64-g++8.a"
+g++ -static -o ss3.exe ss3.obj "C:\ADMB-13.2\lib\libadmb-contrib-mingw64-g++13.a"
dir *.exe
From 0666d6fa5a6386a3b761a55e7aaf9c8842a4b8d2 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Tue, 13 Aug 2024 20:06:06 -0700
Subject: [PATCH 21/29] WIP-more work needed on benchmark calcs
---
SS_benchfore.tpl | 30 ++++++++++++++++++++------
SS_param.tpl | 1 +
SS_popdyn.tpl | 21 +++++++++---------
SS_recruit.tpl | 52 ++++++++++++++++++++++-----------------------
SS_write_report.tpl | 6 +++---
5 files changed, 63 insertions(+), 47 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 13db227f..8bb12412 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -744,6 +744,18 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SR_parm_work(j) = temp / (Bmark_Yr(10) - Bmark_Yr(9) + 1.);
}
}
+
+ if(SR_fxn == 10) // B-H with alpha, beta
+ {
+ alpha = mfexp(SR_parm_work(3));
+ beta = mfexp(SR_parm_work(4));
+ Fishon = 0;
+ Recr_unf = 1.0;
+ Do_Equil_Calc(Recr_unf);
+ SPR_virgin_adj = SSB_equil / 1.0;
+ SR_parm_work(2) = alpha * SPR_virgin_adj / (4. + alpha * SPR_virgin_adj); // steepness
+ SR_parm_work(1) = log(1. / beta * (alpha - (1. / SPR_virgin_adj))); // ln(R0)
+ }
Fishon = 0;
Recr_unf = mfexp(SR_parm_work(1));
Do_Equil_Calc(Recr_unf);
@@ -751,13 +763,17 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SPR_unfished = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
if (show_MSY == 1)
{
- report5 << "SR_parm for benchmark: " << SR_parm_work << endl
+ report5 << "SR_parms for benchmark: " << SR_parm_work << endl
<< "mean from years: " << Bmark_Yr(9) << " " << Bmark_Yr(10) << endl;
// SPR_virgin = SSB_virgin / Recr_virgin; // already defined
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Virgin SPR0, SSB, R: " << SPR_virgin << " " << Equ_SpawnRecr_Result << endl;
+ report5 << " Virgin SSB, R0: " << SSB_virgin << " " << Recr_virgin << " " << SPR_virgin_adj << endl;
+ report5 << " unfished SSB, R0: " << SSB_unf << " " << Recr_unf << " " << SPR_unfished << " with current biology " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_virgin_adj); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Virgin SPR0, result_SSB, R: " << SPR_virgin_adj << " " << Equ_SpawnRecr_Result << endl << endl;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_unfished); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Benchmark SPR0, SSB, R: " << SPR_unfished << " " << Equ_SpawnRecr_Result << endl;
+ report5 << " Benchmark SPR0, result_SSB, R: " << SPR_unfished << " " << Equ_SpawnRecr_Result << endl << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_unfished); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Benchmark SPR0, result_SSB, R: " << SPR_unfished << " " << Equ_SpawnRecr_Result << " with virgin spawn-recr " << endl;
}
SR_parm_work(N_SRparm2 + 1) = SSB_unf;
Mgmt_quant(1) = SSB_unf;
@@ -1930,7 +1946,7 @@ FUNCTION void Get_Forecast()
dvariable OFL_catch;
dvariable Fcast_Crash;
dvariable totcatch;
- dvariable R0_use;
+ dvariable R0_use; // annually updated variable if SR_update_SPR0 == 1
dvariable SSB_use;
dvar_matrix catage_w(1, gmorph, 0, nages);
dvar_vector tempcatch(1, Nfleet);
@@ -2577,7 +2593,7 @@ FUNCTION void Get_Forecast()
}
// SPAWN-RECR: get recruitment in forecast; needs to be area-specific
// SR_fxn
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ if (SR_update_SPR0 == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
@@ -3215,7 +3231,7 @@ FUNCTION void Get_Forecast()
// 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-specific
// SR_fxn
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ if (SR_update_SPR0 == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
diff --git a/SS_param.tpl b/SS_param.tpl
index bd21a12b..3d000222 100644
--- a/SS_param.tpl
+++ b/SS_param.tpl
@@ -166,6 +166,7 @@ PARAMETER_SECTION
number half_sigmaRsq;
number sigmaR;
number SPR_virgin;
+ number SPR_virgin_adj;
number regime_change;
number rho;
number dirichlet_Parm;
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index 3891000a..2f0213c5 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -344,7 +344,7 @@ FUNCTION void get_initial_conditions()
equ_Recr = 1.0;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation. Returns SPR because R = 1.0
SPR_virgin = SSB_equil; // spawners per recruit. Needed for Sr_fxn = 10
-
+ SPR_virgin_adj = SSB_equil;
if(SR_fxn == 10) // B-H with a,b
{
// WHAM based on R = A*S/(1+B*S)
@@ -355,8 +355,8 @@ FUNCTION void get_initial_conditions()
alpha = mfexp(SR_parm(3));
beta = mfexp(SR_parm(4));
- steepness = alpha * SPR_virgin / (4. + alpha * SPR_virgin);
- Recr_virgin = 1. / beta * (alpha - (1. / SPR_virgin));
+ steepness = alpha * SPR_virgin_adj / (4. + alpha * SPR_virgin_adj);
+ Recr_virgin = 1. / beta * (alpha - (1. / SPR_virgin_adj));
// warning << " before AB_calcs " << "parm " << SR_parm(1) << " calc " << log(Recr_virgin) << endl;
SR_parm(1) = log(Recr_virgin);
SR_parm(2) = steepness;
@@ -380,7 +380,7 @@ FUNCTION void get_initial_conditions()
exp_rec(eq_yr, 4) = Recr_virgin;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
SSB_virgin = SSB_equil;
- SPR_virgin = SSB_equil / Recr_virgin; // spawners per recruit
+// SPR_virgin = SSB_equil / Recr_virgin; // spawners per recruit already calculated
if(Do_Benchmark==0)
{
Mgmt_quant(1)=SSB_virgin;
@@ -696,7 +696,7 @@ FUNCTION void get_time_series()
dvariable crashtemp1;
dvariable interim_tot_catch;
dvariable Z_adjuster;
- dvariable R0_use;
+ dvariable R0_use; // annually updated variable if SR_update_SPR0 == 1; gets passed to Spawn_Recr() function
dvariable SSB_use;
if (Do_Morphcomp > 0)
@@ -1022,10 +1022,10 @@ FUNCTION void get_time_series()
}
}
- // SS_Label_Info_24.2.3 #Get the total recruitment produced by this spawning biomass
+ // SS_Label_Info_24.2.3 #Get the total recruitment produced by this spawning biomass at the beginning of the season
// SPAWN-RECR: calc recruitment in time series; need to make this area-specific
// SR_Fxn relevant keyword
- if (SR_parm_timevary(1) == 0) // R0 is not time-varying
+ if (SR_update_SPR0 == 0) // SRparm are not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
@@ -1033,12 +1033,12 @@ FUNCTION void get_time_series()
else
{
R0_use = mfexp(SR_parm_work(1));
- warning << y << " set R0use to SRparm_work " << R0_use << " vir: " << Recr_virgin << endl;
equ_Recr = R0_use;
Fishon = 0;
eq_yr = y;
bio_yr = y;
Do_Equil_Calc(R0_use); // call function to do equilibrium calculation
+ SSB_use = SSB_equil;
if (fishery_on_off == 1)
{
Fishon = 1;
@@ -1047,7 +1047,6 @@ FUNCTION void get_time_series()
{
Fishon = 0;
}
- SSB_use = SSB_equil;
}
Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr
if (SR_fxn != 7) apply_recdev(Recruits, R0_use); // apply recruitment deviation
@@ -1482,7 +1481,7 @@ FUNCTION void get_time_series()
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
// SR_fxn
if (SR_update_SPR0 == 0) // SR parms are not time-varying
@@ -1498,6 +1497,7 @@ FUNCTION void get_time_series()
eq_yr = y;
bio_yr = y;
Do_Equil_Calc(R0_use); // call function to do equilibrium calculation
+ SSB_use = SSB_equil;
if (fishery_on_off == 1)
{
Fishon = 1;
@@ -1506,7 +1506,6 @@ FUNCTION void get_time_series()
{
Fishon = 0;
}
- SSB_use = SSB_equil;
}
Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 2065b277..179ecd6e 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -268,7 +268,7 @@ 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 dvar_vector& SRparm,
- const prevariable& SSB_virgin, const prevariable& Recr_virgin, const prevariable& SPR_temp)
+ const prevariable& SSB_temp, const prevariable& RECR_temp, const prevariable& SPR_temp)
{
RETURN_ARRAYS_INCREMENT();
dvar_vector Equil_Spawn_Recr_Calc(1, 2); // values to return 1 is B_equil, 2 is R_equil
@@ -297,8 +297,8 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// 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_temp * (1. + (log(RECR_temp / SSB_temp) + log(SPR_temp)) / steepness);
+ R_equil = RECR_temp * B_equil / SSB_temp * mfexp(steepness * (1. - B_equil / SSB_temp));
break;
}
@@ -316,14 +316,14 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS3 previously used alternative formulation: R = A*S/(B+S)
// converting SS3 to align with WHAM
- // SPR_virgin = SSB_virgin / Recr_virgin; // this is already defined
- alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
+ // SPR_virgin = SSB_temp / RECR_temp; // this is already defined
+ alpha = 4.0 * steepness / (SPR_virgin_adj * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
// " h " << steepness << " derive " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " " << endl;
- // " R0 " << Recr_virgin << " derive " << 1. / beta * (alpha - 1./SPR_virgin) << endl;
-// report5 <<" SSB_unf "<<SSB_virgin<<" SPR_unf "<<SPR_virgin<<" steep: "<<steepness<<" R0: "<<Recr_virgin << endl;
-// report5 <<" derive_alpha "<<alpha<<" derive_beta "<<beta << endl;
+ // " R0 " << RECR_temp << " derive " << 1. / beta * (alpha - 1./SPR_virgin) << endl;
+// report5 <<" SSB_unf "<<SSB_temp<<" SPR_unf "<<SPR_virgin<<" steep: "<<steepness<<" R0: "<<RECR_temp << endl;
+ report5 <<" derive_alpha "<<alpha<<" derive_beta "<<beta << endl;
// report5 << " deriv_h: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " derive_R0: " << 1. / beta * (alpha - (1. / SPR_virgin))<<endl;
B_equil = (alpha * SPR_temp - 1.0) / beta;
B_equil = posfun(B_equil, 0.0001, temp);
@@ -347,27 +347,27 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS_Label_44.1.4 constant recruitment
case 4: // constant; no bias correction
{
- B_equil = SPR_temp * Recr_virgin;
- R_equil = Recr_virgin;
+ B_equil = SPR_temp * RECR_temp;
+ R_equil = RECR_temp;
break;
}
// SS_Label_44.1.5 Hockey Stick
case 5: // hockey stick
{
- dvariable hockey_min = SRparm(3) * Recr_virgin; // min recruitment level
- // temp=SSB_virgin/R0*steepness; // spawners per recruit at inflection
- dvariable hockey_slope = (Recr_virgin - hockey_min) / (steepness * SSB_virgin); // slope of recruitment on spawners below the inflection
- B_equil = Join_Fxn(0.0 * SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin * steepness, SPR_temp, hockey_min / ((1. / SPR_temp) - hockey_slope), SPR_temp * Recr_virgin);
- R_equil = Join_Fxn(0.0 * SSB_virgin, SSB_virgin, SSB_virgin * steepness, B_equil, hockey_min + hockey_slope * B_equil, Recr_virgin);
+ dvariable hockey_min = SRparm(3) * RECR_temp; // min recruitment level
+ // temp=SSB_temp/R0*steepness; // spawners per recruit at inflection
+ dvariable hockey_slope = (RECR_temp - hockey_min) / (steepness * SSB_temp); // slope of recruitment on spawners below the inflection
+ B_equil = Join_Fxn(0.0 * SSB_temp / RECR_temp, SSB_temp / RECR_temp, SSB_temp / RECR_temp * steepness, SPR_temp, hockey_min / ((1. / SPR_temp) - hockey_slope), SPR_temp * RECR_temp);
+ R_equil = Join_Fxn(0.0 * SSB_temp, SSB_temp, SSB_temp * steepness, B_equil, hockey_min + hockey_slope * B_equil, RECR_temp);
break;
}
// SS_Label_44.1.7 3 parameter survival based
case 7: // survival
{
- SRZ_0 = log(1.0 / (SSB_virgin / Recr_virgin));
+ SRZ_0 = log(1.0 / (SSB_temp / RECR_temp));
srz_min = SRZ_0 * (1.0 - steepness);
- B_equil = SSB_virgin * (1. - (log(1. / SPR_temp) - SRZ_0) / pow((srz_min - SRZ_0), (1. / SRparm(3))));
- SRZ_surv = mfexp((1. - pow((B_equil / SSB_virgin), SRparm(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
+ B_equil = SSB_temp * (1. - (log(1. / SPR_temp) - SRZ_0) / pow((srz_min - SRZ_0), (1. / SRparm(3))));
+ SRZ_surv = mfexp((1. - pow((B_equil / SSB_temp), SRparm(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
R_equil = B_equil * SRZ_surv;
break;
}
@@ -389,10 +389,10 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
Shepherd_c2 = pow(0.2, SRparm(3));
Hupper = 1.0 / (5.0 * Shepherd_c2);
steepness = 0.2 + (SRparm(2) - 0.2) / (0.8) * (Hupper - 0.2);
- Shep_top = 5.0 * steepness * (1.0 - Shepherd_c2) * (SPR_temp * Recr_virgin) / SSB_virgin - (1.0 - 5.0 * steepness * Shepherd_c2);
+ Shep_top = 5.0 * steepness * (1.0 - Shepherd_c2) * (SPR_temp * RECR_temp) / SSB_temp - (1.0 - 5.0 * steepness * Shepherd_c2);
Shep_bot = 5.0 * steepness - 1.0;
Shep_top2 = posfun(Shep_top, 0.001, temp);
- R_equil = (SSB_virgin / SPR_temp) * pow((Shep_top2 / Shep_bot), (1.0 / SRparm(3)));
+ R_equil = (SSB_temp / SPR_temp) * pow((Shep_top2 / Shep_bot), (1.0 / SRparm(3)));
B_equil = R_equil * SPR_temp;
break;
}
@@ -402,10 +402,10 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
{
steepness = SRparm(2);
dvariable RkrPower = SRparm(3);
- temp = SSB_virgin / (SPR_temp * Recr_virgin);
+ temp = SSB_temp / (SPR_temp * RECR_temp);
dvariable RkrTop = pow(0.8, RkrPower) * log(temp) / log(5.0 * steepness);
RkrTop = posfun(RkrTop, 0.000001, CrashPen);
- R_equil = temp * Recr_virgin * (1.0 - pow(RkrTop, 1.0 / RkrPower));
+ R_equil = temp * RECR_temp * (1.0 - pow(RkrTop, 1.0 / RkrPower));
B_equil = R_equil * SPR_temp;
break;
}
@@ -428,10 +428,10 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
Hupper=1.0/(5.0*Shepherd_c2);
steepness=0.20001+((0.8)/(1.0+exp(-SRparm2))-0.2)/(0.8)*(Hupper-0.2);
// steep2=0.20001+(steepness-0.2)/(0.8)*(Hupper-0.2);
- Shep_top=5.0*steepness*(1.0-Shepherd_c2)*(SPR_temp*Recr_virgin)/SSB_virgin-(1.0-5.0*steepness*Shepherd_c2);
+ Shep_top=5.0*steepness*(1.0-Shepherd_c2)*(SPR_temp*RECR_temp)/SSB_temp-(1.0-5.0*steepness*Shepherd_c2);
Shep_bot=5.0*steepness-1.0;
Shep_top2=posfun(Shep_top,0.001,temp);
- R_equil=(SSB_virgin/SPR_temp) * pow((Shep_top2/Shep_bot),(1.0/Shepherd_c));
+ R_equil=(SSB_temp/SPR_temp) * pow((Shep_top2/Shep_bot),(1.0/Shepherd_c));
B_equil=R_equil*SPR_temp;
break;
}
@@ -448,10 +448,10 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// if (Recs < 0) Rec2 = 0; else Rec2 = Recs;
steepness = 0.2 + (10.0 - 0.2)/(1+exp(-SR_parm_work(2)));
dvariable RkrPower=exp(SR_parm_work(3));
- temp=SSB_virgin/(SPR_temp*Recr_virgin);
+ temp=SSB_virgin/(SPR_temp*RECR_temp);
dvariable RkrTop = pow(0.8,RkrPower)*log(temp)/log(5.0*steepness);
RkrTop = posfun(RkrTop,0.000001,CrashPen);
- R_equil = temp *Recr_virgin * (1.0 - pow(RkrTop,1.0/RkrPower));
+ R_equil = temp *RECR_temp * (1.0 - pow(RkrTop,1.0/RkrPower));
B_equil=R_equil*SPR_temp;
break;
}
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index aec89f6d..25e2f936 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1809,7 +1809,7 @@ FUNCTION void write_bigoutput()
{
case 3: // Beverton-Holt with steepness
{
- alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
+ alpha = 4.0 * steepness / (SPR_virgin_adj * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1. - steepness) * SSB_virgin);
SS2out << "Ln(R0): " << SR_parm(1) << endl << "R0: " << mfexp(SR_parm(1)) << endl;
SS2out << "steepness: " << steepness << endl;
@@ -1821,8 +1821,8 @@ FUNCTION void write_bigoutput()
{
SS2out << "Ln(alpha): " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
SS2out << "Ln(beta): " << SR_parm(4) << " beta " << mfexp(SR_parm(4)) << endl;
- SS2out << "ln(R0)_derived: " << log( 1. / beta * (alpha - (1. / SPR_virgin))) << endl; // virgin R0
- SS2out << "steepness_derived: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << endl; // steepness virgin
+ SS2out << "ln(R0)_derived: " << log( 1. / beta * (alpha - (1. / SPR_virgin_adj))) << endl; // virgin R0
+ SS2out << "steepness_derived: " << alpha * SPR_virgin_adj / (4. + alpha * SPR_virgin_adj) << endl; // steepness virgin
break;
}
case 8:
From f2df4f49c177e718ad4ac8aa13aabe4adb26102a Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Thu, 15 Aug 2024 17:02:43 -0700
Subject: [PATCH 22/29] WIP2-need-to-fix-SPR-passing_to_benchmark
---
SS_benchfore.tpl | 36 ++++++++++++++++++------------------
SS_popdyn.tpl | 17 ++++++++---------
2 files changed, 26 insertions(+), 27 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 8bb12412..4bf41b98 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -536,7 +536,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
dvariable last_F1;
dvariable Closer;
dvariable Vbio1_unfished;
- dvariable SPR_unfished;
+ dvariable SPR_unf;
dvariable Vbio_MSY;
dvariable Vbio1_MSY;
dvariable junk;
@@ -751,7 +751,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
beta = mfexp(SR_parm_work(4));
Fishon = 0;
Recr_unf = 1.0;
- Do_Equil_Calc(Recr_unf);
+ Do_Equil_Calc(Recr_unf); //
SPR_virgin_adj = SSB_equil / 1.0;
SR_parm_work(2) = alpha * SPR_virgin_adj / (4. + alpha * SPR_virgin_adj); // steepness
SR_parm_work(1) = log(1. / beta * (alpha - (1. / SPR_virgin_adj))); // ln(R0)
@@ -759,21 +759,21 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Fishon = 0;
Recr_unf = mfexp(SR_parm_work(1));
Do_Equil_Calc(Recr_unf);
- SSB_unf = SSB_equil;
- SPR_unfished = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
+ SSB_unf = SSB_equil; // equilibrium unfished SSB using the benchmark averaged Recr_unf and benchmark averaged biology
+ SPR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
if (show_MSY == 1)
{
report5 << "SR_parms for benchmark: " << SR_parm_work << endl
<< "mean from years: " << Bmark_Yr(9) << " " << Bmark_Yr(10) << endl;
// SPR_virgin = SSB_virgin / Recr_virgin; // already defined
report5 << " Virgin SSB, R0: " << SSB_virgin << " " << Recr_virgin << " " << SPR_virgin_adj << endl;
- report5 << " unfished SSB, R0: " << SSB_unf << " " << Recr_unf << " " << SPR_unfished << " with current biology " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl;
+ report5 << " unfished SSB, R0: " << SSB_unf << " " << Recr_unf << " " << SPR_unf << " with current biology " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_virgin_adj); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
report5 << " Virgin SPR0, result_SSB, R: " << SPR_virgin_adj << " " << Equ_SpawnRecr_Result << endl << endl;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_unfished); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Benchmark SPR0, result_SSB, R: " << SPR_unfished << " " << Equ_SpawnRecr_Result << endl << endl;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_unfished); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Benchmark SPR0, result_SSB, R: " << SPR_unfished << " " << Equ_SpawnRecr_Result << " with virgin spawn-recr " << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Benchmark SPR0, result_SSB, R: " << SPR_unf << " " << Equ_SpawnRecr_Result << endl << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Benchmark SPR0, result_SSB, R: " << SPR_unf << " " << Equ_SpawnRecr_Result << " with virgin spawn-recr " << endl;
}
SR_parm_work(N_SRparm2 + 1) = SSB_unf;
Mgmt_quant(1) = SSB_unf;
@@ -808,7 +808,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SPR_target100 = SPR_target * 100.;
Do_Equil_Calc(equ_Recr);
- SPR_unfished = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
+ SPR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
Vbio1_unfished = smrybio; // gets value from equil_calc
if (show_MSY == 1)
{
@@ -854,7 +854,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Fishon = 1;
Do_Equil_Calc(equ_Recr);
- yld1(ii) = 100. * SSB_equil / SPR_unfished; // spawning potential ratio
+ yld1(ii) = 100. * SSB_equil / SPR_unf; // spawning potential ratio
}
SPR_actual = yld1(1); // spawning potential ratio
@@ -1008,7 +1008,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 / SPR_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;
@@ -1051,7 +1051,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 / SPR_unf;
Vbio_Btgt = totbio;
Vbio1_Btgt = smrybio;
Mgmt_quant(7) = equ_F_std;
@@ -1131,7 +1131,7 @@ 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;
+ SPR_Btgt = SSB_equil / SPR_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, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
@@ -1353,7 +1353,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Do_Equil_Calc(equ_Recr);
// SPAWN-RECR: calc spawn-recr for MSY calcs; need to make area-specific
- MSY_SPR = SSB_equil / SPR_unfished;
+ MSY_SPR = SSB_equil / SPR_unf;
SPR_temp = SSB_equil;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_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
@@ -1446,7 +1446,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
Do_Equil_Calc(equ_Recr);
// SPAWN-RECR: calc spawn-recr for MSY calcs; need to make area-specific
- MSY_SPR = SSB_equil / SPR_unfished;
+ MSY_SPR = SSB_equil / SPR_unf;
SPR_temp = SSB_equil;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Bmsy = Equ_SpawnRecr_Result(1); // MSY is directly estimated
@@ -1739,7 +1739,7 @@ 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;
+ SPR_Btgt2 = SSB_equil / SPR_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, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
@@ -1947,7 +1947,7 @@ FUNCTION void Get_Forecast()
dvariable Fcast_Crash;
dvariable totcatch;
dvariable R0_use; // annually updated variable if SR_update_SPR0 == 1
- dvariable SSB_use;
+ dvariable SSB_use; // selected version of SSB that gets passes to Spawn_Recr
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
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index 2f0213c5..f509cf93 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -344,7 +344,7 @@ FUNCTION void get_initial_conditions()
equ_Recr = 1.0;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation. Returns SPR because R = 1.0
SPR_virgin = SSB_equil; // spawners per recruit. Needed for Sr_fxn = 10
- SPR_virgin_adj = SSB_equil;
+ SPR_virgin_adj = SSB_equil; // also needed for Sr_fxn 10. Will get revised in benchmark to use averaged biology if requested.
if(SR_fxn == 10) // B-H with a,b
{
// WHAM based on R = A*S/(1+B*S)
@@ -381,14 +381,14 @@ FUNCTION void get_initial_conditions()
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
SSB_virgin = SSB_equil;
// SPR_virgin = SSB_equil / Recr_virgin; // spawners per recruit already calculated
- 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;
+ 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;
}
Smry_Table(styr - 2, 1) = totbio; // from equil calcs
Smry_Table(styr - 2, 2) = smrybio; // from equil calcs
@@ -1492,7 +1492,6 @@ FUNCTION void get_time_series()
else // update SSB_use and R0_use where will update SPR0 inside the Spawn_recr fxn
{
R0_use = mfexp(SR_parm_work(1)); // check to be sure this works when R0 is derived from B-H with alpha, beta parameters
- equ_Recr = R0_use;
Fishon = 0;
eq_yr = y;
bio_yr = y;
From a9dacb03d6b28badda5ba3b23dd0b27e89e2fda1 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Fri, 13 Sep 2024 16:27:52 -0700
Subject: [PATCH 23/29] WIP: add switch in benchmark to use correct SPR0
---
SS_benchfore.tpl | 42 +++++++--------
SS_popdyn.tpl | 5 +-
SS_readcontrol_330.tpl | 2 +-
SS_recruit.tpl | 116 ++++++++++++++++++++---------------------
4 files changed, 80 insertions(+), 85 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 4bf41b98..25038380 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -745,35 +745,30 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
}
- if(SR_fxn == 10) // B-H with alpha, beta
+ if (SR_update_SPR0 == 0) // use virgin biology for the spawner-recruitment R0,h calculations
{
- alpha = mfexp(SR_parm_work(3));
- beta = mfexp(SR_parm_work(4));
- Fishon = 0;
- Recr_unf = 1.0;
- Do_Equil_Calc(Recr_unf); //
- SPR_virgin_adj = SSB_equil / 1.0;
- SR_parm_work(2) = alpha * SPR_virgin_adj / (4. + alpha * SPR_virgin_adj); // steepness
- SR_parm_work(1) = log(1. / beta * (alpha - (1. / SPR_virgin_adj))); // ln(R0)
+ Recr_unf = Recr_virgin;
+ SSB_unf = SSB_virgin;
+ SPR_unf = SSB_unf / Recr_unf;
+ }
+ else // use benchmark biology in the spawner-recruitment R0,h calculations
+ {
+ Fishon = 0;
+ Recr_unf = mfexp(SR_parm_work(1));
+ Do_Equil_Calc(Recr_unf);
+ SSB_unf = SSB_equil; // equilibrium unfished SSB using the benchmark averaged Recr_unf and benchmark averaged biology
+ SPR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
}
- Fishon = 0;
- Recr_unf = mfexp(SR_parm_work(1));
- Do_Equil_Calc(Recr_unf);
- SSB_unf = SSB_equil; // equilibrium unfished SSB using the benchmark averaged Recr_unf and benchmark averaged biology
- SPR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
if (show_MSY == 1)
{
report5 << "SR_parms for benchmark: " << SR_parm_work << endl
- << "mean from years: " << Bmark_Yr(9) << " " << Bmark_Yr(10) << endl;
- // SPR_virgin = SSB_virgin / Recr_virgin; // already defined
- report5 << " Virgin SSB, R0: " << SSB_virgin << " " << Recr_virgin << " " << SPR_virgin_adj << endl;
- report5 << " unfished SSB, R0: " << SSB_unf << " " << Recr_unf << " " << SPR_unf << " with current biology " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_virgin_adj); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Virgin SPR0, result_SSB, R: " << SPR_virgin_adj << " " << Equ_SpawnRecr_Result << endl << endl;
+ << "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl;
+ if ( SR_update_SPR0 == 1) report5 << "SPR0 for equilibrium spawner-recruit based on benchmark biology, not virgin biology" << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Virgin SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SPR_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Benchmark SPR0, result_SSB, R: " << SPR_unf << " " << Equ_SpawnRecr_Result << endl << endl;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Benchmark SPR0, result_SSB, R: " << SPR_unf << " " << Equ_SpawnRecr_Result << " with virgin spawn-recr " << endl;
+ report5 << " Benchmark SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " " << SPR_unf << " equil: " << Equ_SpawnRecr_Result << endl;
}
SR_parm_work(N_SRparm2 + 1) = SSB_unf;
Mgmt_quant(1) = SSB_unf;
@@ -925,7 +920,6 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
// SPAWN-RECR: calc equil spawn-recr in YPR; need to make this area-specific
SPR_temp = SSB_equil; // based on most recent call to Do_Equil_Calc
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5<<SPR_temp<<" " << Equ_SpawnRecr_Result<<endl;
Bspr = Equ_SpawnRecr_Result(1);
Bspr_rec = Equ_SpawnRecr_Result(2);
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index f509cf93..a0297e3c 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -535,6 +535,7 @@ FUNCTION void get_initial_conditions()
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
+ // 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(SR_parm_work, SSB_virgin, Recr_virgin, SPR_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;
@@ -1037,7 +1038,7 @@ FUNCTION void get_time_series()
Fishon = 0;
eq_yr = y;
bio_yr = y;
- Do_Equil_Calc(R0_use); // call function to do equilibrium calculation
+ Do_Equil_Calc(R0_use); // call function to do equilibrium calculation with current year's biology and adjusted R0
SSB_use = SSB_equil;
if (fishery_on_off == 1)
{
@@ -1048,7 +1049,7 @@ FUNCTION void get_time_series()
Fishon = 0;
}
}
- Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr
+ Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr using either virgin or adjusted R0 and SSB0
if (SR_fxn != 7) apply_recdev(Recruits, R0_use); // apply recruitment deviation
// distribute Recruitment of age 0 fish among the current and future settlements; and among areas and morphs
// use t offset for each birth event: Settlement_offset(settle)
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index fbf18704..905f9a81 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -6453,7 +6453,7 @@
active_parm(CoVar_Count) = j;
if (y == styr - 2)
{
- ParmLabel += "Recr_Virgin";
+ ParmLabel += "Recr_virgin";
}
else if (y == styr - 1)
{
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 179ecd6e..5766efce 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -6,7 +6,7 @@
//********************************************************************
/* SS_Label_FUNCTION 43 Spawner-recruitment function */
// SPAWN-RECR: function: to calc R from S
-FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariable& Recr_virgin_adj, const prevariable& SSB_current)
+FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_use, const prevariable& Recr_virgin_use, const prevariable& SSB_current)
{
RETURN_ARRAYS_INCREMENT();
dvariable NewRecruits;
@@ -22,7 +22,7 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
dvariable SRZ_0;
dvariable srz_min;
dvariable SRZ_surv;
-// warning << y << " Tester_R0 " << Recr_virgin_adj << " SSB0 " << SSB_virgin_adj << " SSB_curr: " << SSB_current << endl;
+// 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
@@ -30,7 +30,7 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
regime_change = SR_parm_work(N_SRparm2 - 1); // this is a persistent deviation off the S/R curve
// SS_Label_43.3 calculate expected recruitment from the input spawning biomass and the SR curve
- // functions below use Recr_virgin_adj,SSB_virgin_adj which could have been adjusted adjusted above from R0,SSB_virgin
+ // functions below use Recr_virgin_use,SSB_virgin_use which could have been adjusted adjusted above from R0,SSB_virgin
switch (SR_fxn)
{
case 1: // previous placement for B-H constrained
@@ -43,15 +43,15 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 2: // ricker
{
steepness = SR_parm_work(2);
- NewRecruits = Recr_virgin_adj * SSB_curr_adj / SSB_virgin_adj * mfexp(steepness * (1. - SSB_curr_adj / SSB_virgin_adj));
+ NewRecruits = Recr_virgin_use * SSB_curr_adj / SSB_virgin_use * mfexp(steepness * (1. - SSB_curr_adj / SSB_virgin_use));
break;
}
// SS_Label_43.3.3 Beverton-Holt
case 3: // Beverton-Holt
{
steepness = SR_parm_work(2);
- NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_curr_adj) /
- (SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);
+ NewRecruits = (4. * steepness * Recr_virgin_use * SSB_curr_adj) /
+ (SSB_virgin_use * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);
break;
}
@@ -66,7 +66,7 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
// 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
@@ -74,8 +74,8 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
// the 3rd parameter allows for a minimum recruitment level
{
steepness = SR_parm_work(2);
- temp = SR_parm_work(3) * Recr_virgin_adj + SSB_curr_adj / (steepness * SSB_virgin_adj) * (Recr_virgin_adj - SR_parm_work(3) * Recr_virgin_adj); // linear decrease below steepness*SSB_virgin_adj
- NewRecruits = Join_Fxn(0.0 * SSB_virgin_adj, SSB_virgin_adj, steepness * SSB_virgin_adj, SSB_curr_adj, temp, Recr_virgin_adj);
+ temp = SR_parm_work(3) * Recr_virgin_use + SSB_curr_adj / (steepness * SSB_virgin_use) * (Recr_virgin_use - SR_parm_work(3) * Recr_virgin_use); // linear decrease below steepness*SSB_virgin_use
+ NewRecruits = Join_Fxn(0.0 * SSB_virgin_use, SSB_virgin_use, steepness * SSB_virgin_use, SSB_curr_adj, temp, Recr_virgin_use);
break;
}
@@ -83,32 +83,32 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 6: //Beverton-Holt constrained
{
steepness = SR_parm_work(2);
-// dvariable SPR = SSB_virgin_adj / Recr_virgin;
+// 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_adj)
+ if (SSB_curr_adj > SSB_virgin_use)
{
- SSB_BH1 = SSB_virgin_adj;
+ SSB_BH1 = SSB_virgin_use;
}
else
{
SSB_BH1 = SSB_curr_adj;
}
- NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_BH1) / (SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_BH1);
+ NewRecruits = (4. * steepness * Recr_virgin_use * SSB_BH1) / (SSB_virgin_use * (1. - steepness) + (5. * steepness - 1.) * SSB_BH1);
break;
}
// SS_Label_43.3.7 survival based
case 7: // survival based, so constrained such that recruits cannot exceed fecundity
{
- // PPR_0=SSB_virgin_adj/Recr_virgin_adj; // pups per recruit at virgin
+ // PPR_0=SSB_virgin_use/Recr_virgin_use; // pups per recruit at virgin
// Surv_0=1./PPR_0; // recruits per pup at virgin
- // Pups_0=SSB_virgin_adj; // total population fecundity is the number of pups produced
+ // Pups_0=SSB_virgin_use; // total population fecundity is the number of pups produced
// Sfrac=SR_parm(2);
- SRZ_0 = log(1.0 / (SSB_virgin_adj / Recr_virgin_adj));
+ SRZ_0 = log(1.0 / (SSB_virgin_use / Recr_virgin_use));
steepness = SR_parm_work(2);
srz_min = SRZ_0 * (1.0 - steepness);
- SRZ_surv = mfexp((1. - pow((SSB_curr_adj / SSB_virgin_adj), SR_parm_work(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
+ SRZ_surv = mfexp((1. - pow((SSB_curr_adj / SSB_virgin_use), SR_parm_work(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
NewRecruits = SSB_curr_adj * SRZ_surv;
exp_rec(y, 1) = NewRecruits; // expected arithmetic mean recruitment
// SS_Label_43.3.7.1 Do variation in recruitment by adjusting survival
@@ -145,8 +145,8 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
Shepherd_c2 = pow(0.2, SR_parm_work(3));
Hupper = 1.0 / (5.0 * Shepherd_c2);
steepness = 0.2 + (SR_parm_work(2) - 0.2) / (0.8) * (Hupper - 0.2);
- temp = (SSB_curr_adj) / (SSB_virgin_adj);
- NewRecruits = (5. * steepness * Recr_virgin_adj * (1. - Shepherd_c2) * temp) /
+ temp = (SSB_curr_adj) / (SSB_virgin_use);
+ NewRecruits = (5. * steepness * Recr_virgin_use * (1. - Shepherd_c2) * temp) /
(1.0 - 5.0 * steepness * Shepherd_c2 + (5. * steepness - 1.) * pow(temp, Shepherd_c));
break;
}
@@ -156,11 +156,11 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
{
steepness = SR_parm_work(2);
dvariable RkrPower = SR_parm_work(3);
- temp = SSB_curr_adj / SSB_virgin_adj;
+ temp = SSB_curr_adj / SSB_virgin_use;
temp2 = posfun(1.0 - temp, 0.0000001, temp3);
temp = 1.0 - temp2; // Rick's new line to stabilize recruitment at R0 if B>B0
dvariable RkrTop = log(5.0 * steepness) * pow(temp2, RkrPower) / pow(0.8, RkrPower);
- NewRecruits = Recr_virgin_adj * temp * mfexp(RkrTop);
+ NewRecruits = Recr_virgin_use * temp * mfexp(RkrTop);
break;
}
@@ -169,7 +169,7 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
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
@@ -196,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)
{
@@ -228,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
@@ -268,7 +268,7 @@ 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 dvar_vector& SRparm,
- const prevariable& SSB_temp, const prevariable& RECR_temp, const prevariable& SPR_temp)
+ const prevariable& SSB_virgin_use, const prevariable& Recr_virgin_use, const prevariable& SPR_current)
{
RETURN_ARRAYS_INCREMENT();
dvar_vector Equil_Spawn_Recr_Calc(1, 2); // values to return 1 is B_equil, 2 is R_equil
@@ -284,7 +284,7 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
dvariable SRZ_surv;
steepness = SRparm(2); // 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
+ // SS_Label_44.1 calc equilibrium SpawnBio and Recruitment from input SPR_current, which is spawning biomass per recruit at some given F level
switch (SR_fxn)
{
case 1: // previous placement for B-H constrained
@@ -297,8 +297,8 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS_Label_44.1.2 Ricker
case 2: // Ricker
{
- B_equil = SSB_temp * (1. + (log(RECR_temp / SSB_temp) + log(SPR_temp)) / steepness);
- R_equil = RECR_temp * B_equil / SSB_temp * mfexp(steepness * (1. - B_equil / SSB_temp));
+ B_equil = SSB_virgin_use * (1. + (log(Recr_virgin_use / SSB_virgin_use) + log(SPR_current)) / steepness);
+ R_equil = Recr_virgin_use * B_equil / SSB_virgin_use * mfexp(steepness * (1. - B_equil / SSB_virgin_use));
break;
}
@@ -316,19 +316,19 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS3 previously used alternative formulation: R = A*S/(B+S)
// converting SS3 to align with WHAM
- // SPR_virgin = SSB_temp / RECR_temp; // this is already defined
+ // SPR_virgin = SSB_virgin_use / Recr_virgin_use; // this is already defined
alpha = 4.0 * steepness / (SPR_virgin_adj * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
// " h " << steepness << " derive " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " " << endl;
- // " R0 " << RECR_temp << " derive " << 1. / beta * (alpha - 1./SPR_virgin) << endl;
-// report5 <<" SSB_unf "<<SSB_temp<<" SPR_unf "<<SPR_virgin<<" steep: "<<steepness<<" R0: "<<RECR_temp << endl;
- report5 <<" derive_alpha "<<alpha<<" derive_beta "<<beta << endl;
+ // " R0 " << Recr_virgin_use << " derive " << 1. / beta * (alpha - 1./SPR_virgin) << endl;
+// report5 <<" SSB_unf "<<SSB_virgin_use<<" SPR_unf "<<SPR_virgin<<" steep: "<<steepness<<" R0: "<<Recr_virgin_use << endl;
+// report5 <<" derive_alpha "<<alpha<<" derive_beta "<<beta << endl;
// report5 << " deriv_h: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " derive_R0: " << 1. / beta * (alpha - (1. / SPR_virgin))<<endl;
- B_equil = (alpha * SPR_temp - 1.0) / beta;
+ B_equil = (alpha * SPR_current - 1.0) / beta;
B_equil = posfun(B_equil, 0.0001, temp);
R_equil = alpha * B_equil / (1.0 + beta * B_equil);
-// report5 << "SPR_input: " << SPR_temp << " B_equil: " << B_equil << " R_equil: "<<R_equil << endl<<endl;
+// report5 << "SPR_input: " << SPR_current << " B_equil: " << B_equil << " R_equil: "<<R_equil << endl<<endl;
break;
}
@@ -337,37 +337,37 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
{
dvariable alpha = mfexp(SRparm(3));
dvariable beta = mfexp(SRparm(4));
- B_equil = (alpha * SPR_temp - 1.0) / beta;
+ B_equil = (alpha * SPR_current - 1.0) / beta;
B_equil = posfun(B_equil, 0.0001, temp);
R_equil = alpha * B_equil / (1.0 + beta * B_equil);
-// report5<<SPR_temp<<" Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
+// report5<<SPR_current<<" Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
break;
}
// SS_Label_44.1.4 constant recruitment
case 4: // constant; no bias correction
{
- B_equil = SPR_temp * RECR_temp;
- R_equil = RECR_temp;
+ B_equil = SPR_current * Recr_virgin_use;
+ R_equil = Recr_virgin_use;
break;
}
// SS_Label_44.1.5 Hockey Stick
case 5: // hockey stick
{
- dvariable hockey_min = SRparm(3) * RECR_temp; // min recruitment level
- // temp=SSB_temp/R0*steepness; // spawners per recruit at inflection
- dvariable hockey_slope = (RECR_temp - hockey_min) / (steepness * SSB_temp); // slope of recruitment on spawners below the inflection
- B_equil = Join_Fxn(0.0 * SSB_temp / RECR_temp, SSB_temp / RECR_temp, SSB_temp / RECR_temp * steepness, SPR_temp, hockey_min / ((1. / SPR_temp) - hockey_slope), SPR_temp * RECR_temp);
- R_equil = Join_Fxn(0.0 * SSB_temp, SSB_temp, SSB_temp * steepness, B_equil, hockey_min + hockey_slope * B_equil, RECR_temp);
+ dvariable hockey_min = SRparm(3) * Recr_virgin_use; // min recruitment level
+ // temp=SSB_virgin_use/R0*steepness; // spawners per recruit at inflection
+ dvariable hockey_slope = (Recr_virgin_use - hockey_min) / (steepness * SSB_virgin_use); // slope of recruitment on spawners below the inflection
+ B_equil = Join_Fxn(0.0 * SSB_virgin_use / Recr_virgin_use, SSB_virgin_use / Recr_virgin_use, SSB_virgin_use / Recr_virgin_use * steepness, SPR_current, hockey_min / ((1. / SPR_current) - hockey_slope), SPR_current * Recr_virgin_use);
+ R_equil = Join_Fxn(0.0 * SSB_virgin_use, SSB_virgin_use, SSB_virgin_use * steepness, B_equil, hockey_min + hockey_slope * B_equil, Recr_virgin_use);
break;
}
// SS_Label_44.1.7 3 parameter survival based
case 7: // survival
{
- SRZ_0 = log(1.0 / (SSB_temp / RECR_temp));
+ SRZ_0 = log(1.0 / (SSB_virgin_use / Recr_virgin_use));
srz_min = SRZ_0 * (1.0 - steepness);
- B_equil = SSB_temp * (1. - (log(1. / SPR_temp) - SRZ_0) / pow((srz_min - SRZ_0), (1. / SRparm(3))));
- SRZ_surv = mfexp((1. - pow((B_equil / SSB_temp), SRparm(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
+ B_equil = SSB_virgin_use * (1. - (log(1. / SPR_current) - SRZ_0) / pow((srz_min - SRZ_0), (1. / SRparm(3))));
+ SRZ_surv = mfexp((1. - pow((B_equil / SSB_virgin_use), SRparm(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
R_equil = B_equil * SRZ_surv;
break;
}
@@ -389,11 +389,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
Shepherd_c2 = pow(0.2, SRparm(3));
Hupper = 1.0 / (5.0 * Shepherd_c2);
steepness = 0.2 + (SRparm(2) - 0.2) / (0.8) * (Hupper - 0.2);
- Shep_top = 5.0 * steepness * (1.0 - Shepherd_c2) * (SPR_temp * RECR_temp) / SSB_temp - (1.0 - 5.0 * steepness * Shepherd_c2);
+ Shep_top = 5.0 * steepness * (1.0 - Shepherd_c2) * (SPR_current * Recr_virgin_use) / SSB_virgin_use - (1.0 - 5.0 * steepness * Shepherd_c2);
Shep_bot = 5.0 * steepness - 1.0;
Shep_top2 = posfun(Shep_top, 0.001, temp);
- R_equil = (SSB_temp / SPR_temp) * pow((Shep_top2 / Shep_bot), (1.0 / SRparm(3)));
- B_equil = R_equil * SPR_temp;
+ R_equil = (SSB_virgin_use / SPR_current) * pow((Shep_top2 / Shep_bot), (1.0 / SRparm(3)));
+ B_equil = R_equil * SPR_current;
break;
}
@@ -402,11 +402,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
{
steepness = SRparm(2);
dvariable RkrPower = SRparm(3);
- temp = SSB_temp / (SPR_temp * RECR_temp);
+ temp = SSB_virgin_use / (SPR_current * Recr_virgin_use);
dvariable RkrTop = pow(0.8, RkrPower) * log(temp) / log(5.0 * steepness);
RkrTop = posfun(RkrTop, 0.000001, CrashPen);
- R_equil = temp * RECR_temp * (1.0 - pow(RkrTop, 1.0 / RkrPower));
- B_equil = R_equil * SPR_temp;
+ R_equil = temp * Recr_virgin_use * (1.0 - pow(RkrTop, 1.0 / RkrPower));
+ B_equil = R_equil * SPR_current;
break;
}
@@ -428,11 +428,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
Hupper=1.0/(5.0*Shepherd_c2);
steepness=0.20001+((0.8)/(1.0+exp(-SRparm2))-0.2)/(0.8)*(Hupper-0.2);
// steep2=0.20001+(steepness-0.2)/(0.8)*(Hupper-0.2);
- Shep_top=5.0*steepness*(1.0-Shepherd_c2)*(SPR_temp*RECR_temp)/SSB_temp-(1.0-5.0*steepness*Shepherd_c2);
+ Shep_top=5.0*steepness*(1.0-Shepherd_c2)*(SPR_current*Recr_virgin_use)/SSB_virgin_use-(1.0-5.0*steepness*Shepherd_c2);
Shep_bot=5.0*steepness-1.0;
Shep_top2=posfun(Shep_top,0.001,temp);
- R_equil=(SSB_temp/SPR_temp) * pow((Shep_top2/Shep_bot),(1.0/Shepherd_c));
- B_equil=R_equil*SPR_temp;
+ R_equil=(SSB_virgin_use/SPR_current) * pow((Shep_top2/Shep_bot),(1.0/Shepherd_c));
+ B_equil=R_equil*SPR_current;
break;
}
@@ -448,11 +448,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// if (Recs < 0) Rec2 = 0; else Rec2 = Recs;
steepness = 0.2 + (10.0 - 0.2)/(1+exp(-SR_parm_work(2)));
dvariable RkrPower=exp(SR_parm_work(3));
- temp=SSB_virgin/(SPR_temp*RECR_temp);
+ temp=SSB_virgin/(SPR_current*Recr_virgin_use);
dvariable RkrTop = pow(0.8,RkrPower)*log(temp)/log(5.0*steepness);
RkrTop = posfun(RkrTop,0.000001,CrashPen);
- R_equil = temp *RECR_temp * (1.0 - pow(RkrTop,1.0/RkrPower));
- B_equil=R_equil*SPR_temp;
+ R_equil = temp *Recr_virgin_use * (1.0 - pow(RkrTop,1.0/RkrPower));
+ B_equil=R_equil*SPR_current;
break;
}
*/
From d351e200f3961584eef2946b0b1792e1245d1f0c Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Mon, 16 Sep 2024 18:22:06 -0700
Subject: [PATCH 24/29] refactor SPR to SSBpR; clean-up reference to timevary
biology
---
SS_benchfore.tpl | 73 ++++++++++++++++++++++++---------------------
SS_param.tpl | 6 ++--
SS_popdyn.tpl | 18 ++++++-----
SS_recruit.tpl | 54 ++++++++++++++++-----------------
SS_write_report.tpl | 14 ++++-----
5 files changed, 85 insertions(+), 80 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 25038380..79a792fa 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -536,7 +536,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
dvariable last_F1;
dvariable Closer;
dvariable Vbio1_unfished;
- dvariable SPR_unf;
+ dvariable SSBpR_unf;
dvariable Vbio_MSY;
dvariable Vbio1_MSY;
dvariable junk;
@@ -749,7 +749,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
{
Recr_unf = Recr_virgin;
SSB_unf = SSB_virgin;
- SPR_unf = SSB_unf / Recr_unf;
+ SSBpR_unf = SSB_unf / Recr_unf;
}
else // use benchmark biology in the spawner-recruitment R0,h calculations
{
@@ -757,18 +757,21 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Recr_unf = mfexp(SR_parm_work(1));
Do_Equil_Calc(Recr_unf);
SSB_unf = SSB_equil; // equilibrium unfished SSB using the benchmark averaged Recr_unf and benchmark averaged biology
- SPR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
+ SSBpR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
+ SSBpR_virgin_adj = SSB_unf / Recr_unf; // update
}
if (show_MSY == 1)
{
report5 << "SR_parms for benchmark: " << SR_parm_work << endl
<< "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl;
if ( SR_update_SPR0 == 1) report5 << "SPR0 for equilibrium spawner-recruit based on benchmark biology, not virgin biology" << endl;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SPR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Virgin SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SPR_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Virgin SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SSBpR_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
+
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Benchmark SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " " << SSBpR_unf << " equil: " << Equ_SpawnRecr_Result << endl;
+ report5 << "Repro_output_by_age_for_morph_1: " << fec(1) << endl;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Benchmark SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " " << SPR_unf << " equil: " << Equ_SpawnRecr_Result << endl;
}
SR_parm_work(N_SRparm2 + 1) = SSB_unf;
Mgmt_quant(1) = SSB_unf;
@@ -803,7 +806,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SPR_target100 = SPR_target * 100.;
Do_Equil_Calc(equ_Recr);
- SPR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
+ SSBpR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
Vbio1_unfished = smrybio; // gets value from equil_calc
if (show_MSY == 1)
{
@@ -849,7 +852,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Fishon = 1;
Do_Equil_Calc(equ_Recr);
- yld1(ii) = 100. * SSB_equil / SPR_unf; // spawning potential ratio
+ yld1(ii) = 100. * SSB_equil / SSBpR_unf; // spawning potential ratio
}
SPR_actual = yld1(1); // spawning potential ratio
@@ -918,8 +921,8 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
// SPAWN-RECR: calc equil spawn-recr in YPR; need to make this area-specific
- SPR_temp = SSB_equil; // based on most recent call to Do_Equil_Calc
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSBpR_temp = SSB_equil; // based on most recent call to Do_Equil_Calc
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Bspr = Equ_SpawnRecr_Result(1);
Bspr_rec = Equ_SpawnRecr_Result(2);
@@ -1002,7 +1005,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_unf << " " << 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;
@@ -1033,8 +1036,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, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Btgt = Equ_SpawnRecr_Result(1);
Btgt_Rec = Equ_SpawnRecr_Result(2);
YPR_Btgt_enc = YPR_enc; // total encountered yield per recruit
@@ -1045,7 +1048,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_unf;
+ SPR_Btgt = SSB_equil / SSBpR_unf;
Vbio_Btgt = totbio;
Vbio1_Btgt = smrybio;
Mgmt_quant(7) = equ_F_std;
@@ -1060,9 +1063,12 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
// ******************************************************
if (show_MSY == 1)
{
- report5 << "#" << endl
- << "Find_target_SSB/Bzero; where Bzero is for Bmark years, not Virgin" << endl
- << "Iter Fmult ann_F SPR Catch SSB Recruits SSB/Bzero Tot_catch";
+ report5 << "#" << endl;
+ if (SR_update_SPR0 == 0) // use virgin biology for the spawner-recruitment R0,h calculations
+ {report5 << "Find_target_SSB/Bzero; where Bzero is Virgin SSB:" << SSB_unf << " where SSBpR_unf = " << SSBpR_unf << endl;}
+ else
+ {report5 << "Find_target_SSB/Bzero; where Bzero is for Bmark biology and updated SPR0: " << SSB_unf << " where SSBpR_unf = " << SSBpR_unf << endl;}
+ report5 << "Iter Fmult ann_F SPR Catch SSB Recruits SSB/Bzero Tot_catch";
for (p = 1; p <= pop; p++)
for (gp = 1; gp <= N_GP; gp++)
{
@@ -1088,8 +1094,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Nloops = 28;
}
- // Btgttgt=BTGT_target*SSB_virgin; // this is relative to virgin, not to the average biology from benchmark years
- Btgttgt = BTGT_target * SSB_unf; // now relative to Bmark
+ Btgttgt = BTGT_target * SSB_unf;
for (j = 0; j <= Nloops; j++) // loop find Btarget
{
@@ -1124,12 +1129,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_unf;
+ Do_Equil_Calc(equ_Recr); // where equ_Recr=1.0, so returned SSB_equil is in units of SSB/R,
+ SSBpR_temp = SSB_equil;
+ SPR_Btgt = SSBpR_temp / SSBpR_unf; // where SSBpR_unf = SSB_unf / Recr_unf so units of SSB/R; so result is SPR_Btgt = (fished SSB/R) / (unfished SSB/R)
// SPAWN-RECR: calc equil spawn-recr for Btarget calcs; need to make area-specific
- SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5<<SPR_temp<<" " << Equ_SpawnRecr_Result<<endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
yld1(ii) = Equ_SpawnRecr_Result(1);
}
@@ -1169,6 +1173,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
report5 << " " << SSB_equil_pop_gp(p, gp) * Equ_SpawnRecr_Result(2);
}
report5 << endl;
+// " SSB_unf " << SSB_unf << " recr " <<Recr_unf<< " SSB_equil " << SSB_equil << " ssbpr_virginadj: " << SSBpR_virgin_adj << " Btgt " << Btgt << endl;
}
} // end search loop
@@ -1347,9 +1352,9 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Do_Equil_Calc(equ_Recr);
// SPAWN-RECR: calc spawn-recr for MSY calcs; need to make area-specific
- MSY_SPR = SSB_equil / SPR_unf;
- SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ MSY_SPR = SSB_equil / SSBpR_unf;
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
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;
@@ -1440,9 +1445,9 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
Do_Equil_Calc(equ_Recr);
// SPAWN-RECR: calc spawn-recr for MSY calcs; need to make area-specific
- MSY_SPR = SSB_equil / SPR_unf;
- SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ MSY_SPR = SSB_equil / SSBpR_unf;
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Bmsy = Equ_SpawnRecr_Result(1); // MSY is directly estimated
Recr_msy = Equ_SpawnRecr_Result(2);
Profit = (PricePerF * YPR_val_vec) * Recr_msy - Cost;
@@ -1733,10 +1738,10 @@ 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_unf;
+ 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, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
yld1(ii) = Equ_SpawnRecr_Result(1);
}
diff --git a/SS_param.tpl b/SS_param.tpl
index 3d000222..6a454579 100644
--- a/SS_param.tpl
+++ b/SS_param.tpl
@@ -165,8 +165,8 @@ PARAMETER_SECTION
number two_sigmaRsq;
number half_sigmaRsq;
number sigmaR;
- number SPR_virgin;
- number SPR_virgin_adj;
+ number SSBpR_virgin;
+ number SSBpR_virgin_adj;
number regime_change;
number rho;
number dirichlet_Parm;
@@ -249,7 +249,7 @@ PARAMETER_SECTION
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
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index a0297e3c..d475456a 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -343,8 +343,8 @@ FUNCTION void get_initial_conditions()
{
equ_Recr = 1.0;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation. Returns SPR because R = 1.0
- SPR_virgin = SSB_equil; // spawners per recruit. Needed for Sr_fxn = 10
- SPR_virgin_adj = SSB_equil; // also needed for Sr_fxn 10. Will get revised in benchmark to use averaged biology if requested.
+ SSBpR_virgin = SSB_equil; // spawners per recruit. Needed for Sr_fxn = 10
+ SSBpR_virgin_adj = SSB_equil; // also needed for Sr_fxn 10. Will get revised in benchmark to use averaged biology if requested.
if(SR_fxn == 10) // B-H with a,b
{
// WHAM based on R = A*S/(1+B*S)
@@ -355,8 +355,8 @@ FUNCTION void get_initial_conditions()
alpha = mfexp(SR_parm(3));
beta = mfexp(SR_parm(4));
- steepness = alpha * SPR_virgin_adj / (4. + alpha * SPR_virgin_adj);
- Recr_virgin = 1. / beta * (alpha - (1. / SPR_virgin_adj));
+ steepness = alpha * SSBpR_virgin_adj / (4. + alpha * SSBpR_virgin_adj);
+ Recr_virgin = 1. / beta * (alpha - (1. / SSBpR_virgin_adj));
// warning << " before AB_calcs " << "parm " << SR_parm(1) << " calc " << log(Recr_virgin) << endl;
SR_parm(1) = log(Recr_virgin);
SR_parm(2) = steepness;
@@ -380,7 +380,7 @@ FUNCTION void get_initial_conditions()
exp_rec(eq_yr, 4) = Recr_virgin;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
SSB_virgin = SSB_equil;
-// SPR_virgin = SSB_equil / Recr_virgin; // spawners per recruit already calculated
+// SSBpR_virgin = SSB_equil / Recr_virgin; // spawners per recruit already calculated
if(Do_Benchmark==0) // assign values that would be created in benchmark section
{
Mgmt_quant(1) = SSB_virgin;
@@ -533,10 +533,10 @@ FUNCTION void get_initial_conditions()
Do_Equil_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
+ 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(SR_parm_work, SSB_virgin, Recr_virgin, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_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;
@@ -1040,6 +1040,7 @@ FUNCTION void get_time_series()
bio_yr = y;
Do_Equil_Calc(R0_use); // call function to do equilibrium calculation with current year's biology and adjusted R0
SSB_use = SSB_equil;
+ SSBpR_virgin_adj = SSB_use / R0_use; // update
if (fishery_on_off == 1)
{
Fishon = 1;
@@ -1498,6 +1499,7 @@ FUNCTION void get_time_series()
bio_yr = y;
Do_Equil_Calc(R0_use); // call function to do equilibrium calculation
SSB_use = SSB_equil;
+ SSBpR_virgin_adj = SSB_use / R0_use; // update
if (fishery_on_off == 1)
{
Fishon = 1;
diff --git a/SS_recruit.tpl b/SS_recruit.tpl
index 5766efce..645dbcf3 100644
--- a/SS_recruit.tpl
+++ b/SS_recruit.tpl
@@ -268,7 +268,7 @@ 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 dvar_vector& SRparm,
- const prevariable& SSB_virgin_use, const prevariable& Recr_virgin_use, const prevariable& SPR_current)
+ 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
@@ -284,7 +284,7 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
dvariable SRZ_surv;
steepness = SRparm(2); // common usage but some different
- // SS_Label_44.1 calc equilibrium SpawnBio and Recruitment from input SPR_current, which is spawning biomass per recruit at some given F level
+ // 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
@@ -297,7 +297,7 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS_Label_44.1.2 Ricker
case 2: // Ricker
{
- B_equil = SSB_virgin_use * (1. + (log(Recr_virgin_use / SSB_virgin_use) + log(SPR_current)) / steepness);
+ 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;
@@ -316,19 +316,17 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS3 previously used alternative formulation: R = A*S/(B+S)
// converting SS3 to align with WHAM
- // SPR_virgin = SSB_virgin_use / Recr_virgin_use; // this is already defined
- alpha = 4.0 * steepness / (SPR_virgin_adj * (1. - steepness));
- beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
-
- // " h " << steepness << " derive " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " " << endl;
- // " R0 " << Recr_virgin_use << " derive " << 1. / beta * (alpha - 1./SPR_virgin) << endl;
-// report5 <<" SSB_unf "<<SSB_virgin_use<<" SPR_unf "<<SPR_virgin<<" steep: "<<steepness<<" R0: "<<Recr_virgin_use << endl;
+ alpha = 4.0 * steepness / (SSBpR_virgin_adj * (1. - steepness));
+ beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin_use);
+ // " h " << steepness << " derive " << alpha * SSBpR_virgin / (4. + alpha * SSBpR_virgin) << " " << endl;
+ // " R0 " << Recr_virgin_use << " derive " << 1. / beta * (alpha - 1./SSBpR_virgin) << endl;
+// report5 <<" SSB_unf "<<SSB_virgin_use<<" SSBpR_unf "<<SSBpR_virgin<<" steep: "<<steepness<<" R0: "<<Recr_virgin_use << endl;
// report5 <<" derive_alpha "<<alpha<<" derive_beta "<<beta << endl;
-// report5 << " deriv_h: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " derive_R0: " << 1. / beta * (alpha - (1. / SPR_virgin))<<endl;
- B_equil = (alpha * SPR_current - 1.0) / beta;
+// report5 << " deriv_h: " << alpha * SSBpR_virgin / (4. + alpha * SSBpR_virgin) << " derive_R0: " << 1. / beta * (alpha - (1. / SSBpR_virgin))<<endl;
+ B_equil = (alpha * SSBpR_current - 1.0) / beta;
B_equil = posfun(B_equil, 0.0001, temp);
R_equil = alpha * B_equil / (1.0 + beta * B_equil);
-// report5 << "SPR_input: " << SPR_current << " B_equil: " << B_equil << " R_equil: "<<R_equil << endl<<endl;
+// report5 << "SPR_input: " << SSBpR_current << " B_equil: " << B_equil << " R_equil: "<<R_equil << endl<<endl;
break;
}
@@ -337,17 +335,17 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
{
dvariable alpha = mfexp(SRparm(3));
dvariable beta = mfexp(SRparm(4));
- B_equil = (alpha * SPR_current - 1.0) / beta;
+ B_equil = (alpha * SSBpR_current - 1.0) / beta;
B_equil = posfun(B_equil, 0.0001, temp);
R_equil = alpha * B_equil / (1.0 + beta * B_equil);
-// report5<<SPR_current<<" Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
+// report5<<SSBpR_current<<" Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
break;
}
// SS_Label_44.1.4 constant recruitment
case 4: // constant; no bias correction
{
- B_equil = SPR_current * Recr_virgin_use;
+ B_equil = SSBpR_current * Recr_virgin_use;
R_equil = Recr_virgin_use;
break;
}
@@ -357,7 +355,7 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
dvariable hockey_min = SRparm(3) * Recr_virgin_use; // min recruitment level
// temp=SSB_virgin_use/R0*steepness; // spawners per recruit at inflection
dvariable hockey_slope = (Recr_virgin_use - hockey_min) / (steepness * SSB_virgin_use); // slope of recruitment on spawners below the inflection
- B_equil = Join_Fxn(0.0 * SSB_virgin_use / Recr_virgin_use, SSB_virgin_use / Recr_virgin_use, SSB_virgin_use / Recr_virgin_use * steepness, SPR_current, hockey_min / ((1. / SPR_current) - hockey_slope), SPR_current * Recr_virgin_use);
+ B_equil = Join_Fxn(0.0 * SSB_virgin_use / Recr_virgin_use, SSB_virgin_use / Recr_virgin_use, SSB_virgin_use / Recr_virgin_use * steepness, SSBpR_current, hockey_min / ((1. / SSBpR_current) - hockey_slope), SSBpR_current * Recr_virgin_use);
R_equil = Join_Fxn(0.0 * SSB_virgin_use, SSB_virgin_use, SSB_virgin_use * steepness, B_equil, hockey_min + hockey_slope * B_equil, Recr_virgin_use);
break;
}
@@ -366,7 +364,7 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
{
SRZ_0 = log(1.0 / (SSB_virgin_use / Recr_virgin_use));
srz_min = SRZ_0 * (1.0 - steepness);
- B_equil = SSB_virgin_use * (1. - (log(1. / SPR_current) - SRZ_0) / pow((srz_min - SRZ_0), (1. / SRparm(3))));
+ B_equil = SSB_virgin_use * (1. - (log(1. / SSBpR_current) - SRZ_0) / pow((srz_min - SRZ_0), (1. / SRparm(3))));
SRZ_surv = mfexp((1. - pow((B_equil / SSB_virgin_use), SRparm(3))) * (srz_min - SRZ_0) + SRZ_0); // survival
R_equil = B_equil * SRZ_surv;
break;
@@ -389,11 +387,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
Shepherd_c2 = pow(0.2, SRparm(3));
Hupper = 1.0 / (5.0 * Shepherd_c2);
steepness = 0.2 + (SRparm(2) - 0.2) / (0.8) * (Hupper - 0.2);
- Shep_top = 5.0 * steepness * (1.0 - Shepherd_c2) * (SPR_current * Recr_virgin_use) / SSB_virgin_use - (1.0 - 5.0 * steepness * Shepherd_c2);
+ Shep_top = 5.0 * steepness * (1.0 - Shepherd_c2) * (SSBpR_current * Recr_virgin_use) / SSB_virgin_use - (1.0 - 5.0 * steepness * Shepherd_c2);
Shep_bot = 5.0 * steepness - 1.0;
Shep_top2 = posfun(Shep_top, 0.001, temp);
- R_equil = (SSB_virgin_use / SPR_current) * pow((Shep_top2 / Shep_bot), (1.0 / SRparm(3)));
- B_equil = R_equil * SPR_current;
+ R_equil = (SSB_virgin_use / SSBpR_current) * pow((Shep_top2 / Shep_bot), (1.0 / SRparm(3)));
+ B_equil = R_equil * SSBpR_current;
break;
}
@@ -402,11 +400,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
{
steepness = SRparm(2);
dvariable RkrPower = SRparm(3);
- temp = SSB_virgin_use / (SPR_current * Recr_virgin_use);
+ temp = SSB_virgin_use / (SSBpR_current * Recr_virgin_use);
dvariable RkrTop = pow(0.8, RkrPower) * log(temp) / log(5.0 * steepness);
RkrTop = posfun(RkrTop, 0.000001, CrashPen);
R_equil = temp * Recr_virgin_use * (1.0 - pow(RkrTop, 1.0 / RkrPower));
- B_equil = R_equil * SPR_current;
+ B_equil = R_equil * SSBpR_current;
break;
}
@@ -428,11 +426,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
Hupper=1.0/(5.0*Shepherd_c2);
steepness=0.20001+((0.8)/(1.0+exp(-SRparm2))-0.2)/(0.8)*(Hupper-0.2);
// steep2=0.20001+(steepness-0.2)/(0.8)*(Hupper-0.2);
- Shep_top=5.0*steepness*(1.0-Shepherd_c2)*(SPR_current*Recr_virgin_use)/SSB_virgin_use-(1.0-5.0*steepness*Shepherd_c2);
+ Shep_top=5.0*steepness*(1.0-Shepherd_c2)*(SSBpR_current*Recr_virgin_use)/SSB_virgin_use-(1.0-5.0*steepness*Shepherd_c2);
Shep_bot=5.0*steepness-1.0;
Shep_top2=posfun(Shep_top,0.001,temp);
- R_equil=(SSB_virgin_use/SPR_current) * pow((Shep_top2/Shep_bot),(1.0/Shepherd_c));
- B_equil=R_equil*SPR_current;
+ R_equil=(SSB_virgin_use/SSBpR_current) * pow((Shep_top2/Shep_bot),(1.0/Shepherd_c));
+ B_equil=R_equil*SSBpR_current;
break;
}
@@ -448,11 +446,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// if (Recs < 0) Rec2 = 0; else Rec2 = Recs;
steepness = 0.2 + (10.0 - 0.2)/(1+exp(-SR_parm_work(2)));
dvariable RkrPower=exp(SR_parm_work(3));
- temp=SSB_virgin/(SPR_current*Recr_virgin_use);
+ temp=SSB_virgin/(SSBpR_current*Recr_virgin_use);
dvariable RkrTop = pow(0.8,RkrPower)*log(temp)/log(5.0*steepness);
RkrTop = posfun(RkrTop,0.000001,CrashPen);
R_equil = temp *Recr_virgin_use * (1.0 - pow(RkrTop,1.0/RkrPower));
- B_equil=R_equil*SPR_current;
+ B_equil=R_equil*SSBpR_current;
break;
}
*/
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index 25e2f936..4a48fa0f 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1803,13 +1803,13 @@ FUNCTION void write_bigoutput()
SS2out << "# " <<endl;
SS2out << "sigmaR: "<< sigmaR << endl;
- SS2out << "SPR0_(virgin): " << SPR_virgin << " #_uses_biology_from_start_year: " << styr <<endl;
+ SS2out << "SPR0_(virgin): " << SSBpR_virgin << " #_uses_biology_from_start_year: " << styr <<endl;
switch (SR_fxn)
{
case 3: // Beverton-Holt with steepness
{
- alpha = 4.0 * steepness / (SPR_virgin_adj * (1. - steepness));
+ alpha = 4.0 * steepness / (SSBpR_virgin_adj * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1. - steepness) * SSB_virgin);
SS2out << "Ln(R0): " << SR_parm(1) << endl << "R0: " << mfexp(SR_parm(1)) << endl;
SS2out << "steepness: " << steepness << endl;
@@ -1821,8 +1821,8 @@ FUNCTION void write_bigoutput()
{
SS2out << "Ln(alpha): " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
SS2out << "Ln(beta): " << SR_parm(4) << " beta " << mfexp(SR_parm(4)) << endl;
- SS2out << "ln(R0)_derived: " << log( 1. / beta * (alpha - (1. / SPR_virgin_adj))) << endl; // virgin R0
- SS2out << "steepness_derived: " << alpha * SPR_virgin_adj / (4. + alpha * SPR_virgin_adj) << endl; // steepness virgin
+ SS2out << "ln(R0)_derived: " << log( 1. / beta * (alpha - (1. / SSBpR_virgin_adj))) << endl; // virgin R0
+ SS2out << "steepness_derived: " << alpha * SSBpR_virgin_adj / (4. + alpha * SSBpR_virgin_adj) << endl; // steepness virgin
break;
}
case 8:
@@ -1854,7 +1854,7 @@ FUNCTION void write_bigoutput()
{ if (Bmark_Yr(k+1) > Bmark_Yr(k))
{SS2out << "#_range_of_years_is_averaged,_so_reduces_standard_error_of_result;_do_this_only_when_timevarying_makes_necessary: " << k << " "<< k+1 << endl;}
}
- SS2out << "SPR_unfished_benchmark: " << Mgmt_quant(1) / Mgmt_quant(4) << " #_based_on_averaging_biology_over_benchmark_year_range " << endl;
+ SS2out << "SSBpR_unfished_benchmark: " << Mgmt_quant(1) / Mgmt_quant(4) << " #_based_on_averaging_biology_over_benchmark_year_range " << endl;
SS2out << "Bmsy/Bzero: "<< Bmsy / SSB_virgin << " # using styr bio for Bzero" << endl;
SS2out << "Bmsy/Bunf: "<< Bmsy / Mgmt_quant(1) << " # using MSY's averaged bio for Bunf" << endl;
@@ -4862,8 +4862,8 @@ FUNCTION void SPR_profile()
Do_Equil_Calc(equ_Recr);
// SPAWN-RECR: calc equil spawn-recr in the SPR loop
- SPR_temp = SSB_equil;
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ SSBpR_temp = SSB_equil;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Btgt_prof = Equ_SpawnRecr_Result(1);
Btgt_prof_rec = Equ_SpawnRecr_Result(2);
if (Btgt_prof < 0.001 || Btgt_prof_rec < 0.001)
From 7a7ebf7a2a60cbfd90ddeb5f155fa27bb68eaf42 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Fri, 20 Sep 2024 09:58:02 -0700
Subject: [PATCH 25/29] create new switch to control updating of SSBpR0
---
SS_benchfore.tpl | 12 +++++-----
SS_popdyn.tpl | 8 +++----
SS_readcontrol_330.tpl | 53 +++++++++++++++++++++++++++++++++++++-----
SS_write_ssnew.tpl | 6 ++++-
4 files changed, 62 insertions(+), 17 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 79a792fa..d2cee309 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -745,7 +745,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
}
- if (SR_update_SPR0 == 0) // use virgin biology for the spawner-recruitment R0,h calculations
+ if (SR_update_SSBpR0_bmark == 0) // use virgin biology for the spawner-recruitment R0,h calculations in bmark
{
Recr_unf = Recr_virgin;
SSB_unf = SSB_virgin;
@@ -764,7 +764,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
{
report5 << "SR_parms for benchmark: " << SR_parm_work << endl
<< "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl;
- if ( SR_update_SPR0 == 1) report5 << "SPR0 for equilibrium spawner-recruit based on benchmark biology, not virgin biology" << endl;
+ if ( SR_update_SSBpR0_bmark == 1) report5 << "SPR0 for equilibrium spawner-recruit based on benchmark biology, not virgin biology" << endl;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
report5 << " Virgin SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SSBpR_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
@@ -1064,7 +1064,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
if (show_MSY == 1)
{
report5 << "#" << endl;
- if (SR_update_SPR0 == 0) // use virgin biology for the spawner-recruitment R0,h calculations
+ if (SR_update_SSBpR0_bmark == 0) // use virgin biology for the spawner-recruitment R0,h calculations
{report5 << "Find_target_SSB/Bzero; where Bzero is Virgin SSB:" << SSB_unf << " where SSBpR_unf = " << SSBpR_unf << endl;}
else
{report5 << "Find_target_SSB/Bzero; where Bzero is for Bmark biology and updated SPR0: " << SSB_unf << " where SSBpR_unf = " << SSBpR_unf << endl;}
@@ -1945,7 +1945,7 @@ FUNCTION void Get_Forecast()
dvariable OFL_catch;
dvariable Fcast_Crash;
dvariable totcatch;
- dvariable R0_use; // annually updated variable if SR_update_SPR0 == 1
+ dvariable R0_use; // annually updated value if SR_update_SSBpR0_timeseries == 1
dvariable SSB_use; // selected version of SSB that gets passes to Spawn_Recr
dvar_matrix catage_w(1, gmorph, 0, nages);
dvar_vector tempcatch(1, Nfleet);
@@ -2592,7 +2592,7 @@ FUNCTION void Get_Forecast()
}
// SPAWN-RECR: get recruitment in forecast; needs to be area-specific
// SR_fxn
- if (SR_update_SPR0 == 0) // R0 is not time-varying
+ if (SR_update_SSBpR0_timeseries == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
@@ -3230,7 +3230,7 @@ FUNCTION void Get_Forecast()
// 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-specific
// SR_fxn
- if (SR_update_SPR0 == 0) // R0 is not time-varying
+ if (SR_update_SSBpR0_timeseries == 0) // R0 is not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index d475456a..e5fda47c 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -697,7 +697,7 @@ FUNCTION void get_time_series()
dvariable crashtemp1;
dvariable interim_tot_catch;
dvariable Z_adjuster;
- dvariable R0_use; // annually updated variable if SR_update_SPR0 == 1; gets passed to Spawn_Recr() function
+ dvariable R0_use; // annually updated variable if SR_update_SSBpR0_timeseries == 1; gets passed to Spawn_Recr() function
dvariable SSB_use;
if (Do_Morphcomp > 0)
@@ -1026,7 +1026,7 @@ FUNCTION void get_time_series()
// SS_Label_Info_24.2.3 #Get the total recruitment produced by this spawning biomass at the beginning of the season
// SPAWN-RECR: calc recruitment in time series; need to make this area-specific
// SR_Fxn relevant keyword
- if (SR_update_SPR0 == 0) // SRparm are not time-varying
+ if (SR_update_SSBpR0_timeseries == 0) // SRparm are not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
@@ -1486,12 +1486,12 @@ FUNCTION void get_time_series()
// 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
// SR_fxn
- if (SR_update_SPR0 == 0) // SR parms are not time-varying
+ if (SR_update_SSBpR0_timeseries == 0) // SR parms are not time-varying
{
R0_use = Recr_virgin;
SSB_use = SSB_virgin;
}
- else // update SSB_use and R0_use where will update SPR0 inside the Spawn_recr fxn
+ else // update SSB_use and R0_use where will update SSBpR0 inside the Spawn_recr fxn
{
R0_use = mfexp(SR_parm_work(1)); // check to be sure this works when R0 is derived from B-H with alpha, beta parameters
Fishon = 0;
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index 905f9a81..83f2632c 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -1943,8 +1943,21 @@
!!echoinput<<SR_fxn<<" #_SR_function: 1=NA; 2=Ricker(2 parms); 3=BevHolt(2); 4=SCAA(2); 5=Hockey(3); 6=B-H_flattop(2); 7=Survival(3); 8=Shepherd(3); 9=Ricker_Power(3); 10=B-H_a,b(4)"<<endl;
init_int init_equ_steepness;
!!echoinput<<init_equ_steepness<<" # 0/1 to use steepness in initial equ recruitment calculation"<<endl;
- init_int sigmaR_dendep;
-!! echoinput<<sigmaR_dendep<<" # future feature: 0/1 to make realized sigmaR a function of SR curvature"<<endl;
+ init_int SR_update_SSBpR0_rd;
+// SR_update_SSBpR0_rd values
+// 1 best: update SSBpR0 for benchmark and for time series only if SRparm R0 or h (not regime) is set to have time-varying property
+// 2 incorrect, relic: always update SSBpR0 for benchmark's use of spawner-recruitment (old, incorrect SS3 approach), but only for the time series if there is a timevary SR parm
+// 3 option: do not update SSBpR0 (do keep start year SPR0), even if R0 or h is set to have time-varying property
+ int SR_update_SSBpR0_bmark
+ int SR_update_SSBpR0_timeseries
+ LOCAL_CALCS
+ SR_update_SSBpR0_bmark = 0;
+ SR_update_SSBpR0_timeseries = 0;
+ echoinput<<SR_update_SSBpR0_rd<<" # controls effect of time-varying biology on spawner-recruitment updating" << endl;
+ END_CALCS
+
+// echoinput<<sigmaR_dendep<<" # future feature: 0/1 to make realized sigmaR a function of SR curvature"<<endl;
+
ivector N_SRparm(1,10)
!!N_SRparm.fill("{0,2,2,2,3,2,3,3,3,4}");
int N_SRparm2
@@ -1964,7 +1977,7 @@
int timevary_parm_cnt_SR;
ivector timevary_SRparm(styr-3,YrMax+1);
ivector SR_parm_timevary(1,N_SRparm2);
- int SR_update_SPR0 // 0/1 flag to control updating of SPR0 for timevary biology
+ int SR_update_parm // 0/1 flag to control updating of SPR0 for timevary biology
LOCAL_CALCS
// clang-format on
@@ -1976,7 +1989,7 @@
SR_parm_timevary.initialize();
SR_env_link = 0;
SR_env_target = 0;
- SR_update_SPR0 = 0;
+ SR_update_parm = 0;
//#_SR_function: 1=null; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=Survival_3Parm; 10=B-H with a,b "<<endl;
if (SR_fxn == 10)
@@ -2129,8 +2142,8 @@
timevary_SRparm(y) = timevary_pass(y);
} // year vector for this category og MGparm
// special consideration for impact on spawner-recruitment SPR0 calculations
- if (SR_fxn == 10 && (SR_parm_timevary(3) > 0 || SR_parm_timevary(4) > 0) ) {SR_update_SPR0 = 1;} // alpha or beta is time-varying
- else if ((SR_parm_timevary(1) > 0 || SR_parm_timevary(2) > 0) ) {SR_update_SPR0 = 1;} // R0 or steepness is time-varying
+ if (SR_fxn == 10 && (SR_parm_timevary(3) > 0 || SR_parm_timevary(4) > 0) ) {SR_update_parm = 1;} // alpha or beta is time-varying
+ else if ((SR_parm_timevary(1) > 0 || SR_parm_timevary(2) > 0) ) {SR_update_parm = 1;} // R0 or steepness is time-varying
}
}
N_SRparm3 = N_SRparm2;
@@ -2142,6 +2155,34 @@
echoinput << " SR timevary_parm_cnt start and end " << timevary_parm_start_SR << " " << timevary_parm_cnt_SR << endl;
echoinput << "link to timevary parms: " << SR_parm_timevary << endl;
}
+
+// SR_update_SSBpR0_rd values
+// 1 best: update SSBpR0 for benchmark and for time series only if SRparm R0 or h (not regime) is set to have time-varying property
+// 2 incorrect, relic: always update SSBpR0 for benchmark's use of spawner-recruitment (old, incorrect SS3 approach), but only for the time series if there is a timevary SR parm
+// 3 option: do not update SSBpR0 (do keep start year SPR0), even if R0 or h is set to have time-varying property
+
+ switch (SR_update_SSBpR0_rd)
+ {
+ case 0:
+ if(timevary_MG_firstyr < YrMax) // timevary biology exists and SR_update not set
+ {
+ warnstream << "user must select 1, 2, or 3 for updating SPR0 flag (formerly labelled future feature in SR input) because there is time-varying biology";
+ write_message (FATAL, 0); // EXIT!
+ }
+ break;
+ case 1:
+ SR_update_SSBpR0_bmark = 1 * SR_update_parm; // but conditional on SRparm_timevary, so update value there
+ SR_update_SSBpR0_timeseries = 1 * SR_update_parm;
+ break;
+ case 2:
+ SR_update_SSBpR0_bmark = 1;
+ SR_update_SSBpR0_timeseries = 0;
+ break;
+ case 3:
+ SR_update_SSBpR0_bmark = 0;
+ SR_update_SSBpR0_timeseries = 0;
+ break;
+ }
echoinput << "SR_Npar and N_SRparm2 and N_SRparm3: " << N_SRparm(SR_fxn) << " " << N_SRparm2 << " " << N_SRparm3 << endl;
// clang-format off
END_CALCS
diff --git a/SS_write_ssnew.tpl b/SS_write_ssnew.tpl
index ea8b9375..4f692468 100644
--- a/SS_write_ssnew.tpl
+++ b/SS_write_ssnew.tpl
@@ -2125,7 +2125,11 @@ FUNCTION void write_nucontrol()
report4 << "#" << endl;
report4 << SR_fxn << " #_Spawner-Recruitment; Options: 1=NA; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=survival_3Parm; 8=Shepherd_3Parm; 9=RickerPower_3parm" << endl;
report4 << init_equ_steepness << " # 0/1 to use steepness in initial equ recruitment calculation" << endl;
- report4 << sigmaR_dendep << " # future feature: 0/1 to make realized sigmaR a function of SR curvature" << endl;
+ report4 << SR_update_SSBpR0_rd << "# SR_update_SSBpR0" << endl <<
+ "# 0 - OK, but only if no timevary biology or SR parm" << endl <<
+ "# 1 - best: update SSBpR0 for benchmark and for time series only if SRparm R0 or h (not regime) is set to have time-varying property" << endl <<
+ "# 2 - incorrect (old, incorrect SS3 approach): always update SSBpR0 for benchmark's use of spawner-recruitment, but only for the time series if there is a timevary SR parm" << endl <<
+ "# 3 - option: do not update SSBpR0 (do keep start year SSBpR0), even if R0 or h is set to have time-varying property" << endl << "#" << endl;
report4 << "#_ LO HI INIT PRIOR PR_SD PR_type PHASE env-var use_dev dev_mnyr dev_mxyr dev_PH Block Blk_Fxn # parm_name" << endl;
report4.unsetf(std::ios_base::fixed);
report4.unsetf(std::ios_base::floatfield);
From 735fff77e211f11642f9aaf09319543e3c88fb93 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Fri, 4 Oct 2024 14:58:13 -0700
Subject: [PATCH 26/29] too many changes
---
SS_benchfore.tpl | 82 +++++++++++++++++++++++-------------------
SS_global.tpl | 4 +--
SS_param.tpl | 11 +++---
SS_popdyn.tpl | 36 +++++++++----------
SS_proced.tpl | 6 ++--
SS_readcontrol_330.tpl | 21 +++++------
SS_readdata_330.tpl | 4 +--
SS_write.tpl | 2 +-
SS_write_report.tpl | 44 +++++++----------------
SS_write_ssnew.tpl | 6 ++--
10 files changed, 104 insertions(+), 112 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index d2cee309..702beb37 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -4,6 +4,15 @@
// SS_Label_file # * <u>get_forecast()</u> // calculates forecast quantities, includes all popdy characteristics of the time series, writes forecast-report.sso
// SS_Label_file #
+// Terminology
+// SSB refers to spawning stock biomass, calculated from reproductive output at age (fec()) and numbers-at-age at spawn_month in spawn_seas
+// SSBpR refers to SSB per recruit calculated with equilibrium age composition in equil_calc
+// SPR refers to spawner potential ratio which is the ratio of SSBpR at some level of F to SSBpR with F = 0
+
+// SSBpR_virgin calculated and reported, but never used
+// SSBpR_virgin_adj used only in the alpha-beta spawner-recruitment. _adj means could be updated if SR_update_SSBpR0_timeseries == 1. Also used to get alpha in equil_spawn_recr B-H
+// but note that also uses SSB_virgin_use. So need to align _use with _adj
+
FUNCTION void setup_Benchmark() // and forecast
{
// SS_Label_Info_7.5 #Get averages from selected years to use in forecasts
@@ -113,12 +122,7 @@ FUNCTION void setup_Benchmark() // and forecast
}
}
t = styr + (endyr + 1 - styr) * nseas + spawn_seas - 1;
- fec = Wt_Age_t(t, -2);
- // for (g=1;g<=gmorph;g++)
- // if(use_morph(g)>0 && sx(g)==1)
- // {
- // fec(g)=save_sel_num(t,0,g);
- // }
+// fec = Wt_Age_t(t, -2); this will always be overwritten, so deleting
if (Fcast_Loop_Control(3) == 3) // using mean recr_dist from range of years
{
@@ -569,7 +573,6 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
eq_yr = y;
t_base = y + (y - styr) * nseas - 1;
bio_t_base = styr + (bio_yr - styr) * nseas - 1;
-
// set the Hrate for bycatch fleets so not scaled with other fleets
// bycatch_F(f,s) is created here for use in forecast
for (f = 1; f <= Nfleet; f++)
@@ -643,6 +646,8 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Make_AgeLength_Key(s, subseas);
// SPAWN-RECR: call make_fecundity for benchmarks
+ // this means that any calculation of SSB in benchmark will use the updated fec
+ // but SSBpR0 will only use that updated fec if SR_update_SSBpR0_bmark == 1
if (s == spawn_seas)
{
{
@@ -744,7 +749,6 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
SR_parm_work(j) = temp / (Bmark_Yr(10) - Bmark_Yr(9) + 1.);
}
}
-
if (SR_update_SSBpR0_bmark == 0) // use virgin biology for the spawner-recruitment R0,h calculations in bmark
{
Recr_unf = Recr_virgin;
@@ -763,14 +767,17 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
if (show_MSY == 1)
{
report5 << "SR_parms for benchmark: " << SR_parm_work << endl
- << "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl;
- if ( SR_update_SSBpR0_bmark == 1) report5 << "SPR0 for equilibrium spawner-recruit based on benchmark biology, not virgin biology" << endl;
+ << "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << " flag for updating SSBpR0 = " << SR_update_SSBpR0_bmark << endl;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
report5 << " Virgin SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SSBpR_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
-
- Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
- report5 << " Benchmark SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " " << SSBpR_unf << " equil: " << Equ_SpawnRecr_Result << endl;
- report5 << "Repro_output_by_age_for_morph_1: " << fec(1) << endl;
+ if ( SR_update_SSBpR0_bmark == 1)
+ {
+ report5 << "SPR0 for equilibrium spawner-recruit based on benchmark biology, not virgin biology" << endl;
+ Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
+ report5 << " Benchmark SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " " << SSBpR_unf << " equil: " << Equ_SpawnRecr_Result << endl;
+ }
+ report5 << 0 << " y: " << y << " Repro_output_by_age_for_morph_1 bench " << fec(1) << endl;
+// report5 << "Repro_output_by_age_for_morph_1: " << fec(1) << endl;
}
SR_parm_work(N_SRparm2 + 1) = SSB_unf;
@@ -1849,7 +1856,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;
@@ -1869,7 +1876,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;
@@ -1889,7 +1896,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;
@@ -1906,9 +1913,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;
@@ -1934,6 +1941,8 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
<< " " << Vbio1_MSY * Recr_msy << " # " << endl;
}
write_bodywt = write_bodywt_save;
+ report5 << 0 << " y: " << y << " Repro_output_by_age_for_morph_1 end_bench: " << fec(1) << endl;
+// report5 << "Repro_output_by_age_for_morph_1_after_benchmark: " << fec(1) << endl;
} // end benchmarks
FUNCTION void Get_Forecast()
@@ -1945,8 +1954,6 @@ FUNCTION void Get_Forecast()
dvariable OFL_catch;
dvariable Fcast_Crash;
dvariable totcatch;
- dvariable R0_use; // annually updated value if SR_update_SSBpR0_timeseries == 1
- dvariable SSB_use; // selected version of SSB that gets passes to Spawn_Recr
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
@@ -2220,6 +2227,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
@@ -2290,7 +2299,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)
{
@@ -2301,12 +2310,12 @@ FUNCTION void Get_Forecast()
}
else
{
- get_mat_fec();
+ get_mat_fec(); // does spawnseas and stores in wt_Age_t(t, -2)
}
}
}
}
-
+ report5 << Fcast_Loop1 << " y: " << y << " updated_Repro_output_by_age_for_morph_1 endyr: " << fec(1) << endl;
for (y = endyr + 1; y <= YrMax; y++)
{
t_base = styr + (y - styr) * nseas - 1;
@@ -2496,6 +2505,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)
@@ -2575,18 +2585,18 @@ FUNCTION void Get_Forecast()
if (Hermaphro_Option != 0) // get male biomass
{
- MaleSPB(y).initialize();
+ MaleSSB(y).initialize();
for (p = 1; p <= pop; p++)
{
for (g = 1; g <= gmorph; g++)
if (sx(g) == 2 && use_morph(g) > 0) // male; all assumed to be mature
{
- MaleSPB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * natage(t, p, g); // accumulates SSB by area and by growthpattern
+ MaleSSB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * natage(t, p, g); // accumulates SSB by area and by growthpattern
}
}
- if (Hermaphro_maleSPB > 0.0) // add MaleSPB to female SSB
+ if (Hermaphro_maleSSB > 0.0) // add MaleSSB to female SSB
{
- SSB_current += Hermaphro_maleSPB * sum(MaleSPB(y));
+ SSB_current += Hermaphro_maleSSB * sum(MaleSSB(y));
SSB_yr(y) = SSB_current;
}
}
@@ -2615,7 +2625,6 @@ FUNCTION void Get_Forecast()
}
SSB_use = SSB_equil;
}
-
Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr
if (SR_fxn != 7) apply_recdev(Recruits, R0_use); // apply recruitment deviation
if (Fcast_Loop1 < Fcast_Loop_Control(2)) // use expected recruitment this should include environ effect - CHECK THIS
@@ -3212,18 +3221,18 @@ FUNCTION void Get_Forecast()
if (Hermaphro_Option != 0) // get male biomass
{
- MaleSPB(y).initialize();
+ MaleSSB(y).initialize();
for (p = 1; p <= pop; p++)
{
for (g = 1; g <= gmorph; g++)
if (sx(g) == 2 && use_morph(g) > 0) // male; all assumed to be mature
{
- MaleSPB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * elem_prod(natage(t, p, g), mfexp(-Z_rate(t, p, g) * spawn_time_seas)); // accumulates SSB by area and by growthpattern
+ MaleSSB(y, p, GP4(g)) += Wt_Age_t(t, 0, g) * elem_prod(natage(t, p, g), mfexp(-Z_rate(t, p, g) * spawn_time_seas)); // accumulates SSB by area and by growthpattern
}
}
- if (Hermaphro_maleSPB > 0.0) // add MaleSPB to female SSB
+ if (Hermaphro_maleSSB > 0.0) // add MaleSSB to female SSB
{
- SSB_current += Hermaphro_maleSPB * sum(MaleSPB(y));
+ SSB_current += Hermaphro_maleSSB * sum(MaleSSB(y));
SSB_yr(y) = SSB_current;
}
}
@@ -3662,7 +3671,7 @@ 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
@@ -3670,7 +3679,7 @@ FUNCTION void Get_Forecast()
Smry_Table(y, 13) = GenTime;
if( SR_fxn == 10 )
{
- temp = SSB_equil / equ_Recr; // current year's SPB/R with current biology at age
+ temp = SSB_equil / equ_Recr; // current year's SSB/R with current biology at age
alpha = mfexp(SR_parm_work(3));
beta = mfexp(SR_parm_work(4));
SR_parm_byyr(y, 2) = alpha * temp / (4. + alpha * temp); // implied steepness
@@ -3678,6 +3687,7 @@ FUNCTION void Get_Forecast()
}
Fishon = 1;
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation
+// warning <<y<<" SSB_virgin: " << SSB_virgin<< "SSB_use: "<<SSB_use<<" SSB_unf: "<<SSB_unf<<" fec(15): "<<fec(1,15)<<endl;
if (STD_Yr_Reverse_Ofish(y) > 0)
SPR_std(STD_Yr_Reverse_Ofish(y)) = SSB_equil / Smry_Table(y, 11);
Smry_Table(y, 9) = totbio;
diff --git a/SS_global.tpl b/SS_global.tpl
index f2e976e1..aed4dd12 100644
--- a/SS_global.tpl
+++ b/SS_global.tpl
@@ -1270,10 +1270,10 @@ REPORT_SECTION
if (Do_TG > 0)
report << " TG-fleetcomp " << TG_like1 << endl
<< " TG-negbin " << TG_like2 << endl;
- report << " -log(L): " << obj_fun << " Spbio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr)) << endl;
+ report << " -log(L): " << obj_fun << " SSBio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr)) << endl;
report << endl
- << "Year Spbio Recruitment" << endl;
+ << "Year SSBio Recruitment" << endl;
report << "Virg " << SSB_yr(styr - 2) << " " << exp_rec(styr - 2, 4) << endl;
report << "Init " << SSB_yr(styr - 1) << " " << exp_rec(styr - 1, 4) << endl;
for (y = styr; y <= endyr; y++)
diff --git a/SS_param.tpl b/SS_param.tpl
index 6a454579..df290682 100644
--- a/SS_param.tpl
+++ b/SS_param.tpl
@@ -159,7 +159,7 @@ PARAMETER_SECTION
3darray recr_dist_endyr(1,N_GP*gender,1,N_settle_timings,1,pop);
!!// SS_Label_Info_5.1.2 #Create SR_parm vector, recruitment vectors
init_bounded_number_vector SR_parm(1,N_SRparm3,SR_parm_LO,SR_parm_HI,SR_parm_PH)
- matrix SR_parm_byyr(styr-3,YrMax,1,N_SRparm2+1) // R0, steepness, parm3, sigmar, rec_dev_offset, R1, rho, SPB Time_vary implementation of spawner-recruitment
+ matrix SR_parm_byyr(styr-3,YrMax,1,N_SRparm2+1) // R0, steepness, parm3, sigmar, rec_dev_offset, R1, rho, SSB Time_vary implementation of spawner-recruitment
vector SR_parm_virg(1,N_SRparm2+1)
vector SR_parm_work(1,N_SRparm2+1)
number two_sigmaRsq;
@@ -241,9 +241,12 @@ PARAMETER_SECTION
number Recr_virgin
number SSB_vir_LH
- number SSB_unf
+ number SSB_unf // SSB unfished, based on benchmark biology
number Recr_unf
+ number SSB_use
+ number R0_use; // annually updated value if SR_update_SSBpR0_timeseries == 1
+ number SSB_deplete // SSB that will be used as denominator for depletion calculations and as basis for control rule inflection
number SSB_current; // Spawning biomass
number SSB_equil;
@@ -322,7 +325,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);
@@ -602,7 +605,7 @@ PARAMETER_SECTION
!!// SS_Label_Info_5.1.8 #Create matrix called smry to store derived quantities of interest
matrix Smry_Table(styr-3,YrMax,1,20+2*gmorph);
- // 1=totbio, 2=smrybio, 3=smrynum, 4=enc_catch, 5=dead_catch, 6=ret_catch, 7=spbio, 8=recruit,
+ // 1=totbio, 2=smrybio, 3=smrynum, 4=enc_catch, 5=dead_catch, 6=ret_catch, 7=SSBio, 8=recruit,
// 9=equ_totbio, 10=equ_smrybio, 11=equ_SSB_virgin, 12=equ_S1, 13=Gentime, 14=YPR, 15=meanage_spawners, 16=meanage_smrynums, 17=meanage_catch
// 18, 19, 20 not used
// 21+cumF-bymorph, maxF-by morph
diff --git a/SS_popdyn.tpl b/SS_popdyn.tpl
index e5fda47c..e97b9bd1 100644
--- a/SS_popdyn.tpl
+++ b/SS_popdyn.tpl
@@ -383,7 +383,7 @@ FUNCTION void get_initial_conditions()
// SSBpR_virgin = SSB_equil / Recr_virgin; // spawners per recruit already calculated
if(Do_Benchmark==0) // assign values that would be created in benchmark section
{
- Mgmt_quant(1) = SSB_virgin;
+ Mgmt_quant(1) = SSB_virgin; // can be overwritten in benchmark by updated SSB_unf
SSB_unf = SSB_virgin;
Recr_unf = Recr_virgin;
Mgmt_quant(2) = totbio; // from equil calcs
@@ -395,7 +395,7 @@ FUNCTION void get_initial_conditions()
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;
@@ -527,7 +527,7 @@ 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;
@@ -554,7 +554,7 @@ 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;
@@ -697,14 +697,12 @@ FUNCTION void get_time_series()
dvariable crashtemp1;
dvariable interim_tot_catch;
dvariable Z_adjuster;
- dvariable R0_use; // annually updated variable if SR_update_SSBpR0_timeseries == 1; gets passed to Spawn_Recr() function
- 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)
{
@@ -983,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();
@@ -1007,18 +1005,18 @@ 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;
}
}
@@ -1468,18 +1466,18 @@ 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) * 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;
}
}
@@ -1795,7 +1793,7 @@ FUNCTION void get_time_series()
Smry_Table(y, 13) = GenTime;
if( SR_fxn == 10 )
{
- temp = SSB_equil / Recr_virgin; // current year's SPB/R with current biology at age
+ temp = SSB_equil / Recr_virgin; // current year's SSB/R with current biology at age
alpha = mfexp(SR_parm_work(3));
beta = mfexp(SR_parm_work(4));
SR_parm_byyr(y, 2) = alpha * temp / (4. + alpha * temp); // implied steepness
@@ -2302,9 +2300,9 @@ FUNCTION void Do_Equil_Calc(const prevariable& equ_Recr)
GenTime /= SSB_equil;
smryage /= smrynum;
cumF /= (r_ages(nages) - r_ages(Smry_Age) + 1.);
- if (Hermaphro_maleSPB > 0.0) // add MaleSPB to female SSB
+ if (Hermaphro_maleSSB > 0.0) // add MaleSSB to female SSB
{
- SSB_equil += Hermaphro_maleSPB * sum(MaleSSB_equil_pop_gp);
+ SSB_equil += Hermaphro_maleSSB * sum(MaleSSB_equil_pop_gp);
}
} // end equil calcs
diff --git a/SS_proced.tpl b/SS_proced.tpl
index b5fc8025..831e17ed 100644
--- a/SS_proced.tpl
+++ b/SS_proced.tpl
@@ -205,15 +205,15 @@ PROCEDURE_SECTION
temp = sqrt((temp + 0.0000001) / (double(recdev_end - recdev_start + 1)));
if (mcmc_counter == 0 && mceval_counter == 0)
{
- cout << current_phase() << " " << niter << " -log(L): " << obj_fun << " Spbio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
+ cout << current_phase() << " " << niter << " -log(L): " << obj_fun << " SSBio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
}
else if (mcmc_counter > 0)
{
- cout << " MCMC: " << mcmc_counter << " -log(L): " << obj_fun << " Spbio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
+ cout << " MCMC: " << mcmc_counter << " -log(L): " << obj_fun << " SSBio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
}
else if (mceval_counter > 0)
{
- cout << " MCeval: " << mceval_counter << " -log(L): " << obj_fun << " Spbio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
+ cout << " MCeval: " << mceval_counter << " -log(L): " << obj_fun << " SSBio: " << value(SSB_yr(styr)) << " " << value(SSB_yr(endyr));
}
if (F_Method > 1 && sum(catch_like) > 0.01)
{
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index 83f2632c..14bcf479 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -1148,11 +1148,11 @@
int Hermaphro_seas;
int Hermaphro_firstage;
number Hermaphro_seas_rd;
- number Hermaphro_maleSPB;
+ number Hermaphro_maleSSB;
LOCAL_CALCS
// clang-format on
Hermaphro_seas = 0;
- Hermaphro_maleSPB = 0.0;
+ Hermaphro_maleSSB = 0.0;
Hermaphro_firstage = 0;
MGparm_Hermaphro = 0;
@@ -1175,8 +1175,8 @@
// so 2.3 will do switch in season 2 beginning with age 3.
echoinput << Hermaphro_seas << " Hermaphro_season (-1 means all seasons)" << endl;
echoinput << Hermaphro_firstage << " Hermaphro_firstage (from decimal part of seas input; note that firstage can only be a single digit, so 9 is max" << endl;
- *(ad_comm::global_datafile) >> Hermaphro_maleSPB; // read as a fraction (0.0 to 1.0) of the male SSB added into the total SSB
- echoinput << Hermaphro_maleSPB << " Hermaphro_maleSPB " << endl;
+ *(ad_comm::global_datafile) >> Hermaphro_maleSSB; // read as a fraction (0.0 to 1.0) of the male SSB added into the total SSB
+ echoinput << Hermaphro_maleSSB << " Hermaphro_maleSSB " << endl;
}
// clang-format off
END_CALCS
@@ -2183,6 +2183,7 @@
SR_update_SSBpR0_timeseries = 0;
break;
}
+ echoinput << "SRupdate " << SR_update_parm << " use " << SR_update_SSBpR0_bmark << " rd " << SR_update_SSBpR0_rd<< endl;
echoinput << "SR_Npar and N_SRparm2 and N_SRparm3: " << N_SRparm(SR_fxn) << " " << N_SRparm2 << " " << N_SRparm3 << endl;
// clang-format off
END_CALCS
@@ -5795,7 +5796,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;
@@ -6849,23 +6850,23 @@
}
}
- // output ln(SPB) std for selected years
- echoinput << " do ln(SPB) std labels for 3 years" << endl;
+ // output ln(SSB) std for selected years
+ echoinput << " do ln(SSB) std labels for 3 years" << endl;
CoVar_Count++;
j++;
active_parm(CoVar_Count) = j;
sprintf(onenum, "%d", styr);
- ParmLabel += "ln(SPB)_" + onenum + CRLF(1);
+ ParmLabel += "ln(SSB)_" + onenum + CRLF(1);
CoVar_Count++;
j++;
active_parm(CoVar_Count) = j;
sprintf(onenum, "%d", int((endyr + styr) / 2));
- ParmLabel += "ln(SPB)_" + onenum + CRLF(1);
+ ParmLabel += "ln(SSB)_" + onenum + CRLF(1);
CoVar_Count++;
j++;
active_parm(CoVar_Count) = j;
sprintf(onenum, "%d", endyr);
- ParmLabel += "ln(SPB)_" + onenum + CRLF(1);
+ ParmLabel += "ln(SSB)_" + onenum + CRLF(1);
if (Do_se_smrybio > 0)
{
diff --git a/SS_readdata_330.tpl b/SS_readdata_330.tpl
index f10d7f98..b300f874 100644
--- a/SS_readdata_330.tpl
+++ b/SS_readdata_330.tpl
@@ -5034,10 +5034,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_write.tpl b/SS_write.tpl
index 96d516b3..a3664db4 100644
--- a/SS_write.tpl
+++ b/SS_write.tpl
@@ -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 << " ";
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index 4a48fa0f..5a11e6bf 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -1483,8 +1483,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
@@ -1496,7 +1496,7 @@ FUNCTION void write_bigoutput()
{
SS2out << SSB_pop_gp(y, p);
if (Hermaphro_Option != 0)
- SS2out << MaleSPB(y, p);
+ SS2out << MaleSSB(y, p);
}
else
{
@@ -1661,7 +1661,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);
@@ -4683,33 +4683,11 @@ FUNCTION void SPR_profile()
for (s = 1; s <= nseas; s++)
{
t = styr - 3 * nseas + s - 1;
-
- if (MG_active(2) > 0 || MG_active(3) > 0 || save_for_report > 0 || WTage_rd > 0)
- {
- subseas = 1;
- ALK_idx = (s - 1) * N_subseas + subseas; // for midseason
- Make_AgeLength_Key(s, subseas); // for begin season
- subseas = mid_subseas;
- ALK_idx = (s - 1) * N_subseas + subseas; // for midseason
- Make_AgeLength_Key(s, subseas); // for midseason
- if (s == spawn_seas)
- {
- subseas = spawn_subseas;
- if (spawn_subseas != 1 && spawn_subseas != mid_subseas)
- {
- //don't call get_growth3(subseas) because using an average ave_size
- Make_AgeLength_Key(s, subseas); // spawn subseas
- }
- get_mat_fec();
- }
- }
-
- for (g = 1; g <= gmorph; g++)
- if (use_morph(g) > 0)
- {
- ALK_idx = (s - 1) * N_subseas + mid_subseas; // for midseason
- Make_FishSelex();
- }
+ Wt_Age_beg(s) = Wt_Age_t(t, 0); // used for smrybio
+ Wt_Age_mid(s) = Wt_Age_t(t, -1);
+ if (s == spawn_seas)
+ fec = Wt_Age_t(t, -2);
+ report5 << 0 << " y: " << y << " updated_Repro_output spr/ypr: " << fec(1) << endl;
}
SS2out << "SPRloop Iter Bycatch Fmult F_report SPR YPR_dead YPR_dead*Recr YPR_ret*Recr Revenue Cost Profit SSB Recruits SSB/Bzero Tot_Catch ";
@@ -4963,7 +4941,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;
@@ -5042,7 +5020,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 4f692468..b98fefa7 100644
--- a/SS_write_ssnew.tpl
+++ b/SS_write_ssnew.tpl
@@ -1537,7 +1537,7 @@ FUNCTION void write_nucontrol()
NuStart << final_conv << " # final convergence criteria (e.g. 1.0e-04) " << endl;
NuStart << retro_yr - endyr << " # retrospective year relative to end year (e.g. -4)" << endl;
NuStart << Smry_Age << " # min age for calc of summary biomass" << endl;
- NuStart << depletion_basis_rd << " # Depletion basis: denom is: 0=skip; 1=X*SPBvirgin; 2=X*SPBmsy; 3=X*SPB_styr; 4=X*SPB_endyr; 5=X*dyn_Bzero; values>=11 invoke N multiyr with 10s & 100s digit; append .1 to invoke log(ratio); e.g. 122.1 produces log(12 yr trailing average of B/Bmsy)" << endl;
+ NuStart << depletion_basis_rd << " # Depletion basis: denom is: 0=skip; 1=X*SSBvirgin; 2=X*SSBmsy; 3=X*SSB_styr; 4=X*SSB_endyr; 5=X*dyn_Bzero; values>=11 invoke N multiyr with 10s & 100s digit; append .1 to invoke log(ratio); e.g. 122.1 produces log(12 yr trailing average of B/Bmsy)" << endl;
NuStart << depletion_level << " # Fraction (X) for Depletion denominator (e.g. 0.4)" << endl;
NuStart << SPR_reporting << " # SPR_report_basis: 0=skip; 1=(1-SPR)/(1-SPR_tgt); 2=(1-SPR)/(1-SPR_MSY); 3=(1-SPR)/(1-SPR_Btarget); 4=rawSPR" << endl;
NuStart << F_reporting << " # F_std_reporting_units: 0=skip; 1=exploitation(Bio); 2=exploitation(Num); 3=sum(Apical_F's); 4=mean F for range of ages (numbers weighted); 5=unweighted mean F for range of ages" << endl;
@@ -1930,7 +1930,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;
@@ -2125,7 +2125,7 @@ FUNCTION void write_nucontrol()
report4 << "#" << endl;
report4 << SR_fxn << " #_Spawner-Recruitment; Options: 1=NA; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=survival_3Parm; 8=Shepherd_3Parm; 9=RickerPower_3parm" << endl;
report4 << init_equ_steepness << " # 0/1 to use steepness in initial equ recruitment calculation" << endl;
- report4 << SR_update_SSBpR0_rd << "# SR_update_SSBpR0" << endl <<
+ report4 << SR_update_SSBpR0_rd << " # SR_update_SSBpR0" << endl <<
"# 0 - OK, but only if no timevary biology or SR parm" << endl <<
"# 1 - best: update SSBpR0 for benchmark and for time series only if SRparm R0 or h (not regime) is set to have time-varying property" << endl <<
"# 2 - incorrect (old, incorrect SS3 approach): always update SSBpR0 for benchmark's use of spawner-recruitment, but only for the time series if there is a timevary SR parm" << endl <<
From e33a3bdae2a8fa5169a088da661793af8edefcec Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Fri, 4 Oct 2024 16:23:39 -0700
Subject: [PATCH 27/29] add HCR_anchor in forecast.ss
---
SS_benchfore.tpl | 37 ++++++++++++++++++++-----------------
SS_param.tpl | 1 +
SS_readdata_330.tpl | 4 ++--
SS_write_report.tpl | 6 +++---
SS_write_ssnew.tpl | 3 ++-
5 files changed, 28 insertions(+), 23 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index 702beb37..c2293f20 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -776,8 +776,6 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
report5 << " Benchmark SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " " << SSBpR_unf << " equil: " << Equ_SpawnRecr_Result << endl;
}
- report5 << 0 << " y: " << y << " Repro_output_by_age_for_morph_1 bench " << fec(1) << endl;
-// report5 << "Repro_output_by_age_for_morph_1: " << fec(1) << endl;
}
SR_parm_work(N_SRparm2 + 1) = SSB_unf;
@@ -2211,6 +2209,11 @@ FUNCTION void Get_Forecast()
break;
}
}
+ if (Fcast_Loop_Control(5) <= 1)
+ {HCR_anchor = SSB_unf;}
+ else if (Fcast_Loop_Control(5) ==2)
+ {HCR_anchor = SSB_virgin;}
+ report5 << "Control_rule_anchor_approach: " << Fcast_Loop_Control(5) << " HCR_anchor: " << HCR_anchor << endl;
report5 << "#" << endl;
}
@@ -2227,7 +2230,7 @@ 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;
+// 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
{
@@ -2315,7 +2318,7 @@ FUNCTION void Get_Forecast()
}
}
}
- report5 << Fcast_Loop1 << " y: " << y << " updated_Repro_output_by_age_for_morph_1 endyr: " << fec(1) << endl;
+// 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;
@@ -2505,7 +2508,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;
+// 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)
@@ -2659,9 +2662,9 @@ FUNCTION void Get_Forecast()
}
else if (ABC_Loop == 2 && s == 1) // Calc the buffer in season 1, will use last year's spawnbio if multiseas and spawnseas !=1
{
- temp = SSB_unf;
- join1 = 1. / (1. + mfexp(10. * (SSB_current - H4010_bot * temp)));
- join2 = 1. / (1. + mfexp(10. * (SSB_current - H4010_top * temp)));
+
+ join1 = 1. / (1. + mfexp(10. * (SSB_current - H4010_bot * HCR_anchor)));
+ join2 = 1. / (1. + mfexp(10. * (SSB_current - H4010_top * HCR_anchor)));
switch (HarvestPolicy)
{
@@ -2674,8 +2677,8 @@ FUNCTION void Get_Forecast()
// ramp scales catch as f(B) and buffer (H4010_scale) applied to F
{
ABC_buffer(y) = H4010_scale_vec(y) *
- ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
- + (0.0001 + (1.0 - 0.0001) * (H4010_top * temp / SSB_current) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+ ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+ + (0.0001 + (1.0 - 0.0001) * (H4010_top * HCR_anchor / SSB_current) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
) *
(join2) // scale combo
+
@@ -2686,8 +2689,8 @@ FUNCTION void Get_Forecast()
// ramp scales F as f(B) and buffer (H4010_scale) applied to F
{
ABC_buffer(y) = H4010_scale_vec(y) *
- ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
- + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+ ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+ + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
) *
(join2) // scale combo
+
@@ -2698,8 +2701,8 @@ FUNCTION void Get_Forecast()
// ramp scales catch as f(B) and buffer (H4010_scale) applied to catch
{
ABC_buffer(y) = 1.0 *
- ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
- + (0.0001 + (1.0 - 0.0001) * (H4010_top * temp / SSB_current) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+ ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+ + (0.0001 + (1.0 - 0.0001) * (H4010_top * HCR_anchor / SSB_current) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
) *
(join2) // scale combo
+
@@ -2710,8 +2713,8 @@ FUNCTION void Get_Forecast()
// ramp scales F as f(B) and buffer (H4010_scale) applied to catch
{
ABC_buffer(y) = 1.0 *
- ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
- + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+ ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+ + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
) *
(join2) // scale combo
+
@@ -3514,7 +3517,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);
}
diff --git a/SS_param.tpl b/SS_param.tpl
index df290682..3c6147ea 100644
--- a/SS_param.tpl
+++ b/SS_param.tpl
@@ -235,6 +235,7 @@ 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) bor inflection in control rule
// SPAWN-RECR: define some spawning biomass and recruitment entities
number SSB_virgin
diff --git a/SS_readdata_330.tpl b/SS_readdata_330.tpl
index b300f874..4b0a3ea3 100644
--- a/SS_readdata_330.tpl
+++ b/SS_readdata_330.tpl
@@ -4344,8 +4344,8 @@
"even when the base is set to the mean of earlier recruitments" << endl;
}
- echoinput << Fcast_Loop_Control(5) << " #echo: loop control 5 not used" << endl;
-
+ echoinput << Fcast_Loop_Control(5) << " #control rule anchor: 1=unfished_benchmark_SSB(old_approach), 2=virgin_SSB " << endl;
+
echoinput << "#next enter year in which Fcast loop 3 caps and allocations begin to be applied" << endl;
*(ad_comm::global_datafile) >> Fcast_Cap_FirstYear;
echoinput << Fcast_Cap_FirstYear << " # echoed value" << endl;
diff --git a/SS_write_report.tpl b/SS_write_report.tpl
index 5a11e6bf..cf0c4a67 100644
--- a/SS_write_report.tpl
+++ b/SS_write_report.tpl
@@ -4687,7 +4687,7 @@ FUNCTION void SPR_profile()
Wt_Age_mid(s) = Wt_Age_t(t, -1);
if (s == spawn_seas)
fec = Wt_Age_t(t, -2);
- report5 << 0 << " y: " << y << " updated_Repro_output spr/ypr: " << fec(1) << endl;
+// report5 << 0 << " y: " << y << " updated_Repro_output spr/ypr: " << fec(1) << endl;
}
SS2out << "SPRloop Iter Bycatch Fmult F_report SPR YPR_dead YPR_dead*Recr YPR_ret*Recr Revenue Cost Profit SSB Recruits SSB/Bzero Tot_Catch ";
@@ -5020,9 +5020,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;
+// 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;
+// 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 b98fefa7..d8a61230 100644
--- a/SS_write_ssnew.tpl
+++ b/SS_write_ssnew.tpl
@@ -1631,6 +1631,7 @@ 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 << "# Also see HCR_anchor below" << 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;
if (H4010_scale_rd < 0)
@@ -1654,7 +1655,7 @@ FUNCTION void write_nucontrol()
{
NuFore << Fcast_Loop_Control(4) << " # multiplier on base recruitment " << endl;
}
- NuFore << Fcast_Loop_Control(5) << " # not used" << endl << "#" << endl;
+ NuFore << Fcast_Loop_Control(5) << " # HCR_anchor: 0 or 1 uses unfished benchmark SSB (old hardwired approach), 2 = virgin SSB" << endl << "#" << endl;
NuFore << Fcast_Cap_FirstYear << " # FirstYear for caps and allocations (should be after years with fixed inputs) " << endl;
From 46fc4ee4357ef8a42f92783f5fa11944ce8ef630 Mon Sep 17 00:00:00 2001
From: Richard Methot <richard.methot@noaa.gov>
Date: Thu, 17 Oct 2024 16:07:03 -0700
Subject: [PATCH 28/29] latest update
---
SS_benchfore.tpl | 25 +++++++++++++------------
SS_param.tpl | 2 +-
SS_write_ssnew.tpl | 2 +-
3 files changed, 15 insertions(+), 14 deletions(-)
diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl
index c2293f20..08c9c645 100644
--- a/SS_benchfore.tpl
+++ b/SS_benchfore.tpl
@@ -767,7 +767,8 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
if (show_MSY == 1)
{
report5 << "SR_parms for benchmark: " << SR_parm_work << endl
- << "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << " flag for updating SSBpR0 = " << SR_update_SSBpR0_bmark << endl;
+ << "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl <<
+ "input.SR_update_SSBpR0_rd: " << SR_update_SSBpR0_rd << "flag for updating SSBpR0_Bmark: " << SR_update_SSBpR0_bmark << endl;
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
report5 << " Virgin SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SSBpR_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
if ( SR_update_SSBpR0_bmark == 1)
@@ -2176,10 +2177,16 @@ FUNCTION void Get_Forecast()
{
H4010_top = H4010_top_rd;
}
+ if (Fcast_Loop_Control(5) <= 1)
+ {HCR_anchor = SSB_unf;}
+ else if (Fcast_Loop_Control(5) ==2)
+ {HCR_anchor = SSB_virgin;}
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;
+ report5 << "#" << endl;
+
switch (HarvestPolicy)
{
case 0: // none
@@ -2209,13 +2216,7 @@ FUNCTION void Get_Forecast()
break;
}
}
- if (Fcast_Loop_Control(5) <= 1)
- {HCR_anchor = SSB_unf;}
- else if (Fcast_Loop_Control(5) ==2)
- {HCR_anchor = SSB_virgin;}
- report5 << "Control_rule_anchor_approach: " << Fcast_Loop_Control(5) << " HCR_anchor: " << HCR_anchor << endl;
- report5 << "#" << endl;
- }
+ }
int jloop;
if (fishery_on_off == 1 || Do_Dyn_Bzero > 0)
@@ -2264,7 +2265,7 @@ FUNCTION void Get_Forecast()
if (HarvestPolicy == 0)
report5 << "pop year ABC_Loop season No_buffer bio-all bio-Smry SpawnBio Depletion recruit-0 ";
if (HarvestPolicy <= 2)
- report5 << "pop year ABC_Loop season Ramp&Buffer bio-all bio-Smry SpawnBio Depletion recruit-0 ";
+ report5 << "pop year ABC_Loop season Ramp&Buffer Buffer2 bio-all bio-Smry SpawnBio Depletion recruit-0 ";
if (HarvestPolicy >= 3)
report5 << "pop year ABC_Loop season Ramp bio-all bio-Smry SpawnBio Depletion recruit-0 ";
for (int ff = 1; ff <= N_catchfleets(0); ff++)
@@ -3134,7 +3135,7 @@ FUNCTION void Get_Forecast()
}
if (show_MSY == 1)
{
- report5 << p << " " << y << " " << ABC_Loop << " " << s << " " << ABC_buffer(y) << " " << totbio << " " << smrybio << " ";
+ report5 << p << " " << y << " " << ABC_Loop << " " << s << " " << ABC_buffer(y) << " " << H4010_scale_vec(y) << " " << totbio << " " << smrybio << " ";
if (s == spawn_seas)
{
report5 << SSB_current << " ";
diff --git a/SS_param.tpl b/SS_param.tpl
index 3c6147ea..9a86f578 100644
--- a/SS_param.tpl
+++ b/SS_param.tpl
@@ -235,7 +235,7 @@ 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) bor inflection in control rule
+ 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
diff --git a/SS_write_ssnew.tpl b/SS_write_ssnew.tpl
index d8a61230..e7080a9a 100644
--- a/SS_write_ssnew.tpl
+++ b/SS_write_ssnew.tpl
@@ -1631,9 +1631,9 @@ 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 << "# Also see HCR_anchor below" << endl;
NuFore << H4010_bot << " # Control rule cutoff for no F (as frac of Bzero, e.g. 0.10) " << endl;
NuFore << H4010_scale_rd << " # Buffer: enter Control rule target as fraction of Flimit (e.g. 0.75), negative value invokes list of [year, scalar] with filling from year to YrMax " << endl;
+ NuFore << "# Also see HCR_anchor below" << endl;
if (H4010_scale_rd < 0)
{
j = H4010_scale_vec_rd.size() - 1;
From d2acb1e435fbe85d7b9c2561e95267adc386c4d1 Mon Sep 17 00:00:00 2001
From: Ian Taylor <ian.taylor@noaa.gov>
Date: Fri, 18 Oct 2024 09:43:01 -0700
Subject: [PATCH 29/29] revert change in capitalization of Recr_Virgin
---
SS_readcontrol_330.tpl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index 14bcf479..0e4dfab2 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -6495,7 +6495,7 @@
active_parm(CoVar_Count) = j;
if (y == styr - 2)
{
- ParmLabel += "Recr_virgin";
+ ParmLabel += "Recr_Virgin";
}
else if (y == styr - 1)
{