PRVI GMM epistemic uncertainty

src/main/java/gov/usgs/earthquake/nshmp/gmm/UsgsPrviBackbone2025.java

@@ -2,14 +2,18 @@ package gov.usgs.earthquake.nshmp.gmm;
 
 import static gov.usgs.earthquake.nshmp.gmm.Gmm.AG_20_GLOBAL_INTERFACE;
 import static gov.usgs.earthquake.nshmp.gmm.Gmm.AG_20_GLOBAL_SLAB;
-import static gov.usgs.earthquake.nshmp.gmm.Gmm.ASK_14;
-import static gov.usgs.earthquake.nshmp.gmm.Gmm.BSSA_14;
-import static gov.usgs.earthquake.nshmp.gmm.Gmm.CB_14;
-import static gov.usgs.earthquake.nshmp.gmm.Gmm.CY_14;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.ASK_14_BASE;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.BSSA_14_BASE;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.CB_14_BASE;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.CY_14_BASE;
 import static gov.usgs.earthquake.nshmp.gmm.Gmm.KBCG_20_GLOBAL_INTERFACE;
 import static gov.usgs.earthquake.nshmp.gmm.Gmm.KBCG_20_GLOBAL_SLAB;
 import static gov.usgs.earthquake.nshmp.gmm.Gmm.PSBAH_20_GLOBAL_INTERFACE;
 import static gov.usgs.earthquake.nshmp.gmm.Gmm.PSBAH_20_GLOBAL_SLAB;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.Type.ACTIVE_CRUST;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.Type.SUBDUCTION_INTERFACE;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.Type.SUBDUCTION_SLAB;
+import static gov.usgs.earthquake.nshmp.gmm.GroundMotionTables.getPrviEpi;
 import static gov.usgs.earthquake.nshmp.gmm.SubductionZone.Type.INTERFACE;
 import static gov.usgs.earthquake.nshmp.gmm.SubductionZone.Type.SLAB;
 import static java.lang.Math.exp;
@@ -19,6 +23,7 @@ import static java.lang.Math.sqrt;
 import java.util.Map;
 
 import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
+import gov.usgs.earthquake.nshmp.gmm.GroundMotionTables.GroundMotionTable;
 import gov.usgs.earthquake.nshmp.tree.Branch;
 import gov.usgs.earthquake.nshmp.tree.LogicTree;
 
@@ -57,14 +62,12 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
   static final CoefficientContainer COEFFS_DATA_ADJUSTMENT =
       new CoefficientContainer("prvi-25-backbone-adjustments.csv");
 
-  // average bias over all IMTs
-  // private final double ACTIVE_CRUSTAL_MEAN_ADJUSTMENT = -0.371;
-  // private final double SUBDUCTION_INTERFACE_MEAN_ADJUSTMENT = -1.098;
-  // private final double SUBDUCTION_INTRASLAB_MEAN_ADJUSTMENT = -0.405;
+  public static final String MODEL_BASE_ID = "base";
+  public static final String MODEL_ADJUST_ID = "adjust";
 
