Skip to content

Add pred m to z not m #307

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 11 commits into from
May 27, 2022
Merged

Add pred m to z not m #307

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
82ec0b6
convert natM(s to natM(t as first step in tansitional SPR
Rick-Methot-NOAA Apr 27, 2022
af7314c
tune up the natM indexing
Rick-Methot-NOAA Apr 27, 2022
705eda1
add area dim to natM, need testing with predator model
Rick-Methot-NOAA May 3, 2022
f2e213b
fix for 2 seas, 2 area
Rick-Methot-NOAA May 4, 2022
9c19c1c
fix indexing mistake
Rick-Methot-NOAA May 4, 2022
963e8a1
works with vendance_predator, except forecast
Rick-Methot-NOAA May 5, 2022
501b6c9
forecast was correct; main has it wrong with predators
Rick-Methot-NOAA May 5, 2022
0059d34
simplify by removing array natM_M1
Rick-Methot-NOAA May 9, 2022
7719221
style: add spaces
k-doering-NOAA May 11, 2022
187b0f5
Merge branch 'main' into add_predM_to_Z_not_M
k-doering-NOAA May 11, 2022
6787baa
Merge branch 'main' into add_predM_to_Z_not_M
k-doering-NOAA May 27, 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
243 changes: 113 additions & 130 deletions SS_benchfore.tpl
View file
Open in desktop

Large diffs are not rendered by default.

75 changes: 34 additions & 41 deletions SS_biofxn.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -972,7 +972,7 @@ FUNCTION void get_growth3(const int y, const int t, const int s, const int subse

FUNCTION void get_natmort()
{
// SS_Label_Function #17 get_natmort
// SS_Label_Function #17 get_natmort for all seasons given this year's parameters
dvariable Loren_M1;
dvariable Loren_temp;
dvariable Loren_temp2;
Expand All @@ -989,7 +989,7 @@ FUNCTION void get_natmort()
int K_index;
K_index = VBK(1).indexmax();
Do_AveAge = 0;
t_base = styr + (yz - styr) * nseas - 1; // so looping s=1 to nseas; t=t_base+s
t_base = styr + (yz - styr) * nseas - 1; // so looping s=1 to nseas; t=t_base + s
Ip = -N_M_Grow_parms; // start counter for MGparms
// SS_Label_Info_17.1 #loop growth patterns in each gender
gp = 0;
Expand Down Expand Up @@ -1058,22 +1058,16 @@ FUNCTION void get_natmort()
{
for (s = 1; s <= nseas; s++)
{
if (docheckup == 1)
echoinput << "Natmort " << s << " " << gp << " " << gpi << " " << natMparms(1, gp);
natM(s, gpi) = natMparms(1, gp);
surv1(s, gpi) = mfexp(-natMparms(1, gp) * seasdur_half(s)); // refers directly to the constant value
surv2(s, gpi) = square(surv1(s, gpi));
if (docheckup == 1)
echoinput << " surv " << surv1(s, gpi) << endl;
natM(t_base + s, 0, gpi) = natMparms(1, gp);
}
break;
}

// SS_Label_Info_17.1.2.1 #case 1: N breakpoints
case 1: // breakpoints
{
dvariable natM1;
dvariable natM2;
dvariable natM_A;
Copy link
Contributor

@k-doering-NOAA k-doering-NOAA May 6, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why the change to A and B? it would be good to be consistent about non predation natural mortality and predation mortality variable names. For instance, if another variable name is needed, it would be better to name it with how it differs from natM1 (e.g., natM1_all_gps, natM1_user_inputs, etc.);

dvariable natM_B;
for (s = 1; s <= nseas; s++)
{
if (s >= Bseas(g))
Expand All @@ -1087,41 +1081,39 @@ FUNCTION void get_natmort()
t_age = 1.0 + azero_seas(s) - azero_G(g);
}
natM_amax = NatM_break(1);
natM2 = natMparms(1, gp);
natM_B = natMparms(1, gp);
k = a;

for (loop = 1; loop <= N_natMparms + 1; loop++)
{
natM_amin = natM_amax;
natM1 = natM2;
natM_A = natM_B;
if (loop <= N_natMparms)
{
natM_amax = NatM_break(loop);
natM2 = natMparms(loop, gp);
natM_B = natMparms(loop, gp);
}
else
{
natM_amax = r_ages(nages) + 1.;
}
if (natM_amax > natM_amin)
{
temp = (natM2 - natM1) / (natM_amax - natM_amin);
temp = (natM_B - natM_A) / (natM_amax - natM_amin);
} // calc the slope
else
{
temp = 0.0;
}
while (t_age < natM_amax && a <= nages)
{
natM(s, gpi, a) = natM1 + (t_age - natM_amin) * temp;
natM(t_base + s, 0, gpi, a) = natM_A + (t_age - natM_amin) * temp;
t_age += 1.0;
a++;
}
}
if (k == 1)
natM(s, gpi, 0) = natM(s, gpi, 1);
surv1(s, gpi) = mfexp(-natM(s, gpi) * seasdur_half(s));
surv2(s, gpi) = square(surv1(s, gpi));
natM(t_base + s, 0, gpi, 0) = natM(t_base + s, 0, gpi, 1);
} // end season
break;
} // end natM_type==1
Expand All @@ -1135,13 +1127,11 @@ FUNCTION void get_natmort()
for (s = nseas; s >= 1; s--)
{
int Loren_t = styr + (yz - styr) * nseas + s - 1;
natM(s, gpi)(0, nages) = log(
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;
if (s < Bseas(g))
natM(s, gpi, 0) = natM(s + 1, gpi, 0);
surv1(s, gpi) = value(mfexp(-natM(s, gpi) * seasdur_half(s)));
surv2(s, gpi) = value(square(surv1(s, gpi)));
{natM(t_base + s, 0, gpi, 0) = natM(t_base + s + 1, 0, gpi, 0);}
}
break;
}
Expand All @@ -1150,9 +1140,7 @@ FUNCTION void get_natmort()
{
for (s = 1; s <= nseas; s++)
{
natM(s, gpi) = Age_NatMort(gp);
surv1(s, gpi) = value(mfexp(-natM(s, gpi) * seasdur_half(s)));
surv2(s, gpi) = value(square(surv1(s, gpi)));
natM(t_base + s, 0, gpi) = Age_NatMort(gp);
}
break;
}
Expand All @@ -1168,7 +1156,7 @@ FUNCTION void get_natmort()
t_age = azero_seas(s) - azero_G(g);
for (a = k; a <= nages - 1; a++)
{
natM(s, gpi, a) = Age_NatMort(gp, a) + t_age * (Age_NatMort(gp, a + 1) - Age_NatMort(gp, a));
natM(t_base + s, gpi, a) = Age_NatMort(gp, a) + t_age * (Age_NatMort(gp, a + 1) - Age_NatMort(gp, a));
} // end age
}
else
Expand All @@ -1177,13 +1165,11 @@ FUNCTION void get_natmort()
t_age = azero_seas(s) + (1. - azero_G(g));
for (a = k; a <= nages - 1; a++)
{
natM(s, gpi, a) = Age_NatMort(gp, a) + t_age * (Age_NatMort(gp, a + 1) - Age_NatMort(gp, a));
natM(t_base + s, 0, gpi, a) = Age_NatMort(gp, a) + t_age * (Age_NatMort(gp, a + 1) - Age_NatMort(gp, a));
} // end age
natM(s, gpi, 0) = natM(s, gpi, 1);
natM(t_base + s, 0, gpi, 0) = natM(t_base + s, 0, gpi, 1);
}
natM(s, gpi, nages) = Age_NatMort(gp, nages);
surv1(s, gpi) = mfexp(-natM(s, gpi) * seasdur_half(s));
surv2(s, gpi) = square(surv1(s, gpi));
natM(t_base + s, 0, gpi, nages) = Age_NatMort(gp, nages);
} // end season
break;
}
Expand Down Expand Up @@ -1214,10 +1200,10 @@ FUNCTION void get_natmort()
XX_mature(First_Mature_Age, nages) = 1. / (1. + mfexp(Maunder_beta * (Ave_Size(t, mid_subseas, g)(First_Mature_Age, nages) - Maunder_L50)));
{
// original equation had:
// natM(s,gpi,a) = Maunder_Mjuv*pow(Ave_Size(t,ALK_idx,g,a)/Maunder_Lmat,Maunder_lambda) +
// natM(t_base + s,gpi,a) = Maunder_Mjuv*pow(Ave_Size(t,ALK_idx,g,a)/Maunder_Lmat,Maunder_lambda) +
// (Maunder_Mmat-Maunder_Mjuv*pow(Ave_Size(t,ALK_idx,g,a)/Maunder_Lmat,Maunder_lambda))*XXmaturity_Fem(a)XX;
natM(s, gpi) = Maunder_Mjuv * pow((Ave_Size(t, mid_subseas, g) / Maunder_Lmat), Maunder_lambda);
natM(s, gpi) += elem_prod((Maunder_Mmat - natM(s, gpi)), XX_mature);
natM(t_base + s, 0, gpi) = Maunder_Mjuv * pow((Ave_Size(t, mid_subseas, g) / Maunder_Lmat), Maunder_lambda);
natM(t_base + s, 0, gpi) += elem_prod((Maunder_Mmat - natM(t_base + s, 0, gpi)), XX_mature);
}
if (do_once == 1)
{
Expand All @@ -1229,7 +1215,7 @@ FUNCTION void get_natmort()
echoinput << "avesize/Lmat " << Ave_Size(t, mid_subseas, g) / Maunder_Lmat << endl;
echoinput << " natM_juv: " << Maunder_Mjuv * pow((Ave_Size(t, mid_subseas, g) / Maunder_Lmat), Maunder_lambda) << endl;
echoinput << " natM_mat: " << (Maunder_Mmat)*XX_mature << endl;
echoinput << " natM_combined: " << natM(s, gpi) << endl;
echoinput << " natM_combined: " << natM(t_base + s, 0, gpi) << endl;
}
}
break;
Expand All @@ -1240,20 +1226,27 @@ FUNCTION void get_natmort()
if (Do_AveAge == 0)
{
Do_AveAge = 1;
ave_age = 1.0 / natM(1, gpi, nages / 2) - 0.5;
ave_age = 1.0 / natM(t_base+1, 0, gpi, nages / 2) - 0.5;
}

