Skip to content

Commit 6b99f1f

Browse files
Rick-Methot-NOAARick-Methot-NOAA
authored
Merge pull request #391 from nmfs-stock-synthesis/improve-time-varying-parameter-implementation
fix #165 for bounded link; add survey type #41
2 parents f7b1bda + e934175 commit 6b99f1f

File tree

8 files changed

+82
-41
lines changed

8 files changed

+82
-41
lines changed

‎SS_benchfore.tpl‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2204,7 +2204,6 @@ FUNCTION void Get_Forecast()
22042204
}
22052205
env_data(y, -3) = log(smrybio / Smry_Table(styr - 1, 2));
22062206
env_data(y, -4) = log(smrynum / Smry_Table(styr - 1, 3));
2207-
22082207
Smry_Table(y).initialize();
22092208
Smry_Table(y, 2) = smrybio; // in forecast
22102209
Smry_Table(y, 3) = smrynum; //sums to accumulate across platoons and settlements

‎SS_expval.tpl‎

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -334,7 +334,7 @@ FUNCTION void Get_expected_values(const int y, const int t);
334334
}
335335
break;
336336
}
337-
case 31: // recruitment deviation #31
337+
case 31: // exp(recruitment deviation) #31
338338
{
339339
if (y >= recdev_start && y <= recdev_end)
340340
{
@@ -346,6 +346,18 @@ FUNCTION void Get_expected_values(const int y, const int t);
346346
}
347347
break;
348348
}
349+
case 36: // recruitment deviation #36
350+
{
351+
if (y >= recdev_start && y <= recdev_end)
352+
{
353+
vbio = recdev(y);
354+
}
355+
else
356+
{
357+
vbio = 0.0;
358+
}
359+
break;
360+
}
349361
case 32: // recruitment without density-dependence (for pre-recruit survey) #32
350362
{
351363
if (y >= recdev_start && y <= recdev_end)

‎SS_param.tpl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -481,7 +481,7 @@ PARAMETER_SECTION
481481
matrix parm_dev_rwalk(1,N_parm_dev,parm_dev_minyr,parm_dev_maxyr);
482482

483483
init_bounded_number checksum999(998,1000,-999) // set value to 999 to check reading of ss.par
484-
vector timevary_parm(1,timevary_parm_cnt); // will map to the MGparms and selparms that are the actual parameters
484+
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
485485
matrix parm_timevary(1,timevary_cnt,styr-1,YrMax); // time series of adjusted parm values for block and trend
486486

487487
LOCAL_CALCS

‎SS_readcontrol_330.tpl‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1627,7 +1627,7 @@
16271627
// timevary_setup(6)=env link type
16281628
// timevary_setup(7)=env variable
16291629
// timevary_setup(8)=dev vector used
1630-
// timevary_setup(9)=dev link type
1630+
// timevary_setup(9)=dev link type used in SS_timevarmparm
16311631
// timevary_setup(10)=dev min year
16321632
// timevary_setup(11)=dev maxyear
16331633
// timevary_setup(12)=dev phase
@@ -4480,6 +4480,10 @@
44804480
// clang-format off
44814481
END_CALCS
44824482
4483+
!!// SS_Label_Info_4.9.xx #Create arrays needed for timevary_parameters
4484+
vector baseparm_min(1,timevary_parm_cnt)
4485+
vector baseparm_max(1,timevary_parm_cnt)
4486+
44834487
!!// SS_Label_Info_4.9.9 #Create arrays for the total set of selex parameters
44844488
vector selparm_LO(1,N_selparm2)
44854489
vector selparm_HI(1,N_selparm2)
@@ -4778,7 +4782,7 @@
47784782
47794783
ivector parm_dev_type(1,N_parm_dev); // distinguish parameter dev vectors from 2DAR devs
47804784
ivector parm_dev_use_rho(1,N_parm_dev); // uses rho parameter, or not
4781-
ivector parm_dev_info(1,N_parm_dev); // pointer from list of devvectorsto 2DAR list
4785+
ivector parm_dev_info(1,N_parm_dev); // pointer from list of devvectors to 2DAR list
47824786
ivector TwoD_AR_ymin(1,TwoD_AR_cnt)
47834787
ivector TwoD_AR_ymax(1,TwoD_AR_cnt)
47844788
ivector TwoD_AR_amin(1,TwoD_AR_cnt)
@@ -4806,7 +4810,8 @@
48064810
echoinput << " dev vector #: " << k << " setup: " << timevary_setup << " phase: " << parm_dev_PH(k) << endl;
48074811
f = timevary_setup(13); // index of base parameter
48084812
int picker = timevary_setup(9);
4809-
parm_dev_type(k) = 1; // so P'=P+dev*se with objfun using -log(1)
4813+
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
4814+
// parm_dev_type is used in SS_objfunc.tpl
48104815
if (picker > 20)
48114816
{
48124817
picker -= 20;
@@ -4821,7 +4826,7 @@
48214826
}
48224827
if (picker > 10)
48234828
{
4824-
parm_dev_type(k) = 3; // use objfun using -log(se) to match 3.30.12 and earlier
4829+
parm_dev_type(k) = 3; // P'=P+dev; objfun using -log(se) to match 3.30.12 and earlier
48254830
picker -= 10;
48264831
}
48274832
if (picker == 6) parm_dev_type(k) = 4; // add penalty to keep rmse near 1. Needs to estimate stddev factor
@@ -4830,6 +4835,7 @@
48304835
timevary_def[j](9) = picker; // save in array also
48314836
48324837
parm_dev_use_rho(k) = 0;
4838+
// require rho to be used for some dev approaches
48334839
if (picker == 4 || picker == 5 || picker == 6) parm_dev_use_rho(k) = 1;
48344840
for (y = parm_dev_minyr(k); y <= parm_dev_maxyr(k); y++)
48354841
{

‎SS_readdata_330.tpl‎

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -652,13 +652,14 @@
652652
int Svy_N_rd;
653653
int Svy_N;
654654
init_imatrix Svy_units_rd(1,Nfleet,1,4);
655-
ivector Svy_units(1,Nfleet); //0=num; 1=bio; 2=F; >=30 for special patterns
655+
ivector Svy_units(1,Nfleet); // 0=num; 1=bio; 2=F; >=30 for special patterns
656656
ivector Svy_errtype(1,Nfleet); // -1=normal / 0=lognormal / >0=T
657657
ivector Svy_sdreport(1,Nfleet); // 0=no sdreport; 1=enable sdreport
658658
int Svy_N_sdreport;
659659
660660
LOCAL_CALCS
661661
// clang-format on
662+
662663
data_type = 1; // for surveys
663664
echoinput << "Units: 0=numbers; 1=biomass; 2=F; >=30 for special patterns" << endl;
664665
echoinput << "Errtype: -1=normal; 0=lognormal; >0=T" << endl;
@@ -669,6 +670,16 @@
669670
Svy_errtype = column(Svy_units_rd, 3);
670671
Svy_sdreport = column(Svy_units_rd, 4);
671672

673+
for (f = 1; f<=Nfleet; f++)
674+
{
675+
if (Svy_units(f) >= 35 && Svy_errtype(f) == 0)
676+
{
677+
warnstream << " survey error type must not be lognormal for surveys of deviations for fleet: " << f << fleetname(f) << endl;
678+
write_message(FATAL, 1);
679+
}
680+
}
681+
682+
// read survey data
672683
ender = 0;
673684
do
674685
{
@@ -4099,6 +4110,10 @@
40994110
<< "# next read flag for specifying selectivity used in forecasts; 0 is value that mimics 3.24, 1 is experimental" << endl;
41004111
*(ad_comm::global_datafile) >> Fcast_Specify_Selex;
41014112
echoinput << Fcast_Specify_Selex << " # echoed Fcast_Specify_Selex value" << endl;
4113+
if (do_densitydependent == 1 && Fcast_Specify_Selex == 0) {
4114+
warnstream << "Fcast_specify_selex is 0 but user should change to 1 if density dependence affects a selectivity parameter or growth "<<endl;
4115+
write_message(WARN, 0);
4116+
}
41024117
41034118
echoinput << endl
41044119
<< "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;

‎SS_selex.tpl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@ FUNCTION void get_selectivity()
5050
else // time-varying
5151
{
5252
sp(j) = parm_timevary(selparm_timevary(Ip + j), y);
53+
5354
if (parm_adjust_method == 1 && (save_for_report > 0 || do_once == 1))
5455
{
5556
if (sp(j) > -999 && (sp(j) < selparm_1(Ip + j, 1) || sp(j) > selparm_1(Ip + j, 2)))

‎SS_timevaryparm.tpl‎

Lines changed: 38 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,8 @@
88
FUNCTION void make_timevaryparm()
99
{
1010
dvariable baseparm;
11-
double baseparm_min = -999.;
12-
double baseparm_max = 999;
11+
baseparm_min = -999.; // fill array with default
12+
baseparm_max = 999; // fill array with default
1313
dvariable endtrend;
1414
dvariable infl_year;
1515
dvariable slope;
@@ -36,8 +36,8 @@ FUNCTION void make_timevaryparm()
3636
case 1: // MG
3737
{
3838
baseparm = MGparm(timevary_setup(2)); // index of base parm
39-
baseparm_min = MGparm_LO(timevary_setup(2));
40-
baseparm_max = MGparm_HI(timevary_setup(2));
39+
baseparm_min(tvary) = MGparm_LO(timevary_setup(2));
40+
baseparm_max(tvary) = MGparm_HI(timevary_setup(2));
4141
if (do_once == 1)
4242
echoinput << "base MGparm " << baseparm << endl;
4343
for (j = timevary_setup(3); j < timevary_def[tvary + 1](3); j++)
@@ -53,8 +53,8 @@ FUNCTION void make_timevaryparm()
5353
case 2: // SR
5454
{
5555
baseparm = SR_parm(timevary_setup(2)); // index of base parm
56-
baseparm_min = SR_parm_LO(timevary_setup(2));
57-
baseparm_max = SR_parm_HI(timevary_setup(2));
56+
baseparm_min(tvary) = SR_parm_LO(timevary_setup(2));
57+
baseparm_max(tvary) = SR_parm_HI(timevary_setup(2));
5858
if (do_once == 1)
5959
echoinput << "base SR_parm " << baseparm << endl;
6060
for (j = timevary_setup(3); j < timevary_def[tvary + 1](3); j++)
@@ -70,8 +70,8 @@ FUNCTION void make_timevaryparm()
7070
case 3: // Q
7171
{
7272
baseparm = Q_parm(timevary_setup(2)); // index of base parm
73-
baseparm_min = Q_parm_LO(timevary_setup(2));
74-
baseparm_max = Q_parm_HI(timevary_setup(2));
73+
baseparm_min(tvary) = Q_parm_LO(timevary_setup(2));
74+
baseparm_max(tvary) = Q_parm_HI(timevary_setup(2));
7575
if (do_once == 1)
7676
echoinput << "base Qparm " << baseparm << endl;
7777
for (j = timevary_setup(3); j < timevary_def[tvary + 1](3); j++)
@@ -87,8 +87,8 @@ FUNCTION void make_timevaryparm()
8787
case 5: // selex
8888
{
8989
baseparm = selparm(timevary_setup(2)); // index of base parm
90-
baseparm_min = selparm_LO(timevary_setup(2));
91-
baseparm_max = selparm_HI(timevary_setup(2));
90+
baseparm_min(tvary) = selparm_LO(timevary_setup(2));
91+
baseparm_max(tvary) = selparm_HI(timevary_setup(2));
9292
if (do_once == 1)
9393
echoinput << "base selparm " << baseparm << endl;
9494
for (j = timevary_setup(3); j < timevary_def[tvary + 1](3); j++)
@@ -160,9 +160,9 @@ FUNCTION void make_timevaryparm()
160160
echoinput << "logistic trend over time " << endl;
161161
if (timevary_setup(4) == -1) // use logistic transform to keep with bounds of the base parameter
162162
{
163-
endtrend = log((baseparm_max - baseparm_min + 0.0000002) / (baseparm - baseparm_min + 0.0000001) - 1.) / (-2.); // transform the base parameter
163+
endtrend = log((baseparm_max(tvary) - baseparm_min(tvary) + 0.0000002) / (baseparm - baseparm_min(tvary) + 0.0000001) - 1.) / (-2.); // transform the base parameter
164164
endtrend += timevary_parm(timevary_parm_cnt); // add the offset Note that offset value is in the transform space
165-
endtrend = baseparm_min + (baseparm_max - baseparm_min) / (1. + mfexp(-2. * endtrend)); // backtransform
165+
endtrend = baseparm_min(tvary) + (baseparm_max(tvary) - baseparm_min(tvary)) / (1. + mfexp(-2. * endtrend)); // backtransform
166166
infl_year = log(0.5) / (-2.); // transform the base parameter
167167
infl_year += timevary_parm(timevary_parm_cnt + 1); // add the offset Note that offset value is in the transform space
168168
infl_year = r_years(styr) + (r_years(endyr) - r_years(styr)) / (1. + mfexp(-2. * infl_year)); // backtransform
@@ -174,7 +174,7 @@ FUNCTION void make_timevaryparm()
174174
}
175175
else if (timevary_setup(4) == -3) // use parm as fraction of way between bounds
176176
{
177-
endtrend = baseparm_min + (baseparm_max - baseparm_min) * timevary_parm(timevary_parm_cnt);
177+
endtrend = baseparm_min(tvary) + (baseparm_max(tvary) - baseparm_min(tvary)) * timevary_parm(timevary_parm_cnt);
178178
infl_year = r_years(styr) + (r_years(endyr) - r_years(styr)) * timevary_parm(timevary_parm_cnt + 1);
179179
}
180180
slope = timevary_parm(timevary_parm_cnt + 2);
@@ -197,7 +197,7 @@ FUNCTION void make_timevaryparm()
197197
parm_timevary(tvary, styr - 1) = baseparm;
198198
}
199199

200-
if (timevary_setup(7) > 0) // env link
200+
if (timevary_setup(7) > 0) // env link (negative value indicates density-dependence which is calculated year-by-year in different function)
201201
{
202202
if (do_once == 1)
203203
echoinput << "env_link to env_variable: " << timevary_setup(7) << " using link_type " << timevary_setup(6) << endl;
@@ -224,13 +224,13 @@ FUNCTION void make_timevaryparm()
224224
case 3: // result constrained by baseparm_min-max; input values are unit normal
225225
{
226226
dvariable temp;
227-
double p_range = baseparm_max - baseparm_min;
227+
double p_range = baseparm_max(tvary) - baseparm_min(tvary);
228228
229229
for (int y1 = env_data_minyr(timevary_setup(7)); y1 <= env_data_maxyr(timevary_setup(7)); y1++)
230230
{
231-
temp = log((parm_timevary(tvary, y1) - baseparm_min + 1.0e-7) / (baseparm_max - parm_timevary(tvary, y1) + 1.0e-7));
231+
temp = log((parm_timevary(tvary, y1) - baseparm_min(tvary) + 1.0e-7) / (baseparm_max(tvary) - parm_timevary(tvary, y1) + 1.0e-7));
232232
temp += timevary_parm(timevary_parm_cnt) * env_data(y1, timevary_setup(7));
233-
parm_timevary(tvary, y1) = baseparm_min + p_range / (1.0 + exp(-temp));
233+
parm_timevary(tvary, y1) = baseparm_min(tvary) + p_range / (1.0 + exp(-temp));
234234
}
235235
timevary_parm_cnt++;
236236
break;
@@ -255,7 +255,7 @@ FUNCTION void make_timevaryparm()
255255
echoinput << "dev vector #: " << k << endl;
256256
parm_dev_stddev(k) = timevary_parm(timevary_parm_cnt);
257257
parm_dev_rho(k) = timevary_parm(timevary_parm_cnt + 1);
258-
int picker = timevary_setup(9);
258+
int picker = timevary_setup(9); // selects the method for creating time-vary parameter from dev vector
259259
260260
switch (picker)
261261
{
@@ -314,19 +314,19 @@ FUNCTION void make_timevaryparm()
314314
{
315315
// NOTE: if the stddev parameter is greater than 1.8, the distribution of adjusted parameters will become U-shaped
316316
dvariable temp;
317-
double p_range = baseparm_max - baseparm_min;
317+
double p_range = baseparm_max(tvary) - baseparm_min(tvary);
318318
int j = timevary_setup(10);
319319
parm_dev_rwalk(k, j) = parm_dev(k, j) * parm_dev_stddev(k); // 1st yr dev
320-
// p_base=(parm_timevary(tvary,j)-baseparm_min)/(baseparm_max-baseparm_min); // convert parm to (0,1) scale
320+
// p_base=(parm_timevary(tvary,j)-baseparm_min(tvary))/(baseparm_max(tvary)-baseparm_min(tvary)); // convert parm to (0,1) scale
321321
// 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
322-
temp = log((parm_timevary(tvary, j) - baseparm_min + 1.0e-7) / (baseparm_max - parm_timevary(tvary, j) + 1.0e-7));
323-
parm_timevary(tvary, j) = baseparm_min + p_range / (1.0 + exp(-temp - parm_dev_rwalk(k, j)));
322+
temp = log((parm_timevary(tvary, j) - baseparm_min(tvary) + 1.0e-7) / (baseparm_max(tvary) - parm_timevary(tvary, j) + 1.0e-7));
323+
parm_timevary(tvary, j) = baseparm_min(tvary) + p_range / (1.0 + exp(-temp - parm_dev_rwalk(k, j)));
324324
for (j = timevary_setup(10) + 1; j <= timevary_setup(11); j++)
325325
{
326326
// =(1-rho)*mean + rho*prevval + dev // where mean = 0.0
327327
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
328-
temp = log((parm_timevary(tvary, j) - baseparm_min + 1.0e-7) / (baseparm_max - parm_timevary(tvary, j) + 1.0e-7));
329-
parm_timevary(tvary, j) = baseparm_min + p_range / (1.0 + exp(-temp - parm_dev_rwalk(k, j)));
328+
temp = log((parm_timevary(tvary, j) - baseparm_min(tvary) + 1.0e-7) / (baseparm_max(tvary) - parm_timevary(tvary, j) + 1.0e-7));
329+
parm_timevary(tvary, j) = baseparm_min(tvary) + p_range / (1.0 + exp(-temp - parm_dev_rwalk(k, j)));
330330
}
331331
break;
332332
}
@@ -351,29 +351,36 @@ FUNCTION void make_densitydependent_parm(int const y1)
351351
timevary_setup(1, 13) = timevary_def[tvary](1, 13);
352352
if (timevary_setup(7) < 0) // density-dependent
353353
{
354+
int env_var = timevary_setup(7);
354355
timevary_parm_cnt = timevary_setup(3); // link parameter index
355356
if (do_once == 1)
356-
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;
357+
echoinput << y1 << " density-dependent to env_variable: " << env_var << " using link_type "
358+
<< timevary_setup(6) << " env: " << env_data(y1, env_var) << " parm: " << timevary_parm(timevary_parm_cnt) << endl;
357359
switch (int(timevary_setup(6)))
358360
{
359361
case 1: // exponential env link
360362
{
361-
parm_timevary(tvary, y1) *= mfexp(timevary_parm(timevary_parm_cnt) * (env_data(y1, timevary_setup(7))));
363+
parm_timevary(tvary, y1) *= mfexp(timevary_parm(timevary_parm_cnt) * env_data(y1, env_var));
362364
break;
363365
}
364366
case 2: // linear env link
365367
{
366-
parm_timevary(tvary, y1) += timevary_parm(timevary_parm_cnt) * env_data(y1, timevary_setup(7));
368+
parm_timevary(tvary, y1) += timevary_parm(timevary_parm_cnt) * env_data(y1, env_var);
367369
break;
368370
}
369-
case 3:
371+
case 3: // result constrained by baseparm_min-max; input values are unit normal
370372
{
371-
// not implemented
373+
dvariable temp;
374+
double p_range = baseparm_max(tvary) - baseparm_min(tvary);
375+
temp = log((parm_timevary(tvary, y1) - baseparm_min(tvary) + 1.0e-7) / (baseparm_max(tvary) - parm_timevary(tvary, y1) + 1.0e-7));
376+
temp += timevary_parm(timevary_parm_cnt) * env_data(y1, env_var);
377+
parm_timevary(tvary, y1) = baseparm_min(tvary) + p_range / (1.0 + exp(-temp));
378+
break;
372379
}
373-
case 4: // logistic MGparm env link
380+
case 4: // logistic env link
374381
{
375382
// first parm is offset ; second is slope
376-
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))));
383+
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))));
377384
break;
378385
}
379386
}

‎SS_write_ssnew.tpl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -143,10 +143,11 @@ FUNCTION void write_nudata()
143143
report1 << "-9999 0 0 0 0" << endl
144144
<< "#" << endl;
145145

146-
report1 << " #_CPUE_and_surveyabundance_observations" << endl;
147-
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;
146+
report1 << "#_CPUE_and_surveyabundance_and_index_observations" << endl;
147+
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;
148148
report1 << "#_Errtype: -1=normal; 0=lognormal; >0=T" << endl;
149149
report1 << "#_SD_Report: 0=no sdreport; 1=enable sdreport" << endl;
150+
report1 << "#_note that link functions are specified in Q_setup section of control file" << endl;
150151
report1 << "#_Fleet Units Errtype SD_Report" << endl;
151152
for (f = 1; f <= Nfleet; f++)
152153
report1 << f << " " << Svy_units(f) << " " << Svy_errtype(f) << " " << Svy_sdreport(f) << " # " << fleetname(f) << endl;

0 commit comments

Comments
 (0)