-  public static final String SIGMA_REF_ID = "σ-model1";
-  public static final String SIGMA_PRVI_ID = "σ-model2";
-  private static final String[] SIGMA_IDS = { SIGMA_REF_ID, SIGMA_PRVI_ID };
+  public static final String SIGMA_NGA_ID = "σ-nga";
+  public static final String SIGMA_PRVI_ID = "σ-prvi";
+  private static final String[] SIGMA_IDS = { SIGMA_NGA_ID, SIGMA_PRVI_ID };
   private static final double[] SIGMA_WTS = { 0.5, 0.5 };
 
   private static final class Coefficients {
@@ -74,11 +77,6 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
     final double φs2sPrvi, φs2sWus;
     final double τSub, φSub, φSubSS, φSubS2S;
 
-    final double εSubLo = 0.3;
-    final double εSubHi = 0.3;
-    final double εCrustalLo = 0.3;
-    final double εCrustalHi = 0.3;
-
     Coefficients(
         Imt imt,
         CoefficientContainer ccTau,
@@ -111,17 +109,18 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
 
   private final Coefficients coeffs;
 
-  private final Map<Gmm, Double> gmms;
   private final LogicTree<GroundMotionModel> tree;
   private final Imt imt;
+  private final GroundMotionTable[] epiTables; // 0=lower, 1=upper
 
   /* Supply map of ground motion models initialized to the required IMT. */
   UsgsPrviBackbone2025(
       Imt imt,
+      String name,
       Map<Gmm, Double> gmms,
-      String name) {
-    this.gmms = gmms;
+      GroundMotionTable[] epiTables) {
     this.imt = imt;
+    this.epiTables = epiTables;
     LogicTree.Builder<GroundMotionModel> b = LogicTree.builder(name);
     gmms.entrySet().stream().forEach(e -> b.addBranch(
         e.getKey().name(),
@@ -142,15 +141,11 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
     return imt;
   }
 
-  abstract double[] calcSigmas(GmmInput in);
-
-  abstract double[] applyEpistemic(double mu);
-
   @Override
   public LogicTree<GroundMotion> calc(GmmInput in) {
 
     double μ = calcMean(in);
-    double[] μs = applyEpistemic(μ);
+    double[] μs = applyEpistemic(μ, in.Mw, in.rRup);
     double[] σs = calcSigmas(in);
 
     return GroundMotions.createTree(
@@ -169,13 +164,17 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
           : gmTree.get(1).value();
       median += exp(gm.mean()) * b.weight();
     }
+    return log(median);
+  }
 
-    // epistemic uncertainty -
-    // double[] μs = GroundMotions.create5th50th95th(log(median), median);
-    double m2 = log(median);
-    return m2;
+  double[] applyEpistemic(double μ, double Mw, double rRup) {
+    double εLo = epiTables[0].get(rRup, Mw);
+    double εHi = epiTables[1].get(rRup, Mw);
+    return new double[] { μ + εLo, μ, μ + εHi };
   }
 
+  abstract double[] calcSigmas(GmmInput in);
+
   /*
    * Implementations. For each, coefficients are only used to get the set of
    * supported IMTs and therefore reference a model that supports the
@@ -185,21 +184,21 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
   /* PRVI 2025 Active Crust */
   static class ActiveCrust extends UsgsPrviBackbone2025 {
 
-    static final String NAME = UsgsPrviBackbone2025.NAME + ": Active Crust";
+    static final String NAME = UsgsPrviBackbone2025.NAME + " : Active Crust";
     static final Constraints CONSTRAINTS = AbrahamsonEtAl_2014.CONSTRAINTS;
 
     static final Map<Gmm, Double> GMM_MAP = Map.of(
-        ASK_14, 0.25,
-        BSSA_14, 0.25,
-        CB_14, 0.25,
-        CY_14, 0.25);
+        ASK_14_BASE, 0.25,
+        BSSA_14_BASE, 0.25,
+        CB_14_BASE, 0.25,
+        CY_14_BASE, 0.25);
 
     ActiveCrust(Imt imt) {
-      super(imt, GMM_MAP, NAME);
+      this(imt, NAME);
     }
 
-    ActiveCrust(Imt imt, Map<Gmm, Double> gmms, String name) {
-      super(imt, gmms, name);
+    ActiveCrust(Imt imt, String name) {
+      super(imt, name, GMM_MAP, getPrviEpi(ACTIVE_CRUST, imt));
     }
 
     @Override
@@ -231,12 +230,6 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
     static double calcPhiSS(double Mw, Coefficients c) {
       return Mw <= 5.0 ? c.a : Mw > 6.5 ? c.b : c.a + (Mw - 5.0) * ((c.b - c.a) / 1.5);
     }
-
-    @Override
-    double[] applyEpistemic(double μ) {
-      Coefficients c = super.coeffs;
-      return new double[] { μ - c.εCrustalLo, μ, μ + c.εCrustalHi };
-    }
   }
 
   /* PRVI 2025 Active Crust Data Adjusted */
@@ -247,7 +240,7 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
     final double bias;
 
     ActiveCrustAdjusted(Imt imt) {
-      super(imt, ActiveCrust.GMM_MAP, NAME);
+      super(imt, NAME);
       bias = COEFFS_DATA_ADJUSTMENT.get(imt, "active_crust");
     }
 
@@ -257,10 +250,38 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
     }
   }
 
+  /* TODO Deprecate */
+  static final class ActiveCrustTotal implements GroundMotionModel {
+    static final String NAME = ActiveCrust.NAME + " (total tree)";
+    static final Constraints CONSTRAINTS = AbrahamsonEtAl_2014.CONSTRAINTS;
+
+    private final Imt imt;
+    private final ActiveCrust gmm1;
+    private final ActiveCrustAdjusted gmm2;
+
+    ActiveCrustTotal(Imt imt) {
+      this.imt = imt;
+      gmm1 = new ActiveCrust(imt);
+      gmm2 = new ActiveCrustAdjusted(imt);
+    }
+
+    @Override
+    public Imt imt() {
+      return imt;
+    }
+
+    @Override
+    public LogicTree<GroundMotion> calc(GmmInput in) {
+      LogicTree<GroundMotion> treeBase = gmm1.calc(in);
+      LogicTree<GroundMotion> treeAdj = gmm2.calc(in);
+      return combine(treeBase, treeAdj, "active-crust");
+    }
+  }
+
   /* PRVI 2025 Subduction Interface */
   static class Interface extends UsgsPrviBackbone2025 {
 
-    static final String NAME = UsgsPrviBackbone2025.NAME + ":" + INTERFACE;
+    static final String NAME = UsgsPrviBackbone2025.NAME + " : " + INTERFACE;
     static final Constraints CONSTRAINTS = AbrahamsonGulerce_2020.CONSTRAINTS_INTERFACE;
 
     static final Map<Gmm, Double> GMM_MAP = Map.of(
@@ -269,11 +290,12 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
         PSBAH_20_GLOBAL_INTERFACE, 0.3334);
 
     Interface(Imt imt) {
-      super(imt, GMM_MAP, NAME);
+      this(imt, NAME);
     }
 
-    Interface(Imt imt, Map<Gmm, Double> gmms, String name) {
-      super(imt, gmms, name);
+    Interface(Imt imt, String name) {
+      super(imt, name, GMM_MAP, getPrviEpi(SUBDUCTION_INTERFACE, imt));
+
     }
 
     @Override
@@ -283,12 +305,6 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
       double σPrvi = sqrt(c.τSub * c.τSub + c.φSubS2S * c.φSubS2S + c.φSubSS * c.φSubSS);
       return new double[] { σ, σPrvi };
     }
-
-    @Override
-    double[] applyEpistemic(double μ) {
-      Coefficients c = super.coeffs;
-      return new double[] { μ - c.εSubLo, μ, μ + c.εSubHi };
-    }
   }
 
   /* PRVI 2025 Interface Data Adjusted */
@@ -299,7 +315,7 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
     final double bias;
 
     InterfaceAdjusted(Imt imt) {
-      super(imt, Interface.GMM_MAP, NAME);
+      super(imt, NAME);
       bias = COEFFS_DATA_ADJUSTMENT.get(imt, "interface");
     }
 
@@ -309,10 +325,38 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
     }
   }
 
