From aa858f3f45d0b268d6f88496759fab0e5ec1c345 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 12 Nov 2021 10:32:41 -0700
Subject: [PATCH] Per #1809, saving off the current state of development before
 proceeding with changes to simplify the matching logic.

---
 met/src/libcode/vx_tc_util/prob_gen_info.cc |  24 ++++
 met/src/libcode/vx_tc_util/prob_gen_info.h  |   8 +-
 met/src/tools/tc_utils/tc_gen/tc_gen.cc     | 128 +++++++++++++++++++-
 met/src/tools/tc_utils/tc_pairs/Makefile.am |   1 +
 4 files changed, 154 insertions(+), 7 deletions(-)

diff --git a/met/src/libcode/vx_tc_util/prob_gen_info.cc b/met/src/libcode/vx_tc_util/prob_gen_info.cc
index 863fc90bc1..e292ce897a 100644
--- a/met/src/libcode/vx_tc_util/prob_gen_info.cc
+++ b/met/src/libcode/vx_tc_util/prob_gen_info.cc
@@ -17,6 +17,8 @@ using namespace std;
 #include <cstdio>
 #include <cmath>
 
+#include "nav.h"
+
 #include "prob_gen_info.h"
 #include "atcf_offsets.h"
 
@@ -201,3 +203,25 @@ bool ProbGenInfo::add(const ATCFProbLine &l, double dland, bool check_dup) {
 }
 
 ////////////////////////////////////////////////////////////////////////
