Skip to content

Alt lorenzen #329

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rick-Methot-NOAA merged 7 commits into from
Jun 8, 2022
Merged

Alt lorenzen #329

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
a036ed6
Alternate Lorenzen M
nathanvaughan-NOAA Jun 6, 2022
7ad19b6
Merge remote-tracking branch 'origin/main' into Alt_Lorenzen
nathanvaughan-NOAA Jun 6, 2022
34c2700
natM update
nathanvaughan-NOAA Jun 6, 2022
bc7b2d0
warning fix
nathanvaughan-NOAA Jun 7, 2022
500f874
Survivorship Lorenzen case 6
nathanvaughan-NOAA Jun 7, 2022
c1b4714
Lorenzen M average M change
nathanvaughan-NOAA Jun 7, 2022
bee54eb
style Lorenzen M changes
k-doering-NOAA Jun 8, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 56 additions & 2 deletions SS_biofxn.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1128,13 +1128,14 @@ FUNCTION void get_natmort()
{
int Loren_t = styr + (yz - styr) * nseas + s - 1;
natM(t_base + s, 0, gpi)(0, nages) = log(
elem_div(Ave_Size(Loren_t, mid_subseas, g)(0, nages), (Ave_Size(Loren_t, mid_subseas, g)(0, nages) + Loren_temp2))) *
Loren_M1;
elem_div(Ave_Size(Loren_t, mid_subseas, g)(0, nages), (Ave_Size(Loren_t, mid_subseas, g)(0, nages) + Loren_temp2))) *
Loren_M1;
if (s < Bseas(g))
{natM(t_base + s, 0, gpi, 0) = natM(t_base + s + 1, 0, gpi, 0);}
}
break;
}

// SS_Label_Info_17.1.2.3 #case 3: set to empirical M as read from file, no seasonal interpolation
case (3): // read age_natmort as constant
{
Expand Down Expand Up @@ -1220,6 +1221,59 @@ FUNCTION void get_natmort()
}
break;
}
// SS_Label_Info_17.1.2.6 #case 6: Calculate lorenzen M from survivorship over fixed age range
case 6: // Survivorship based Lorenzen M
{
Loren_temp2 = L_inf(gp) * (mfexp(-VBK(gp, K_index) * VBK_seas(0)) - 1.); // need to verify use of VBK_seas here
Loren_M1 = (natMparms(1, gp)); //This is the user specified average M over the input range of ages.
for (s = nseas ; s >= 1; s--)
{
int Loren_t = styr + (yz - styr) * nseas + s - 1;
dvariable loren_scale_extra = 0; //start with no extra scaler. This will be used if the maximum reference age is greater than nages.
int ref_age = int(natM_amax); //start with reference age equal to the input maximum age. This will be adjusted below to equal nages if the maximum age is greater than nages.
if (ref_age > nages)//if reference age is greater than accumulator age need math to approximate the unknown size/age bins
{
int extra_years = ref_age - nages;//determine how many extra ages will be included between accumulator age and reference age

//The following code is a simple difference approach to approximate the first and second rate of change in relative M to estimate approximate M for ages older than nages
//calculate proportional change in lorenzen M between second to last and last age group
dvariable d1 = 1 + (log((Ave_Size(Loren_t, mid_subseas, g)(nages)) / (Ave_Size(Loren_t, mid_subseas, g)(nages) + Loren_temp2)) -
log((Ave_Size(Loren_t, mid_subseas, g)(nages - 1)) / (Ave_Size(Loren_t, mid_subseas, g)(nages-1) + Loren_temp2))) /
log((Ave_Size(Loren_t, mid_subseas, g)(nages)) / (Ave_Size(Loren_t, mid_subseas, g)(nages) + Loren_temp2));

//calculate proportional change in lorenzen M between third to last and second to last age group
dvariable d2 = 1 + (log((Ave_Size(Loren_t, mid_subseas, g)(nages - 1))/(Ave_Size(Loren_t, mid_subseas, g)(nages - 1) + Loren_temp2)) -
log((Ave_Size(Loren_t, mid_subseas, g)(nages - 2)) / (Ave_Size(Loren_t, mid_subseas, g)(nages - 2) + Loren_temp2))) /
log((Ave_Size(Loren_t, mid_subseas, g)(nages - 1)) / (Ave_Size(Loren_t, mid_subseas, g)(nages - 1) + Loren_temp2));

//calculate the second order proportional change in proportional changes during the last two age pairs
dvariable d3 = 1 + (d1 - d2) / d1;

//project total proportion of last years M that will occur in all ages older than nages
for (int ey = 1; ey <= extra_years; ey++)
{
d1 = d1 * d3;//each year adjust the first order proportion by the second order proportion
loren_scale_extra += d1;//add that proportion to a scaler that will be multiplied by the nages M value
}
ref_age = nages; //set reference age to nages to use all available Ave_Size values
}

//Calculate loren_temp multiplier that achieves target average M
Loren_temp = (Loren_M1 * (natM_amax - natM_amin + 1)) / (sum(log(
elem_div(Ave_Size(Loren_t, mid_subseas, g)(natM_amin, ref_age), (Ave_Size(Loren_t, mid_subseas, g)(natM_amin, ref_age) + Loren_temp2))
)) + loren_scale_extra * log((Ave_Size(Loren_t, mid_subseas, g)(ref_age)) / (Ave_Size(Loren_t, mid_subseas, g)(ref_age) + Loren_temp2)));

natM(t_base + s, 0, gpi)(0, nages) = log(
elem_div(Ave_Size(Loren_t, mid_subseas, g)(0, nages),
(Ave_Size(Loren_t, mid_subseas, g)(0, nages) + Loren_temp2)))
* Loren_temp;
if (s < Bseas(g))
{
natM(t_base + s, 0, gpi, 0) = natM(t_base + s + 1, 0, gpi, 0);
}
}
break;
}
} // end natM_type switch

