diff --git a/src/main/java/gov/usgs/earthquake/nshmp/calc/ExceedanceModel.java b/src/main/java/gov/usgs/earthquake/nshmp/calc/ExceedanceModel.java index 30d385a93195ab9c443f89577260cbeb0c645e45..7fadca4c5f8da488c93431cb3e82b0f4767aa6d4 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/calc/ExceedanceModel.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/calc/ExceedanceModel.java @@ -11,6 +11,7 @@ import java.util.List; import gov.usgs.earthquake.nshmp.Maths; import gov.usgs.earthquake.nshmp.data.MutableXySequence; +import gov.usgs.earthquake.nshmp.data.XyPoint; import gov.usgs.earthquake.nshmp.data.XySequence; import gov.usgs.earthquake.nshmp.gmm.Imt; import gov.usgs.earthquake.nshmp.gmm.MultiScalarGroundMotion; @@ -37,6 +38,10 @@ public enum ExceedanceModel { * TODO We probably want to refactor this to probability model and provide * 'occurrence' in addition to exceedence. See commented distribution function * at eof. + * + * Developer note: Exceedance methods can be called millions of times in a + * hazard calculation; do not be tempted to use streams to update XySequence + * values. */ /** @@ -52,8 +57,11 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { - sequence.stream().forEach(p -> p.set(Maths.stepFunction(μ, p.x()))); + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { + for (XyPoint p : sequence) { + p.set(Maths.stepFunction(μ, p.x())); + } return sequence; } }, @@ -70,7 +78,8 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { return boundedCcdFn(μ, σ, sequence, 0.0, 1.0); } }, @@ -87,7 +96,8 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { return boundedCcdFn(μ, σ, sequence, prob(μ, σ, n), 1.0); } }, @@ -105,7 +115,8 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { double pHi = prob(μ, σ, n); return boundedCcdFn(μ, σ, sequence, pHi, 1.0 - pHi); } @@ -123,7 +134,8 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { return Ccdfs.UPPER_3SIGMA.get(μ, σ, sequence); } }, @@ -141,7 +153,8 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { return boundedCcdFn(μ, 0.65, sequence, 0.0, 1.0); } }, @@ -164,8 +177,11 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { - sequence.stream().forEach(p -> p.set(exceedance(μ, σ, n, imt, p.x()))); + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { + for (XyPoint p : sequence) { + p.set(exceedance(μ, σ, n, imt, p.x())); + } return sequence; } }, @@ -185,7 +201,8 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { double pHi = prob(μ, σ, n, Math.log(ceusMaxValue(imt))); return boundedCcdFn(μ, σ, sequence, pHi, 1.0); } @@ -209,7 +226,8 @@ public enum ExceedanceModel { } @Override - XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence) { + MutableXySequence exceedance(double μ, double σ, double n, Imt imt, + MutableXySequence sequence) { double lnMaxGm = Math.log(ceusMaxValue(imt)); double ln3σGm = μ + 3.0 * σ; if (ln3σGm < lnMaxGm) { @@ -249,7 +267,12 @@ public enum ExceedanceModel { * {@link #NSHM_CEUS_MAX_INTENSITY} * @param value for which to compute the exceedance probability */ - abstract double exceedance(double μ, double σ, double n, Imt imt, double value); + abstract double exceedance( + double μ, + double σ, + double n, + Imt imt, + double value); /** * Compute the probability of exceeding a sequence of x-values. @@ -262,7 +285,12 @@ public enum ExceedanceModel { * @param sequence the x-values for which to compute exceedance probabilities * @return the supplied {@code sequence} */ - abstract XySequence exceedance(double μ, double σ, double n, Imt imt, XySequence sequence); + abstract MutableXySequence exceedance( + double μ, + double σ, + double n, + Imt imt, + MutableXySequence sequence); /* * Return a list of exceedance curves, one for each tree branch in the @@ -285,8 +313,7 @@ public enum ExceedanceModel { for (int i = 0; i < means.length; i++) { double μ = means[i]; for (int j = 0; j < sigmas.length; j++) { - curves.add(MutableXySequence.copyOf( - exceedance(μ, sigmas[j], n, imt, MutableXySequence.emptyCopyOf(sequence)))); + curves.add(exceedance(μ, sigmas[j], n, imt, MutableXySequence.emptyCopyOf(sequence))); } } return curves; @@ -343,14 +370,16 @@ public enum ExceedanceModel { * lower probability limits. Return the supplied {@code XySequence} populated * with probabilities. */ - private static XySequence boundedCcdFn( + private static MutableXySequence boundedCcdFn( double μ, double σ, - XySequence sequence, + MutableXySequence sequence, double pHi, double pLo) { - sequence.stream().forEach(p -> p.set(boundedCcdFn(μ, σ, p.x(), pHi, pLo))); + for (XyPoint p : sequence) { + p.set(boundedCcdFn(μ, σ, p.x(), pHi, pLo)); + } return sequence; } @@ -462,8 +491,10 @@ public enum ExceedanceModel { return 0.0; } - XySequence get(double μ, double σ, XySequence sequence) { - sequence.stream().forEach(p -> p.set(get(μ, σ, p.x()))); + MutableXySequence get(double μ, double σ, MutableXySequence sequence) { + for (XyPoint p : sequence) { + p.set(get(μ, σ, p.x())); + } return sequence; } }