fixed issues with median gm calc, used expected results from official python...

fixed issues with median gm calc, used expected results from official python implementation (modified to use updated coeffs). uncertainty not yet implemented

fixed issues with median gm calc, used expected results from official python...
999dba56 · Jason Altekruse · 3e08c026 · 999dba56 · 999dba56 · 999dba56
Commit 999dba56 authored 4 years ago by Jason Altekruse
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
@@ -21,6 +21,7 @@ import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;

 /**
 * Implementation of the PEER NGA-Subduction GMM, Parker et al., 2020 GMM
+ * (PSHAB20)
 *
 * Notes:
 *
@@ -28,6 +29,11 @@ import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
 * analysis is planned to determine if there are attenuation differences between
 * forearc and backarc zones.
 *
+ * PSHAB20 uses zHyp (randomized source location on fault), but this can be
+ * replaced this event-specific zHyp with the mean depth expected for the fault
+ * plane, zHypBar, in Eq. 4.6 without modification to coefficients or
+ * uncertainty.
+ *
 * <p><b>Reference:</b> Parker, G.A, Stewart, J.P., Boore, D.M., Atkinson, G.M.,
 * Hassani, B. (2020) NGA-Subduction global ground-motion models with regional
 * adjustment factors. Report no. 2020/03. Berkeley, CA: PEER, 131 p.
@@ -37,14 +43,6 @@ import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
 * @author U.S. Geological Survey
 */

