From e880c0cc2335fa92ec1fc9c681290670a577d4bf Mon Sep 17 00:00:00 2001 From: Rick-Methot-NOAA Date: Fri, 14 Aug 2026 16:29:18 -0700 Subject: [PATCH 1/8] first commit of rebin --- SS_miscfxn.tpl | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/SS_miscfxn.tpl b/SS_miscfxn.tpl index 76964cb3..79ec66b4 100644 --- a/SS_miscfxn.tpl +++ b/SS_miscfxn.tpl @@ -120,3 +120,42 @@ FUNCTION dvariable Comp_logL_Dirichlet(const double& Nsamp, const dvariable& dir logL = sum(gammln(Nsamp * obs_comp + dirichlet_Parm * exp_comp)) - sum(gammln(dirichlet_Parm * exp_comp)); return (logL); } + + +FUNCTION dvector rebin(const dvector& src_edges, const dvector& src_counts, const dvector& dest_edges) + { + /* + This implementation takes two vectors representing the boundaries (edges) of the source + and destination bins, and one vector for the source counts. + + Rebins frequency data from one set of boundaries to another. + @param src_edges Boundaries of the original bins (size N+1). + @param src_counts Frequency/counts in original bins (size N). + @param dest_edges Boundaries of the new bins (size M+1). + @return Vector of rebinned frequency data (size M). + */ + + dvector dest_counts(dest_edges.size() - 1, 0.0); + for (int i = 0; i < dest_counts.size(); ++i) { + double d_low = dest_edges[i]; + double d_high = dest_edges[i + 1]; + + for (int j = 0; j < src_counts.size(); ++j) { + double s_low = src_edges[j]; + double s_high = src_edges[j + 1]; + + // Calculate the overlap between [d_low, d_high] and [s_low, s_high] + double overlap_low = max(d_low, s_low); + double overlap_high = min(d_high, s_high); + + if (overlap_low < overlap_high) { + double overlap_width = overlap_high - overlap_low; + double src_bin_width = s_high - s_low; + + // Distribute source count proportionally to the overlap area + dest_counts[i] += src_counts[j] * (overlap_width / src_bin_width); + } + } + } + return (dest_counts); + } \ No newline at end of file From f2fe2b733d6bd1b6c43f4c774ad8d4d25d744a56 Mon Sep 17 00:00:00 2001 From: Rick-Methot-NOAA Date: Wed, 26 Aug 2026 16:55:17 -0700 Subject: [PATCH 2/8] rebin method matches legacy for one configuration --- SS_expval.tpl | 390 +++++++++++++++++++++++++++++++++++++++++++- SS_miscfxn.tpl | 36 ++-- SS_readdata_330.tpl | 8 +- 3 files changed, 414 insertions(+), 20 deletions(-) diff --git a/SS_expval.tpl b/SS_expval.tpl index df6fb5eb..d5b7d798 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -666,13 +666,15 @@ FUNCTION void Get_expected_values(const int y, const int t); { if (SzFreq_Nmeth > 0) // have some sizefreq data { - if (have_data(ALK_time, f, data_type, 0) > 0) { for (j = 1; j <= have_data(ALK_time, f, data_type, 0); j++) // loop all obs of this type { iobs = have_data(ALK_time, f, data_type, j); // observation index SzFreqMethod = SzFreq_obs_hdr(iobs, 6); + + { + echoinput<<" begin szfreq legacy method"< 0) // first occurrence of this method at this time is with fleet = f { @@ -699,7 +701,328 @@ FUNCTION void Get_expected_values(const int y, const int t); // 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 + echoinput<<" ready "<= z2) + { + warnstream << "max population size " << wt_len_low(s, 1, z1) << " is less than first data bin " << SzFreq_bins(SzFreqMethod, 1) << " for SzFreqMethod " << SzFreqMethod; + write_message (FATAL, 0); // EXIT! + } + + if (wt_len_low(s, 1, nlength2) < SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod))) + { + warnstream << "max population size " << wt_len_low(s, 1, nlength2) << " is less than max data bin " << SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod)) << " for SzFreqMethod " << SzFreqMethod; + write_message (FATAL, 0); // EXIT! + } + for (z = z1; z <= z2; z++) + { + if (ibin == SzFreq_Nbins(SzFreqMethod)) + { + SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); + } + else + { + if (wt_len_low(s, 1, z) >= topbin) + { + ibin++; + ibinsave++; + } + if (ibin > 1) + { + botbin = SzFreq_bins2(SzFreqMethod, ibin); + } + if (ibin == SzFreq_Nbins(SzFreqMethod)) + { + SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); + topbin = 99999.; + } + else + { + topbin = SzFreq_bins2(SzFreqMethod, ibin + 1); + if (wt_len_low(s, 1, z) >= botbin && wt_len_low(s, 1, z + 1) <= topbin) + { + SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); + } + else + { + temp = (wt_len_low(s, 1, z + 1) - topbin) / wt_len_fd(s, 1, z); // frac in pop bin above (data bin +1) + temp1 = wt_len_low(s, 1, z) + (1. - temp * 0.5) * wt_len_fd(s, 1, z); // approx body wt for these fish + temp2 = wt_len_low(s, 1, z) + (1. - temp) * 0.5 * wt_len_fd(s, 1, z); // approx body wt for fish below + SzFreqTrans(SzFreqMethod_seas, z, ibinsave + 1) = temp * temp1; + SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = (1. - temp) * temp2; + } + } + } + } + if (SzFreq_scale(SzFreqMethod) == 2 && gg == gender) // convert to pounds + { + SzFreqTrans(SzFreqMethod_seas) /= 0.4536; + } + break; + } // end of units in biomass + // NOTE: even though the transition matrix is currently in units of biomass distribution, there is no need to + // normalize to sum to 1.0 here because the normalization will occur after it gets used to create SzFreq_exp + case (2): // units are numbers + { + echoinput<<" ibin "<= z2) + { + warnstream << "max population size " << wt_len_low(s, 1, z1) << " is less than first data bin " << SzFreq_bins(SzFreqMethod, 1) << " for SzFreqMethod " << SzFreqMethod; + write_message (FATAL, 0); // EXIT! + } + if (wt_len_low(s, 1, nlength2) < SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod))) + { + warnstream << "max population size " << wt_len_low(s, 1, nlength2) << " is less than max data bin " << SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod)) << " for SzFreqMethod " << SzFreqMethod; + write_message (FATAL, 0); // EXIT! + } + for (z = z1; z <= z2; z++) + { +// echoinput<<"ibin "<= topbin) + { + ibin++; + ibinsave++; + } + if (ibin > 1) + { + botbin = SzFreq_bins2(SzFreqMethod, ibin); + } + if (ibin == SzFreq_Nbins(SzFreqMethod)) // checkup<<" got to last ibin, so put rest of popbins here"<= botbin && wt_len_low(s, 1, z + 1) <= topbin) //checkup<<" pop inside dat, put here"<= topbin) + { + ibin++; + ibinsave++; + } //checkup<<" incr ibin "< 1) + { + botbin = SzFreq_bins2(SzFreqMethod, ibin); + } + if (ibin == SzFreq_Nbins(SzFreqMethod)) // checkup<<" got to last ibin, so put rest of popbins here"<= botbin && len_bins2(z + 1) <= topbin) //checkup<<" pop inside dat, put here"< 0.0) + { + SzFreq_exp(iobs)(1, k) += SzFreq_mincomp(SzFreqMethod); + SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); + } + break; + } + case (1): // female only + { + SzFreq_exp(iobs)(k + 1, 2 * k) = 0.0; // zero out the males so will not interfere with data generation + SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); + if (SzFreq_mincomp(SzFreqMethod) > 0.0) + { + SzFreq_exp(iobs)(1, k) += SzFreq_mincomp(SzFreqMethod); + SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); + } + break; + } + case (2): // male only + { + ibin = SzFreq_obs_hdr(iobs, 7); + SzFreq_exp(iobs)(1, ibin - 1) = 0.0; // zero out the females so will not interfere with data generation + SzFreq_exp(iobs)(ibin, k) /= sum(SzFreq_exp(iobs)(ibin, k)); + if (SzFreq_mincomp(SzFreqMethod) > 0.0) + { + SzFreq_exp(iobs)(ibin, k) += SzFreq_mincomp(SzFreqMethod); + SzFreq_exp(iobs)(ibin, k) /= sum(SzFreq_exp(iobs)(ibin, k)); + } + break; + } + case (3): // each gender + { + SzFreq_exp(iobs) /= sum(SzFreq_exp(iobs)); + if (SzFreq_mincomp(SzFreqMethod) > 0.0) + { + SzFreq_exp(iobs) += SzFreq_mincomp(SzFreqMethod); + SzFreq_exp(iobs) /= sum(SzFreq_exp(iobs)); + } + break; + } + } // end gender switch + } // end have 2 genders + else + { + k = SzFreq_obs_hdr(iobs, 8); // max bins for this method + SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); + if (SzFreq_mincomp(SzFreqMethod) > 0.0) + { + SzFreq_exp(iobs)(1, k) += SzFreq_mincomp(SzFreqMethod); + SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); + } + } + #ifdef DO_ONCE + if (do_once == 1) + echoinput << y << " " << f << " szfreq_exp_final " << SzFreq_exp(iobs) << endl; + #endif + + } // end szfreqtrans method + + { // begin rebin method + echoinput<<"begin rebin method"< 0) // first occurrence of this method at this time is with fleet = f + { + if (do_once == 1 || (MG_active(3) > 0 && (timevary_MG(y, 3) > 0))) // calc matrix because wtlen parameters have changed + { + for (gg = 1; gg <= gender; gg++) + { + if (gg == 1) + { + z1 = 1; + z2 = nlength; + ibin = 0; + ibinsave = 0; + } // female + else + { + z1 = nlength1; + z2 = nlength2; + ibin = 0; + ibinsave = SzFreq_Nbins(SzFreqMethod); + } // male + topbin = 0.; + botbin = 0.; + + // 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 + /* switch (SzFreq_units(SzFreqMethod)) // biomass vs. numbers are accumulated in the bins { case (1): // units are biomass, so accumulate body weight into the bins; Assume that bin demarcations are also in biomass @@ -894,10 +1217,68 @@ FUNCTION void Get_expected_values(const int y, const int t); if (docheckup == 1 && gg == gender) echoinput << " sizefreq trans_matrix: method/season " << SzFreqMethod << " / " << s << endl << trans(SzFreqTrans(SzFreqMethod_seas)) << endl; - } // end gender loop - } // end needing to calc the matrix because it may have changed - } // done calculating the SzFreqTransition matrix for this method + */ + 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)); + echoinput << "src_edges "< d_low) overlap_low = s_low; + dvariable overlap_high = d_high; + if(s_high < d_high) overlap_high = s_high; +// = dmin(d_high, s_high); if (overlap_low < overlap_high) { - double overlap_width = overlap_high - overlap_low; - double src_bin_width = s_high - s_low; + dvariable overlap_width = overlap_high - overlap_low; + dvariable src_bin_width = s_high - s_low; // Distribute source count proportionally to the overlap area +// echoinput << " count: "< 0) echoinput << SzFreq_nobs << " Sizefreq N obs per method" << endl; ivector SzFreq_Nbins_seas_g(1,SzFreq_Nmeth*nseas); // array dimensioner used only for the SzFreqTrans array - ivector SzFreq_Nbins3(1,SzFreq_Nmeth); // doubles the Nbins if gender==2 + ivector SzFreq_Nbins3(1,SzFreq_Nmeth); // values will be doubled if gender==2 int SzFreqMethod_seas; ivector Comp_Err_Sz(1,SzFreq_Nmeth); ivector Comp_Err_Sz2(1,SzFreq_Nmeth); @@ -3196,7 +3196,6 @@ } } // clang-format off - echoinput<<"here"< Date: Thu, 27 Aug 2026 15:21:56 -0700 Subject: [PATCH 3/8] revise rebin for efficiency --- SS_miscfxn.tpl | 64 +++++++++++++++++++++++++++------------------ SS_write_report.tpl | 2 +- 2 files changed, 39 insertions(+), 27 deletions(-) diff --git a/SS_miscfxn.tpl b/SS_miscfxn.tpl index 2082c3f1..f67c059d 100644 --- a/SS_miscfxn.tpl +++ b/SS_miscfxn.tpl @@ -133,37 +133,49 @@ FUNCTION dvar_vector rebin(const dvector& src_edges, const dvar_vector& src_coun @param src_counts Frequency/counts in original bins (size N). @param dest_edges Boundaries of the new bins (size M+1). @return Vector of rebinned frequency data (size M). + the original code searched all source bins for each destination bins. + Here the ordered characteristic of the bins allows for the search for bin (i+1) to continue from search for bin(i) */ dvar_vector dest_counts(1, dest_edges.size() - 1); // size to leave off the topbin bounary dest_counts.initialize(); echoinput<<"in rebin: "< d_low) overlap_low = s_low; - dvariable overlap_high = d_high; - if(s_high < d_high) overlap_high = s_high; -// = dmin(d_high, s_high); - - if (overlap_low < overlap_high) { - dvariable overlap_width = overlap_high - overlap_low; - dvariable src_bin_width = s_high - s_low; - - // Distribute source count proportionally to the overlap area -// echoinput << " count: "< dest_edges[i]) overlap_low = src_edges[j]; + if(src_edges[j + 1] < dest_edges[i + 1]) overlap_high = src_edges[j + 1]; +// echoinput<<" "< Date: Fri, 28 Aug 2026 15:55:30 -0700 Subject: [PATCH 4/8] legacy code gone; echos first call only; ready for testing; needs indent clean-up --- SS_biofxn.tpl | 3 +- SS_expval.tpl | 381 +++++++++++++------------------------------------ SS_miscfxn.tpl | 11 -- 3 files changed, 102 insertions(+), 293 deletions(-) diff --git a/SS_biofxn.tpl b/SS_biofxn.tpl index 00acbfc8..e0217d73 100644 --- a/SS_biofxn.tpl +++ b/SS_biofxn.tpl @@ -1532,7 +1532,8 @@ FUNCTION void get_wtlen() wt_len(s, gp).shift(1); } - // SS_Label_Info_19.2.3 #calculate first diff of wt_len for use in generalized sizp comp bin calculations + // SS_Label_Info_19.2.3 #calculate first diff of wt_len for use in generalized size comp bin calculations + // this may no longer be needed with the rebin method if (gg == gender) { wt_len2_sq(s, GPat) = elem_prod(wt_len2(s, GPat), wt_len2(s, GPat)); diff --git a/SS_expval.tpl b/SS_expval.tpl index d5b7d798..4a2474a4 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -672,7 +672,7 @@ FUNCTION void Get_expected_values(const int y, const int t); { iobs = have_data(ALK_time, f, data_type, j); // observation index SzFreqMethod = SzFreq_obs_hdr(iobs, 6); - + /* { echoinput<<" begin szfreq legacy method"< 0) // first occurrence of this method at this time is with fleet = f - { - if (do_once == 1 || (MG_active(3) > 0 && (timevary_MG(y, 3) > 0))) // calc matrix because wtlen parameters have changed - { - for (gg = 1; gg <= gender; gg++) - { - if (gg == 1) - { - z1 = 1; - z2 = nlength; - ibin = 0; - ibinsave = 0; - } // female - else - { - z1 = nlength1; - z2 = nlength2; - ibin = 0; - ibinsave = SzFreq_Nbins(SzFreqMethod); - } // male - topbin = 0.; - botbin = 0.; - - // 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 - /* - switch (SzFreq_units(SzFreqMethod)) // biomass vs. numbers are accumulated in the bins - { - case (1): // units are biomass, so accumulate body weight into the bins; Assume that bin demarcations are also in biomass - { - if (SzFreq_Omit_Small(SzFreqMethod) == 1) - { - while (wt_len_low(s, 1, z1 + 1) < SzFreq_bins(SzFreqMethod, 1) && z1 < z2) - { - z1++; - } - } // ignore tiny fish - if (z1 + 1 >= z2) - { - warnstream << "max population size " << wt_len_low(s, 1, z1) << " is less than first data bin " << SzFreq_bins(SzFreqMethod, 1) << " for SzFreqMethod " << SzFreqMethod; - write_message (FATAL, 0); // EXIT! - } - - if (wt_len_low(s, 1, nlength2) < SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod))) - { - warnstream << "max population size " << wt_len_low(s, 1, nlength2) << " is less than max data bin " << SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod)) << " for SzFreqMethod " << SzFreqMethod; - write_message (FATAL, 0); // EXIT! - } - - for (z = z1; z <= z2; z++) - { - if (ibin == SzFreq_Nbins(SzFreqMethod)) - { - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); - } - else - { - if (wt_len_low(s, 1, z) >= topbin) - { - ibin++; - ibinsave++; - } - if (ibin > 1) - { - botbin = SzFreq_bins2(SzFreqMethod, ibin); - } - if (ibin == SzFreq_Nbins(SzFreqMethod)) - { - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); - topbin = 99999.; - } - else - { - topbin = SzFreq_bins2(SzFreqMethod, ibin + 1); - if (wt_len_low(s, 1, z) >= botbin && wt_len_low(s, 1, z + 1) <= topbin) - { - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); - } - else - { - temp = (wt_len_low(s, 1, z + 1) - topbin) / wt_len_fd(s, 1, z); // frac in pop bin above (data bin +1) - temp1 = wt_len_low(s, 1, z) + (1. - temp * 0.5) * wt_len_fd(s, 1, z); // approx body wt for these fish - temp2 = wt_len_low(s, 1, z) + (1. - temp) * 0.5 * wt_len_fd(s, 1, z); // approx body wt for fish below - SzFreqTrans(SzFreqMethod_seas, z, ibinsave + 1) = temp * temp1; - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = (1. - temp) * temp2; - } - } - } - } - if (SzFreq_scale(SzFreqMethod) == 2 && gg == gender) // convert to pounds - { - SzFreqTrans(SzFreqMethod_seas) /= 0.4536; - } - break; - } // end of units in biomass - // NOTE: even though the transition matrix is currently in units of biomass distribution, there is no need to - // normalize to sum to 1.0 here because the normalization will occur after it gets used to create SzFreq_exp - - case (2): // units are numbers - { - if (SzFreq_scale(SzFreqMethod) <= 2) // bin demarcations are in weight units (1=kg, 2=lbs), so uses wt_len to compare to bins - { - if (SzFreq_Omit_Small(SzFreqMethod) == 1) - { - while (wt_len_low(s, 1, z1 + 1) < SzFreq_bins(SzFreqMethod, 1) && z1 < z2) - { - z1++; - } - } // ignore tiny fish - if (z1 + 1 >= z2) - { - warnstream << "max population size " << wt_len_low(s, 1, z1) << " is less than first data bin " << SzFreq_bins(SzFreqMethod, 1) << " for SzFreqMethod " << SzFreqMethod; - write_message (FATAL, 0); // EXIT! - } - if (wt_len_low(s, 1, nlength2) < SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod))) - { - warnstream << "max population size " << wt_len_low(s, 1, nlength2) << " is less than max data bin " << SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod)) << " for SzFreqMethod " << SzFreqMethod; - write_message (FATAL, 0); // EXIT! - } - - for (z = z1; z <= z2; z++) - { - if (ibin == SzFreq_Nbins(SzFreqMethod)) - { - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = 1.; - } //checkup<<" got to last ibin, so put rest of popbins here"<= topbin) - { - ibin++; - ibinsave++; - } //checkup<<" incr ibin "< 1) - { - botbin = SzFreq_bins2(SzFreqMethod, ibin); - } - if (ibin == SzFreq_Nbins(SzFreqMethod)) // checkup<<" got to last ibin, so put rest of popbins here"<= botbin && wt_len_low(s, 1, z + 1) <= topbin) //checkup<<" pop inside dat, put here"< 0) // first occurrence of this method at this time is with fleet = f + // check to see if multiple fleets using szfreq will still get included + { + for (gg = 1; gg <= gender; gg++) + { + if (gg == 1) + { + z1 = 1; + z2 = nlength; + } // female + else + { + z1 = nlength1; + z2 = nlength2; + } // male - for (z = z1; z <= z2; z++) - { - if (ibin == SzFreq_Nbins(SzFreqMethod)) - { - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = 1.; - } //checkup<<" got to last ibin, so put rest of popbins here"<= topbin) - { - ibin++; - ibinsave++; - } //checkup<<" incr ibin "< 1) - { - botbin = SzFreq_bins2(SzFreqMethod, ibin); - } - if (ibin == SzFreq_Nbins(SzFreqMethod)) // checkup<<" got to last ibin, so put rest of popbins here"<= botbin && len_bins2(z + 1) <= topbin) //checkup<<" pop inside dat, put here"< 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); + break; + } + case (3): // scale is in cm + { + break; + } + case 4: // scale is in inches. Already converted to cm in readdata. Relic feature for some ancient recreational fishery data + { + break; + } + } // end calc of new dest_edges + if (do_once == 1) echoinput << "sex: " << gg << " dest_edges_in_cm "< dest_edges[i]) overlap_low = src_edges[j]; if(src_edges[j + 1] < dest_edges[i + 1]) overlap_high = src_edges[j + 1]; -// echoinput<<" "< Date: Tue, 1 Sep 2026 15:58:22 -0700 Subject: [PATCH 5/8] Gemini's version was better than mine --- SS_expval.tpl | 359 +++---------------------------------------------- SS_miscfxn.tpl | 62 ++++++++- 2 files changed, 73 insertions(+), 348 deletions(-) diff --git a/SS_expval.tpl b/SS_expval.tpl index 4a2474a4..46be4329 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -672,332 +672,8 @@ FUNCTION void Get_expected_values(const int y, const int t); { iobs = have_data(ALK_time, f, data_type, j); // observation index SzFreqMethod = SzFreq_obs_hdr(iobs, 6); - /* - { - echoinput<<" begin szfreq legacy method"< 0) // first occurrence of this method at this time is with fleet = f - { - if (do_once == 1 || (MG_active(3) > 0 && (timevary_MG(y, 3) > 0))) // calc matrix because wtlen parameters have changed - { - for (gg = 1; gg <= gender; gg++) - { - if (gg == 1) - { - z1 = 1; - z2 = nlength; - ibin = 0; - ibinsave = 0; - } // female - else - { - z1 = nlength1; - z2 = nlength2; - ibin = 0; - ibinsave = SzFreq_Nbins(SzFreqMethod); - } // male - topbin = 0.; - botbin = 0.; - - // 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 - echoinput<<" ready "<= z2) - { - warnstream << "max population size " << wt_len_low(s, 1, z1) << " is less than first data bin " << SzFreq_bins(SzFreqMethod, 1) << " for SzFreqMethod " << SzFreqMethod; - write_message (FATAL, 0); // EXIT! - } - - if (wt_len_low(s, 1, nlength2) < SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod))) - { - warnstream << "max population size " << wt_len_low(s, 1, nlength2) << " is less than max data bin " << SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod)) << " for SzFreqMethod " << SzFreqMethod; - write_message (FATAL, 0); // EXIT! - } - for (z = z1; z <= z2; z++) - { - if (ibin == SzFreq_Nbins(SzFreqMethod)) - { - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); - } - else - { - if (wt_len_low(s, 1, z) >= topbin) - { - ibin++; - ibinsave++; - } - if (ibin > 1) - { - botbin = SzFreq_bins2(SzFreqMethod, ibin); - } - if (ibin == SzFreq_Nbins(SzFreqMethod)) - { - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); - topbin = 99999.; - } - else - { - topbin = SzFreq_bins2(SzFreqMethod, ibin + 1); - if (wt_len_low(s, 1, z) >= botbin && wt_len_low(s, 1, z + 1) <= topbin) - { - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = wt_len2(s, 1, z); - } - else - { - temp = (wt_len_low(s, 1, z + 1) - topbin) / wt_len_fd(s, 1, z); // frac in pop bin above (data bin +1) - temp1 = wt_len_low(s, 1, z) + (1. - temp * 0.5) * wt_len_fd(s, 1, z); // approx body wt for these fish - temp2 = wt_len_low(s, 1, z) + (1. - temp) * 0.5 * wt_len_fd(s, 1, z); // approx body wt for fish below - SzFreqTrans(SzFreqMethod_seas, z, ibinsave + 1) = temp * temp1; - SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = (1. - temp) * temp2; - } - } - } - } - if (SzFreq_scale(SzFreqMethod) == 2 && gg == gender) // convert to pounds - { - SzFreqTrans(SzFreqMethod_seas) /= 0.4536; - } - break; - } // end of units in biomass - // NOTE: even though the transition matrix is currently in units of biomass distribution, there is no need to - // normalize to sum to 1.0 here because the normalization will occur after it gets used to create SzFreq_exp - - case (2): // units are numbers - { - echoinput<<" ibin "<= z2) - { - warnstream << "max population size " << wt_len_low(s, 1, z1) << " is less than first data bin " << SzFreq_bins(SzFreqMethod, 1) << " for SzFreqMethod " << SzFreqMethod; - write_message (FATAL, 0); // EXIT! - } - if (wt_len_low(s, 1, nlength2) < SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod))) - { - warnstream << "max population size " << wt_len_low(s, 1, nlength2) << " is less than max data bin " << SzFreq_bins(SzFreqMethod, SzFreq_Nbins(SzFreqMethod)) << " for SzFreqMethod " << SzFreqMethod; - write_message (FATAL, 0); // EXIT! - } - for (z = z1; z <= z2; z++) - { -// echoinput<<"ibin "<= topbin) - { - ibin++; - ibinsave++; - } - if (ibin > 1) - { - botbin = SzFreq_bins2(SzFreqMethod, ibin); - } - if (ibin == SzFreq_Nbins(SzFreqMethod)) // checkup<<" got to last ibin, so put rest of popbins here"<= botbin && wt_len_low(s, 1, z + 1) <= topbin) //checkup<<" pop inside dat, put here"<= topbin) - { - ibin++; - ibinsave++; - } //checkup<<" incr ibin "< 1) - { - botbin = SzFreq_bins2(SzFreqMethod, ibin); - } - if (ibin == SzFreq_Nbins(SzFreqMethod)) // checkup<<" got to last ibin, so put rest of popbins here"<= botbin && len_bins2(z + 1) <= topbin) //checkup<<" pop inside dat, put here"< 0.0) - { - SzFreq_exp(iobs)(1, k) += SzFreq_mincomp(SzFreqMethod); - SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); - } - break; - } - case (1): // female only - { - SzFreq_exp(iobs)(k + 1, 2 * k) = 0.0; // zero out the males so will not interfere with data generation - SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); - if (SzFreq_mincomp(SzFreqMethod) > 0.0) - { - SzFreq_exp(iobs)(1, k) += SzFreq_mincomp(SzFreqMethod); - SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); - } - break; - } - case (2): // male only - { - ibin = SzFreq_obs_hdr(iobs, 7); - SzFreq_exp(iobs)(1, ibin - 1) = 0.0; // zero out the females so will not interfere with data generation - SzFreq_exp(iobs)(ibin, k) /= sum(SzFreq_exp(iobs)(ibin, k)); - if (SzFreq_mincomp(SzFreqMethod) > 0.0) - { - SzFreq_exp(iobs)(ibin, k) += SzFreq_mincomp(SzFreqMethod); - SzFreq_exp(iobs)(ibin, k) /= sum(SzFreq_exp(iobs)(ibin, k)); - } - break; - } - case (3): // each gender - { - SzFreq_exp(iobs) /= sum(SzFreq_exp(iobs)); - if (SzFreq_mincomp(SzFreqMethod) > 0.0) - { - SzFreq_exp(iobs) += SzFreq_mincomp(SzFreqMethod); - SzFreq_exp(iobs) /= sum(SzFreq_exp(iobs)); - } - break; - } - } // end gender switch - } // end have 2 genders - else - { - k = SzFreq_obs_hdr(iobs, 8); // max bins for this method - SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); - if (SzFreq_mincomp(SzFreqMethod) > 0.0) - { - SzFreq_exp(iobs)(1, k) += SzFreq_mincomp(SzFreqMethod); - SzFreq_exp(iobs)(1, k) /= sum(SzFreq_exp(iobs)(1, k)); - } - } - #ifdef DO_ONCE - if (do_once == 1) - echoinput << y << " " << f << " szfreq_exp_final " << SzFreq_exp(iobs) << endl; - #endif - - } // end szfreqtrans method - - - */ - { // begin rebin method - SzFreqMethod_seas = nseas * (SzFreqMethod - 1) + s; // index that combines sizefreqmethod and season and used in SzFreqTrans + SzFreqMethod_seas = nseas * (SzFreqMethod - 1) + s; // index that combines sizefreqmethod and season and used in SzFreqTrans if (SzFreq_obs_hdr(iobs, 9) > 0) // first occurrence of this method at this time is with fleet = f // check to see if multiple fleets using szfreq will still get included { @@ -1025,8 +701,11 @@ FUNCTION void Get_expected_values(const int y, const int t); dest_edges(ibin) = SzFreq_bins(SzFreqMethod, ibin); } if(SzFreq_Omit_Small(SzFreqMethod) < 0.0) - {dest_edges(1) = 0.;} // so all small fish go into this bin - dest_edges(SzFreq_Nbins(SzFreqMethod) + 1) = 99999.; // dest_edges(SzFreq_Nbins(SzFreqMethod)) + (dest_edges(SzFreq_Nbins(SzFreqMethod)) - dest_edges(SzFreq_Nbins(SzFreqMethod) - 1)); + { + dest_edges(1) = 0.; // so src_bins smaller than first actual dest_bin are included in first bin + } + + dest_edges(SzFreq_Nbins(SzFreqMethod) + 1) = 99999.; // dest_edges(SzFreq_Nbins(SzFreqMethod)) + (dest_edges(SzFreq_Nbins(SzFreqMethod)) - dest_edges(SzFreq_Nbins(SzFreqMethod) - 1)); dvar_vector dest_comp(1,SzFreq_Nbins(SzFreqMethod)); // destination composition. one-sex at a time // convert dest_edges to cm @@ -1067,14 +746,11 @@ FUNCTION void Get_expected_values(const int y, const int t); break; } } // end calc of new dest_edges - if (do_once == 1) echoinput << "sex: " << gg << " dest_edges_in_cm "< dest_edges[i]) overlap_low = src_edges[j]; if(src_edges[j + 1] < dest_edges[i + 1]) overlap_high = src_edges[j + 1]; - - if (overlap_low < overlap_high) + echoinput<<"rebin: dest: "<= src_edges[j_start + 1]); +// echoinput<<"jstart-after-omit "< d_low) overlap_low = s_low; + + dvariable overlap_high = d_high; + if (s_high < d_high) overlap_high = s_high; +// echoinput<<"rebin: dest: "< 1e-8) { + dest_counts[i] += src_counts[j] * (overlap_width / src_bin_width); +// echoinput<<"add src to dest: "< Date: Wed, 2 Sep 2026 10:04:46 -0700 Subject: [PATCH 6/8] fix issue affecting simple_with_DM_sizefreq --- SS_biofxn.tpl | 8 +------- SS_expval.tpl | 10 +++------- SS_param.tpl | 1 - SS_prelim.tpl | 1 - 4 files changed, 4 insertions(+), 16 deletions(-) diff --git a/SS_biofxn.tpl b/SS_biofxn.tpl index e0217d73..ecf9aafe 100644 --- a/SS_biofxn.tpl +++ b/SS_biofxn.tpl @@ -1532,20 +1532,14 @@ FUNCTION void get_wtlen() wt_len(s, gp).shift(1); } - // SS_Label_Info_19.2.3 #calculate first diff of wt_len for use in generalized size comp bin calculations - // this may no longer be needed with the rebin method if (gg == gender) { wt_len2_sq(s, GPat) = elem_prod(wt_len2(s, GPat), wt_len2(s, GPat)); - wt_len_fd(s, GPat) = first_difference(wt_len_low(s, GPat)); - if (gender == 2) - wt_len_fd(s, GPat, nlength) = wt_len_fd(s, GPat, nlength - 1); #ifdef DO_ONCE if (do_once == 1) echoinput << "wtlen2 " << endl << wt_len2 << endl - << "wtlen2^2 " << wt_len2_sq << endl - << "wtlen2:firstdiff " << wt_len_fd << endl; + << "wtlen2^2 " << wt_len2_sq << endl; #endif } } diff --git a/SS_expval.tpl b/SS_expval.tpl index 46be4329..ba6493ef 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -672,12 +672,8 @@ FUNCTION void Get_expected_values(const int y, const int t); { iobs = have_data(ALK_time, f, data_type, j); // observation index SzFreqMethod = SzFreq_obs_hdr(iobs, 6); - - SzFreqMethod_seas = nseas * (SzFreqMethod - 1) + s; // index that combines sizefreqmethod and season and used in SzFreqTrans - if (SzFreq_obs_hdr(iobs, 9) > 0) // first occurrence of this method at this time is with fleet = f - // check to see if multiple fleets using szfreq will still get included - { - for (gg = 1; gg <= gender; gg++) + SzFreqMethod_seas = nseas * (SzFreqMethod - 1) + s; // index that combines sizefreqmethod and season and used in SzFreqTrans + for (gg = 1; gg <= gender; gg++) { if (gg == 1) { @@ -729,6 +725,7 @@ FUNCTION void Get_expected_values(const int y, const int t); int parm_loc = 1 + (gg - 1) * 6; wtlen_a = wtlen_p(GPat, parm_loc); wtlen_b = wtlen_p(GPat, parm_loc + 1); +// if (do_once==1) echoinput< 0 || MGparm_seas_effects(2) > 0) // get seasonal effect on wtlen parameters { wtlen_a *= wtlen_seas(s, GPat, parm_loc); @@ -863,7 +860,6 @@ FUNCTION void Get_expected_values(const int y, const int t); } } - } // done calculating the SzFreq obs if (do_once == 1) echoinput << y << " " << f << " method: " << SzFreqMethod << " szfreq_exp_after_add_mincomp " << SzFreq_exp(iobs) << endl; } // end loop of obs for fleet = f diff --git a/SS_param.tpl b/SS_param.tpl index feab63af..e602f2f0 100644 --- a/SS_param.tpl +++ b/SS_param.tpl @@ -97,7 +97,6 @@ PARAMETER_SECTION 3darray wt_len2(1,nseas,1,N_GP,1,nlength2) // stores wt at midbin; stacked genders 3darray wt_len2_sq(1,nseas,1,N_GP,1,nlength2) // stores wt at midbin^2; stacked genders 3darray wt_len_low(1,nseas,1,N_GP,1,nlength2) // wt at lower edge of size bin - 3darray wt_len_fd(1,nseas,1,N_GP,1,nlength2-1) // first diff of wt_len_low matrix mat_len(1,N_GP,1,nlength) matrix fec_len(1,N_GP,1,nlength) // fecundity at length diff --git a/SS_prelim.tpl b/SS_prelim.tpl index 633a5f42..f520d9eb 100644 --- a/SS_prelim.tpl +++ b/SS_prelim.tpl @@ -1353,7 +1353,6 @@ get_mat_fec(); wt_len = value(wt_len); wt_len2 = value(wt_len2); - wt_len_fd = value(wt_len_fd); mat_len = value(mat_len); mat_fec_len = value(mat_fec_len); mat_age = value(mat_age); From e5569451fece219bee0c16249b83511139628fb2 Mon Sep 17 00:00:00 2001 From: Rick-Methot-NOAA Date: Wed, 2 Sep 2026 16:13:18 -0700 Subject: [PATCH 7/8] simplify treatment of small bins --- SS_expval.tpl | 7 +++---- SS_miscfxn.tpl | 51 ++++++-------------------------------------------- 2 files changed, 9 insertions(+), 49 deletions(-) diff --git a/SS_expval.tpl b/SS_expval.tpl index ba6493ef..83a02640 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -779,9 +779,9 @@ FUNCTION void Get_expected_values(const int y, const int t); if (do_once == 1) echoinput<<"pass_comp_biomass: "< dest_edges[i]) overlap_low = src_edges[j]; - if(src_edges[j + 1] < dest_edges[i + 1]) overlap_high = src_edges[j + 1]; - echoinput<<"rebin: dest: "<= src_edges[j_start + 1]); -// echoinput<<"jstart-after-omit "< 1e-8) { - dest_counts[i] += src_counts[j] * (overlap_width / src_bin_width); -// echoinput<<"add src to dest: "< Date: Fri, 4 Sep 2026 13:11:55 -0700 Subject: [PATCH 8/8] rebin now working for biomass composition; improved over legacy --- SS_expval.tpl | 235 ++++++++++++++++++++++++------------------------- SS_miscfxn.tpl | 73 ++++++++++++++- 2 files changed, 188 insertions(+), 120 deletions(-) diff --git a/SS_expval.tpl b/SS_expval.tpl index 83a02640..21c213e0 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -674,126 +674,125 @@ FUNCTION void Get_expected_values(const int y, const int t); SzFreqMethod = SzFreq_obs_hdr(iobs, 6); SzFreqMethod_seas = nseas * (SzFreqMethod - 1) + s; // index that combines sizefreqmethod and season and used in SzFreqTrans for (gg = 1; gg <= gender; gg++) - { - if (gg == 1) - { - z1 = 1; - z2 = nlength; - } // female - else - { - z1 = nlength1; - z2 = nlength2; - } // male - - // 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 - dest_edges.initialize(); - for (ibin = 1; ibin <= SzFreq_Nbins(SzFreqMethod); ibin++) { - dest_edges(ibin) = SzFreq_bins(SzFreqMethod, ibin); - } - 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 - } - - dest_edges(SzFreq_Nbins(SzFreqMethod) + 1) = 99999.; // dest_edges(SzFreq_Nbins(SzFreqMethod)) + (dest_edges(SzFreq_Nbins(SzFreqMethod)) - dest_edges(SzFreq_Nbins(SzFreqMethod) - 1)); - 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 s==1 and Gpat==1 will be used. It will be sex-specific -// wt_len_low(s, GPat)(1, nlength) = wtlen_p(GPat, 1) * pow(len_bins2(1, nlength), wtlen_p(GPat, 2)); -// W=a*L^b; so L = exp((ln(W)-ln(a))/B) - 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 (do_once==1) echoinput< 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); - break; - } - case (3): // scale is in cm - { - break; - } - case 4: // scale is in inches. Already converted to cm in readdata. Relic feature for some ancient recreational fishery data - { - break; - } - } // end calc of new dest_edges - if (do_once == 1) echoinput << "Szfreq_method: " << SzFreqMethod << " sex: " << gg << " dest_edges_in_cm "< 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); + break; + } + case (3): // scale is in cm + { + break; + } + case 4: // scale is in inches. Already converted to cm in readdata. Relic feature for some ancient recreational fishery data + { + break; + } + } // end calc of new dest_edges + if (do_once == 1) echoinput << "Szfreq_method: " << SzFreqMethod << " sex: " << gg << " dest_edges_converted_to_cm "< 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; + + // get mean weight of fish in range of src being assigned to dest + // d_low and d_high are the original weight bin boundaries (assumes the weight comp data always uses weight bin boundaries) + // s_low and s_high have weights stored in wt_len_low(s,GPat) + // src_edges_wt, dest_edges_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_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;} + + dest_counts[i] += src_counts[j] * mean_szwt * (overlap_width / src_bin_width); // add overlap fraction of biomass to dest + // echoinput<<" src: "<