Merge branch 'dsf-edits' into 'main'

DSF edits

See merge request !414

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

@@ -2,20 +2,25 @@ package gov.usgs.earthquake.nshmp.gmm;
 
 import static com.google.common.base.Preconditions.checkArgument;
 import static gov.usgs.earthquake.nshmp.data.DoubleData.checkInRange;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.Type.ACTIVE_CRUST;
+import static gov.usgs.earthquake.nshmp.gmm.Gmm.Type.STABLE_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 java.lang.Math.abs;
 import static java.lang.Math.log;
 import static java.lang.Math.sqrt;
+import static java.util.stream.Collectors.toList;
 
+import java.util.Arrays;
 import java.util.EnumMap;
 import java.util.List;
 import java.util.Map;
 import java.util.Optional;
+import java.util.Set;
 
 import com.google.common.annotations.Beta;
-import com.google.common.collect.FluentIterable;
 import com.google.common.collect.Lists;
 import com.google.common.collect.Range;
-import com.google.common.primitives.Doubles;
 
 import gov.usgs.earthquake.nshmp.data.Interpolator;
 import gov.usgs.earthquake.nshmp.data.XySequence;
@@ -32,8 +37,7 @@ import gov.usgs.earthquake.nshmp.tree.Trees;
  *
  * <p><b>Implementation notes:</b>
  *
- * <ul><li>Active crust model uses coefficients for RotD50, Table 1 in Rezaeian
- * et al. (2014)
+ * <ul><li>Active crust model is used for stable crust GMMs.
  *
  * </ul>
  *
