diff --git a/SS_biofxn.tpl b/SS_biofxn.tpl index 00acbfc8..ecf9aafe 100644 --- a/SS_biofxn.tpl +++ b/SS_biofxn.tpl @@ -1532,19 +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 sizp comp bin calculations 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 df6fb5eb..08f17e3c 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -666,260 +666,129 @@ 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); - 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 + SzFreqMethod_seas = nseas * (SzFreqMethod - 1) + s; // index that combines sizefreqmethod and season + for (gg = 1; gg <= gender; gg++) { - if (do_once == 1 || (MG_active(3) > 0 && (timevary_MG(y, 3) > 0))) // calc matrix because wtlen parameters have changed + if (gg == 1) { - 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"<= 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); // W=a*L^b; so L = exp((ln(W)-ln(a))/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; +// echoinput<<"rebin: dest: "< 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: "< 0) echoinput << SzFreq_mincomp << " Sizefreq: add small constant to comps, per method " << endl; init_ivector SzFreq_nobs(1,SzFreq_Nmeth); !!if (SzFreq_Nmeth > 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); @@ -3188,15 +3187,9 @@ write_message(FATAL, 0); } SzFreq_Nbins3(k) = gender * SzFreq_Nbins(k); - for (s = 1; s <= nseas; s++) - { - g++; - SzFreq_Nbins_seas_g(g) = SzFreq_Nbins(k) * gender; - } } } // clang-format off - echoinput<<"here"< 0) - { - SS2out << endl - << pick_report_name(12) << endl; - SS2out << "#_NOTE: rows_are_population_length_bins;_columns_are_recipient_size_bins_according_to_the_specified_method" << endl; - for (SzFreqMethod = 1; SzFreqMethod <= SzFreq_Nmeth; SzFreqMethod++) - { - SS2out << SzFreqMethod << " gp seas len mid-len "; - if (SzFreq_scale(SzFreqMethod) == 1) - { - SS2out << " mid-kg "; - } - else if (SzFreq_scale(SzFreqMethod) == 2) - { - SS2out << " mid-lbs "; - } - else if (SzFreq_scale(SzFreqMethod) == 3) - { - SS2out << " mid-cm "; - } - else - { - SS2out << " mid-inch "; - } - SS2out << SzFreq_bins1(SzFreqMethod); - if (gender == 2) - SS2out << SzFreq_bins1(SzFreqMethod); - SS2out << endl - << SzFreqMethod << " gp seas len mid-len metric " << SzFreq_bins(SzFreqMethod) << endl; - ; - for (gp = 1; gp <= N_GP; gp++) - for (s = 1; s <= nseas; s++) - { - SzFreqMethod_seas = nseas * (SzFreqMethod - 1) + s; // index that combines sizefreqmethod and season and used in SzFreqTrans - for (z = 1; z <= nlength2; z++) - { - SS2out << SzFreqMethod << " " << gp << " " << s << " " << len_bins2(z) << " " << len_bins_m2(z) << " "; - if (SzFreq_scale(SzFreqMethod) == 1) - { - SS2out << wt_len2(s, gp, z) << " "; - } - else if (SzFreq_scale(SzFreqMethod) == 2) - { - SS2out << wt_len2(s, gp, z) / 0.4536 << " "; - } - else if (SzFreq_scale(SzFreqMethod) == 3) - { - SS2out << len_bins_m2(z) << " "; - } - else - { - SS2out << len_bins_m2(z) / 2.54 << " "; - } - for (j = 1; j <= gender * SzFreq_Nbins(SzFreqMethod); j++) - { - SS2out << SzFreqTrans(SzFreqMethod_seas, z, j) << " "; - if (SzFreqTrans(SzFreqMethod_seas, z, j) < 0.0) - { - warnstream << "Bin widths narrower than pop len bins caused negative allocation in sizefreq method:"; - warnstream << " method, season, size, bin: " << SzFreqMethod << " " << s << " " << len_bins2(z) << " " << j; - write_message (FATAL, 0); // EXIT! - } - } - SS2out << endl; - } - } - } - } + // deprecated table when method transitioned to rebin // REPORT_KEYWORD 13 MOVEMENT - if (pick_report_use(13) == "Y" && do_migration > 0) { SS2out << endl @@ -2705,7 +2634,7 @@ FUNCTION void write_bigoutput() dvector tempvec_l(1, SzFreq_exp(iobs).size()); tempvec_l = value(SzFreq_exp(iobs)); more_comp_info = process_comps(gender, gg, SzFreq_bins(sz_method), SzFreq_means(sz_method), sz_tails, SzFreq_obs(iobs), tempvec_l); - Nsamp_DM = SzFreq_sampleN(iobs); // Will remain this if not used; there is no "adjusted" sample size for sizwfreq + Nsamp_DM = SzFreq_sampleN(iobs); // Will remain this if not used; there is no "adjusted" sample size for sizefreq if (Comp_Err_Sz(sz_method) == 1) // Dirichlet #1 { dirichlet_Parm = mfexp(selparm(Comp_Err_parmloc(Comp_Err_Sz2(sz_method),1))); // Thorson's theta from eq 10 @@ -3624,7 +3553,7 @@ FUNCTION void write_bigoutput() for (f1 = 1; f1 <= N_pred; f1++) { f = predator(f1); - SS2out << fleetname(f) << "_M2 comsume_Bio consume_Num"; + SS2out << fleetname(f) << "_M2 consume_Bio consume_Num"; } SS2out << endl; for (y = styr - 2; y <= YrMax; y++)