tapered GR mfd temp fix

src/main/java/gov/usgs/earthquake/nshmp/mfd/Mfd.java

@@ -27,6 +27,7 @@ import gov.usgs.earthquake.nshmp.data.MutableXySequence;
 import gov.usgs.earthquake.nshmp.data.Sequences;
 import gov.usgs.earthquake.nshmp.data.XyPoint;
 import gov.usgs.earthquake.nshmp.data.XySequence;
+import gov.usgs.earthquake.nshmp.mfd.Mfd.Properties.TaperedGr;
 
 /**
  * Entry point for creating a magnitude frequency distribution (MFD). An MFD
@@ -429,6 +430,10 @@ public final class Mfd {
      * @return this {@code Builder} object
      */
     public Builder scaleToIncrementalRate(double incrementalRate) {
+      // TODO remove (needs sensitivity testing)
+      if (props.type == Type.GR_TAPER) {
+        return TaperedGr.taperCorrection(this, incrementalRate);
+      }
       return scale(checkRate(incrementalRate) / mfd.y(0));
     }
 
@@ -453,6 +458,7 @@ public final class Mfd {
       mfd.transform(checkNotNull(action));
       return this;
     }
+
   }
 
   /**
@@ -779,6 +785,13 @@ public final class Mfd {
          * directly with a given a-value and by creating a generic taperd MFD
          * that we then scaleToIncrementalRate() using the given a-value, we
          * correct the tapered MFD by the offset.
+         *
+         * Now, this is how we wished the nshmp-haz implementation behaved, but
+         * since it doesn't we need an intermediate solution to match nshmp-haz
+         * until further sensitivity testing through to hazard can be
+         * undertaken. What we want: after calling scaleToIncrementalRate(), the
+         * tapered GR MFD rates should match those that would've resulted from
+         * the direct calculation.
          */
         double grRate = Mfds.gutenbergRichterRate(a(), b(), builder.mfd.x(0));
         double tgrRate = builder.mfd.y(0);
@@ -798,6 +811,24 @@ public final class Mfd {
             type(),
             super.propsString() + ", mc=" + mc);
       }
+
+      /*
+       * TODO clean after sensitivity testing for tapered rates when calling
+       * scale toIncrementalRate() tapered GR values become exact matches to the
+       * target rate as opposed to the slightly scaled rate resulting from
+       * creating the tapered GR directly with a given a-value. We want to
+       * replicate nshmp-haz behavior for now but get rid of it in the not too
+       * distant future. This is basically reverting the scaling correction
+       * applied in the builder.
+       */
+      private static Builder taperCorrection(Builder b, double rate) {
+        TaperedGr p = b.props.getAsGrTaper();
+        TaperFunction taperFn = new TaperFunction(p.Δm(), p.b(), p.mc);
+        double mMin = b.mfd.x(0);
+        double taperRate = checkRate(rate) * taperFn.scale(mMin);
+        return b.scale(taperRate / b.mfd.y(0));
+      }
+
     }
 
     private static final double M_MIN_MOMENT = Earthquakes.magToMoment(4.0);

src/test/java/gov/usgs/earthquake/nshmp/mfd/MfdTests.java

@@ -10,6 +10,7 @@ import java.io.IOException;
 import java.nio.file.Files;
 import java.nio.file.Path;
 import java.nio.file.Paths;
+import java.util.Arrays;
 
 import org.junit.jupiter.api.Test;
 
@@ -202,6 +203,66 @@ class MfdTests {
         grProps.toString());
   }
 
