From 9a82a4f7b5fce5785e23fa2fee05dc8ba354957b Mon Sep 17 00:00:00 2001
From: Peter Powers <pmpowers@usgs.gov>
Date: Fri, 25 Feb 2022 11:46:39 -0700
Subject: [PATCH] cs edits
---
.../nshmp/calc/ConditionalSpectrum.java | 178 +++++++++++-------
1 file changed, 108 insertions(+), 70 deletions(-)
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/calc/ConditionalSpectrum.java b/src/main/java/gov/usgs/earthquake/nshmp/calc/ConditionalSpectrum.java
index dbe85e40..df8329cf 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/calc/ConditionalSpectrum.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/calc/ConditionalSpectrum.java
@@ -8,16 +8,12 @@ import static java.lang.Math.log;
import static java.lang.Math.max;
import static java.lang.Math.min;
import static java.lang.Math.sqrt;
+import static java.util.stream.Collectors.toMap;
-import java.util.EnumSet;
+import java.util.EnumMap;
+import java.util.Map;
import java.util.Set;
-import java.util.function.ToDoubleBiFunction;
-import com.google.common.collect.ArrayTable;
-import com.google.common.collect.ImmutableTable;
-import com.google.common.collect.Table;
-
-import gov.usgs.earthquake.nshmp.Maths;
import gov.usgs.earthquake.nshmp.gmm.Imt;
/**
@@ -27,7 +23,7 @@ import gov.usgs.earthquake.nshmp.gmm.Imt;
*/
class ConditionalSpectrum {
- interface ImtCorrelation extends ToDoubleBiFunction<Imt, Imt> {}
+ // interface ImtCorrelation extends ToDoubleBiFunction<Imt, Imt> {}
/*
* Method 4: "Exact" CS with disaggregation weights.
@@ -39,94 +35,136 @@ class ConditionalSpectrum {
* For each source/gmm, we need epsilon at the target IMT, Maths.epsilon(...)
* which we can obtain from the disaggregation
*
- * Primary loop is target, iterate over all other SAs including Tj = T* (shoud
- * just yield μCond = μ + ε and σ = 0.0
+ * Primary loop is target, iterate over all other SAs including Tj = T*
+ * (should just yield μCond = μ + ε and σ = 0.0
*/
- static final Table<Imt, Imt, Double> ÏTable = initCorrelations();
+ // static final Table<Imt, Imt, Double> ÏTable = initCorrelations();
/*
* Create a table of SA IMT correlations. Table contains duplicate values as
* all possible ordered pairs (e.g. Imt1:Imt2 and Imt2:Imt1) are included.
*/
- static Table<Imt, Imt, Double> initCorrelations() {
- BakerJayaram08 bj08 = new BakerJayaram08();
- Set<Imt> saImts = EnumSet.copyOf(Imt.mprsImts());
- saImts.remove(Imt.PGA);
- Table<Imt, Imt, Double> table = ArrayTable.create(saImts, saImts);
- for (Imt imt1 : saImts) {
- saImts.stream().forEach(imt2 -> table.put(
- imt1,
- imt2,
- Maths.round(bj08.applyAsDouble(imt1, imt2), 8)));
- }
- return ImmutableTable.copyOf(table);
+ // static Table<Imt, Imt, Double> initCorrelations() {
+ //
+ // Set<Imt> saImts = EnumSet.copyOf(Imt.mprsImts());
+ // saImts.remove(Imt.PGA);
+ // Table<Imt, Imt, Double> table = ArrayTable.create(saImts, saImts);
+ // for (Imt imt1 : saImts) {
+ // saImts.stream().forEach(imt2 -> table.put(
+ // imt1,
+ // imt2,
+ // Maths.round(CorrelationModel.BAKER_JAYARAM_08.value(imt1, imt2), 8)));
+ // }
+ // return ImmutableTable.copyOf(table);
+ // }
+
+ static double totalMean() {
+
+ // Summation over all μ * disaggWt
+ // double disaggWt =
+ return 0;
}
- /*
- * Compute the mean at the conditional IMT from the conditional SA mean, μ,
- * the conditional SA sigma, σ, the IMT correlation coefficient, Ï, and
- * epsilon, ε, at the target IMT.
+ /**
+ * @param μt conditional mean ground motion over all ruptures at a target
+ * period
+ * @param μc conditional mean ground motion of a rupture at a target period
+ * @param cσ conditional sigma of ground motion of a rupture at a target
+ * period
+ * @param p
+ * @return
*/
- private static double conditionalMean(double μ, double σ, double Ï, double ε) {
- return μ + Ï * ε * σ;
+ static double totalSigma(double μt, double μc, double σc, double p) {
+ double sigma = p * ((σc * σc) + (μc - μt) * (μc - μt));
+ // TODO summation of scaled sigmas over all ruptures
+ return sqrt(sigma);
}
/**
- * Baker and Jayaram (2008) intensity measure correlation model for spectral
- * accelerations between 0.01 and 10s.
+ * Compute the conditional mean ground motion of a source at some target
+ * period.
*
- * <p><b>Reference:</b> J. W. Baker and N. Jayaram. (2008). Correlation of
- * spectral acceleration values from NGA ground motion models, Earthquake
- * Spectra, 24: 1, 299-317.
+ * @param μ mean ground motion of a rupture at a target period
+ * @param σ sigma of ground motion of a rupture at a target period
+ * @param Ï correlation coefficient between reference and target period
+ * epsilons
+ * @param ε epsilon of a ground motion of a rupture at the reference period
*/
- static class BakerJayaram08 implements ImtCorrelation {
+ private static double mean(double μ, double σ, double Ï, double ε) {
+ return μ + Ï * ε * σ;
+ }
+
+ /**
+ * @param σ sigma of ground motion of a rupture at a target period
+ * @param Ï correlation coefficient between reference and target period
+ * epsilons
+ * @return
+ */
+ static double sigma(double σ, double Ï) {
+ return σ * sqrt(1 - Ï * Ï);
+ }
+
+ public static enum CorrelationModel {
+
+ BAKER_JAYARAM_08() {
+
+ @Override
+ public double value(Imt imt1, Imt imt2) {
+ checkArgument(imt1.isSA());
+ checkArgument(imt2.isSA());
+
+ double t1 = imt1.period();
+ double t2 = imt2.period();
+
+ double tMin = min(t1, t2);
+ double tMax = max(t1, t2);
+
+ double c1 = (1.0 - cos(PI / 2.0 - log(tMax / max(tMin, 0.109)) * 0.366));
+ double c2 = (tMax < 0.2)
+ ? 1.0 - 0.105 * (1.0 - 1.0 / (1.0 + exp(100.0 * tMax - 5.0))) *
+ (tMax - tMin) / (tMax - 0.0099)
+ : Double.NaN;
+ double c3 = (tMax < 0.109) ? c2 : c1;
+ double c4 = c1 + 0.5 * (sqrt(c3) - c3) * (1.0 + cos(PI * tMin / 0.109));
+
+ if (tMax <= 0.109) {
+ return c2;
+ } else if (tMin > 0.109) {
+ return c1;
+ } else if (tMax < 0.2) {
+ return min(c2, c4);
+ } else {
+ return c4;
+ }
+ }
+ };
/**
- * Computes the correlation coefficient between two spectral acceleration
- * IMTs. The order of the two IMTs does not matter.
+ * Compute the correlation coefficient between two IMTs. The order of the
+ * two IMTs does not matter. See individual models for supported IMTs.
*
* @param imt1 the first IMT
* @param imt2 the second IMT
- * @return pearson correlation coefficient between {@code imt1} and
+ * @return Pearson correlation coefficient between {@code imt1} and
* {@code imt2}
*/
- @Override
- public double applyAsDouble(Imt imt1, Imt imt2) {
-
- checkArgument(imt1.isSA());
- checkArgument(imt2.isSA());
-
- double t1 = imt1.period();
- double t2 = imt2.period();
-
- double tMin = min(t1, t2);
- double tMax = max(t1, t2);
-
- double c1 = (1.0 - cos(PI / 2.0 - log(tMax / max(tMin, 0.109)) * 0.366));
- double c2 = (tMax < 0.2)
- ? 1.0 - 0.105 * (1.0 - 1.0 / (1.0 + exp(100.0 * tMax - 5.0))) *
- (tMax - tMin) / (tMax - 0.0099)
- : Double.NaN;
- double c3 = (tMax < 0.109) ? c2 : c1;
- double c4 = c1 + 0.5 * (sqrt(c3) - c3) * (1.0 + cos(PI * tMin / 0.109));
-
- if (tMax <= 0.109) {
- return c2;
- } else if (tMin > 0.109) {
- return c1;
- } else if (tMax < 0.2) {
- return min(c2, c4);
- } else {
- return c4;
- }
+ public abstract double value(Imt imt1, Imt imt2);
+
+ /**
+ * Create a map of correlation coefficients between the supplied {@code imt}
+ * and {@code imts}.
+ */
+ public Map<Imt, Double> valueMap(Set<Imt> imts, Imt targetImt) {
+ Map<Imt, Double> coeffs = imts.stream().collect(toMap(
+ imt -> imt,
+ imt -> value(imt, targetImt)));
+ return new EnumMap<>(coeffs);
}
}
public static void main(String[] args) {
- BakerJayaram08 bj = new BakerJayaram08();
- // System.out.println(ÏTable);
- System.out.println(bj.applyAsDouble(Imt.SA0P02, Imt.SA0P02));
+ System.out.println(CorrelationModel.BAKER_JAYARAM_08.value(Imt.SA0P02, Imt.SA0P02));
}
}
--
GitLab