From e33a3bdae2a8fa5169a088da661793af8edefcec Mon Sep 17 00:00:00 2001 From: Richard Methot Date: Fri, 4 Oct 2024 16:23:39 -0700 Subject: [PATCH] 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;