From b4da4755d0f44a05331660cf8af981959fdf4b36 Mon Sep 17 00:00:00 2001 From: RickMethot Date: Wed, 21 Dec 2022 16:53:45 -0800 Subject: [PATCH 1/2] fix #165 for bounded link enable use of density-dependence with a bounded link env. --- SS_benchfore.tpl | 1 - SS_readcontrol_330.tpl | 14 +++++-- SS_readdata_330.tpl | 4 ++ SS_selex.tpl | 1 + SS_timevaryparm.tpl | 69 +++++++++++++++++++---------------- StockSynthesis.code-workspace | 4 +- 6 files changed, 56 insertions(+), 37 deletions(-) diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl index 8521d215..1aae3cad 100644 --- a/SS_benchfore.tpl +++ b/SS_benchfore.tpl @@ -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 diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl index 588a198a..f4f081ce 100644 --- a/SS_readcontrol_330.tpl +++ b/SS_readcontrol_330.tpl @@ -1617,7 +1617,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 @@ -4468,6 +4468,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) @@ -4766,7 +4770,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) @@ -4794,7 +4798,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; @@ -4809,7 +4814,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 @@ -4818,6 +4823,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++) { diff --git a/SS_readdata_330.tpl b/SS_readdata_330.tpl index 9f895a4a..7bc8aede 100644 --- a/SS_readdata_330.tpl +++ b/SS_readdata_330.tpl @@ -4029,6 +4029,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 "< 0 || do_once == 1)) { if (sp(j) > -999 && (sp(j) < selparm_1(Ip + j, 1) || sp(j) > selparm_1(Ip + j, 2))) diff --git a/SS_timevaryparm.tpl b/SS_timevaryparm.tpl index fdf6c696..3c721455 100644 --- a/SS_timevaryparm.tpl +++ b/SS_timevaryparm.tpl @@ -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; @@ -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++) @@ -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++) @@ -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++) @@ -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++) @@ -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 @@ -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); @@ -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; @@ -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; @@ -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) { @@ -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; } @@ -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; } } diff --git a/StockSynthesis.code-workspace b/StockSynthesis.code-workspace index 22aacea5..a10c413f 100644 --- a/StockSynthesis.code-workspace +++ b/StockSynthesis.code-workspace @@ -8,7 +8,9 @@ "files.associations": { "*.tpl": "c", "\"*.extension\":": "\"tpl\"", - "*.htp": "c" + "*.htp": "c", + "cmath": "c", + "ostream": "c" } } } \ No newline at end of file From 46e958757060dff31be8b7e71def3ab41bcf7331 Mon Sep 17 00:00:00 2001 From: RickMethot Date: Thu, 22 Dec 2022 14:03:40 -0800 Subject: [PATCH 2/2] add survey type 36 add survey type 36 as recruitment dev --- SS_expval.tpl | 14 +++++++++++++- SS_param.tpl | 2 +- SS_readdata_330.tpl | 13 ++++++++++++- SS_write_ssnew.tpl | 5 +++-- 4 files changed, 29 insertions(+), 5 deletions(-) diff --git a/SS_expval.tpl b/SS_expval.tpl index 4c788c63..43388fde 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -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) { @@ -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) diff --git a/SS_param.tpl b/SS_param.tpl index 2dc1dbbf..119862c2 100644 --- a/SS_param.tpl +++ b/SS_param.tpl @@ -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 diff --git a/SS_readdata_330.tpl b/SS_readdata_330.tpl index 7bc8aede..90c4e890 100644 --- a/SS_readdata_330.tpl +++ b/SS_readdata_330.tpl @@ -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; @@ -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 { diff --git a/SS_write_ssnew.tpl b/SS_write_ssnew.tpl index 08a86de6..7ed38ddd 100644 --- a/SS_write_ssnew.tpl +++ b/SS_write_ssnew.tpl @@ -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;