diff --git a/.gitignore b/.gitignore
index 7933bef31eee328b596e12080caafe952c7e5e12..1d4aedc826db62eaeae637253d8d871be6c84dbe 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,6 +9,7 @@ tmp
.DS_Store
libs
HazardCalc.java
+DisaggCalc.java
node_modules
package-lock.json
package.json
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 dbe85e40e3579496b883f5cb3d91ca2f74b0541e..df8329cf53633ecf46519a3e6eea82c583e45fbb 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));
}
}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/calc/DisaggContributor.java b/src/main/java/gov/usgs/earthquake/nshmp/calc/DisaggContributor.java
index a97b29fde4904279934b0764e32cdffe4aa82059..928482fef82e593c16b15f15d730dda7db1a5577 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/calc/DisaggContributor.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/calc/DisaggContributor.java
@@ -3,14 +3,13 @@ package gov.usgs.earthquake.nshmp.calc;
import static com.google.common.base.Preconditions.checkState;
import static gov.usgs.earthquake.nshmp.Text.NEWLINE;
import static gov.usgs.earthquake.nshmp.model.SourceType.FAULT_SYSTEM;
+import static java.util.Comparator.comparingDouble;
import static java.util.stream.Collectors.toList;
import java.util.ArrayList;
import java.util.Collection;
-import java.util.Comparator;
import java.util.Iterator;
import java.util.List;
-import java.util.function.Function;
import gov.usgs.earthquake.nshmp.Maths;
import gov.usgs.earthquake.nshmp.calc.DisaggExport.ContributionFilter;
@@ -778,24 +777,9 @@ abstract class DisaggContributor {
private static List<DisaggContributor> buildAndSort(
Collection<DisaggContributor.Builder> builders) {
return builders.stream()
- .map(BUILDER::apply)
- .sorted(SORTER)
+ .map(Builder::build)
+ .sorted(comparingDouble(DisaggContributor::total).reversed())
.collect(toList());
- }
-
- static final Function<DisaggContributor.Builder, DisaggContributor> BUILDER =
- new Function<DisaggContributor.Builder, DisaggContributor>() {
- @Override
- public DisaggContributor apply(Builder builder) {
- return builder.build();
- }
- };
-
- static final Comparator<DisaggContributor> SORTER = new Comparator<DisaggContributor>() {
- @Override
- public int compare(DisaggContributor left, DisaggContributor right) {
- return Double.compare(right.total(), left.total());
- }
- };
+ }
}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/calc/DisaggDataset.java b/src/main/java/gov/usgs/earthquake/nshmp/calc/DisaggDataset.java
index be992d77bd6757636200657cdc1b7359249a6bc0..7edf1079bd066fa33a0508949007882356b66bf4 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/calc/DisaggDataset.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/calc/DisaggDataset.java
@@ -2,6 +2,7 @@ package gov.usgs.earthquake.nshmp.calc;
import static com.google.common.base.Preconditions.checkState;
import static gov.usgs.earthquake.nshmp.data.DoubleData.checkInRange;
+import static java.util.Comparator.comparingDouble;
import static java.util.stream.Collectors.toList;
import java.util.ArrayList;
@@ -438,7 +439,7 @@ final class DisaggDataset {
DisaggDataset build() {
List<DisaggContributor> sortedList = contributors.stream()
- .sorted(DisaggContributor.SORTER)
+ .sorted(comparingDouble(DisaggContributor::total).reversed())
.collect(toList());
return super.build(sortedList);
}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/calc/HazardExport.java b/src/main/java/gov/usgs/earthquake/nshmp/calc/HazardExport.java
index 47ccd92e82dc3bc96579438e7745172c9763a5d1..8a62bce2c709d3a6813632596cedba6a4118b184 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/calc/HazardExport.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/calc/HazardExport.java
@@ -69,7 +69,7 @@ public final class HazardExport {
private HazardExport(
HazardModel model,
CalcConfig config,
- List<Site> sites,
+ boolean namedSites,
Path out) throws IOException {
this.out = out;
@@ -83,10 +83,8 @@ public final class HazardExport {
? Maths.annualRateToProbabilityConverter().andThen(formatter)
: formatter;
- Site demoSite = sites.get(0);
- this.namedSites = demoSite.name() != Site.NO_NAME;
-
- init(sites);
+ this.namedSites = namedSites;
+ init();
}
/**
@@ -94,21 +92,22 @@ public final class HazardExport {
*
* @param model being run
* @param config that specifies output options and formats
- * @param sites reference to the sites to be processed (not retained)
+ * @param namedSites whether the sites being processed contain a 'name' column
+ * that should be preserved in output files
* @param out results directory; this may have been been modified from
* {@code config.output.directory} by calling program
*/
public static HazardExport create(
HazardModel model,
CalcConfig config,
- List<Site> sites,
+ boolean namedSites,
Path out) throws IOException {
- return new HazardExport(model, config, sites, out);
+ return new HazardExport(model, config, namedSites, out);
}
/* Initialize output directories. */
- private void init(List<Site> sites) throws IOException {
+ private void init() throws IOException {
for (Imt imt : config.hazard.imts) {
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/GroundMotion.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/GroundMotion.java
index 2647d3d86dbfebc5bbfe252c10832dfca6d7a67d..76dff8be232ceee20d4bfd71e081457474133876 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/GroundMotion.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/GroundMotion.java
@@ -19,6 +19,12 @@ public interface GroundMotion {
*/
public double sigma();
+ // public double rate();
+ //
+ // public double epsilon();
+ //
+ // public void setDisaggData(double rate, double epsilon);
+
/**
* Create a new ground motion container.
*
@@ -29,6 +35,9 @@ public interface GroundMotion {
public static GroundMotion create(double mean, double sigma) {
return new GroundMotion() {
+ // private double rate;
+ // private double epsilon;
+
@Override
public double mean() {
return mean;
@@ -39,6 +48,22 @@ public interface GroundMotion {
return sigma;
}
+ // @Override
+ // public double rate() {
+ // return rate;
+ // }
+ //
+ // @Override
+ // public double epsilon() {
+ // return epsilon;
+ // }
+ //
+ // @Override
+ // public void setDisaggData(double rate, double epsilon) {
+ // this.rate = rate;
+ // this.epsilon = epsilon;
+ // }
+
@Override
public String toString() {
return new StringBuilder(getClass().getSimpleName())
@@ -49,6 +74,7 @@ public interface GroundMotion {
.append("]")
.toString();
}
+
};
}
}
diff --git a/src/test/java/gov/usgs/earthquake/nshmp/model/LoaderTests.java b/src/test/java/gov/usgs/earthquake/nshmp/model/LoaderTests.java
index 0b79509351ee2bb23c1705a18ec4e8230811faa3..aca6c3d41cd0b1e20becc6f8b3082d367c8f019b 100644
--- a/src/test/java/gov/usgs/earthquake/nshmp/model/LoaderTests.java
+++ b/src/test/java/gov/usgs/earthquake/nshmp/model/LoaderTests.java
@@ -214,7 +214,8 @@ class LoaderTests {
List<Site> sites) throws IOException {
ExecutorService exec = initExecutor(config.performance.threadCount);
- HazardExport handler = HazardExport.create(model, config, sites, config.output.directory);
+ boolean namedSites = sites.get(0).name() != Site.NO_NAME;
+ HazardExport handler = HazardExport.create(model, config, namedSites, config.output.directory);
for (Site site : sites) {
Hazard hazard = HazardCalcs.hazard(model, config, site, exec);
handler.write(hazard);