@@ -62,9 +66,8 @@ public class UsgsDampingScaling {
    * Current set up is to use instances by tectonic setting.
    *
    * The calcDsf() model returns a value much closer to zero for a damping ratio
-   * of 5.027590868... The assumption is that only coarser damping ratios will
-   * be used (e.g. 3.5, 4.0) such that this inconsistency will not make a
-   * difference. Are we going to support a range of fixed values?
+   * of 5.027590868. By supporting only fixed, coarser damping ratios (e.g. 3,
+   * 4) this inconsistency does not make a difference.
    *
    * For subduction, (1) special handling for periods > 5s, (2) correlation
    * coefficient table missing 5% column (is 0.0 in crustal), added column for
@@ -79,24 +82,26 @@ public class UsgsDampingScaling {
   /** The range of supported damping ratios: {@code [0.5..30]}. */
   public static final Range<Double> DAMPING_RATIO_RANGE = Range.closed(0.05, 30.0);
 
-  private static final CoefficientContainer COEFFS_ACTIVE_CRUST;
-  private static final CoefficientContainer COEFFS_INTERFACE;
-  private static final CoefficientContainer COEFFS_SLAB;
-
-  private static final String[] corrCoeffKeys =
+  private static final String[] DSF_KEYS =
       { "0.5", "1", "2", "3", "5", "7", "10", "15", "20", "25", "30" };
 
-  private static final List<Double> CORR_COEFF_DSFS =
-      FluentIterable.from(corrCoeffKeys)
-          .transform(Doubles.stringConverter())
-          .toList();
+  private static final List<Double> DSF_VALUES = Arrays.stream(DSF_KEYS)
+      .map(Double::valueOf)
+      .collect(toList());
 
   private static final Interpolator CC_INTERP = Interpolator.builder().build();
 
+  private static final Map<Gmm.Type, UsgsDampingScaling> INSTANCE_MAP;
+
   static {
-    COEFFS_ACTIVE_CRUST = new CoefficientContainer("Rezaeian14-active-crust.csv");
-    COEFFS_INTERFACE = new CoefficientContainer("Rezaeian21-interface.csv");
-    COEFFS_SLAB = new CoefficientContainer("Rezaeian21-slab.csv");
+    var dsfActive = new UsgsDampingScaling("Rezaeian14-active-crust.csv");
+    var dsfInterface = new UsgsDampingScaling("Rezaeian21-interface.csv");
+    var dsfSlab = new UsgsDampingScaling("Rezaeian21-slab.csv");
+    INSTANCE_MAP = Map.of(
+        ACTIVE_CRUST, dsfActive,
+        STABLE_CRUST, dsfActive,
+        SUBDUCTION_INTERFACE, dsfInterface,
+        SUBDUCTION_SLAB, dsfSlab);
   }
 
   private static final class Coefficients {
@@ -122,39 +127,8 @@ public class UsgsDampingScaling {
   private final Map<Imt, Coefficients> coeffs = new EnumMap<>(Imt.class);
   private final Map<Imt, XySequence> corrCoeffs = new EnumMap<>(Imt.class);
 
-  // private final SourceType type;
-
-  /**
-   * Create a new damping scaling factor (DSF) model instance for the specified
-   * tectonic settign and source type.
-   *
-   * @param setting tectonic setting for the DSF
-   * @param type source type for the DSF; this argument is only relevant if
-   *        {@code TectonicSetting == SUBDUCTION}
-   */
-  private UsgsDampingScaling(Gmm.Type type) {
-    // Dallin- Modified the function to only use source type to determine the
-    // correct coeffs
-    // this.type = type;
-    CoefficientContainer cc = (type == Gmm.Type.SUBDUCTION_INTERFACE)
-        ? COEFFS_INTERFACE
-        : (type == Gmm.Type.SUBDUCTION_SLAB)
-            ? COEFFS_SLAB
-            : COEFFS_ACTIVE_CRUST;
-
-    // COEFFS_ACTIVE_CRUST;
-    // cc = ()
-    // if (type == Gmm.Type.SUBDUCTION_INTERFACE)
-    // if (type == INTERFACE || type == INTERFACE_CLUSTER || type ==
-    // INTERFACE_SYSTEM ||
-    // type == INTERFACE_GRID) {
-    // cc = COEFFS_INTERFACE;
-    // } else if (type == SLAB || type == SLAB_GRID) {
-    // cc = COEFFS_SLAB;
-    // } else {
-    // cc = COEFFS_ACTIVE_CRUST;
-    // }
-
+  private UsgsDampingScaling(String coeffsFile) {
+    CoefficientContainer cc = new CoefficientContainer(coeffsFile);
     for (Imt imt : cc.imts()) {
       coeffs.put(imt, new Coefficients(imt, cc));
       corrCoeffs.put(imt, initCorrCoeff(imt, cc));
@@ -168,10 +142,15 @@ public class UsgsDampingScaling {
   private static XySequence initCorrCoeff(Imt imt, CoefficientContainer cc) {
     Map<String, Double> coeffs = cc.get(imt);
     List<Double> corrCoeffs = Lists.newArrayList();
-    for (String key : corrCoeffKeys) {
+    for (String key : DSF_KEYS) {
       corrCoeffs.add(coeffs.get(key));
     }
-    return XySequence.create(CORR_COEFF_DSFS, corrCoeffs);
+    return XySequence.create(DSF_VALUES, corrCoeffs);
+  }
+
+  /** The set of IMTs supported by this class. */
+  public static final Set<Imt> supportedImts() {
+    return INSTANCE_MAP.get(ACTIVE_CRUST).coeffs.keySet();
   }
 
   /**
@@ -186,21 +165,9 @@ public class UsgsDampingScaling {
   public double factor(Imt imt, double Mw, double rRup, double dampingRatio) {
     checkArgument(coeffs.containsKey(imt), "Unsupported IMT: %s", imt.name());
     checkInRange(DAMPING_RATIO_RANGE, "Damping ratio", dampingRatio);
-    return (dampingRatio != 5.0) ? calcFactor(imt, Mw, rRup, dampingRatio) : 0;
-    // double lnDSF = 0.0;
-    // if (dampingRatio != 5.0) {
-    // if (imt.period() > 5.0 &&
-    // (type == INTERFACE || type == INTERFACE_CLUSTER || type ==
-    // INTERFACE_SYSTEM ||
-    // type == INTERFACE_GRID || type == SLAB || type == SLAB_GRID)) {
-    // // System.out.println("LONGPERIOD");
-    // lnDSF = calclnDsflongPeriod(imt, Mw, rRup, dampingRatio);
-    // } else {
-    // // System.out.println("normal");
-    // lnDSF = calclnDsf(imt, Mw, rRup, dampingRatio);
-    // }
-    // }
-    // return lnDSF;
+    return (dampingRatio != 5.0)
+        ? calcFactor(imt, Mw, rRup, dampingRatio)
+        : 0;
   }
 
   double calcFactor(Imt imt, double Mw, double rRup, double dampingRatio) {
@@ -281,13 +248,6 @@ public class UsgsDampingScaling {
     return sqrt(σPsa * σPsa + σDsf * σDsf + 2.0 * σPsa * σDsf * ρ);
   }
 
-  private static final UsgsDampingScaling INSTANCE_DEFAULT =
-      new UsgsDampingScaling(Gmm.Type.ACTIVE_CRUST);
-  private static final UsgsDampingScaling INSTANCE_INTERFACE =
-      new UsgsDampingScaling(Gmm.Type.SUBDUCTION_INTERFACE);
-  private static final UsgsDampingScaling INSTANCE_SLAB =
-      new UsgsDampingScaling(Gmm.Type.SUBDUCTION_SLAB);
-
   public static LogicTree<GroundMotion> apply(
       Imt imt,
       GmmInput in,
@@ -299,20 +259,8 @@ public class UsgsDampingScaling {
       return tree;
     }
 
-    /*
-     * Changed method signature
-     *
-     * Switch to set instance to use.
-     *
-     * factor() sigma() and rho() don't depend on mean and sigma values; only be
-     * called once
-     */
-
-    UsgsDampingScaling instance = (type == Gmm.Type.SUBDUCTION_INTERFACE)
-        ? INSTANCE_INTERFACE
-        : (type == Gmm.Type.SUBDUCTION_SLAB)
-            ? INSTANCE_SLAB
-            : INSTANCE_DEFAULT;
+    UsgsDampingScaling instance = INSTANCE_MAP.get(type);
+    checkArgument(instance.coeffs.containsKey(imt), "Unsupported IMT: %s", imt.name());
 
     double factor = instance.calcFactor(imt, in.Mw, in.rRup, dampingRatio);
     double sigma = instance.calcSigma(imt, dampingRatio);
@@ -324,52 +272,17 @@ public class UsgsDampingScaling {
             gm.mean() + factor,
             updateSigma(gm.sigma(), sigma, rho)),
         Optional.empty());
-
-    // if (type == Gmm.Type.SUBDUCTION_INTERFACE) {
-    //
-    // // if (gmmInput.srcType == INTERFACE || gmmInput.srcType ==
-    // // INTERFACE_CLUSTER ||
-    // // gmmInput.srcType == INTERFACE_SYSTEM || gmmInput.srcType ==
-    // // INTERFACE_GRID) {
-    //
-    // return Trees.transform(tree,
-    // gm -> GroundMotion.create(
-    // gm.mean() +
-    // INSTANCE_INTERFACE.factor(imt, gmmInput.Mw, gmmInput.rRup, dampingRatio),
-    // updateSigma(gm.sigma(), INSTANCE_INTERFACE.sigma(imt, dampingRatio),
-    // INSTANCE_INTERFACE.rho(imt, dampingRatio))),
-    // Optional.empty());
-    //
-    // } else if (type == Gmm.Type.SUBDUCTION_SLAB) {
-    // // if (gmmInput.srcType == SLAB || gmmInput.srcType == SLAB_GRID) {
-    // return Trees.transform(tree,
-    // gm -> GroundMotion.create(
-    // gm.mean() +
-    // INSTANCE_SLAB.factor(imt, gmmInput.Mw, gmmInput.rRup, dampingRatio),
-    // updateSigma(gm.sigma(), INSTANCE_SLAB.sigma(imt, dampingRatio),
-    // INSTANCE_SLAB.rho(imt, dampingRatio))),
-    // Optional.empty());
-    //
-    // } else {
-    // // active and stable crust
-    // return Trees.transform(tree,
-    // gm -> GroundMotion.create(
-    // gm.mean() +
-    // INSTANCE_DEFAULT.factor(imt, gmmInput.Mw, gmmInput.rRup, dampingRatio),
-    // updateSigma(gm.sigma(), INSTANCE_DEFAULT.sigma(imt, dampingRatio),
-    // INSTANCE_DEFAULT.rho(imt, dampingRatio))),
-    // Optional.empty());
-    // }
   }
 
+  // TODO clean
   public static void main(String[] args) {
 
-    UsgsDampingScaling model = new UsgsDampingScaling(Gmm.Type.SUBDUCTION_INTERFACE);
+    UsgsDampingScaling model = INSTANCE_MAP.get(SUBDUCTION_INTERFACE);
 
-    Imt imt = Imt.SA5P0;
+    Imt imt = Imt.SA1P0;
     double Mw = 7.0;
     double rRup = 120;
-    double ratio = 5.0;
+    double ratio = 3.0;
 
     double scaleFactor = Math.exp(model.factor(imt, Mw, rRup, ratio));
     double dsfSigma = model.sigma(imt, ratio);