From 79ef4a9cd6e875b88cfc027226d13c3f12ea0c6b Mon Sep 17 00:00:00 2001
From: Peter Powers <pmpowers@usgs.gov>
Date: Mon, 3 Apr 2017 10:50:55 -0600
Subject: [PATCH] added EqRate calculator
---
src/org/opensha2/calc/EqRate.java | 213 +++++++++++++++++-
.../opensha2/eq/model/SystemSourceSet.java | 61 ++++-
2 files changed, 263 insertions(+), 11 deletions(-)
diff --git a/src/org/opensha2/calc/EqRate.java b/src/org/opensha2/calc/EqRate.java
index 1a38e1444..c6be87234 100644
--- a/src/org/opensha2/calc/EqRate.java
+++ b/src/org/opensha2/calc/EqRate.java
@@ -1,24 +1,225 @@
package org.opensha2.calc;
+import org.opensha2.data.IntervalArray;
import org.opensha2.data.XySequence;
+import org.opensha2.eq.model.ClusterSource;
+import org.opensha2.eq.model.ClusterSourceSet;
+import org.opensha2.eq.model.Distance;
+import org.opensha2.eq.model.HazardModel;
+import org.opensha2.eq.model.Rupture;
+import org.opensha2.eq.model.Source;
+import org.opensha2.eq.model.SourceSet;
import org.opensha2.eq.model.SourceType;
+import org.opensha2.eq.model.SystemSourceSet;
+import org.opensha2.geo.Location;
+import org.opensha2.mfd.Mfds;
+import com.google.common.collect.ImmutableMap;
+import com.google.common.collect.Maps;
+
+import java.util.EnumMap;
import java.util.Map;
+import java.util.Map.Entry;
/**
- * Magnitude-frequency distribution data container. The class may be used for a
- * single source or a region.
+ * General purpose magnitude-frequency distribution data container. CUrrent;y
+ * implemented for annual rate only.
*
* @author Peter Powers
*/
public class EqRate {
+ final XySequence totalMfd;
+ final Map<SourceType, XySequence> typeMfds;
+
+ private EqRate(XySequence totalMfd, Map<SourceType, XySequence> typeMfds) {
+ this.totalMfd = totalMfd;
+ this.typeMfds = typeMfds;
+ }
+
+ /*
+ * Developer notes:
+ *
+ * TODO the receiver MFD needs to be built by querying hazard model for min-max
+ * magnitudes; currently it is built using fixed, known values from the 2014
+ * COUS model.
+ */
+
+ public static EqRate createIncremental(
+ HazardModel model,
+ Location location,
+ double distance) {
+
+ IntervalArray modelMfd = IntervalArray.Builder
+ .withRows(4.7, 9.4, 0.1)
+ .build();
+
+ /* Initialize SourceType mfd builders. */
+ Map<SourceType, IntervalArray.Builder> typeMfdBuilders = new EnumMap<>(SourceType.class);
+ for (SourceType type : model.types()) {
+ typeMfdBuilders.put(type, IntervalArray.Builder.fromModel(modelMfd));
+ }
+
+ /* Populate builders. */
+ for (SourceSet<? extends Source> sourceSet : model) {
+ IntervalArray sourceSetMfd = mfd(sourceSet, location, distance, modelMfd);
+ typeMfdBuilders.get(sourceSet.type()).add(sourceSetMfd);
+ }
+
+ /* Compute total and convert to sequences. */
+ IntervalArray.Builder totalMfd = IntervalArray.Builder.fromModel(modelMfd);
+ ImmutableMap.Builder<SourceType, XySequence> typeMfds = ImmutableMap.builder();
+ for (Entry<SourceType, IntervalArray.Builder> entry : typeMfdBuilders.entrySet()) {
+ IntervalArray typeMfd = entry.getValue().build();
+ totalMfd.add(typeMfd);
+ }
+
+ return new EqRate(
+ totalMfd.build().values(),
+ typeMfds.build());
+ }
+
+ /**
+ * Get the total
+ * @param model
+ * @param location
+ * @param distance
+ * @return
+ */
+ public static EqRate createCumulative(
+ HazardModel model,
+ Location location,
+ double distance) {
+
+ EqRate incremental = createIncremental(model, location, distance);
+ XySequence cumulativeTotal = Mfds.toCumulative(incremental.totalMfd);
+ ImmutableMap.Builder<SourceType, XySequence> cumulativeTypes = ImmutableMap.builder();
+ for (Entry<SourceType, XySequence> entry : incremental.typeMfds.entrySet()) {
+ cumulativeTypes.put(
+ entry.getKey(),
+ Mfds.toCumulative(entry.getValue()));
+ }
+ return new EqRate(
+ cumulativeTotal,
+ cumulativeTypes.build());
+ }
+
+ static EqRate combine(EqRate... eqRates) {
+ // validate xs; or not if only for internal use where we know the receiver
+ // mfds are all sources from the same model
+
+ XySequence totalMfd = XySequence.emptyCopyOf(eqRates[0].totalMfd);
+ EnumMap<SourceType, XySequence> typeMfds = new EnumMap<>(SourceType.class);
+ for (EqRate eqRate : eqRates) {
+ totalMfd.add(eqRate.totalMfd);
+ for (Entry<SourceType, XySequence> entry : eqRate.typeMfds.entrySet()) {
+ SourceType type = entry.getKey();
+ if (!typeMfds.containsKey(type)) {
+ XySequence typeMfd = XySequence.emptyCopyOf(totalMfd);
+ typeMfds.put(type, typeMfd);
+ }
+ typeMfds.get(type).add(entry.getValue());
+ }
+ }
+ return new EqRate(
+ XySequence.immutableCopyOf(totalMfd),
+ Maps.immutableEnumMap(typeMfds));
+ }
+
+ private static IntervalArray mfd(
+ SourceSet<? extends Source> sourceSet,
+ Location location,
+ double distance,
+ IntervalArray modelMfd) {
+
+ switch (sourceSet.type()) {
+ case GRID:
+ return gridMfd(sourceSet, location, distance, modelMfd);
+ case SLAB:
+ return gridMfd(sourceSet, location, distance, modelMfd);
+ case CLUSTER:
+ return clusterMfd((ClusterSourceSet) sourceSet, location, distance, modelMfd);
+ case SYSTEM:
+ return systemMfd((SystemSourceSet) sourceSet, location, distance, modelMfd);
+ default:
+ return faultMfd(sourceSet, location, distance, modelMfd);
+ // TODO AREA?
+ }
+ }
+
+ /*
+ * Short-circuit GridSourceSet by summing the relevant node mfds. Handles both
+ * GRID and SLAB types.
+ */
+ private static IntervalArray gridMfd(
+ SourceSet<? extends Source> sourceSet,
+ Location location,
+ double distance,
+ IntervalArray modelMfd) {
+
+ IntervalArray.Builder sourceSetMfd = IntervalArray.Builder.fromModel(modelMfd);
+ for (Source source : sourceSet.iterableForLocation(location, distance)) {
+ for (XySequence mfd : source.mfds()) {
+ sourceSetMfd.addEach(mfd);
+ }
+ }
+ return sourceSetMfd.build();
+ }
+
+ /*
+ * Special case ClusterSourceSet.
+ */
+ private static IntervalArray clusterMfd(
+ ClusterSourceSet sourceSet,
+ Location location,
+ double distance,
+ IntervalArray modelMfd) {
+
+ IntervalArray.Builder sourceSetMfd = IntervalArray.Builder.fromModel(modelMfd);
+ for (Source source : sourceSet.iterableForLocation(location, distance)) {
+ ClusterSource clusterSource = (ClusterSource) source;
+ sourceSetMfd.add(faultMfd(
+ clusterSource.faults(),
+ location,
+ distance,
+ modelMfd));
+ }
+ return sourceSetMfd.build();
+ }
+
+ /*
+ * Special case and delegate to SystemSourceSet.
+ */
+ static IntervalArray systemMfd(
+ SystemSourceSet sourceSet,
+ Location location,
+ double distance,
+ IntervalArray modelMfd) {
+
+ return SystemSourceSet
+ .toRatesFunction(location, distance, modelMfd)
+ .apply(sourceSet);
+ }
+
/*
- * Depending on operation, this could be the raw combination of multiple MFDs,
- * or resampled to some uniform spacing.
+ * Default approach: distance filter on ruptures.
*/
- XySequence totalMfd;
+ private static IntervalArray faultMfd(
+ SourceSet<? extends Source> sourceSet,
+ Location location,
+ double distance,
+ IntervalArray modelMfd) {
- Map<SourceType, XySequence> sourceMfds;
+ IntervalArray.Builder sourceSetMfd = IntervalArray.Builder.fromModel(modelMfd);
+ for (Source source : sourceSet.iterableForLocation(location, distance)) {
+ for (Rupture rupture : source) {
+ Distance d = rupture.surface().distanceTo(location);
+ if (d.rJB <= distance) {
+ sourceSetMfd.add(rupture.mag(), rupture.rate());
+ }
+ }
+ }
+ return sourceSetMfd.build();
+ }
}
diff --git a/src/org/opensha2/eq/model/SystemSourceSet.java b/src/org/opensha2/eq/model/SystemSourceSet.java
index 0431df6b1..7545f89dc 100644
--- a/src/org/opensha2/eq/model/SystemSourceSet.java
+++ b/src/org/opensha2/eq/model/SystemSourceSet.java
@@ -18,6 +18,7 @@ import org.opensha2.calc.InputList;
import org.opensha2.calc.Site;
import org.opensha2.calc.SystemInputList;
import org.opensha2.data.Indexing;
+import org.opensha2.data.IntervalArray;
import org.opensha2.data.XySequence;
import org.opensha2.eq.fault.Faults;
import org.opensha2.eq.fault.surface.GriddedSurface;
@@ -125,8 +126,7 @@ public final class SystemSourceSet extends AbstractSourceSet<SystemSourceSet.Sys
}
@Override
- public Predicate<SystemSource> distanceFilter(Location loc,
- double distance) {
+ public Predicate<SystemSource> distanceFilter(Location loc, double distance) {
BitSet siteBitset = bitsetForLocation(loc, distance);
return new BitsetFilter(siteBitset);
}
@@ -395,6 +395,59 @@ public final class SystemSourceSet extends AbstractSourceSet<SystemSourceSet.Sys
stats);
}
}
+
+ /*
+ * Handle rate calculations internally as SystemSource is not fully implemented.
+ * If/when it is, this should be removed in favor using iterableForLocation
+ * and getRupture(0).
+ */
+
+ /**
+ * Return an instance of a {@code Function} that converts a
+ * {@code SystemSourceSet} to a ground motion model {@code InputList}.
+ *
+ * @param location with which to initialize instance.
+ * @param distance if interest (relevant source radius)
+ * @param modelMfd MFD to populate
+ */
+ public static Function<SystemSourceSet, IntervalArray> toRatesFunction(
+ Location location,
+ double distance,
+ IntervalArray modelMfd) {
+ return new ToRates(location, distance, modelMfd);
+ }
+
+ private static final class ToRates implements Function<SystemSourceSet, IntervalArray> {
+
+ private final Location location;
+ private final double distance;
+ private final IntervalArray modelMfd;
+
+ ToRates(
+ final Location location,
+ final double distance,
+ final IntervalArray modelMfd) {
+
+ this.location = location;
+ this.distance = distance;
+ this.modelMfd = modelMfd;
+ }
+
+ @Override
+ public IntervalArray apply(final SystemSourceSet sourceSet) {
+ IntervalArray.Builder mfdForLocation = IntervalArray.Builder.fromModel(modelMfd);
+ BitSet bitsetForLocation = sourceSet.bitsetForLocation(location, distance);
+ if (bitsetForLocation.isEmpty()) {
+ return modelMfd;
+ }
+ int[] sourceIndices = Indexing.bitsToIndices(bitsetForLocation);
+ for (int i : sourceIndices) {
+ mfdForLocation.add(sourceSet.mags[i], sourceSet.rates[i]);
+ }
+ return mfdForLocation.build();
+ }
+ }
+
/*
* System source calculation pipeline.
@@ -506,9 +559,7 @@ public final class SystemSourceSet extends AbstractSourceSet<SystemSourceSet.Sys
/* Create inputs. */
Map<Integer, double[]> rMap = rMapBuilder.build();
- Function<SystemSource, HazardInput> inputGenerator = new InputGenerator(
- rMap,
- site);
+ Function<SystemSource, HazardInput> inputGenerator = new InputGenerator(rMap, site);
Predicate<SystemSource> rFilter = new BitsetFilter(siteBitset);
Iterable<SystemSource> sources = Iterables.filter(sourceSet, rFilter);
--
GitLab