diff --git a/src/main/java/scratch/kevin/nga/Idriss2013Debug.java b/src/main/java/scratch/kevin/nga/Idriss2013Debug.java index dc8c2358..3c31e34b 100644 --- a/src/main/java/scratch/kevin/nga/Idriss2013Debug.java +++ b/src/main/java/scratch/kevin/nga/Idriss2013Debug.java @@ -10,10 +10,11 @@ import org.opensha.sha.earthquake.param.BackgroundRupType; import org.opensha.sha.earthquake.param.ProbabilityModelOptions; import org.opensha.sha.earthquake.rupForecastImpl.FaultRuptureSource; -import org.opensha.sha.earthquake.rupForecastImpl.PointSource13b; +import org.opensha.sha.earthquake.rupForecastImpl.PointSourceNshm; import org.opensha.sha.earthquake.rupForecastImpl.WGCEP_UCERF_2_Final.griddedSeis.Point2Vert_FaultPoisSource; import org.opensha.sha.faultSurface.EvenlyGriddedSurface; import org.opensha.sha.faultSurface.FourPointEvenlyGriddedSurface; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrections; import org.opensha.sha.imr.attenRelImpl.ngaw2.NGAW2_Wrappers.Idriss_2014_Wrapper; import org.opensha.sha.imr.param.IntensityMeasureParams.PGA_Param; import org.opensha.sha.imr.param.SiteParams.DepthTo1pt0kmPerSecParam; @@ -66,7 +67,6 @@ public static void main(String[] args) { double fracReverse = 0d; WC1994_MagLengthRelationship magLenRel = new WC1994_MagLengthRelationship(); double ptSrcCutoff = 6.0; - double[] DEPTHS = new double[] {5.0, 1.0}; ProbEqkSource pointSource; @@ -74,19 +74,19 @@ public static void main(String[] args) { case CROSSHAIR: pointSource = new Point2Vert_FaultPoisSource(faultLoc, mfd, magLenRel, duration, ptSrcCutoff, fracStrikeSlip, fracNormal, - fracReverse, true); + fracReverse, true, null); break; case FINITE: pointSource = new Point2Vert_FaultPoisSource(faultLoc, mfd, magLenRel, duration, ptSrcCutoff, fracStrikeSlip, fracNormal, - fracReverse, false); + fracReverse, false, null); break; case POINT: Map mechMap = Maps.newHashMap(); mechMap.put(FocalMech.STRIKE_SLIP, fracStrikeSlip); mechMap.put(FocalMech.REVERSE, fracReverse); mechMap.put(FocalMech.NORMAL, fracNormal); - pointSource = new PointSource13b(faultLoc, mfd, duration, DEPTHS, mechMap); + pointSource = new PointSourceNshm(faultLoc, mfd, duration, mechMap, PointSourceDistanceCorrections.NONE.get()); break; default: throw new IllegalStateException("Unknown Background Rup Type: "+bgRupType); diff --git a/src/main/java/scratch/nshm23/BValSweepHazardComparison.java b/src/main/java/scratch/kevin/nshm23/BValSweepHazardComparison.java similarity index 99% rename from src/main/java/scratch/nshm23/BValSweepHazardComparison.java rename to src/main/java/scratch/kevin/nshm23/BValSweepHazardComparison.java index cc4ac896..0203593a 100644 --- a/src/main/java/scratch/nshm23/BValSweepHazardComparison.java +++ b/src/main/java/scratch/kevin/nshm23/BValSweepHazardComparison.java @@ -1,4 +1,4 @@ -package scratch.nshm23; +package scratch.kevin.nshm23; import java.awt.Color; import java.io.File; diff --git a/src/main/java/scratch/kevin/nshm23/BranchAveragedHazardScriptWriter.java b/src/main/java/scratch/kevin/nshm23/BranchAveragedHazardScriptWriter.java index 33e94cdd..0a9bad8a 100644 --- a/src/main/java/scratch/kevin/nshm23/BranchAveragedHazardScriptWriter.java +++ b/src/main/java/scratch/kevin/nshm23/BranchAveragedHazardScriptWriter.java @@ -33,6 +33,7 @@ public static void main(String[] args) throws IOException { boolean linkFromBase = true; Double vs30 = null; double gridSpacing = 0.1; + boolean supersample = false; double[] periods = { 0d, 0.2d, 1d, 5d }; AttenRelRef[] gmms = null; @@ -77,7 +78,9 @@ public static void main(String[] args) throws IOException { // String solFileName = "results_PRVI_SUB_FM_LARGE_branch_averaged_gridded.zip"; String baseDirName = "2024_11_19-prvi25_crustal_subduction_combined_branches"; - String suffix = "ba_only"; +// String suffix = "ba_only"; +// String suffix = "ba_only-test_grid_branch"; + String suffix = "ba_only-test_grid_branch-supersample"; supersample = true; String solFileName = "combined_branch_averaged_solution.zip"; // String baseDirName = "2024_11_19-prvi25_crustal_branches-dmSample5x"; @@ -158,6 +161,9 @@ public static void main(String[] args) throws IOException { // remoteTotalMemGB*1024, null); // BatchScriptWriter pbsWrite = new HovenweepScriptWriter(); + if (parallelMPJWrite instanceof FastMPJShellScriptWriter) + ((FastMPJShellScriptWriter)parallelMPJWrite).setUseLaunchWrapper(true); + parallelMPJWrite.setEnvVar("MAIN_DIR", remoteMainDir.getAbsolutePath()); singleMPJWrite.setEnvVar("MAIN_DIR", remoteMainDir.getAbsolutePath()); String mainDirPath = "$MAIN_DIR"; @@ -239,6 +245,8 @@ else if (gridReg.getNodeCount() > 5000) argz += (float)periods[p]; } } + if (supersample) + argz += " --supersample"; argz += " "+dispatchArgs; File jobFile = new File(localDir, "batch_hazard_"+bgOp.name()+".slurm"); diff --git a/src/main/java/scratch/kevin/nshm23/SectRIDistCSVWriter.java b/src/main/java/scratch/kevin/nshm23/SectRIDistCSVWriter.java new file mode 100644 index 00000000..305bdd66 --- /dev/null +++ b/src/main/java/scratch/kevin/nshm23/SectRIDistCSVWriter.java @@ -0,0 +1,123 @@ +package scratch.kevin.nshm23; + +import java.awt.geom.Point2D; +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; +import java.util.List; + +import org.opensha.commons.data.CSVFile; +import org.opensha.commons.data.function.ArbDiscrEmpiricalDistFunc; +import org.opensha.commons.data.function.EvenlyDiscretizedFunc; +import org.opensha.sha.earthquake.faultSysSolution.FaultSystemRupSet; +import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; +import org.opensha.sha.earthquake.faultSysSolution.modules.BranchSectParticMFDs; +import org.opensha.sha.earthquake.faultSysSolution.util.FaultSysTools; +import org.opensha.sha.faultSurface.FaultSection; +import org.opensha.sha.magdist.IncrementalMagFreqDist; + +public class SectRIDistCSVWriter { + + public static void main(String[] args) throws IOException { + FaultSystemSolution sol = FaultSystemSolution.load(new File("/home/kevin/OpenSHA/nshm23/batch_inversions/" + + "2024_02_02-nshm23_branches-WUS_FM_v3/results_WUS_FM_v3_branch_averaged.zip")); + FaultSystemRupSet rupSet = sol.getRupSet(); + + BranchSectParticMFDs sectMFDs = sol.requireModule(BranchSectParticMFDs.class); + + double[] mags = {0d, 6d, 6.5d, 7d, 7.5d}; + + double[] fractiles = {0.025,0.16, 0.84, 0.975}; + + List riHeader = List.of("Subsection Index", "Parent Section ID", "Subsection Name" + ,"Mean RI", "Min RI", "Max RI", "Std. Dev", "p2.5", "p16.0", "p84.0", "p97.5"); + List rateHeader = new ArrayList<>(riHeader); + for (int i=0; i> riCSVs = new ArrayList<>(mags.length); + List> rateCSVs = new ArrayList<>(mags.length); + for (int m=0; m riCSV = new CSVFile<>(true); + riCSV.addLine(riHeader); + riCSVs.add(riCSV); + CSVFile rateCSV = new CSVFile<>(true); + rateCSV.addLine(rateHeader); + rateCSVs.add(rateCSV); + } + + EvenlyDiscretizedFunc refMFD = FaultSysTools.initEmptyMFD(5.05, 8.45); + List mfdHeader = new ArrayList<>(); + mfdHeader.addAll(riHeader.subList(0, 3)); + for (int i=0; i mfdParticCSV = new CSVFile<>(true); + mfdParticCSV.addLine(mfdHeader); + CSVFile mfdNuclCSV = new CSVFile<>(true); + mfdNuclCSV.addLine(mfdHeader); + + for (int s=0; s common = List.of(s+"", sect.getParentSectionId()+"", sect.getSectionName()); + List partic = new ArrayList<>(common); + List nucl = new ArrayList<>(common); + for (int i=0; i= (float)mags[m]) + rates[m] += pt.getY(); + if (!rate) + for (int m=0; m line = new ArrayList<>(common); + line.add((float)dists[m].getMean()+""); + line.add((float)dists[m].getMinX()+""); + line.add((float)dists[m].getMinY()+""); + line.add((float)dists[m].getStdDev()+""); + for (double fract : fractiles) + line.add((float)dists[m].getInterpolatedFractile(fract)+""); + if (rate) + rateCSVs.get(m).addLine(line); + else + riCSVs.get(m).addLine(line); + } + } + } + + File outputDir = new File("/tmp"); + String prefix = "nshm23_wus"; + mfdParticCSV.writeToFile(new File(outputDir, prefix+"_partic_mfds.csv")); + mfdNuclCSV.writeToFile(new File(outputDir, prefix+"_nucl_mfds.csv")); + for (int m=0; m trts = EnumSet.of(TectonicRegionType.ACTIVE_SHALLOW); if (subduction) { trts.add(TectonicRegionType.SUBDUCTION_INTERFACE); @@ -144,7 +149,7 @@ public static void main(String[] args) throws IOException { calc.getHazardCurve(logXVals, site, gmpe, erf); DiscretizedFunc wrapperFaultCurve = toLinear(logXVals, xVals); - ProbEqkSource modelTestSrc = gridProv.getSource(testIndex, 1d, null, BackgroundRupType.POINT); + ProbEqkSource modelTestSrc = gridProv.getSource(testIndex, 1d, null, gridSettings); FaultSystemSolutionERF modelERF = new FaultSystemSolutionERF(baSol); modelERF.setParameter(ApplyGardnerKnopoffAftershockFilterParam.NAME, false); modelERF.setParameter(ProbabilityModelParam.NAME, ProbabilityModelOptions.POISSON); diff --git a/src/main/java/scratch/kevin/pointSources/PointSourceHazardComparison.java b/src/main/java/scratch/kevin/pointSources/PointSourceHazardComparison.java index 456946d5..1d460c8e 100644 --- a/src/main/java/scratch/kevin/pointSources/PointSourceHazardComparison.java +++ b/src/main/java/scratch/kevin/pointSources/PointSourceHazardComparison.java @@ -31,6 +31,8 @@ import org.opensha.commons.calc.magScalingRelations.magScalingRelImpl.WC1994_MagLengthRelationship; import org.opensha.commons.data.CSVFile; import org.opensha.commons.data.Site; +import org.opensha.commons.data.WeightedList; +import org.opensha.commons.data.WeightedValue; import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc; import org.opensha.commons.data.function.DefaultXY_DataSet; import org.opensha.commons.data.function.DiscretizedFunc; @@ -59,6 +61,7 @@ import org.opensha.sha.calc.HazardCurveCalculator; import org.opensha.sha.earthquake.AbstractERF; import org.opensha.sha.earthquake.FocalMechanism; +import org.opensha.sha.earthquake.PointSource; import org.opensha.sha.earthquake.ProbEqkRupture; import org.opensha.sha.earthquake.ProbEqkSource; import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; @@ -66,20 +69,16 @@ import org.opensha.sha.earthquake.faultSysSolution.modules.MFDGridSourceProvider; import org.opensha.sha.earthquake.faultSysSolution.util.FaultSysTools; import org.opensha.sha.earthquake.faultSysSolution.util.SolHazardMapCalc.ReturnPeriods; -import org.opensha.sha.earthquake.rupForecastImpl.PointEqkSource; -import org.opensha.sha.earthquake.rupForecastImpl.PointSource13b; import org.opensha.sha.earthquake.rupForecastImpl.PointSourceNshm; -import org.opensha.sha.earthquake.rupForecastImpl.PointSourceNshm.PointSurfaceNshm; import org.opensha.sha.earthquake.rupForecastImpl.PointToFiniteSource; import org.opensha.sha.earthquake.rupForecastImpl.nshm23.gridded.NSHM23_AbstractGridSourceProvider; -import org.opensha.sha.earthquake.rupForecastImpl.nshm23.logicTree.NSHM23_SingleStates; -import org.opensha.sha.faultSurface.EvenlyGriddedSurface; import org.opensha.sha.faultSurface.FiniteApproxPointSurface; import org.opensha.sha.faultSurface.PointSurface; import org.opensha.sha.faultSurface.QuadSurface; import org.opensha.sha.faultSurface.RuptureSurface; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrection; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrections; import org.opensha.sha.faultSurface.utils.PointSurfaceBuilder; -import org.opensha.sha.faultSurface.utils.PtSrcDistCorr.Type; import org.opensha.sha.gui.infoTools.IMT_Info; import org.opensha.sha.imr.AttenRelRef; import org.opensha.sha.imr.ScalarIMR; @@ -106,36 +105,26 @@ public class PointSourceHazardComparison { enum PointSourceType { - TRUE_POINT_SOURCE("True Point Source", false, Type.NONE) { + TRUE_POINT_SOURCE("True Point Source", false, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return new PointEqkSource(centerLoc, mfd, 1d, aveRake, aveDip); +// return new PointEqkSource(centerLoc, mfd, 1d, aveRake, aveDip); + return PointSource.poissonBuilder(centerLoc) + .truePointSources() + .forMFDAndFocalMech(mfd, new FocalMechanism(Double.NaN, aveDip, aveRake)) + .duration(1d) + .build(); } }, - POINT_SOURCE_13b_NO_CORR("PointSource13b (no distance correction)", false, Type.NONE) { + POINT_SOURCE_NSHM("PointSourceNshm", false, PointSourceDistanceCorrections.NSHM_2013) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return new PointSource13b(centerLoc, mfd, 1d, new double[] {5.0, 1.0}, mechWtMapForRake(aveRake)); + return new PointSourceNshm(centerLoc, mfd, 1d, mechWtMapForRake(aveRake), getDistCorrs()); } }, - POINT_SOURCE_13b_NSHMP_CORR("PointSource13b (NSHMP08 distance correction)", false, Type.NSHMP08) { - @Override - public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, - double aveRake, double aveDip, boolean isSupersample, Random r) { - return new CorrOverrideProbEqkSource(POINT_SOURCE_13b_NO_CORR.buildSource( - centerLoc, mfd, aveRake, aveDip, isSupersample, r), Type.NSHMP08); - } - }, - POINT_SOURCE_NSHM("PointSourceNshm", false, Type.NONE) { // none is fine here, it's baked in - @Override - public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, - double aveRake, double aveDip, boolean isSupersample, Random r) { - return new PointSourceNshm(centerLoc, mfd, 1d, mechWtMapForRake(aveRake)); - } - }, - APROX_SUPERSAMPLE_POINT_SOURCE_NSHM("Approx. Supersampled PointSourceNshm", false, Type.NONE) { // none is fine here, it's baked in + APROX_SUPERSAMPLE_POINT_SOURCE_NSHM("Approx. Supersampled PointSourceNshm", false, PointSourceDistanceCorrections.NONE) { // none is fine here, it's baked in @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -144,7 +133,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, // not externally supersampled, approximate it here List modRups = new ArrayList<>(); for (ProbEqkRupture rup : source) { - PointSurfaceNshm ptSurf = (PointSurfaceNshm)rup.getRuptureSurface(); + FiniteApproxPointSurface ptSurf = (FiniteApproxPointSurface)rup.getRuptureSurface(); rup.setRuptureSurface(new SupersampledApproxPointSurfaceNshm(ptSurf, rup.getMag(), 0.1, SUPERSAMPLE_SPACING)); modRups.add(rup); } @@ -153,7 +142,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return source; } }, - POINT_TO_FINITE("PointToFiniteSource", true, Type.NONE) { + POINT_TO_FINITE("PointToFiniteSource", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -181,21 +170,21 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return new RupListSource(rups, source0); } }, - SIMPLE_APPROX_FINITE_POINT_NO_CORR("Simple Approx. Finite, No Dist. Corr.", true, Type.NONE) { + SIMPLE_APPROX_FINITE_POINT_NO_CORR("Simple Approx. Finite, No Dist. Corr.", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildPointSource(centerLoc, null, mfd, aveRake, aveDip, r, getDistCorrType()); + return buildPointSource(centerLoc, null, mfd, aveRake, aveDip, r, getDistCorrs()); } }, - SIMPLE_APPROX_FINITE_POINT_NSHMP_CORR("Simple Approx. Finite, NSHMP08 Dist. Corr.", true, Type.NSHMP08) { + SIMPLE_APPROX_FINITE_POINT_NSHMP_CORR("Simple Approx. Finite, NSHMP08 Dist. Corr.", true, PointSourceDistanceCorrections.NSHM_2008) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildPointSource(centerLoc, null, mfd, aveRake, aveDip, r, getDistCorrType()); + return buildPointSource(centerLoc, null, mfd, aveRake, aveDip, r, getDistCorrs()); } }, - SINGLE_QUAD_TRACE_CENTERED("Single Quad Surface, Trace Centered", true, Type.NONE) { + SINGLE_QUAD_TRACE_CENTERED("Single Quad Surface, Trace Centered", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -206,7 +195,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return buildQuadSource(builder, mfd, aveRake, aveDip, 1); } }, - SINGLE_QUAD_SURF_CENTER("Single Quad Surface", true, Type.NONE) { + SINGLE_QUAD_SURF_CENTER("Single Quad Surface", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -217,7 +206,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return buildQuadSource(builder, mfd, aveRake, aveDip, 1); } }, - CROSSHAIR_QUAD("Crosshair Quad Surface", true, Type.NONE) { + CROSSHAIR_QUAD("Crosshair Quad Surface", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -228,7 +217,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return buildQuadSource(builder, mfd, aveRake, aveDip, 2); } }, - QUAD_QUAD("4x Quad Surface", true, Type.NONE) { + QUAD_QUAD("4x Quad Surface", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -239,7 +228,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 4); } }, - OCT_QUAD("8x Quad Surface", true, Type.NONE) { + OCT_QUAD("8x Quad Surface", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -250,7 +239,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 8); } }, - OCT_QUAD_RAND_DAS_DD("8x Quad Surface, Random DAS & DD", true, Type.NONE) { + OCT_QUAD_RAND_DAS_DD("8x Quad Surface, Random DAS & DD", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -261,7 +250,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 8); } }, - OCT_QUAD_RAND_CELL("8x Quad Surface, Random Cell Locations", true, Type.NONE) { + OCT_QUAD_RAND_CELL("8x Quad Surface, Random Cell Locations", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -274,7 +263,7 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 8); } }, - OCT_QUAD_RAND_DAS_DD_CELL("8x Quad Surface, Random DAS, DD, & Cell Locations", true, Type.NONE) { + OCT_QUAD_RAND_DAS_DD_CELL("8x Quad Surface, Random DAS, DD, & Cell Locations", true, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -287,35 +276,35 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 8); } }, - TABLE_FINITE_APPROX_POINT_SOURCE("Table-Based Finite Approx", false, Type.NONE) { + TABLE_FINITE_APPROX_POINT_SOURCE("Table-Based Finite Approx", false, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { return buildTableBackedSource(centerLoc, mfd, aveRake, aveDip, false, r); } }, - TABLE_FINITE_APPROX_POINT_SOURCE_SUPERSAMPLE("Table-Based Finite Supersample Approx", false, Type.NONE) { + TABLE_FINITE_APPROX_POINT_SOURCE_SUPERSAMPLE("Table-Based Finite Supersample Approx", false, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { return buildTableBackedSource(centerLoc, mfd, aveRake, aveDip, isSupersample ? false : true, r); } }, - HAZ_EQUIV_APPROX_POINT_SOURCE("Hazard Equivalent Finite Approx", false, Type.NONE) { + HAZ_EQUIV_APPROX_POINT_SOURCE("Hazard Equivalent Finite Approx", false, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { return buildHazEquivSource(centerLoc, mfd, aveRake, aveDip, this); } }, - CDF_BASED_RJB_CORR_APPROX_POINT_SOURCE("RJB Distribution Approx Finite Point Source", false, Type.NONE) { + CDF_BASED_RJB_CORR_APPROX_POINT_SOURCE("RJB Distribution Approx Finite Point Source", false, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { return buildRJBCorrSource(centerLoc, mfd, aveRake, aveDip, OCT_QUAD, 5); } }, - CDF_BASED_RJB_CORR_APPROX_SUPERSAMPLED_POINT_SOURCE("RJB Distribution Approx Supersampled Finite Point Source", false, Type.NONE) { + CDF_BASED_RJB_CORR_APPROX_SUPERSAMPLED_POINT_SOURCE("RJB Distribution Approx Supersampled Finite Point Source", false, PointSourceDistanceCorrections.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { @@ -325,17 +314,21 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, private String label; private boolean stochastic; - private Type distCorrType; + private WeightedList distCorrs; + + private PointSourceType(String label, boolean stochastic, PointSourceDistanceCorrections distCorrType) { + this(label, stochastic, distCorrType.get()); + } - private PointSourceType(String label, boolean stochastic, Type distCorrType) { + private PointSourceType(String label, boolean stochastic, WeightedList distCorrs) { this.label = label; this.stochastic = stochastic; - Preconditions.checkNotNull(distCorrType); - this.distCorrType = distCorrType; + Preconditions.checkNotNull(distCorrs); + this.distCorrs = distCorrs; } - public Type getDistCorrType() { - return distCorrType; + public WeightedList getDistCorrs() { + return distCorrs; } public abstract ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, @@ -405,14 +398,18 @@ private static ProbEqkSource buildQuadSource(PointSurfaceBuilder builder, Increm } private static ProbEqkSource buildPointSource(Location centerLoc, Region cell, IncrementalMagFreqDist mfd, - double aveRake, double aveDip, Random r, Type corrType) { + double aveRake, double aveDip, Random r, WeightedList distCorrs) { PointSurfaceBuilder builder = new PointSurfaceBuilder(centerLoc); builder.dip(aveDip); builder.random(r); builder.sampleFromCell(cell); double dipRad = Math.toRadians(aveDip); - boolean[] footwalls = aveDip == 90 ? new boolean[] {true} : new boolean[] {false, true}; - List rups = new ArrayList<>(mfd.size()*footwalls.length); + int expectedSize = mfd.size(); + if (aveDip != 90d) + expectedSize *= 2; + if (distCorrs != null) + expectedSize *= distCorrs.size(); + List rups = new ArrayList<>(expectedSize); for (int i=0; i surfs = builder.buildFiniteApproxPointSurfaces(distCorrs); - double rateEach = mfd.getY(i)/(double)footwalls.length; - double probEach = 1-Math.exp(-rateEach); - for (boolean footwall : footwalls) { - builder.footwall(footwall); - FiniteApproxPointSurface surf = builder.buildFiniteApproxPointSurface(); - surf.setDistCorrMagAndType(mag, corrType); - rups.add(new ProbEqkRupture(mag, aveRake, probEach, surf, null)); + for (WeightedValue surfWeight : surfs) { + double myRate = totRate*surfWeight.weight; + double myProb = 1-Math.exp(-myRate); + ProbEqkRupture rup = new ProbEqkRupture(mag, aveRake, myProb, surfWeight.value, null); + rups.add(rup); } } return new RupListSource(rups, null); @@ -671,11 +668,12 @@ public ProbEqkRupture getRupture(int nRupture) { private static class CorrOverrideProbEqkSource extends ProbEqkSource { private ProbEqkSource source; - private Type corrType; + private PointSourceDistanceCorrection distCorr; - public CorrOverrideProbEqkSource(ProbEqkSource source, Type corrType) { + public CorrOverrideProbEqkSource(ProbEqkSource source, WeightedList distCorrs) { this.source = source; - this.corrType = corrType; + Preconditions.checkState(distCorrs.size() == 1); + this.distCorr = distCorrs.getValue(0); } @Override @@ -702,7 +700,7 @@ public int getNumRuptures() { public ProbEqkRupture getRupture(int nRupture) { ProbEqkRupture rup = source.getRupture(nRupture); Preconditions.checkState(rup.getRuptureSurface() instanceof PointSurface); - ((PointSurface)rup.getRuptureSurface()).setDistCorrMagAndType(rup.getMag(), corrType); + ((PointSurface)rup.getRuptureSurface()).setDistanceCorrection(distCorr, rup); return rup; } @@ -710,13 +708,11 @@ public ProbEqkRupture getRupture(int nRupture) { static class PointSourceCalcERF extends AbstractERF { - private Type distCorrType; private List sources; private List> sourceCalls; public PointSourceCalcERF(PointSourceType type, Location loc, IncrementalMagFreqDist mfd, double rake, double dip, int numPerStochastic, Random r) { - this.distCorrType = type.distCorrType; if (type.stochastic) { sources = new ArrayList<>(numPerStochastic); IncrementalMagFreqDist mfdEach = mfd.deepClone(); @@ -730,7 +726,6 @@ public PointSourceCalcERF(PointSourceType type, Location loc, IncrementalMagFreq public PointSourceCalcERF(PointSourceType type, GriddedRegion locs, IncrementalMagFreqDist mfd, double rake, double dip, int numPerStochastic, Random r) { - this.distCorrType = type.distCorrType; IncrementalMagFreqDist mfdEach = mfd.deepClone(); sources = new ArrayList<>(type.stochastic ? numPerStochastic*locs.getNodeCount() : locs.getNodeCount()); mfdEach.scale(1d/locs.getNodeCount()); @@ -749,7 +744,6 @@ public PointSourceCalcERF(PointSourceType type, GriddedRegion locs, IncrementalM public PointSourceCalcERF(PointSourceType type, GridSourceProvider gridProv, double supersampleDisc, int numPerStochastic, Random r, Region calcRegion, double maxDist) { sourceCalls = new ArrayList<>(); - this.distCorrType = type.distCorrType; GriddedRegion gridReg = gridProv.getGriddedRegion(); for (int i=0; i testCases = new ArrayList<>(); + + GriddedRegion superDuperSampled = new GriddedRegion(cell, 0.0001, new Location(0.00005, 0.00005)); + System.out.println("SuperDuperSampled has "+superDuperSampled.getNodeCount()+" nodes"); + TestCase refCase = new LocListTestCase(superDuperSampled.getNodeList(), "Ideal sample"); + + testCases.add(new GriddedRegionTestCase(cell, 21)); + testCases.add(new GriddedRegionTestCase(cell, 20)); + testCases.add(new GriddedRegionTestCase(cell, 12)); + testCases.add(new GriddedRegionTestCase(cell, 11)); + testCases.add(new GriddedRegionTestCase(cell, 10)); + testCases.add(new GriddedRegionTestCase(cell, 6)); + testCases.add(new GriddedRegionTestCase(cell, 5)); + + testCases.add(new ExteriorNodesTestCase(cell, 12)); + testCases.add(new ExteriorNodesTestCase(cell, 11)); + testCases.add(new ExteriorNodesTestCase(cell, 10)); + testCases.add(new ExteriorNodesTestCase(cell, 6)); + testCases.add(new ExteriorNodesTestCase(cell, 5)); + + LocationList border = cell.getBorder(); + if (LocationUtils.areSimilar(border.first(), border.last())) { + // remove the closing point (don't want any duplicates) + border = new LocationList(border); + border.remove(border.size()-1); + } + + testCases.add(new LocListTestCase(border, "Corners")); + + testCases.add(new ResampledBorderTestCase(border, 1d)); + testCases.add(new ResampledBorderTestCase(border, 2d)); + testCases.add(new ResampledBorderTestCase(border, 5d)); + + LocationList sideCenterBorder = new LocationList(); + sideCenterBorder.add(new Location(centerLoc.getLatitude(), cell.getMaxLon())); + sideCenterBorder.add(new Location(centerLoc.getLatitude(), cell.getMinLon())); + sideCenterBorder.add(new Location(cell.getMaxLat(), centerLoc.getLongitude())); + sideCenterBorder.add(new Location(cell.getMinLat(), centerLoc.getLongitude())); + testCases.add(new LocListTestCase(sideCenterBorder, "Side centers")); + LocationList centersAndCorners = new LocationList(); + centersAndCorners.addAll(sideCenterBorder); + centersAndCorners.addAll(border); + testCases.add(new LocListTestCase(centersAndCorners, "Corners and centers")); + + LocationList centerAndCorners = new LocationList(); + centerAndCorners.add(centerLoc); + centerAndCorners.addAll(border); + testCases.add(new LocListTestCase(centerAndCorners, "Corners and center")); + + LocationList centerOnly = new LocationList(); + centerOnly.add(centerLoc); + testCases.add(new LocListTestCase(centerOnly, "Center only")); + + double[] dists = { + 0d, 10d, 20d, 30d, 40d, 50d, 75d, 100d, 150d, 200d, 300d + }; + + int longestName = 0; + for (TestCase test : testCases) + longestName = Integer.max(longestName, test.getName().length()); + + for (double dist : dists) { + System.out.println("Center Distance: "+(float)dist); + int numTests = dist == 0d ? 1 : 50; + Location[] testLocs = new Location[numTests]; + if (numTests == 1) { + testLocs[0] = centerLoc; + } else { + for (int i=0; i= 0d) + ret += "+"; + ret += pDF.format(diff/ref.getAverage())+")"; + } + return ret; + } + + private static class LocListTestCase implements TestCase { + + protected LocationList locs; + private String name; + + public LocListTestCase(LocationList locs, String name) { + this.locs = locs; + this.name = name; + } + + @Override + public String getName() { + return name+" ["+locs.size()+"]"; + } + + @Override + public MinMaxAveTracker run(Location centerLoc, Location[] testLocs) { + MinMaxAveTracker track = new MinMaxAveTracker(); + for (Location testLoc : testLocs) + for (Location loc : locs) + track.addValue(LocationUtils.linearDistanceFast(testLoc, loc)); + return track; + } + + } + + private static final DecimalFormat oDF = new DecimalFormat("0.##"); + + private static class ResampledBorderTestCase extends LocListTestCase { + + public ResampledBorderTestCase(LocationList border, double discrKM) { + super(null, oDF.format(discrKM)+"km border"); + if (LocationUtils.areSimilar(border.first(), border.last())) { + // remove the closing point (don't want any duplicates) + border = new LocationList(border); + border.remove(border.size()-1); + } + FaultTrace borderAsTrace = new FaultTrace(null); + borderAsTrace.addAll(border); + borderAsTrace.add(border.first()); // need to close it, will remove at the end + int samples = (int)Math.round(borderAsTrace.getTraceLength()/discrKM); + locs = FaultUtils.resampleTrace(borderAsTrace, samples); + locs.remove(locs.size()-1); // last is a duplicate + } + + } + + private static class GriddedRegionTestCase extends LocListTestCase { + + public GriddedRegionTestCase(Region cell, int numSamples) { + super(getSampled(cell, numSamples).getNodeList(), numSamples+"x sampled"); + } + + } + + private static class ExteriorNodesTestCase extends LocListTestCase { + + public ExteriorNodesTestCase(Region cell, int numSamples) { + super(null, numSamples+"x exterior"); + GriddedRegion sampled = getSampled(cell, numSamples); + int numLats = sampled.getNumLatNodes(); + int[] minLonIndexes = new int[numLats]; + int[] maxLonIndexes = new int[numLats]; + for (int i=0; i= 0 && lonIndex >= 0); + // overall min + minLatIndex = Integer.min(latIndex, minLatIndex); + maxLatIndex = Integer.max(latIndex, maxLatIndex); + // min/max lon for this lat + minLonIndexes[latIndex] = Integer.min(minLonIndexes[latIndex], lonIndex); + maxLonIndexes[latIndex] = Integer.max(maxLonIndexes[latIndex], lonIndex); + } + System.out.println("LatIndex range: ["+minLatIndex+", "+maxLatIndex+"] for numLat="+numLats); + this.locs = new LocationList(); + for (Location loc : nodes) { + int latIndex = sampled.getLatIndex(loc); + int lonIndex = sampled.getLonIndex(loc); + if (latIndex == minLatIndex || latIndex == maxLatIndex + || lonIndex == minLonIndexes[latIndex] || lonIndex == maxLonIndexes[latIndex]) + locs.add(loc); + } + System.out.println("Kept "+locs.size()+" exterior nodes (of "+sampled.getNodeCount()+" nodes) for "+numSamples+"x sampled"); + } + + } + + private static GriddedRegion getSampled(Region cell, int numSamples) { + double discrX = cell.getMaxLon() - cell.getMinLon(); + double discrY = cell.getMaxLat() - cell.getMinLat(); + return new GriddedRegion(cell, discrY/(double)numSamples, discrX/(double)numSamples, + new Location(discrY/(2d*numSamples), discrX/(2d*numSamples))); + } + + private static interface TestCase extends Named { + + public MinMaxAveTracker run(Location centerLoc, Location[] testLocs); + } + +} diff --git a/src/main/java/scratch/kevin/pointSources/SuperSamplingResolutionTests.java b/src/main/java/scratch/kevin/pointSources/SuperSamplingResolutionTests.java new file mode 100644 index 00000000..caf19154 --- /dev/null +++ b/src/main/java/scratch/kevin/pointSources/SuperSamplingResolutionTests.java @@ -0,0 +1,459 @@ +package scratch.kevin.pointSources; + +import java.io.File; +import java.io.IOException; +import java.text.DecimalFormat; +import java.util.ArrayDeque; +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; +import java.util.Random; +import java.util.concurrent.Callable; +import java.util.concurrent.ExecutionException; +import java.util.concurrent.ExecutorService; +import java.util.concurrent.Executors; +import java.util.concurrent.Future; +import java.util.concurrent.TimeUnit; + +import org.apache.commons.lang3.exception.ExceptionUtils; +import org.opensha.commons.data.CSVFile; +import org.opensha.commons.data.Site; +import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc; +import org.opensha.commons.data.function.DiscretizedFunc; +import org.opensha.commons.data.function.LightFixedXFunc; +import org.opensha.commons.geo.GriddedRegion; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.LocationList; +import org.opensha.commons.geo.Region; +import org.opensha.commons.util.ComparablePairing; +import org.opensha.commons.util.DataUtils.MinMaxAveTracker; +import org.opensha.sha.calc.HazardCurveCalculator; +import org.opensha.sha.calc.params.filters.SourceFilterManager; +import org.opensha.sha.calc.params.filters.SourceFilters; +import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; +import org.opensha.sha.earthquake.faultSysSolution.modules.GridSourceList; +import org.opensha.sha.earthquake.faultSysSolution.util.FaultSysTools; +import org.opensha.sha.earthquake.faultSysSolution.util.SolHazardMapCalc; +import org.opensha.sha.earthquake.faultSysSolution.util.SolHazardMapCalc.ReturnPeriods; +import org.opensha.sha.earthquake.param.IncludeBackgroundOption; +import org.opensha.sha.earthquake.param.IncludeBackgroundParam; +import org.opensha.sha.earthquake.rupForecastImpl.nshm23.util.NSHM23_RegionLoader; +import org.opensha.sha.earthquake.util.GridCellSupersamplingSettings; +import org.opensha.sha.earthquake.util.GriddedSeismicitySettings; +import org.opensha.sha.gui.infoTools.IMT_Info; +import org.opensha.sha.imr.AttenRelRef; +import org.opensha.sha.imr.ScalarIMR; +import org.opensha.sha.imr.param.IntensityMeasureParams.PGA_Param; +import org.opensha.sha.imr.param.IntensityMeasureParams.SA_Param; + +import com.google.common.base.Preconditions; +import com.google.common.base.Stopwatch; + +import scratch.UCERF3.erf.FaultSystemSolutionERF; + +public class SuperSamplingResolutionTests { + + public static void main(String[] args) throws IOException { + Region testReg = NSHM23_RegionLoader.loadFullConterminousWUS(); +// double testSpacing = 0.33d; + double testSpacing = 3d; +// double testSpacing = 1; + GriddedRegion testGridded = new GriddedRegion(testReg, testSpacing, GriddedRegion.ANCHOR_0_0); + System.out.println("Have "+testGridded.getNodeCount()+" test sites"); + + LocationList siteLocs = new LocationList(testGridded.getNodeCount()); + // randomly shift so that they're not perfectly co-located + Random rand = new Random(testGridded.getNodeCount()); + for (Location loc : testGridded.getNodeList()) + siteLocs.add(new Location(loc.lat + testSpacing*(rand.nextDouble() - 0.5d), + loc.lon + testSpacing*(rand.nextDouble() - 0.5d))); + + FaultSystemSolution sol = FaultSystemSolution.load( + new File("/home/kevin/OpenSHA/nshm23/batch_inversions/" + + "2024_02_02-nshm23_branches-WUS_FM_v3/results_WUS_FM_v3_branch_averaged_gridded.zip")); + + File outputDir = new File("/tmp"); + String csvPrefix = "super_sampling_results_"+siteLocs.size()+"sites"; + + ReturnPeriods[] testRPs = ReturnPeriods.values(); + +// AttenRelRef gmmRef = AttenRelRef.ASK_2014; + AttenRelRef gmmRef = AttenRelRef.WRAPPED_ASK_2014; + + double period = 0d; + DiscretizedFunc xVals = new IMT_Info().getDefaultHazardCurve(PGA_Param.NAME); + +// double period = 1d; +// DiscretizedFunc xVals = new IMT_Info().getDefaultHazardCurve(SA_Param.NAME); + + boolean cacheGridSources = true; + + int threads = FaultSysTools.defaultNumThreads(); + + SourceFilterManager sourceFilters = new SourceFilterManager(SourceFilters.TRT_DIST_CUTOFFS); + + ArrayDeque gmmDeque = new ArrayDeque<>(threads); + ArrayDeque calcDeque = new ArrayDeque<>(threads); + + ExecutorService exec = Executors.newFixedThreadPool(threads); + + FaultSystemSolutionERF erf = new FaultSystemSolutionERF(sol); + erf.setParameter(IncludeBackgroundParam.NAME, IncludeBackgroundOption.EXCLUDE); + erf.getTimeSpan().setDuration(1d); + erf.setCacheGridSources(cacheGridSources); + erf.updateForecast(); + + System.out.println("Will calculate "+siteLocs.size()+" curves per test"); + + System.out.println("Calculating fualt-only curves"); + DiscretizedFunc[] faultCurves = calcCurves(siteLocs, erf, period, gmmRef, gmmDeque, calcDeque, xVals, exec, sourceFilters, null); + System.out.println(); + + erf.setParameter(IncludeBackgroundParam.NAME, IncludeBackgroundOption.ONLY); + erf.updateForecast(); + + List samplingParams = new ArrayList<>(); + + int refIndex = 0; // whatever is first is the reference for comparing hazard; put your best foot forward + + samplingParams.add(new double[] {0.5d, 300d, 0d, 0d}); // double resolution to 300 km + samplingParams.add(new double[] {1d, 300d, 0d, 0d}); // full to 300 km + + int benchmarkIndex = samplingParams.size(); // no supersampling is time benchmark + samplingParams.add(null); // no supersampling + + double[] fixedDists = { 20, 30, 40, 60 }; + double[] borderDists = { 0d, 40d, 60d, 80d, 100d }; + double[] cornerDists = { 0d, 60d, 100d, 200d, 300d }; + for (double fixedDist : fixedDists) { + for (double borderDist : borderDists) { + if (borderDist > 0d && borderDist <= fixedDist) + continue; + for (double cornerDist : cornerDists) { + if (cornerDist > 0d && (cornerDist <= fixedDist || cornerDist <= borderDist)) + continue; + samplingParams.add(new double[] {1d, fixedDist, borderDist, cornerDist}); + } + } + } +// samplingParams.add(new double[] {1d, 20d, 0d, 0d}); // full to 20 km +// samplingParams.add(new double[] {1d, 30d, 0d, 0d}); // full to 30 km +// samplingParams.add(new double[] {1d, 40d, 0d, 0d}); // full to 40 km +// samplingParams.add(new double[] {1d, 40d, 0d, 120d}); // variable +// samplingParams.add(new double[] {1d, 30d, 60d, 120d}); // variable +// samplingParams.add(new double[] {0.5d, 30d, 60d, 120d}); // variable +// samplingParams.add(new double[] {1d, 30d, 100d, 200d}); // variable +// samplingParams.add(new double[] {1d, 40d, 100d, 200d}); // variable +// samplingParams.add(new double[] {1d, 20d, 40d, 120d}); // variable +// samplingParams.add(new double[] {1d, 20d, 40d, 0d}); // variable +// samplingParams.add(new double[] {1d, 20d, 60d, 0d}); // variable +// samplingParams.add(new double[] {1d, 30d, 60d, 0d}); // variable +// samplingParams.add(new double[] {1d, 40d, 80d, 0d}); // variable + + List samplingCurves = new ArrayList<>(samplingParams.size()); + List times = new ArrayList<>(samplingParams.size()); + List indexes = new ArrayList<>(samplingCurves.size()); + + GriddedSeismicitySettings gridSettings = GriddedSeismicitySettings.DEFAULT; + + DiscretizedFunc[] refCurves = null; + + for (int s=0; s> csvs = new ArrayList<>(testRPs.length+1); + for (int r=0; r csv = new CSVFile<>(true); + + csv.addLine("Full Sample Distance (km)", "Perimeter Sample Distance (km)", "Corner Distance (km)", + "Time (s)", "Time Factor", "Minimum", "Maximum", "Average", "MAD", "SRSS", "|Min|+|Max|+MAD"); + + // placeholder + List empty = new ArrayList<>(csv.getNumCols()); + for (int i=0; i csv = csvs.get(r); + String rpPrefix = r < testRPs.length ? testRPs[r].name() : "COMBINED"; + File outputFile = new File(outputDir, csvPrefix+"_"+rpPrefix+".csv"); + csv.writeToFile(outputFile); + } + } + + private static final DecimalFormat twoDF = new DecimalFormat("0.00"); + + private static String pFormat(double percent) { + return pFormat(percent, false); + } + + private static String pFormat(double percent, boolean plus) { + String ret = twoDF.format(percent)+"%"; + if (plus && !ret.startsWith("-")) + ret = "+"+ret; + return ret; + } + + private static double[] calcMap(DiscretizedFunc[] curves, ReturnPeriods rp) { + double[] ret = new double[curves.length]; + + double level = rp.oneYearProb; + + for (int i=0; i curves[i].getMaxY()) + val = 0d; + else if (level < curves[i].getMinY()) + // saturated + val = curves[i].getMaxX(); + else + val = curves[i].getFirstInterpolatedX_inLogXLogYDomain(level); + ret[i] = val; + } + + return ret; + } + + private static final DecimalFormat oDF = new DecimalFormat("0.##"); + + private static String samplingName(double[] params) { + if (params == null) + return "No supersampling"; + if (params[2] == 0d && params[3] == 0d) + return "Full up to "+oDF.format(params[1])+" km ("+oDF.format(params[0])+" km spacing)"; + return "Full up to "+oDF.format(params[1])+"-"+oDF.format(params[2])+"-"+oDF.format(params[3])+" km ("+oDF.format(params[0])+" km spacing)"; + } + + private static DiscretizedFunc[] calcCurves(LocationList siteLocs, FaultSystemSolutionERF erf, double period, + AttenRelRef gmmRef, ArrayDeque gmmDeque, ArrayDeque calcDeque, + DiscretizedFunc xVals, ExecutorService exec, SourceFilterManager sourceFilters, + DiscretizedFunc[] combineWith) { + DiscretizedFunc logXVals = new ArbitrarilyDiscretizedFunc(); + for (int i=0; i> futures = new ArrayList<>(siteLocs.size()); + +// System.out.println("Calculating "+siteLocs.size()+" hazard curves"); + + int printMod; + if (siteLocs.size() > 500) + printMod = 100; + else if (siteLocs.size() > 250) + printMod = 50; + else + printMod = 20; + + for (Location siteLoc : siteLocs) { + futures.add(exec.submit(new Callable() { + + @Override + public DiscretizedFunc call() throws Exception { + HazardCurveCalculator calc = null; + ScalarIMR gmm = null; + synchronized (SuperSamplingResolutionTests.class) { + if (!calcDeque.isEmpty()) + calc = calcDeque.pop(); + if (!gmmDeque.isEmpty()) + gmm = gmmDeque.pop(); + } + if (calc == null) { + calc = new HazardCurveCalculator(sourceFilters); + calc.setTrackProgress(false); + } + if (gmm == null) + gmm = gmmRef.get(); + else + gmm.setParamDefaults(); + + if (period == 0d) { + gmm.setIntensityMeasure(PGA_Param.NAME); + } else { + Preconditions.checkState(period > 0d); + gmm.setIntensityMeasure(SA_Param.NAME); + SA_Param.setPeriodInSA_Param(gmm.getIntensityMeasure(), period); + } + + DiscretizedFunc logCurve = new LightFixedXFunc(logXVals); + + Site site = new Site(siteLoc); + site.addParameterList(gmm.getSiteParams()); + + calc.getHazardCurve(logCurve, site, gmm, erf); + + synchronized (SuperSamplingResolutionTests.class) { + calcDeque.push(calc); + gmmDeque.push(gmm); + } + + DiscretizedFunc ret = new ArbitrarilyDiscretizedFunc(); + for (int i=0; i gmms) { diff --git a/src/main/java/scratch/kevin/prvi25/InterfacePointMagCutTests.java b/src/main/java/scratch/kevin/prvi25/InterfacePointMagCutTests.java new file mode 100644 index 00000000..44aaf7f2 --- /dev/null +++ b/src/main/java/scratch/kevin/prvi25/InterfacePointMagCutTests.java @@ -0,0 +1,196 @@ +package scratch.kevin.prvi25; + +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; +import java.util.List; + +import org.opensha.commons.data.Site; +import org.opensha.commons.data.function.EvenlyDiscretizedFunc; +import org.opensha.commons.geo.Location; +import org.opensha.commons.util.DataUtils.MinMaxAveTracker; +import org.opensha.sha.earthquake.ProbEqkRupture; +import org.opensha.sha.earthquake.ProbEqkSource; +import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; +import org.opensha.sha.earthquake.faultSysSolution.modules.GridSourceProvider; +import org.opensha.sha.earthquake.faultSysSolution.modules.GridSourceList.GriddedRupture; +import org.opensha.sha.earthquake.faultSysSolution.util.FaultSysTools; +import org.opensha.sha.earthquake.util.GriddedSeismicitySettings; +import org.opensha.sha.faultSurface.RuptureSurface; +import org.opensha.sha.imr.AttenRelRef; +import org.opensha.sha.imr.ScalarIMR; +import org.opensha.sha.imr.attenRelImpl.nshmp.NSHMP_GMM_Wrapper; +import org.opensha.sha.imr.param.IntensityMeasureParams.PGA_Param; +import org.opensha.sha.util.TectonicRegionType; + +import com.google.common.base.Preconditions; + +import gov.usgs.earthquake.nshmp.gmm.Gmm; + +public class InterfacePointMagCutTests { + + public static void main(String[] args) throws IOException { + FaultSystemSolution sol = FaultSystemSolution.load(new File("/home/kevin/OpenSHA/nshm23/batch_inversions/" + + "2024_11_19-prvi25_subduction_branches/results_PRVI_INTERFACE_ONLY_branch_averaged_gridded.zip")); + + GridSourceProvider gridProv = sol.getGridSourceProvider(); + + GriddedSeismicitySettings fullFiniteSettings = GriddedSeismicitySettings.DEFAULT.forPointSourceMagCutoff(5d); + GriddedSeismicitySettings magCut6FiniteSettings = GriddedSeismicitySettings.DEFAULT.forPointSourceMagCutoff(6d); + + System.out.println("Full finite settings: "+fullFiniteSettings); + System.out.println("Mag cut settings: "+magCut6FiniteSettings); + +// double maxMag = 8.55; +// double maxMag = 7d; + double maxMag = 6d; + EvenlyDiscretizedFunc refMFD = FaultSysTools.initEmptyMFD(5.01, maxMag); + Quantity[] quants = Quantity.values(); + MinMaxAveTracker[][] fullTracks = new MinMaxAveTracker[quants.length][refMFD.size()]; + MinMaxAveTracker[][] cutTracks = new MinMaxAveTracker[quants.length][refMFD.size()]; + MinMaxAveTracker[][] diffTracks = new MinMaxAveTracker[quants.length][refMFD.size()]; + MinMaxAveTracker[][] pDiffTracks = new MinMaxAveTracker[quants.length][refMFD.size()]; + +// ScalarIMR gmm = AttenRelRef.USGS_PRVI_INTERFACE.get(); + +// Gmm gmmRef = Gmm.AG_20_GLOBAL_INTERFACE; +// Gmm gmmRef = Gmm.KBCG_20_GLOBAL_INTERFACE; + Gmm gmmRef = Gmm.PSBAH_20_GLOBAL_INTERFACE; + ScalarIMR gmm = new NSHMP_GMM_Wrapper(gmmRef, gmmRef.toString(), gmmRef.name(), false, null); + gmm.setIntensityMeasure(PGA_Param.NAME); + + for (int q=0; q locs; +// locs = new ArrayList<>(gridProv.getNumLocations()); +// for (int l=0; l> fullRupsBinned = new ArrayList<>(refMFD.size()); + List> cutRupsBinned = new ArrayList<>(refMFD.size()); + for (int i=0; i()); + cutRupsBinned.add(new ArrayList<>()); + } + double rateSumFull = 0d; + for (ProbEqkRupture rup : sourceFull) { + if (rup.getMag() > maxMag) + continue; + fullRupsBinned.get(refMFD.getClosestXIndex(rup.getMag())).add(rup); + rateSumFull += rup.getMeanAnnualRate(1d); + } + double rateSumCut = 0d; + for (ProbEqkRupture rup : sourceCut) { + if (rup.getMag() > maxMag) + continue; + cutRupsBinned.get(refMFD.getClosestXIndex(rup.getMag())).add(rup); + rateSumCut += rup.getMeanAnnualRate(1d); + } + Preconditions.checkState((float)rateSumFull == (float)rateSumCut, "Rate mismatch: %s != %s", (float)rateSumFull, (float)rateSumCut); + + for (int i=0; i rupsFull = fullRupsBinned.get(i); + List rupsCut = cutRupsBinned.get(i); + if (rupsFull.isEmpty() || rupsCut.isEmpty()) { + Preconditions.checkArgument(rupsFull.isEmpty() && rupsCut.isEmpty()); + continue; + } + + for (int q=0; q fullVals = new ArrayList<>(); + for (ProbEqkRupture rup : rupsFull) { + double val = quants[q].calc(rup, testLoc, gmm); + fullVals.add(val); + fullTracks[q][i].addValue(val); + } + List cutVals = new ArrayList<>(); + for (ProbEqkRupture rup : rupsCut) { + double val = quants[q].calc(rup, testLoc, gmm); + cutVals.add(val); + cutTracks[q][i].addValue(val); + } + + for (double fullVal : fullVals) { + for (double cutVal : cutVals) { + double diff = cutVal - fullVal; + double pDiff = 100d*diff/fullVal; + + diffTracks[q][i].addValue(diff); + pDiffTracks[q][i].addValue(pDiff); + } + } + } + } + } + + for (int m=0; m crustalSects = crustalSol.getRupSet().getFaultSectionDataList(); + List subSmallSects = subSmallSol.getRupSet().getFaultSectionDataList(); + List subLargeSects = subLargeSol.getRupSet().getFaultSectionDataList(); + System.out.println("SubSmall has "+subSmallSects.size()+" sects"); + System.out.println("SubLarge has "+subLargeSects.size()+" sects"); + + boolean allMatch = subSmallSects.size() == subLargeSects.size(); + for (int i=0; i= subSmallSects.size()) { + name1 = "MISSING"; + } else { + FaultSection sect = subSmallSects.get(i); + name1 = sect.getSectionName()+" ("+sect.getParentSectionId()+")"; + } + String name2; + if (i >= subLargeSects.size()) { + name2 = "MISSING"; + } else { + FaultSection sect = subLargeSects.get(i); + name2 = sect.getSectionName()+" ("+sect.getParentSectionId()+")"; + } + boolean match = name1.equals(name2); + System.out.println(i+". "+name1+"\t"+name2+(match ? "" : "\tMISMATCH")); + allMatch &= match; + } + Preconditions.checkState(allMatch, "Sub FM sects don't match"); + + double[] minMags = { 0d, 6d, 6.5, 7d, 7.5d, 8d }; + double[] gridMinMags = { 5d, 6d, 6.5, 7d, 7.5d, 8d }; + + List combSects = new ArrayList<>(); + List> combSectRates = new ArrayList<>(); + List> combSectSortables = new ArrayList<>(); + for (int m=0; m()); + combSectSortables.add(new ArrayList<>()); + } + + // first crustal + for (FaultSection sect : crustalSects) { + for (int m=0; m smallGriddedRates = null; + List largeGriddedRates = null; + if (includeInterfaceGridded) { + for (boolean large : new boolean[] {false,true}) { + List griddedRates = new ArrayList<>(subSmallSects.size()); + for (int s=0; s= minMags[m]) + rates[m] += rup.associatedSectionFracts[i] * rup.rate; + } + } + } + if (large) + largeGriddedRates = griddedRates; + else + smallGriddedRates = griddedRates; + } + } + + // now subduction combined + for (int s=0; s 0d) + rates[m] += gridRates[m]; + } + + if (large) + largeRates = rates; + else + smallRates = rates; + } + + if (smallSect.getParentSectionId() == muertosParent || smallSect.getParentSectionId() == hispaniolaParent) { + // simple, just weight-average them and represent once + for (int m=0; m rates : combSectRates) + for (int r=0; r rates : combSectRates) + for (int r=0; r 0d) { + label += "M>"+oDF.format(minMags[m]); + prefix += "m"+(float)minMags[m]; + noLegendTitle = "Fault, M>"+oDF.format(minMags[m]); + } else { + label += "Supra-Seismogenic"; + prefix += "supra_seis"; + noLegendTitle = "Fault, Supra-Seismogenic"; + } + if (includeInterfaceGridded) + prefix += "_incl_interface_gridded"; + label += " Fault Participation Rate (/yr)"; + + mapMaker.plotSectScalars(combSectRates.get(m), combSectSortables.get(m), cpt, label); + mapMaker.plot(outputDir, prefix, " "); + + // now again without a CPT legend + mapMaker.plotSectScalars(combSectRates.get(m), combSectSortables.get(m), cpt, null); + mapMaker.plot(outputDir, prefix+"_no_cpt", noLegendTitle); + } + + // now gridded + TectonicRegionType[] trts = { + null, + TectonicRegionType.ACTIVE_SHALLOW, + TectonicRegionType.SUBDUCTION_INTERFACE, + TectonicRegionType.SUBDUCTION_SLAB + }; + + mapMaker.clearSectScalars(); + mapMaker.setSectTraceChar(new PlotCurveCharacterstics(PlotLineType.SOLID, 1f, Color.DARK_GRAY)); + GridSourceList gridProv = combSol.requireModule(GridSourceList.class); + GriddedRegion gridReg = gridProv.getGriddedRegion(); + cpt.setNanColor(Color.WHITE); + double minGridMag = Doubles.min(gridMinMags); + int numRandPartic = 1000; + BitSet nodeBits = null; + Table particsCache = HashBasedTable.create(); + Random rand = new Random(123456789l); + for (boolean psuedoPartic : new boolean[] {false, true}) { + if (psuedoPartic) { + // precache + ExecutorService exec = Executors.newFixedThreadPool(FaultSysTools.defaultNumThreads()); + Table> particsCacheFutures = HashBasedTable.create(); + System.out.println("Precaching gridded participation maps"); + for (int l=0; l future = particsCacheFutures.get(gridLoc.lat, rup.properties); + if (future == null) { + long seed = rand.nextLong(); + future = exec.submit(new Callable() { + + @Override + public CachedGriddedParticResult call() throws Exception { + return new CachedGriddedParticResult(rup.properties, gridReg, gridLoc.lat, numRandPartic, new Random(seed)); + } + }); + particsCacheFutures.put(gridLoc.lat, rup.properties, future); + } + } + } + System.out.println("Waiting on "+particsCacheFutures.size()+" futures"); + for (Cell> cell : particsCacheFutures.cellSet()) { + try { + CachedGriddedParticResult result = cell.getValue().get(); + particsCache.put(cell.getRowKey(), cell.getColumnKey(), result); + } catch (InterruptedException | ExecutionException e) { + throw ExceptionUtils.asRuntimeException(e); + } + } + + exec.shutdown(); + } + for (TectonicRegionType trt : trts) { + GriddedGeoDataSet[] xyzs = new GriddedGeoDataSet[gridMinMags.length]; + for (int m=0; m 0d); + + double regWidth = props.length; + if (props.dip != 90d) { + // sin(dip) = horz / vert + // horz = vert * sin(dip) + double ddwProj = Math.sin(Math.toRadians(props.dip)) * (props.lowerDepth - props.upperDepth); + regWidth += ddwProj; // rotations get funky, buffer by some extra + } + regWidth += 50; + + Location refLoc = new Location(lat, 0.5*(gridReg.getMinLon() + gridReg.getMaxLon())); + Location origRefLoc = refLoc; + refLoc = gridReg.getLocation(gridReg.indexForLocation(refLoc)); // snap to grid + Preconditions.checkNotNull(refLoc, "Couldn't snap middle loc to grid? orig=[%s]", origRefLoc); + + Preconditions.checkState(gridReg.getLatSpacing() == gridReg.getLonSpacing()); + double spacing = gridReg.getSpacing(); + + Location bottom = LocationUtils.location(refLoc, Math.PI, 0.5*regWidth); + Location bottomLeft = LocationUtils.location(bottom, 1.5*Math.PI, 0.5*regWidth); + Location top = LocationUtils.location(refLoc, 0d, 0.5*regWidth); + Location topRight = LocationUtils.location(top, 0.5*Math.PI, 0.5*regWidth); + Region surroundingRegion = new Region(bottomLeft, topRight); + + surroundingGrid = new GriddedRegion(surroundingRegion, spacing, gridReg.getLocation(0)); + xyz = new EvenlyDiscrXYZ_DataSet( + surroundingGrid.getNumLonNodes(), surroundingGrid.getNumLatNodes(), + surroundingGrid.getMinGridLon()-refLoc.lon, surroundingGrid.getMinGridLat()-refLoc.lat, spacing); + Preconditions.checkState(xyz.size() == surroundingGrid.getNodeCount()); + + GriddedRupture fakeRup = new GriddedRupture(0, refLoc, props, 1d); + + PointSurfaceBuilder builder = GridSourceList.surfBuilderForRup(fakeRup); + + builder.random(rand); + + EvenlyGriddedSurface[] surfs; + if (props.dip == 90d) + surfs = builder.buildRandLineSurfaces(numSurfs); + else + surfs = builder.buildRandGriddedSurfaces(numSurfs); + + BitSet mappings = new BitSet(xyz.size()); + double rateEach = 1d/surfs.length; + for (int s=0; s 0) + mappings.clear(); + for (Location loc : surfs[s]) { + int mappedIndex = surroundingGrid.indexForLocation(loc); + Preconditions.checkState(mappedIndex >= 0, + "Buffer wasn't big enough? Surf location (%s) didn't map to index. Len=%s, regWidth=%s" + + "\n\trefLoc=[%s], botLeft=[%s], topRight=[%s]" + + "\n\tgridMin=[%s, %s], gridMax=[%s, %s]" + + "\n\tprops=%s", + loc, props.length, regWidth, refLoc, bottomLeft, topRight, + surroundingGrid.getMinGridLat(), surroundingGrid.getMinGridLon(), + surroundingGrid.getMaxGridLat(), surroundingGrid.getMaxGridLon(), + props); + mappings.set(mappedIndex); + } + + for (int i=0; i elements, File outputDir if (fssForComparison != null) { IncrementalMagFreqDist compSupraSeis = fssForComparison.calcTotalNucleationMFD( - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF, 8.55, 0.1); + 5.05, 8.55, 0.1); compSupraSeis.setName("UCERF3 Supra-Seis"); SummedMagFreqDist compSubSeis = null; SummedMagFreqDist compTotal = null; diff --git a/src/main/java/scratch/kevin/spatialVar/U3CatalogVarCalcDemo.java b/src/main/java/scratch/kevin/spatialVar/U3CatalogVarCalcDemo.java index 77b4b87a..8fef30a5 100644 --- a/src/main/java/scratch/kevin/spatialVar/U3CatalogVarCalcDemo.java +++ b/src/main/java/scratch/kevin/spatialVar/U3CatalogVarCalcDemo.java @@ -14,10 +14,14 @@ import org.opensha.sha.earthquake.ProbEqkSource; import org.opensha.sha.earthquake.faultSysSolution.FaultSystemRupSet; import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; +import org.opensha.sha.earthquake.faultSysSolution.erf.BaseFaultSystemSolutionERF; import org.opensha.sha.earthquake.faultSysSolution.modules.GridSourceProvider; import org.opensha.sha.earthquake.param.BackgroundRupType; +import org.opensha.sha.earthquake.util.GriddedSeismicitySettings; import org.opensha.sha.faultSurface.PointSurface; import org.opensha.sha.faultSurface.RuptureSurface; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrection; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrections; import org.opensha.sha.imr.attenRelImpl.ngaw2.NGAW2_Wrappers.ASK_2014_Wrapper; import org.opensha.sha.imr.attenRelImpl.ngaw2.ScalarGroundMotion; import org.opensha.sha.imr.param.IntensityMeasureParams.SA_Param; @@ -82,6 +86,8 @@ public static void main(String[] args) throws IOException { GridSourceProvider gridProv = fss.getGridSourceProvider(); GriddedRegion gridReg = gridProv.getGriddedRegion(); + GriddedSeismicitySettings gridSettings = GriddedSeismicitySettings.DEFAULT; + System.out.println("Calculating ground motions for catalog"); for (int i=0; i= 0) { // match - ProbEqkSource gridSource = gridProv.getSource(gridIndex, 1d, null, BackgroundRupType.POINT); + ProbEqkSource gridSource = gridProv.getSource(gridIndex, 1d, null, gridSettings); // there can be multiple ruptures with that magnitude, possibly with different rakes. randomly choose one double magTol = 0.05; diff --git a/src/main/java/scratch/kevin/ucerf3/PureScratch.java b/src/main/java/scratch/kevin/ucerf3/PureScratch.java index c8be0db2..d1ec73cb 100644 --- a/src/main/java/scratch/kevin/ucerf3/PureScratch.java +++ b/src/main/java/scratch/kevin/ucerf3/PureScratch.java @@ -156,8 +156,6 @@ import org.opensha.sha.earthquake.param.ProbabilityModelOptions; import org.opensha.sha.earthquake.param.ProbabilityModelParam; import org.opensha.sha.earthquake.param.UseRupMFDsParam; -import org.opensha.sha.earthquake.rupForecastImpl.PointSource13b.PointSurface13b; -import org.opensha.sha.earthquake.rupForecastImpl.PointSourceNshm.PointSurfaceNshm; import org.opensha.sha.earthquake.rupForecastImpl.WGCEP_UCERF_2_Final.FaultSegmentData; import org.opensha.sha.earthquake.rupForecastImpl.WGCEP_UCERF_2_Final.MeanUCERF2.MeanUCERF2; import org.opensha.sha.earthquake.rupForecastImpl.WGCEP_UCERF_2_Final.data.A_FaultsFetcher; @@ -1492,7 +1490,6 @@ private static void test290() throws IOException { } private static void test291() throws IOException { - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = 2.55; FaultSystemSolution sol = FaultSystemSolution.load(new File("/data/kevin/git/ucerf3-etas-launcher/inputs/" // FaultSystemSolution sol = FaultSystemSolution.load(new File("/home/kevin/OpenSHA/UCERF3/rup_sets/orig/" + "2013_05_10-ucerf3p3-production-10runs_COMPOUND_SOL_FM3_1_SpatSeisU3_MEAN_BRANCH_AVG_SOL.zip")); @@ -1727,77 +1724,77 @@ private static void test296() throws IOException { System.out.println("All nth tests passed"); } - private static void test297() throws IOException { - // nth rup from catalog and ETAS_Launcher.buildERF - File etasDir = new File("/home/kevin/OpenSHA/UCERF3/etas/simulations/" - + "2024_05_24-Start2012_500yr_kCOV1p5_Spontaneous_HistCatalog/"); - File configFile = new File(etasDir, "config.json"); - File binFile = new File(etasDir, "results_m5_preserve_chain.bin"); - - ETAS_Config config = ETAS_Config.readJSON(configFile); - - ETAS_Launcher launcher = new ETAS_Launcher(config, false); - FaultSystemSolutionERF erf = (FaultSystemSolutionERF) launcher.checkOutERF(); - MFDGridSourceProvider gridProv = erf.getSolution().requireModule(MFDGridSourceProvider.class); - - int randSrcIndex = 252798; - int gridRegionIndex = randSrcIndex - erf.getNumFaultSystemSources(); - int isSubSeismo = 1; - int filteredR = 2; - Location translatedParLoc = new Location(35, -118); - ProbEqkSource src = erf.getSource(randSrcIndex); - ProbEqkSource filteredSrc; - if(isSubSeismo == 1) - filteredSrc = gridProv.getSourceSubSeisOnFault(gridRegionIndex, - erf.getTimeSpan().getDuration(), null, BackgroundRupType.POINT); - else - filteredSrc = gridProv.getSourceUnassociated(gridRegionIndex, - erf.getTimeSpan().getDuration(), null, BackgroundRupType.POINT); - - ProbEqkRupture filteredRup = filteredSrc.getRupture(filteredR); - System.out.println("Filtered r="+filteredR+", M"+filteredRup.getMag()+", rake="+filteredRup.getAveRake()); - - int r = -1; - for (int r2=0; r2= 0, "No rupture mapping found sourceID=%s, subSeismo=%s," - + "filteredR=%s, mag=%s, rake=%s, sourceType=%s, surfType=%s", - randSrcIndex, isSubSeismo, filteredR, filteredRup.getMag(), filteredRup.getAveRake(), - filteredSrc.getClass().getName(), filteredRup.getRuptureSurface().getClass().getName()); - ProbEqkRupture rup = src.getRupture(r); - System.out.println("Matched with r="+r+", M"+rup.getMag()+", rake="+rup.getAveRake()); - - } +// private static void test297() throws IOException { +// // nth rup from catalog and ETAS_Launcher.buildERF +// File etasDir = new File("/home/kevin/OpenSHA/UCERF3/etas/simulations/" +// + "2024_05_24-Start2012_500yr_kCOV1p5_Spontaneous_HistCatalog/"); +// File configFile = new File(etasDir, "config.json"); +// File binFile = new File(etasDir, "results_m5_preserve_chain.bin"); +// +// ETAS_Config config = ETAS_Config.readJSON(configFile); +// +// ETAS_Launcher launcher = new ETAS_Launcher(config, false); +// FaultSystemSolutionERF erf = (FaultSystemSolutionERF) launcher.checkOutERF(); +// MFDGridSourceProvider gridProv = erf.getSolution().requireModule(MFDGridSourceProvider.class); +// +// int randSrcIndex = 252798; +// int gridRegionIndex = randSrcIndex - erf.getNumFaultSystemSources(); +// int isSubSeismo = 1; +// int filteredR = 2; +// Location translatedParLoc = new Location(35, -118); +// ProbEqkSource src = erf.getSource(randSrcIndex); +// ProbEqkSource filteredSrc; +// if(isSubSeismo == 1) +// filteredSrc = gridProv.getSourceSubSeisOnFault(gridRegionIndex, +// erf.getTimeSpan().getDuration(), null, BackgroundRupType.POINT); +// else +// filteredSrc = gridProv.getSourceUnassociated(gridRegionIndex, +// erf.getTimeSpan().getDuration(), null, BackgroundRupType.POINT); +// +// ProbEqkRupture filteredRup = filteredSrc.getRupture(filteredR); +// System.out.println("Filtered r="+filteredR+", M"+filteredRup.getMag()+", rake="+filteredRup.getAveRake()); +// +// int r = -1; +// for (int r2=0; r2= 0, "No rupture mapping found sourceID=%s, subSeismo=%s," +// + "filteredR=%s, mag=%s, rake=%s, sourceType=%s, surfType=%s", +// randSrcIndex, isSubSeismo, filteredR, filteredRup.getMag(), filteredRup.getAveRake(), +// filteredSrc.getClass().getName(), filteredRup.getRuptureSurface().getClass().getName()); +// ProbEqkRupture rup = src.getRupture(r); +// System.out.println("Matched with r="+r+", M"+rup.getMag()+", rake="+rup.getAveRake()); +// +// } private static void test298() throws IOException { File inDir = new File("/data/kevin/nshm23/batch_inversions/2024_02_02-nshm23_branches-WUS_FM_v3/node_branch_averaged"); @@ -2211,25 +2208,25 @@ private static void test309() throws IOException { } } - private static void test310() throws IOException { - FaultSystemSolution sol = FaultSystemSolution.load(new File("/tmp/PRVI_SUB_FM_LARGE_with_gridded.zip")); - GridSourceList gridSources = sol.requireModule(GridSourceList.class); - for (TectonicRegionType trt : gridSources.getTectonicRegionTypes()) { - System.out.println(trt); - for (int gridIndex=0; gridIndex etasCatalog; private boolean fullTraces; + private double minMag; public UCERF3_IM_Calculator(CommandLine cmd) throws IOException, DocumentException { if (cmd.hasOption("min-mag")) { - double minMag = Double.parseDouble(cmd.getOptionValue("min-mag")); + minMag = Double.parseDouble(cmd.getOptionValue("min-mag")); System.out.println("Setting UCERF3 minimum magnitude to: "+(float)minMag); - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = minMag; + } else { + minMag = BaseFaultSystemSolutionERF.GRID_SETTINGS_DEFAULT.minimumMagnitude; } File fssFile = new File(cmd.getOptionValue("sol-file")); FaultSystemSolution sol = U3FaultSystemIO.loadSol(fssFile); @@ -139,10 +142,12 @@ private void init(FaultSystemSolution sol, Site site, AttenRelRef gmpe, boolean public void calculate() throws IOException { FaultSystemSolutionERF erf = new FaultSystemSolutionERF(sol); - if (doGridded) + if (doGridded) { erf.setParameter(IncludeBackgroundParam.NAME, IncludeBackgroundOption.INCLUDE); - else + erf.setGriddedSeismicitySettings(erf.getGriddedSeismicitySettings().forMinimumMagnitude(minMag)); + } else { erf.setParameter(IncludeBackgroundParam.NAME, IncludeBackgroundOption.EXCLUDE); + } erf.getTimeSpan().setDuration(ERF_DURATION); erf.updateForecast(); @@ -225,7 +230,7 @@ public void calculate() throws IOException { csv.addLine(processRupture(rup, site, gmpe, etasRup, fssIndex, gridNodeIndex, hasInterStdDev, hasIntraStdDev, hasRrup, hasRJB)); } else if (doGridded) { - if (etasRup.getMag() < AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF) + if (etasRup.getMag() < minMag) continue; int sourceID = gridNodeIndex + numSourcesFSS; Preconditions.checkState(sourceID >= numSourcesFSS && sourceID < erf.getNumSources(), diff --git a/src/main/java/scratch/kevin/ucerf3/eal/UCERF3_BranchAvgLossFetcher.java b/src/main/java/scratch/kevin/ucerf3/eal/UCERF3_BranchAvgLossFetcher.java index 6cbe8f6e..a349aa93 100644 --- a/src/main/java/scratch/kevin/ucerf3/eal/UCERF3_BranchAvgLossFetcher.java +++ b/src/main/java/scratch/kevin/ucerf3/eal/UCERF3_BranchAvgLossFetcher.java @@ -11,6 +11,7 @@ import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc; import org.opensha.commons.data.function.DiscretizedFunc; import org.opensha.commons.geo.GriddedRegion; +import org.opensha.sha.earthquake.faultSysSolution.erf.BaseFaultSystemSolutionERF; import org.opensha.sha.earthquake.faultSysSolution.modules.GridSourceProvider; import org.opensha.sha.imr.AttenRelRef; import org.opensha.sha.magdist.IncrementalMagFreqDist; @@ -260,7 +261,7 @@ public static void main(String[] args) throws IOException, DocumentException { watch.reset(); double totGriddedLosses = 0; for (int i=0; i> loadCatalogs(File resultsBinFile) throws IOException { Preconditions.checkArgument(resultsBinFile.exists(), "catalog file doesn't exist"); - return loadCatalogs(resultsBinFile, AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF-0.05); + return loadCatalogs(resultsBinFile, LOSS_CALC_MIN_MAG-0.05); } public ETAS_CatalogEALCalculator(UCERF3_BranchAvgLossFetcher fetcher, U3FaultSystemSolution meanSol, @@ -148,11 +146,8 @@ public ETAS_CatalogEALCalculator(UCERF3_BranchAvgLossFetcher fetcher, U3FaultSys LastEventData.populateSubSects(meanSol.getRupSet().getFaultSectionDataList(), LastEventData.load()); this.meanSol = meanSol; System.out.println("Loading ERF"); - double origMinMag = AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF; - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = 2.55; erf = ETAS_Launcher.buildERF(meanSol, false, 1d, 2014); erf.updateForecast(); - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = origMinMag; System.out.println("Done loading ERF"); } @@ -519,7 +514,7 @@ static int calcNodeIndex(ETAS_EqkRupture rup, GriddedRegion region) { double calcGridSourceLoss(ETAS_EqkRupture rup, GriddedRegion region, DiscretizedFunc[] griddedMagLossDists, String catName) { double mag = rup.getMag(); - if ((float)mag < (float)AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF) + if ((float)mag < (float)LOSS_CALC_MIN_MAG) // below our min mag return 0; Location loc = rup.getHypocenterLocation(); @@ -565,8 +560,8 @@ static int calcNodeIndex(ETAS_EqkRupture rup, GriddedRegion region) { Preconditions.checkState(sourceID >= 0); ProbEqkSource source = erf.getSource(sourceID); Location sourceLoc = null; - if (source instanceof PointSource13b) - sourceLoc = ((PointSource13b)source).getLocation(); + if (source instanceof PointSource) + sourceLoc = ((PointSource)source).getLocation(); else if (source instanceof Point2Vert_FaultPoisSource) sourceLoc = ((Point2Vert_FaultPoisSource)source).getLocation(); else @@ -1407,7 +1402,7 @@ public void writeLossesToCSV(File csvFile, List lossDists) thro CSVFile csv = new CSVFile(true); - double cutoffMag = AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF; + double cutoffMag = LOSS_CALC_MIN_MAG; csv.addLine("Index", "Total Mean Loss", "# FSS Ruptures", "# M>="+(float)cutoffMag, "Max Mag"); diff --git a/src/main/java/scratch/kevin/ucerf3/etas/ETAS_CatalogGridSourceProvider.java b/src/main/java/scratch/kevin/ucerf3/etas/ETAS_CatalogGridSourceProvider.java index a60383d7..73bdb15c 100644 --- a/src/main/java/scratch/kevin/ucerf3/etas/ETAS_CatalogGridSourceProvider.java +++ b/src/main/java/scratch/kevin/ucerf3/etas/ETAS_CatalogGridSourceProvider.java @@ -14,9 +14,13 @@ import org.opensha.sha.earthquake.AbstractERF; import org.opensha.sha.earthquake.ProbEqkRupture; import org.opensha.sha.earthquake.ProbEqkSource; +import org.opensha.sha.earthquake.faultSysSolution.erf.BaseFaultSystemSolutionERF; import org.opensha.sha.earthquake.faultSysSolution.modules.GridSourceProvider; import org.opensha.sha.earthquake.faultSysSolution.modules.MFDGridSourceProvider; import org.opensha.sha.earthquake.param.BackgroundRupType; +import org.opensha.sha.earthquake.util.GriddedSeismicitySettings; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrection; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrections; import org.opensha.sha.magdist.IncrementalMagFreqDist; import org.opensha.sha.util.TectonicRegionType; @@ -49,6 +53,8 @@ public class ETAS_CatalogGridSourceProvider extends AbstractGridSourceProvider { private Map nodeMFDs; + private GriddedSeismicitySettings gridSeisSettings = GriddedSeismicitySettings.DEFAULT; + /** * If true, then this class is just used to create ruptures with correct focal mechanisms and such, and * ProbEqkRup objects are created with relative rates for each focal mech conditioned on the occurance @@ -244,7 +250,7 @@ public Iterable getConditionalRuptures(ETAS_EqkRupture etasRup) IncrementalMagFreqDist nodeMFD = getMFD_Unassociated(node); Preconditions.checkNotNull(nodeMFD, "Rupture maps to uninitialized node!"); Preconditions.checkState(nodeMFD.getY(mfdIndex) > 0, "Mag uninitialized in MFD node!"); - IncrementalMagFreqDist trimmedMFD = getMFD(node, SOURCE_MIN_MAG_CUTOFF); + IncrementalMagFreqDist trimmedMFD = getMFD(node, gridSeisSettings.minimumMagnitude); // Preconditions.checkState(trimmedMFD.size() > mfdIndex, // "Trimmed MFD cuts it off! WTF?\n\nOrig MFD:\n%s\n\nTrimmedMFD\n%s", nodeMFD, trimmedMFD); if (trimmedMFD.size() <= mfdIndex) { @@ -259,7 +265,7 @@ public Iterable getConditionalRuptures(ETAS_EqkRupture etasRup) // this generates a new source instance, but the ruptures in that source // reuse properties. so we can only iterate over it, thus the custom iterable - ProbEqkSource source = getSource(node, 1d, null, BackgroundRupType.POINT); + ProbEqkSource source = getSource(node, 1d, null, gridSeisSettings); SubsetIterable iterable = new SubsetIterable(source); List rupMags = Lists.newArrayList(); @@ -336,7 +342,7 @@ public void updateForecast() {} @Override public ProbEqkSource getSource(int idx) { - return gridProv.getSource(sourceIndexes.get(idx), 1d, null, BackgroundRupType.POINT); + return gridProv.getSource(sourceIndexes.get(idx), 1d, null, gridSeisSettings); } @Override diff --git a/src/main/java/scratch/kevin/ucerf3/etas/MPJ_ETAS_CharFactorCalc.java b/src/main/java/scratch/kevin/ucerf3/etas/MPJ_ETAS_CharFactorCalc.java index 505b9b5f..3a937947 100644 --- a/src/main/java/scratch/kevin/ucerf3/etas/MPJ_ETAS_CharFactorCalc.java +++ b/src/main/java/scratch/kevin/ucerf3/etas/MPJ_ETAS_CharFactorCalc.java @@ -67,7 +67,6 @@ public class MPJ_ETAS_CharFactorCalc extends MPJTaskCalculator { public MPJ_ETAS_CharFactorCalc(CommandLine cmd, File compoundFile, File outputDir) throws ZipException, IOException { super(cmd); - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = 2.55; weightProv = new U3APrioriBranchWeightProvider(); cfss = U3CompoundFaultSystemSolution.fromZipFile(compoundFile); diff --git a/src/main/java/scratch/kevin/ucerf3/etas/PostEventETASPlotter.java b/src/main/java/scratch/kevin/ucerf3/etas/PostEventETASPlotter.java index 8f5ea878..2d8fab1d 100644 --- a/src/main/java/scratch/kevin/ucerf3/etas/PostEventETASPlotter.java +++ b/src/main/java/scratch/kevin/ucerf3/etas/PostEventETASPlotter.java @@ -59,7 +59,6 @@ public static void main(String[] args) throws IOException, DocumentException { ETAS_MultiSimAnalysisTools.plotMagNum(catalogs, plotDir, name, "one_week", null, 0, null, oneWeekOT); if (scenario != null) { - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = 2.55; File fssFile = new File("/home/kevin/workspace/OpenSHA/dev/scratch/UCERF3/data/scratch/InversionSolutions/" + "2013_05_10-ucerf3p3-production-10runs_COMPOUND_SOL_FM3_1_MEAN_BRANCH_AVG_SOL.zip"); FaultSystemSolution fss = U3FaultSystemIO.loadSol(fssFile); diff --git a/src/main/java/scratch/kevin/ucerf3/inversion/CompoundGriddedRecalc.java b/src/main/java/scratch/kevin/ucerf3/inversion/CompoundGriddedRecalc.java deleted file mode 100644 index 860d7b4d..00000000 --- a/src/main/java/scratch/kevin/ucerf3/inversion/CompoundGriddedRecalc.java +++ /dev/null @@ -1,110 +0,0 @@ -package scratch.kevin.ucerf3.inversion; - -import java.io.File; -import java.io.IOException; -import java.util.ArrayDeque; -import java.util.List; -import java.util.NoSuchElementException; -import java.util.zip.ZipException; - -import org.opensha.commons.util.ExceptionUtils; -import org.opensha.sha.earthquake.faultSysSolution.modules.GridSourceProvider; -import org.opensha.sha.earthquake.faultSysSolution.modules.MFDGridSourceProvider; - -import com.google.common.collect.Lists; - -import scratch.UCERF3.U3CompoundFaultSystemSolution; -import scratch.UCERF3.griddedSeismicity.AbstractGridSourceProvider; -import scratch.UCERF3.griddedSeismicity.GridSourceFileReader; -import scratch.UCERF3.inversion.InversionFaultSystemSolution; -import scratch.UCERF3.logicTree.U3LogicTreeBranch; - -public class CompoundGriddedRecalc { - - public static void main(String[] args) throws ZipException, IOException { - File compoundFile = new File("/home/kevin/workspace/OpenSHA/dev/scratch/UCERF3/data/scratch/" - + "InversionSolutions/2013_05_10-ucerf3p3-production-10runs_COMPOUND_SOL.zip"); - U3CompoundFaultSystemSolution cfss = U3CompoundFaultSystemSolution.fromZipFile(compoundFile); - - File writeDir = new File("/tmp/cfss_branches"); - if (!writeDir.exists()) - writeDir.mkdir(); - - ArrayDeque branchQueue = new ArrayDeque(); - branchQueue.addAll(cfss.getBranches()); - - int numThreads = 4; - - List threads = Lists.newArrayList(); - - for (int i=0; i branches) { - if (branches.size() % 50 == 0) - System.out.println(branches.size()+" braches left"); - return branches.pop(); - } - - private static boolean region_written = false; - - private static synchronized boolean isRegionWritten() { - if (!region_written) { - region_written = true; - return false; - } - return true; - } - - private static class CalcThread extends Thread { - private ArrayDeque branches; - private File dir; - private U3CompoundFaultSystemSolution cfss; - - public CalcThread(U3CompoundFaultSystemSolution cfss, File dir, ArrayDeque branches) { - this.branches = branches; - this.dir = dir; - this.cfss = cfss; - } - - @Override - public void run() { - while (true) { - U3LogicTreeBranch branch; - try { - branch = popBranch(branches); - } catch (NoSuchElementException e) { - // done - return; - } - - InversionFaultSystemSolution sol = cfss.getSolution(branch); - MFDGridSourceProvider gridSources = sol.requireModule(MFDGridSourceProvider.class); - - File regXMLFile = null; - if (!isRegionWritten()) - regXMLFile = new File(dir, U3CompoundFaultSystemSolution.getRemappedName("grid_sources_reg.xml", branch)); - File binFile = new File(dir, U3CompoundFaultSystemSolution.getRemappedName("grid_sources.bin", branch)); - try { -// GridSourceFileReader.writeGriddedSeisFile(file, gridSources); - GridSourceFileReader.writeGriddedSeisBinFile(binFile, regXMLFile, gridSources, - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF); - } catch (IOException e) { - ExceptionUtils.throwAsRuntimeException(e); - } - } - } - } - -} diff --git a/src/main/java/scratch/ned/ETAS_ERF/testModels/TestModel1_ERF.java b/src/main/java/scratch/ned/ETAS_ERF/testModels/TestModel1_ERF.java index 46c9550d..a0435237 100644 --- a/src/main/java/scratch/ned/ETAS_ERF/testModels/TestModel1_ERF.java +++ b/src/main/java/scratch/ned/ETAS_ERF/testModels/TestModel1_ERF.java @@ -157,7 +157,8 @@ protected ProbEqkSource getOtherSource(int iSource) { double magCutOff =8.0; // all rups below are treated a point sources boolean isCrossHair=true; return new Point2Vert_FaultPoisSource(griddedRegion.getLocation(regionIndex), mfd, - magLengthRel,timeSpan.getDuration(), magCutOff ,1.0, 0.0,0.0, isCrossHair); + magLengthRel,timeSpan.getDuration(), magCutOff ,1.0, 0.0,0.0, isCrossHair, + /* no dist corrs */ null); } diff --git a/src/main/java/scratch/ned/ETAS_Tests/Tests.java b/src/main/java/scratch/ned/ETAS_Tests/Tests.java index a4e7a10e..c3b6b390 100644 --- a/src/main/java/scratch/ned/ETAS_Tests/Tests.java +++ b/src/main/java/scratch/ned/ETAS_Tests/Tests.java @@ -134,7 +134,7 @@ public void mkProjectPlanETAS_FigureData() { Location[] locs; int numLocs; - NSHMP_GridSourceGenerator nshmpGen = new NSHMP_GridSourceGenerator(); + NSHMP_GridSourceGenerator nshmpGen = new NSHMP_GridSourceGenerator(null); GriddedRegion region = nshmpGen.getGriddedRegion(); // GriddedRegion region = new GriddedRegion(new Location(32.5,-118), new Location(35.5,-114), 0.1, null); numLocs = region.getNodeCount(); diff --git a/src/main/java/scratch/ned/GK_Declustering/U3ETAS_LossSimulationAnalysis.java b/src/main/java/scratch/ned/GK_Declustering/U3ETAS_LossSimulationAnalysis.java index 3a50591c..c1014b66 100644 --- a/src/main/java/scratch/ned/GK_Declustering/U3ETAS_LossSimulationAnalysis.java +++ b/src/main/java/scratch/ned/GK_Declustering/U3ETAS_LossSimulationAnalysis.java @@ -4,68 +4,47 @@ import java.io.File; import java.io.IOException; import java.util.ArrayList; -import java.util.GregorianCalendar; import java.util.HashMap; import java.util.List; -import java.util.Random; import org.apache.commons.math3.distribution.LogNormalDistribution; import org.apache.commons.math3.distribution.PoissonDistribution; import org.apache.commons.math3.random.JDKRandomGenerator; -import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics; import org.dom4j.DocumentException; import org.jfree.data.Range; import org.opensha.commons.data.CSVFile; -import org.opensha.commons.data.Site; import org.opensha.commons.data.TimeSpan; -import org.opensha.commons.data.function.AbstractDiscretizedFunc; import org.opensha.commons.data.function.ArbDiscrEmpiricalDistFunc; import org.opensha.commons.data.function.ArbDiscrEmpiricalDistFunc_3D; import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc; -import org.opensha.commons.data.function.DefaultXY_DataSet; import org.opensha.commons.data.function.DiscretizedFunc; import org.opensha.commons.data.function.EvenlyDiscretizedFunc; import org.opensha.commons.data.function.HistogramFunction; import org.opensha.commons.data.function.XY_DataSet; -import org.opensha.commons.data.region.CaliforniaRegions; import org.opensha.commons.data.uncertainty.UncertainArbDiscFunc; -import org.opensha.commons.geo.GriddedRegion; -import org.opensha.commons.geo.Location; import org.opensha.commons.geo.LocationUtils; import org.opensha.commons.gui.plot.PlotCurveCharacterstics; import org.opensha.commons.gui.plot.PlotLineType; -import org.opensha.commons.gui.plot.PlotSymbol; -import org.opensha.commons.param.Parameter; -import org.opensha.commons.param.ParameterList; import org.opensha.commons.util.ExceptionUtils; -import org.opensha.sha.calc.HazardCurveCalculator; -import org.opensha.sha.calc.params.PtSrcDistanceCorrectionParam; -import org.opensha.sha.earthquake.ERF; import org.opensha.sha.earthquake.EqkRupture; import org.opensha.sha.earthquake.ProbEqkRupture; -import org.opensha.sha.earthquake.calc.ERF_Calculator; import org.opensha.sha.earthquake.calc.recurInterval.LognormalDistCalc; import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; -import org.opensha.sha.earthquake.faultSysSolution.modules.RupMFDsModule; import org.opensha.sha.earthquake.observedEarthquake.ObsEqkRupList; import org.opensha.sha.earthquake.observedEarthquake.ObsEqkRupture; import org.opensha.sha.earthquake.observedEarthquake.Declustering.GardnerKnopoffDeclustering; -import org.opensha.sha.earthquake.param.MagDependentAperiodicityParam; import org.opensha.sha.earthquake.param.ProbabilityModelOptions; import org.opensha.sha.earthquake.param.ProbabilityModelParam; import org.opensha.sha.magdist.IncrementalMagFreqDist; import org.opensha.sha.magdist.SummedMagFreqDist; -import com.google.common.base.Preconditions; -import com.google.common.collect.HashBasedTable; import com.google.common.math.Quantiles; -import scratch.UCERF3.erf.ETAS.ETAS_Utils; import scratch.UCERF3.erf.FaultSystemSolutionERF; import scratch.UCERF3.erf.ETAS.ETAS_CatalogIO; import scratch.UCERF3.erf.ETAS.ETAS_CatalogIO.ETAS_Catalog; import scratch.UCERF3.erf.ETAS.ETAS_EqkRupture; -import scratch.UCERF3.erf.ETAS.ETAS_Simulator; +import scratch.UCERF3.erf.ETAS.ETAS_Utils; import scratch.UCERF3.erf.ETAS.FaultSystemSolutionERF_ETAS; import scratch.UCERF3.erf.ETAS.launcher.ETAS_Launcher; import scratch.UCERF3.griddedSeismicity.AbstractGridSourceProvider; @@ -1242,8 +1221,6 @@ public ArrayList old_makeCatalogMFDs(boolean mkPopUpPlots) { public static FaultSystemSolutionERF_ETAS getERF() { - // temporary hack (this needs to be done first) - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = 2.55; FaultSystemSolution sol=null; try { sol = FaultSystemSolution.load(new File(fssFileName)); diff --git a/src/main/java/scratch/ned/GK_Declustering/U3ETAS_SimulationAnalysis.java b/src/main/java/scratch/ned/GK_Declustering/U3ETAS_SimulationAnalysis.java index f7992d0b..6b1973a2 100644 --- a/src/main/java/scratch/ned/GK_Declustering/U3ETAS_SimulationAnalysis.java +++ b/src/main/java/scratch/ned/GK_Declustering/U3ETAS_SimulationAnalysis.java @@ -1,31 +1,21 @@ package scratch.ned.GK_Declustering; import java.awt.Color; -import java.io.BufferedInputStream; -import java.io.DataInputStream; import java.io.File; import java.io.IOException; -import java.io.InputStream; import java.util.ArrayList; -import java.util.Enumeration; -import java.util.GregorianCalendar; import java.util.List; import java.util.Random; -import java.util.TimeZone; -import java.util.zip.GZIPInputStream; -import java.util.zip.ZipEntry; -import java.util.zip.ZipFile; import org.apache.commons.math3.distribution.PoissonDistribution; import org.dom4j.DocumentException; import org.jfree.data.Range; -import org.opensha.commons.data.CSVFile; import org.opensha.commons.data.Site; +import org.opensha.commons.data.WeightedList; import org.opensha.commons.data.function.AbstractDiscretizedFunc; import org.opensha.commons.data.function.ArbDiscrEmpiricalDistFunc; import org.opensha.commons.data.function.ArbDiscrEmpiricalDistFunc_3D; import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc; -import org.opensha.commons.data.function.DefaultXY_DataSet; import org.opensha.commons.data.function.DiscretizedFunc; import org.opensha.commons.data.function.EvenlyDiscretizedFunc; import org.opensha.commons.data.function.HistogramFunction; @@ -36,33 +26,26 @@ import org.opensha.commons.geo.Location; import org.opensha.commons.gui.plot.PlotCurveCharacterstics; import org.opensha.commons.gui.plot.PlotLineType; -import org.opensha.commons.gui.plot.PlotSymbol; import org.opensha.commons.param.Parameter; import org.opensha.commons.param.ParameterList; import org.opensha.commons.util.ExceptionUtils; import org.opensha.sha.calc.HazardCurveCalculator; -import org.opensha.sha.calc.params.PtSrcDistanceCorrectionParam; import org.opensha.sha.earthquake.EqkRupture; -import org.opensha.sha.earthquake.calc.ERF_Calculator; import org.opensha.sha.earthquake.faultSysSolution.FaultSystemSolution; import org.opensha.sha.earthquake.observedEarthquake.ObsEqkRupList; import org.opensha.sha.earthquake.observedEarthquake.ObsEqkRupture; import org.opensha.sha.earthquake.observedEarthquake.Declustering.GardnerKnopoffDeclustering; -import org.opensha.sha.earthquake.observedEarthquake.parsers.USGS_NSHMP_CatalogParser; import org.opensha.sha.earthquake.param.AleatoryMagAreaStdDevParam; import org.opensha.sha.earthquake.param.ApplyGardnerKnopoffAftershockFilterParam; -import org.opensha.sha.earthquake.param.BPTAveragingTypeOptions; -import org.opensha.sha.earthquake.param.BPTAveragingTypeParam; import org.opensha.sha.earthquake.param.BackgroundRupParam; import org.opensha.sha.earthquake.param.BackgroundRupType; -import org.opensha.sha.earthquake.param.HistoricOpenIntervalParam; import org.opensha.sha.earthquake.param.IncludeBackgroundOption; import org.opensha.sha.earthquake.param.IncludeBackgroundParam; -import org.opensha.sha.earthquake.param.MagDependentAperiodicityOptions; -import org.opensha.sha.earthquake.param.MagDependentAperiodicityParam; import org.opensha.sha.earthquake.param.ProbabilityModelOptions; import org.opensha.sha.earthquake.param.ProbabilityModelParam; -import org.opensha.sha.faultSurface.utils.PtSrcDistCorr; +import org.opensha.sha.faultSurface.PointSurface; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrection; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrections; import org.opensha.sha.imr.AttenRelRef; import org.opensha.sha.imr.ScalarIMR; import org.opensha.sha.imr.param.IntensityMeasureParams.PGA_Param; @@ -70,7 +53,6 @@ import org.opensha.sha.imr.param.OtherParams.SigmaTruncLevelParam; import org.opensha.sha.imr.param.OtherParams.SigmaTruncTypeParam; import org.opensha.sha.magdist.IncrementalMagFreqDist; -import org.opensha.sha.magdist.SummedMagFreqDist; import com.google.common.base.Preconditions; import com.google.common.math.Quantiles; @@ -79,16 +61,12 @@ import scratch.UCERF3.erf.ETAS.ETAS_CatalogIO; import scratch.UCERF3.erf.ETAS.ETAS_CatalogIO.ETAS_Catalog; import scratch.UCERF3.erf.ETAS.ETAS_EqkRupture; -import scratch.UCERF3.erf.ETAS.ETAS_SimAnalysisTools; import scratch.UCERF3.erf.ETAS.FaultSystemSolutionERF_ETAS; import scratch.UCERF3.erf.ETAS.launcher.ETAS_Launcher; -import scratch.UCERF3.erf.utils.ProbabilityModelsCalc; import scratch.UCERF3.griddedSeismicity.AbstractGridSourceProvider; -import scratch.UCERF3.utils.U3FaultSystemIO; import scratch.UCERF3.utils.GardnerKnopoffAftershockFilter; -import scratch.ned.FSS_Inversion2019.FaultSystemRuptureRateInversion; +import scratch.UCERF3.utils.U3FaultSystemIO; import scratch.ned.FSS_Inversion2019.PlottingUtils; -import scratch.ned.GK_Declustering.MakeFigures; public class U3ETAS_SimulationAnalysis { @@ -361,9 +339,6 @@ public ArrayList makeMFDs(boolean mkPopUpPlots) { public static ArrayList loadCatalogs(File fssFile, File catalogsFile, double minMag) throws IOException, DocumentException { - - // temporary hack - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = 2.55; FaultSystemSolution sol = U3FaultSystemIO.loadSol(fssFile);; List catalogs = ETAS_CatalogIO.loadCatalogsBinary(catalogsFile, minMag); @@ -1214,9 +1189,14 @@ private static EvenlyDiscretizedFunc computeHazardCurveLnX(ObsEqkRupList obsQkLi HazardCurveCalculator calc = getHazardCurveCalculator(); + WeightedList distCorrs = PointSourceDistanceCorrections.NSHM_2013.get(); + ArrayList eqkRupList = new ArrayList(); - for(ObsEqkRupture rup: obsQkList) + for(ObsEqkRupture rup: obsQkList) { + if (rup.getRuptureSurface() instanceof PointSurface) + ((PointSurface)rup.getRuptureSurface()).setDistanceCorrection(distCorrs.sample(), rup); eqkRupList.add(rup); + } if(randomIML) calc.getEventSetHazardCurveRandomIML(curveLogXvalues, site, imr, eqkRupList, false, random); @@ -1235,7 +1215,6 @@ private static EvenlyDiscretizedFunc computeHazardCurveLnX(ObsEqkRupList obsQkLi */ private static HazardCurveCalculator getHazardCurveCalculator() { HazardCurveCalculator calc = new HazardCurveCalculator(); - calc.setPtSrcDistCorrType(PtSrcDistCorr.Type.NSHMP08); calc.setMinMagnitude(5.0); return calc; } @@ -1271,16 +1250,20 @@ private static EvenlyDiscretizedFunc computeExpNumExceedCurveLnX(ObsEqkRupList o HazardCurveCalculator calc = getHazardCurveCalculator(); + WeightedList distCorr = PointSourceDistanceCorrections.NSHM_2013.get(); + ArrayList eqkRupList = new ArrayList(); - for(ObsEqkRupture rup: obsQkList) + for(ObsEqkRupture rup: obsQkList) { + if (rup.getRuptureSurface() instanceof PointSurface) + ((PointSurface)rup.getRuptureSurface()).setDistanceCorrection(distCorr.sample(), rup); eqkRupList.add(rup); + } calc.getEventSetExpNumExceedCurve(curveLogXvalues, site, imr, eqkRupList, false); return curveLogXvalues; } - /** */ private EvenlyDiscretizedFunc computeNumExceedCurveRandomIML_LnX(ObsEqkRupList obsQkList, Location location, @@ -1311,9 +1294,14 @@ private EvenlyDiscretizedFunc computeNumExceedCurveRandomIML_LnX(ObsEqkRupList o HazardCurveCalculator calc = getHazardCurveCalculator(); + WeightedList distCorr = PointSourceDistanceCorrections.NSHM_2013.get(); + ArrayList eqkRupList = new ArrayList(); - for(ObsEqkRupture rup: obsQkList) + for(ObsEqkRupture rup: obsQkList) { + if (rup.getRuptureSurface() instanceof PointSurface) + ((PointSurface)rup.getRuptureSurface()).setDistanceCorrection(distCorr.sample(), rup); eqkRupList.add(rup); + } calc.getEventSetNumExceedCurveRandomIML(curveLogXvalues, site, imr, eqkRupList, false, random); diff --git a/src/test/java/scratch/kevin/ucerf3/eal/UCER3_EAL_CombinerTest.java b/src/test/java/scratch/kevin/ucerf3/eal/UCER3_EAL_CombinerTest.java index ae5721a7..2cc0fce1 100644 --- a/src/test/java/scratch/kevin/ucerf3/eal/UCER3_EAL_CombinerTest.java +++ b/src/test/java/scratch/kevin/ucerf3/eal/UCER3_EAL_CombinerTest.java @@ -31,7 +31,11 @@ import org.opensha.sha.earthquake.param.BackgroundRupType; import org.opensha.sha.earthquake.param.IncludeBackgroundOption; import org.opensha.sha.earthquake.param.IncludeBackgroundParam; +import org.opensha.sha.earthquake.param.PointSourceDistanceCorrectionParam; +import org.opensha.sha.earthquake.util.GriddedSeismicitySettings; import org.opensha.sha.faultSurface.FaultSection; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrection; +import org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrections; import org.opensha.sha.magdist.IncrementalMagFreqDist; import org.opensha.sha.magdist.SummedMagFreqDist; import org.opensha.sha.util.FocalMech; @@ -228,9 +232,9 @@ public static void setUpBeforeClass() throws Exception { ((SummedMagFreqDist)meanSubMFD).addIncrementalMagFreqDist(subMFD); } - if (gridProv.getMFD(n, AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF).calcSumOfY_Vals()>0) + if (gridProv.getMFD(n, 5d).calcSumOfY_Vals()>0) numAboveZero++; - Preconditions.checkState(gridProv.getMFD(n, AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF).size() > 0); + Preconditions.checkState(gridProv.getMFD(n, 5d).size() > 0); } Preconditions.checkState(numAboveZero>0, "Sol "+i+" has all zero mfd nodes!"); @@ -425,6 +429,8 @@ public void doTestGridded(BackgroundRupType bgType) throws DocumentException, IO double[] griddedEALs = comb.getGriddedEALs(); List branches = comb.getBranches(); + GriddedSeismicitySettings gridSettings = erf.getGriddedSeismicitySettings(); + assertEquals(branches.size(), numSols); assertEquals(branches.size(), griddedEALs.length); @@ -437,11 +443,11 @@ public void doTestGridded(BackgroundRupType bgType) throws DocumentException, IO int numNonZeroLoss = 0; int numMFDNonZero = 0; for (int n=0; n 0) + if (gridProv.getMFD(n, gridSettings.minimumMagnitude).calcSumOfY_Vals() > 0) numMFDNonZero++; for (ProbEqkRupture rup : src) { double loss = lossForGridRup(rup, n);