+  /* TODO Deprecate */
+  static final class InterfaceTotal implements GroundMotionModel {
+    static final String NAME = Interface.NAME + " (total tree)";
+    static final Constraints CONSTRAINTS = AbrahamsonGulerce_2020.CONSTRAINTS_INTERFACE;
+
+    private final Imt imt;
+    private final Interface gmm1;
+    private final InterfaceAdjusted gmm2;
+
+    InterfaceTotal(Imt imt) {
+      this.imt = imt;
+      gmm1 = new Interface(imt);
+      gmm2 = new InterfaceAdjusted(imt);
+    }
+
+    @Override
+    public Imt imt() {
+      return imt;
+    }
+
+    @Override
+    public LogicTree<GroundMotion> calc(GmmInput in) {
+      LogicTree<GroundMotion> treeBase = gmm1.calc(in);
+      LogicTree<GroundMotion> treeAdj = gmm2.calc(in);
+      return combine(treeBase, treeAdj, "interface");
+    }
+  }
+
   /* PRVI 2025 Subduction Intraslab */
   static class Slab extends UsgsPrviBackbone2025 {
 
-    static final String NAME = UsgsPrviBackbone2025.NAME + ":" + SLAB;
+    static final String NAME = UsgsPrviBackbone2025.NAME + " : " + SLAB;
     static final Constraints CONSTRAINTS = AbrahamsonGulerce_2020.CONSTRAINTS_SLAB;
 
     static final Map<Gmm, Double> GMM_MAP = Map.of(
@@ -321,11 +365,11 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
         PSBAH_20_GLOBAL_SLAB, 0.3334);
 
     Slab(Imt imt) {
-      super(imt, GMM_MAP, NAME);
+      this(imt, NAME);
     }
 
-    Slab(Imt imt, Map<Gmm, Double> gmms, String name) {
-      super(imt, gmms, name);
+    Slab(Imt imt, String name) {
+      super(imt, name, GMM_MAP, getPrviEpi(SUBDUCTION_SLAB, imt));
     }
 
     @Override
@@ -335,12 +379,6 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
       double σPrvi = sqrt(c.τSub * c.τSub + c.φSubS2S * c.φSubS2S + c.φSubSS * c.φSubSS);
       return new double[] { σ, σPrvi };
     }
-
-    @Override
-    double[] applyEpistemic(double μ) {
-      Coefficients c = super.coeffs;
-      return new double[] { μ - c.εSubLo, μ, μ + c.εSubHi };
-    }
   }
 
   /* PRVI 2025 Interface Data Adjusted */
