diff --git a/src/main/java/gov/usgs/earthquake/nshmp/Faults.java b/src/main/java/gov/usgs/earthquake/nshmp/Faults.java
index 48122228992ba0e08201792f6de9675aad07d2b2..d1ede73fe9bb49465691b638008a5a5c42614425 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/Faults.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/Faults.java
@@ -138,15 +138,14 @@ public final class Faults {
return zTor + Math.sin(dip * Maths.TO_RADIANS) * width / 2.0;
}
- public static double width(double dip, double zUpper, double zLower) {
- return sinProjection(dip, zLower - zUpper);
- }
-
- public static double projectVerticalSlip(double dip, double slip) {
- return sinProjection(dip, slip);
- }
-
- private static double sinProjection(double dip, double value) {
+ /**
+ * Project a vertically measured value (e.g. slip rate, depth) onto a dipping
+ * surface.
+ *
+ * @param dip of surface (i.e. a fault)
+ * @param value to project (vertical)
+ */
+ public static double projectVertical(double dip, double value) {
return value / Math.sin(dip * Maths.TO_RADIANS);
}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/GmmInput.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/GmmInput.java
index d7b443f286da917e12d8b0ee2e9963a80fe368ea..e612c0c3b90eab5806259cb125acfdd3519011bc 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/GmmInput.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/GmmInput.java
@@ -568,14 +568,14 @@ public class GmmInput {
@Override
public String toString() {
return MoreObjects.toStringHelper(this)
- .add(MW.toString(), String.format("%.2f", Mw))
- .add(RJB.toString(), String.format("%.3f", rJB))
- .add(RRUP.toString(), String.format("%.3f", rRup))
- .add(RX.toString(), String.format("%.3f", rX))
+ .add(MW.toString(), String.format("%.4f", Mw))
+ .add(RJB.toString(), String.format("%.6f", rJB))
+ .add(RRUP.toString(), String.format("%.6f", rRup))
+ .add(RX.toString(), String.format("%.6f", rX))
.add(DIP.toString(), dip)
- .add(WIDTH.toString(), String.format("%.3f", width))
+ .add(WIDTH.toString(), String.format("%.6f", width))
.add(ZTOP.toString(), zTop)
- .add(ZHYP.toString(), String.format("%.3f", zHyp))
+ .add(ZHYP.toString(), String.format("%.6f", zHyp))
.add(RAKE.toString(), rake)
.add(VS30.toString(), vs30)
.add(VSINF.toString(), vsInf)
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/mfd/Mfd.java b/src/main/java/gov/usgs/earthquake/nshmp/mfd/Mfd.java
index 3378ad66b672cd06e7fac9af1d9fb73cb330513c..4696a965ac8cf107bf9b7d5dd22bd90d1dbbbae0 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/mfd/Mfd.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/mfd/Mfd.java
@@ -650,15 +650,27 @@ public final class Mfd {
}
/**
- * Create a new properties object for a Gutenberg-Richter MFD with a copy
- * of the supplied properties and an updated maximum magnitude.
+ * Create a new Gutenberg-Richter MFD properties object from a copy of the
+ * supplied properties and an updated maximum magnitude.
*
* @param gr properties to copy
- * @param mMax updated maximum magnitude
+ * @param mMax maximum magnitude
*/
- public GutenbergRichter(GutenbergRichter gr, double mMax) {
- /* Supports common magnitude arrays for variable mMax MFDs */
- this(gr.type(), gr.a, gr.b, gr.Δm, gr.mMin, mMax);
+ public static GutenbergRichter copyOfWithMax(GutenbergRichter gr, double mMax) {
+ /* Supports common magnitude arrays for variable mMax MFDs. */
+ return new GutenbergRichter(gr.type(), gr.a, gr.b, gr.Δm, gr.mMin, mMax);
+ }
+
+ /**
+ * Create a new Gutenberg-Richter MFD properties object from a copy of the
+ * supplied properties and an updated a-value.
+ *
+ * @param gr properties to copy
+ * @param a value
+ */
+ public static GutenbergRichter copyOfWithRate(GutenbergRichter gr, double a) {
+ /* Supports prescriptive a-values when mo-rate balancing. */
+ return new GutenbergRichter(gr.type(), a, gr.b, gr.Δm, gr.mMin, gr.mMax);
}
/** The Gutenberg-Richter {@code a}-value. */
@@ -699,7 +711,7 @@ public final class Mfd {
double mMin = checkMagnitude(this.mMin);
double mMax = checkMagnitude(this.mMax);
double[] magnitudes = Sequences.arrayBuilder(mMin, mMax, this.Δm)
- .scale(4)
+ .scale(5)
.build();
Builder builder = new Builder(this, magnitudes);
builder.transform(p -> p.set(Mfds.gutenbergRichterRate(a, b, p.x())));
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/Deserialize.java b/src/main/java/gov/usgs/earthquake/nshmp/model/Deserialize.java
index 9f3243c56b468bad872723f834abc4b9449b25bf..84b8b0145bfd73757ed0b44a82085eb347aecc2e 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/Deserialize.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/Deserialize.java
@@ -359,6 +359,11 @@ class Deserialize {
return clusterSet.build();
}
+ /* Convenience method to use Deserialize.GSON instance external to class. */
+ static <T> T fromJson(JsonElement json, Type typeOfT) {
+ return GSON.fromJson(json, typeOfT);
+ }
+
/* Create an mfd-tree from a JSON element or via map lookup. */
private static LogicTree<Mfd.Properties> mfdTree(
JsonElement mfdElement,
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/FaultRuptureSet.java b/src/main/java/gov/usgs/earthquake/nshmp/model/FaultRuptureSet.java
index 28fa26e0f32cf64650b45465c7a9172cdeec5dda..f2d829c0b8a7f3a5e591c6f1796b447435a7a1f1 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/FaultRuptureSet.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/FaultRuptureSet.java
@@ -14,6 +14,7 @@ import static java.util.stream.Collectors.toList;
import java.util.List;
import java.util.Map;
import java.util.Optional;
+import java.util.OptionalDouble;
import java.util.stream.Stream;
import gov.usgs.earthquake.nshmp.Earthquakes;
@@ -27,9 +28,12 @@ import gov.usgs.earthquake.nshmp.geo.Locations;
import gov.usgs.earthquake.nshmp.geo.json.Feature;
import gov.usgs.earthquake.nshmp.geo.json.Properties;
import gov.usgs.earthquake.nshmp.mfd.Mfd;
+import gov.usgs.earthquake.nshmp.mfd.Mfd.Properties.GutenbergRichter;
import gov.usgs.earthquake.nshmp.mfd.Mfd.Properties.Single;
import gov.usgs.earthquake.nshmp.mfd.Mfds;
+import gov.usgs.earthquake.nshmp.model.MfdConfig.AleatoryProperties;
import gov.usgs.earthquake.nshmp.model.SourceFeature.Key;
+import gov.usgs.earthquake.nshmp.model.SourceFeature.NshmFault.Rates;
import gov.usgs.earthquake.nshmp.tree.Branch;
import gov.usgs.earthquake.nshmp.tree.LogicTree;
@@ -210,10 +214,12 @@ public class FaultRuptureSet implements RuptureSet {
/*
* Fault MFDs are constructed in a variety of ways. Single fault MFDs may
* be created from a logic tree of slip rates.
+ *
+ *
*/
MfdConfig mfdConfig = data.mfdConfig().orElseThrow();
- Optional<LogicTree<Double>> moRateTree = Optional.empty();
+ Optional<LogicTree<MoRate>> moRateTree = Optional.empty();
/* Build mfd-tree from fault slip or recurrence rates. */
if (data.faultRateTree().isPresent()) {
@@ -221,9 +227,7 @@ public class FaultRuptureSet implements RuptureSet {
// System.out.println(data.faultRateTree().get());
/* Create Mo-rate tree. */
- moRateTree = Optional.of(createMoRateTree(
- feature,
- data.faultRateTree().orElseThrow()));
+ moRateTree = Optional.of(createMoRateTree(feature, data));
/* Update MFD properties tree based on source magnitude. */
mfdPropsTree = updateMfdsForSlip();
@@ -261,10 +265,7 @@ public class FaultRuptureSet implements RuptureSet {
.trace(feature.trace)
.depth(feature.upperDepth)
.dip(feature.dip)
- .width(Faults.width(
- feature.dip,
- feature.upperDepth,
- feature.lowerDepth))
+ .width(feature.width)
.spacing(config.surfaceSpacing)
.build();
@@ -436,12 +437,12 @@ public class FaultRuptureSet implements RuptureSet {
/* Logic tree of moment rates. */
private static LogicTree<Mfd> createMfdTree(
LogicTree<Mfd.Properties> mfdPropsTree,
- Optional<LogicTree<Double>> moRateTree,
+ Optional<LogicTree<MoRate>> moRateTree,
MfdConfig mfdConfig) {
LogicTree.Builder<Mfd> mfdTree = LogicTree.builder(MFD_TREE);
if (moRateTree.isPresent()) {
- LogicTree<Double> moTree = moRateTree.orElseThrow();
+ LogicTree<MoRate> moTree = moRateTree.orElseThrow();
moTree.forEach(branch -> addMfds(mfdTree, mfdPropsTree, Optional.of(branch), mfdConfig));
} else {
addMfds(mfdTree, mfdPropsTree, Optional.empty(), mfdConfig);
@@ -452,16 +453,31 @@ public class FaultRuptureSet implements RuptureSet {
private static void addMfds(
LogicTree.Builder<Mfd> mfdTree, // receiver tree
LogicTree<Mfd.Properties> mfdPropsTree,
- Optional<Branch<Double>> moBranch,
+ Optional<Branch<MoRate>> moBranch,
MfdConfig mfdConfig) {
- for (Branch<Mfd.Properties> mfdBranch : mfdPropsTree) {
+ /*
+ * Working with branch list because weights may not sum to 1.0 as we loop
+ * over moment rate branches to build receiver mfdTree. It's simpler to
+ * update branch names and weights here than do it in both Single and GR MFD
+ * factories below.
+ *
+ * Update properties tree with mo rate branch name and weight.
+ */
+ List<Branch<Mfd.Properties>> propsBranchList = moBranch.isPresent()
+ ? MfdTrees.toBranchList(mfdPropsTree, moBranch.orElseThrow())
+ : mfdPropsTree;
+ Optional<MoRate> moRateIn = moBranch.map(Branch::value);
+
+ for (Branch<Mfd.Properties> mfdBranch : propsBranchList) {
+ double moRate = computeMoRate(mfdBranch.value(), moRateIn);
+
switch (mfdBranch.value().type()) {
case SINGLE:
- addSingleMfds(mfdTree, mfdBranch, moBranch, mfdConfig);
+ addSingleMfds(mfdTree, mfdBranch, moRate, mfdConfig);
break;
case GR:
- addGrMfds(mfdTree, mfdBranch, moBranch, mfdConfig);
+ addGrMfds(mfdTree, mfdBranch, moRate, mfdConfig);
break;
default:
throw new UnsupportedOperationException(
@@ -469,37 +485,108 @@ public class FaultRuptureSet implements RuptureSet {
}
}
}
+ //
+
+ /*
+ * Container class for slip-rate based moRate.
+ *
+ * Optional fields, if present, are used downstream to override slip-based
+ * moRates to yield better agreement with legacy codes.
+ */
+ private static class MoRate {
+
+ final double moRate;
+ final OptionalDouble graValue;
+ final OptionalDouble singleRate;
+
+ MoRate(
+ double moRate,
+ OptionalDouble graValue,
+ OptionalDouble singleRate) {
+
+ this.moRate = moRate;
+ this.graValue = graValue;
+ this.singleRate = singleRate;
+ }
+ }
+
+ private static double computeMoRate(Mfd.Properties props, Optional<MoRate> moRateData) {
+
+ /* No mo-rate data supplied, compute from properties */
+ if (moRateData.isEmpty()) {
+ switch (props.type()) {
+ case SINGLE:
+ Mfd.Properties.Single single = props.getAsSingle();
+ return single.rate() * Earthquakes.magToMoment(single.magnitude());
+ case GR:
+ Mfd.Properties.GutenbergRichter gr = props.getAsGr();
+ return grMomentRate(gr);
+ default:
+ throw new UnsupportedOperationException(
+ props.type() + " MFD not supported");
+ }
+ }
+
+ /* Return supplied mo-rate or compute from updated properties. */
+ MoRate rateData = moRateData.orElseThrow();
+ switch (props.type()) {
+ case SINGLE:
+ if (rateData.singleRate.isPresent()) {
+ return rateData.singleRate.orElseThrow() *
+ Earthquakes.magToMoment(props.getAsSingle().magnitude());
+ }
+ return rateData.moRate;
+ case GR:
+ if (rateData.graValue.isPresent()) {
+ Mfd.Properties.GutenbergRichter grNew =
+ GutenbergRichter.copyOfWithRate(props.getAsGr(), rateData.graValue.orElseThrow());
+ return grMomentRate(grNew);
+ }
+ return rateData.moRate;
+ default:
+ throw new UnsupportedOperationException(
+ props.type() + " MFD not supported");
+ }
+ }
/*
* Create a logic tree of moment rates for a fault source given a feature and
* rate data. If present, a rate-dip may be used to compute a constraining
* moment rate (e.g. WUS normal faults)
*/
- private static LogicTree<Double> createMoRateTree(
+ private static LogicTree<MoRate> createMoRateTree(
SourceFeature.NshmFault feature,
- LogicTree<String> faultRateTree) {
+ ModelData data) {
+
+ double dip = feature.dip;
+ double width = feature.width;
+
+ // TODO this needs revisiting for 2008 and/or
+ // future models if slip rate is to be used
+ // rateDip should not be present for 2014/2018 models
- double dip = feature.rateDip.isPresent()
- ? feature.rateDip.getAsDouble()
- : feature.dip;
- double width = Faults.width(dip, feature.upperDepth, feature.lowerDepth);
- double length = feature.length.orElseThrow();
+ /* Dip slip model dependent parameters. */
+ if (feature.rateDip.isPresent()) {
+ dip = feature.rateDip.orElseThrow();
+ double height = feature.lowerDepth - feature.upperDepth;
+ width = Faults.projectVertical(dip, height);
+ }
+
+ double length = feature.length;
double area = length * width * 1e6;
- double Mw = feature.magnitude.getAsDouble();
+ double Mw = feature.magnitude.orElseThrow();
double mo = Earthquakes.magToMoment(Mw);
- Map<String, Double> rateMap = feature.rateMap;
+ Map<String, Rates> rateMap = feature.rateMap;
RateType rateType = feature.rateType;
- LogicTree.Builder<Double> moRateTree = LogicTree.builder(RATE_TREE);
+ LogicTree.Builder<MoRate> moRateTree = LogicTree.builder(RATE_TREE);
// TODO clean
- // System.out.println(feature.name);
- // System.out.println(rateMap);
if (rateMap.containsValue(null)) {
System.out.println(" " + rateMap);
}
- for (Branch<String> branch : faultRateTree) {
+ for (Branch<String> branch : data.faultRateTree().orElseThrow()) {
/*
* TODO Verify missing fault rates and imbalanced logic trees. In some
@@ -509,30 +596,21 @@ public class FaultRuptureSet implements RuptureSet {
* because GEO weight is 0.8 and imbalanced; maybe we just handle this
* with setting rate to 0, as below.
*/
- Double rateCheck = rateMap.get(branch.id());
- double rate = (rateCheck != null) ? rateCheck : 0.0;
+ Rates rates = rateMap.get(branch.id());
+
+ double rate = (rates != null) ? rates.rate : 0.0;
- /* CONUS: 2008 are vertical, 2014 are on-plane. */
if (rateType == RateType.VERTICAL_SLIP) {
- rate = Faults.projectVerticalSlip(dip, rate);
+ rate = Faults.projectVertical(dip, rate);
}
- double moRate = (rateType == RECURRENCE)
+ double moRateValue = (rateType == RECURRENCE)
? mo * rate
: Earthquakes.moment(area, rate / 1000.0);
- /*
- * 2014 model consistency. Convert to annual rate and set scale to seven
- * decimal places to match the MFD values in fortran and XML input files
- * and then convert back to moment rate. The steps of going from slip rate
- * to explicit MFD parameters and then back to build MFDs is no longer
- * necessary, but the cumulative effects of rounding at various steps in
- * legacy codes are difficult to replicate without doing things like this.
- * TODO consider describing different rounding/precision characteristics
- * in a model defined in a source config.
- */
- double eqRate = Maths.round(moRate / mo, 7);
- moRate = eqRate * mo;
+ MoRate moRate = (rates != null)
+ ? new MoRate(moRateValue, rates.graValue, rates.singleRate)
+ : new MoRate(moRateValue, OptionalDouble.empty(), OptionalDouble.empty());
moRateTree.addBranch(branch.id(), moRate, branch.weight());
}
return moRateTree.build();
@@ -541,8 +619,8 @@ public class FaultRuptureSet implements RuptureSet {
/* Fault GR */
private static void addGrMfds(
LogicTree.Builder<Mfd> mfdTree, // receiver tree
- Branch<Mfd.Properties> mfdBranch, // properties may include rate
- Optional<Branch<Double>> moBranch, // or rate supplied separately
+ Branch<Mfd.Properties> mfdBranch,
+ double moRate,
MfdConfig mfdConfig) {
Mfd.Properties.GutenbergRichter gr = mfdBranch.value().getAsGr();
@@ -550,19 +628,10 @@ public class FaultRuptureSet implements RuptureSet {
String mfdId = mfdBranch.id();
double mfdWt = mfdBranch.weight();
- if (moBranch.isPresent()) {
- mfdId = mfdBranchId(moBranch.orElseThrow().id(), mfdId);
- mfdWt *= moBranch.orElseThrow().weight();
- }
-
// TODO does this ever throw?
int nMag = magCount(gr.mMin(), gr.mMax(), gr.Δm());
checkState(nMag > 0, "GR MFD with no magnitudes [%s]", mfdBranch.id());
- double moRate = moBranch.isPresent()
- ? moBranch.orElseThrow().value()
- : grMomentRate(gr);
-
/*
* Optional epistemic uncertainty; no aleatory with GR. This allows
* epistemic uncertainty because the cutoff is based on mMin (6.5), but in
@@ -606,10 +675,15 @@ public class FaultRuptureSet implements RuptureSet {
System.out.println("GR MFD epistemic branch with no magnitudes");
}
- String id = mfdBranchId(mfdId, epiBranch.id());
- Mfd mfd = newGrBuilder(gr, mMaxEpi)
+ String id = MfdTrees.combineBranchIds(mfdId, epiBranch.id());
+ Mfd mfd = GutenbergRichter.copyOfWithMax(gr, mMaxEpi)
+ .toBuilder()
.scaleToMomentRate(epiOk ? moRate : 0)
.build();
+
+ // Mfd mfd = newGrBuilder(gr, mMaxEpi)
+ // .scaleToMomentRate(epiOk ? moRate : 0)
+ // .build();
mfdTree.addBranch(id, mfd, weightEpi);
}
} else {
@@ -646,80 +720,98 @@ public class FaultRuptureSet implements RuptureSet {
/* Fault SINGLE */
private static void addSingleMfds(
LogicTree.Builder<Mfd> mfdTree, // receiver tree
- Branch<Mfd.Properties> mfdBranch, // properties may include rate
- Optional<Branch<Double>> moBranch, // or rate supplied separately
+ Branch<Mfd.Properties> mfdBranch,
+ double moRate,
MfdConfig mfdConfig) {
- Mfd.Properties.Single single = mfdBranch.value().getAsSingle();
-
- String mfdId = mfdBranch.id();
- double mfdWt = mfdBranch.weight();
-
- if (moBranch.isPresent()) {
- mfdId = mfdBranchId(moBranch.orElseThrow().id(), mfdId);
- mfdWt *= moBranch.orElseThrow().weight();
- }
+ double magnitude = mfdBranch.value().getAsSingle().magnitude();
- double moRate = moBranch.isPresent()
- ? moBranch.orElseThrow().value()
- : single.rate() * Earthquakes.magToMoment(single.magnitude());
+ // String mfdId = mfdBranch.id(); // TODO clean
+ // double mfdWt = mfdBranch.weight();
/* Optional epistemic uncertainty. */
- boolean epistemic = hasEpistemic(mfdConfig, single.magnitude());
+ boolean epistemic = hasEpistemic(mfdConfig, magnitude);
/* Optional aleatory variability. */
- boolean aleatory = mfdConfig.aleatoryProperties.isPresent();
+ Optional<AleatoryProperties> aleaProps = mfdConfig.aleatoryProperties;
/* mMax epistemic uncertainty branches */
if (epistemic) {
for (Branch<Double> epiBranch : mfdConfig.epistemicTree.orElseThrow()) {
- double mEpi = single.magnitude() + epiBranch.value();
- double weightEpi = mfdWt * epiBranch.weight();
- String id = mfdBranchId(mfdId, epiBranch.id());
+ double mEpi = magnitude + epiBranch.value();
+ double wtEpi = mfdBranch.weight() * epiBranch.weight();
+ String id = MfdTrees.combineBranchIds(mfdBranch.id(), epiBranch.id());
/* Possibly apply aleatory variability */
- if (aleatory) {
-
- MfdConfig.AleatoryProperties aleaProps = mfdConfig.aleatoryProperties.orElseThrow();
- Mfd mfd = newGaussianBuilder(mEpi, aleaProps)
- .scaleToMomentRate(moRate)
- .build();
- mfdTree.addBranch(id, mfd, weightEpi);
+ // Mfd.Builder builder = (aleaProps.isPresent())
+ // ? newGaussianBuilder(mEpi, aleaProps.orElseThrow())
+ // : Mfd.newSingleBuilder(mEpi);
+ // Mfd mfd = builder.scaleToMomentRate(moRate).build();
- } else {
+ Mfd mfd = createSingleMfd(mEpi, moRate, aleaProps);
+ mfdTree.addBranch(id, mfd, wtEpi);
- Mfd mfd = Mfd.newSingleBuilder(mEpi)
- .scaleToMomentRate(moRate)
- .build();
- mfdTree.addBranch(id, mfd, weightEpi);
- }
+ /* Possibly apply aleatory variability */
+ // TODO clean
+ // if (aleatory) {
+ //
+ // MfdConfig.AleatoryProperties aleaProps =
+ // mfdConfig.aleatoryProperties.orElseThrow();
+ // Mfd mfd = newGaussianBuilder(mEpi, aleaProps)
+ // .scaleToMomentRate(moRate)
+ // .build();
+ // mfdTree.addBranch(id, mfd, weightEpi);
+ //
+ // } else {
+ //
+ // Mfd mfd = Mfd.newSingleBuilder(mEpi)
+ // .scaleToMomentRate(moRate)
+ // .build();
+ // mfdTree.addBranch(id, mfd, weightEpi);
+ // }
}
} else { /* No epistemic uncertainty */
- /* Possibly apply aleatory variability */
- if (aleatory) {
-
- MfdConfig.AleatoryProperties aleaProps = mfdConfig.aleatoryProperties.orElseThrow();
- Mfd mfd = newGaussianBuilder(single.magnitude(), aleaProps)
- .scaleToMomentRate(moRate)
- .build();
- mfdTree.addBranch(mfdId, mfd, mfdWt);
+ // Mfd.Builder builder = (aleaProps.isPresent())
+ // ? newGaussianBuilder(magnitude, aleaProps.orElseThrow())
+ // : Mfd.newSingleBuilder(magnitude);
+ // Mfd mfd = builder.scaleToMomentRate(moRate).build();
- } else {
+ Mfd mfd = createSingleMfd(magnitude, moRate, aleaProps);
+ mfdTree.addBranch(mfdBranch.id(), mfd, mfdBranch.weight());
- Mfd mfd = Mfd.newSingleBuilder(single.magnitude())
- .scaleToMomentRate(moRate)
- .build();
- mfdTree.addBranch(mfdId, mfd, mfdWt);
- }
+ /* Possibly apply aleatory variability */
+ // TODO clean
+ // if (aleaProps.isPresent()) {
+ //
+ // MfdConfig.AleatoryProperties aleaProps =
+ // mfdConfig.aleatoryProperties.orElseThrow();
+ // Mfd mfd = newGaussianBuilder(magnitude, aleaProps)
+ // .scaleToMomentRate(moRate)
+ // .build();
+ // mfdTree.addBranch(mfdId, mfd, mfdWt);
+ //
+ // } else {
+ //
+ // Mfd mfd = Mfd.newSingleBuilder(magnitude)
+ // .scaleToMomentRate(moRate)
+ // .build();
+ // mfdTree.addBranch(mfdId, mfd, mfdWt);
+ // }
}
}
- private static String mfdBranchId(String... ids) {
- return String.join(" : ", ids);
+ private static Mfd createSingleMfd(
+ double magnitude, double moRate,
+ Optional<AleatoryProperties> aleaProps) {
+
+ Mfd.Builder builder = (aleaProps.isPresent())
+ ? newGaussianBuilder(magnitude, aleaProps.orElseThrow())
+ : Mfd.newSingleBuilder(magnitude);
+ return builder.scaleToMomentRate(moRate).build();
}
private static boolean hasEpistemic(MfdConfig mfdConfig, double magnitude) {
@@ -727,8 +819,8 @@ public class FaultRuptureSet implements RuptureSet {
/* No epistemic uncertainty if: m - Δm < mfdConfig.minMagnitude */
if (epistemic) {
- double mEpiCutoff = mfdConfig.minMagnitude.getAsDouble();
- double mEpiOffset = mfdConfig.minEpiOffset.getAsDouble();
+ double mEpiCutoff = mfdConfig.minMagnitude.orElseThrow();
+ double mEpiOffset = mfdConfig.minEpiOffset.orElseThrow();
epistemic = (magnitude + mEpiOffset) > mEpiCutoff;
}
return epistemic;
@@ -767,11 +859,6 @@ public class FaultRuptureSet implements RuptureSet {
aleaProps.nσ);
}
- /* GR MFD builder with mMax override. */
- static Mfd.Builder newGrBuilder(Mfd.Properties.GutenbergRichter gr, double mMax) {
- return Mfd.newGutenbergRichterBuilder(gr.b(), gr.Δm(), gr.mMin(), mMax);
- }
-
/* Expected moment rate of a GR MFD. */
static double grMomentRate(Mfd.Properties.GutenbergRichter gr) {
double mMin = gr.mMin() + gr.Δm() / 2;
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/FaultSource.java b/src/main/java/gov/usgs/earthquake/nshmp/model/FaultSource.java
index c1047556f4e2ca6ed47b13b035159c4b0b4090c2..acfa83cf6d513f90420cd0e7010f2da72162df33 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/FaultSource.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/FaultSource.java
@@ -103,6 +103,15 @@ public class FaultSource implements Source {
checkState(Iterables.size(Iterables.concat(ruptureLists)) > 0,
"FaultSource has no ruptures");
+ // System.out.println(this);
+
+ // System.out.println(scaledMfds);
+ // long nRuptures = ruptureLists.stream().flatMap(List::stream).count();
+ // System.out.println(nRuptures);
+
+ // System.out.println(Mfds.momentRate(scaledMfds.get(0)));
+ // System.out.println(trace);
+ // System.out.println(surface);
}
@Override
@@ -208,7 +217,10 @@ public class FaultSource implements Source {
// Mfd.Type type = props.type();
double length = trace.length();
- // boolean floats = (type == Type.GR) || floatRuptures;
+ Mfd.Type mfdType = mfd.properties().type();
+
+ boolean floatSingle = mfdType == Mfd.Type.SINGLE && floatSingleRuptures;
+ boolean floats = (mfdType == Mfd.Type.GR) || floatSingle;
for (XyPoint xy : mfd.data()) {
double mag = xy.x();
@@ -231,7 +243,7 @@ public class FaultSource implements Source {
// TODO we want to get the 'floats' attribute out of MFDs
// the only reason it is there is to allow SINGLE to flip-flop
// it should just be a SourceProperty
- if (floatSingleRuptures && msrLength < length) {
+ if (floats && msrLength < length) {
// TODO change to using RuptureFloating.OFF
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java b/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java
index 0751493ae567f1ffa5484db8e48e20962e8c6884..c4738f2bc7ffc2d5ed397fd71dda920af210975e 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java
@@ -284,7 +284,7 @@ class GridLoader {
/* Create total MFDs and transpose mfd tree lists. */
for (FeatureData featureData : dataMap.values()) {
featureData.mfds = reduceMfds(featureData.mfdTrees);
- featureData.mfdsTree = MfdTrees.transposeTree(featureData.mfdTrees);
+ featureData.mfdsTree = MfdTrees.transpose(featureData.mfdTrees);
}
}
@@ -399,7 +399,7 @@ class GridLoader {
double mMax) {
Mfd.Type type = grProps.type();
if (type == Type.GR_MMAX_GR) {
- builder.transform(xy -> zeroAboveM(xy, 6.5));
+ builder.transform(xy -> zeroAboveM(xy, Math.min(mMax, 6.5)));
} else if (type == Type.GR_MMAX_SINGLE) {
/*
* Keeps the magnitude bin just above mMax if mMax > m_gr + Δm/2. Rounding
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/HazardModel.java b/src/main/java/gov/usgs/earthquake/nshmp/model/HazardModel.java
index e85733b6a7e3cb4f03356ad80d1264fcfd1884a3..6cf6ee8c073ec0e12eb95adf0f3a04ef77e426bd 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/HazardModel.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/HazardModel.java
@@ -17,7 +17,6 @@ import com.google.common.collect.SetMultimap;
import com.google.common.collect.Sets;
import com.google.common.collect.Streams;
-import gov.usgs.earthquake.nshmp.Faults;
import gov.usgs.earthquake.nshmp.Maths;
import gov.usgs.earthquake.nshmp.calc.CalcConfig;
import gov.usgs.earthquake.nshmp.geo.LocationList;
@@ -365,10 +364,7 @@ public final class HazardModel implements Iterable<SourceSet<? extends Source>>
.id(feature.id)
.trace(feature.trace)
.dip(feature.dip)
- .width(Faults.width(
- feature.dip,
- feature.upperDepth,
- feature.lowerDepth))
+ .width(feature.width)
.depth(feature.upperDepth)
.rake(feature.rake)
.mfdTree(frs.mfdTree)
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/MfdTrees.java b/src/main/java/gov/usgs/earthquake/nshmp/model/MfdTrees.java
index 2d5995293343365dd827865bcef50c9bb2c00f46..c0c2a76c47dad31334e6385938a7e520a74a2568 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/MfdTrees.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/MfdTrees.java
@@ -28,12 +28,16 @@ import gov.usgs.earthquake.nshmp.tree.Branch;
import gov.usgs.earthquake.nshmp.tree.LogicTree;
/**
- * Factory class for manipulating and querying MFD logic trees.
+ * Factory class for manipulating and querying logic trees. Most, but not all,
+ * methods pertain to MFD logic trees.
+ *
*
* @author U.S. Geological Survey
*/
class MfdTrees {
+ // TODO rename to Trees
+
/* All Gutenberg-Richter MFD types. */
private static final Set<Type> GR_TYPES = EnumSet.of(
GR,
@@ -99,7 +103,7 @@ class MfdTrees {
private static Mfd mMaxGrMfd(Mfd.Properties props, double mMax) {
Mfd.Properties newProps = (props.type() == Type.GR_TAPER)
? new TaperedGr(props.getAsGrTaper(), mMax)
- : new GutenbergRichter(props.getAsGr(), mMax);
+ : GutenbergRichter.copyOfWithMax(props.getAsGr(), mMax);
return newProps.toBuilder()
.transform(p -> p.set((p.x() > props.getAsGr().mMax()) ? 0.0 : p.y()))
.build();
@@ -248,7 +252,7 @@ class MfdTrees {
* Transpose a list of logic trees to a logic tree of immutable lists.
* Supplied logic trees are assumed to have same branch names and IDs.
*/
- static <T> LogicTree<List<T>> transposeTree(List<LogicTree<T>> treeList) {
+ static <T> LogicTree<List<T>> transpose(List<LogicTree<T>> treeList) {
LogicTree<T> model = treeList.get(0);
/* Init branch value lists. */
@@ -273,4 +277,54 @@ class MfdTrees {
}
return treeOfLists.build();
}
+
+ /* Join the supplied IDs on ' : ' */
+ static String combineBranchIds(String... ids) {
+ return String.join(" : ", ids);
+ }
+
+ /*
+ * Combine IDs and weights of a tree and branch, preserving tree values.
+ * Branches in resultant list do not likely sum to one.
+ */
+ static <T, R> List<Branch<T>> toBranchList(
+ LogicTree<T> tree,
+ Branch<R> branch) {
+
+ return tree.stream()
+ .map(treeBranch -> new IdWeightBranch<>(
+ combineBranchIds(branch.id(), treeBranch.id()),
+ treeBranch.value(),
+ branch.weight() * treeBranch.weight()))
+ .collect(toUnmodifiableList());
+ }
+
+ private static class IdWeightBranch<T> implements Branch<T> {
+
+ final String id;
+ final double weight;
+ final T value;
+
+ @Override
+ public String id() {
+ return id;
+ }
+
+ @Override
+ public double weight() {
+ return weight;
+ }
+
+ @Override
+ public T value() {
+ return value;
+ }
+
+ IdWeightBranch(String id, T value, double weight) {
+ this.id = id;
+ this.weight = weight;
+ this.value = value;
+ }
+ }
+
}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/ModelLoader.java b/src/main/java/gov/usgs/earthquake/nshmp/model/ModelLoader.java
index 7096ceeef5248f6b4a26b165f7363be5e861c291..316a4986a6bc5454d9b4d21209ad17f76e95b56b 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/ModelLoader.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/ModelLoader.java
@@ -37,6 +37,7 @@ import static gov.usgs.earthquake.nshmp.model.SourceType.ZONE;
import static java.nio.file.Files.newDirectoryStream;
import java.io.IOException;
+import java.math.RoundingMode;
import java.nio.file.DirectoryStream;
import java.nio.file.Files;
import java.nio.file.Path;
@@ -51,6 +52,7 @@ import java.util.logging.Logger;
import com.google.common.collect.Iterables;
import com.google.gson.JsonElement;
+import gov.usgs.earthquake.nshmp.Maths;
import gov.usgs.earthquake.nshmp.fault.FocalMech;
import gov.usgs.earthquake.nshmp.geo.json.Feature;
import gov.usgs.earthquake.nshmp.geo.json.Properties;
@@ -180,12 +182,12 @@ abstract class ModelLoader {
this.root = root;
}
- /* Load source type specific rate tree, if present. */
- abstract void loadRateTree(Path dir, ModelData data);
-
/* Load source type specific config, if present. */
abstract void loadSourceConfig(Path dir, ModelData data);
+ /* Load source type specific rate tree, if present. */
+ abstract void loadRateTree(Path dir, ModelData data);
+
/* Load single (branch) type specific source. */
abstract SourceTree loadSingleSource(
Path geojson,
@@ -277,12 +279,13 @@ abstract class ModelLoader {
System.out.println(" source: " + root.relativize(geojson));
/* Create new feature with updated properties. */
- double length = feature.trace.length();
- double Mw = FaultRuptureSet.magnitudeWc94(length);
+ double Mw = FaultRuptureSet.magnitudeWc94(feature.length);
+ double Δz = feature.lowerDepth - feature.upperDepth;
+ double width = computeWidth(Δz, feature.dip);
Feature nshmFeature = Feature.copyOf(feature.source)
.properties(Properties.fromFeature(feature.source)
- .put(Key.LENGTH, length)
.put(Key.MAGNITUDE, Mw)
+ .put(Key.WIDTH, width)
.build())
.build();
feature = new SourceFeature.NshmFault(nshmFeature);
@@ -321,21 +324,21 @@ abstract class ModelLoader {
.gmms(data.gmms())
.root(dipSourceTree);
- double length = feature.trace.length();
- double Mw = FaultRuptureSet.magnitudeWc94(length);
-
+ double Mw = FaultRuptureSet.magnitudeWc94(feature.length);
for (int i = 0; i < dipTree.size(); i++) {
Branch<Double> dipBranch = dipTree.get(i);
int branchDip = (int) (feature.dip + dipBranch.value());
int rateDip = (int) ((dipSlipModel == FIXED) ? feature.dip : branchDip);
String branchName = String.format("%s [%s°]", feature.name, branchDip);
+ double Δz = feature.lowerDepth - feature.upperDepth;
+ double width = computeWidth(Δz, branchDip);
Feature dipFeature = Feature.copyOf(feature.source)
.properties(Properties.fromFeature(feature.source)
.put(Key.NAME, branchName)
.put(Key.DIP, branchDip)
.put(Key.RATE_DIP, rateDip)
- .put(Key.LENGTH, length)
.put(Key.MAGNITUDE, Mw)
+ .put(Key.WIDTH, width)
.build())
.build();
@@ -348,6 +351,15 @@ abstract class ModelLoader {
return treeBuilder.build();
}
+ /* This exists for consistency with legacy codes */
+ @Deprecated
+ private static final double TO_RAD = 3.14159 / 180.0;
+
+ @Deprecated
+ private static double computeWidth(double Δz, double dip) {
+ return Maths.round(Δz / Math.sin(dip * TO_RAD), 5, RoundingMode.FLOOR);
+ }
+
private void loadSingleRuptureSet(
Branch<Path> branch,
SourceTree.Builder treeBuilder,
@@ -607,13 +619,13 @@ abstract class ModelLoader {
}
@Override
- void loadRateTree(Path dir, ModelData data) {
- readGridRateTree(dir).ifPresent(data::gridRateTree);
+ void loadSourceConfig(Path dir, ModelData data) {
+ readGridConfig(dir).ifPresent(data::sourceConfig);
}
@Override
- void loadSourceConfig(Path dir, ModelData data) {
- readGridConfig(dir).ifPresent(data::sourceConfig);
+ void loadRateTree(Path dir, ModelData data) {
+ readGridRateTree(dir).ifPresent(data::gridRateTree);
}
@Override
@@ -714,7 +726,7 @@ abstract class ModelLoader {
/* Attach singleton as leaf. */
if (ruptureSets.size() == 1) {
treeBuilder.addLeaf(rateBranch, ruptureSets.get(0));
- return;
+ continue;
}
/* Otherwise create logic group of rupture set pseudo paths. */
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/SourceFeature.java b/src/main/java/gov/usgs/earthquake/nshmp/model/SourceFeature.java
index b2fd609fc7a2b032f0b15d3eefb62b18d4a5df39..51d35515fd138d442edcdc6d4f2a0eddb2a91187 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/SourceFeature.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/SourceFeature.java
@@ -1,6 +1,7 @@
package gov.usgs.earthquake.nshmp.model;
import static com.google.common.base.Preconditions.checkArgument;
+import static com.google.common.base.Preconditions.checkNotNull;
import static gov.usgs.earthquake.nshmp.Earthquakes.checkCrustalDepth;
import static gov.usgs.earthquake.nshmp.Faults.checkAseismicityFactor;
import static gov.usgs.earthquake.nshmp.Faults.checkDip;
@@ -8,6 +9,7 @@ import static gov.usgs.earthquake.nshmp.Faults.checkDipDirection;
import static gov.usgs.earthquake.nshmp.Faults.checkRake;
import static gov.usgs.earthquake.nshmp.Faults.checkSlipRate;
import static gov.usgs.earthquake.nshmp.Text.checkName;
+import static java.util.stream.Collectors.toMap;
import java.awt.geom.Area;
import java.lang.reflect.Type;
@@ -16,13 +18,16 @@ import java.nio.file.Path;
import java.util.ArrayList;
import java.util.List;
import java.util.Map;
+import java.util.Map.Entry;
import java.util.NoSuchElementException;
import java.util.Optional;
import java.util.OptionalDouble;
import com.google.common.reflect.TypeToken;
+import com.google.gson.JsonElement;
import gov.usgs.earthquake.nshmp.Earthquakes;
+import gov.usgs.earthquake.nshmp.Faults;
import gov.usgs.earthquake.nshmp.geo.Location;
import gov.usgs.earthquake.nshmp.geo.LocationList;
import gov.usgs.earthquake.nshmp.geo.Locations;
@@ -199,6 +204,8 @@ public abstract class SourceFeature {
static final String RATE = "rate";
static final String RATE_MAP = "rate-map";
static final String RATE_TYPE = "rate-type";
+ static final String GR_A_VALUE = "gr-a-value";
+ static final String SINGLE_RATE = "single-rate";
static final String STRIKE = "strike";
static final String MFD_TREE = "mfd-tree";
static final String ASEISMICITY = "aseismicity";
@@ -273,11 +280,17 @@ public abstract class SourceFeature {
/** The fault section rake. */
final double rake;
+ /** The fault section length. */
+ final double length;
+
+ /** The fault section width. */
+ final double width;
+
/** The fault section {@code RateType}. */
final RateType rateType;
/** The mapping of rates from different deformation models. */
- final Map<String, Double> rateMap;
+ final Map<String, Rates> rateMap;
/*
* Fields for use in FaultRuptureSet.Builder. length, width, and magnitude
@@ -291,14 +304,14 @@ public abstract class SourceFeature {
* parameters as used in hazard calculations.
*/
final OptionalDouble rateDip;
- final OptionalDouble length;
- final OptionalDouble width;
final OptionalDouble magnitude;
NshmFault(Feature section) {
super(section);
trace = section.asLineString();
Properties props = section.properties();
+
+ /* Required properties. */
dip = props.getDouble(Key.DIP).orElseThrow();
checkDip(dip);
upperDepth = props.getDouble(Key.UPPER_DEPTH).orElseThrow();
@@ -310,31 +323,86 @@ public abstract class SourceFeature {
rateType = RateType.valueOf(props.getString(Key.RATE_TYPE).orElseThrow());
rateMap = nshmRateMap(props);
props.getString(Key.STATE).orElseThrow();
+
+ /* Properties computed if absent. */
+ length = props.getDouble(Key.LENGTH).orElseGet(trace::length);
+ width = props.getDouble(Key.WIDTH).orElseGet(this::computeWidth);
+
+ /* Optional properties. */
rateDip = props.getDouble(Key.RATE_DIP);
- length = props.getDouble(Key.LENGTH);
- width = props.getDouble(Key.WIDTH);
magnitude = props.getDouble(Key.MAGNITUDE);
+ }
+ private double computeWidth() {
+ return Faults.projectVertical(dip, lowerDepth - upperDepth);
}
- }
- /*
- * A NSHM feature can have either single rate (maps to 'GEO'), no rate because
- * type is RECURRENCE (maps to 'GEO' with rate of 0), a rate-map with multiple
- * rates whose type can be SLIP or RECURRENCE.
- */
- private static Map<String, Double> nshmRateMap(Properties props) {
- OptionalDouble rate = props.getDouble(Key.RATE);
- RateType type = props.get(Key.RATE_TYPE, RateType.class).orElseThrow();
- if (rate.isPresent()) {
- return Map.of("GEO", rate.orElseThrow());
+ /* Utility container for rate values. */
+ static class Rates {
+
+ /* Slip or recurrence rate */
+ final double rate;
+
+ /* Pre-computed Gutenberg-Richter MFD a-value. */
+ final OptionalDouble graValue;
+
+ /* Pre-computed single MFD rate. */
+ final OptionalDouble singleRate;
+
+ Rates(double rate) {
+ this.rate = rate;
+ this.graValue = OptionalDouble.empty();
+ this.singleRate = OptionalDouble.empty();
+ }
+
+ Rates(JsonElement json) {
+ Map<String, Double> rateMap = Deserialize.fromJson(json, RATE_DATA_TOKEN);
+ this.rate = checkNotNull(rateMap.get(Key.RATE));
+ Double graValue = rateMap.get(Key.GR_A_VALUE);
+ this.graValue = (graValue != null)
+ ? OptionalDouble.of(graValue)
+ : OptionalDouble.empty();
+ Double singleRate = rateMap.get(Key.SINGLE_RATE);
+ this.singleRate = (singleRate != null)
+ ? OptionalDouble.of(singleRate)
+ : OptionalDouble.empty();
+ }
+ }
+
+ private static Type RATE_MAP_TOKEN = new TypeToken<Map<String, JsonElement>>() {}.getType();
+ private static Type RATE_DATA_TOKEN = new TypeToken<Map<String, Double>>() {}.getType();
+
+ /*
+ * A NSHM feature can have either single rate (maps to 'GEO'), no rate
+ * because type is RECURRENCE (maps to 'GEO' with rate of 0), a rate-map
+ * with multiple rates whose type can be SLIP or RECURRENCE.
+ */
+ private static Map<String, Rates> nshmRateMap(Properties props) {
+ OptionalDouble rate = props.getDouble(Key.RATE);
+ RateType type = props.get(Key.RATE_TYPE, RateType.class).orElseThrow();
+ if (rate.isPresent()) {
+ return Map.of("GEO", new Rates(rate.orElseThrow()));
+ }
+
+ Optional<Map<String, Rates>> rateMap = readRateMap(props);
+ if (rateMap.isEmpty() && type == RateType.RECURRENCE) {
+ return Map.of("GEO", new Rates(0.0));
+ }
+ return rateMap.orElseThrow();
}
- Type mapType = new TypeToken<Map<String, Double>>() {}.getType();
- Optional<Map<String, Double>> rateMap = props.get(Key.RATE_MAP, mapType);
- if (rateMap.isEmpty() && type == RateType.RECURRENCE) {
- return Map.of("GEO", 0.0);
+
+ private static Optional<Map<String, Rates>> readRateMap(Properties props) {
+ Optional<Map<String, JsonElement>> rateElements = props.get(Key.RATE_MAP, RATE_MAP_TOKEN);
+ if (rateElements.isEmpty()) {
+ return Optional.empty();
+ }
+ Map<String, Rates> rateMap = rateElements.orElseThrow().entrySet().stream()
+ .filter(e -> !e.getValue().isJsonNull())
+ .collect(toMap(
+ Entry::getKey,
+ e -> new Rates(e.getValue())));
+ return Optional.of(rateMap);
}
- return rateMap.orElseThrow();
}
/** A sub-section of a crustal fault in a NSHM fault-system. */
diff --git a/src/test/java/gov/usgs/earthquake/nshmp/mfd/MfdTests.java b/src/test/java/gov/usgs/earthquake/nshmp/mfd/MfdTests.java
index 6270622bdf0bdace3ad89574687e77dea941a3fa..7d9ab76955b0119734a67e0dfaa2dac8d7b8b70e 100644
--- a/src/test/java/gov/usgs/earthquake/nshmp/mfd/MfdTests.java
+++ b/src/test/java/gov/usgs/earthquake/nshmp/mfd/MfdTests.java
@@ -186,7 +186,7 @@ class MfdTests {
assertEquals(5.0, grProps.mMin());
assertEquals(5.5, grProps.mMax());
- grProps = new GutenbergRichter(grProps, 7.5);
+ grProps = GutenbergRichter.copyOfWithMax(grProps, 7.5);
assertEquals(1.0, grProps.a());
assertEquals(1.0, grProps.b());
assertEquals(0.1, grProps.Δm());