From 2a212aca66d109932dabae5ce085533ffe7b3e8e Mon Sep 17 00:00:00 2001
From: Peter Powers <pmpowers@usgs.gov>
Date: Wed, 21 Feb 2024 08:56:55 -0700
Subject: [PATCH] fixes to individual seeds to handle PGV

---
 .../usgs/earthquake/nshmp/gmm/NgaEast.java    | 57 ++++++++++++++-----
 1 file changed, 44 insertions(+), 13 deletions(-)

diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/NgaEast.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/NgaEast.java
index c68317b8..2b7bedf7 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/NgaEast.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/NgaEast.java
@@ -13,7 +13,7 @@ import static java.lang.Math.exp;
 import static java.lang.Math.log;
 import static java.lang.Math.min;
 import static java.lang.Math.sqrt;
-import static java.util.stream.Collectors.toMap;
+import static java.util.stream.Collectors.toUnmodifiableMap;
 
 import java.io.IOException;
 import java.net.URL;
@@ -197,11 +197,11 @@ public abstract class NgaEast implements GroundMotionModel {
           .stream()
           .skip(1)
           .map(line -> Text.splitToList(line, Delimiter.COMMA))
-          .collect(toMap(
-              entry -> entry.get(0),
-              entry -> Double.valueOf(entry.get(1))));
+          .collect(toUnmodifiableMap(
+              e -> e.get(0),
+              e -> Double.valueOf(e.get(1))));
       checkWeights(wtMap.values());
-      return Map.copyOf(wtMap);
+      return wtMap;
 
     } catch (IOException ioe) {
       throw new RuntimeException(ioe);
@@ -556,15 +556,26 @@ public abstract class NgaEast implements GroundMotionModel {
     static final String NAME = NgaEast.NAME + " Seed Tree (2018)";
     static final String SEED_PREFIX = "NGA_EAST_SEED_";
     static final String SP16_ID = "SP16";
-    static final Set<Gmm> noPgvSeeds = EnumSet.of(
-        Gmm.NGA_EAST_SEED_PEER_GP,
+    static final Set<Gmm> NO_PGV_SEEDS = EnumSet.of(
+        Gmm.NGA_EAST_SEED_GRAIZER, // not used
         Gmm.NGA_EAST_SEED_GRAIZER16,
         Gmm.NGA_EAST_SEED_GRAIZER17,
+        Gmm.NGA_EAST_SEED_PEER_EX, // not used
+        Gmm.NGA_EAST_SEED_PEER_GP,
         Gmm.NGA_EAST_SEED_PZCT15_M1SS,
         Gmm.NGA_EAST_SEED_PZCT15_M2ES);
+    static final Set<String> NO_PGV_IDS = Set.of(
+        "Graizer",
+        "Graizer16",
+        "Graizer17",
+        "PEER_EX",
+        "PEER_GP",
+        "PZCT15_M1SS",
+        "PZCT15_M2ES");
 
     /* Valid Seed IDs for table based models; skips SP16 */
     static final List<String> SEED_TABLE_IDS; // size = 13
+    static final List<String> SEED_TABLE_IDS_PGV; // size = 8
 
     /* All Gmms */
     static final List<Gmm> GMMS; // pgv helper, size = 14
@@ -581,6 +592,7 @@ public abstract class NgaEast implements GroundMotionModel {
 
     static {
       List<String> tableIds = new ArrayList<>(); // valid IDs for tables
+      List<String> tableIdsPgv = new ArrayList<>(); // valid IDs for PGV tables
       List<Gmm> gmms = new ArrayList<>();
       List<String> μIds = new ArrayList<>();
       List<Double> μWts = new ArrayList<>();
@@ -588,12 +600,16 @@ public abstract class NgaEast implements GroundMotionModel {
         String seedId = e.getKey();
         if (!seedId.equals(SP16_ID)) {
           tableIds.add(seedId);
+          if (!NO_PGV_IDS.contains(seedId)) {
+            tableIdsPgv.add(seedId);
+          }
         }
         gmms.add(Gmm.valueOf(SEED_PREFIX + e.getKey().toUpperCase()));
         μIds.add("Seed (" + seedId + ")");
         μWts.add(e.getValue());
       });
       SEED_TABLE_IDS = List.copyOf(tableIds);
+      SEED_TABLE_IDS_PGV = List.copyOf(tableIdsPgv);
       GMMS = List.copyOf(gmms);
       MEAN_IDS = μIds.toArray(new String[μIds.size()]);
       MEAN_WTS = μWts.stream()
@@ -615,9 +631,10 @@ public abstract class NgaEast implements GroundMotionModel {
     }
 
     private static Map<Gmm, GroundMotionTable> getImtTables(Imt imt) {
-      return GroundMotionTables.getNgaEastSeeds(SEED_TABLE_IDS, imt)
+      List<String> ids = (imt == Imt.PGV) ? SEED_TABLE_IDS_PGV : SEED_TABLE_IDS;
+      return GroundMotionTables.getNgaEastSeeds(ids, imt)
           .entrySet().stream()
-          .collect(toMap(
+          .collect(toUnmodifiableMap(
               e -> Gmm.valueOf(SEED_PREFIX + e.getKey().toUpperCase()),
               Entry::getValue));
     }
@@ -633,7 +650,7 @@ public abstract class NgaEast implements GroundMotionModel {
           double μPga = exp(sp16pga.calcMeanRock(in.Mw, in.rJB));
           SiteTerm fSite = siteAmp.calc(μPga, in.vs30);
           μs[i] = fSite.apply(μ);
-        } else if (imt == Imt.PGV && noPgvSeeds.contains(seed)) {
+        } else if (imt == Imt.PGV && NO_PGV_SEEDS.contains(seed)) {
           // site is included when calling individual seed
           μs[i] = UsgsPgvSupport.calcAB20Pgv(seed, in).mean();
         } else {
@@ -688,7 +705,7 @@ public abstract class NgaEast implements GroundMotionModel {
         if (seed == Gmm.NGA_EAST_SEED_SP16) {
           μ = sp16.calcMeanRock(in.Mw, in.rJB);
           μPga = sp16pga.calcMeanRock(in.Mw, in.rJB);
-        } else if (imt == Imt.PGV && noPgvSeeds.contains(seed)) {
+        } else if (imt == Imt.PGV && NO_PGV_SEEDS.contains(seed)) {
           // note rock input; TODO this doesn't really work as the PGV
           // is conditioned on vs30 so site should already be applied
           μ = UsgsPgvSupport.calcAB20Pgv(seed, inRock).mean();
@@ -766,10 +783,24 @@ public abstract class NgaEast implements GroundMotionModel {
     public LogicTree<GroundMotion> calc(GmmInput in) {
       double σEpri = sigmaEpri(in.Mw);
       double σPanel = sigmaPanel(in.Mw, in.vs30);
-      Position p = table.position(in.rRup, in.Mw);
+      // PGV fails if 'table' used; all seeds support PGA natively
+      Position p = pgaTable.position(in.rRup, in.Mw);
       double μPga = exp(pgaTable.get(p));
       SiteTerm fSite = siteAmp.calc(μPga, in.vs30);
-      double μ = fSite.apply(table.get(p));
+
+      Gmm seedGmm = Gmm.valueOf(UsgsSeeds_2018.SEED_PREFIX + id.toUpperCase());
+      GmmInput inRock = GmmInput.builder().fromCopy(in).vs30(3000).build();
+      double μ = (imt == Imt.PGV && UsgsSeeds_2018.NO_PGV_SEEDS.contains(seedGmm))
+          ? UsgsPgvSupport.calcAB20Pgv(seedGmm, inRock).mean()
+          : table.get(p);
+      μ = fSite.apply(μ);
+
+      /*
+       * TODO the order of operations for PGV for individual seeds is consistent
+       * with seed 2023 but not seed 2018 and needs to be rechecked. Conditional
+       * PGV is dependent on vs30; should we realy be using the consitional
+       * model with vs30=3000
+       */
       double σ = Maths.srssWeighted(
           new double[] { σEpri, σPanel },
           SIGMA_WTS);
-- 
GitLab