@@ -351,7 +389,7 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
     final double bias;
 
     SlabAdjusted(Imt imt) {
-      super(imt, Slab.GMM_MAP, NAME);
+      super(imt, NAME);
       bias = COEFFS_DATA_ADJUSTMENT.get(imt, "intraslab");
     }
 
@@ -360,4 +398,57 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel {
       return super.calcMean(in) + bias;
     }
   }
+
+  /* TODO Deprecate */
+  static final class SlabTotal implements GroundMotionModel {
+    static final String NAME = Slab.NAME + " (total tree)";
+    static final Constraints CONSTRAINTS = AbrahamsonGulerce_2020.CONSTRAINTS_SLAB;
+
+    private final Imt imt;
+    private final Slab gmm1;
+    private final SlabAdjusted gmm2;
+
+    SlabTotal(Imt imt) {
+      this.imt = imt;
+      gmm1 = new Slab(imt);
+      gmm2 = new SlabAdjusted(imt);
+    }
+
+    @Override
+    public Imt imt() {
+      return imt;
+    }
+
+    @Override
+    public LogicTree<GroundMotion> calc(GmmInput in) {
+      LogicTree<GroundMotion> treeBase = gmm1.calc(in);
+      LogicTree<GroundMotion> treeAdj = gmm2.calc(in);
+      return combine(treeBase, treeAdj, "intraslab");
+    }
+  }
+
+  /*
+   * Hopefully temporary support for collapsed/combined GMM logic trees. For
+   * PRVI, kmilner requested access to a single logic tree object for each
+   * tectonic setting.
+   */
+
+  static final double WT_SCALE = 0.5;
+
+  /* Merge base and adjusted model logic trees */
+  static LogicTree<GroundMotion> combine(
+      LogicTree<GroundMotion> base,
+      LogicTree<GroundMotion> adj,
+      String name) {
+
+    LogicTree.Builder<GroundMotion> tree = LogicTree.builder("gmm-tree-prvi-" + name);
+    for (Branch<GroundMotion> b : base) {
+      tree.addBranch(MODEL_BASE_ID + "-" + b.id(), b.value(), b.weight() * 0.5);
+    }
+    for (Branch<GroundMotion> b : adj) {
+      tree.addBranch(MODEL_ADJUST_ID + "-" + b.id(), b.value(), b.weight() * 0.5);
+    }
+    return tree.build();
+  }
+
 }