+
+bool ProbGenInfo::is_match(const TrackPoint &p, const double rad,
+                           const int beg, const int end) const {
+
+   // Check for matching in time and space
+   return(p.valid() >= (GenTime + beg) &&
+          p.valid() <= (GenTime + end) &&
+          gc_dist(Lat, Lon, p.lat(), p.lon()) <= rad);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+bool ProbGenInfo::is_match(const GenesisInfo &gi, const double rad,
+                           const int beg, const int end) const {
+
+   // Input genesis point
+   const TrackPoint *p = gi.genesis();
+
+   return(p ? is_match(*p, rad, beg, end) : false);
+}
+
+////////////////////////////////////////////////////////////////////////
diff --git a/met/src/libcode/vx_tc_util/prob_gen_info.h b/met/src/libcode/vx_tc_util/prob_gen_info.h
index 1ce3dc10d5..49a6534ea1 100644
--- a/met/src/libcode/vx_tc_util/prob_gen_info.h
+++ b/met/src/libcode/vx_tc_util/prob_gen_info.h
@@ -16,6 +16,7 @@
 #include <iostream>
 
 #include "prob_info_base.h"
+#include "genesis_info.h"
 
 #include "vx_util.h"
 
@@ -69,7 +70,12 @@ class ProbGenInfo : public ProbInfoBase {
       void initialize(const ATCFProbLine &, double);
       bool is_match  (const ATCFProbLine &) const;
       bool add       (const ATCFProbLine &, double, bool check_dup = false);
-      
+
+      bool is_match  (const TrackPoint &,
+                      const double, const int, const int) const;
+      bool is_match  (const GenesisInfo &,
+                      const double, const int, const int) const;
+
 };
 
 ////////////////////////////////////////////////////////////////////////
diff --git a/met/src/tools/tc_utils/tc_gen/tc_gen.cc b/met/src/tools/tc_utils/tc_gen/tc_gen.cc
index dcbd87f7c6..eb65f75ec5 100644
--- a/met/src/tools/tc_utils/tc_gen/tc_gen.cc
+++ b/met/src/tools/tc_utils/tc_gen/tc_gen.cc
@@ -95,6 +95,10 @@ static int    find_genesis_match   (const GenesisInfo &,
                                     const GenesisInfoArray &,
                                     const TrackInfoArray &,
                                     bool, double, int, int);
+static int    find_probgen_match   (const ProbGenInfo &,
+                                    const GenesisInfoArray &,
+                                    const TrackInfoArray &,
+                                    bool, double, int, int);
 
 static void   setup_txt_files      (int, int);
 static void   setup_table          (AsciiTable &);
@@ -725,8 +729,9 @@ void do_probgen_pct(const TCGenVxOpt &vx_opt,
                     const GenesisInfoArray &best_ga,
                     const TrackInfoArray &oper_ta,
                     ProbGenPCTInfo &pgi) {
-   int i, j, prob_lead;
+   int i, j, i_bga, prob_lead;
    double prob_value;
+   bool dev_match, ops_match;
 
    // Initialize
    pgi.clear();
@@ -738,11 +743,12 @@ void do_probgen_pct(const TCGenVxOpt &vx_opt,
       // Store info about each probability of genesis forecast
       pgi.add(model_pa.prob_gen(i));
 
-      // JHG search for a match!
-      // For now, just assume it's a hit.
-      // Also need to define SHC values.
-      bool dev_match = true;
-      bool ops_match = true;
+      // Search for a BEST track match
+      i_bga = find_probgen_match(model_pa.prob_gen(i), best_ga, oper_ta,
+                                 vx_opt.GenesisMatchPointTrack,
+                                 vx_opt.GenesisMatchRadius,
+                                 vx_opt.GenesisMatchBeg,
+                                 vx_opt.GenesisMatchEnd);
 
       // Loop over the individual probabilities
       for(j=0; j<model_pa.prob_gen(i).n_prob(); j++) {
@@ -750,6 +756,15 @@ void do_probgen_pct(const TCGenVxOpt &vx_opt,
          prob_lead  = nint(model_pa.prob_gen(i).prob_item(j));
          prob_value = model_pa.prob_gen(i).prob(j) / 100.0;
 
+         // Score no matching BEST track as a non-event
+         if(!is_bad_data(i_bga)) {
+            dev_match = ops_match = false;
+         }
+         // Otherwise, check the time window
+         else {
+            dev_match = ops_match = true;
+         }
+
          // Increment counts
          if(dev_match) pgi.PCTDev[prob_lead].pct.inc_event   (prob_value);
          else          pgi.PCTDev[prob_lead].pct.inc_nonevent(prob_value);
@@ -851,6 +866,94 @@ int find_genesis_match(const GenesisInfo      &fcst_gi,
    return(i_best);
 }
 
+////////////////////////////////////////////////////////////////////////
+
+int find_probgen_match(const ProbGenInfo      &prob_gi,
+                       const GenesisInfoArray &bga,
+                       const TrackInfoArray   &ota,
+                       bool point2track, double rad,
+                       int beg, int end) {
+   int i, j, i_best, i_oper;
+
+   ConcatString case_cs;
+   case_cs << prob_gi.technique() << " "
+           << unix_to_yyyymmdd_hhmmss(prob_gi.init())
+           << " initialization, "
+           << unix_to_yyyymmdd_hhmmss(prob_gi.gen_time())
+           << " forecast genesis at (" << prob_gi.lat() << ", "
+           << prob_gi.lon() << ")";
+
+   // Search for a BEST track genesis match
+   for(i=0, i_best=bad_data_int;
+       i<bga.n() && is_bad_data(i_best);
+       i++) {
+
+      // Check all BEST track points
+      if(point2track) {
+
+         for(j=0; j<bga[i].n_points(); j++) {
+            if(prob_gi.is_match(bga[i][j], rad, beg, end)) {
+               i_best = i;
+               mlog << Debug(4) << case_cs
+                    << " MATCHES BEST genesis track "
+                    << bga[i].storm_id() << ".\n";
+               break;
+            }
+         }
+      }
+      // Check only the BEST genesis points
+      else {
+
+         if(prob_gi.is_match(bga[i], rad, beg, end)) {
+            i_best = i;
+            mlog << Debug(4) << case_cs
+                 << " MATCHES BEST genesis point "
+                 << bga[i].storm_id() << ".\n";
+            break;
+         }
+      }
+   } // end for bga
+
+   // If no BEST track match was found, search the operational tracks
+   if(is_bad_data(i_best)) {
+
+      for(i=0, i_oper=bad_data_int;
+          i<ota.n() && is_bad_data(i_oper);
+          i++) {
+
+         // Each operational track contains only lead time 0
+         if(ota[i].n_points() == 0) continue;
+
+         if(prob_gi.is_match(ota[i][0], rad, beg, end)) {
+            i_oper = i;
+            mlog << Debug(4) << case_cs
+                 << " MATCHES operational " << ota[i].technique()
+                 << " genesis track " << ota[i].storm_id() << ".\n";
+            break;
+         }
+      } // end for ota
+
+      // Find BEST track for this operational track
+      if(!is_bad_data(i_oper)) {
+         for(i=0; i<bga.n(); i++) {
+            if(bga[i].storm_id() == ota[i_oper].storm_id()) {
+               i_best = i;
+               break;
+            }
+         }
+      }
+   }
+
+   // Check for no match
+   if(is_bad_data(i_best)) {
+      mlog << Debug(4) << case_cs
+           << " has NO MATCH in the BEST or operational tracks.\n";
+   }
+
+   return(i_best);
+}
+
+
 ////////////////////////////////////////////////////////////////////////
 
 void get_atcf_files(const StringArray &source,
@@ -1163,6 +1266,19 @@ void process_tracks(const StringArray &files,
    mlog << Debug(2) << "Found " << best_ga.n()
         << " BEST genesis events.\n";
 
+   // Dump out very verbose output
+   if(mlog.verbosity_level() > 6) {
+      mlog << Debug(6) << best_ga.serialize_r() << "\n";
+   }
+   // Dump out track info
+   else {
+      for(i=0; i<best_ga.n(); i++) {
+         mlog << Debug(6) << "[Genesis " << i+1 << " of "
+              << best_ga.n() << "] " << best_ga[i].serialize()
+              << "\n";
+      }
+   }
+
    return;
 }
 
diff --git a/met/src/tools/tc_utils/tc_pairs/Makefile.am b/met/src/tools/tc_utils/tc_pairs/Makefile.am
index d5d4436698..1dd9e1275f 100644
--- a/met/src/tools/tc_utils/tc_pairs/Makefile.am
+++ b/met/src/tools/tc_utils/tc_pairs/Makefile.am
@@ -38,6 +38,7 @@ tc_pairs_LDADD = -lvx_stat_out \
 	-lvx_gsl_prob \
 	-lvx_pb_util \
 	-lvx_cal \
+	-lvx_nav \
 	-lvx_util \
 	-lvx_math \
 	-lvx_color \