first commit of rebin

[Rick-Methot-NOAA]

Concisely describe what has been changed/addressed in the pull request.

This PR adds a FUNCTION "rebin" in miscfxn.tpl then uses it as a superior, streamlined way to implement generalized sizecomp data. It corrects the legacy problem with narrow data bins, and corrects a flaw in the accumulation of small population bins.

• Resolves new approach for weight composition - superceded #752

What tests have been done?

Where are the relevant files?

-- - [x] Test files are in the issue. -->

What tests/review still need to be done?

Is there an input change for users to Stock Synthesis?

Additional information (optional).

[Rick-Methot-NOAA]

please do not link copilot until I am ready

[Rick-Methot-NOAA]

get a match to the legacy method using this setup:

5 #_Sizefreq N bins
 2 #_Sizefreq units(1=bio/2=num)
 1 #_Sizefreq scale(1=kg/2=lbs/3=cm/4=inches)
 0.001 #_Sizefreq:  small constant to add to comps
 1 #_Sizefreq number of obs per method
#_Sizefreq bins. one row for each method
#Note: negative value for first bin makes it accumulate all smaller fish vs. truncate small fish
-0.1 1 3 5 7
#_method year month fleet sex part Nsamp <data> 
 1 1971 7 1 3 0 125 1 1 1 1 1 1 1 1 1 1

Beware: lots of temporary echoinput: statements still in the output.

• Next step is to extend capability to all scale and units options
• then test with intentional examples that break the legacy method

[e-perl-NOAA]

@Rick-Methot-NOAA while on my own I don't fully understand what each line of code is doing, I did have Gemini explain it to me and I asked it if it had any suggestions for improvement and it suggested two things:

1. A safety check after dvariable src_bin_width for the situation that if a user provides faulty data where two adjacent source bin edges are identical (making the bin width 0), this will cause a division by zero, resulting in a fatal NaN (Not a Number) error that propagates through the model.

dvariable src_bin_width = s_high - s_low;
if (src_bin_width > 1e-8) { 
    dest_counts[i] += src_counts[j] * (overlap_width / src_bin_width);
}

2. The current rebin function uses a nested loop that compares every destination bin to every source bin. If there are $M$ destination bins and $N$ source bins, this requires $O(M \times N)$ operations. Because size frequency bins are strictly increasing, it is inefficient to check if the 1st source bin overlaps with the 100th destination bin. A fix would be to use a "two-pointer" approach. Since the arrays are sorted, you can track the current source bin and only advance it when you've moved past its upper edge, reducing the time complexity to $O(M + N)$. The following was provided as an example:

FUNCTION dvar_vector rebin(const dvector& src_edges, const dvar_vector& src_counts, const dvar_vector& dest_edges)
{
    dvar_vector dest_counts(1, dest_edges.size() - 1);
    dest_counts.initialize();

    int j_start = 1; // Pointer to track the lowest relevant source bin

    for (int i = 1; i <= dest_counts.size(); i++) {
        dvariable d_low = dest_edges[i];
        dvariable d_high = dest_edges[i + 1];

        // Advance j_start if the source bin is entirely below the current destination bin.
        // Because d_low increases with 'i', j_start only ever moves forward.
        while (j_start <= src_counts.size() && src_edges[j_start + 1] <= d_low) {
            j_start++;
        }

        // Iterate through source bins starting from j_start, but stop as soon 
        // as the source bin is completely above the current destination bin.
        for (int j = j_start; j <= src_counts.size() && src_edges[j] < d_high; j++) {
            dvariable s_low = src_edges[j];
            dvariable s_high = src_edges[j + 1];

            // Calculate the overlap bounds
            dvariable overlap_low = d_low;
            if (s_low > d_low) overlap_low = s_low;
            
            dvariable overlap_high = d_high;
            if (s_high < d_high) overlap_high = s_high;

            // If there is valid overlap, distribute the counts
            if (overlap_low < overlap_high) {
                dvariable overlap_width = overlap_high - overlap_low;
                dvariable src_bin_width = s_high - s_low;
                
                // Safety check: Prevent division by zero if source bin edges are identical
                if (src_bin_width > 1e-8) {
                    dest_counts[i] += src_counts[j] * (overlap_width / src_bin_width);
                }
            }
        }
    }
    
    return (dest_counts);
}

