Skip to content

Commit

Permalink
Merge pull request #204 from nmfs-stock-synthesis/predator_fleet
Browse files Browse the repository at this point in the history 
first commit for predator fleets
  • Loading branch information
@Rick-Methot-NOAA
Rick-Methot-NOAA authored Sep 21, 2021
2 parents cfbe093 + 48a5543 commit 568b406
Show file tree
Hide file tree
Showing 13 changed files with 292 additions and 59 deletions.
Binary file modified Change_log_for_SS_3.30.xlsx
View file
Binary file not shown.
22 changes: 22 additions & 0 deletions SS_benchfore.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1651,6 +1651,28 @@ FUNCTION void Get_Forecast()
{
Make_FishSelex(); // calcs fishery selex by current season, all fleets, current gmorph
}
if(N_pred>0)
{
for(f1=1;f1<=N_pred;f1++)
{
f=predator(f1);
pred_M2(f1,t)=mgp_adj(predparm_pointer(f1)); // base with no seasonal effect
if(nseas>1) pred_M2(f1,t)*=mgp_adj(predparm_pointer(f1)+s);

for (gp=1;gp<=N_GP*gender;gp++)
{g=g_Start(gp); // base platoon
for (settle=1;settle<=N_settle_timings;settle++)
{
g+=N_platoon;
int gpi=GP3(g); // GP*gender*settlement
natM(s,gpi)+=pred_M2(f1,t)*sel_al_3(s,g,f);
surv1(s,gpi)=mfexp(-natM(s,gpi)*seasdur_half(s));
surv2(s,gpi)=square(surv1(s,gpi));
}
}
}
}

} // end of seasonal biology

