Skip to content

Commit

Permalink
Merge pull request #391 from nmfs-stock-synthesis/improve-time-varyin…
Browse files Browse the repository at this point in the history 
…g-parameter-implementation

fix #165 for bounded link; add survey type #41
  • Loading branch information
@Rick-Methot-NOAA
Rick-Methot-NOAA committed Feb 9, 2023
2 parents f7b1bda + e934175 commit 6b99f1f
Show file tree
Hide file tree
Showing 8 changed files with 82 additions and 41 deletions.
1 change: 0 additions & 1 deletion SS_benchfore.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2204,7 +2204,6 @@ FUNCTION void Get_Forecast()
}
env_data(y, -3) = log(smrybio / Smry_Table(styr - 1, 2));
env_data(y, -4) = log(smrynum / Smry_Table(styr - 1, 3));

Smry_Table(y).initialize();
Smry_Table(y, 2) = smrybio; // in forecast
Smry_Table(y, 3) = smrynum; //sums to accumulate across platoons and settlements
Expand Down
14 changes: 13 additions & 1 deletion SS_expval.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -334,7 +334,7 @@ FUNCTION void Get_expected_values(const int y, const int t);
}
break;
}
case 31: // recruitment deviation #31
case 31: // exp(recruitment deviation) #31
{
if (y >= recdev_start && y <= recdev_end)
{
Expand All @@ -346,6 +346,18 @@ FUNCTION void Get_expected_values(const int y, const int t);
}
break;
}
case 36: // recruitment deviation #36
{
if (y >= recdev_start && y <= recdev_end)
{
vbio = recdev(y);
}
else
{
vbio = 0.0;
}
break;
}
case 32: // recruitment without density-dependence (for pre-recruit survey) #32
{
if (y >= recdev_start && y <= recdev_end)
Expand Down
2 changes: 1 addition & 1 deletion SS_param.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -481,7 +481,7 @@ PARAMETER_SECTION
matrix parm_dev_rwalk(1,N_parm_dev,parm_dev_minyr,parm_dev_maxyr);

init_bounded_number checksum999(998,1000,-999) // set value to 999 to check reading of ss.par
vector timevary_parm(1,timevary_parm_cnt); // will map to the MGparms and selparms that are the actual parameters
vector timevary_parm(1,timevary_parm_cnt); // holds the link parameters; in SS_timevaryparm these are set to actual parms in MGparms, SRparms, Qparms, selparms
matrix parm_timevary(1,timevary_cnt,styr-1,YrMax); // time series of adjusted parm values for block and trend

LOCAL_CALCS
Expand Down
14 changes: 10 additions & 4 deletions SS_readcontrol_330.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1627,7 +1627,7 @@
// timevary_setup(6)=env link type
// timevary_setup(7)=env variable
// timevary_setup(8)=dev vector used
// timevary_setup(9)=dev link type
// timevary_setup(9)=dev link type used in SS_timevarmparm
// timevary_setup(10)=dev min year
// timevary_setup(11)=dev maxyear
// timevary_setup(12)=dev phase
Expand Down Expand Up @@ -4480,6 +4480,10 @@
// clang-format off
END_CALCS

!!// SS_Label_Info_4.9.xx #Create arrays needed for timevary_parameters
vector baseparm_min(1,timevary_parm_cnt)
vector baseparm_max(1,timevary_parm_cnt)

!!// SS_Label_Info_4.9.9 #Create arrays for the total set of selex parameters
vector selparm_LO(1,N_selparm2)
vector selparm_HI(1,N_selparm2)
Expand Down Expand Up @@ -4778,7 +4782,7 @@

ivector parm_dev_type(1,N_parm_dev); // distinguish parameter dev vectors from 2DAR devs
ivector parm_dev_use_rho(1,N_parm_dev); // uses rho parameter, or not
ivector parm_dev_info(1,N_parm_dev); // pointer from list of devvectorsto 2DAR list
ivector parm_dev_info(1,N_parm_dev); // pointer from list of devvectors to 2DAR list
ivector TwoD_AR_ymin(1,TwoD_AR_cnt)
ivector TwoD_AR_ymax(1,TwoD_AR_cnt)
ivector TwoD_AR_amin(1,TwoD_AR_cnt)
Expand Down Expand Up @@ -4806,7 +4810,8 @@
echoinput << " dev vector #: " << k << " setup: " << timevary_setup << " phase: " << parm_dev_PH(k) << endl;
f = timevary_setup(13); // index of base parameter
int picker = timevary_setup(9);
parm_dev_type(k) = 1; // so P'=P+dev*se with objfun using -log(1)
parm_dev_type(k) = 1; // so P'=P+dev*se with objfun using -log(1); so expects se of devs to be approx unit normal
// parm_dev_type is used in SS_objfunc.tpl
if (picker > 20)
{
picker -= 20;
Expand All @@ -4821,7 +4826,7 @@
}
if (picker > 10)
{
parm_dev_type(k) = 3; // use objfun using -log(se) to match 3.30.12 and earlier
parm_dev_type(k) = 3; // P'=P+dev; objfun using -log(se) to match 3.30.12 and earlier
picker -= 10;
}
if (picker == 6) parm_dev_type(k) = 4; // add penalty to keep rmse near 1. Needs to estimate stddev factor
Expand All @@ -4830,6 +4835,7 @@
timevary_def[j](9) = picker; // save in array also

parm_dev_use_rho(k) = 0;
// require rho to be used for some dev approaches
if (picker == 4 || picker == 5 || picker == 6) parm_dev_use_rho(k) = 1;
for (y = parm_dev_minyr(k); y <= parm_dev_maxyr(k); y++)
{
Expand Down
17 changes: 16 additions & 1 deletion SS_readdata_330.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -652,13 +652,14 @@
int Svy_N_rd;
int Svy_N;
init_imatrix Svy_units_rd(1,Nfleet,1,4);
ivector Svy_units(1,Nfleet); //0=num; 1=bio; 2=F; >=30 for special patterns
ivector Svy_units(1,Nfleet); // 0=num; 1=bio; 2=F; >=30 for special patterns
ivector Svy_errtype(1,Nfleet); // -1=normal / 0=lognormal / >0=T
ivector Svy_sdreport(1,Nfleet); // 0=no sdreport; 1=enable sdreport
int Svy_N_sdreport;

LOCAL_CALCS
// clang-format on

data_type = 1; // for surveys
echoinput << "Units: 0=numbers; 1=biomass; 2=F; >=30 for special patterns" << endl;
echoinput << "Errtype: -1=normal; 0=lognormal; >0=T" << endl;
Expand All @@ -669,6 +670,16 @@
Svy_errtype = column(Svy_units_rd, 3);
Svy_sdreport = column(Svy_units_rd, 4);

for (f = 1; f<=Nfleet; f++)
{
if (Svy_units(f) >= 35 && Svy_errtype(f) == 0)
{
warnstream << " survey error type must not be lognormal for surveys of deviations for fleet: " << f << fleetname(f) << endl;
write_message(FATAL, 1);
}
}

// read survey data
ender = 0;
do
{
Expand Down Expand Up @@ -4099,6 +4110,10 @@
<< "# next read flag for specifying selectivity used in forecasts; 0 is value that mimics 3.24, 1 is experimental" << endl;
*(ad_comm::global_datafile) >> Fcast_Specify_Selex;
echoinput << Fcast_Specify_Selex << " # echoed Fcast_Specify_Selex value" << endl;
if (do_densitydependent == 1 && Fcast_Specify_Selex == 0) {
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);
}

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;
Expand Down
1 change: 1 addition & 0 deletions SS_selex.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ FUNCTION void get_selectivity()
else // time-varying
{
sp(j) = parm_timevary(selparm_timevary(Ip + j), y);

if (parm_adjust_method == 1 && (save_for_report > 0 || do_once == 1))
{
if (sp(j) > -999 && (sp(j) < selparm_1(Ip + j, 1) || sp(j) > selparm_1(Ip + j, 2)))
Expand Down
69 changes: 38 additions & 31 deletions SS_timevaryparm.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@
FUNCTION void make_timevaryparm()
{
dvariable baseparm;
double baseparm_min = -999.;
double baseparm_max = 999;
baseparm_min = -999.; // fill array with default
baseparm_max = 999; // fill array with default
dvariable endtrend;
dvariable infl_year;
dvariable slope;
Expand All @@ -36,8 +36,8 @@ FUNCTION void make_timevaryparm()
case 1: // MG
{
baseparm = MGparm(timevary_setup(2)); // index of base parm
baseparm_min = MGparm_LO(timevary_setup(2));
baseparm_max = MGparm_HI(timevary_setup(2));
baseparm_min(tvary) = MGparm_LO(timevary_setup(2));
baseparm_max(tvary) = MGparm_HI(timevary_setup(2));
if (do_once == 1)
echoinput << "base MGparm " << baseparm << endl;
for (j = timevary_setup(3); j < timevary_def[tvary + 1](3); j++)
Expand All @@ -53,8 +53,8 @@ FUNCTION void make_timevaryparm()
case 2: // SR
{
baseparm = SR_parm(timevary_setup(2)); // index of base parm
baseparm_min = SR_parm_LO(timevary_setup(2));
baseparm_max = SR_parm_HI(timevary_setup(2));
baseparm_min(tvary) = SR_parm_LO(timevary_setup(2));
baseparm_max(tvary) = SR_parm_HI(timevary_setup(2));
if (do_once == 1)
echoinput << "base SR_parm " << baseparm << endl;
for (j = timevary_setup(3); j < timevary_def[tvary + 1](3); j++)
Expand All @@ -70,8 +70,8 @@ FUNCTION void make_timevaryparm()
case 3: // Q
{
baseparm = Q_parm(timevary_setup(2)); // index of base parm
baseparm_min = Q_parm_LO(timevary_setup(2));
baseparm_max = Q_parm_HI(timevary_setup(2));
baseparm_min(tvary) = Q_parm_LO(timevary_setup(2));
baseparm_max(tvary) = Q_parm_HI(timevary_setup(2));
if (do_once == 1)
echoinput << "base Qparm " << baseparm << endl;
for (j = timevary_setup(3); j < timevary_def[tvary + 1](3); j++)
Expand All @@ -87,8 +87,8 @@ FUNCTION void make_timevaryparm()
case 5: // selex
{
baseparm = selparm(timevary_setup(2)); // index of base parm
baseparm_min = selparm_LO(timevary_setup(2));
baseparm_max = selparm_HI(timevary_setup(2));
baseparm_min(tvary) = selparm_LO(timevary_setup(2));
baseparm_max(tvary) = selparm_HI(timevary_setup(2));
if (do_once == 1)
echoinput << "base selparm " << baseparm << endl;
for (j = timevary_setup(3); j < timevary_def[tvary + 1](3); j++)
Expand Down Expand Up @@ -160,9 +160,9 @@ FUNCTION void make_timevaryparm()
echoinput << "logistic trend over time " << endl;
if (timevary_setup(4) == -1) // use logistic transform to keep with bounds of the base parameter
{
endtrend = log((baseparm_max - baseparm_min + 0.0000002) / (baseparm - baseparm_min + 0.0000001) - 1.) / (-2.); // transform the base parameter
endtrend = log((baseparm_max(tvary) - baseparm_min(tvary) + 0.0000002) / (baseparm - baseparm_min(tvary) + 0.0000001) - 1.) / (-2.); // transform the base parameter
endtrend += timevary_parm(timevary_parm_cnt); // add the offset Note that offset value is in the transform space
endtrend = baseparm_min + (baseparm_max - baseparm_min) / (1. + mfexp(-2. * endtrend)); // backtransform
endtrend = baseparm_min(tvary) + (baseparm_max(tvary) - baseparm_min(tvary)) / (1. + mfexp(-2. * endtrend)); // backtransform
infl_year = log(0.5) / (-2.); // transform the base parameter
infl_year += timevary_parm(timevary_parm_cnt + 1); // add the offset Note that offset value is in the transform space
infl_year = r_years(styr) + (r_years(endyr) - r_years(styr)) / (1. + mfexp(-2. * infl_year)); // backtransform
Expand All @@ -174,7 +174,7 @@ FUNCTION void make_timevaryparm()
}
else if (timevary_setup(4) == -3) // use parm as fraction of way between bounds
{
endtrend = baseparm_min + (baseparm_max - baseparm_min) * timevary_parm(timevary_parm_cnt);
endtrend = baseparm_min(tvary) + (baseparm_max(tvary) - baseparm_min(tvary)) * timevary_parm(timevary_parm_cnt);
infl_year = r_years(styr) + (r_years(endyr) - r_years(styr)) * timevary_parm(timevary_parm_cnt + 1);
}
slope = timevary_parm(timevary_parm_cnt + 2);
Expand All @@ -197,7 +197,7 @@ FUNCTION void make_timevaryparm()
parm_timevary(tvary, styr - 1) = baseparm;
}

if (timevary_setup(7) > 0) // env link
if (timevary_setup(7) > 0) // env link (negative value indicates density-dependence which is calculated year-by-year in different function)
{
if (do_once == 1)
echoinput << "env_link to env_variable: " << timevary_setup(7) << " using link_type " << timevary_setup(6) << endl;
Expand All @@ -224,13 +224,13 @@ FUNCTION void make_timevaryparm()
case 3: // result constrained by baseparm_min-max; input values are unit normal
{
dvariable temp;
double p_range = baseparm_max - baseparm_min;
double p_range = baseparm_max(tvary) - baseparm_min(tvary);

for (int y1 = env_data_minyr(timevary_setup(7)); y1 <= env_data_maxyr(timevary_setup(7)); y1++)
{
temp = log((parm_timevary(tvary, y1) - baseparm_min + 1.0e-7) / (baseparm_max - parm_timevary(tvary, y1) + 1.0e-7));
temp = log((parm_timevary(tvary, y1) - baseparm_min(tvary) + 1.0e-7) / (baseparm_max(tvary) - parm_timevary(tvary, y1) + 1.0e-7));
temp += timevary_parm(timevary_parm_cnt) * env_data(y1, timevary_setup(7));
parm_timevary(tvary, y1) = baseparm_min + p_range / (1.0 + exp(-temp));
parm_timevary(tvary, y1) = baseparm_min(tvary) + p_range / (1.0 + exp(-temp));
}
timevary_parm_cnt++;
break;
Expand All @@ -255,7 +255,7 @@ FUNCTION void make_timevaryparm()
echoinput << "dev vector #: " << k << endl;
parm_dev_stddev(k) = timevary_parm(timevary_parm_cnt);
parm_dev_rho(k) = timevary_parm(timevary_parm_cnt + 1);
int picker = timevary_setup(9);
int picker = timevary_setup(9); // selects the method for creating time-vary parameter from dev vector

switch (picker)
{
Expand Down Expand Up @@ -314,19 +314,19 @@ FUNCTION void make_timevaryparm()
{
// NOTE: if the stddev parameter is greater than 1.8, the distribution of adjusted parameters will become U-shaped
dvariable temp;
double p_range = baseparm_max - baseparm_min;
double p_range = baseparm_max(tvary) - baseparm_min(tvary);
int j = timevary_setup(10);
parm_dev_rwalk(k, j) = parm_dev(k, j) * parm_dev_stddev(k); // 1st yr dev
// p_base=(parm_timevary(tvary,j)-baseparm_min)/(baseparm_max-baseparm_min); // convert parm to (0,1) scale
// p_base=(parm_timevary(tvary,j)-baseparm_min(tvary))/(baseparm_max(tvary)-baseparm_min(tvary)); // convert parm to (0,1) scale
// temp=log(p_base/(1.-p_base)) + parm_dev_rwalk(k,j); // convert to logit and add dev; so dev must be in units of the logit
temp = log((parm_timevary(tvary, j) - baseparm_min + 1.0e-7) / (baseparm_max - parm_timevary(tvary, j) + 1.0e-7));
parm_timevary(tvary, j) = baseparm_min + p_range / (1.0 + exp(-temp - parm_dev_rwalk(k, j)));
temp = log((parm_timevary(tvary, j) - baseparm_min(tvary) + 1.0e-7) / (baseparm_max(tvary) - parm_timevary(tvary, j) + 1.0e-7));
parm_timevary(tvary, j) = baseparm_min(tvary) + p_range / (1.0 + exp(-temp - parm_dev_rwalk(k, j)));
for (j = timevary_setup(10) + 1; j <= timevary_setup(11); j++)
{
// =(1-rho)*mean + rho*prevval + dev // where mean = 0.0
parm_dev_rwalk(k, j) = parm_dev_rho(k) * parm_dev_rwalk(k, j - 1) + parm_dev(k, j) * parm_dev_stddev(k); // update MRRW using annual dev
temp = log((parm_timevary(tvary, j) - baseparm_min + 1.0e-7) / (baseparm_max - parm_timevary(tvary, j) + 1.0e-7));
parm_timevary(tvary, j) = baseparm_min + p_range / (1.0 + exp(-temp - parm_dev_rwalk(k, j)));
temp = log((parm_timevary(tvary, j) - baseparm_min(tvary) + 1.0e-7) / (baseparm_max(tvary) - parm_timevary(tvary, j) + 1.0e-7));
parm_timevary(tvary, j) = baseparm_min(tvary) + p_range / (1.0 + exp(-temp - parm_dev_rwalk(k, j)));
}
break;
}
Expand All @@ -351,29 +351,36 @@ FUNCTION void make_densitydependent_parm(int const y1)
timevary_setup(1, 13) = timevary_def[tvary](1, 13);
if (timevary_setup(7) < 0) // density-dependent
{
int env_var = timevary_setup(7);
timevary_parm_cnt = timevary_setup(3); // link parameter index
if (do_once == 1)
echoinput << y1 << " density-dependent to env_variable: " << timevary_setup(7) << " using link_type " << timevary_setup(6) << " env: " << env_data(y1, timevary_setup(7)) << " parm: " << timevary_parm(timevary_parm_cnt) << endl;
echoinput << y1 << " density-dependent to env_variable: " << env_var << " using link_type "
<< timevary_setup(6) << " env: " << env_data(y1, env_var) << " parm: " << timevary_parm(timevary_parm_cnt) << endl;
switch (int(timevary_setup(6)))
{
case 1: // exponential env link
{
parm_timevary(tvary, y1) *= mfexp(timevary_parm(timevary_parm_cnt) * (env_data(y1, timevary_setup(7))));
parm_timevary(tvary, y1) *= mfexp(timevary_parm(timevary_parm_cnt) * env_data(y1, env_var));
break;
}
case 2: // linear env link
{
parm_timevary(tvary, y1) += timevary_parm(timevary_parm_cnt) * env_data(y1, timevary_setup(7));
parm_timevary(tvary, y1) += timevary_parm(timevary_parm_cnt) * env_data(y1, env_var);
break;
}
case 3:
case 3: // result constrained by baseparm_min-max; input values are unit normal
{
// not implemented
dvariable temp;
double p_range = baseparm_max(tvary) - baseparm_min(tvary);
temp = log((parm_timevary(tvary, y1) - baseparm_min(tvary) + 1.0e-7) / (baseparm_max(tvary) - parm_timevary(tvary, y1) + 1.0e-7));
temp += timevary_parm(timevary_parm_cnt) * env_data(y1, env_var);
parm_timevary(tvary, y1) = baseparm_min(tvary) + p_range / (1.0 + exp(-temp));
break;
}
case 4: // logistic MGparm env link
case 4: // logistic env link
{
// first parm is offset ; second is slope
parm_timevary(tvary, y1) = 2.00000 / (1.00000 + mfexp(-timevary_parm(timevary_parm_cnt + 1) * (env_data(y1, timevary_setup(7)) - timevary_parm(timevary_parm_cnt))));
parm_timevary(tvary, y1) = 2.00000 / (1.00000 + mfexp(-timevary_parm(timevary_parm_cnt + 1) * (env_data(y1, env_var) - timevary_parm(timevary_parm_cnt))));
break;
}
}
Expand Down
5 changes: 3 additions & 2 deletions SS_write_ssnew.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,10 +143,11 @@ FUNCTION void write_nudata()
report1 << "-9999 0 0 0 0" << endl
<< "#" << endl;

report1 << " #_CPUE_and_surveyabundance_observations" << endl;
report1 << "#_Units: 0=numbers; 1=biomass; 2=F; 30=spawnbio; 31=recdev; 32=spawnbio*recdev; 33=recruitment; 34=depletion(&see Qsetup); 35=parm_dev(&see Qsetup)" << endl;
report1 << "#_CPUE_and_surveyabundance_and_index_observations" << endl;
report1 << "#_Units: 0=numbers; 1=biomass; 2=F; 30=spawnbio; 31=exp(recdev); 36=recdev; 32=spawnbio*recdev; 33=recruitment; 34=depletion(&see Qsetup); 35=parm_dev(&see Qsetup)" << endl;
report1 << "#_Errtype: -1=normal; 0=lognormal; >0=T" << endl;
report1 << "#_SD_Report: 0=no sdreport; 1=enable sdreport" << endl;
report1 << "#_note that link functions are specified in Q_setup section of control file" << endl;
report1 << "#_Fleet Units Errtype SD_Report" << endl;
for (f = 1; f <= Nfleet; f++)
report1 << f << " " << Svy_units(f) << " " << Svy_errtype(f) << " " << Svy_sdreport(f) << " # " << fleetname(f) << endl;
Expand Down

0 comments on commit 6b99f1f

Please sign in to comment.