Skip to content

Commit

Permalink
working version with a,b
Browse files Browse the repository at this point in the history
  • Loading branch information
@Rick-Methot-NOAA
Rick-Methot-NOAA committed May 22, 2024
1 parent 4de6922 commit e166f2e
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 26 deletions.
41 changes: 19 additions & 22 deletions SS_recruit.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,15 +49,19 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
case 3: // Beverton-Holt
{
steepness = SR_parm_work(2);
// dvariable SPR = SSB_virgin_adj / Recr_virgin;
// dvariable alpha = (4.0 * steepness) / ( SPR * (1. - steepness)) ;
// dvariable beta = (5.0 * steepness - 1.0 )/((1.0 - steepness ) * Recr_virgin * SPR );
// beta = (1.0 / Recr_virgin) * (alpha - (1.0 / SPR));
NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_curr_adj) /
(SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);
break;
}

case 10: // Beverton-Holt with alpha beta per WHAM: R = A*S/(1+B*S)
{
alpha = mfexp(SR_parm_work(3));
beta = mfexp(SR_parm_work(4));
NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
break;
}

// SS_Label_43.3.4 constant expected recruitment
case 4: // none
{
Expand Down Expand Up @@ -159,14 +163,6 @@ FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariab
break;
}

case 10: // Beverton-Holt with alpha beta
{
// WHAM based on R = A*S/(1+B*S)
alpha = mfexp(SR_parm_work(3));
beta = mfexp(SR_parm_work(4));
NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
break;
}
}
RETURN_ARRAYS_DECREMENT();
return NewRecruits;
Expand Down Expand Up @@ -318,18 +314,19 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
SPR_virgin = SSB_virgin / Recr_virgin;
alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
report5<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_virgin<<" SPR "<<SPR_virgin<<" steep "<<steepness<<endl;

// report5 <<" SSB_unf "<<SSB_virgin<<" SPR_unf "<<SPR_virgin<<" steep: "<<steepness<<" R0: "<<Recr_virgin << endl;
// report5 <<" derive_alpha "<<alpha<<" derive_beta "<<beta << endl;
// report5 << " deriv_h: " << alpha * SPR_virgin / (4. + alpha * SPR_virgin) << " derive_R0: " << 1. / beta * (alpha - (1. / SPR_virgin))<<endl;
B_equil = (alpha * SPR_temp - 1.0) / beta;
// report5<<" spr "<<SPR_temp<<" Beq "<<B_equil<<endl;
B_equil = posfun(B_equil, 0.0001, temp);
R_equil = alpha * B_equil / (1.0 + beta * B_equil);
// report5<<" Beq "<<B_equil<<" Req "<<R_equil<<endl;
// R_equil = (4. * steepness * Recr_virgin * B_equil) / (SSB_virgin * (1. - steepness) + (5. * steepness - 1.) * B_equil); //Beverton-Holt
// report5<<" SS3 Beq "<<B_equil<<" Req "<<R_equil<<" alpha "<<alpha<<" beta "<<beta<<" SSB_unf "<<SSB_unf<<endl;
// report5 << "SPR_input: " << SPR_temp << " B_equil: " << B_equil << " R_equil: "<<R_equil << endl<<endl;

break;
}

case 10: // Beverton-Holt with alpha and beta parameterization using R = A*S/(B+S) approach
case 10: // Beverton-Holt with alpha and beta parameterization using R = A*S/(1+B*S) approach; same as WHAM
{
alpha = mfexp(SRparm(3));
beta = mfexp(SRparm(4));
Expand All @@ -350,11 +347,11 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
// SS_Label_44.1.5 Hockey Stick
case 5: // hockey stick
{
alpha = SRparm(3) * Recr_virgin; // min recruitment level
dvariable hockey_min = SRparm(3) * Recr_virgin; // min recruitment level
// temp=SSB_virgin/R0*steepness; // spawners per recruit at inflection
beta = (Recr_virgin - alpha) / (steepness * SSB_virgin); // slope of recruitment on spawners below the inflection
B_equil = Join_Fxn(0.0 * SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin * steepness, SPR_temp, alpha / ((1. / SPR_temp) - beta), SPR_temp * Recr_virgin);
R_equil = Join_Fxn(0.0 * SSB_virgin, SSB_virgin, SSB_virgin * steepness, B_equil, alpha + beta * B_equil, Recr_virgin);
dvariable hockey_slope = (Recr_virgin - hockey_min) / (steepness * SSB_virgin); // slope of recruitment on spawners below the inflection
B_equil = Join_Fxn(0.0 * SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin, SSB_virgin / Recr_virgin * steepness, SPR_temp, hockey_min / ((1. / SPR_temp) - hockey_slope), SPR_temp * Recr_virgin);
R_equil = Join_Fxn(0.0 * SSB_virgin, SSB_virgin, SSB_virgin * steepness, B_equil, hockey_min + hockey_slope * B_equil, Recr_virgin);
break;
}
// SS_Label_44.1.7 3 parameter survival based
Expand Down
8 changes: 4 additions & 4 deletions SS_write_report.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1739,7 +1739,7 @@ FUNCTION void write_bigoutput()
else if (SR_fxn == 3)
{
alpha = 4.0 * steepness / (SPR_virgin * (1. - steepness));
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin);
beta = (5.0 * steepness - 1.0) / ((1. - steepness) * SSB_virgin);
SS2out << " Ln_alpha_derived: " << log(alpha) << " alpha " << alpha << endl;
SS2out << " Ln_beta_derived: " << log(beta) << " beta " << beta;
}
Expand Down Expand Up @@ -1791,7 +1791,7 @@ FUNCTION void write_bigoutput()
}
SS2out << endl;

SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SPR0 h R0" << endl;
SS2out << "Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SSBe/R0 h R0" << endl;
SS2out << "S/Rcurve " << SSB_virgin << " " << Recr_virgin << endl;
y = styr - 2;
SS2out << "Virg " << SSB_yr(y) << " " << exp_rec(y) << " - " << 0.0 << " Virg " << SSB_B_yr(y) << " " << SSB_N_yr(y) << " 0.0 " << endl;
Expand Down Expand Up @@ -1842,10 +1842,10 @@ FUNCTION void write_bigoutput()
{
SS2out << " _ _ Fixed";
}
temp = Smry_Table(y,11) / Recr_virgin;
dvariable SPR = Smry_Table(y, 11) / Recr_virgin;
alpha = mfexp(SR_parm_byyr(y,3));
beta = mfexp(SR_parm_byyr(y,4));
SS2out << " " << temp << " " << alpha * temp / (4. + alpha * temp) << " " << 1. / beta * (alpha - (1. / temp));
SS2out << " " << SPR << " " << alpha * SPR / (4. + alpha * SPR) << " " << 1. / beta * (alpha - (1. / SPR));
SS2out << SR_parm_byyr(y)(1,4) << endl;
}

Expand Down

0 comments on commit e166f2e

Please sign in to comment.