if(s==nseas) {adv_age=1;} else {adv_age=0;} // advance age or not when doing survivorship
Expand Down
54 changes: 37 additions & 17 deletions SS_expval.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,10 @@ FUNCTION void Get_expected_values(const int y,const int t);
}
case 2: // F rate
{
vbio=Hrate(f,t);
if(fleet_type(f)<4)
{vbio=Hrate(f,t);}
else
{vbio=pred_M2(predator_rev(f),t);}
break;
}
case 30: // spawning biomass #30
Expand Down Expand Up @@ -448,31 +451,48 @@ FUNCTION void Get_expected_values(const int y,const int t);
case(2): // DISCARD_OUTPUT
/* SS_Label_46.2 expected discard amount */
{
j=have_data(ALK_time,f,data_type,0); // number of observations
if(j>0)
if(have_data(ALK_time,f,data_type,0)>0) // number of observations
{
j=have_data(ALK_time,f,data_type,1); // only getting first observation for now
if(catch_ret_obs(f,t)>0.0 || y>endyr)
if(fleet_type(f)<=2)
{
if(disc_units(f)==3) // numbers regardless of catchunits for retained catch
if(catch_ret_obs(f,t)>0.0 || y>endyr)
{
if(disc_units(f)==3) // numbers regardless of catchunits for retained catch
{
exp_disc(f,j)=catch_fleet(t,f,4)-catch_fleet(t,f,6);
}
else if(catchunits(f)==1) // biomass units for retained and discarded catch
{
exp_disc(f,j)=catch_fleet(t,f,1)-catch_fleet(t,f,3); // discard in biomass
if(disc_units(f)==2) exp_disc(f,j) /= (catch_fleet(t,f,1) + 0.0000001);
}
else // numbers for retained and discarded catch
{
exp_disc(f,j)=catch_fleet(t,f,4)-catch_fleet(t,f,6); // discard in numbers
if(disc_units(f)==2) exp_disc(f,j) /= (catch_fleet(t,f,4) + 0.0000001);
}
if(exp_disc(f,j)<0.0) {
N_warn++;
warning<<N_warn<<" "<<f<<" "<<j<<" "<<exp_disc(f,j)<<" catches "<<catch_fleet(t,f)<<endl;
}
}
else
{
exp_disc(f,j)=catch_fleet(t,f,4)-catch_fleet(t,f,6);
exp_disc(f,j)=-1.;
}
else if(catchunits(f)==1) // biomass units for retained and discarded catch
}
else if (fleet_type(f)==4) // predator consumption stored in discard
{
if(disc_units(f)==3) // numbers
{
exp_disc(f,j)=catch_fleet(t,f,1)-catch_fleet(t,f,3); // discard in biomass
if(disc_units(f)==2) exp_disc(f,j) /= (catch_fleet(t,f,1) + 0.0000001);
exp_disc(f,j)=catch_fleet(t,f,4); // consumption in numbers
}
else // numbers for retained and discarded catch
else
{
exp_disc(f,j)=catch_fleet(t,f,4)-catch_fleet(t,f,6); // discard in numbers
if(disc_units(f)==2) exp_disc(f,j) /= (catch_fleet(t,f,4) + 0.0000001);
exp_disc(f,j)=catch_fleet(t,f,1); // consumption in biomass
}
if(exp_disc(f,j)<0.0) {
N_warn++;
warning<<N_warn<<" "<<f<<" "<<j<<" "<<exp_disc(f,j)<<" catches "<<catch_fleet(t,f)<<endl;
}
}
}
else
{
exp_disc(f,j)=-1.;
Expand Down
1 change: 1 addition & 0 deletions SS_global.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,7 @@ GLOBALS_SECTION
adstring F_report_label;
adstring SPR_report_label;
adstring onenum(4);
adstring onenum2(4);
adstring anystring;
adstring anystring2;
adstring_array version_info;
Expand Down
2 changes: 2 additions & 0 deletions SS_param.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,8 @@ PARAMETER_SECTION
// 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)
Expand Down
126 changes: 102 additions & 24 deletions SS_popdyn.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,20 +131,14 @@ FUNCTION void get_initial_conditions()
get_mat_fec(); // does just spawn season and subseason using ALK calculated just above
if(Hermaphro_Option!=0) get_Hermaphro();

if(MG_active(1)>0) get_natmort();
if(MG_active(1)>0)
{
get_natmort();
natM_M1=natM;
}
#ifdef DO_ONCE
if(do_once==1) cout<<" natmort OK"<<endl;
#endif
if(y>=Bmark_Yr(1)&&y<=Bmark_Yr(2))
{
for (s=1;s<=nseas;s++)
for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++)
{
natM_unf(s,gp)+=natM(s,gp); // need nseas to capture differences due to settlement
surv1_unf(s,gp)+=surv1(s,gp); // need nseas to capture differences due to settlement
surv2_unf(s,gp)+=surv2(s,gp); // need nseas to capture differences due to settlement
}
}

if(MG_active(4)>0) get_recr_distribution();
if(y>=Bmark_Yr(7)&&y<=Bmark_Yr(8))
Expand Down Expand Up @@ -240,6 +234,41 @@ FUNCTION void get_initial_conditions()
// SS_Label_Info_23.3.3 #for each platoon, combine size_at_age distribution with length selectivity and weight-at-length to get combined selectivity vectors
Make_FishSelex();
}

// SS_Label_Info_23.3.4 #add predator M2 to M1 to update seasonal M in styr
if(N_pred>0)
{
for(f1=1;f1<=N_pred;f1++)
{
f=predator(f1);
pred_M2(f1,t)=mgp_adj(predparm_pointer(f1)); // base with no seasonal effect
if(nseas>1) pred_M2(f1,t)*=mgp_adj(predparm_pointer(f1)+s);

for (gp=1;gp<=N_GP*gender;gp++)
{g=g_Start(gp); // base platoon
for (settle=1;settle<=N_settle_timings;settle++)
{
g+=N_platoon;
int gpi=GP3(g); // GP*gender*settlement
natM(s,gpi)+=pred_M2(f1,t)*sel_al_3(s,g,f);
surv1(s,gpi)=mfexp(-natM(s,gpi)*seasdur_half(s));
surv2(s,gpi)=square(surv1(s,gpi));
if(do_once==1) echoinput<<y<<" s "<<s<<" t "<<t<<" gp "<<gpi<<" M1: "<<natM_M1(s,gpi)<<endl;
if(do_once==1) echoinput<<y<<" s "<<s<<" t "<<t<<" gp "<<gpi<<" M1+M2: "<<natM(s,gpi)<<endl;
}
}
}
}

