Per #2339, refine logic for -job aggregate -line_type SEEPS.

src/libcode/vx_seeps/seeps.cc

@@ -41,6 +41,8 @@ static const char *var_name_elv       = "elv";
 static const char *dim_name_lat       = "latitude";
 static const char *dim_name_lon       = "longitude";
 
+double weighted_average(double, double, double, double);
+
 ////////////////////////////////////////////////////////////////////////
 
 SeepsClimo *get_seeps_climo() {
@@ -546,6 +548,65 @@ void SeepsAggScore::clear() {
 ////////////////////////////////////////////////////////////////////////
 
 
+SeepsAggScore & SeepsAggScore::operator+=(const SeepsAggScore &c) {
+
+   // Check for degenerate case
+   if(n_obs == 0 && c.n_obs == 0) return(*this);
+
+   // Compute weights
+   double w1 = (double)   n_obs / (n_obs + c.n_obs);
+   double w2 = (double) c.n_obs / (n_obs + c.n_obs);
+
+   // Increment number of obs
+   n_obs += c.n_obs;
+
+   // Increment counts
+   c12 += c.c12;
+   c13 += c.c13;
+   c21 += c.c21;
+   c23 += c.c23;
+   c31 += c.c31;
+   c32 += c.c32;
+
+   // Compute weighted averages
+   s12 = weighted_average(s12, w1, c.s12, w2);
+   s13 = weighted_average(s13, w1, c.s13, w2);
+   s21 = weighted_average(s21, w1, c.s21, w2);
+   s23 = weighted_average(s23, w1, c.s23, w2);
+   s31 = weighted_average(s31, w1, c.s31, w2);
+   s32 = weighted_average(s32, w1, c.s32, w2);
+
+   pv1 = weighted_average(pv1, w1, c.pv1, w2);
+   pv2 = weighted_average(pv2, w1, c.pv2, w2);
+   pv3 = weighted_average(pv3, w1, c.pv3, w2);
+
+   pf1 = weighted_average(pf1, w1, c.pf1, w2);
+   pf2 = weighted_average(pf2, w1, c.pf2, w2);
+   pf3 = weighted_average(pf3, w1, c.pf3, w2);
+
+   mean_fcst = weighted_average(mean_fcst, w1, c.mean_fcst, w2);
+   mean_obs  = weighted_average(mean_obs,  w1, c.mean_obs,  w2);
+
+   score          = weighted_average(score,          w1, c.score,          w2);
+   weighted_score = weighted_average(weighted_score, w1, c.weighted_score, w2);
+
+   return(*this);
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+
+double weighted_average(double v1, double w1, double v2, double w2) {
+   return(is_bad_data(v1) || is_bad_data(v2) ?
+          bad_data_double :
+          v1 * w1 + v2 * w2);
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+
 SeepsClimoGrid::SeepsClimoGrid(int month, int hour) : month{month}, hour{hour}
 {
 

src/libcode/vx_seeps/seeps.h

@@ -74,7 +74,8 @@ struct SeepsScore { // For SEEPS_MPR
 struct SeepsAggScore {  // For SEEPS
 
    void clear();
-
+   SeepsAggScore & operator+=(const SeepsAggScore &);
+
    int   n_obs;
    int   c12;
    int   c13;

src/tools/core/stat_analysis/aggr_stat_line.cc

@@ -35,6 +35,8 @@
 //   015    01/24/20  Halley Gotway   Add aggregate RPS lines.
 //   016    04/12/21  Halley Gotway   MET #1735 Support multiple
 //                    -out_thresh and -out_line_type options.
+//   017    11/10/22  Halley Gotway   MET #2339 Add SEEPS and SEEPS_MPR
+//                                      line types.
 //
 ////////////////////////////////////////////////////////////////////////
 
@@ -3356,6 +3358,89 @@ void aggr_ssvar_lines(LineDataFile &f, STATAnalysisJob &job,
 
 ////////////////////////////////////////////////////////////////////////
 
+void aggr_seeps_lines(LineDataFile &f, STATAnalysisJob &job,
+                      map<ConcatString, AggrSEEPSInfo> &m,
+                      int &n_in, int &n_out) {
+   STATLine line;
+   AggrSEEPSInfo aggr;
+   SeepsAggScore cur;
+   ConcatString key, fcst_var, obs_var;
+
+   //
+   // Process the STAT lines
+   //
+   while(f >> line) {
+
+      if(line.is_header()) continue;
+
+      n_in++;
+
+      if(job.is_keeper(line)) {
+
+         job.dump_stat_line(line);
+
+         if(line.type() != stat_seeps) {
+            mlog << Error << "\naggr_seeps_lines() -> "
+                 << "should only encounter SEEPS line types.\n"
+                 << "ERROR occurred on STAT line:\n" << line << "\n\n";
+            throw(1);
+         }
+
+         //
+         // Only aggregate consistent variable names
+         //
+         if(fcst_var.empty()) fcst_var = line.fcst_var();
+         if(obs_var.empty())  obs_var  = line.obs_var();
+
+         if(fcst_var != line.fcst_var() ||
+            obs_var  != line.obs_var()) {
+            mlog << Error << "\naggr_seeps_lines() -> "
+                 << "both the forecast and observation variable types must "
+                 << "remain constant.  Try setting \"-fcst_var\" and/or "
+                 << "\"-obs_var\".\n"
+                 << "ERROR occurred on STAT line:\n" << line << "\n\n";
+            throw(1);
+         }
+
+         //
+         // Parse the current SSVAR line
+         //
+         parse_seeps_line(line, cur);
+
+         //
+         // Build the map key for the current line
+         //
+         key = job.get_case_info(line);
+
+         //
+         // Add a new map entry, if necessary
+         //
+         if(m.count(key) == 0) {
+            aggr.agg_score = cur;
+            aggr.hdr.clear();
+            m[key] = aggr;
+         }
+         //
+         // Increment counts in the existing map entry
+         //
+         else {
+            m[key].agg_score += cur;
+         }
+
+         //
+         // Keep track of the unique header column entries
+         //
+         m[key].hdr.add(line);
+
+         n_out++;
+      }
+   } // end while
+
+   return;
+}
+
+////////////////////////////////////////////////////////////////////////
+
 void aggr_time_series_lines(LineDataFile &f, STATAnalysisJob &job,
                             map<ConcatString, AggrTimeSeriesInfo> &m,
                             int &n_in, int &n_out) {

src/tools/core/stat_analysis/aggr_stat_line.h

@@ -192,6 +192,11 @@ struct AggrSSVARInfo {
    std::map<ConcatString, SSVARInfo, ssvar_bin_cmp> ssvar_bins;
 };
 
+struct AggrSEEPSInfo {
+   StatHdrInfo hdr;
+   SeepsAggScore agg_score;
+};
+
 struct AggrTimeSeriesInfo {
    StatHdrInfo hdr;
    ConcatString fcst_var, obs_var;
@@ -290,6 +295,11 @@ extern void aggr_ssvar_lines(
                std::map<ConcatString, AggrSSVARInfo> &,
                int &, int &);
 
+extern void aggr_seeps_lines(
+               LineDataFile &, STATAnalysisJob &,
+               std::map<ConcatString, AggrSEEPSInfo> &,
+               int &, int &);
+
 extern void aggr_time_series_lines(
                LineDataFile &, STATAnalysisJob &,
                std::map<ConcatString, AggrTimeSeriesInfo> &,

src/tools/core/stat_analysis/stat_analysis_job.cc

@@ -43,6 +43,8 @@
 //   021    04/12/21  Halley Gotway   MET #1735 Support multiple 
 //                    -out_thresh and -out_line_type options.
 //   022    05/28/21  Halley Gotway   Add MCTS HSS_EC output.
+//   023    11/10/22  Halley Gotway   MET #2339 Add SEEPS and SEEPS_MPR
+//                                      line types.
 //
 ////////////////////////////////////////////////////////////////////////
 
@@ -421,6 +423,7 @@ void do_job_aggr(const ConcatString &jobstring, LineDataFile &f,
    map<ConcatString, AggrENSInfo>   ens_map;
    map<ConcatString, AggrRPSInfo>   rps_map;
    map<ConcatString, AggrSSVARInfo> ssvar_map;
+   map<ConcatString, AggrSEEPSInfo> seeps_map;
 
    //
    // Check that the -line_type option has been supplied exactly once
@@ -449,13 +452,14 @@ void do_job_aggr(const ConcatString &jobstring, LineDataFile &f,
       lt != stat_grad   && lt != stat_ecnt   &&
       lt != stat_rps    && lt != stat_rhist  &&
       lt != stat_phist  && lt != stat_relp   &&
-      lt != stat_ssvar  && lt != stat_isc) {
+      lt != stat_ssvar  && lt != stat_isc    &&
+      lt != stat_seeps) {
       mlog << Error << "\ndo_job_aggr() -> "
            << "the \"-line_type\" option must be set to one of:\n"
            << "\tFHO, CTC, MCTC,\n"
            << "\tSL1L2, SAL1L2, VL1L2, VAL1L2,\n"
            << "\tPCT, NBRCTC, NBRCNT, GRAD, ISC,\n"
-           << "\tECNT, RPS, RHIST, PHIST, RELP, SSVAR\n\n";
+           << "\tECNT, RPS, RHIST, PHIST, RELP, SSVAR, SEEPS\n\n";
       throw(1);
    }
 
@@ -571,6 +575,14 @@ void do_job_aggr(const ConcatString &jobstring, LineDataFile &f,
       write_job_aggr_ssvar(job, lt, ssvar_map, out_at);
    }
 
+   //
+   // Sum the SEEPS line types
+   //
+   else if(lt == stat_seeps) {
+      aggr_seeps_lines(f, job, seeps_map, n_in, n_out);
+      write_job_aggr_seeps(job, lt, seeps_map, out_at);
+   }
+
    //
    // Check for no matching STAT lines
    //
@@ -2677,6 +2689,70 @@ void write_job_aggr_ssvar(STATAnalysisJob &job, STATLineType lt,
 
 ////////////////////////////////////////////////////////////////////////
 
+void write_job_aggr_seeps(STATAnalysisJob &job, STATLineType lt,
+                         map<ConcatString, AggrSEEPSInfo> &m,
+                         AsciiTable &at) {
+   map<ConcatString, AggrSEEPSInfo>::iterator it;
+   int n, n_row, n_col, r, c;
+   StatHdrColumns shc;
+
+   //
+   // Setup the output table
+   //
+   n_row = 1 + m.size();
+   n_col = 1 + job.by_column.n() + n_seeps_columns;
+   write_job_aggr_hdr(job, n_row, n_col, at);
+
+   //
+   // Write the rest of the header row
+   //
+   c = 1 + job.by_column.n();
+   write_header_row(seeps_columns, n_seeps_columns, 0, at, 0, c);
+
+   //
+   // Setup the output STAT file
+   //
+   job.setup_stat_file(n_row, n);
+
+   mlog << Debug(2) << "Computing output for "
+        << (int) m.size() << " case(s).\n";
+
+   //
+   // Loop through the map
+   //
+   for(it = m.begin(), r=1; it != m.end(); it++) {
+
+      //
+      // Write the output STAT header columns
+      //
+      shc = it->second.hdr.get_shc(it->first, job.by_column,
+                                   job.hdr_name, job.hdr_value, lt);
+
+      //
+      // Initialize
+      //
+      c = 0;
+
+      //
+      // SEEPS output line
+      //
+      if(job.stat_out) {
+         write_header_cols(shc, job.stat_at, job.stat_row);
+         write_seeps_cols(&it->second.agg_score, job.stat_at,
+                          job.stat_row++, n_header_columns);
+      }
+      else {
+         at.set_entry(r, c++, "SEEPS:");
+         write_case_cols(it->first, at, r, c);
+         write_seeps_cols(&it->second.agg_score, at, r++, c);
+      }
+   } // end for it
+
+   return;
+}
+
+////////////////////////////////////////////////////////////////////////
+
 void write_job_aggr_orank(STATAnalysisJob &job, STATLineType lt,
                           map<ConcatString, AggrENSInfo> &m,
                           AsciiTable &at, gsl_rng *rng_ptr) {

src/tools/core/stat_analysis/stat_analysis_job.h

@@ -118,6 +118,9 @@ extern void write_job_aggr_relp(STATAnalysisJob &, STATLineType,
 extern void write_job_aggr_ssvar(STATAnalysisJob &, STATLineType,
                std::map<ConcatString, AggrSSVARInfo> &, AsciiTable &);
 
+extern void write_job_aggr_seeps(STATAnalysisJob &, STATLineType,
+               std::map<ConcatString, AggrSEEPSInfo> &, AsciiTable &);
+
 extern void write_job_aggr_orank(STATAnalysisJob &, STATLineType,
                std::map<ConcatString, AggrENSInfo> &, AsciiTable &,
                gsl_rng *);