#ifdef DO_ONCE
if (do_once == 1)
{
for (s = 1; s <= nseas; s++)
echoinput << "Natmort seas:" << s << " sex:" << gg << " Gpat:" << GPat << " sex*Gpat:" << gp << " settlement:" << settle << " gpi:" << gpi << endl
<< " M: " << natM(s, gpi) << endl;
<< " M: " << natM(t_base + s, 0, gpi) << endl;
}
#endif
} // end use of this morph
} // end settlement
} // end growth pattern x gender loop
natM_M1 = natM; // set M1 equal to M; M2 can be added later if predators are used
for (s = 1; s <= nseas; s++)
for (p = 1; p <= pop; p++)
{
natM(t_base + s, p) = natM(t_base + s, 0); // copy M1 to eack area's M;
// p=0 holds that M1 as the base M with no predators
// pred_M2 will be added later on area-specific basis
}
} // end nat mort

FUNCTION void get_recr_distribution()
Expand Down Expand Up @@ -1942,8 +1935,8 @@ FUNCTION void get_saveGparm()
save_G_parm(save_gparm, 12) = value(CVLmax(gp));
save_G_parm(save_gparm, 13) = natM_amin;
save_G_parm(save_gparm, 14) = natM_amax;
save_G_parm(save_gparm, 15) = value(natM(1, GP3(g), 0));
save_G_parm(save_gparm, 16) = value(natM(1, GP3(g), nages));
save_G_parm(save_gparm, 15) = value(natM(t_base+1, 0, GP3(g), 0));
save_G_parm(save_gparm, 16) = value(natM(t_base+1, 0, GP3(g), nages));
if (gg == 1)
{
for (k = 1; k <= 6; k++)
Expand Down
6 changes: 4 additions & 2 deletions SS_objfunc.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1320,17 +1320,19 @@ FUNCTION void Process_STDquant()
}

