Skip to content

Commit

Permalink
new fcast year approach
Browse files Browse the repository at this point in the history
  • Loading branch information
@Rick-Methot-NOAA
Rick-Methot-NOAA committed Aug 15, 2023
1 parent 1e126bd commit ad5dc10
Show file tree
Hide file tree
Showing 5 changed files with 200 additions and 156 deletions.
9 changes: 7 additions & 2 deletions SS_benchfore.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,9 +147,9 @@ FUNCTION void setup_Benchmark()
if (Fcast_Loop_Control(5) == 1) //
{

for (int parm_type = 1; parm_type <= 8; parm_type++)
for (int parm_type = 1; parm_type <= 12; parm_type++)
{
if(Fcast_MGparm_ave(parm_type, 2) == 1) // do averaging of derived biology
if(Fcast_MGparm_ave(parm_type, 2) == 1) // do averaging of derived factor
{
double ave_styr = Fcast_MGparm_ave(parm_type,3);
double ave_endyr = Fcast_MGparm_ave(parm_type,4);
Expand Down Expand Up @@ -248,6 +248,11 @@ FUNCTION void setup_Benchmark()
warnstream << "Maturity & fecundity params averaging is not implemented, execution continues. " ;
write_message (WARN, 1);
break;

case 10: // 9=selectivity
tempvec_a.initialize();
break;

}
}
}
Expand Down
20 changes: 17 additions & 3 deletions SS_readcontrol_330.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1805,6 +1805,20 @@

echoinput << y << " timevary_MG: " << timevary_MG(y) << endl;
}

for (y = endyr + 1; y <= YrMax; y++)
{
for (f = 1; f <= 7; f++)
{
if (timevary_MG(y, f) > 0 && Fcast_MGparm_ave(f,2) > 0)
{
warnstream << "mean MGparm for forecast is incompatible with timevary parm in forecast yr: " << y << "; for type: " << f << "; SS3 will disable time-vary";
write_message(WARN, 0);
timevary_MG(y, f) = 0;
}
}
}

// clang-format off
END_CALCS

Expand Down Expand Up @@ -2251,12 +2265,12 @@
echoinput << Fcast_recr_lambda << " #_lambda for Fcast_recr_like occurring before endyr+1" << endl;
if (Fcast_Loop_Control(3) >= 3 && Fcast_recr_PH_rd >= 0)
{
warnstream << "Mean recruitment for forecast is incompatible with pos. phase for forecast rec_devs; set phase to neg. unless using late rec_devs";
write_message (WARN, 0);
warnstream << "Forecast devs will be applied to mean base recruitment over range of historical years in forecast.ss";
write_message (NOTE, 0);
}
if (Do_Impl_Error > 0 && Fcast_recr_PH_rd < 0)
{
warnstream << "Implementation error incompatible with neg. phase for forecast rec_devs; SS3 will run without active impl error";
warnstream << "Implementation error has null effect unless Fcast_recr_PH is >=0";
write_message (WARN, 0);
}
echoinput << recdev_adj(1) << " #_last_early_yr_nobias_adj_in_MPD" << endl;
Expand Down
238 changes: 137 additions & 101 deletions SS_readdata_330.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4032,13 +4032,14 @@
Do_Rebuilder = 0;
// clang-format off
END_CALCS
matrix Fcast_MGparm_ave_rd(1,8,1,4) // for the 8 MG_types, method, st_year, end_year
matrix Fcast_MGparm_ave(1,8,1,4) // for the 8 MG_types, method, st_year, end_year (real years)
matrix Fcast_MGparm_ave_rd(1,12,1,4) // for the 8 MG_types plus, method, st_year, end_year
matrix Fcast_MGparm_ave(1,12,1,4) // for the 8 MG_types plus, method, st_year, end_year (real years)

LOCAL_CALCS
// clang-format on
Fcast_Specify_Selex = 0; // default