The first suggestion seems logical and harmless enough. The second suggestion, I don't fully understand but is only a change for computational savings.

[Rick-Methot-NOAA]

Thanks Elizabeth. The second suggestion is something I was already thinking about and is what the legacy code was already capable of.

I implemented an efficient search so that each source bin is only addressed once or twice for each destination bin. Sequence from some temporary echoinputs looks like this:
sum_src: 18008.2
in rebin:
dest_edges: 0 47.5267 65.9921 76.8734 85.0034 1481.97 (this is the weight bins converted to cm)
dest_bin 1 src_start 1
1 10 12 store
2 12 14 store
3 14 16 store
4 16 18 store
5 18 20 store
6 20 22 store
7 22 24 store
8 24 26 store
9 26 28 store
10 28 30 store
11 30 32 store
12 32 34 store
13 34 36 store
14 36 38 store
15 38 40 store
16 40 42 store
17 42 44 store
18 44 46 store
19 46 47.5267 store
20 48 47.5267 stop_here: 20
dest_bin 2 src_start 19
19 47.5267 48 store
20 48 50 store
21 50 52 store
22 52 54 store
23 54 56 store
24 56 58 store
25 58 60 store
26 60 62 store
27 62 64 store
28 64 65.9921 store
29 66 65.9921 stop_here: 29
dest_bin 3 src_start 28
28 65.9921 66 store
29 66 68 store
30 68 70 store
31 70 72 store
32 72 74 store
33 74 76 store
34 76 76.8734 store
35 78 76.8734 stop_here: 35
dest_bin 4 src_start 34
34 76.8734 78 store
35 78 80 store
36 80 82 store
37 82 84 store
38 84 85.0034 store
39 86 85.0034 stop_here: 39
dest_bin 5 src_start 38
38 85.0034 86 store
39 86 88 store
40 88 90 store
41 90 92 store
42 92 94 store
43 94 96 store
dest_comps 1206.04 8567.37 6438.2 1575.4 221.159
sum_dest: 18008.2 (which matches sum_src, so all source bins have been allocated to the destination bins).

Next push will be for the streamlined search algorithm.

looping in @N-DucharmeBarth-NOAA who prompted this work long ago.

[iantaylor-NOAA]

@Rick-Methot-NOAA, thanks for moving this forward.

I didn't try to fully understand the algorithm but ran it with the Simple_with_DM_sizefreq model and the results are identical to 3.30.25.1, so that's good enough for me in terms of the calculations.

Regarding efficiency, I ran the model before and after the latest commit and unfortunately things are running slower (based on sample size of 1 with each version).

• 3.30.25.1 model (run with -nohess): 373 iterations, 35 seconds
• before latest commit: 373 iterations, 49 seconds
• after latest commit: 373 iterations, 1 minutes, 33 seconds

Could the echoinput.sso be slowing things down? The rebin versions have 94,000+ lines vs ~2800 in the original (although that doesn't explain the difference caused by the latest commit). In the final version of this PR, it would be good to reduce the echoinput and remove the legacy rebin calculations entirely (perhaps that's already noted somewhere as the plan).

[iantaylor-NOAA]

comment here seems unfinished

[iantaylor-NOAA]

One more comment. I modified the "Simple_with_DM_sizefreq" model to turn off the movement and the time-varying allocation of recruitment among areas which have different growth patterns: control.ss.txt. This causes the Change? column in MGparm_By_Year_after_adjustments to be 0 for all years, but the echoinput suggests that the rebinning is occurring in all years. How hard would it be to rebin only once for a model without time-varying quantities? Or maybe I'm misunderstanding the interaction between changes in size structure due to fishing and the rebinning process.