if(y>=Bmark_Yr(1)&&y<=Bmark_Yr(2))
{
for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++)
{
natM_unf(s,gp)+=natM(s,gp); // need nseas to capture differences due to settlement
surv1_unf(s,gp)+=surv1(s,gp); // need nseas to capture differences due to settlement
surv2_unf(s,gp)+=surv2(s,gp); // need nseas to capture differences due to settlement
}
}
}

#ifdef DO_ONCE
Expand Down Expand Up @@ -630,17 +659,12 @@ FUNCTION void get_time_series()
// but Maunder's M in get_natmort() may use the fecundity vector, so would be using the most recently calculated Problem??
if(Hermaphro_Option!=0) get_Hermaphro();
}
if(timevary_MG(y,1)>0) get_natmort();
if(y>=Bmark_Yr(1)&&y<=Bmark_Yr(2))
if(timevary_MG(y,1)>0)
{
for (s=1;s<=nseas;s++)
for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++)
{
natM_unf(s,gp)+=natM(s,gp);
surv1_unf(s,gp)+=surv1(s,gp);
surv2_unf(s,gp)+=surv2(s,gp);
}
get_natmort();
natM_M1=natM;
}

if(timevary_MG(y,4)>0) get_recr_distribution();
if(y>=Bmark_Yr(7)&&y<=Bmark_Yr(8))
{
Expand Down Expand Up @@ -752,8 +776,42 @@ FUNCTION void get_time_series()
for (g=1;g<=gmorph;g++)
if(use_morph(g)>0)
{Make_FishSelex();}

// SS_Label_Info_24.x.x #add predator M2 inside the yr,seas loop
if(N_pred>0)
{
for(f1=1;f1<=N_pred;f1++)
{
f=predator(f1);
pred_M2(f1,t)=mgp_adj(predparm_pointer(f1)); // base with no seasonal effect
if(nseas>1) pred_M2(f1,t)*=mgp_adj(predparm_pointer(f1)+s);
for (gp=1;gp<=N_GP*gender;gp++)
{g=g_Start(gp); // base platoon
for (settle=1;settle<=N_settle_timings;settle++)
{
g+=N_platoon;
int gpi=GP3(g); // GP*gender*settlement
natM(s,gpi)+=pred_M2(f1,t)*sel_al_3(s,g,f);
surv1(s,gpi)=mfexp(-natM(s,gpi)*seasdur_half(s));
surv2(s,gpi)=square(surv1(s,gpi));
if(do_once==1) echoinput<<y<<" s "<<s<<" t "<<t<<" gp "<<gpi<<" M1: "<<natM_M1(s,gpi)<<endl;
if(do_once==1) echoinput<<y<<" s "<<s<<" t "<<t<<" gp "<<gpi<<" M1+M2: "<<natM(s,gpi)<<endl;
}
}
}
}

if(y>=Bmark_Yr(1)&&y<=Bmark_Yr(2))
{
for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++)
{
natM_unf(s,gp)+=natM(s,gp);
surv1_unf(s,gp)+=surv1(s,gp);
surv2_unf(s,gp)+=surv2(s,gp);
}
}

}
// SS_Label_Info_24.2.2 #Compute spawning biomass if this is spawning season so recruits could occur later this season
// SPAWN-RECR: calc SPB in time series if spawning is at beginning of the season
if(s==spawn_seas && spawn_time_seas<0.0001) // compute spawning biomass if spawning at beginning of season so recruits could occur later this season
Expand Down Expand Up @@ -892,7 +950,7 @@ FUNCTION void get_time_series()
// SS_Label_Info_24.3.3 #Do fishing mortality
catage_tot.initialize();

