Alt lorenzen

[nathanvaughan-NOAA]

What issue(s) does this PR address? Describe and add issue numbers, if applicable.

This pull request addresses issue #328 to create a survivorship-based Lorenzen M formulation. This formulation facilitates the existing workflow for SEDAR stock assessments.

What tests have been done? Upload any model input files created for testing in a zip file, if possible.

This modification has been compiled locally and tested on queen triggerfish and scamp with results matching our externally calculated expectations.
Google drive link to scamp example (this example didn't have the final pred M refactoring but I did test that locally)
https://drive.google.com/drive/folders/1U6hY1t32pm8AjxP-aZKYVd-DNYZMslcH?usp=sharing

What tests/review still need to be done? Who can do it, and by when is it needed (ideally)?

I don't believe this is significantly different than the existing formulation in terms of potential errors.
This change does not impact the current formulation so shouldn't cause errors in existing models.

Check which is true. This PR requires:

• no further changes to r4ss
• no further changes to the manual
• no further changes to SSI (the SS3 GUI)
• no further changes to the stock synthesis change log (new features, bug reports)

Describe any changes in r4ss/SS3 manual/SSI/change log that are needed (if not checked):

I will also make a pull request to the doc that takes a first pass at adding a description of this feature but I think it needs discussion.

Additional information (optional):

[Rick-Methot-NOAA]

Hi Nathan.
Given that Lorenzen M is invoked inside of a switch(), perhaps this would be cleaner to do as another block within the switch rather than a IF clause inside the switch.
Rick

[nathanvaughan-NOAA]

Hey Rick,
I'm happy to change this to run via a new case=5 option if you think that would be cleaner, I was just trying to minimize the amount of code I modified. In that case, I will also need to modify the r4ss read/write control options. Do you know of any other places in the code I will need to modify if a new switch option is added?
Cheers,
Nathan

[Rick-Methot-NOAA]

Thanks. I do think it will be cleaner as the code is much different.
related modifications will just be in readcontrol and write_ssnew

[nathanvaughan-NOAA]

Thanks, working on it now.

[nathanvaughan-NOAA]

Code is now updated to include a case 6 option for survivorship-based Lorenzen. Now that it is a new case I modified it to have 2 conditional input values, one for the minimum age and one for the maximum age, this should be much simpler for users to understand. I have also submitted pull requests to update the SS manual with the new option and to r4ss to update the read/write control functions. I have not made any update to the SSGUI so this option will not work with that for now.

[Rick-Methot-NOAA]

So rather than the parameter being M at a specified age, the parameter is now survivorship that is averaged over a range of input reference ages. Correct? So this equation:
Loren_M1=(-log(natMparms(1,gp)));
converts the average annual survivorship during that age range to a M. Correct? Unclear to me why the parameter cannot just be M as it is for other natM options.

[nathanvaughan-NOAA]

The new function takes the total survivorship over the age range and then calculates Loren_M1 as the total mortality over that time (i.e. if 1% survivorship is input as 0.01 then Loren_M1=4.60517). You are correct that this could be specified also in terms of average M over that age range which we did discuss as an option. The rationale for using survivorship as the input was that it was considered more intuitive to discuss with reviewers/SSC/etc. I can always change that though if you think a broader audience will find it more confusing than average M.

[Rick-Methot-NOAA]

I do not understand "total survivorship over the age range"
Is it survivorship from Refage_min to ref_age_max? or survivorship from age 0 to refage_mid, or what.

[nathanvaughan-NOAA]

Yes, survivorship from age_min to age_max.

[Rick-Methot-NOAA]

then:
Loren_M1=(-log(natMparms(1,gp)) / (age_max - age_min));

[nathanvaughan-NOAA]

It doesn't need that division because I need to sum over the Lorenzen scaled values as well on lines 1261:1264. I would just have to do the division by ages again there so it would be a wasted computation.

[Rick-Methot-NOAA]

OK. So adding about 6 comments for the various steps of the computation would help clarity (not that I am always good at that). So first comment would be:
// get M * (delta_age) over the age range for survivorship.

Note that by the parameter being the total survivorship over the age range, if the user changes the age range, then the parameter value also needs rescaled.

Makes me think that greatest clarity is if the M input is the average M over the specified age range. In this way, if age_max = age_min, this Lorenzen implementation provides identical result as the original.

[nathanvaughan-NOAA]

No worries, I did actually add comments in most of the code just not that first line (and the line 1233 comment is now incorrect as I changed that code for the new case 6 version). I can expand the comments though and modify the code to be an average M over the age range.

[Rick-Methot-NOAA]

Thanks.

[Rick-Methot-NOAA]

This looks good Nathan.
Good if you could attach a brief demo showing that this produces same result as previous method.
@nschindler-noaa or @k-doering-NOAA can you adjust the styling? I see missing spaces.
I approve merging after those actions.

[nathanvaughan-NOAA]

Excellent, thanks Rick. I will coordinate with Kathryn tomorrow to adjust the styling and add the demo example.

[k-doering-NOAA]

Styling changes were made - Also, it looks like the checks are passing, but some of them only show up on Nathan's branch but not on the PR: https://github.com/nathanvaughan-NOAA/stock-synthesis/tree/Alt_Lorenzen

The with est run was still going when I typed this, but I think once this passes, this should be good to merge in (if Rick agrees).

[nathanvaughan-NOAA]

Thanks for touching up the code Kathryn. Here is the google drive link where I have added three example cases
1: New case 6 option with reference ages of 3 and 28 (expected use example).
2: New case 6 option with reference ages of 11 and 11 (replicate of existing case 2 use example).
3: Existing case 2 option with reference age of 11 (control to compare with new case 6).

examples 2 and 3 show essentially identical results as expected.

https://drive.google.com/drive/folders/1U6hY1t32pm8AjxP-aZKYVd-DNYZMslcH?usp=sharing

[Rick-Methot-NOAA]

Excellent work Nathan and glad to see the full write-up and examples.
We'll look forward to your next contribution to SS3.
Rick

[Rick-Methot-NOAA]

like the name