[Rick-Methot-NOAA]

Thanks for taking a look Ian, but this is still not ready for any real testing. The expected values are still coming from the old method and all my testing has just been unit testing within that code segment. I should have just looped you in with an @, rather that letting it seem like it is ready for review.

[iantaylor-NOAA]

@Rick-Methot-NOAA no problem. It's good to start understanding how this is working even if it's not ready for prime time yet.

[Rick-Methot-NOAA]

The logic of the new rebin method is quite different from the logic of the legacy method. Several internal arrays with the legacy method can be removed before this is merged.
Legacy method:

1. for each poplenbin, calculate the proportion that goes to each of the Szfreq bins and store in a transition matrix.
2. that matrix only needed to be calculated once, unless there was time-varying wt-at-len
3. the complications of changing units and scale was dealt with in the calculation of that transition matrix.
   that matrix creation failed with some narrow bin situations
4. It is a full matrix, so all cells get accessed when it later gets used to convert a lengthcomp expected vector into a szfreq vector.

new rebin method is much more compact:

1. convert the szfreq bin boundaries into cm. This happens for every observation, but could be modified to only occur for the first, or time-varying, of each szfreq method. But it is much simpler than the legacy method that embedded spanning bin calculations into the creation of the transition matrix.
2. if the szfreq obs is in units of weight, then do elem_prod of the length expected value vector and the wt@len vector to get biomass at length. Uses the current wt@len if time-varying.
3. send to rebin the poplen bin boundaries, the szfreq boundaries, and the expected value comp from step 2. Returns the comp in the new bins. Rebin is smart enough to not have to process all poplen bins for each szfreq bin.

[Rick-Methot-NOAA]

obviously not working yet. I'll let you know when ready for real testing

[Rick-Methot-NOAA]

While working on the new approach, it is not clear that the legacy approach was properly implementing the omit small feature correctly. Watch for this in the testing.

[Rick-Methot-NOAA]

Here is my test file which includes 5 permutations on size composition inputs. Included excel sheet shows close match of expected value between legacy and rebin.
szcomp_test.zip
However, the GHA test files include "simple_with_DM_sizefreq" which works for early years, then two of the fleets fail beginning in 1977. Not sure what is going on there.
Any clues Ian? @iantaylor-NOAA
more: The simple with DM_sizefreq example has a lot of complexity, particularly having 2 Gpatterns. I suspect that is the cause of the problem. I'll investigate,

[Rick-Methot-NOAA]

The "simple_with_DM_sizefreq" example is showing a small difference. It is entirely possible that the new, simpler, code is more correct than the legacy custom code. I will do a detailed compare of the expected values of the legacy and rebin outputs to see if I see a clue as to the reason for the difference.

[iantaylor-NOAA]

Thanks for your work on this @Rick-Methot-NOAA. I've been distracted by other projects so haven't looked at the latest changes, but please let you me know if you want help after your detailed comparison. We could also look for other models with size comps to see how big the scale of the differences are elsewhere.

[Rick-Methot-NOAA]

Found it. The results are identical if the first bin is set to -1 to accumulate the tiny fish into first bin. Results differ when that is set to truncate those fish. I suspect there is an errant <= vs < in the code.

[Rick-Methot-NOAA]

The gha test is showing a failure, but this is because the legacy method is slightly wrong.

The attached file shows that the new rebin method does very well at assigning pop len bins to szfreq data bins.
demonstration_szfreq_bins.txt

However, close examination of the szfreq_translation table in the legacy model run shows that both the 24 and 26 cm poplen bins are assigned to the first 26cm data bin. Then examination of the expected szfreq output in compreport.sso shows that this results in the legacy method having a small surplus of fish assigned to the first szfreq data bin (26 cm).

Remaining issues:

• clean-up the code to remove reporting of the deprecated szfreq_translation table. This is not needed because the efficient rebin method does not need it.
• consider if a fix is feasible to allow the rebin method to more precisely assign biomass to bins. The legacy method accounted for fact that fish at smaller edge of a poplen bin have lower body weight than fish at upper portion of the poplen bin. This matters when there is overlap, but the legacy code was rather convoluted to account for the phenomenon.

[Rick-Methot-NOAA]

@iantaylor-NOAA Does r4ss read this output table? I am preparing to deprecate it as it is not populated with the rebin method.
SIZEFREQ_TRANSLATION report:12

[iantaylor-NOAA]

@iantaylor-NOAA Does r4ss read this output table? I am preparing to deprecate it as it is not populated with the rebin method. SIZEFREQ_TRANSLATION report:12

No, it does not read that table. I thought that I had added that functionality at some point, but II don't even see it on the issues list, let alone in the code. So no problem removing the table.

[Rick-Methot-NOAA]

@N-DucharmeBarth-NOAA Can you give this a try with your model configuration?

[N-DucharmeBarth-NOAA]

@N-DucharmeBarth-NOAA Can you give this a try with your model configuration?

Hi @Rick-Methot-NOAA please forgive my ignorance but where can I grab the executable that has this feature implemented?

[iantaylor-NOAA]

@N-DucharmeBarth-NOAA Can you give this a try with your model configuration?

Hi @Rick-Methot-NOAA please forgive my ignorance but where can I grab the executable that has this feature implemented?

For some reason the build-ss3 github action failed to compile the windows version, so I just restarted it.
Once it's complete, hoping that it passes, you should see artifacts at the bottom of the github action page which for the current build is at https://github.com/nmfs-ost/ss3-source-code/actions/runs/35132811441.

[iantaylor-NOAA]

The Windows github action was successful on the 2nd try and here's the direct link to the zip with the executable: https://github.com/nmfs-ost/ss3-source-code/actions/runs/35132811441/artifacts/10473437088.
However, somehow re-running only the windows action has caused the artifacts for the other operating systems to be excluded. Let us know if you need something other than widows for testing this out.

[N-DucharmeBarth-NOAA]

Hi @Rick-Methot-NOAA and @iantaylor-NOAA,

Apologies for the delay in getting back to you on this. I used the attached files to test the new version. The model appears to execute properly now with the new version that you shared with me. Using my previous version 3.30.24.1/ss3_win.exe on the exact same files I received the below screen output. Since the model no longer crashes out at this stage it appears that the fix was successful. For reference, I believe that the source of the initial failure was a mismatch between the population length bins 1cm from 5cm-220cm and the generalized size comp bins 1kg from 1kg to 200kg.
Best,
Nicholas

 reading from starter.ss
 reading from data file
Data read successful 999

 reading forecast file 
 reading from control file
End of control file successful! 999
Begin setting up parameters ...  end of parameter section 
Evaluating biology calculations once ... done
Starting optimization of 'ss3_win' in phase 1 of 7 at Fri Sep 18 07:54:52 2026
matrix column bound exceeded -- column index 1 is higher than 216 in linad99\fvma_acc.cpp:49.
[1] "model run failed"

starter.txt
control.txt
data.txt
forecast.txt

[Rick-Methot-NOAA]

Just reviewed the test example and the match to legacy is not as good as it should be. Hit pause on merging until I investigate.
The difference is only when body weight was being accumulated into the bins. A detailed comparison of the code and intermediate values showed that the old code used the body weight at the mid-length of a length bin, but the new code was using the midweight between the body weights at the bin edges. This can be easily rectified for the length bins that fall entirely within the length range of a sizefreq weight bin. However, the complex logic that I concocted years ago for the overlapping bins will require some further thought as to how to replicate.

[e-perl-NOAA]

@Rick-Methot-NOAA, copilot suggests changes the overlap cases to estimate weight from the midpoint of the overlap in source length space, while still preserving src_mid_wt[j] when the whole source bin lands inside one destination bin and it provided the following code to do so:

@@
-FUNCTION dvar_vector rebin_bio(const dvector& src_edges, const dvar_vector& src_counts, const dvar_vector& dest_edges, const dvar_vector& src_edges_wt, const dvar_vector& dest_edges_wt, const dvar_vector& src_mid_wt)
+FUNCTION dvar_vector rebin_bio(const dvector& src_edges, const dvar_vector& src_counts, const dvar_vector& dest_edges, const dvar_vector& src_edges_wt, const dvar_vector& dest_edges_wt, const dvar_vector& src_mid_wt)
   {
@@
-            dvariable mean_szwt;  // body weight of fish to be allocated from src to dest
-            if (s_low >= d_low)
-            {
-              if (s_high <= d_high)  // src entirely in dest range range
-                {mean_szwt = src_mid_wt[j];}  //  (src_edges_wt[j] + src_edges_wt[j+1]) * 0.5;}
-              else //  s_high > d_high; calc mean size from s_low and d_high
-                {mean_szwt = (src_edges_wt[j] + dest_edges_wt[i+1]) * 0.5;}
-            }
-            else if (s_high <= d_high )  // overlap with s_low < d_low; calc the mean size from d_low and s_high 
-              {mean_szwt = (dest_edges_wt[i] + src_edges_wt[j+1]) * 0.5;}
-            else  // s_high > d_high, so dest bin is internal to src_bin
-              {mean_szwt = (dest_edges_wt[i] + dest_edges_wt[i+1]) * 0.5;}
+            dvariable mean_szwt;  // body weight of fish allocated from src to dest
+            if (s_low >= d_low && s_high <= d_high)
+            {
+              // full source bin goes into this destination bin:
+              // match legacy intent by using body weight at the source-bin mid-length
+              mean_szwt = src_mid_wt[j];
+            }
+            else
+            {
+              // overlap case:
+              // approximate mean biomass for the allocated fraction using the midpoint
+              // of the overlap interval in source-length space, then linearly interpolate
+              // weight between source-bin edge weights
+              dvariable overlap_mid = 0.5 * (overlap_low + overlap_high);
+              dvariable frac_mid = (overlap_mid - s_low) / src_bin_width;
+              if (frac_mid < 0.0) frac_mid = 0.0;
+              if (frac_mid > 1.0) frac_mid = 1.0;
+              mean_szwt = src_edges_wt[j] + frac_mid * (src_edges_wt[j + 1] - src_edges_wt[j]);
+            }
 
-              dest_counts[i] += src_counts[j] * mean_szwt * (overlap_width / src_bin_width);  // add overlap fraction of biomass to dest
+              dest_counts[i] += src_counts[j] * mean_szwt * (overlap_width / src_bin_width);  // add overlap fraction of biomass to dest

You obviously are better equipped to assess if this code is useful or not.

Given the two test models that are "failing" (different enough from legacy to be flagged) use multiple growth patterns, I'm also curious if the only using GP 1 biology for calculating the sizefreq transition matrix (

ss3-source-code/SS_expval.tpl

Lines 689 to 733 in 9e708cb

	// NOTE: wt_len_low is calculated separately for each growth pattern (GPat)
	// but the code below still just uses GPat=1 for calculation of the sizefreq transition matrix
	dvector src_edges(1, nlength + 1); // +1 to create a dummy bin holding upper edge of the actual last bin
	src_edges(1, nlength) = len_bins(1, nlength);
	src_edges(nlength + 1) = src_edges(nlength) + (src_edges(nlength) - src_edges(nlength - 1));
	dvar_vector dest_edges(1,SzFreq_Nbins(SzFreqMethod) + 1); // +1 to create a dummy bin holding upper edge of the actual last bin
	dvar_vector dest_edges_wt(1,SzFreq_Nbins(SzFreqMethod) + 1);
	dest_edges.initialize();
	for (ibin = 1; ibin <= SzFreq_Nbins(SzFreqMethod); ibin++) {
	dest_edges(ibin) = SzFreq_bins(SzFreqMethod, ibin);
	dest_edges_wt(ibin) = dest_edges(ibin); // save in raw units
	}
	dest_edges(SzFreq_Nbins(SzFreqMethod) + 1) = 99999.;
	dest_edges_wt(ibin) = dest_edges_wt(ibin - 1) + (dest_edges_wt(ibin - 1) - dest_edges_wt(ibin - 2));
	if(SzFreq_Omit_Small(SzFreqMethod) < 0.0)
	{
	dest_edges(1) = 0.; // so src_bins smaller than first actual dest_bin are included in first bin
	}
	dvar_vector dest_comp(1,SzFreq_Nbins(SzFreqMethod)); // destination composition. one-sex at a time
	// convert dest_edges to cm
	// as in the legacy method, this conversion is sex and season specific, but only uses Gpat = 1
	switch (SzFreq_scale(SzFreqMethod)) // biomass vs. numbers are accumulated in the bins
	{
	case (1): // scale is kg
	{
	// fall thru to calcs in the lb case
	}
	case (2): // scale is lb, or fall thru from kg
	{
	// bin boundary in lb have already been converted to kg
	// wt_len(s, gp) is potentially by season and by Gpat within sex. Here, only Gpat==1 will be used. It will be sex-specific
	dvariable wtlen_a;
	dvariable wtlen_b;
	GPat = 1; // this means that any (unusual) use of multiple Gpatterns will still use biology of GPat 1 for szfreq purposes
	// female parms are at 1, 2; males at 7,8
	int parm_loc = 1 + (gg - 1) * 6;
	wtlen_a = wtlen_p(GPat, parm_loc);
	wtlen_b = wtlen_p(GPat, parm_loc + 1);
	if (MGparm_seas_effects(1) > 0 || MGparm_seas_effects(2) > 0) // get seasonal effect on wtlen parameters
	{
	wtlen_a *= wtlen_seas(s, GPat, parm_loc);
	wtlen_b *= wtlen_seas(s, GPat, parm_loc + 1);
	}
	dest_edges = exp((log(dest_edges)-log(wtlen_a))/wtlen_b); // W=a*L^b; so L = exp((ln(W)-ln(a))/B)

) is the real culprit for why the two models are being flagged.

[Rick-Methot-NOAA]

Thanks Elizabeth; I was hoping you'd loop in CoPilot. :)

I thought the failing models were only due to the legacy being wrong regarding initial bin calculation, but I will now look deeper into it.

Rick

[Rick-Methot-NOAA]

the GHA run seems to be failing because the model: simple_with_DM_sizefreq is getting a 0.0 for SR_ln(R0)_se.
When I run that same test model locally, I get a match to the base model.
Nevertheless, the CoPilot code introduced some IF() statements dependent on calculated quantities, and these can mess with gradients which can cause Hessian to fail and result in se=0.0.
I will now try to replace those IF() with logistic joiners.
RIck

[e-perl-NOAA]

I think this needs to be changed to the following to get rid of the build warning:

FUNCTION dvar_vector rebin_bio(
  const dvector& src_edges,
  const dvar_vector& src_counts,
  const dvar_vector& dest_edges,
  const dvar_vector& src_edges_wt,
  const dvar_vector& src_mid_wt)

[e-perl-NOAA]

I also think this needs to be changed to the following to get rid of the build warning:

dest_comp = rebin_bio(src_edges, pass_comp, dest_edges, pass_src_wt, wt_len(s, GPat));

[e-perl-NOAA]

In my conversations with copilot, aside from the changes that I suggested about the build warnings, the rest of the changes that are causing the tests to fail, are all to be expected based on the intent of the changes in this PR. If you agree with the changes I proposed to address the build warnings and implement those, I will add copilot as a reviewer to do just one more spot check, and then I think we will be good to go.

[Rick-Methot-NOAA]

@e-perl-NOAA Wht is the build warning that your changes would resolve? I am much more worried about the se=0 for the models that use sizefreq data. I do not get that when I run locally.