first commit of rebin

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
         }
       }

SS_miscfxn.tpl

@@ -120,3 +120,129 @@ 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 dvar_vector rebin(const dvector& src_edges, const dvar_vector& src_counts, const dvar_vector& 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).
+    src_counts need not be counts; works for real number of fish, or for biomass of fish.
+    However, if src_counts are in biomass units, this rebin method does not account for fact that fish in lower portion of a src bin
+       will weigh less than fish in upper portion.
+    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();
+    int j_start = 1;
+    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;
+//            echoinput<<"rebin: dest: "<<i<<" "<<dest_edges[i]<<" src: "<<j<<" "<<src_edges[j]<<" overlap_lo "<<overlap_low <<" overlap_hi "<<overlap_high <<  endl;
+
+            // 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;
+                dest_counts[i] += src_counts[j] * (overlap_width / src_bin_width);
+//                echoinput<<"add src to dest: "<<dest_edges[i]<<" "<<src_edges[j]<<" result: "<<(overlap_width / src_bin_width)<<endl;
+            }
+        }
+      }
+    return (dest_counts);
+  }
+
+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 )
+  {
+    /*
+    This modification of rebin is used when biomass is accumulated into the bins
+    for example, with catch weight composition
+    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
+    NOTE:  need to undo the conversion of numbers to biomass in SS_expval.  It needs to occur here.
+    legacy szfreq method used this approach:
+      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;
+    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
+    */
+    dvar_vector dest_counts(1, dest_edges.size() - 1);  // sized to leave off the topbin boundary
+    dest_counts.initialize();
+    int j_start = 1;
+    for (int i = 1; i <= dest_counts.size(); i++) {  // loop the destination bins
+        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;
+
+            // 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: "<<j<<" lo_hi "<<s_low << " "<<s_high<<" dest: "<<i<<" lo_hi "<<d_low<<" "<<d_high<<" overlap: "<<overlap_low <<" "<<overlap_high
+            //  <<" fraction: "<<(overlap_width / src_bin_width)<<" mean: "<<mean_szwt<<endl;
+          }  // end having overlap to be allocated
+        }  // end j loop of source bins
+      }
+    return (dest_counts);
+  }

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

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);

SS_readdata_330.tpl

@@ -3136,7 +3136,7 @@
   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);
@@ -3196,7 +3196,6 @@
     }
   }
   // clang-format off
-  echoinput<<"here"<<endl;
  END_CALCS
 
  !!echoinput<<"bins "<<SzFreq_Nbins<<endl;
@@ -3266,9 +3265,14 @@
           SzFreq_means(k, z) = SzFreq_means(k, z - 1) + (SzFreq_bins2(k, z) - SzFreq_bins2(k, z - 1));
         }
         if (gender == 2)
+        {
           SzFreq_means(k, z + SzFreq_Nbins(k)) = SzFreq_means(k, z);
+        }
       }
+//      SzFreq_bins2(k, SzFreq_Nbins(k)) = 99999.;
+//      if (gender == 2) SzFreq_bins2(k,SzFreq_Nbins(k) + SzFreq_Nbins(k) + 1) = 99999.;
       echoinput << "Processed_SizeFreqMethod_bins for method: " << k << endl
+                << "bins2: " << SzFreq_bins2(k) << endl
                 << "low: " << SzFreq_bins(k) << endl
                 << "mean: " << SzFreq_means(k) << endl;
     }

SS_write_report.tpl

@@ -3624,7 +3624,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++)