+  /*
+   * TODO clean
+   *
+   * TODO something in the TaperedGR moment based scaling leads to differences
+   * in the 4th to 5th significant figure that should be better understood, and
+   * then this main method cleaned out. The results are still pretty close
+   * overall.
+   */
+  public static void main(String[] args) {
+    double bVal = 0.8;
+    double aVal = Math.log10(4.0);
+    double Δm = 0.1;
+    double mMin = 5.0;
+    double mMax = 8.0;
+    double mc = 7.5;
+
+    // /* GR direct vs scaled builder comparison */
+    // System.out.println("Gutenberg Richter");
+    // Mfd mfd = new GutenbergRichter(aVal, bVal, Δm, mMin, mMax)
+    // .toBuilder()
+    // .build();
+    //
+    // System.out.println(Arrays.toString(mfd.data().xValues().toArray()));
+    // System.out.println(Arrays.toString(mfd.data().yValues().toArray()));
+    //
+    // double incrRate = Mfds.gutenbergRichterRate(aVal, bVal, 5.05);
+    // mfd = Mfd.newGutenbergRichterBuilder(bVal, Δm, mMin, mMax)
+    // .scaleToIncrementalRate(incrRate)
+    // .build();
+    //
+    // System.out.println(Arrays.toString(mfd.data().yValues().toArray()));
+
+    /* GR direct vs scaled builder comparison */
+    System.out.println("Tapered Gutenberg Richter");
+    Mfd mfd = new TaperedGr(aVal, bVal, Δm, mMin, mMax, mc)
+        .toBuilder()
+        .build();
+
+    double[] directRates = mfd.data().yValues().toArray();
+    System.out.println(Arrays.toString(mfd.data().xValues().toArray()));
+    System.out.println(Arrays.toString(directRates));
+
+    double incrRate = Mfds.gutenbergRichterRate(aVal, bVal, 5.05);
+    mfd = Mfd.newTaperedGutenbergRichterBuilder(bVal, Δm, mMin, mMax, mc)
+        .scaleToIncrementalRate(incrRate)
+        .build();
+
+    double[] scaledRates = mfd.data().yValues().toArray();
+    System.out.println(Arrays.toString(scaledRates));
+
+    int size = directRates.length;
+    double[] ratios = new double[size];
+    for (int i = 0; i < size; i++) {
+      double ratio = 100.0 * (directRates[i] / scaledRates[i]);
+      ratios[i] = ratio;
+    }
+    System.out.println(Arrays.toString(ratios));
+
+  }
+
   private static final double[] TAPERED_GR_M = {
       5.05, 5.15, 5.25, 5.35, 5.45,
       5.55, 5.65, 5.75, 5.85, 5.95,
@@ -243,6 +304,39 @@ class MfdTests {
       5.437829885686743E-7,
       1.5949675112240326E-7 };
 
+  /* TODO These match nshmp-haz and should be removed */
+  private static final double[] TAPERED_GR_R_TEMP = {
+      3.648734771264746E-4,
+      3.0351155247723867E-4,
+      2.5247694248332143E-4,
+      2.1003291504182058E-4,
+      1.7473503715378214E-4,
+      1.4538201763727166E-4,
+      1.2097482132130438E-4,
+      1.0068266229609009E-4,
+      8.381471718000333E-5,
+      6.979659287661517E-5,
+      5.815074326255645E-5,
+      4.848016071724236E-5,
+      4.0454775273297684E-5,
+      3.3800079282566184E-5,
+      2.8287560844165223E-5,
+      2.3726576420233305E-5,
+      1.9957324528955872E-5,
+      1.6844604598702632E-5,
+      1.4272075425536471E-5,
+      1.213680849238659E-5,
+      1.0344152445466782E-5,
+      8.803445832444013E-6,
+      7.426137715686031E-6,
+      6.1294171417451564E-6,
+      4.849632254896957E-6,
+      3.5675342487878286E-6,
+      2.3362943546550987E-6,
+      1.2825840184413834E-6,
+      5.438860517858617E-7,
+      1.5952698054966957E-7 };
+
   @Test
   void testTaperedGr() {
 
@@ -262,7 +356,7 @@ class MfdTests {
     assertTrue(xy.size() == 30);
     assertArrayEquals(TAPERED_GR_M, xy.xValues().toArray());
     /* These end up being very slightly different due to scaling */
-    assertArrayEquals(TAPERED_GR_R, xy.yValues().toArray(), 1e-17);
+    assertArrayEquals(TAPERED_GR_R_TEMP, xy.yValues().toArray(), 1e-17);
 
     /* Properties */
     Properties props = new TaperedGr(aVal, bVal, Δm, mMin, mMax, mc);