Fcast_MGparm_ave.initialize();
Fcast_MGparm_ave_rd.initialize();
if (Do_Forecast_rd > 0)
{
// Fcast_Input(1,k)=Fcast_Input_rd(1,k);
Expand Down Expand Up @@ -4067,17 +4068,15 @@
*(ad_comm::global_datafile) >> Fcast_Flevel;
echoinput << Fcast_Flevel << " # echoed Fmult value" << endl;

N_Fcast_parm_aves = 0;
echoinput << endl
<< "# next enter Fcast_years: beg_selex, end_selex, beg_relF, end_relF, beg_recruits, end_recruits" << endl
<< "# enter actual year, or values of 0 or -integer to be relative to endyr)" << endl;
*(ad_comm::global_datafile) >> Fcast_yr_rd(1, 6);
// k++; Fcast_yr(1) = int(Fcast_Input(k));
// k++; Fcast_yr(2) = int(Fcast_Input(k));
// k++; Fcast_yr(3) = int(Fcast_Input(k));
// k++; Fcast_yr(4) = int(Fcast_Input(k));
// k++; Fcast_yr(5) = int(Fcast_Input(k));
// k++; Fcast_yr(6) = int(Fcast_Input(k));

*(ad_comm::global_datafile) >> Fcast_yr_rd(1);

if(Fcast_yr_rd(1) != -12345) // continue with old approach
{
*(ad_comm::global_datafile) >> Fcast_yr_rd(2,6);
echoinput << Fcast_yr_rd << " # echoed Fcast years as read" << endl;
Fcast_yr = Fcast_yr_rd;
for (i = 1; i <= 6; i++)
Expand Down Expand Up @@ -4118,100 +4117,49 @@
warnstream << "Fcast_specify_selex is 0 but user should change to 1 if density dependence affects a selectivity parameter or growth "<<endl;
write_message(WARN, 0);
}
// set equivalent values using new approach
Fcast_MGparm_ave(10,1) = 10;
Fcast_MGparm_ave(10,2) = 1;
Fcast_MGparm_ave(10,3) = Fcast_yr_rd(1);
Fcast_MGparm_ave(10,4) = Fcast_yr_rd(2);
Fcast_MGparm_ave(11,1) = 11;
Fcast_MGparm_ave(11,2) = 1;
Fcast_MGparm_ave(11,3) = Fcast_yr_rd(3);
Fcast_MGparm_ave(11,4) = Fcast_yr_rd(4);
Fcast_MGparm_ave(12,1) = 12;
Fcast_MGparm_ave(12,2) = 1;
Fcast_MGparm_ave(12,3) = Fcast_yr_rd(5);
Fcast_MGparm_ave(12,4) = Fcast_yr_rd(6);

} // end old approach for Fcast years
else
{ // read fcast year ranges in new list-based format

echoinput << endl
<< "next read 4 values for: control rule shape(0, 1, 2, 3 or 4), inflection (like 0.40), cutoff(like 0.10), scale(like 0.75)" << endl;
*(ad_comm::global_datafile) >> HarvestPolicy;
if (HarvestPolicy == 0)
echoinput << "HarvestPolicy=0, so values for top, bottom, buffer will be ignored" << endl;

echoinput << HarvestPolicy << " # echoed HarvestPolicy " << endl;
*(ad_comm::global_datafile) >> H4010_top_rd; // use -1 to set H4010_top to Bmsy/SSB_unf
echoinput << H4010_top_rd << " # echoed control rule inflection (-1 to set to Bmsy/SSB_unf)" << endl;
*(ad_comm::global_datafile) >> H4010_bot;
echoinput << H4010_bot << " # echoed control rule cutoff " << endl;
*(ad_comm::global_datafile) >> H4010_scale_rd;
H4010_scale = H4010_scale_rd;
echoinput << H4010_scale << " # echoed control rule scalar " << endl;
if (H4010_top_rd > 0.0 && H4010_top_rd <= H4010_bot)
{
warnstream << "control rule inflection: " << H4010_top_rd << " must be > control rule cutoff " << H4010_bot;
write_message(FATAL, 0);
}
if (H4010_scale > 1.0)
{
warnstream << "Sure you want control rule scalar > 1.0? " << H4010_scale;
write_message(WARN, 0);
}

if (H4010_scale < 0.0)
{
echoinput << "# now read pairs of year,H4010scale; each read fills from that year to YrMax; end with year<0.0 " << endl;
ender = 0;
do
{
dvector tempvec(1, 2);
*(ad_comm::global_datafile) >> tempvec(1, 2);
if (tempvec(1) < 0.0)
ender = 1;
H4010_scale_vec_rd.push_back(tempvec(1, 2));
echoinput << " H4010 read: " << tempvec(1, 2) << endl;
} while (ender == 0);
}

echoinput << endl
<< "# next enter 2 values that control looping through the forecast (see manual), then 3 additional controls" << endl;
echoinput << "# first does F_msy or proxy; 2nd applies control rule; 3rd applies caps and allocations" << endl;
*(ad_comm::global_datafile) >> Fcast_Loop_Control(1, 5);
echoinput << Fcast_Loop_Control(1) << " #echoed N forecast loops (1-3) (recommend 3 to get full variance for short-term forecasts)" << endl;
echoinput << Fcast_Loop_Control(2) << " #echoed First forecast loop with stochastic recruitment (recommend 3)" << endl;
echoinput << Fcast_Loop_Control(3) << " #echoed Forecast recruitment: 0=spawn_recr; 1=mult*spawn_recr; 2=mult*VirginRecr; 3=mean from yr range; 4=mult*mean from yr range" << endl;
if (Fcast_Loop_Control(3) == 0)
{
echoinput << Fcast_Loop_Control(4) << " Forecast control #4 (not used) " << endl;
}
else if (Fcast_Loop_Control(3) == 1)
{
echoinput << Fcast_Loop_Control(4) << "multiplier on spawn_recr" << endl;
}
else if (Fcast_Loop_Control(3) == 2)
{
echoinput << Fcast_Loop_Control(4) << "multiplier on virgin recr" << endl;
}
else if (Fcast_Loop_Control(3) == 3)
{
echoinput << "mean recruitment and recr_dist from years: " << Fcast_Rec_yr1 << " to " << Fcast_Rec_yr2 << endl;
}
else if (Fcast_Loop_Control(3) == 4)
{
echoinput << "mult * mean recruitment from years: " << Fcast_Rec_yr1 << " to " << Fcast_Rec_yr2 << " recrdist from parameters, or average using control_5" << endl;
echoinput << Fcast_Loop_Control(4) << " multiplier on recent mean recr" << endl;
}
else // input probably was a -1 from pre 3.30.15, so convert to 0
{
Fcast_Loop_Control(3) = 0;
Fcast_Loop_Control(4) = 1.0;
echoinput << Fcast_Loop_Control(4) << " #echoed Forecast control #4: multiplier on recruitment from spawn_recr" << endl;
}
// set defaults, but each can be overridden
Fcast_yr(1,6) = endyr;
Fcast_Sel_yr1 = Fcast_yr(1);
Fcast_Sel_yr2 = Fcast_yr(2);
Fcast_RelF_yr1 = Fcast_yr(3);
Fcast_RelF_yr2 = Fcast_yr(4);
Fcast_Rec_yr1 = Fcast_yr(5);
Fcast_Rec_yr2 = Fcast_yr(6);
Fcast_Specify_Selex = 1;

echoinput << Fcast_Loop_Control(5) << " #echoed Forecast biology averaging: enter 1 to start list input" << endl;
echoinput << " #_Forecast factors averaging" << endl;
// Fcast_MGparm_ave_rd: read MG_type, method, start year, end year
// terminate with MG_type = -9999
// MG_type: 1=M, 2=growth 3=wtlen, 4=recr_dist&femfrac, 5=migration, 6=ageerror, 7=catchmult 8=hermaphroditism
N_Fcast_parm_aves = 0;
if (Fcast_Loop_Control(5) == 1)
{
echoinput << "Fcast Control(5)==1, so now read list of MG_type, method (1, 2), start year, end year" << endl
echoinput << "read list of MG_type, method (0,1), start year, end year" << endl
<< "Terminate with -9999 for MG_type" << endl
<< "MG_type: 1=M, 2=growth, 3=wtlen, 4=recr_dist&femfrac, 5=migration, 6=ageerror, 7=catchmult, 8=hermaphroditism" << endl
<< "Method = 0 (default) means continue time_vary parms; 1 means to use average of derived biology; 2 (future) means average parameter then apply as if time-vary"<<endl;
<< "Method = 0, or omitted, means continue time_vary parms; 1 means to use average of derived biology; 2 (future) means average parameter then apply as if time-vary"<<endl;
ender = 0;
do
{
dvector tempvec(1, 4);
*(ad_comm::global_datafile) >> tempvec(1, 4);
echoinput << tempvec << endl;
if (tempvec(1) == -9999. || tempvec(1) > 8)
if (tempvec(1) == -9999. || tempvec(1) > 12)
ender = 1;
else
{
Expand All @@ -4223,7 +4171,7 @@

// Adjusting Fcast_MGparm_ave_rd minyear and maxyear values
// for Fcast_MGparm_ave
for (i = 1; i <= 8; i++)
for (i = 1; i <= 12; i++)
{
if (Fcast_MGparm_ave_rd(i,1) > 0)
{
Expand Down Expand Up @@ -4267,16 +4215,97 @@
warnstream << "Fcast_MGparm_ave maxyear exceeds endyr, setting to: " << endyr;
write_message(ADJUST, 0);
}
switch (i)
{
case 10: // 10=selectivity
Fcast_Sel_yr1 = Fcast_MGparm_ave(i,3);
Fcast_Sel_yr2 = Fcast_MGparm_ave(i,4);
Fcast_Specify_Selex = 0; // tells SS3 to use averages, not time-vary parms
break;

case 11: // 11=relative F
Fcast_RelF_yr1 = Fcast_MGparm_ave(i,3);
Fcast_RelF_yr2 = Fcast_MGparm_ave(i,4);
// only year range read here; invocation will be read later: Fcast_RelF_Basis; // tells SS3 to use averages, not time-vary parms
break;

}
}
}
echoinput << "Forecast MGparm averaging: " << endl << Fcast_MGparm_ave << endl;
echoinput << "Forecast factor averaging: " << endl << Fcast_MGparm_ave << endl;
echoinput << "operational values will be calculated or assigned in benchmark_forecast setup" << endl;
}
// clang-format off
END_CALCS

LOCAL_CALCS
// clang-format on
echoinput << "#next enter year in which Fcast loop 3 caps and allocations begin to be applied" << endl;

echoinput << endl
<< "next read 4 values for: control rule shape(0, 1, 2, 3 or 4), inflection (like 0.40), cutoff(like 0.10), scale(like 0.75)" << endl;
*(ad_comm::global_datafile) >> HarvestPolicy;
if (HarvestPolicy == 0)
echoinput << "HarvestPolicy=0, so values for top, bottom, buffer will be ignored" << endl;

echoinput << HarvestPolicy << " # echoed HarvestPolicy " << endl;
*(ad_comm::global_datafile) >> H4010_top_rd; // use -1 to set H4010_top to Bmsy/SSB_unf
echoinput << H4010_top_rd << " # echoed control rule inflection (-1 to set to Bmsy/SSB_unf)" << endl;
*(ad_comm::global_datafile) >> H4010_bot;
echoinput << H4010_bot << " # echoed control rule cutoff " << endl;
*(ad_comm::global_datafile) >> H4010_scale_rd;
H4010_scale = H4010_scale_rd;
echoinput << H4010_scale << " # echoed control rule scalar " << endl;
if (H4010_top_rd > 0.0 && H4010_top_rd <= H4010_bot)
{
warnstream << "control rule inflection: " << H4010_top_rd << " must be > control rule cutoff " << H4010_bot;
write_message(FATAL, 0);
}
if (H4010_scale > 1.0)
{
warnstream << "Sure you want control rule scalar > 1.0? " << H4010_scale;
write_message(WARN, 0);
}

if (H4010_scale < 0.0)
{
echoinput << "# now read pairs of year,H4010scale; each read fills from that year to YrMax; end with year<0.0 " << endl;
ender = 0;
do
{
dvector tempvec(1, 2);
*(ad_comm::global_datafile) >> tempvec(1, 2);
if (tempvec(1) < 0.0)
ender = 1;
H4010_scale_vec_rd.push_back(tempvec(1, 2));
echoinput << " H4010 read: " << tempvec(1, 2) << endl;
} while (ender == 0);
}

echoinput << endl
<< "# next enter 2 values that control looping through the forecast (see manual), then 3 additional controls" << endl;
echoinput << "# first does F_msy or proxy; 2nd applies control rule; 3rd applies caps and allocations" << endl;
*(ad_comm::global_datafile) >> Fcast_Loop_Control(1, 5);
echoinput << Fcast_Loop_Control(1) << " #echo: N forecast loops (1-3) (recommend 3 to get full variance for short-term forecasts)" << endl;
echoinput << Fcast_Loop_Control(2) << " #echo: First forecast loop with stochastic recruitment (recommend 3)" << endl;
echoinput << Fcast_Loop_Control(3) << " #echo: Forecast base recruitment: 0=spawn_recr; 1=mult*spawn_recr; 2=mult*VirginRecr; 3=mean from yr range; 4=mult*mean from yr range" << endl;
if (Fcast_Loop_Control(3) == 3)
{
echoinput << "Forecast base recruitment AND recr_dist is mean from years: " << Fcast_Rec_yr1 << " to " << Fcast_Rec_yr2 << endl;
}
else if (Fcast_Loop_Control(3) == 4)
{
echoinput << "Forecast base recruitment is mean from years: " << Fcast_Rec_yr1 << " to " << Fcast_Rec_yr2 << " recrdist from parameters, or average using control_5" << endl;
}
else if (Fcast_Loop_Control(3) < 0) // input probably was a -1 from pre 3.30.15, so convert to 0
{
Fcast_Loop_Control(3) = 0;
Fcast_Loop_Control(4) = 1.0;
}
if (Fcast_Loop_Control(3) > 0)
{
echoinput << Fcast_Loop_Control(4) << "#echo: multiplier on forecast base recruitment" << endl;
echoinput << "forecast devs will be applied after the multiplier," << endl <<
"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 << "#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;

Expand Down Expand Up @@ -4316,11 +4345,18 @@
<< "#next select fleet relative F: 1=use first-last alloc year read above; 2=read list of seas, fleet, relF below" << endl;
echoinput << "# Note that fleet allocation is used directly as average F if Do_Forecast=4 " << endl;
*(ad_comm::global_datafile) >> Fcast_RelF_Basis;
echoinput << Fcast_RelF_Basis << " # echoed value" << endl;
echoinput << Fcast_RelF_Basis << " # echoed value" << Fcast_MGparm_ave(11,2) << endl;
if (Fcast_RelF_Basis < 1 || Fcast_RelF_Basis > 2) {
warnstream << "Fcast_relF_Basis value must be 1 or 2" << endl;
write_message(FATAL, 1);
}
if (Fcast_RelF_Basis == 1 && Fcast_MGparm_ave(11,2) == 0) {
echoinput << "Fcast_relF_Basis requires that year range is set above" << endl;
warnstream << "Fcast_relF_Basis requires that year range is set above" << endl;
write_message(FATAL, 1);

}

if (Do_Forecast_rd == 4 && Fcast_RelF_Basis == 2) {
warnstream << "Cannot specify forecast fleet relative F because Do_Forecast==4 specifies relative F directly as F;" << endl
<< " need to align choice of forecast basis and forecast relative F basis";
Expand Down
2 changes: 2 additions & 0 deletions SS_recruit.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -258,6 +258,8 @@ FUNCTION void apply_recdev(prevariable& NewRecruits, const prevariable& Recr_vir
break;
}
}
// note that if user requests "mean" as base forecast recr, then devs are still applied
// so, phase for forecast recdevs must be <0 to assure that forecast recr do not get added variability
if (do_recdev > 0)
NewRecruits *= mfexp(Fcast_recruitments(y)); // recruitment deviation
}
Expand Down
Loading

0 comments on commit ad5dc10

Please sign in to comment.