-/*
- * Regions: global, Alaska, Cascadia. Unused Regions: CAM (Central America),
- * Japan, SA (South America), Taiwan
- *
- * Saturation Regions: global, Aleutian, Alaska, Cascadia. Unused Regions:
- * Central America S, Central America N, Japan Pacific, Japan Phillipines, South
- * America N, South America S, Taiwan W, Taiwan E
- */
 public abstract class ParkerEtAl_2020 implements GroundMotionModel {

  static final String NAME = "Parker et al. (2020) NGA-Subduction";
@@ -62,6 +60,14 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
      .set(VS30, Range.closed(150.0, 2000.0))
      .build();

+  /*
+   * Regions: global, Alaska, Cascadia, [CAM (Central America), Japan, SA (South
+   * America), Taiwan]
+   *
+   * Saturation Regions: global, Aleutian, Alaska, Cascadia, [Central America S,
+   * Central America N, Japan Pacific, Japan Phillipines, South America N, South
+   * America S, Taiwan W, Taiwan E]
+   */
  static final class Region {
    static final String GLOBAL = "Global";
    static final String ALASKA = "Alaska";
@@ -253,11 +259,21 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
      final double Mw,
      final double rRup,
      final double zHyp) {
-    double h = nearSourceSaturationTerm(Mw, c.mc, isSlab);
-    return c.c0 +
-        getPathTerm(c, rRup, Mw, h) +
-        getMagnitudeTerm(c, Mw) +
-        getDepthTerm(c, isSlab, zHyp);
+    double h = nearSourceSaturation(Mw, c.mc, isSlab);
+    double Fp = getPathTermPy(c, rRup, Mw, h);
+    double Fm = getMagnitudeTerm(c, Mw);
+    double Fd = getDepthTerm(c, isSlab, zHyp);
+
+    // System.out.println("calcMean Reference for T = " + c.imt.toString());
+    // System.out.printf(" mc = %.2f\n", c.mc);
+    // System.out.printf(" h = %.6f\n", h);
+    // System.out.printf(" Fp = %.6f\n", Fp);
+    System.out.printf("  Fm(ref) = %.6f\n", Fm);
+    // System.out.printf(" Fd = %.6f\n", Fd);
+    // System.out.printf(" Fb = %.6f (reference calc)\n", 0.0);
+    // System.out.printf(" mu = %.6f\n", c.c0 + Fp + Fm + Fd);
+
+    return c.c0 + Fp + Fm + Fd;
  }

  /* Equation 4.1 */
@@ -272,28 +288,57 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
      final double vs30,
      final double z2p5) {

-    /* Equation 4.7: Fs = Flin + Fnl + Fb */
-    double siteTerm = getSiteTermLinear(c, vs30) + getSiteTermNonLinear(c, pgaRef, vs30) +
-        getSiteBasinTerm(c, basin, vs30, z2p5);
-
-    double h = nearSourceSaturationTerm(Mw, c.mc, isSlab);
-
-    double μ = c.c0 +
-        getPathTerm(c, rRup, Mw, h) +
-        getMagnitudeTerm(c, Mw) +
-        getDepthTerm(c, isSlab, zHyp) +
-        siteTerm;
+    double h = nearSourceSaturation(Mw, c.mc, isSlab);
+    double Fp = getPathTermPy(c, rRup, Mw, h);
+    double Fm = getMagnitudeTerm(c, Mw);
+    double Fd = getDepthTerm(c, isSlab, zHyp);

-    return μ;
+    /* Equation 4.7: Fs = Flin + Fnl + Fb */
+    double Fslin = getSiteTermLinear(c, vs30);
+    double Fsnl = getSiteTermNonLinear(c, pgaRef, vs30);
+    double Fsb = getSiteBasinTerm(c, basin, vs30, z2p5);
+    double Fs = Fslin + Fsnl + Fsb;
+
+    System.out.println("calcMean for T = " + c.imt.toString());
+    System.out.printf("  Vs30 = %.2f\n", vs30);
+    System.out.printf("  h = %.6f\n", h);
+    System.out.printf("  Fp = %.6f\n", Fp);
+    System.out.printf("  Fm = %.6f\n", Fm);
+    System.out.printf("  Fd = %.6f\n", Fd);
+    System.out.printf("  PGAref = %.6f\n", pgaRef);
+    System.out.printf("  Fslin = %.6e\n", Fslin);
+    System.out.printf("  Fsnl = %.6e\n", Fsnl);
+    System.out.printf("  Fssb = %.6e\n", Fsb);
+    System.out.printf("  Fs = %.6e\n", Fs);
+    System.out.printf("  mu = %.6f\n", c.c0 + Fp + Fm + Fd + Fs);
+    System.out.printf("  y = %.6f\n", exp(c.c0 + Fp + Fm + Fd + Fs));
+
+    return c.c0 + Fp + Fm + Fd + Fs;
  }

  /* Equations 4.2 and 4.3 */
  private static final double getPathTerm(final Coefficients c, double rRup, double Mw, double h) {
    double R = Maths.hypot(rRup, h);
-    return c.c1 * log(R) + c.b4 * Mw * log(R) + c.a0 * R;
+    return c.c1 * log(R) + c.b4 * Mw * log(rRup / R) + c.a0 * R;
+  }
+
+  /* Equations 4.2 and 4.3, following official Python implementation */
+  private static final double getPathTermPy(final Coefficients c, double rRup, double Mw,
+      double h) {
+    double Rref = Maths.hypot(1, h);
+    double R = Maths.hypot(rRup, h);
+    double R_Rref = R / Rref;
+    // System.out.printf(" c1 = %.6e\n", c.c1);
+    // System.out.printf(" b4 = %.6e\n", c.b4);
+    // System.out.printf(" a0 = %.6e\n", c.a0);
+    // System.out.printf(" Rref = %.6e\n", R_Rref);
+    // System.out.printf(" R = %.6e\n", R);
+    // System.out.printf(" Mw = %.6e\n", Mw);
+
+    return c.c1 * log(R) + c.b4 * Mw * log(R_Rref) + c.a0 * R;
  }