// NatM
// shortcut natM( doing area 1 only
if (NatM_Std_Cnt > 0)
{
t = styr + (endyr - styr) * nseas; // first season of selected year
for (i = 1; i <= NatM_Std_Cnt; i++)
{
j = NatM_Std_Pick(i); // selected age
k = g_finder(Do_NatM_Std, 1); // selected GP and gender gp3
Extra_Std(gender * (Selex_Std_Cnt + Growth_Std_Cnt + NatAge_Std_Cnt) + i) = natM(1, k, j);
Extra_Std(gender * (Selex_Std_Cnt + Growth_Std_Cnt + NatAge_Std_Cnt) + i) = natM(t, 1, k, j);
if (gender == 2)
{
k = g_finder(Do_NatM_Std, 2); // selected GP and gender gp3
Extra_Std(gender * (Selex_Std_Cnt + Growth_Std_Cnt + NatAge_Std_Cnt) + NatM_Std_Cnt + i) = natM(1, k, j);
Extra_Std(gender * (Selex_Std_Cnt + Growth_Std_Cnt + NatAge_Std_Cnt) + NatM_Std_Cnt + i) = natM(t, 1, k, j);
}
}
}
Expand Down
34 changes: 15 additions & 19 deletions SS_param.tpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ PARAMETER_SECTION
init_bounded_number dummy_parm(0,2,dummy_phase) // estimate in phase 0

!! // SS_Label_Info_5.1.1 #Create MGparm vector and associated arrays
// natural mortality and growth
// growth
init_bounded_number_vector MGparm(1,N_MGparm2,MGparm_LO,MGparm_HI,MGparm_PH)
vector femfrac(1,N_GP*gender);
vector L_inf(1,N_GP*gender);
Expand All @@ -52,20 +52,7 @@ PARAMETER_SECTION