// SS_Label_Info_17.2 #calc an ave_age for the first gp as a scaling factor in logL for initial recruitment (R1) deviation
Expand Down
24 changes: 22 additions & 2 deletions SS_readcontrol_330.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -628,7 +628,7 @@
if(nseas>1) N_predparms+=N_pred*nseas;
natM_5_opt=0;
MGparm_point.initialize();
// 0=1Parm; 1=segmented; 2=Lorenzen; 3=agespecific; 4=agespec with seas interpolate; 5=Maunder_M
// 0=1Parm; 1=segmented; 2=Lorenzen; 3=agespecific; 4=agespec with seas interpolate; 5=Maunder_M; 6=Lorenzen range
*(ad_comm::global_datafile) >> natM_type;
echoinput<<natM_type<<" natM_type"<<endl;
switch(natM_type)
Expand Down Expand Up @@ -688,6 +688,15 @@
// Maunder_beta = natMparms(6,gp);
break;
}
case 6:
{
N_natMparms = 1;
*(ad_comm::global_datafile) >> natM_amin;
echoinput << natM_amin << " natM_minage for Lorenzen" << endl;
*(ad_comm::global_datafile) >> natM_amax;
echoinput << natM_amax << " natM_maxage for Lorenzen" << endl;
break;
}
}
END_CALCS

Expand Down Expand Up @@ -1047,6 +1056,14 @@
}
break;
}
case 6:
{
ParCount++;
ParmLabel += "NatM_Lorenzen_average" + GenderLbl(gg) + GP_Lbl(gp);
Parm_info += "val";
Parm_minmax.push_back(3);
break;
}
default:
{
break;
Expand Down Expand Up @@ -1580,7 +1597,10 @@
if(MGparm_PH(f)>0)
{MG_active(mgp_type(f))=1;}
}
if(natM_type==2 && MG_active(2)>0) MG_active(1)=1; // lorenzen M depends on growth
if ((natM_type == 2 || natM_type == 6) && MG_active(2) > 0)
{
MG_active(1) = 1; // lorenzen M depends on growth
}

j=N_MGparm;
if(timevary_parm_cnt_MG>0)
Expand Down
7 changes: 6 additions & 1 deletion SS_write_ssnew.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1708,7 +1708,7 @@ FUNCTION void write_nucontrol()
report4 << "#" << endl
<< "# setup for M, growth, wt-len, maturity, fecundity, (hermaphro), recr_distr, cohort_grow, (movement), (age error), (catch_mult), sex ratio " << endl;
report4 << "#_NATMORT" << endl
<< natM_type << " #_natM_type:_0=1Parm; 1=N_breakpoints;_2=Lorenzen;_3=agespecific;_4=agespec_withseasinterpolate;_5=BETA:_Maunder_link_to_maturity" << endl;
<< natM_type << " #_natM_type:_0=1Parm; 1=N_breakpoints;_2=Lorenzen;_3=agespecific;_4=agespec_withseasinterpolate;_5=BETA:_Maunder_link_to_maturity;_6=Lorenzen_range" << endl;
if (natM_type == 0)
{
report4 << " #_no additional input for selected M option; read 1P per morph" << endl;
Expand All @@ -1722,6 +1722,11 @@ FUNCTION void write_nucontrol()
{
report4 << natM_amin << " #_reference age for Lorenzen M; read 1P per morph" << endl;
}
else if (natM_type == 6)
{
report4 << natM_amin << " #_minimum age for Lorenzen" << endl
<< natM_amax << " #_maximum age for Lorenzen; read 1P per morph" << endl;
}
else if (natM_type >= 3 && natM_type < 5)
{
report4 << " #_Age_natmort_by sex x growthpattern (nest GP in sex)" << endl
Expand Down