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 a09ed62251d82e6788f04c05072f15de33fe6145..b0340803b76435a7233f48f5b85d37160370462e 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/mfd/Mfd.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/mfd/Mfd.java
@@ -4,6 +4,8 @@ import static com.google.common.base.Preconditions.checkNotNull;
import static com.google.common.base.Preconditions.checkState;
import static gov.usgs.earthquake.nshmp.Earthquakes.checkMagnitude;
import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR;
+import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR_MMAX_GR;
+import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR_MMAX_SINGLE;
import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR_TAPER;
import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.SINGLE;
import static gov.usgs.earthquake.nshmp.mfd.Mfds.checkGaussianMagnitudeCount;
@@ -324,9 +326,6 @@ public final class Mfd {
* double, double)
*/
GR_TAPER;
-
- /** A set containing all Gutenberg-Richter MFD types. */
- public static final Set<Type> GR_TYPES = EnumSet.of(GR, GR_MMAX_GR, GR_MMAX_SINGLE, GR_TAPER);
}
/**
@@ -644,6 +643,18 @@ public final class Mfd {
this(GR, a, b, Δm, mMin, mMax);
}
+ /**
+ * Create a new properties object for a Gutenberg-Richter MFD with a copy
+ * of the supplied properties and an updated maximum magnitude.
+ *
+ * @param gr properties to copy
+ * @param mMax updated 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);
+ }
+
/** The Gutenberg-Richter {@code a}-value. */
public double a() {
return a;
@@ -669,9 +680,14 @@ public final class Mfd {
return mMax;
};
+ private static final Set<Type> GR_TYPES = EnumSet.of(
+ GR,
+ GR_MMAX_GR,
+ GR_MMAX_SINGLE);
+
@Override
public Builder toBuilder() {
- checkState(type().equals(GR));
+ checkState(GR_TYPES.contains(type()));
checkGrBvalue(this.b);
/* Pre-check magnitudes; final arrays are checked in builder */
double mMin = checkMagnitude(this.mMin);
@@ -680,7 +696,7 @@ public final class Mfd {
.scale(4)
.build();
Builder builder = new Builder(this, magnitudes);
- builder.mfd.forEach(p -> p.set(Mfds.gutenbergRichterRate(this.a, this.b, p.x())));
+ builder.transform(p -> p.set(Mfds.gutenbergRichterRate(a, b, p.x())));
return builder;
}
@@ -711,8 +727,7 @@ public final class Mfd {
}
/**
- * Create a new properties object for a tapered Gutenberg-Richter MFD with
- * an initial {@code a}-value of one.
+ * Create a new properties object for a tapered Gutenberg-Richter MFD.
*
* @param a value of the distribution
* @param b value of the distribution
@@ -725,6 +740,18 @@ public final class Mfd {
this.mc = mc;
}
+ /**
+ * Create a new properties object for a tapered Gutenberg-Richter MFD with
+ * a copy of the supplied properties and an updated maximum magnitude.
+ *
+ * @param gr properties to copy
+ * @param mMax updated maximum magnitude
+ */
+ public TaperedGr(TaperedGr gr, double mMax) {
+ /* Supports common magnitude arrays for variable mMax MFDs */
+ this(gr.a(), gr.b(), gr.Δm(), gr.mMin(), mMax, gr.mc);
+ }
+
/** The corner magnitude. */
public final double mc() {
return mc;
@@ -743,11 +770,28 @@ public final class Mfd {
.build();
Builder builder = new Builder(this, magnitudes);
TaperFunction grTaper = new TaperFunction(Δm(), b(), mc);
- builder.mfd.forEach(p -> p.set(
- Mfds.gutenbergRichterRate(a(), b(), p.x()) * grTaper.scale(p.x())));
+ builder.transform(p -> p.set(taperedRate(p, grTaper)));
+
+ /*
+ * The tapering function actually does not recover the un-tapered rates;
+ * it's close, but not exact. The difference, however, is uniform across
+ * all m-bins. In order to get the same result by creating a tapered GR
+ * 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.
+ */
+ double grRate = Mfds.gutenbergRichterRate(a(), b(), builder.mfd.x(0));
+ double tgrRate = builder.mfd.y(0);
+ double aCorrection = grRate / tgrRate;
+ builder.scale(aCorrection);
+
return builder;
}
+ private double taperedRate(XyPoint p, TaperFunction fn) {
+ return Mfds.gutenbergRichterRate(a(), b(), p.x()) * fn.scale(p.x());
+ }
+
@Override
public String toString() {
return Mfds.propsToString(
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/GmmSet.java b/src/main/java/gov/usgs/earthquake/nshmp/model/GmmSet.java
index b2081955556ef0e907fec593fef44e4f020534db..d0128579ba98ef0208638c9ec5b67ea4dcb1562e 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/GmmSet.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/GmmSet.java
@@ -8,14 +8,11 @@ import static gov.usgs.earthquake.nshmp.data.DoubleData.checkWeights;
import java.util.Map;
import java.util.Objects;
-import java.util.Optional;
import java.util.Set;
-import java.util.stream.StreamSupport;
import com.google.common.collect.Maps;
import com.google.common.collect.Range;
-import gov.usgs.earthquake.nshmp.data.DoubleData;
import gov.usgs.earthquake.nshmp.gmm.Gmm;
import gov.usgs.earthquake.nshmp.gmm.GroundMotionModel;
import gov.usgs.earthquake.nshmp.gmm.UncertaintyModel;
@@ -51,47 +48,44 @@ public final class GmmSet {
private final int hashCode;
private final UncertType uncertainty;
- private final double epiValue;
- private final double[][] epiValues;
- private final double[] epiWeights; // TODO DataArray (vs. Table, Volume)
-
- private final Optional<UncertaintyModel> epiModel;
-
- GmmSet(
- Map<Gmm, Double> weightMapLo,
- double maxDistLo,
- Map<Gmm, Double> weightMapHi,
- double maxDistHi,
- double[] epiValues,
- double[] epiWeights,
- Optional<UncertaintyModel> epiModel) {
-
- this.weightMapLo = weightMapLo;
- this.maxDistLo = maxDistLo;
- this.weightMapHi = weightMapHi;
- this.maxDistHi = (weightMapHi != null) ? maxDistHi : maxDistLo;
- this.singular = weightMapHi == null;
+ private final double[] epiWeights;
+
+ private final UncertaintyModel epiModel;
+
+ GmmSet(Builder builder) {
+
+ this.weightMapLo = builder.gmmWtMapLo;
+ this.maxDistLo = builder.maxDistanceLo;
+ this.weightMapHi = builder.gmmWtMapHi;
+ this.maxDistHi = (this.weightMapHi != null)
+ ? builder.maxDistanceHi
+ : maxDistLo;
+ this.singular = this.weightMapHi == null;
this.hashCode = Objects.hash(
this.weightMapLo, this.maxDistLo,
this.weightMapHi, this.maxDistHi);
- uncertainty = (epiValues == null) ? UncertType.NONE : (epiValues.length == 1)
- ? UncertType.SINGLE : UncertType.MULTI;
-
- this.epiWeights = epiWeights;
- if (uncertainty == UncertType.NONE) {
- this.epiValue = Double.NaN;
- this.epiValues = null;
- } else if (uncertainty == UncertType.SINGLE) {
- this.epiValue = epiValues[0];
- this.epiValues = null;
- } else {
- this.epiValue = Double.NaN;
- this.epiValues = initEpiValues(epiValues);
- }
+ this.uncertainty = (builder.epiModel == null)
+ ? UncertType.NONE
+ : UncertType.MULTI;
+
+ this.epiWeights = builder.uncWeights;
+
+ // if (uncertainty == UncertType.NONE) {
+ // this.epiValue = Double.NaN;
+ // this.epiValues = null;
+ // } else {
+ // this.epiValue = Double.NaN;
+ // this.epiValues = initEpiValues(builder.uncValues);
+ // }
- this.epiModel = epiModel;
+ this.epiModel = builder.epiModel;
+
+ // System.out.println(epiModel);
+ // System.out.println(epiValues);
+ // System.out.println(epiWeights);
+ // System.out.println(uncertainty);
}
/**
@@ -142,13 +136,13 @@ public final class GmmSet {
this.maxDistHi == that.maxDistHi;
}
- private static double[][] initEpiValues(double[] v) {
- return new double[][] {
- { v[0], v[1], v[2] },
- { v[3], v[4], v[5] },
- { v[6], v[7], v[8] }
- };
- }
+ // private static double[][] initEpiValues(double[] v) {
+ // return new double[][] {
+ // { v[0], v[1], v[2] },
+ // { v[3], v[4], v[5] },
+ // { v[6], v[7], v[8] }
+ // };
+ // }
public boolean epiUncertainty() {
return uncertainty != UncertType.NONE;
@@ -163,7 +157,7 @@ public final class GmmSet {
public double epiValue(double m, double r) {
switch (uncertainty) {
case MULTI:
- return epiModel.orElseThrow().value(m, r);
+ return epiModel.value(m, r);
default:
return 0.0;
}
@@ -189,9 +183,9 @@ public final class GmmSet {
return epiWeights;
}
+ // TODO put NONE in UncertaintyModel
private static enum UncertType {
NONE,
- SINGLE,
MULTI;
}
@@ -209,63 +203,69 @@ public final class GmmSet {
private double maxDistanceHi;
// optional
- private double[] uncValues;
+ // private double[] uncValues;
+ private UncertaintyModel epiModel;
private double[] uncWeights;
- private Optional<UncertaintyModel> epiModel;
-
// leave maxDistanceHi as primitive unless validation required
// at some later date; GmmSet throws NPE if Double useds
// refactor uncertainty only set with primary gmms
Builder primaryModels(LogicTree<Gmm> tree, GmmConfig config) {
- // checkArgument(checkNotNull(tree, "Map is null").size() > 0, "Map is
- // empty");
+ checkNotNull(tree);
+ checkNotNull(config);
// setBuilder.uncertainty(uncValues, uncWeights); TODO from gmm-config
- gmmWtMapLo = StreamSupport.stream(tree.spliterator(), false)
+ gmmWtMapLo = tree.stream()
.collect(Maps.toImmutableEnumMap(
Branch::value,
Branch::weight));
- primaryMaxDistance(config.maxDistance);
-
- epiModel = config.epistemicModel;
+ maxDistanceLo = config.maxDistance;
+ // TODO use tree instead of extracting weights
+ if (config.epistemicModel.isPresent()) {
+ epiModel = config.epistemicModel.orElseThrow();
+ uncWeights = config.epistemicTree.orElseThrow().stream()
+ .mapToDouble(Branch::weight)
+ .toArray();
+ }
return this;
}
- Builder primaryModels(Map<Gmm, Double> gmmWtMap, GmmConfig config) {
- checkArgument(checkNotNull(gmmWtMap, "Map is null").size() > 0, "Map is empty");
- checkWeights(gmmWtMap.values());
- gmmWtMapLo = Maps.immutableEnumMap(gmmWtMap);
-
- primaryMaxDistance(config.maxDistance);
-
- // TODO just pass config to GmmSet
- // if (config.epistemicTree.isPresent()) {
- // double[] values = new double[3];
- // double[] weights = new double[3];
- // List<Branch<Double>> epiBranches =
- // config.epistemicTree.get().branches();
- // checkState(epiBranches.size() == 3);
- // for (int i = 0; i < 3; i++) {
- // Branch<Double> branch = epiBranches.get(0);
- // values[i] = branch.value();
- // weights[i] = branch.weight();
- // }
- // uncertainty(values, weights);
- // }
-
- return this;
- }
+ // Builder primaryModels(Map<Gmm, Double> gmmWtMap, GmmConfig config) {
+ // checkArgument(checkNotNull(gmmWtMap, "Map is null").size() > 0, "Map is
+ // empty");
+ // checkWeights(gmmWtMap.values());
+ // gmmWtMapLo = Maps.immutableEnumMap(gmmWtMap);
+ //
+ // primaryMaxDistance(config.maxDistance);
+ //
+ // // TODO just pass config to GmmSet
+ // // if (config.epistemicTree.isPresent()) {
+ // // double[] values = new double[3];
+ // // double[] weights = new double[3];
+ // // List<Branch<Double>> epiBranches =
+ // // config.epistemicTree.get().branches();
+ // // checkState(epiBranches.size() == 3);
+ // // for (int i = 0; i < 3; i++) {
+ // // Branch<Double> branch = epiBranches.get(0);
+ // // values[i] = branch.value();
+ // // weights[i] = branch.weight();
+ // // }
+ // // uncertainty(values, weights);
+ // // }
+ //
+ // return this;
+ // }
- @Deprecated
- Builder primaryMaxDistance(double maxDistance) {
- maxDistanceLo = checkInRange(MAX_DIST_RANGE, "Max distance", maxDistance);
- return this;
- }
+ // @Deprecated
+ // Builder primaryMaxDistance(double maxDistance) {
+ // maxDistanceLo = checkInRange(MAX_DIST_RANGE, "Max distance",
+ // maxDistance);
+ // return this;
+ // }
Builder secondaryModels(Map<Gmm, Double> gmmWtMap, GmmConfig config) {
checkArgument(checkNotNull(gmmWtMap, "Map is null").size() > 0, "Map is empty");
@@ -283,17 +283,20 @@ public final class GmmSet {
return this;
}
- Builder uncertainty(double[] values, double[] weights) {
- checkNotNull(values, "Values is null");
- checkArgument(values.length == 9 || values.length == 1,
- "Values must contain 1 or 9 values");
- checkArgument(checkNotNull(weights, "Weights are null").length == 3,
- "Weights must contain 3 values");
- checkArgument(DoubleData.sum(weights) == 1.0, "%s uncertainty weights must sum to 1");
- uncValues = values;
- uncWeights = weights;
- return this;
- }
+ // TODO need to checkArgument that values is 9 long
+
+ // Builder uncertainty(double[] values, double[] weights) {
+ // checkNotNull(values, "Values is null");
+ // checkArgument(values.length == 9 || values.length == 1,
+ // "Values must contain 1 or 9 values");
+ // checkArgument(checkNotNull(weights, "Weights are null").length == 3,
+ // "Weights must contain 3 values");
+ // checkArgument(DoubleData.sum(weights) == 1.0, "%s uncertainty weights
+ // must sum to 1");
+ // uncValues = values;
+ // uncWeights = weights;
+ // return this;
+ // }
void validateState(String id) {
@@ -314,11 +317,11 @@ public final class GmmSet {
maxDistanceLo);
}
- if (uncValues != null) {
+ if (epiModel != null) {
checkNotNull(uncWeights, "%s uncertainty weights not set", id);
}
if (uncWeights != null) {
- checkNotNull(uncValues, "%s uncertainty values not set", id);
+ checkNotNull(epiModel, "%s uncertainty values not set", id);
}
built = true;
@@ -328,11 +331,7 @@ public final class GmmSet {
validateState(ID);
try {
- GmmSet gmmSet = new GmmSet(
- gmmWtMapLo, maxDistanceLo,
- gmmWtMapHi, maxDistanceHi,
- uncValues, uncWeights,
- epiModel);
+ GmmSet gmmSet = new GmmSet(this);
return gmmSet;
} catch (Exception e) {
e.printStackTrace();
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 37c48725cba9294a2e8f5f7cc1d026dd549269ca..0751493ae567f1ffa5484db8e48e20962e8c6884 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/GridLoader.java
@@ -16,6 +16,7 @@ import java.util.Optional;
import java.util.OptionalDouble;
import java.util.stream.Stream;
+import gov.usgs.earthquake.nshmp.Maths;
import gov.usgs.earthquake.nshmp.data.DelimitedData;
import gov.usgs.earthquake.nshmp.data.DelimitedData.Record;
import gov.usgs.earthquake.nshmp.data.XyPoint;
@@ -24,6 +25,7 @@ import gov.usgs.earthquake.nshmp.geo.Location;
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.Type;
import gov.usgs.earthquake.nshmp.mfd.Mfds;
import gov.usgs.earthquake.nshmp.tree.Branch;
@@ -292,6 +294,8 @@ class GridLoader {
.collect(toUnmodifiableList());
}
+ // private boolean temp = true;
+
private void processRecord(Record record, boolean mMaxPresent) {
Location loc = Location.create(
@@ -307,9 +311,19 @@ class GridLoader {
FeatureData featureData = dataMap.get(gridId);
featureData.locations.add(loc);
- LogicTree<Mfd> modelMfd = modelMfdTrees.get(gridId);
- LogicTree<Mfd> nodeMfd = createNodeMfdTree(modelMfd, a, mMax);
- featureData.mfdTrees.add(nodeMfd);
+ LogicTree<Mfd> modelMfds = modelMfdTrees.get(gridId);
+ LogicTree<Mfd> nodeMfds = createNodeMfdTree(modelMfds, a, mMax);
+
+ // if (temp) {
+ // System.out.println(nodeMfds);
+ // System.out.println(a + " " + mMax);
+ // for (Branch<Mfd> branch : nodeMfds) {
+ // System.out.println(branch.value().properties().getAsGr().mMin());
+ // System.out.println(branch.value());
+ // }
+ // temp = false;
+ // }
+ featureData.mfdTrees.add(nodeMfds);
}
private LogicTree<Mfd> createNodeMfdTree(
@@ -321,10 +335,11 @@ class GridLoader {
for (Branch<Mfd> branch : modelTree) {
Mfd.Properties.GutenbergRichter grProps = branch.value().properties().getAsGr();
- double rate = Mfds.gutenbergRichterRate(a, grProps.b(), grProps.mMin());
+ double mMin = grProps.mMin() + grProps.Δm() * 0.5;
+ double rate = Mfds.gutenbergRichterRate(a, grProps.b(), mMin);
Mfd.Builder nodeMfd = Mfd.Builder.from(branch.value())
.scaleToIncrementalRate(rate);
- mMax.ifPresent(m -> truncate(grProps.type(), nodeMfd, m));
+ mMax.ifPresent(m -> truncate(grProps, nodeMfd, m));
nodeMfdTree.addBranch(
branch.id(),
nodeMfd.build(),
@@ -345,8 +360,7 @@ class GridLoader {
LogicTree<Mfd.Properties> propsTree = Deserialize.mfdTree(
grid.source.properties(),
data);
- MfdTrees.checkGrMinAndDelta(propsTree);
- LogicTree<Mfd> mfdTree = MfdTrees.propsTreeToMmaxMfdTree(propsTree);
+ LogicTree<Mfd> mfdTree = MfdTrees.grPropsTreeToMmaxMfdTree(propsTree);
modelMfds.put(entry.getKey(), mfdTree);
}
return Map.copyOf(modelMfds);
@@ -379,10 +393,26 @@ class GridLoader {
}
/* Truncate special GR cases; e.g. WUS double counting. */
- private static void truncate(Mfd.Type type, Mfd.Builder builder, double mMax) {
- // TODO 6.5 should be obtained from config
- double m = (type == Type.GR_MMAX_GR) ? 6.5 : mMax;
- builder.transform(xy -> zeroAboveM(xy, m));
+ private static void truncate(
+ GutenbergRichter grProps,
+ Mfd.Builder builder,
+ double mMax) {
+ Mfd.Type type = grProps.type();
+ if (type == Type.GR_MMAX_GR) {
+ builder.transform(xy -> zeroAboveM(xy, 6.5));
+ } else if (type == Type.GR_MMAX_SINGLE) {
+ /*
+ * Keeps the magnitude bin just above mMax if mMax > m_gr + Δm/2. Rounding
+ * to 5 decimal places ensures that mMax/Δm doesn't equal *.00...001,
+ * which ceil() would take up to *.1
+ *
+ * We really should just be cutting out the uppermost bin with overlap but
+ * this is consistent with what was being done in legacy fortran codes.
+ */
+ double Δm = grProps.Δm();
+ double m = Math.ceil(Maths.round(mMax / Δm, 5)) * Δm;
+ builder.transform(xy -> zeroAboveM(xy, m));
+ }
}
static GridSourceSet createGrid(
@@ -405,7 +435,7 @@ class GridLoader {
builder.gridConfig(config);
props.getDouble(Key.STRIKE).ifPresent(builder::strike);
- System.out.println(mfds.get(0).properties().type());
+ // System.out.println(mfds.get(0).properties().type());
builder.locations(locations, mfds, mfdsTree, focalMechMaps);
return builder.build();
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/GridRuptureSet.java b/src/main/java/gov/usgs/earthquake/nshmp/model/GridRuptureSet.java
index 6fd95b1baf66397329d9f1b8ea990777abd486ef..1a263b5b046ce60b1d8160b794ca04df4f7155e9 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/GridRuptureSet.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/GridRuptureSet.java
@@ -48,6 +48,17 @@ class GridRuptureSet implements RuptureSet {
if (gss == null) {
gss = createSourceSet(weight);
}
+
+ // System.out.println("--- GSS ---");
+ // for (PointSource ptSrc : gss) {
+ // System.out.println(ptSrc);
+ // for (Rupture rup : ptSrc) {
+ // System.out.println(rup);
+ // break;
+ // }
+ // }
+ // System.out.println(gss.mechMaps);
+
return gss;
}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/GridSourceSet.java b/src/main/java/gov/usgs/earthquake/nshmp/model/GridSourceSet.java
index 6056c84553311d4b705a2d1b1c27683e54d82a18..ab326580b2bd1419472557f8b0fa8a483ab1775b 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/GridSourceSet.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/GridSourceSet.java
@@ -58,7 +58,7 @@ public class GridSourceSet extends AbstractSourceSet<PointSource> {
final List<Mfd> mfds;
final LogicTree<List<Mfd>> mfdsTree;
// private final Map<FocalMech, Double> mechMap; // default, used for copyOf
- private final List<Map<FocalMech, Double>> mechMaps; // may be nCopies
+ final List<Map<FocalMech, Double>> mechMaps; // may be nCopies
private final boolean singularMechs;
private final NavigableMap<Double, Map<Double, Double>> magDepthMap;
private final OptionalDouble strike;
@@ -107,17 +107,25 @@ public class GridSourceSet extends AbstractSourceSet<PointSource> {
// this.optimizable = (sourceType() != FIXED_STRIKE) && !Double.isNaN(Δm);
this.optimizable = this.magMaster.isPresent();
- System.out.println(Arrays.toString(magMaster.orElseThrow()));
+ // System.out.println(Arrays.toString(magMaster.orElseThrow()));
double[] depthMags = this.mfds.get(0).data().xValues().toArray();
// System.out.println(mfdsTree);
- System.out.println(Arrays.toString(depthMags));
+ // System.out.println(Arrays.toString(depthMags));
this.depthModel = DepthModel.create(
magDepthMap,
Doubles.asList(depthMags),
maxDepth);
+
+ // System.out.println("singulearMechs: " + singularMechs);
+ // System.out.println("ptSrcType: " + sourceType);
+ // for (int i = 0; i < locations.size(); i++) {
+ // System.out.println(locations.get(i));
+ // System.out.println(mfds.get(i));
+ // }
+ // System.out.println(mfds);
}
@Override
@@ -266,8 +274,10 @@ public class GridSourceSet extends AbstractSourceSet<PointSource> {
// this.spacing = gridConfig.spacing;
this.maxDepth = gridConfig.maxDepth;
+
// TODO makes bad assumption of GRID
- validateDepth(this.maxDepth, SourceType.GRID);
+ // TODO reenable but need slab handling
+ // validateDepth(this.maxDepth, SourceType.GRID);
this.rupScaling = gridConfig.ruptureScaling;
this.sourceType = gridConfig.pointSourceType;
@@ -281,7 +291,8 @@ public class GridSourceSet extends AbstractSourceSet<PointSource> {
validateMagCutoffs(this.magDepthMap);
// TODO makes bad assumption of GRID
- validateDepthMap(this.magDepthMap, SourceType.GRID);
+ // TODO reenable but need slab handling
+ // validateDepthMap(this.magDepthMap, SourceType.GRID);
return this;
}
@@ -718,6 +729,8 @@ public class GridSourceSet extends AbstractSourceSet<PointSource> {
Mfd.Properties props = modelMfd.properties(); // probably INCR
double[] mags = modelMfd.data().xValues().toArray();
+ // TODO can we get this from the mfd.properties?
+ // the original min max won't have offsets
double Δm = mags[1] - mags[0];
double ΔmBy2 = Δm / 2.0;
double mMin = mags[0] - ΔmBy2;
@@ -749,6 +762,8 @@ public class GridSourceSet extends AbstractSourceSet<PointSource> {
// System.out.println(parent.name());
// System.out.println(mfdTable);
+ // System.out.println("TableSum: " + mfdTable.collapse().sum());
+ // System.out.println(mfdTable.rows());
List<Double> distances = mfdTable.rows();
maximumSize = distances.size();
@@ -759,7 +774,6 @@ public class GridSourceSet extends AbstractSourceSet<PointSource> {
continue;
}
Location loc = Locations.location(origin, SRC_TO_SITE_AZIMUTH, r);
-
b.add(PointSources.pointSource(
parent.type(),
parent.sourceType,
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 441ce5fb23982c8017b3b1c319d2b3a5f8e2abbc..ad3421e52a8ca8ae6baa963138b9e141042e1d02 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/MfdTrees.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/MfdTrees.java
@@ -1,13 +1,17 @@
package gov.usgs.earthquake.nshmp.model;
import static com.google.common.base.Preconditions.checkArgument;
-import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR_TYPES;
+import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR;
+import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR_MMAX_GR;
+import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR_MMAX_SINGLE;
+import static gov.usgs.earthquake.nshmp.mfd.Mfd.Type.GR_TAPER;
import static java.util.stream.Collectors.toMap;
import static java.util.stream.Collectors.toSet;
import static java.util.stream.Collectors.toUnmodifiableList;
import java.nio.file.Path;
import java.util.ArrayList;
+import java.util.EnumSet;
import java.util.List;
import java.util.Map;
import java.util.Set;
@@ -17,6 +21,8 @@ import gov.usgs.earthquake.nshmp.data.XySequence;
import gov.usgs.earthquake.nshmp.mfd.Mfd;
import gov.usgs.earthquake.nshmp.mfd.Mfd.Properties;
import gov.usgs.earthquake.nshmp.mfd.Mfd.Properties.GutenbergRichter;
+import gov.usgs.earthquake.nshmp.mfd.Mfd.Properties.TaperedGr;
+import gov.usgs.earthquake.nshmp.mfd.Mfd.Type;
import gov.usgs.earthquake.nshmp.mfd.Mfds;
import gov.usgs.earthquake.nshmp.tree.Branch;
import gov.usgs.earthquake.nshmp.tree.LogicTree;
@@ -28,7 +34,15 @@ import gov.usgs.earthquake.nshmp.tree.LogicTree;
*/
class MfdTrees {
+ /* All Gutenberg-Richter MFD types. */
+ private static final Set<Type> GR_TYPES = EnumSet.of(
+ GR,
+ GR_MMAX_GR,
+ GR_MMAX_SINGLE,
+ GR_TAPER);
+
/* Convert a logic tree of mfd properties to builders. */
+ @Deprecated // TODO clean
static LogicTree<Mfd.Builder> mfdPropsToBuilders(LogicTree<Mfd.Properties> tree) {
LogicTree.Builder<Mfd.Builder> mfdTree = LogicTree.builder(tree.name());
for (Branch<Mfd.Properties> branch : tree) {
@@ -42,25 +56,55 @@ class MfdTrees {
* b-values possible) to a tree of model MFDs with identical x-values. This
* supports optimizations in XySequence.combine().
*/
- static LogicTree<Mfd> propsTreeToMmaxMfdTree(LogicTree<Mfd.Properties> tree) {
+ static LogicTree<Mfd> grPropsTreeToMmaxMfdTree(LogicTree<Mfd.Properties> tree) {
+ checkGrMinAndDelta(tree);
double mMax = mfdTreeMaxMagnitude(tree);
LogicTree.Builder<Mfd> mfdTree = LogicTree.builder(tree.name());
- tree.stream()
- .forEach(branch -> mfdTree.addBranch(
- branch.id(),
- mMaxGrMfd(branch.value(), mMax),
- branch.weight()));
+ tree.stream().forEach(branch -> mfdTree.addBranch(
+ branch.id(),
+ mMaxGrMfd(branch.value(), mMax),
+ branch.weight()));
return mfdTree.build();
}
/* Find the maximum mMax, checking that tree is all GR MFDs. */
- static double mfdTreeMaxMagnitude(LogicTree<Mfd.Properties> tree) {
+ private static double mfdTreeMaxMagnitude(LogicTree<Mfd.Properties> tree) {
return tree.stream()
.mapToDouble(b -> b.value().getAsGr().mMax())
.max()
.getAsDouble();
}
+ /*
+ * TODO clean/consolidate these notes
+ *
+ * Note: if we create custom MFDs for each M* branch, we'll have the same
+ * reference x-values for each model, but when added, a new set of x-vlues
+ * will be created each time by XySequence.combine(). We therefore create MFDs
+ * with common x-values so that XySequence.combine() just copies those of the
+ * first M* sequence.
+ *
+ * Create a GR mfd with possibly zero rate bins above gr.mMax. The supplied
+ * mMax sets the upper limit of the actual MFD. Zero rate bins will not scale
+ * in new builder.from(); method is for the combination of GR logic trees of
+ * MFDs with similar binning.
+ *
+ * The rate zeroing done below permits combining of GR_TAPER and GR branches,
+ * where GR_TAPER might go as high as M8. It should not be confused with the
+ * rate zeroing that occurrs later on as individual node MFDs are being
+ * created for GR_MMAX_GR and GR_MMAX_SINGLE to avoid double counting.
+ *
+ * once transforming, is OK to cast GR_TAPER as GR to get original mMax
+ */
+ 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);
+ return newProps.toBuilder()
+ .transform(p -> p.set((p.x() > props.getAsGr().mMax()) ? 0.0 : p.y()))
+ .build();
+ }
+
/* Ensure all MFDs in a tree are the same as specified type. */
static void checkType(LogicTree<Mfd.Properties> tree, Mfd.Type type) {
Mfd.Type refType = tree.get(0).value().type();
@@ -78,7 +122,7 @@ class MfdTrees {
}
/* Ensure all MFDs in a tree are some Gutenberg-Richter type. */
- static void checkGrTypes(LogicTree<Mfd.Properties> tree) {
+ private static void checkGrTypes(LogicTree<Mfd.Properties> tree) {
Set<Mfd.Type> types = tree.stream()
.map(Branch::value)
.map(Mfd.Properties::type)
@@ -89,7 +133,14 @@ class MfdTrees {
}
/* Ensure all MFDs are Gutenberg-Richter and that mMin and Δm match. */
- static void checkGrMinAndDelta(LogicTree<Mfd.Properties> tree) {
+ private static void checkGrMinAndDelta(LogicTree<Mfd.Properties> tree) {
+
+ // TODO consider condensing
+ // Properties::getAsGr will fail as checkGrTypes would
+ // (GR_TAPER. is subclass so getAsGr works)
+ // checkMinAndDelta as predicate could throw exception
+ // and then map to max(mMax)
+
checkGrTypes(tree);
GutenbergRichter refProps = tree.get(0).value().getAsGr();
boolean minAndDeltaEqual = tree.stream()
@@ -108,24 +159,6 @@ class MfdTrees {
return (props1.mMin() == props2.mMin() && props1.Δm() == props2.Δm());
}
- private static double checkMinMagnitude(LogicTree<Mfd.Properties> mfdTree) {
- Set<Double> mins = mfdTree.stream()
- .map(Branch::value)
- .map(p -> p.getAsGr().mMin())
- .collect(toSet());
- checkArgument(mins.size() == 1, "Grid mfd-tree has different mfd mMin: %s", mins);
- return mins.iterator().next();
- }
-
- private static double checkMagnitudeDelta(LogicTree<Mfd.Properties> mfdTree) {
- Set<Double> deltas = mfdTree.stream()
- .map(Branch::value)
- .map(p -> p.getAsGr().Δm())
- .collect(toSet());
- checkArgument(deltas.size() == 1, "Grid mfd-tree has different mfd Δm: %s", deltas);
- return deltas.iterator().next();
- }
-
/* Check if the IDs and weights of two mfd-trees are the same. */
static void checkTreeIdsAndWeights(LogicTree<Mfd> tree1, LogicTree<Mfd> tree2) {
checkArgument(
@@ -139,25 +172,6 @@ class MfdTrees {
return tree.stream().collect(toMap(Branch::id, Branch::weight));
}
- /*
- * Note: if we create custom MFDs for each M* branch, we'll have the same
- * reference x-values for each model, but when added, a new set of x-vlues
- * will be created each time by XySequence.combine(). We therefore create MFDs
- * with common x-values so that XySequence.combine() just copies those of the
- * first M* sequence.
- *
- * Create a GR mfd with possibly zero rate bins above gr.mMax. The supplied
- * mMax sets the upper limit of the actual MFD. Zero rate bins will not scale
- * in new builder.from(); method is for the combination of GR logic trees of
- * MFDs with similar binning.
- */
- private static Mfd mMaxGrMfd(Mfd.Properties props, double mMax) {
- GutenbergRichter gr = props.getAsGr();
- return Mfd.newGutenbergRichterBuilder(gr.b(), gr.Δm(), gr.mMin(), mMax)
- .transform(p -> p.set(p.x() > gr.mMax() ? 0.0 : p.y()))
- .build();
- }
-
/* Get a list of branch IDs from a tree. */
static List<String> treeIds(LogicTree<?> tree) {
return tree.stream()
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 cbd3a01f0e7a22a495e35ca039b4333caea0fd6c..7096ceeef5248f6b4a26b165f7363be5e861c291 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/model/ModelLoader.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/ModelLoader.java
@@ -73,7 +73,7 @@ import gov.usgs.earthquake.nshmp.tree.LogicTree;
abstract class ModelLoader {
public static void main(String[] args) {
- Path testModel = Paths.get("../nshm-conus-2018-tmp");
+ Path testModel = Paths.get("../nshm-conus-2018");
HazardModel model = ModelLoader.load(testModel);
System.out.println(model);
}
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 24b371b0a7c18b78dadf74ae1cb9a0b2906ac9ac..6c501842245f6e20e32fca3c0a8a9032c78bac6d 100644
--- a/src/test/java/gov/usgs/earthquake/nshmp/mfd/MfdTests.java
+++ b/src/test/java/gov/usgs/earthquake/nshmp/mfd/MfdTests.java
@@ -10,7 +10,6 @@ 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;
@@ -30,6 +29,9 @@ import gov.usgs.earthquake.nshmp.mfd.Mfd.Type;
class MfdTests {
+ // TODO add test to see if we can slip in negative a value
+ // via GR properties
+
static final double[] M = { 5.05, 5.15, 5.25, 5.35, 5.45 };
static final double[] R = { 0.1, 0.08, 0.06, 0.04, 0.02 };
@@ -183,6 +185,13 @@ class MfdTests {
assertEquals(5.0, grProps.mMin());
assertEquals(5.5, grProps.mMax());
+ grProps = new GutenbergRichter(grProps, 7.5);
+ assertEquals(1.0, grProps.a());
+ assertEquals(1.0, grProps.b());
+ assertEquals(0.1, grProps.Δm());
+ assertEquals(5.0, grProps.mMin());
+ assertEquals(7.5, grProps.mMax());
+
/* props.toString() */
assertEquals(
grProps.type() + " {a=" + grProps.a() +
@@ -193,56 +202,6 @@ 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 = 1.0;
- double aVal = Math.log10(150.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 = new TaperedGr(aVal, bVal, Δm, mMin, mMax, mc)
- .toBuilder()
- .build();
-
- System.out.println(Arrays.toString(mfd.data().xValues().toArray()));
- System.out.println(Arrays.toString(mfd.data().yValues().toArray()));
-
- incrRate = Mfds.gutenbergRichterRate(aVal, bVal, 5.05);
- mfd = Mfd.newTaperedGutenbergRichterBuilder(bVal, Δm, mMin, mMax, mc)
- .scaleToIncrementalRate(incrRate)
- .build();
-
- System.out.println(Arrays.toString(mfd.data().yValues().toArray()));
-
- }
-
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,
@@ -251,100 +210,82 @@ class MfdTests {
7.05, 7.15, 7.25, 7.35, 7.45,
7.55, 7.65, 7.75, 7.85, 7.95 };
- /* using scaleToInrementalRate */
- private static final double[] TAPERED_GR_R_SCALED = {
+ /* Generated using props.toBuilder().build() */
+ private static final double[] TAPERED_GR_R = {
3.6480433574236396E-4,
3.0345403881799134E-4,
2.5242909958337484E-4,
2.099931150361607E-4,
1.7470192588897718E-4,
1.4535446859752395E-4,
- 1.2095189730212225E-4,
+ 1.2095189730212227E-4,
1.0066358352203966E-4,
8.37988347826301E-5,
6.978336683158799E-5,
5.813972403773392E-5,
- 4.847097401109892E-5,
- 4.044710933063091E-5,
+ 4.8470974011098916E-5,
+ 4.04471093306309E-5,
3.37936743657631E-5,
- 2.8282200517277984E-5,
- 2.372208037314646E-5,
+ 2.8282200517277988E-5,
+ 2.3722080373146463E-5,
1.9953542732997602E-5,
- 1.684141264491656E-5,
+ 1.6841412644916558E-5,
1.4269370950942086E-5,
1.213450863835374E-5,
- 1.034219228924177E-5,
- 8.801777631086678E-6,
+ 1.0342192289241769E-5,
+ 8.80177763108668E-6,
7.4247305061395814E-6,
- 6.128255653142835E-6,
- 4.84871327802538E-6,
- 3.5668582219695296E-6,
+ 6.1282556531428345E-6,
+ 4.848713278025381E-6,
+ 3.56685822196953E-6,
2.3358516405761174E-6,
- 1.2823409762916182E-6,
- 5.437829885686744E-7,
+ 1.2823409762916184E-6,
+ 5.437829885686743E-7,
1.5949675112240326E-7 };
- /* using a value directly from properties */
- private static final double[] TAPERED_GR_R_DIRECT = {
- 3.648734771264746E-4,
- 3.0351155247723867E-4,
- 2.5247694248332143E-4,
- 2.100329150418206E-4,
- 1.7473503715378214E-4,
- 1.4538201763727166E-4,
- 1.2097482132130438E-4,
- 1.006826622960901E-4,
- 8.381471718000331E-5,
- 6.979659287661517E-5,
- 5.815074326255646E-5,
- 4.848016071724236E-5,
- 4.0454775273297684E-5,
- 3.3800079282566184E-5,
- 2.8287560844165223E-5,
- 2.3726576420233308E-5,
- 1.9957324528955875E-5,
- 1.684460459870263E-5,
- 1.4272075425536471E-5,
- 1.213680849238659E-5,
- 1.034415244546678E-5,
- 8.803445832444015E-6,
- 7.426137715686032E-6,
- 6.1294171417451564E-6,
- 4.849632254896958E-6,
- 3.5675342487878286E-6,
- 2.3362943546550987E-6,
- 1.2825840184413838E-6,
- 5.438860517858616E-7,
- 1.5952698054966957E-7 };
-
@Test
void testTaperedGr() {
+ double aVal = Math.log10(4.0);
+ double bVal = 0.8;
+ double Δm = 0.1;
+ double mMin = 5.0;
+ double mMax = 8.0;
+ double mc = 7.5;
+
/* Factory builder; covers props.toBuilder() */
- double incrRate = Mfds.gutenbergRichterRate(Math.log10(4.0), 0.8, 5.05);
- Mfd mfd = Mfd.newTaperedGutenbergRichterBuilder(0.8, 0.1, 5.0, 8.0, 7.5)
+ double incrRate = Mfds.gutenbergRichterRate(aVal, bVal, 5.05);
+ Mfd mfd = Mfd.newTaperedGutenbergRichterBuilder(bVal, Δm, mMin, mMax, mc)
.scaleToIncrementalRate(incrRate)
.build();
XySequence xy = mfd.data();
assertTrue(xy.size() == 30);
assertArrayEquals(TAPERED_GR_M, xy.xValues().toArray());
- assertArrayEquals(TAPERED_GR_R_SCALED, xy.yValues().toArray());
+ /* These end up being very slightly different due to scaling */
+ assertArrayEquals(TAPERED_GR_R, xy.yValues().toArray(), 1e-17);
/* Properties */
- double aValue = Math.log10(4.0);
- Properties props = new TaperedGr(aValue, 0.8, 0.1, 5.0, 8.0, 7.5);
+ Properties props = new TaperedGr(aVal, bVal, Δm, mMin, mMax, mc);
TaperedGr grProps = props.getAsGrTaper(); // cover getAs
- assertEquals(aValue, grProps.a());
- assertEquals(0.8, grProps.b());
- assertEquals(0.1, grProps.Δm());
- assertEquals(5.0, grProps.mMin());
- assertEquals(8.0, grProps.mMax());
- assertEquals(7.5, grProps.mc());
+ assertEquals(aVal, grProps.a());
+ assertEquals(bVal, grProps.b());
+ assertEquals(Δm, grProps.Δm());
+ assertEquals(mMin, grProps.mMin());
+ assertEquals(mMax, grProps.mMax());
+ assertEquals(mc, grProps.mc());
assertTrue(xy.size() == 30);
mfd = props.toBuilder().build();
xy = mfd.data();
assertArrayEquals(TAPERED_GR_M, xy.xValues().toArray());
- assertArrayEquals(TAPERED_GR_R_DIRECT, xy.yValues().toArray());
+ assertArrayEquals(TAPERED_GR_R, xy.yValues().toArray());
+
+ grProps = new TaperedGr(grProps, 8.5);
+ assertEquals(aVal, grProps.a());
+ assertEquals(bVal, grProps.b());
+ assertEquals(Δm, grProps.Δm());
+ assertEquals(mMin, grProps.mMin());
+ assertEquals(8.5, grProps.mMax());
+ assertEquals(mc, grProps.mc());
/* props.toString() */
assertEquals(
@@ -357,40 +298,6 @@ class MfdTests {
grProps.toString());
}
- private static final double[] GR_M_TRANDSORM = TAPERED_GR_M;
-
- private static final double[] GR_R_TRANSFORM = {
- 0.0013368764072006192,
- 0.0010619186765762063,
- 8.435119877855239E-4,
- 6.700253882264454E-4,
- 5.322200838503631E-4,
- 4.2275743968966836E-4,
- 3.3580817078525074E-4,
- 2.667419115058385E-4,
- 2.1188063169341338E-4,
- 1.683027681452945E-4,
- 1.3368764072006192E-4,
- 1.0619186765762062E-4,
- 8.435119877855238E-5,
- 6.700253882264454E-5,
- 5.322200838503631E-5,
- 4.227574396896683E-5,
- 3.358081707852507E-5,
- 2.667419115058385E-5,
- 2.1188063169341338E-5,
- 1.683027681452945E-5,
- 1.3368764072006191E-5,
- 1.0619186765762062E-5,
- 8.435119877855238E-6,
- 6.700253882264454E-6,
- 5.322200838503631E-6,
- 4.227574396896683E-6,
- 3.3580817078525076E-6,
- 2.6674191150583848E-6,
- 2.1188063169341338E-6,
- 1.683027681452945E-6 };
-
@Test
void testBuilderTransforms() {