nmfs-ost
191 spawner recruit rebased
continuation of PR #442 after updating from 3.30.23 main on 11/06/2024 This implements SRR function 10 for Bev-Holt with alpha, beta also adds two new controls to give user more control of how time-varying biology affects benchmark and control rule implementation.
Resolves issue #191, #341, #625
See the checkboxes at top of issue #191 for the list of major changes with this branch.
@e-perl-NOAA, do you know why most of the github actions for this Pull Request were all skipped? I was able to compile on my computer for testing, it would be helpful to know which tests are passing.
It is because it's a draft PR, once it's converted from a draft PR they will all run.
@Rick-Methot-NOAA The test is now failing because the following values have changed in the test models that have time-varying biology:
| quantity | new value | ref_value | diff | perc_change | ratio | mod_name |
|---|---|---|---|---|---|---|
| maxgrad | 0.00E+00 | 6.22E-05 | -6.22E-05 | -100 | 0 | growth_timevary |
| SR_LN(R0)_se | 0.00E+00 | 1.44E-01 | -1.44E-01 | -100 | 0 | growth_timevary |
| SSB_unfished_se | 0.00E+00 | 3.29E+05 | -3.29E+05 | -100 | 0 | growth_timevary |
| ForeCatch_2012_se | 0.00E+00 | 8.81E+04 | -8.81E+04 | -100 | 0 | growth_timevary |
| Catch_like | 5.78E-13 | 5.78E-13 | 1.00E-18 | 1.73E-04 | 1 | platoons |
| Sum_recdevs_like | -9.51E-16 | 1.49E-16 | -1.10E-15 | -7.37E+02 | -6.3721 | platoons |
| maxgrad | 0.00E+00 | 8.30E-05 | -8.30E-05 | -1.00E+02 | 0 | platoons |
| SR_LN(R0)_se | 0.00E+00 | 9.39E-02 | -9.39E-02 | -1.00E+02 | 0 | platoons |
| SSB_unfished_se | 0.00E+00 | 5.71E+03 | -5.71E+03 | -1.00E+02 | 0 | platoons |
| ForeCatch_2004_val | NaN | 3.55E+03 | NaN | NaN | NaN | platoons |
| ForeCatch_2004_se | 0.00E+00 | 1.58E+03 | -1.58E+03 | -1.00E+02 | 0 | platoons |
| Sum_recdevs_like | -1.78E-15 | 2.66E-15 | -4.44E-15 | -166.667 | -0.66667 | three_area_nomove |
| maxgrad | 0.00E+00 | 2.97E-04 | -2.97E-04 | -100 | 0 | three_area_nomove |
| SR_LN(R0)_se | 0.00E+00 | 3.31E-02 | -3.31E-02 | -100 | 0 | three_area_nomove |
| SSB_unfished_se | 0.00E+00 | 5.19E+03 | -5.19E+03 | -100 | 0 | three_area_nomove |
| ForeCatch_2010_val | NaN | 1.06E+04 | NaN | NaN | NaN | three_area_nomove |
| ForeCatch_2010_se | 0.00E+00 | 1.28E+03 | -1.28E+03 | -100 | 0 | three_area_nomove |
Looks like some issue with the Hessian matrix not inverting to provide uncertainty for those models.
I think the problem is the NAN for forecatch, which then causes hessian failure.
I am ready to convert this to a real PR. OK?
@Rick-Methot-NOAA, that sounds good to me. @e-perl-NOAA and I will get to work on the r4ss side soon.
for the time-vary growth model, the blocks end in 2009, so the parameter values revert to 0 in 2010 for the forecast. Not sure why this is only a problem now. the above seems not the problem since selex parameters are using average approach in forecast. Trial and error with number of forecast loops finds that the problem occurs in the second forecast loop. Examination of forecast-report.sso shows problem to be in the calculation of the ABC "Ramp&buffer"
@Rick-Methot-NOAA, I'm working on sorting out the r4ss read of the revised output and I see see that the GLOBAL_MSY outputs are no longer filled in at the end of Report.sso (comparison in screenshot below). Those tables aren't actually read by r4ss::SS_output(), so it's not an issue for what I'm working on, but they seem like useful quantities and could be added to r4ss in the future if they get output.
@Rick-Methot-NOAA I updated the two time-varying models in the test models repo to use the new format. The no-est action is now "failing" (the results are different for one model). The changes I made to those models can be seen here. Please make sure that I updated them correctly.
| quantity | value | ref_value | diff | perc_change | ratio | mod_name |
|---|---|---|---|---|---|---|
| InitEQ_regime_like | 1.46E-30 | 0.00E+00 | 1.46E-30 | 0 | NA | tagging_mirrored_sel |
| Sum_recdevs_like | -2.61E-16 | 4.42E-17 | -3.05E-16 | -690.197 | -5.90197 | tagging_mirrored_sel |
| Forecast_Recruitment_like | 0.00E+00 | 0.00E+00 | 0.00E+00 | 0 | NA | tagging_mirrored_sel |
| Parm_devs_like | 0.00E+00 | 0.00E+00 | 0.00E+00 | 0 | NA | tagging_mirrored_sel |
| F_Ballpark_like | 0.00E+00 | 0.00E+00 | 0.00E+00 | 0 | NA | tagging_mirrored_sel |
| Crash_Pen_like | 0.00E+00 | 0.00E+00 | 0.00E+00 | 0 | NA | tagging_mirrored_sel |
| maxgrad | 0.00E+00 | 3.89E-03 | -3.89E-03 | -100 | 0 | tagging_mirrored_sel |
| SR_LN(R0)_se | 0.00E+00 | 3.01E-02 | -3.01E-02 | -100 | 0 | tagging_mirrored_sel |
| SSB_unfished_se | 0.00E+00 | 5.24E+02 | 5.24E+02 | -100 | 0 | tagging_mirrored_sel |
| ForeCatch_2008_val | 5.72E+02 | 4.14E+02 | 1.58E+02 | 38.3118 | 1.38312 | tagging_mirrored_sel |
| ForeCatch_2008_se | 0.00E+00 | 5.72E+01 | 5.72E+01 | -100 | 0 | tagging_mirrored_sel |
In order to match an old model, the compatibility flag needs to be set to 0 and the HCR anchor also set to 0.
Yes, but to test the models with the new way is what I was trying to do. So it's okay if it "fails" since the values would be different. I just want to make sure that what we are seeing when it uses the new way looks correct.
got it. I will run locally and look at the details
my local run with compatibility = 1 and HCR = 3 ran normally and did not produce 0.0 for SR_LN(R0)_se or for fore_catch_2008, So, I cannot tell what went awry for the gha. the HCR=3 did produce a different value for forecast_2008 as expected. I think that for testing, a better value to use for HCR anchor would be 2 to maintain old result
After changes to the r4ss package, the test-r4ss-with-ss3 github action is now passing when I trigger it manually to run on this branch: https://github.com/nmfs-ost/ss3-source-code/actions/workflows/test-r4ss-with-ss3.yml. And I'm also not seeing any issues with the test-ss3-no-est check, which I think is failing because of the HCR anchor 3 vs 2 causing the ForeCatch_2008 to change as you noted elsewhere.
One question on the change to the SS3 code is whether we need to keep the 18., ... 20. numbers preceding the new management quantities added in
https://github.com/nmfs-ost/ss3-source-code/blob/191-spawner-recruit-rebased/SS_readcontrol_330.tpl#L6670-L6685.
There's also an "ignore" quantity that appears in the derived quantities with a estimate and variance but would be nice to either remove or provide a more informative name.
That’s good news. I will remove both next week.
Richard D. Methot Jr.
Stock Assessment Research Scientist (ST)
Northwest Fisheries Science Center
NOAA Fisheries | U.S. Department of Commerce
Office: (425) 666-9893
Mobile: (301) 830-2454
www.fisheries.noaa.gov
On Fri, Jun 27, 2025 at 1:20 PM Ian Taylor @.***> wrote:
iantaylor-NOAA left a comment (nmfs-ost/ss3-source-code#640) https://github.com/nmfs-ost/ss3-source-code/pull/640#issuecomment-3014268249
After changes to the r4ss package, the test-r4ss-with-ss3 github action is now passing when I trigger it manually to run on this branch: https://github.com/nmfs-ost/ss3-source-code/actions/workflows/test-r4ss-with-ss3.yml. And I'm also not seeing any issues with the test-ss3-no-est check, which I think is failing because of the HCR anchor 3 vs 2 causing the ForeCatch_2008 to change as you noted elsewhere.
One question on the change to the SS3 code is whether we need to keep the 18., ... 20. numbers preceding the new management quantities added in
https://github.com/nmfs-ost/ss3-source-code/blob/191-spawner-recruit-rebased/SS_readcontrol_330.tpl#L6670-L6685 . There's also an "ignore" quantity that appears in the derived quantities with a estimate and variance but would be nice to either remove or provide a more informative name.
— Reply to this email directly, view it on GitHub https://github.com/nmfs-ost/ss3-source-code/pull/640#issuecomment-3014268249, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABPV4IC7WPFO7BQ4W7VP57L3FWRRNAVCNFSM6AAAAABRJKPHDCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTAMJUGI3DQMRUHE . You are receiving this because you were mentioned.Message ID: @.***>
@Rick-Methot-NOAA, here's one more thing to look at next week: Thanks for bringing back the KNIFE_AGE_SELECTIVITY_MSY and SLOT_AGE_SELECTIVITY_MSY, tables at the end of Report.sso, but I notice that there are some columns where non-zero values in earlier versions are now zero. The columns are "Exploit", "Y_MSY", "Y_dead", and "Y_ret" in all of the sections (SPR, BTGT and MSY).
@Rick-Methot-NOAA Now that it is running with estimation there are some values that are different than the previous run's values which is making the action fail. Some of these are expected or very small but the Simple_Lorenzen_tv_trend has some big changes that look like they may need attention.
| quantity | value | ref_value | diff | perc_change | ratio | mod_name |
|---|---|---|---|---|---|---|
| Sum_recdevs_like | 4.14E-12 | 4.47E-09 | -4.46E-09 | -9.99E+01 | 0.000926 | growth_timevary |
| Forecast_Recruitment_like | 5.07E-24 | 5.91E-18 | -5.91E-18 | -1.00E+02 | 0.000001 | growth_timevary |
| maxgrad | 4.46E-05 | 9.27E-05 | -4.80E-05 | -5.18E+01 | 0.481532 | growth_timevary |
| SSB_unfished_se | 6.00E+05 | 3.29E+05 | 2.71E+05 | 8.22E+01 | 1.821990 | growth_timevary |
| Catch_like | 5.78E-13 | 5.78E-13 | 1.00E-18 | 1.73E-04 | 1.000002 | platoons |
| Sum_recdevs_like | -4.89E-16 | 1.49E-16 | -6.38E-16 | -4.28E+02 | -3.279074 | platoons |
| maxgrad | 8.59E-05 | 8.30E-05 | 2.86E-06 | 3.44E+00 | 1.034426 | platoons |
| SSB_unfished_se | 4.78E+03 | 5.71E+03 | -9.39E+02 | -1.64E+01 | 0.835735 | platoons |
| SSB_unfished_se | 1.41E+03 | 1.39E+03 | 16.38 | 1.17651 | 1.011770 | Simple |
| SSB_unfished_se | 0.00E+00 | 1.70E+03 | -1699.49 | -100 | 0.000000 | Simple_Lorenzen_tv_trend |
| ForeCatch_2011_se | 0.00E+00 | 4.51E+02 | -450.739 | -100 | 0.000000 | Simple_Lorenzen_tv_trend |
| SSB_unfished_se | 1.41E+03 | 1.39E+03 | 19.58 | 1.40988 | 1.014100 | Simple_NoCPUE |
| SSB_unfished_se | 1.48E+03 | 1.46E+03 | 18.03 | 1.2335 | 1.012330 | Simple_with_Discard |
respond to Ian's request for label changes, but reworded rather that delete one label because indexing of labels is complex
My local run of simple_lorensen gave: Convergence_Level: 6.84275e-06 is_final_gradient Hessian: 154.951 is ln(determinant). Final_phase: 5 N_iterations: 479 total_LogL: 794.96
Here are a few notes on my looking into the differences found in the tests:
-
for likelihoods components that are tiny (e.g. 4.14E-12), it would be good to ignore differences even if the percent change or ratio is above the threshold. I would suggest that any likelihood contribution with absolute value less than 0.001 could be ignored as long as the reference and new values are both below that value.
-
For "Simple_Lorenzen_tv_trend" I get the same estimates and gradient that @Rick-Methot-NOAA reports above, and all of the parameters have reasonable uncertainty estimates. However, the the uncertainty in the derived quantities isn't working, where this warning appears after estimation
Warning: Non-positive variance of sdreport variables: 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, ...and the report files shows 0 in the StdDev column of all the derived quantities:
Label Value StdDev (Val-1.0)/Stddev CumNorm
SSB_Virgin 45831.8 0
SSB_Initial 45831.8 0
SSB_1971 45831.8 0
- For the "platoons" model, the uncertainty in many of the management quantities has decreased by up to about 20% (screenshot below). The new values may well be better than the old values, but it would be good to understand why these quantities would have changed for the platoons model and not the others.
In order to try to narrow down the issue, I changed from N_platoons = 5 to 1 and reran with both 3.30.23.2 and latest build (from today) and I see a similar difference. This model also happens to have implementation error parameters turned on in the forecast,
0.1 # stddev of log(realized catch/target catch) in forecast (set value>0.0 to cause active impl_error)but changing that to 0 also didn't have any impact on the management quantities, so I'm not sure what it could be about this model which is different.
When I turn off dynamic Bzero in the lorenzen control file, then the derived variances work
Okay, after making the change to the lorenzen to not include dynamic B0, the values that have changed are the following:
| quantity | value | ref_value | diff | perc_change | ratio | mod_name |
|---|---|---|---|---|---|---|
| Sum_recdevs_like | 4.14E-12 | 4.47E-09 | -0.000000004 | -99.9074 | 0.00093 | growth_timevary |
| Forecast_Recruitment_like | 5.07E-24 | 5.91E-18 | -5.91E-18 | -99.9999 | 0.0000009 | growth_timevary |
| Parm_priors_like | 0.2915 | 0.2915 | 2.00E-06 | 0.0007 | 1.00001 | growth_timevary |
| maxgrad | 0.0000 | 0.0001 | -0.000048050 | -51.8468 | 0.48153 | growth_timevary |
| SSB_unfished_se | 600018.0000 | 329321.0000 | 270697.000000000 | 82.1985 | 1.82199 | growth_timevary |
| Sum_recdevs_like | -4.89E-16 | 1.49E-16 | -6.38E-16 | -427.9070 | -3.27907 | platoons |
| maxgrad | 0.0001 | 0.0001 | 0.000002858 | 3.4426 | 1.03443 | platoons |
| SSB_unfished_se | 4775.4400 | 5714.0600 | -938.620000000 | -16.4265 | 0.83574 | platoons |
| SSB_unfished_se | 1408.6300 | 1392.2500 | 16.380000000 | 1.1765 | 1.01177 | Simple |
| SSB_unfished_se | 1408.3500 | 1388.7700 | 19.580000000 | 1.4099 | 1.01410 | Simple_NoCPUE |
| SSB_unfished_se | 1479.7300 | 1461.7000 | 18.030000000 | 1.2335 | 1.01233 | Simple_with_Discard |
The growth_timevary we expect, the 3 simple ones have very small changes so I think we are okay there, but the platoons one still has some differences that we don't understand yet. @Rick-Methot-NOAA, could you dig into that model to see if you can figure out why we are seeing these differences?
I suspect that the changes to SSB_unfished_se are because the calculation of that quantity now can involves a call to Equ_SpawnRecr_Result function, rather than the simpler approach. But the simpler existing approach is still used in compatibility mode, so this might require a deep dive to understand it.
@Rick-Methot-NOAA Would you like to proceed with the pre-release while you do that deep dive into that?
Probably yes with the pre-release; that will help find other differences. Please send me the ss_summary.sso file from the simple model run in the gha. Seems different from my local.
The simple model in the gha produces a different result than what I get locally.
gha simple is using recdev option 1 (zero-centered), but it produces a sum_recdevs of 1.37009. How can that happen?
The compare picture attached here has the gha run on the left, my local run in center, and an old run using 3.30.20 on the right.
@Rick-Methot-NOAA, did you add the compatibility lines in the simple model or did you run it just as is?