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 attempts to use it as a better way to implement generalized sizecomp data

• 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]

I think commented-out lines like these should either be removed or have a comment above them noting what they would do if they were added back

[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.