vector Lmin(1,N_GP*gender);
vector Lmin_last(1,N_GP*gender);
// vector natM1(1,N_GP*gender)
// vector natM2(1,N_GP*gender)
matrix natMparms(1,N_natMparms,1,N_GP*gender)
3darray natM(1,nseas,1,N_GP*gender*N_settle_timings,0,nages) // need nseas to capture differences due to settlement
3darray natM_M1(1,nseas,1,N_GP*gender*N_settle_timings,0,nages) // need nseas to capture differences due to settlement
matrix pred_M2(1,N_pred,styr-3*nseas,TimeMax_Fcast_std+nseas); // index by t
3darray natM_unf(1,nseas,1,N_GP*gender*N_settle_timings,0,nages) // need nseas to capture differences due to settlement
3darray natM_endyr(1,nseas,1,N_GP*gender*N_settle_timings,0,nages) // need nseas to capture differences due to settlement
3darray surv1(1,nseas,1,N_GP*gender*N_settle_timings,0,nages)
3darray surv1_unf(1,nseas,1,N_GP*gender*N_settle_timings,0,nages)
3darray surv1_endyr(1,nseas,1,N_GP*gender*N_settle_timings,0,nages)
3darray surv2(1,nseas,1,N_GP*gender*N_settle_timings,0,nages)
3darray surv2_unf(1,nseas,1,N_GP*gender*N_settle_timings,0,nages)
3darray surv2_endyr(1,nseas,1,N_GP*gender*N_settle_timings,0,nages)

vector CVLmin(1,N_GP*gender)
vector CVLmax(1,N_GP*gender)
vector CV_const(1,N_GP*gender)
Expand Down Expand Up @@ -278,19 +265,30 @@ PARAMETER_SECTION
3darray natage_temp(1,pop,1,gmorph,0,nages)
number ave_age // average age of fish in unfished population; used to weight R1

!!// SS_Label_Info_5.1.3 #Create F parameters and associated arrays and constants
!!// SS_Label_Info_5.1.3 #Create M, F, and Z parameters and associated arrays and constants
init_bounded_number_vector init_F(1,N_init_F,init_F_LO,init_F_HI,init_F_PH)
matrix est_equ_catch(1,nseas,1,Nfleet)

// natural, predation and fishing mortality
matrix natMparms(1,N_natMparms,1,N_GP*gender) // will be derived from the MGparms
!!if(Do_Forecast>0) {k=TimeMax_Fcast_std+nseas;} else {k=TimeMax+nseas;}
4darray natM(styr-3*nseas,k,0,pop,1,N_GP*gender*N_settle_timings,0,nages) // M1 + pred_M2, see desc. in biofxn.tpl
// 3darray natM_M1(1,nseas,1,N_GP*gender*N_settle_timings,0,nages) // base M, biology only
matrix pred_M2(1,N_pred,styr-3*nseas,TimeMax_Fcast_std+nseas); // predator M2

// add area (pop) dimension to same dimension as season; use s1=(p-1)*pop + s
3darray surv1(1,nseas*pop,1,N_GP*gender*N_settle_timings,0,nages)
3darray surv2(1,nseas*pop,1,N_GP*gender*N_settle_timings,0,nages)
4darray Z_rate(styr-3*nseas,k,1,pop,1,gmorph,0,nages)
3darray Zrate2(1,pop,1,gmorph,0,nages)
matrix Hrate(1,Nfleet,styr-3*nseas,k) //Harvest Rate for each fleet; this is F
4darray natage(styr-3*nseas,k,1,pop,1,gmorph,0,nages) // add +1 year
4darray catage(styr-nseas,k,1,Nfleet,1,gmorph,0,nages)
4darray disc_age(styr-3*nseas,TimeMax_Fcast_std+nseas,1,2*N_retain_fleets,1,gmorph,0,nages);
4darray equ_catage(1,nseas,1,Nfleet,1,gmorph,0,nages)
4darray equ_numbers(1,nseas,1,pop,1,gmorph,0,3*nages)
4darray equ_Z(1,nseas,1,pop,1,gmorph,0,nages)
matrix catage_tot(1,gmorph,0,nages)//sum the catches for all fleets, reuse matrix each year
matrix Hrate(1,Nfleet,styr-3*nseas,k) //Harvest Rate for each fleet
matrix bycatch_F(1,Nfleet,1,nseas)
3darray catch_fleet(styr-3*nseas,k,1,Nfleet,1,6) // 1=sel_bio, 2=kill_bio; 3=ret_bio; 4=sel_num; 5=kill_num; 6=ret_num
matrix annual_catch(styr-1,YrMax,1,6) // same six as above
Expand Down Expand Up @@ -318,8 +316,6 @@ PARAMETER_SECTION
number harvest_rate; // Harvest rate
number maxpossF;

4darray Z_rate(styr-3*nseas,k,1,pop,1,gmorph,0,nages)
3darray Zrate2(1,pop,1,gmorph,0,nages)

LOCAL_CALCS
if(N_Fparm>0) // continuous F
Expand Down
Loading