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/3] 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/3] 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/3] 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<<" "<