if(catch_seas_area(t,p,0)==1 && fishery_on_off==1)
if( (catch_seas_area(t,p,0)==1 && fishery_on_off==1))
{
if(F_Method>1) // not Pope's
{
Expand Down Expand Up @@ -933,7 +991,7 @@ FUNCTION void get_time_series()
for (g=1;g<=gmorph;g++)
if(use_morph(g)>0)
{
Z_rate(t,p,g)=natM(s,GP3(g));
Z_rate(t,p,g)=natM(s,GP3(g)); // where natM already includes M2
for (f=1;f<=Nfleet;f++) //loop over fishing fleets to get Z
if (catch_seas_area(t,p,f)!=0)
{
Expand Down Expand Up @@ -1041,7 +1099,7 @@ FUNCTION void get_time_series()
}
} //close gmorph loop
} // close fishery
} // end continuous F method
} // end continuous F method

}
else
Expand Down Expand Up @@ -1129,8 +1187,28 @@ FUNCTION void get_time_series()
// SS_Label_Info_24.3.3.4 #No catch or fishery turned off, so set Z=M
for (g=1;g<=gmorph;g++)
if(use_morph(g)>0)
{Z_rate(t,p,g)=natM(s,GP3(g));}
{Z_rate(t,p,g)=natM(s,GP3(g));
Zrate2(p,g)=elem_div( (1.-mfexp(-seasdur(s)*Z_rate(t,p,g))), Z_rate(t,p,g));
}
}
if(N_predparms>0)
{
for(f1=1;f1<=N_pred;f1++)
{
f=predator(f1);
for (g=1;g<=gmorph;g++)
if(use_morph(g)>0)
{
catch_fleet(t,f,1)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),sel_al_1(s,g,f))*Zrate2(p,g);
catch_fleet(t,f,2)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),deadfish_B(s,g,f))*Zrate2(p,g);
catch_fleet(t,f,3)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),sel_al_2(s,g,f))*Zrate2(p,g); // retained bio
catch_fleet(t,f,4)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),sel_al_3(s,g,f))*Zrate2(p,g);
catch_fleet(t,f,5)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),deadfish(s,g,f))*Zrate2(p,g);
catch_fleet(t,f,6)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),sel_al_4(s,g,f))*Zrate2(p,g); // retained numbers
catage(t,f,g)=pred_M2(f1,t)*elem_prod(elem_prod(natage(t,p,g),deadfish(s,g,f)),Zrate2(p,g));
} //close gmorph loop
}
}
} //close area loop
if(s==1 && save_for_report==1){
Smry_Table(y,1)=totbio;
Expand Down
24 changes: 24 additions & 0 deletions SS_prelim.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -951,6 +951,30 @@ PRELIMINARY_CALCS_SECTION
// SS_Label_Info_6.8.4 #Call fxn get_natmort()
echoinput<<"ready to do natmort "<<endl;
get_natmort();
natM_M1=natM; // base M1 to which M2 is added

// SS_Label_Info_6.8.5 #add M2 for predator mortality
if(N_pred>0)
{
for(f1=1;f1<=N_pred;f1++)
{
pred_M2(f1,styr)=MGparm(predparm_pointer(f1)); // base with no seasonal effect
if(nseas==1)
{
for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++) {natM(s,gp)+=pred_M2(f1,styr);}
}
else
{
t=styr-1; // resets for each predator
for (s=1;s<=nseas;s++)
{
t++;
for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++) {natM(s,gp)+=pred_M2(f1,styr)*MGparm(predparm_pointer(f1)+s);}
}
}
}
}

natM = value(natM);
surv1 = value(surv1);
surv2 = value(surv2);
Expand Down
Loading

0 comments on commit 568b406

Please sign in to comment.