diff --git a/SS_biofxn.tpl b/SS_biofxn.tpl index b9781597..0839e29e 100644 --- a/SS_biofxn.tpl +++ b/SS_biofxn.tpl @@ -94,15 +94,14 @@ FUNCTION void get_growth1() { VBK_seas = sum(seasdur); // set vector to null effect } + if (do_once == 1) + echoinput << "season_duration_as_used_in_growth_calculations: " << VBK_seas(1, nseas) << endl << "summed_for_annual_calcs: " << VBK_seas(0)< 0) { wtlen_seas(0, gp, j) = 0.0; @@ -127,7 +126,6 @@ FUNCTION void get_growth1() for (s = 1; s <= nseas; s++) wtlen_seas(s) = 1.0; // set vector to null effect } - // SS_Label_Info_15.2 #create variability of size-at-age factors using direct assignment or offset approaches gp = 0; for (gg = 1; gg <= gender; gg++) @@ -206,19 +204,25 @@ FUNCTION void get_growth2(const int y) // called at beginning of each year, so y is known // if y=styr, then does equilibrium size-at-age according to start year growth parameters // for any year, calculates for each season the size at the beginning of the next season, with growth increment calculated according to that year's parameters + // VBK(gp,a) is the growth parameter for growth type gp and age. VBK for age 0 is used for all ages. set to negative of parameter + // VBK_seas(0, nseas) is optional seasonal adjustment to VBK + // VBK_by_seas is VBK adjusted by VBK_seas + // VBK_work is current value of VBK for the selected model condition + //Growth Cessation Model code added by Mark Maunder October 2018 //The growth cessation model is described in //Maunder, M.N., Deriso, R.B., Schaefer, K.M., Fuller, D.W., Aires-da-Silva, A.M., Minte‑Vera, C.V., Campana, S.E. 2018. The growth cessation model: a growth model for species showing a near cessation in growth with application to bigeye tuna (Thunnus obesus). Marine Biology (2018) 165:76. - //Ian Taylor derived the formula for Linf + //Ian Taylor derived the formula for L_inf int k2; int add_age; int ALK_idx2; // beginning of first subseas of next season + int L_inf_plus = 0; dvariable plusgroupsize; dvariable current_size; - dvariable VBK_temp; - dvariable VBK_temp2; // with VBKseas(s) multiplied + dvariable VBK_work; // working version of VBK + dvariable VBK_by_seas; // with VBKseas(s) multiplied dvariable LminR; dvariable LmaxR; dvariable LinfR; @@ -240,7 +244,7 @@ FUNCTION void get_growth2(const int y) #ifdef DO_ONCE { if (do_once == 1) - echoinput << "GROWTH, yr= " << y << endl; + echoinput << "GROWTH2 by season in year " << y << endl; } #endif for (gg = 1; gg <= gender; gg++) @@ -248,32 +252,36 @@ FUNCTION void get_growth2(const int y) { gp++; Ip = MGparm_point(gg, GPat) + N_natMparms; - switch (Grow_type) // create specific growth parameters from the mgp_adj list of current MGparms + // SS_Label_Info_16.2.1 #set Lmin, Lmax, VBK, Richards to this year's values for mgp_adj + // SS_Label_Info_16.2.2 #Set up age specific k + if (Grow_type != 7) { - case 7: // empirical length + // get basic parameters used by all grow_type + if (MGparm_def > 1 && gp > 1) // do offset approach { - break; + Lmin(gp) = Lmin(1) * mfexp(mgp_adj(Ip)); + Lmax_temp(gp) = Lmax_temp(1) * mfexp(mgp_adj(Ip + 1)); + VBK(gp) = VBK(1) * mfexp(mgp_adj(Ip + 2)); // assigns to all ages for which VBK is defined } + else + { + Lmin(gp) = mgp_adj(Ip); + Lmax_temp(gp) = mgp_adj(Ip + 1); // size at A2; could be 999 to indicate L_inf + VBK(gp) = -mgp_adj(Ip + 2); // because always used as negative; assigns to all ages for which VBK is defined + } + L_inf(gp) = Lmax_temp(gp); // because this is the typical assignment + VBK_work = VBK(gp, 0); // will be reset to VBK(gp,nages) if using age-specific K; note that this inherits the negative sign - default: // process parameters for all other grow_type + // get specific growth parameters: Richards/Gompertz, growth cessation, and age-specific VBK + switch (Grow_type) { - // SS_Label_Info_16.2.1 #set Lmin, Lmax, VBK, Richards to this year's values for mgp_adj - if (MGparm_def > 1 && gp > 1) // do offset approach - { - Lmin(gp) = Lmin(1) * mfexp(mgp_adj(Ip)); - Lmax_temp(gp) = Lmax_temp(1) * mfexp(mgp_adj(Ip + 1)); - VBK(gp) = VBK(1) * mfexp(mgp_adj(Ip + 2)); // assigns to all ages for which VBK is defined - } - else + case 1: // standard von Bertallanfy no action needed { - Lmin(gp) = mgp_adj(Ip); - Lmax_temp(gp) = mgp_adj(Ip + 1); // size at A2; could be 999 to indicate Linf - VBK(gp) = -mgp_adj(Ip + 2); // because always used as negative; assigns to all ages for which VBK is defined + if (AFIX2 != 999) L_inf_plus = 1; + break; } - VBK_temp = VBK(gp, 0); // will be reset to VBK(gp,nages) if using age-specific K - // SS_Label_Info_16.2.2 #Set up age specific k - if (Grow_type == 3) // age specific k + case 3: // age specific k ascending sequence { j = 1; for (a = 1; a <= nages; a++) @@ -289,12 +297,15 @@ FUNCTION void get_growth2(const int y) VBK(gp, a) = VBK(gp, a - 1); } } - VBK_temp = VBK(gp, nages); + VBK_work = VBK(gp, nages); + if (AFIX2 != 999) L_inf_plus = 1; + break; } - else if (Grow_type == 4) // age specific k reverse order, so age_k_points need to be descending + + case 4: // age specific k reverse order, so age_k_points need to be descending { j = 1; - VBK(gp, nages) = VBK_temp; + VBK(gp, nages) = VBK_work; // oldest age has the base parameter value for (a = nages - 1; a >= 0; a--) { if (a == Age_K_points(j)) @@ -308,29 +319,32 @@ FUNCTION void get_growth2(const int y) VBK(gp, a) = VBK(gp, a + 1); } } + if (AFIX2 != 999) L_inf_plus = 1; + break; } - else if (Grow_type == 5) // age specific k replacement, so age_k_points need to be descending + + case 5: // age specific k replacement, age_k_points need to be descending { + // all ages start at base parameter value, which is VBK_work j = 1; for (a = nages; a >= 0; a--) { if (a == Age_K_points(j)) { - VBK(gp, a) = mgp_adj(Ip + 2 + j) * VBK_temp; + VBK(gp, a) = mgp_adj(Ip + 2 + j) * VBK_work; if (j < Age_K_count) j++; } else { - VBK(gp, a) = VBK_temp; + VBK(gp, a) = VBK_work; } } - VBK_temp = VBK(gp, nages); + if (AFIX2 != 999) L_inf_plus = 1; + break; } - // get Linf from Lmax - // get Richards or growth cessation parameter if appropriate - if (Grow_type == 2) // Richards + case 2: // Richards or Gompertz growth { if (MGparm_def > 1 && gp > 1) { @@ -344,21 +358,19 @@ FUNCTION void get_growth2(const int y) inv_Richards = 1.0 / Richards(gp); if (AFIX2 == 999) { - L_inf(gp) = Lmax_temp(gp); + // L_inf(gp) = Lmax_temp(gp); // already set LinfR = pow(L_inf(gp), Richards(gp)); } else { LmaxR = pow(Lmax_temp(gp), Richards(gp)); - LinfR = LminR + (LmaxR - LminR) / (1. - mfexp(VBK_temp * VBK_seas(0) * (AFIX_delta))); + LinfR = LminR + (LmaxR - LminR) / (1. - mfexp(VBK_work * VBK_seas(0) * (AFIX_delta))); L_inf(gp) = pow(LinfR, inv_Richards); } - #ifdef DO_ONCE - if (do_once == 1) - echoinput << " linf " << L_inf(gp) << " VBK: " << VBK_temp << endl; - #endif + break; } - else if (Grow_type == 8) + + case 8: { if (MGparm_def > 1 && gp > 1) { @@ -369,51 +381,68 @@ FUNCTION void get_growth2(const int y) Richards(gp) = mgp_adj(Ip + 3); } L_inf(gp) = Lmax_temp(gp); - VBK_temp = -VBK(gp, 0) * VBK_seas(0); // t50 is the calculated inflection age for the decline in K - t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / VBK_temp) - 1.0) / Richards(gp); - } - else - { - if (AFIX2 == 999) - { - L_inf(gp) = Lmax_temp(gp); - } - else - { - L_inf(gp) = Lmin(gp) + (Lmax_temp(gp) - Lmin(gp)) / (1. - mfexp(VBK_temp * VBK_seas(0) * (AFIX_delta))); - #ifdef DO_ONCE - if (do_once == 1) - echoinput << VBK_temp << " " << VBK_seas(0) << " " << VBK_temp * VBK_seas(0) << " " << Lmax_temp(gp) << " " << L_inf(gp) << endl; - #endif - } + // note use of -VBK_work because VBK_work has negative sign + t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / (-VBK_work * VBK_seas(0))) - 1.0) / Richards(gp); + break; } + } - // SS_Label_Info_16.2.3 #Set up Lmin and Lmax in Start Year - if (y == styr) - { - Cohort_Lmin(gp) = Lmin(gp); // sets for all years and ages - } - else if (timevary_MG(y, 2) > 0) // using time-vary growth + if (L_inf_plus == 1) + { + L_inf(gp) = Lmin(gp) + (Lmax_temp(gp) - Lmin(gp)) / (1. - mfexp(VBK_work * VBK_seas(0) * (AFIX_delta))); + } + + // SS_Label_Info_16.2.3 #Set up Cohort_Lmin in Start Year + if (y == styr) + { + Cohort_Lmin(gp) = Lmin(gp); // sets for all years and ages + } + else if (timevary_MG(y, 2) > 0) // using time-vary growth + { + k = min(nages, (YrMax - y)); + for (a = 0; a <= k; a++) { - k = min(nages, (YrMax - y)); - for (a = 0; a <= k; a++) - { - Cohort_Lmin(gp, y + a, a) = Lmin(gp); - } // sets for future years so cohort remembers its size at birth; with Lmin(gp) being size at birth this year - } - } // end setup of parametric growth parameters - } // end switch between parametric and non-parametric growth + Cohort_Lmin(gp, y + a, a) = Lmin(gp); + } // sets for future years so cohort remembers its size at birth; with Lmin(gp) being size at birth this year + } + } // end parametric growth parameters + else + { + // grow_type == 7} // not implemented placeholder for empirical length entry + } #ifdef DO_ONCE if (do_once == 1) { - echoinput << "sex: " << gg << " GP: " << gp << " Lmin: " << Lmin(gp) << " Linf: " << L_inf(gp) << " VBK_temp: " << VBK_temp << " VBK@age: " << -VBK(gp) << endl; + echoinput << "sex: " << gg << " GP: " << gp << " Lmin: " << Lmin(gp) << " L_max: " << Lmax_temp(gp) << " Linf: " << L_inf(gp) << " VBK@age: " << -VBK(gp) << endl; if (Grow_type == 2) echoinput << " Richards: " << Richards(gp) << endl; if (Grow_type == 8) - echoinput << " Cessation_decay: " << Richards(gp) << endl; + echoinput << "t50 for growth cessation: " << t50 << " using decay_rate: " << Richards(gp) << endl; } #endif + if (y == styr) + { + // SS_Label_Info_16.2.4.1 #set up the delta in growth variability across ages if needed + // get more CV factors that may depend upon growth parameters, but cannot be time-varying + if (CV_const(gp) > 0) + { + if (CV_depvar_a == 0) + { + CV_delta(gp) = (CVLmax(gp) - CVLmin(gp)) / (Lmax_temp(gp) - Lmin(gp)); + } + else + { + CV_delta(gp) = (CVLmax(gp) - CVLmin(gp)) / (AFIX2_forCV - AFIX); + } + } + else + { + CV_delta(gp) = 0.0; + CV_G(gp) = CVLmin(gp); // sets all seasons and whole age range + } + } // end preliminary calcs + // SS_Label_Info_16.2.4 #Loop settlement events because growth starts at time of settlement g = g_Start(gp); // base platoon for (settle = 1; settle <= N_settle_timings; settle++) @@ -421,56 +450,28 @@ FUNCTION void get_growth2(const int y) g += N_platoon; // increment by N_platoon because only middle platoon has growth modeled if (use_morph(g) > 0) { - if (y == styr) - { - switch (Grow_type) - { - case 7: // non-parametric - { - break; - } - default: - { - // SS_Label_Info_16.2.4.1 #set up the delta in growth variability across ages if needed - if (CV_const(gp) > 0) - { - if (CV_depvar_a == 0) - { - CV_delta(gp) = (CVLmax(gp) - CVLmin(gp)) / (Lmax_temp(gp) - Lmin(gp)); - } - else - { - CV_delta(gp) = (CVLmax(gp) - CVLmin(gp)) / (AFIX2_forCV - AFIX); - } - } - else - { - CV_delta(gp) = 0.0; - CV_G(gp) = CVLmin(gp); // sets all seasons and whole age range - } - } - } - - // SS_Label_Info_16.2.4.1.1 #if y=styr, get size-at-age in first subseason of first season of this first year + // SS_Label_Info_16.2.4.1.1 #if y=styr, get size-at-age in first subseason of first season of this first year + // VBK_seas(0) is the sum of all the season durations * seasonal K adjustments. Normally it has a value of 1.0 + // but for seasons-as-years, it will be the value equal to the selected duration of those seasons +// Ave_Size(styr, 1, g)(0, first_grow_age(g)) = Lmin(gp); + if (y == styr) + { + switch (Grow_type) { case 1: { - VBK_temp2 = VBK_temp * VBK_seas(0); for (a = 0; a <= nages; a++) { - // Ave_Size(styr,1,g,a) = Lmin(gp) + (Lmin(gp)-L_inf(gp))* (mfexp(VBK_temp2*(real_age(g,1,a)-AFIX))-1.0); - Ave_Size(styr, 1, g, a) = L_inf(gp) + (Lmin(gp) - L_inf(gp)) * mfexp(VBK_temp2 * (real_age(g, 1, a) - AFIX)); + Ave_Size(styr, 1, g, a) = L_inf(gp) + (Lmin(gp) - L_inf(gp)) * mfexp(VBK_work * VBK_seas(0) * (real_age(g, 1, a) - AFIX)); } // done ageloop break; } case 2: // Richards { - Ave_Size(styr, 1, g)(0, first_grow_age(g)) = Lmin(gp); - VBK_temp2 = VBK_temp * VBK_seas(0); - for (a = first_grow_age(g); a <= nages; a++) + for (a = 0; a <= nages; a++) { - temp = LinfR + (LminR - LinfR) * mfexp(VBK_temp2 * (real_age(g, 1, a) - AFIX)); + temp = LinfR + (LminR - LinfR) * mfexp(VBK_work * VBK_seas(0) * (real_age(g, 1, a) - AFIX)); Ave_Size(styr, 1, g, a) = pow(temp, inv_Richards); } // done ageloop break; @@ -484,19 +485,18 @@ FUNCTION void get_growth2(const int y) case 3: // age-specific K, so need age-by-age calculations { ALK_idx = 1; - //VBK_seas(0) accounts for season duration + // NOTE: VBK_seas(0) accounts for season duration for (a = 0; a <= nages; a++) { - k2 = a - 1; + k2 = max(0, a - 1); if (lin_grow(g, ALK_idx, a) >= -1.0) // linear segment, or first time point beyond AFIX; { - Ave_Size(styr, 1, g, a) = Lmin(gp) + (Lmin(gp) - L_inf(gp)) * (mfexp(VBK(gp, 0) * VBK_seas(0) * (real_age(g, 1, a) - AFIX)) - 1.0); + Ave_Size(styr, 1, g, a) = Lmin(gp) + (Lmin(gp) - L_inf(gp)) * (mfexp(VBK(gp, k2) * VBK_seas(0) * (real_age(g, 1, a) - AFIX)) - 1.0); } else { Ave_Size(styr, 1, g, a) = Ave_Size(styr, 1, g, k2) + (mfexp(VBK(gp, k2) * VBK_seas(0)) - 1.0) * (Ave_Size(styr, 1, g, k2) - L_inf(gp)); } - // echoinput< -997.) // decay rate has been read; uses same code for Richards and standard + if (Linf_decay > 0.0001) // decay rate has been read. Should be an estimate of initial Z { temp1 = 1.0; temp4 = 1.0; temp = current_size; temp2 = mfexp(-Linf_decay); // cannot use natM or Z because growth is calculated first - #ifdef DO_ONCE - if (do_once == 1) - echoinput << " L_inf " << L_inf(gp) << " size@exactly maxage " << current_size << endl; - #endif - if (Grow_type < 3) - { - VBK_temp2 = VBK(gp, 0); - } - else - { - VBK_temp2 = VBK(gp, nages); - } - VBK_temp2 = (1.0 - mfexp(VBK_temp2 * VBK_seas(0))); + + dvariable VBK_temp = (1.0 - mfexp(VBK_work * VBK_seas(0))); for (a = nages + 1; a <= 3 * nages; a++) { temp4 *= temp2; // decay numbers at age by exp(-0.xxx) - current_size += (L_inf(gp) - current_size) * VBK_temp2; + current_size += (L_inf(gp) - current_size) * VBK_temp; temp += temp4 * current_size; temp1 += temp4; // accumulate numbers to create denominator for mean size calculation } @@ -567,15 +556,15 @@ FUNCTION void get_growth2(const int y) } Ave_Size(styr, 1, g, nages) = temp / temp1; // this is weighted mean size at nages } - else + else if (Linf_decay == -998.) // decay code set to -998 which omits growth in the plus group { // no adjustment } #ifdef DO_ONCE if (do_once == 1) - echoinput << " adjusted size at maxage " << Ave_Size(styr, 1, g, nages) << " using decay of: " << Linf_decay << endl; + echoinput << " adjusted size in plusgroup: " << Ave_Size(styr, 1, g, nages) << " using decay of: " << Linf_decay << endl; #endif - } // end initial year calcs + } // end initial year calcs // SS_Label_Info_16.2.4.2 #loop seasons for growth calculation for (s = 1; s <= nseas; s++) @@ -595,12 +584,13 @@ FUNCTION void get_growth2(const int y) else add_age = 0; // advance age or not // growth to next season - VBK_temp2 = (mfexp(VBK_temp * seasdur(s) * VBK_seas(s)) - 1.0); - // warning< styr && s == nseas) // do plus group mean size update { - if (y > styr && Linf_decay != -998.) + if (Linf_decay != -998.) { - + if (do_once == 1 && timevary_MG_firstyr == y) + { + warnstream << "plus group mean size is updated in years with time-vary growth beginning in: " << y << "; can turn this off with Linf_decay = -998"; + write_message (NOTE, 0); + } // 3.24 code - #ifdef DO_ONCE if (do_once == 1) - echoinput << " plus group calc: " - << " N _entering: " << natage(t, 1, g, nages - 1) << " N_inplus: " << natage(t, 1, g, nages) << " size in: " << Ave_Size(t + 1, 1, g, nages) << " old size: " << plusgroupsize << " "; - #endif + echoinput << niter << " " << y << " plus group calc: " + << " N _entering: " << natage(t, 1, g, nages - 1) << " N_inplus: " << natage(t, 1, g, nages) << " size in: " << Ave_Size(t + 1, 1, g, nages) << " old size: " << plusgroupsize << " "; temp = ((natage(t, 1, g, nages - 1) + 0.01) * Ave_Size(t + 1, 1, g, nages) + (natage(t, 1, g, nages) + 0.01) * plusgroupsize) / (natage(t, 1, g, nages - 1) + natage(t, 1, g, nages) + 0.02); Ave_Size(t + 1, 1, g, nages) = temp; - #ifdef DO_ONCE + + #ifdef DO_ONCE if (do_once == 1 && g == 1) echoinput << " final_val " << Ave_Size(t + 1, 1, g, nages) << endl; #endif @@ -845,7 +820,6 @@ FUNCTION void get_growth3(const int y, const int t, const int s, const int subse dvariable LminR; dvariable inv_Richards; dvariable t50; - dvariable VBK_temp2; ALK_idx = (s - 1) * N_subseas + subseas; // note that this changes a global value for (g = g_Start(1) + N_platoon; g <= gmorph; g += N_platoon) // looping the middle platoons for each sex*gp @@ -919,10 +893,9 @@ FUNCTION void get_growth3(const int y, const int t, const int s, const int subse } // done Richards case 8: // Cessation { - // VBK_temp2=-VBK_temp*seasdur(s); // negative to restore positive // t50 is the calculated inflection age for the decline in K - dvariable VBK_temp = -VBK(gp, 0); - t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / (VBK_temp)) - 1.0) / Richards(gp); + dvariable VBK_temp = -VBK(gp, 0) * VBK_seas(0); + t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / VBK_temp) - 1.0) / Richards(gp); for (a = 0; a <= nages; a++) { // calculate a full year's growth increment, then multiple by seasdur(s) @@ -968,6 +941,10 @@ FUNCTION void get_growth3(const int y, const int t, const int s, const int subse break; } } // done switch + #ifdef DO_ONCE + if (do_once == 1) + echoinput << " seas.subseas: " << s << "." << subseas <<" g:" << g <<" size: " << Ave_Size(t, subseas, g)(0, min(6, nages)) << " plusgroup: " << Ave_Size(t + 1, 1, g, nages) << endl; + #endif } // end need this platoon } // done platoon } // end calc size-at-age at a particular subseason diff --git a/SS_prelim.tpl b/SS_prelim.tpl index f02b9636..d270d012 100644 --- a/SS_prelim.tpl +++ b/SS_prelim.tpl @@ -1319,6 +1319,16 @@ write_message (WARN, 0); } } + // check for vulnerability in plusgroup size + warning << Ave_Size(styr, 1, g, nages - 1) << " " << 0.99 * Ave_Size(styr, 1, g, nages) << endl; + if (Ave_Size(styr, 1, g, nages - 1) < 0.99 * Ave_Size(styr, 1, g, nages)) + { + if ( Linf_decay == -998 && MG_active_firstyr(2) > 0 ) + { + warnstream << "There is greater than 1% growth in plusgroup, timevary growth starts in : " << MG_active_firstyr(2) << " and Linf_decay = -998 disables updating plusgroup size; suggest setting a Linf_decay value like initial Z"; + write_message (WARN, 0); + } + } } } } diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl index 77df33bd..e5ad3679 100644 --- a/SS_readcontrol_330.tpl +++ b/SS_readcontrol_330.tpl @@ -1689,6 +1689,7 @@ // stores years to calc non-constant MG parms (1=natmort; 2=growth; 3=wtlen & fec; 4=recr_dist&femfrac; 5=movement; 6=ageerrorkey; 7=catchmult) ivector timevary_pass(styr-3,YrMax+1) // extracted column ivector MG_active(0,7) // 0=all, 1=M, 2=growth 3=wtlen, 4=recr_dist&femfrac, 5=migration, 6=ageerror, 7=catchmult + ivector MG_active_firstyr(0,7) // first year in which a MGparm is timevarying vector env_data_pass(1,2) // holds min-max year with env data int do_densitydependent; @@ -1724,6 +1725,7 @@ echoinput << "Now read env, block/trend, and dev adjustments to MGparms " << endl; timevary_MG.initialize(); // stores years to calc non-constant MG parms (1=natmort; 2=growth; 3=wtlen & fec; 4=recr_dist; 5=movement) MG_active.initialize(); + MG_active_firstyr.initialize(); CGD_onoff = 0; timevary_parm_start_MG = 0; @@ -1829,6 +1831,7 @@ if (timevary_MG(y, f) > 0) { MG_active(f) = 1; + if (MG_active_firstyr(2) == 0 ) MG_active_firstyr(f) = y; timevary_MG(y, 0) = 1; // tracks active status for all MG types if(timevary_MG_firstyr == YrMax) timevary_MG_firstyr = y; // save for reporting in MSY and spawn_recruit output }