-  private static final double nearSourceSaturationTerm(double Mw, double mc, boolean slab) {
+  private static final double nearSourceSaturation(double Mw, double mc, boolean slab) {
    if (slab) {
      /* Equation 4.4b */
      return (Mw <= mc) ? pow(10, (1.050 / (mc - 4.0)) * (Mw - mc) + 1.544) : 35;

--- a/src/main/resources/gmm/coeffs/PSHAB20_slab.csv
+++ b/src/main/resources/gmm/coeffs/PSHAB20_slab.csv
 T,Global_c0,Alaska_c0,Aleutian_c0,Cascadia_c0,c1,b4,Global_a0,Alaska_a0,Cascadia_a0,c4,c5,c6,d,m,db,V1,V2,Vref,Global_s2,Alaska_s2,Cascadia_s2,f4,f5,C_e1,C_e2,C_e3,del_None,del_Seattle,Tau,phi21,phi22,phi2V,VM,phi2S2S0,a1,phi2SS1,phi2SS2,a2
 PGV,13.194,12.79,13.6,12.874,-2.422,0.1,-1.90E-03,-2.38E-03,-1.09E-03,1.84,-0.05,0.8,0.2693,0.0252,67,270,850,760,-0.601,-1.031,-0.671,-0.31763,-0.0052,0,0.115,0.068,-0.115,0,0.477,0.348,0.288,-0.179,423,0.142,0.047,0.153,0.166,0.011
-PGV,9.907,9.404,9.912,9.6,-2.543,0.1,-2.55E-03,-2.27E-03,-3.54E-03,1.84,-0.05,0.4,0.3004,0.0314,67,270,1350,760,-0.498,-0.785,-0.572,-0.44169,-0.0052,0,0,1,0,0,0.48,0.396,0.565,-0.18,423,0.221,0.093,0.149,0.327,0.068
+PGA,9.907,9.404,9.912,9.6,-2.543,0.1,-2.55E-03,-2.27E-03,-3.54E-03,1.84,-0.05,0.4,0.3004,0.0314,67,270,1350,760,-0.498,-0.785,-0.572,-0.44169,-0.0052,0,0,1,0,0,0.48,0.396,0.565,-0.18,423,0.221,0.093,0.149,0.327,0.068
 0.01,9.962,9.451,9.954,9.802,-2.554,0.1,-2.55E-03,-2.19E-03,-4.01E-03,1.84,-0.05,0.4,0.2839,0.0296,67,270,1300,760,-0.498,-0.803,-0.571,-0.4859,-0.0052,0,0,1,0,0,0.476,0.397,0.56,-0.18,423,0.223,0.098,0.148,0.294,0.071
 0.02,10.099,9.587,10.086,9.933,-2.566,0.1,-2.55E-03,-2.19E-03,-4.01E-03,1.884,-0.05,0.415,0.2854,0.0298,67,270,1225,760,-0.478,-0.785,-0.575,-0.4859,-0.00518,0,0,1,0,0,0.482,0.401,0.563,-0.181,423,0.227,0.105,0.149,0.294,0.073
 0.025,10.181,9.667,10.172,10.009,-2.578,0.1,-2.55E-03,-2.19E-03,-4.01E-03,1.884,-0.05,0.43,0.2891,0.0302,67,270,1200,760,-0.464,-0.745,-0.573,-0.4859,-0.00515,0,0,1,0,0,0.49,0.405,0.575,-0.183,423,0.231,0.12,0.15,0.31,0.076

--- a/src/test/java/gov/usgs/earthquake/nshmp/gmm/PSHAB20.java
+++ b/src/test/java/gov/usgs/earthquake/nshmp/gmm/PSHAB20.java
 package gov.usgs.earthquake.nshmp.gmm;

 import static gov.usgs.earthquake.nshmp.gmm.Gmm.PSHAB20_INTERFACE;
+import static gov.usgs.earthquake.nshmp.gmm.Imt.PGA;

 import java.io.IOException;
 import java.util.EnumSet;
@@ -33,7 +34,7 @@ class PSHAB20 implements ArgumentsProvider {

  /* Result generation sets */
  private static Set<Gmm> gmms = EnumSet.of(PSHAB20_INTERFACE);
-  private static Set<Imt> imts = PSHAB20_INTERFACE.supportedIMTs();
+  private static Set<Imt> imts = EnumSet.of(PGA);// PSHAB20_INTERFACE.supportedIMTs();

  public static void main(String[] args) throws IOException {
    GmmTest.generateResults(gmms, imts, GMM_INPUTS, GMM_RESULTS);

--- a/src/test/resources/gmm/PSHAB20_results.csv
+++ b/src/test/resources/gmm/PSHAB20_results.csv
-0-PSHAB20_INTERFACE-SA0P01,0.020897,0.885198
-0-PSHAB20_INTERFACE-SA0P02,0.021238,0.890687
-0-PSHAB20_INTERFACE-SA0P03,0.022458,0.907193
-0-PSHAB20_INTERFACE-SA0P05,0.025828,0.954874
-0-PSHAB20_INTERFACE-SA0P075,0.029636,0.984835
-0-PSHAB20_INTERFACE-SA0P1,0.032537,0.975488
-0-PSHAB20_INTERFACE-SA0P15,0.035565,0.936536
-0-PSHAB20_INTERFACE-SA0P2,0.036517,0.918151
-0-PSHAB20_INTERFACE-SA0P25,0.037127,0.923073
-0-PSHAB20_INTERFACE-SA0P3,0.037353,0.915458
-0-PSHAB20_INTERFACE-SA0P4,0.034904,0.913818
-0-PSHAB20_INTERFACE-SA0P5,0.031459,0.911078
-0-PSHAB20_INTERFACE-SA0P75,0.025667,0.921989
-0-PSHAB20_INTERFACE-SA1P0,0.020612,0.912175
-0-PSHAB20_INTERFACE-SA1P5,0.014567,0.898924
-0-PSHAB20_INTERFACE-SA2P0,0.010897,0.891664
-0-PSHAB20_INTERFACE-SA3P0,0.006269,0.876963
-0-PSHAB20_INTERFACE-SA4P0,0.004365,0.855607
-0-PSHAB20_INTERFACE-SA5P0,0.003441,0.836698
-0-PSHAB20_INTERFACE-SA7P5,0.001955,0.809978
-0-PSHAB20_INTERFACE-SA10P0,0.001252,0.795653
-0-PSHAB20_INTERFACE-PGA,0.020287,0.886792
-0-PSHAB20_INTERFACE-PGV,2.053897,0.857630
+0-PSHAB20_INTERFACE-SA0P01,2.08965604936e-02,0.0
+0-PSHAB20_INTERFACE-SA0P02,2.12376645069e-02,0.0
+0-PSHAB20_INTERFACE-SA0P03,2.24575034172e-02,0.0
+0-PSHAB20_INTERFACE-SA0P05,2.58276475982e-02,0.0
+0-PSHAB20_INTERFACE-SA0P075,2.96364966819e-02,0.0
+0-PSHAB20_INTERFACE-SA0P1,3.25369383278e-02,0.0
+0-PSHAB20_INTERFACE-SA0P15,3.55651022874e-02,0.0
+0-PSHAB20_INTERFACE-SA0P2,3.64296484198e-02,0.0
+0-PSHAB20_INTERFACE-SA0P25,3.69487943052e-02,0.0
+0-PSHAB20_INTERFACE-SA0P3,3.71144622014e-02,0.0
+0-PSHAB20_INTERFACE-SA0P4,3.46117876964e-02,0.0
+0-PSHAB20_INTERFACE-SA0P5,3.11738351662e-02,0.0
+0-PSHAB20_INTERFACE-SA0P75,2.56085017514e-02,0.0
+0-PSHAB20_INTERFACE-SA1P0,2.06123859520e-02,0.0
+0-PSHAB20_INTERFACE-SA1P5,1.45673139070e-02,0.0
+0-PSHAB20_INTERFACE-SA2P0,1.08967064936e-02,0.0
+0-PSHAB20_INTERFACE-SA3P0,6.26946254992e-03,0.0
+0-PSHAB20_INTERFACE-SA4P0,4.36498415871e-03,0.0
+0-PSHAB20_INTERFACE-SA5P0,3.44133448941e-03,0.0
+0-PSHAB20_INTERFACE-SA7P5,1.95537566087e-03,0.0
+0-PSHAB20_INTERFACE-SA10P0,1.25157789513e-03,0.0
+0-PSHAB20_INTERFACE-PGA,2.02870001012e-02,0.0
+0-PSHAB20_INTERFACE-PGV,1.95216532021e+00,0.0