Skip to content

Commit 2fa8db0

Browse files
Rick-Methot-NOAARick-Methot-NOAA
committed
rebin now working for biomass composition; improved over legacy
1 parent e556945 commit 2fa8db0

2 files changed

Lines changed: 188 additions & 120 deletions

File tree

‎SS_expval.tpl‎

Lines changed: 117 additions & 118 deletions
Original file line numberDiff line numberDiff line change
@@ -674,126 +674,125 @@ FUNCTION void Get_expected_values(const int y, const int t);
674674
SzFreqMethod = SzFreq_obs_hdr(iobs, 6);
675675
SzFreqMethod_seas = nseas * (SzFreqMethod - 1) + s; // index that combines sizefreqmethod and season and used in SzFreqTrans
676676
for (gg = 1; gg <= gender; gg++)
677-
{
678-
if (gg == 1)
679-
{
680-
z1 = 1;
681-
z2 = nlength;
682-
} // female
683-
else
684-
{
685-
z1 = nlength1;
686-
z2 = nlength2;
687-
} // male
688-
689-
// NOTE: wt_len_low is calculated separately for each growth pattern (GPat)
690-
// but the code below still just uses GPat=1 for calculation of the sizefreq transition matrix
691-
dvector src_edges(1, nlength + 1); // +1 to create a dummy bin holding upper edge of the actual last bin
692-
src_edges(1, nlength) = len_bins(1, nlength);
693-
src_edges(nlength + 1) = src_edges(nlength) + (src_edges(nlength) - src_edges(nlength - 1));
694-
dvar_vector dest_edges(1,SzFreq_Nbins(SzFreqMethod) + 1); // +1 to create a dummy bin holding upper edge of the actual last bin
695-
dest_edges.initialize();
696-
for (ibin = 1; ibin <= SzFreq_Nbins(SzFreqMethod); ibin++) {
697-
dest_edges(ibin) = SzFreq_bins(SzFreqMethod, ibin);
698-
}
699-
if(SzFreq_Omit_Small(SzFreqMethod) < 0.0)
700-
{
701-
dest_edges(1) = 0.; // so src_bins smaller than first actual dest_bin are included in first bin
702-
}
703-
704-
dest_edges(SzFreq_Nbins(SzFreqMethod) + 1) = 99999.; // dest_edges(SzFreq_Nbins(SzFreqMethod)) + (dest_edges(SzFreq_Nbins(SzFreqMethod)) - dest_edges(SzFreq_Nbins(SzFreqMethod) - 1));
705-
dvar_vector dest_comp(1,SzFreq_Nbins(SzFreqMethod)); // destination composition. one-sex at a time
706-
707-
// convert dest_edges to cm
708-
// as in the legacy method, this conversion is sex and season specific, but only uses Gpat = 1
709-
switch (SzFreq_scale(SzFreqMethod)) // biomass vs. numbers are accumulated in the bins
710-
{
711-
case (1): // scale is kg
712-
{
713-
// fall thru to calcs in the lb case
714-
}
715-
case (2): // scale is lb, or fall thru from kg
716-
{
717-
// bin boundary in lb have already been converted to kg
718-
// 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
719-
// wt_len_low(s, GPat)(1, nlength) = wtlen_p(GPat, 1) * pow(len_bins2(1, nlength), wtlen_p(GPat, 2));
720-
// W=a*L^b; so L = exp((ln(W)-ln(a))/B)
721-
dvariable wtlen_a;
722-
dvariable wtlen_b;
723-
GPat = 1; // this means that any (unusual) use of multiple Gpatterns will still use biology of GPat 1 for szfreq purposes
724-
// female parms are at 1, 2; males at 7,8
725-
int parm_loc = 1 + (gg - 1) * 6;
726-
wtlen_a = wtlen_p(GPat, parm_loc);
727-
wtlen_b = wtlen_p(GPat, parm_loc + 1);
728-
// if (do_once==1) echoinput<<gg<<" parmloc "<<parm_loc<<" a, b "<<wtlen_a<<" " << wtlen_b<<endl;
729-
if (MGparm_seas_effects(1) > 0 || MGparm_seas_effects(2) > 0) // get seasonal effect on wtlen parameters
730-
{
731-
wtlen_a *= wtlen_seas(s, GPat, parm_loc);
732-
wtlen_b *= wtlen_seas(s, GPat, parm_loc + 1);
733-
}
734-
dest_edges = exp((log(dest_edges)-log(wtlen_a))/wtlen_b);
735-
break;
736-
}
737-
case (3): // scale is in cm
738-
{
739-
break;
740-
}
741-
case 4: // scale is in inches. Already converted to cm in readdata. Relic feature for some ancient recreational fishery data
742-
{
743-
break;
744-
}
745-
} // end calc of new dest_edges
746-
if (do_once == 1) echoinput << "Szfreq_method: " << SzFreqMethod << " sex: " << gg << " dest_edges_in_cm "<<dest_edges<<endl;
747-
748-
dvar_vector pass_comp(z1, z2); // this gets males or females range according to z1, z2
749-
pass_comp.initialize();
750-
// if (do_once == 1) echoinput << "retain/discard obs type: "<<SzFreq_obs_hdr(iobs, 5)<<endl;
751-
switch (SzFreq_obs_hdr(iobs, 5)) // discard/retained partition
752-
{
753-
case (0):
754-
{
755-
pass_comp = exp_l_temp(z1, z2);
756-
// SzFreq_exp(iobs) = trans(SzFreqTrans(SzFreqMethod_seas)) * exp_l_temp;
757-
break;
758-
}
759-
case (1):
760-
{
761-
pass_comp = exp_l_temp(z1, z2) - exp_l_temp_ret(z1, z2);
762-
// SzFreq_exp(iobs) = trans(SzFreqTrans(SzFreqMethod_seas)) * (exp_l_temp - exp_l_temp_ret);
763-
break;
764-
}
765-
case (2):
766-
{
767-
pass_comp = exp_l_temp_ret(z1, z2);
768-
// SzFreq_exp(iobs) = trans(SzFreqTrans(SzFreqMethod_seas)) * exp_l_temp_ret;
769-
break;
770-
}
771-
}
772-
int use_GP;
773-
use_GP = 1;
774-
if (gg == 2) {use_GP += N_GP;} // because GPat is nested in sex
775-
pass_comp.shift(1); // change index to 1
776-
if (SzFreq_units(SzFreqMethod) == 1) // biomass is accumulated in the bins
777-
{
778-
pass_comp = elem_prod(pass_comp, wt_len(s, use_GP)); // where wtlen has been calculated using mid length of length bins
779-
if (do_once == 1) echoinput<<"pass_comp_biomass: "<<pass_comp<<endl;
780-
// wt_len being referenced here by gg (sex), but actually is stored by gp which is sex and GPat. So, only the first Gpat gets used
781-
}
782-
// if (do_once == 1) echoinput<<"ready to call rebin, pass_comp: "<<pass_comp<<endl<<" sum_pass: " << sum(pass_comp)<<endl;
783-
dest_comp = rebin(SzFreq_Omit_Small(SzFreqMethod), src_edges, pass_comp, dest_edges);
784-
// if (do_once == 1) echoinput << "dest_comp: "<<dest_comp<<endl<<" sum_dest: " << sum(dest_comp)<<endl;
677+
{
678+
if (gg == 1)
679+
{
680+
z1 = 1;
681+
z2 = nlength;
682+
} // female
683+
else
684+
{
685+
z1 = nlength1;
686+
z2 = nlength2;
687+
} // male
688+
689+
// NOTE: wt_len_low is calculated separately for each growth pattern (GPat)
690+
// but the code below still just uses GPat=1 for calculation of the sizefreq transition matrix
691+
dvector src_edges(1, nlength + 1); // +1 to create a dummy bin holding upper edge of the actual last bin
692+
src_edges(1, nlength) = len_bins(1, nlength);
693+
src_edges(nlength + 1) = src_edges(nlength) + (src_edges(nlength) - src_edges(nlength - 1));
694+
dvar_vector dest_edges(1,SzFreq_Nbins(SzFreqMethod) + 1); // +1 to create a dummy bin holding upper edge of the actual last bin
695+
dvar_vector dest_edges_wt(1,SzFreq_Nbins(SzFreqMethod) + 1);
696+
dest_edges.initialize();
697+
for (ibin = 1; ibin <= SzFreq_Nbins(SzFreqMethod); ibin++) {
698+
dest_edges(ibin) = SzFreq_bins(SzFreqMethod, ibin);
699+
dest_edges_wt(ibin) = dest_edges(ibin); // save in raw units
700+
}
701+
dest_edges(SzFreq_Nbins(SzFreqMethod) + 1) = 99999.;
702+
dest_edges_wt(ibin) = dest_edges_wt(ibin - 1) + (dest_edges_wt(ibin - 1) - dest_edges_wt(ibin - 2));
703+
if(SzFreq_Omit_Small(SzFreqMethod) < 0.0)
704+
{
705+
dest_edges(1) = 0.; // so src_bins smaller than first actual dest_bin are included in first bin
706+
}
785707
786-
if( gg == 1)
787-
{
788-
SzFreq_exp(iobs)(1, SzFreq_Nbins(SzFreqMethod)) = dest_comp;
789-
}
790-
else
791-
{
792-
dest_comp.shift(SzFreq_Nbins(SzFreqMethod) + 1);
793-
SzFreq_exp(iobs)(SzFreq_Nbins(SzFreqMethod) + 1, 2 * SzFreq_Nbins(SzFreqMethod)) = dest_comp;
794-
}
708+
dvar_vector dest_comp(1,SzFreq_Nbins(SzFreqMethod)); // destination composition. one-sex at a time
709+
// convert dest_edges to cm
710+
// as in the legacy method, this conversion is sex and season specific, but only uses Gpat = 1
711+
switch (SzFreq_scale(SzFreqMethod)) // biomass vs. numbers are accumulated in the bins
712+
{
713+
case (1): // scale is kg
714+
{
715+
// fall thru to calcs in the lb case
716+
}
717+
case (2): // scale is lb, or fall thru from kg
718+
{
719+
// bin boundary in lb have already been converted to kg
720+
// 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
721+
// W=a*L^b; so L = exp((ln(W)-ln(a))/B)
722+
dvariable wtlen_a;
723+
dvariable wtlen_b;
724+
GPat = 1; // this means that any (unusual) use of multiple Gpatterns will still use biology of GPat 1 for szfreq purposes
725+
// female parms are at 1, 2; males at 7,8
726+
int parm_loc = 1 + (gg - 1) * 6;
727+
wtlen_a = wtlen_p(GPat, parm_loc);
728+
wtlen_b = wtlen_p(GPat, parm_loc + 1);
729+
if (MGparm_seas_effects(1) > 0 || MGparm_seas_effects(2) > 0) // get seasonal effect on wtlen parameters
730+
{
731+
wtlen_a *= wtlen_seas(s, GPat, parm_loc);
732+
wtlen_b *= wtlen_seas(s, GPat, parm_loc + 1);
733+
}
734+
dest_edges = exp((log(dest_edges)-log(wtlen_a))/wtlen_b);
735+
break;
736+
}
737+
case (3): // scale is in cm
738+
{
739+
break;
740+
}
741+
case 4: // scale is in inches. Already converted to cm in readdata. Relic feature for some ancient recreational fishery data
742+
{
743+
break;
744+
}
745+
} // end calc of new dest_edges
746+
if (do_once == 1) echoinput << "Szfreq_method: " << SzFreqMethod << " sex: " << gg << " dest_edges_converted_to_cm "<<dest_edges<<endl;
795747
796-
} // end gender loop
748+
dvar_vector pass_comp(z1, z2); // this gets males or females range according to z1, z2
749+
pass_comp.initialize();
750+
switch (SzFreq_obs_hdr(iobs, 5)) // discard/retained partition
751+
{
752+
case (0):
753+
{
754+
pass_comp = exp_l_temp(z1, z2);
755+
// SzFreq_exp(iobs) = trans(SzFreqTrans(SzFreqMethod_seas)) * exp_l_temp;
756+
break;
757+
}
758+
case (1):
759+
{
760+
pass_comp = exp_l_temp(z1, z2) - exp_l_temp_ret(z1, z2);
761+
// SzFreq_exp(iobs) = trans(SzFreqTrans(SzFreqMethod_seas)) * (exp_l_temp - exp_l_temp_ret);
762+
break;
763+
}
764+
case (2):
765+
{
766+
pass_comp = exp_l_temp_ret(z1, z2);
767+
// SzFreq_exp(iobs) = trans(SzFreqTrans(SzFreqMethod_seas)) * exp_l_temp_ret;
768+
break;
769+
}
770+
}
771+
int use_GP;
772+
use_GP = 1;
773+
if (gg == 2) {use_GP += N_GP;} // because GPat is nested in sex; only GPat 1 bio gets used for szfreq data
774+
pass_comp.shift(1); // change index to 1
775+
if (SzFreq_units(SzFreqMethod) == 1) // biomass is accumulated in the bins
776+
{
777+
dvar_vector pass_src_wt(z1, z2+1);
778+
pass_src_wt(z1,z2) = wt_len_low(s, 1)(z1,z2); // just one sex's range is passed
779+
pass_src_wt.shift(1); // change index to 1
780+
dest_comp = rebin_bio(src_edges, pass_comp, dest_edges, pass_src_wt, dest_edges_wt);
781+
}
782+
else
783+
{
784+
dest_comp = rebin(src_edges, pass_comp, dest_edges);
785+
}
786+
if( gg == 1)
787+
{
788+
SzFreq_exp(iobs)(1, SzFreq_Nbins(SzFreqMethod)) = dest_comp;
789+
}
790+
else
791+
{
792+
dest_comp.shift(SzFreq_Nbins(SzFreqMethod) + 1);
793+
SzFreq_exp(iobs)(SzFreq_Nbins(SzFreqMethod) + 1, 2 * SzFreq_Nbins(SzFreqMethod)) = dest_comp;
794+
}
795+
} // end gender loop
797796
798797
if (gender == 2)
799798
{

‎SS_miscfxn.tpl‎

Lines changed: 71 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@ FUNCTION dvariable Comp_logL_Dirichlet(const double& Nsamp, const dvariable& dir
122122
}
123123

124124

125-
FUNCTION dvar_vector rebin(const int omit_small, const dvector& src_edges, const dvar_vector& src_counts, const dvar_vector& dest_edges)
125+
FUNCTION dvar_vector rebin(const dvector& src_edges, const dvar_vector& src_counts, const dvar_vector& dest_edges)
126126
{
127127
/*
128128
This implementation takes two vectors representing the boundaries (edges) of the source
@@ -171,9 +171,78 @@ FUNCTION dvar_vector rebin(const int omit_small, const dvector& src_edges, const
171171
dvariable overlap_width = overlap_high - overlap_low;
172172
dvariable src_bin_width = s_high - s_low;
173173
dest_counts[i] += src_counts[j] * (overlap_width / src_bin_width);
174-
// echoinput<<"add src to dest: "<<dest_edges[i]<<" "<<src_edges[j]<<" result: "<<(overlap_width / src_bin_width)<<endl;
174+
// echoinput<<"add src to dest: "<<dest_edges[i]<<" "<<src_edges[j]<<" result: "<<(overlap_width / src_bin_width)<<endl;
175175
}
176176
}
177177
}
178178
return (dest_counts);
179+
}
180+
181+
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 )
182+
{
183+
/*
184+
This modification of rebin is used when biomass is accumulated into the bins
185+
for example, with catch weight composition
186+
it takes into account the fact that fish in the lower portion of a length bin have less body weight than fish in the opper portion of the length bin
187+
NOTE: need to undo the conversion of numbers to biomass in SS_expval. It needs to occur here.
188+
legacy szfreq method used this approach:
189+
temp = (wt_len_low(s, 1, z + 1) - topbin) / wt_len_fd(s, 1, z); // frac in pop bin above (data bin +1)
190+
temp1 = wt_len_low(s, 1, z) + (1. - temp * 0.5) * wt_len_fd(s, 1, z); // approx body wt for these fish
191+
temp2 = wt_len_low(s, 1, z) + (1. - temp) * 0.5 * wt_len_fd(s, 1, z); // approx body wt for fish below
192+
SzFreqTrans(SzFreqMethod_seas, z, ibinsave + 1) = temp * temp1;
193+
SzFreqTrans(SzFreqMethod_seas, z, ibinsave) = (1. - temp) * temp2;
194+
new approach has access to the body wt at boundaries of size range of fish getting rebinned, so will use that to get more exact body weights
195+
*/
196+
dvar_vector dest_counts(1, dest_edges.size() - 1); // sized to leave off the topbin boundary
197+
dest_counts.initialize();
198+
int j_start = 1;
199+
for (int i = 1; i <= dest_counts.size(); i++) { // loop the destination bins
200+
dvariable d_low = dest_edges[i];
201+
dvariable d_high = dest_edges[i + 1];
202+
// Advance j_start if the source bin is entirely below the current destination bin.
203+
// Because d_low increases with 'i', j_start only ever moves forward.
204+
while (j_start <= src_counts.size() && src_edges[j_start + 1] <= d_low) {
205+
j_start++;
206+
}
207+
// Iterate through source bins starting from j_start, but stop as soon
208+
// as the source bin is completely above the current destination bin.
209+
for (int j = j_start; j <= src_counts.size() && src_edges[j] < d_high; j++) {
210+
dvariable s_low = src_edges[j];
211+
dvariable s_high = src_edges[j + 1];
212+
// Calculate the overlap bounds
213+
dvariable overlap_low = d_low;
214+
if (s_low > d_low) overlap_low = s_low;
215+
216+
dvariable overlap_high = d_high;
217+
if (s_high < d_high) overlap_high = s_high;
218+
219+
// If there is valid overlap, distribute the counts
220+
if (overlap_low < overlap_high) {
221+
dvariable overlap_width = overlap_high - overlap_low;
222+
dvariable src_bin_width = s_high - s_low;
223+
224+
// get mean weight of fish in range of src being assigned to dest
225+
// d_low and d_high are the original weight bin boundaries (assumes the weight comp data always uses weight bin boundaries)
226+
// s_low and s_high have weights stored in wt_len_low(s,GPat)
227+
// src_edges_wt, dest_edges_wt
228+
dvariable mean_szwt; // body weight of fish to be allocated from src to dest
229+
if (s_low >= d_low)
230+
{
231+
if (s_high <= d_high) // src entirely in dest range range
232+
{mean_szwt = (src_edges_wt[j] + src_edges_wt[j+1]) * 0.5;}
233+
else // s_high > d_high; calc mean size from s_low and d_high
234+
{mean_szwt = (src_edges_wt[j] + dest_edges_wt[i+1]) * 0.5;}
235+
}
236+
else if (s_high <= d_high ) // overlap with s_low < d_low; calc the mean size from d_low and s_high
237+
{mean_szwt = (dest_edges_wt[i] + src_edges_wt[j+1]) * 0.5;}
238+
else // s_high > d_high, so dest bin is internal to src_bin
239+
{mean_szwt = (dest_edges_wt[i] + dest_edges_wt[i+1]) * 0.5;}
240+
241+
dest_counts[i] += src_counts[j] * mean_szwt * (overlap_width / src_bin_width); // add overlap fraction of biomass to dest
242+
// echoinput<<" src: "<<j<<" lo_hi "<<s_low << " "<<s_high<<" dest: "<<i<<" lo_hi "<<d_low<<" "<<d_high<<" overlap: "<<overlap_low <<" "<<overlap_high
243+
// <<" fraction: "<<(overlap_width / src_bin_width)<<" mean: "<<mean_szwt<<endl;
244+
} // end having overlap to be allocated
245+
} // end j loop of source bins
246+
}
247+
return (dest_counts);
179248
}

0 commit comments

Comments
 (0)