From 23ad85a37e31c4df69f2136e0753a153ef40d392 Mon Sep 17 00:00:00 2001 From: Kevin Milner Date: Wed, 29 Nov 2023 16:26:04 -0800 Subject: [PATCH] simulator experiments separating crustal and subduction rups --- src/main/java/scratch/kevin/MemTest.java | 36 + .../java/scratch/kevin/bbp/BBP_Module.java | 4 +- .../MultiRupGMPE_ComparePageGen.java | 22 +- .../simCompare/SimulationHazardPlotter.java | 9 +- .../kevin/simCompare/ZScoreHistPlot.java | 7 +- .../simulators/ruptures/ASK_EventData.java | 74 ++ .../ruptures/CatalogGMPE_Compare.java | 54 +- .../ruptures/MPJ_BBP_CatalogSimScriptGen.java | 51 +- .../ruptures/RSQSimGeographicMapMaker.java | 244 +++++ .../subduction/BruceSubductionFilter.java | 4 +- .../CrustalAndSubductionCSVWriter.java | 440 -------- .../CrustalAndSubductionPageGen.java | 963 ++++++++++++++++++ .../scratch/kevin/ucerf3/PureScratch.java | 84 +- 13 files changed, 1485 insertions(+), 507 deletions(-) create mode 100644 src/main/java/scratch/kevin/MemTest.java create mode 100644 src/main/java/scratch/kevin/simulators/ruptures/RSQSimGeographicMapMaker.java delete mode 100644 src/main/java/scratch/kevin/simulators/ruptures/subduction/CrustalAndSubductionCSVWriter.java create mode 100644 src/main/java/scratch/kevin/simulators/ruptures/subduction/CrustalAndSubductionPageGen.java diff --git a/src/main/java/scratch/kevin/MemTest.java b/src/main/java/scratch/kevin/MemTest.java new file mode 100644 index 00000000..f5e09b36 --- /dev/null +++ b/src/main/java/scratch/kevin/MemTest.java @@ -0,0 +1,36 @@ +package scratch.kevin; + +import java.text.DecimalFormat; +import java.util.Random; + +import com.google.common.base.Preconditions; + +public class MemTest { + + public static void main(String[] args) { + long maxMemBytes = Runtime.getRuntime().maxMemory(); + long gbToBytes = 1024l*1024l*1024l; + double maxMemGB = (double)maxMemBytes/(double)gbToBytes; + long curMemBytes = Runtime.getRuntime().totalMemory(); + double curMemGB = (double)curMemBytes/(double)gbToBytes; + DecimalFormat gbDF = new DecimalFormat("0.00"); + System.gc(); + System.out.println("Max VM memory: "+gbDF.format(maxMemGB)+" GB"); + System.out.println("Starting VM memory: "+gbDF.format(curMemGB)+" GB"); + int targetGB = Integer.parseInt(args[0]); + long targetBytes = gbToBytes*targetGB; + long numLongs = targetBytes / 8l; + System.out.println("Will try to allocate "+targetGB+" GB = "+targetBytes+" bytes = "+numLongs+" longs"); + Preconditions.checkState(numLongs < (long)Integer.MAX_VALUE); + long[] hugeArray = new long[(int)numLongs]; + System.out.println("Done initializing array"); + curMemBytes = Runtime.getRuntime().totalMemory(); + curMemGB = (double)curMemBytes/(double)gbToBytes; + System.out.println("Current VM memory: "+gbDF.format(curMemGB)); + System.gc(); + Random r = new Random(); + while (true) + hugeArray[r.nextInt(hugeArray.length)] = r.nextLong(); + } + +} diff --git a/src/main/java/scratch/kevin/bbp/BBP_Module.java b/src/main/java/scratch/kevin/bbp/BBP_Module.java index a27eb18c..cd0f6da0 100644 --- a/src/main/java/scratch/kevin/bbp/BBP_Module.java +++ b/src/main/java/scratch/kevin/bbp/BBP_Module.java @@ -19,7 +19,9 @@ public static enum VelocityModel { LA_BASIN_863("LA Basin 863 (m/s)", "LABasin863", "$BBP_INSTALL_GF/LABasin863/gp/genslip_nr_generic1d-gp01.vmod", 863, 0.02, 3), LA_BASIN_500("LA Basin 500 (m/s)", "LABasin500", - "$BBP_INSTALL_GF/LABasin500/gp/nr02-vs500.fk1d", 500, 0.2, 2.5); + "$BBP_INSTALL_GF/LABasin500/gp/nr02-vs500.fk1d", 500, 0.2, 2.5), + CENTRAL_JAPAN("Central Japan", "CentralJapan500", + "$BBP_INSTALL_GF/CentralJapan500/gp/niigata04-vs500.fk1d", 500, 0.1, 0.9); private String name; private String xmlName; diff --git a/src/main/java/scratch/kevin/simCompare/MultiRupGMPE_ComparePageGen.java b/src/main/java/scratch/kevin/simCompare/MultiRupGMPE_ComparePageGen.java index 3028fb77..5ac2ff05 100644 --- a/src/main/java/scratch/kevin/simCompare/MultiRupGMPE_ComparePageGen.java +++ b/src/main/java/scratch/kevin/simCompare/MultiRupGMPE_ComparePageGen.java @@ -19,6 +19,7 @@ import java.util.concurrent.ExecutorService; import java.util.concurrent.Executors; import java.util.concurrent.Future; +import java.util.function.Supplier; import org.apache.commons.math3.stat.descriptive.moment.StandardDeviation; import org.jfree.chart.annotations.XYPolygonAnnotation; @@ -55,6 +56,7 @@ import org.opensha.sha.earthquake.faultSysSolution.ruptures.util.RupSetMapMaker; import org.opensha.sha.faultSurface.FaultSection; import org.opensha.sha.imr.AttenRelRef; +import org.opensha.sha.imr.AttenRelSupplier; import org.opensha.sha.imr.ScalarIMR; import org.opensha.sha.imr.attenRelImpl.MultiIMR_Averaged_AttenRel; import org.opensha.sha.imr.param.EqkRuptureParams.MagParam; @@ -108,7 +110,7 @@ public abstract class MultiRupGMPE_ComparePageGen { protected ExecutorService exec; - private Map> gmpesInstancesCache; + private Map, LinkedList> gmpesInstancesCache; private boolean replotScatters = true; private boolean replotZScores = true; @@ -335,7 +337,7 @@ public boolean matches(RuptureComparison comp, Site site) { private static final int MAX_SCATTER_NUM_POINTS = 100000; public boolean plotScatter(Collection> eventComps, Collection sites, IMT imt, - AttenRelRef gmpe, RuptureComparisonFilter filter, List binDescriptions, File outputDir, String prefix) + AttenRelSupplier gmpe, RuptureComparisonFilter filter, List binDescriptions, File outputDir, String prefix) throws IOException { DefaultXY_DataSet xy = new DefaultXY_DataSet(); @@ -451,7 +453,7 @@ public File call() throws Exception { } - protected ScalarIMR checkOutGMPE(AttenRelRef gmpeRef) { + protected ScalarIMR checkOutGMPE(Supplier gmpeRef) { synchronized (gmpesInstancesCache) { LinkedList gmpes = gmpesInstancesCache.get(gmpeRef); if (gmpes == null) { @@ -461,12 +463,12 @@ protected ScalarIMR checkOutGMPE(AttenRelRef gmpeRef) { if (!gmpes.isEmpty()) return gmpes.pop(); } - ScalarIMR gmpe = gmpeRef.instance(null); + ScalarIMR gmpe = gmpeRef.get(); gmpe.setParamDefaults(); return gmpe; } - protected void checkInGMPE(AttenRelRef gmpeRef, ScalarIMR gmpe) { + protected void checkInGMPE(Supplier gmpeRef, ScalarIMR gmpe) { synchronized (gmpesInstancesCache) { gmpesInstancesCache.get(gmpeRef).push(gmpe); } @@ -507,7 +509,11 @@ private static List getZscoreTableLine(String siteName, ZScoreResult[] s return line; } - protected void checkTransSiteParams(AttenRelRef gmpeRef) { + protected void checkTransSiteParams(Supplier gmpeRef) { + checkTransSiteParams(gmpeRef, sites); + } + + public static void checkTransSiteParams(Supplier gmpeRef, List sites) { SiteTranslator trans = null; ScalarIMR gmmInstance = gmpeRef.get(); Site site0 = sites.get(0); @@ -543,7 +549,7 @@ protected void checkTransSiteParams(AttenRelRef gmpeRef) { } } - public void generateGMPE_Page(File outputDir, List headerLines, AttenRelRef gmpeRef, IMT[] imts, + public void generateGMPE_Page(File outputDir, List headerLines, AttenRelSupplier gmpeRef, IMT[] imts, List> comps, List highlightSites) throws IOException { File resourcesDir = new File(outputDir, "resources"); Preconditions.checkState(resourcesDir.exists() || resourcesDir.mkdir()); @@ -1806,7 +1812,7 @@ private static enum HazardRatioType { } public void generateNonErgodicMapPage(File outputDir, List headerLines, List subSects, - Region siteRegion, Region sourceRegion, AttenRelRef gmpeRef, List> comps, + Region siteRegion, Region sourceRegion, AttenRelSupplier gmpeRef, List> comps, Map> rupSectMappings, Map rupNuclSects, IMT[] imts, List highlightSites) throws IOException { File resourcesDir = new File(outputDir, "resources"); diff --git a/src/main/java/scratch/kevin/simCompare/SimulationHazardPlotter.java b/src/main/java/scratch/kevin/simCompare/SimulationHazardPlotter.java index 69a0bf88..535c8e27 100644 --- a/src/main/java/scratch/kevin/simCompare/SimulationHazardPlotter.java +++ b/src/main/java/scratch/kevin/simCompare/SimulationHazardPlotter.java @@ -42,6 +42,7 @@ import org.opensha.commons.util.FileUtils; import org.opensha.commons.util.cpt.CPT; import org.opensha.sha.imr.AttenRelRef; +import org.opensha.sha.imr.AttenRelSupplier; import org.opensha.sha.imr.param.IntensityMeasureParams.PGA_Param; import org.opensha.sha.imr.param.IntensityMeasureParams.PGV_Param; import org.opensha.sha.imr.param.IntensityMeasureParams.SA_Param; @@ -100,7 +101,7 @@ public class SimulationHazardPlotter { private List> comps; private Site site; private double curveDuration; - private AttenRelRef gmpeRef; + private AttenRelSupplier gmpeRef; private Map xValsMap; private Table sourceRupContribFracts; @@ -121,12 +122,12 @@ public class SimulationHazardPlotter { private Table gmpeSourceCurves = HashBasedTable.create(); public SimulationHazardPlotter(SimulationHazardCurveCalc simCalc, List> comps, - Site site, double curveDuration, AttenRelRef gmpeRef) { + Site site, double curveDuration, AttenRelSupplier gmpeRef) { this(simCalc, null, comps, site, curveDuration, gmpeRef); } public SimulationHazardPlotter(SimulationHazardCurveCalc simCalc, List> compCalcs, - List> comps, Site site, double curveDuration, AttenRelRef gmpeRef) { + List> comps, Site site, double curveDuration, AttenRelSupplier gmpeRef) { this.simCalc = simCalc; if (compCalcs == null) compCalcs = new ArrayList<>(); @@ -148,7 +149,7 @@ public Site getSite() { return site; } - public AttenRelRef getGMPE() { + public AttenRelSupplier getGMPE() { return gmpeRef; } diff --git a/src/main/java/scratch/kevin/simCompare/ZScoreHistPlot.java b/src/main/java/scratch/kevin/simCompare/ZScoreHistPlot.java index 99d01fde..6e9ab773 100644 --- a/src/main/java/scratch/kevin/simCompare/ZScoreHistPlot.java +++ b/src/main/java/scratch/kevin/simCompare/ZScoreHistPlot.java @@ -36,6 +36,7 @@ import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; import org.opensha.commons.util.cpt.CPT; import org.opensha.sha.imr.AttenRelRef; +import org.opensha.sha.imr.AttenRelSupplier; import com.google.common.base.Preconditions; import com.google.common.collect.Table; @@ -188,7 +189,7 @@ else if (numMatches < 500) } public static boolean plotStandardNormal(SimulationRotDProvider simProv, Collection> eventComps, - List sites, IMT[] imts, AttenRelRef gmpe, RuptureComparisonFilter filter, List binDescriptions, + List sites, IMT[] imts, AttenRelSupplier gmpe, RuptureComparisonFilter filter, List binDescriptions, File outputDir, String prefix) throws IOException { return plotStandardNormal(simProv, eventComps, sites, imts, gmpe, filter, binDescriptions, outputDir, prefix, null, 0, false); } @@ -197,7 +198,7 @@ public static boolean plotStandardNormal(SimulationRotDProvider simProv, public static final double numStdDev = 3.75; public static boolean plotStandardNormal(SimulationRotDProvider simProv, Collection> eventComps, - List sites, IMT[] imts, AttenRelRef gmpe, RuptureComparisonFilter filter, List binDescriptions, + List sites, IMT[] imts, AttenRelSupplier gmpe, RuptureComparisonFilter filter, List binDescriptions, File outputDir, String prefix, Table sourceRupContribFracts, int maxNumSourceContribs, boolean pub) throws IOException { @@ -208,7 +209,7 @@ public static boolean plotStandardNormal(SimulationRotDProvider simProv, return plotStandardNormal(scores, simProv.getName(), imts, gmpe, binDescriptions, outputDir, prefix, maxNumSourceContribs, pub, maxY); } - public static boolean plotStandardNormal(ZScoreResult[] scores, String name, IMT[] imts, AttenRelRef gmpe, + public static boolean plotStandardNormal(ZScoreResult[] scores, String name, IMT[] imts, AttenRelSupplier gmpe, List binDescriptions, File outputDir, String prefix, int maxNumSourceContribs, boolean pub, double maxY) throws IOException { diff --git a/src/main/java/scratch/kevin/simulators/ruptures/ASK_EventData.java b/src/main/java/scratch/kevin/simulators/ruptures/ASK_EventData.java index e286d99c..54f169cb 100644 --- a/src/main/java/scratch/kevin/simulators/ruptures/ASK_EventData.java +++ b/src/main/java/scratch/kevin/simulators/ruptures/ASK_EventData.java @@ -18,11 +18,16 @@ import org.opensha.commons.data.CSVFile; import org.opensha.commons.data.function.DiscretizedFunc; import org.opensha.commons.data.function.HistogramFunction; +import org.opensha.commons.data.xyz.EvenlyDiscrXYZ_DataSet; import org.opensha.commons.gui.plot.HeadlessGraphPanel; import org.opensha.commons.gui.plot.PlotCurveCharacterstics; import org.opensha.commons.gui.plot.PlotLineType; import org.opensha.commons.gui.plot.PlotSpec; +import org.opensha.commons.gui.plot.PlotUtils; +import org.opensha.commons.gui.plot.jfreechart.xyzPlot.XYZPlotSpec; +import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; import org.opensha.commons.util.DataUtils.MinMaxAveTracker; +import org.opensha.commons.util.cpt.CPT; import org.opensha.sha.imr.attenRelImpl.ngaw2.FaultStyle; import com.google.common.base.Preconditions; @@ -315,6 +320,64 @@ public static void plotMultiCountHist(File outputDir, String prefix, String titl gp.saveAsPDF(file.getAbsolutePath()+".pdf"); } + public static void plotMagDistXYZ(File outputDir, String prefix, String title, + Collection> data, double minMag, double maxMag, double minDist, double maxDist, + boolean rJB, boolean eventCount, boolean logZ) throws IOException { + int numDist = 20; + int numMag = 15; + double distSpacing = (Math.log10(maxDist) - Math.log10(minDist)) / (double)numDist; + double magSpacing = (maxMag - minMag) / (double)numMag; + EvenlyDiscrXYZ_DataSet xyz = new EvenlyDiscrXYZ_DataSet(numDist, numMag, Math.log10(minDist)+0.5*distSpacing, + minMag+0.5*magSpacing, distSpacing, magSpacing); + + for (List eventData : data) { + double mag = eventData.get(0).mag; + if (mag < minMag || mag > maxMag) + continue; + int magBin = xyz.getYIndex(mag); + int[] distCounts = new int[xyz.getNumX()]; + for (ASK_EventData event : eventData) { + double dist = rJB ? event.rJB : event.rRup; + if (dist < minDist || dist > maxDist) + continue; + double logDist = Math.log10(dist); + int distBin = xyz.getXIndex(logDist); + if (eventCount) + distCounts[distBin] = 1; + else + distCounts[distBin]++; + } + for (int x=0; x 0) + xyz.set(x, magBin, xyz.get(x, magBin)+distCounts[x]); + } + + CPT cpt = GMT_CPT_Files.BLACK_RED_YELLOW_UNIFORM.instance().reverse(); + if (logZ) { + xyz.log10(); + cpt = cpt.rescale(0d, Math.ceil(xyz.getMaxZ())); + cpt.setNanColor(Color.WHITE); + } else { + cpt = cpt.rescale(1d, xyz.getMaxZ()); + } + cpt.setBelowMinColor(Color.WHITE); + + String xLabel = rJB ? "Log10 Distance JB" : "Log10 Distance Rup"; + String yLabel = "Magnitude"; + String zLabel = eventCount ? "Num Events" : "Num Recordings"; + if (logZ) + zLabel = "Lot10 "+zLabel; + + XYZPlotSpec xyzSpec = new XYZPlotSpec(xyz, cpt, title, xLabel, yLabel, zLabel); + xyzSpec.setCPTTickUnit(0.5d); + HeadlessGraphPanel xyzGP = PlotUtils.initHeadless(); + xyzGP.drawGraphPanel(xyzSpec, false, false, new org.jfree.data.Range(Math.log10(minDist), Math.log10(maxDist)), + new org.jfree.data.Range(minMag, maxMag)); +// xyzGP.getChartPanel().getChart().setBackgroundPaint(Color.WHITE); + xyzGP.getChartPanel().setSize(700, 550); + PlotUtils.writePlots(outputDir, prefix, xyzGP, 700, 550, true, true, false); + } + private static final DecimalFormat twoDigit = new DecimalFormat("0.00"); public static void main(String[] args) throws IOException { @@ -361,6 +424,17 @@ public static void main(String[] args) throws IOException { System.out.println("Loaded "+totNum+" total records"); System.out.println(); periodData.put(period, data); + + double minDist = 1; + double maxDist = 200; + double minMag = 6.5d; + double maxMag = 8.2d; + + String perStr = new DecimalFormat("0.#").format(period); + String title = "ASK (2014) Data, "+perStr+"s"; + + plotMagDistXYZ(new File("/tmp"), "ask_mag_dist_"+perStr+"s", title, data.values(), + minMag, maxMag, minDist, maxDist, false, false, true); } for (Range magRange : magRanges) { diff --git a/src/main/java/scratch/kevin/simulators/ruptures/CatalogGMPE_Compare.java b/src/main/java/scratch/kevin/simulators/ruptures/CatalogGMPE_Compare.java index f034a861..aee0ff8f 100644 --- a/src/main/java/scratch/kevin/simulators/ruptures/CatalogGMPE_Compare.java +++ b/src/main/java/scratch/kevin/simulators/ruptures/CatalogGMPE_Compare.java @@ -15,6 +15,7 @@ import java.util.Set; import java.util.concurrent.ExecutionException; import java.util.concurrent.Future; +import java.util.function.Supplier; import java.util.zip.ZipException; import java.util.zip.ZipFile; @@ -42,6 +43,7 @@ import org.opensha.sha.earthquake.EqkRupture; import org.opensha.sha.faultSurface.FaultSection; import org.opensha.sha.imr.AttenRelRef; +import org.opensha.sha.imr.AttenRelSupplier; import org.opensha.sha.imr.IntensityMeasureRelationship; import org.opensha.sha.imr.ScalarIMR; import org.opensha.sha.imr.param.IntensityMeasureParams.PGV_Param; @@ -216,7 +218,7 @@ private EqkRupture getGMPE_Rup(RSQSimEvent event) { } } - private int loadCache(File file, Site site, Map eventComps) throws IOException { + private static int loadCache(File file, Site site, Map eventComps) throws IOException { CSVFile csv = CSVFile.readFile(file, true); int count = 0; @@ -241,7 +243,7 @@ private int loadCache(File file, Site site, Map eventC return count; } - private void writeCache(File file, Site site, Map comps) throws IOException { + private static void writeCache(File file, Site site, Map comps) throws IOException { CSVFile csv = new CSVFile<>(true); csv.addLine("Event ID", "Period (s)", "Ln(Mean)", "Std Dev", "rRup (km)", "rJB (km)"); @@ -257,19 +259,20 @@ private void writeCache(File file, Site site, Map comp csv.writeToFile(file); } - private > Iterable sorted(Set vals) { + private static > Iterable sorted(Set vals) { List list = new ArrayList<>(vals); Collections.sort(list); return list; } - public synchronized void calcGMPE(Map eventComps, Site site, AttenRelRef gmpeRef, IMT... imts) { + public synchronized void calcGMPE(Map eventComps, Site site, Supplier gmpeRef, IMT... imts) { List> futures = new ArrayList<>(); File gmpeCacheFile = null; if (gmpeCacheDir != null) { String siteName = sitesGMPEtoBBP.get(site).getName(); - gmpeCacheFile = new File(gmpeCacheDir, siteName+"_"+gmpeRef.getShortName()+".csv"); + String shortName = gmpeRef instanceof AttenRelRef ? ((AttenRelRef)gmpeRef).getShortName() : gmpeRef.get().getShortName(); + gmpeCacheFile = new File(gmpeCacheDir, siteName+"_"+shortName+".csv"); System.out.println(gmpeCacheFile.getAbsolutePath()+" exits ? "+gmpeCacheFile.exists()); if (gmpeCacheFile.exists()) { try { @@ -332,11 +335,11 @@ public synchronized void calcGMPE(Map eventComps, Site private class EventComparisonCalc implements Runnable { private EventComparison eventComp; - private AttenRelRef gmpeRef; + private Supplier gmpeRef; private Site site; private List imts; - public EventComparisonCalc(EventComparison eventComp, AttenRelRef gmpeRef, Site site, List imts) { + public EventComparisonCalc(EventComparison eventComp, Supplier gmpeRef, Site site, List imts) { this.eventComp = eventComp; this.gmpeRef = gmpeRef; this.site = site; @@ -398,7 +401,7 @@ public double getRuptureTimeYears() { } } - public List loadCalcComps(AttenRelRef gmpeRef, IMT[] imts) { + public List loadCalcComps(Supplier gmpeRef, IMT[] imts) { Map compsMap = new HashMap<>(); for (BBP_Site site : sites) { @@ -418,7 +421,7 @@ public List loadCalcComps(AttenRelRef gmpeRef, IMT[] imts) { return new ArrayList<>(compsMap.values()); } - public void generateGMPE_Page(File outputDir, AttenRelRef gmpeRef, IMT[] imts, List comps) + public void generateGMPE_Page(File outputDir, AttenRelSupplier gmpeRef, IMT[] imts, List comps) throws IOException { LinkedList lines = new LinkedList<>(); @@ -566,9 +569,9 @@ public static void main(String[] args) throws ZipException, IOException { // RSQSimCatalog catalog = Catalogs.BRUCE_2585_1MYR.instance(); // RSQSimCatalog catalog = Catalogs.BRUCE_4983_STITCHED.instance(); // RSQSimCatalog catalog = Catalogs.BRUCE_5413.instance(); - RSQSimCatalog catalog = Catalogs.BRUCE_5652.instance(); +// RSQSimCatalog catalog = Catalogs.BRUCE_5652.instance(); // RSQSimCatalog catalog = Catalogs.BRUCE_5566_CRUSTAL.instance(); -// RSQSimCatalog catalog = Catalogs.BRUCE_5585_SUB.instance(); + RSQSimCatalog catalog = Catalogs.BRUCE_5566_SUB.instance(); boolean doGMPE = true; boolean doRotD = false; @@ -586,15 +589,15 @@ public static void main(String[] args) throws ZipException, IOException { // VelocityModel forceVM = VelocityModel.LA_BASIN_863; VelocityModel forceVM = null; - AttenRelRef[] gmpeRefs = { AttenRelRef.NGAWest_2014_AVG_NOIDRISS, AttenRelRef.ASK_2014, - AttenRelRef.BSSA_2014, AttenRelRef.CB_2014, AttenRelRef.CY_2014 }; -// AttenRelRef[] gmpeRefs = { AttenRelRef.NGAWest_2014_AVG_NOIDRISS, AttenRelRef.ASK_2014 }; -// AttenRelRef[] gmpeRefs = { AttenRelRef.NGAWest_2014_AVG_NOIDRISS }; -// AttenRelRef[] gmpeRefs = { AttenRelRef.BSSA_2014, AttenRelRef.CB_2014, AttenRelRef.CY_2014 }; -// IMT[] imts = { IMT.SA3P0 }; -// AttenRelRef[] gmpeRefs = { AttenRelRef.ASK_2014 }; - IMT[] imts = { IMT.PGV, IMT.SA2P0, IMT.SA3P0, IMT.SA5P0, IMT.SA10P0 }; - AttenRelRef rotDGMPE = AttenRelRef.ASK_2014; +// AttenRelSupplier[] gmpeRefs = { AttenRelRef.NGAWest_2014_AVG_NOIDRISS, AttenRelRef.ASK_2014, +// AttenRelRef.BSSA_2014, AttenRelRef.CB_2014, AttenRelRef.CY_2014 }; +//// AttenRelRef[] gmpeRefs = { AttenRelRef.NGAWest_2014_AVG_NOIDRISS, AttenRelRef.ASK_2014 }; +//// AttenRelRef[] gmpeRefs = { AttenRelRef.NGAWest_2014_AVG_NOIDRISS }; +//// AttenRelRef[] gmpeRefs = { AttenRelRef.BSSA_2014, AttenRelRef.CB_2014, AttenRelRef.CY_2014 }; +//// IMT[] imts = { IMT.SA3P0 }; +//// AttenRelRef[] gmpeRefs = { AttenRelRef.ASK_2014 }; +// IMT[] imts = { IMT.PGV, IMT.SA2P0, IMT.SA3P0, IMT.SA5P0, IMT.SA10P0 }; +// AttenRelRef rotDGMPE = AttenRelRef.ASK_2014; // AttenRelRef[] gmpeRefs = { AttenRelRef.AFSHARI_STEWART_2016 }; // IMT[] imts = { IMT.DUR_5_75, IMT.DUR_5_95, IMT.DUR_20_80 }; @@ -604,6 +607,13 @@ public static void main(String[] args) throws ZipException, IOException { // IMT[] imts = { IMT.SA2P0, IMT.SA3P0, IMT.SA5P0 }; // AttenRelRef rotDGMPE = null; + AttenRelSupplier[] gmpeRefs = { + new org.opensha.sha.imr.attenRelImpl.nshmp.NSHMP_AttenRelSupplier( + gov.usgs.earthquake.nshmp.gmm.Gmm.AG_20_GLOBAL_INTERFACE, "AG2020_Global") + }; + IMT[] imts = { IMT.SA2P0, IMT.SA3P0, IMT.SA5P0, IMT.SA10P0 }; + AttenRelRef rotDGMPE = null; + String[] highlightNames; if (doGridded) highlightNames = new String[0]; @@ -748,7 +758,7 @@ public static void main(String[] args) throws ZipException, IOException { System.out.println("Has RotD100? "+comp.bbpZipFile.hasRotD100()); doRotD = doRotD && comp.bbpZipFile.hasRotD100(); - for (AttenRelRef ref : gmpeRefs) + for (Supplier ref : gmpeRefs) comp.checkTransSiteParams(ref); // comp.testPlotRupAzDiffs(); @@ -766,7 +776,7 @@ public static void main(String[] args) throws ZipException, IOException { List rotD_GMPEcomps = null; if (doGMPE || doNonErgodicMaps) { if (!hasRG || !rgOnlyIfPossible) { - for (AttenRelRef gmpeRef : gmpeRefs) { + for (AttenRelSupplier gmpeRef : gmpeRefs) { if (highlightNames != null) { if (gmpeRef == rotDGMPE) comp.setHighlightSites(highlightNames); diff --git a/src/main/java/scratch/kevin/simulators/ruptures/MPJ_BBP_CatalogSimScriptGen.java b/src/main/java/scratch/kevin/simulators/ruptures/MPJ_BBP_CatalogSimScriptGen.java index 16703bee..6b3838bc 100644 --- a/src/main/java/scratch/kevin/simulators/ruptures/MPJ_BBP_CatalogSimScriptGen.java +++ b/src/main/java/scratch/kevin/simulators/ruptures/MPJ_BBP_CatalogSimScriptGen.java @@ -34,10 +34,10 @@ public static void main(String[] args) throws IOException { // String catalogDirName = "rundir5450"; // String catalogDirName = "rundir4983_stitched"; // String catalogDirName = "rundir5566"; - String catalogDirName = "rundir5672"; -// String catalogDirName = "rundir5566_subduction_corupture"; +// String catalogDirName = "rundir5597"; // String catalogDirName = "rundir5413_multifault_separate"; -// String catalogDirName = "rundir5566_crustal_corupture"; +// String catalogDirName = "rundir5597_subduction_corupture"; + String catalogDirName = "rundir5597_crustal_corupture"; // int skipYears = 20000; // int skipYears = 5000; @@ -51,29 +51,33 @@ public static void main(String[] args) throws IOException { double griddedSpacing = 1d; // CA - Integer utmZone = null; - Character utmBand = null; - boolean standardSites = false; - boolean csInitialLASites = false; - boolean cs500LASites = false; - boolean csLAMapSites = false; - boolean griddedCASites = true; - boolean griddedSoCalSites = false; - boolean griddedNZSites = false; - boolean nzStandardSites = false; - - // NZ -// Integer utmZone = 59; -// Character utmBand = 'G'; -// System.out.println("New Zealand!"); +// Integer utmZone = null; +// Character utmBand = null; // boolean standardSites = false; // boolean csInitialLASites = false; // boolean cs500LASites = false; // boolean csLAMapSites = false; -// boolean griddedCASites = false; +// boolean griddedCASites = true; // boolean griddedSoCalSites = false; -// boolean griddedNZSites = true; -// boolean nzStandardSites = true; +// boolean griddedNZSites = false; +// boolean nzStandardSites = false; +//// VelocityModel vm = VelocityModel.LA_BASIN_863; // uncomment only if you need the old 863 +// VelocityModel vm = VelocityModel.LA_BASIN_500; + + // NZ + Integer utmZone = 59; + Character utmBand = 'G'; + System.out.println("New Zealand!"); + boolean standardSites = false; + boolean csInitialLASites = false; + boolean cs500LASites = false; + boolean csLAMapSites = false; + boolean griddedCASites = false; + boolean griddedSoCalSites = false; + boolean griddedNZSites = true; + boolean nzStandardSites = true; + VelocityModel vm = VelocityModel.CENTRAL_JAPAN; +// VelocityModel vm = VelocityModel.LA_BASIN_500; // double minMag = 0d; // double minMag = 5; @@ -84,11 +88,6 @@ public static void main(String[] args) throws IOException { // double minMag = 7; // int numRG = 20; -// RSQSimBBP_Config.VM = VelocityModel.LA_BASIN_863; - if (cs500LASites) - RSQSimBBP_Config.VM = VelocityModel.LA_BASIN_500; - VelocityModel vm = RSQSimBBP_Config.VM; - double timeScalar = 1d; boolean scaleVelocities = true; diff --git a/src/main/java/scratch/kevin/simulators/ruptures/RSQSimGeographicMapMaker.java b/src/main/java/scratch/kevin/simulators/ruptures/RSQSimGeographicMapMaker.java new file mode 100644 index 00000000..2a29b4c6 --- /dev/null +++ b/src/main/java/scratch/kevin/simulators/ruptures/RSQSimGeographicMapMaker.java @@ -0,0 +1,244 @@ +package scratch.kevin.simulators.ruptures; + +import java.awt.BasicStroke; +import java.awt.Color; +import java.util.ArrayList; +import java.util.List; + +import org.jfree.chart.annotations.XYPolygonAnnotation; +import org.jfree.data.Range; +import org.opensha.commons.data.function.DefaultXY_DataSet; +import org.opensha.commons.data.function.XY_DataSet; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.Region; +import org.opensha.commons.gui.plot.GeographicMapMaker; +import org.opensha.commons.gui.plot.PlotCurveCharacterstics; +import org.opensha.commons.gui.plot.PlotLineType; +import org.opensha.commons.gui.plot.PlotUtils; +import org.opensha.commons.util.cpt.CPT; +import org.opensha.sha.simulators.RSQSimEvent; +import org.opensha.sha.simulators.SimulatorElement; +import org.opensha.sha.simulators.Vertex; +import org.opensha.sha.simulators.utils.RSQSimUtils; + +import com.google.common.base.Preconditions; + +public class RSQSimGeographicMapMaker extends GeographicMapMaker { + + private RSQSimEvent event; + private Color elemFillColor; + private List elemScalars; + private CPT elemCPT; + private String elemScalarLabel; + private boolean fillElemScalars; + private Color elemOutlineColor; + private float elemOutlineThickness; + + private Color eventHypoColor = null; + private double eventHypoRadius = 0.2; + + + public RSQSimGeographicMapMaker(Region region) { + super(region); + setWriteGeoJSON(false); + setWritePDFs(false); + } + + public RSQSimGeographicMapMaker(Region region, XY_DataSet[] politicalBoundaries) { + super(region, politicalBoundaries); + setWriteGeoJSON(false); + setWritePDFs(false); + } + + public void plotEvent(RSQSimEvent event, Color fillColor, Color outlineColor, float outlineThickness) { + clearEvent(); + Preconditions.checkState(fillColor != null || outlineColor != null); + this.event = event; + this.elemFillColor = fillColor; + this.elemOutlineColor = outlineColor; + this.elemOutlineThickness = outlineThickness; + } + + public void plotEventFillScalars(RSQSimEvent event, List scalars, CPT cpt, String label) { + plotEventFillScalars(event, scalars, cpt, null, Float.NaN, label); + } + + public void plotEventFillScalars(RSQSimEvent event, List scalars, CPT cpt, Color outlineColor, float outlineThickness, String label) { + clearEvent(); + Preconditions.checkState(scalars.size() == event.getNumElements()); + Preconditions.checkNotNull(cpt); + this.event = event; + this.elemScalars = scalars; + this.elemCPT = cpt; + this.elemScalarLabel = label; + this.elemOutlineColor = outlineColor; + this.elemOutlineThickness = outlineThickness; + this.fillElemScalars = true; + } + + public void plotEventOutlineScalars(RSQSimEvent event, List scalars, CPT cpt, float outlineThickness, String label) { + clearEvent(); + Preconditions.checkState(scalars.size() == event.getNumElements()); + Preconditions.checkNotNull(cpt); + this.event = event; + this.elemScalars = scalars; + this.elemScalarLabel = label; + this.elemCPT = cpt; + this.elemOutlineThickness = outlineThickness; + this.fillElemScalars = false; + } + + private void clearEvent() { + this.event = null; + this.elemFillColor = null; + this.elemScalars = null; + this.elemCPT = null; + this.elemScalarLabel = null; + this.elemOutlineColor = null; + this.elemOutlineThickness = Float.NaN; + } + + public void plotEventHypocenter(Color color) { + plotEventHypocenter(color, eventHypoRadius); + } + + public void plotEventHypocenter(Color color, double radius) { + eventHypoColor = color; + eventHypoRadius = radius; + } + + @Override + protected PlotBuilder initPlotBuilder() { + return new RSQSimPlotBuilder(); + } + + protected class RSQSimPlotBuilder extends PlotBuilder { + + @Override + protected void plotFirst() { + super.plotFirst(); + + if (event != null) { + // plot event + + ArrayList elems = event.getAllElements(); + + XY_DataSet[] elemXYs = buildElemXYs(elems); + + if (elemCPT != null) { + Preconditions.checkState(elems.size() == elemScalars.size()); + + if (fillElemScalars) { + // fill them + for (int i=0; i sites = BBP_Site.readFile(crustalBBPdir); - - HashBiMap sitesBBPtoGMPE = HashBiMap.create(); - List gmpeSites = new ArrayList<>(); - for (BBP_Site site : sites) { - Site gmpeSite = site.buildGMPE_Site(vm); - gmpeSite.setName(RSQSimBBP_Config.siteCleanName(site)); - gmpeSites.add(gmpeSite); - sitesBBPtoGMPE.put(site, gmpeSite); - } - - List events = fullCatalog.loader().minMag(minMag).load(); - Map eventsMap = new HashMap<>(); - for (RSQSimEvent event : events) - eventsMap.put(event.getID(), event); - - Map crustalMap = crustalCatalog.loader().minMag(minMag).loadMap(); - Map subductionMap = subductionCatalog.loader().minMag(minMag).loadMap(); - - BBP_CatalogSimZipLoader fullLoader = new BBP_CatalogSimZipLoader(fullZip, sites, sitesBBPtoGMPE, eventsMap); - BBP_CatalogSimZipLoader crustalLoader = new BBP_CatalogSimZipLoader(curstalZip, sites, sitesBBPtoGMPE, eventsMap); - BBP_CatalogSimZipLoader subductionLoader = new BBP_CatalogSimZipLoader(subductionZip, sites, sitesBBPtoGMPE, eventsMap); - - double[] periods = { 2d, 3d, 5d, 10d }; - - double[] minFracts = { 0d, 0.2, 0.5 }; - - Quantity[] quantities = Quantity.values(); - - DefaultXY_DataSet[][][] scatters = new DefaultXY_DataSet[minFracts.length][quantities.length][periods.length]; - - for (int i=0; i combCSV = null; - - HashSet uniqueEvents = new HashSet<>(); - - List lines = new ArrayList<>(); - - lines.add("# Crustal and Subduction Co-Rupture Ground Motions"); - lines.add(""); - - int tocIndex = lines.size(); - String topLink = "*[(top)](#table-of-contents)*"; - - lines.add("## Sites and CSV Files"); - lines.add(topLink); lines.add(""); - - TableBuilder table = MarkdownUtils.tableBuilder(); - - table.addLine("Site Name", "Latitude", "Longitude", "Record Count", "CSV File"); - - for (BBP_Site site : sites) { - CSVFile csv = new CSVFile<>(true); - - List header = new ArrayList<>(); - header.add("Site Name"); - header.add("Event ID"); - - for (Type type : types) - header.add(type.prefix+" Moment"); - for (Type type : types) - header.add(type.prefix+" Magnitude"); - for (Type type : types) - header.add(type.prefix+" rRup"); - for (Type type : types) - header.add(type.prefix+" rJB"); - for (double period : periods) - for (Type type : types) - header.add(type.prefix+" "+oDF.format(period)+"s RotD50"); - - csv.addLine(header); - - if (combCSV == null) { - combCSV = new CSVFile<>(true); - combCSV.addLine(header); - } - - int matches = 0; - - for (RSQSimEvent fullEvent : events) { - int eventID = fullEvent.getID(); - RSQSimEvent crustalEvent = crustalMap.get(eventID); - RSQSimEvent subductionEvent = subductionMap.get(eventID); - - if (crustalEvent != null && crustalLoader.contains(site, eventID) - && subductionEvent != null && subductionLoader.contains(site, eventID)) { - // exists for both - matches++; - uniqueEvents.add(eventID); - - Map typeEventMap = new EnumMap<>(Type.class); - Map typeEqkRupMap = new EnumMap<>(Type.class); - Map typeRDMap = new EnumMap<>(Type.class); - for (Type type : types) { - RSQSimEvent event; - RSQSimCatalog catalog; - BBP_CatalogSimZipLoader loader; - switch (type) { - case FULL: - event = fullEvent; - catalog = fullCatalog; - loader = fullLoader; - break; - case CRUSTAL: - event = crustalMap.get(eventID); - catalog = crustalCatalog; - loader = crustalLoader; - break; - case SUBDUCTION: - event = subductionMap.get(eventID); - catalog = subductionCatalog; - loader = subductionLoader; - break; - - default: - throw new IllegalStateException(); - } - typeEventMap.put(type, event); - typeEqkRupMap.put(type, catalog.getEqkRupture(event)); - typeRDMap.put(type, loader.readRotD50(site, eventID)); - - } - - List line = new ArrayList<>(header.size()); - line.add(site.getName()); - line.add(eventID+""); - for (Type type : types) - line.add((float)moment(typeEventMap.get(type))+""); - for (Type type : types) - line.add((float)typeEventMap.get(type).getMagnitude()+""); - for (Type type : types) - line.add((float)typeEqkRupMap.get(type).getRuptureSurface().getDistanceRup(site.getLoc())+""); - for (Type type : types) - line.add((float)typeEqkRupMap.get(type).getRuptureSurface().getDistanceJB(site.getLoc())+""); - for (double period : periods) - for (Type type : types) - line.add(typeRDMap.get(type).getInterpolatedY(period)+""); - - csv.addLine(line); - combCSV.addLine(line); - - for (int p=0; p 0, "Bad fullGM=%s for event %s, p=%s, site %s", - fullGM, eventID, periods[p], site.getName()); - Preconditions.checkState(crustalGM > 0, "Bad crustalGM=%s for event %s, p=%s, site %s", - crustalGM, eventID, periods[p], site.getName()); - Preconditions.checkState(subductionGM > 0, "Bad subductionGM=%s for event %s, p=%s, site %s", - subductionGM, eventID, periods[p], site.getName()); - - double fract = crustalGM > subductionGM ? subductionGM/crustalGM : crustalGM/subductionGM; - - for (int i=0; i= minFracts[i]) { - for (int q=0; q 0) { - File csvFile = new File(resourcesDir, site.getName()+".csv"); - csv.writeToFile(csvFile); - table.addLine(site.getName(), (float)site.getLoc().lat, (float)site.getLoc().lon, matches+" events", - "["+csvFile.getName()+"]("+resourcesDir.getName()+"/"+csvFile.getName()+")"); - } - } - - System.out.println("Found "+(combCSV.getNumRows()-1)+" records for "+uniqueEvents.size()+" events"); - File csvFile = new File(resourcesDir, "all_sites.csv"); - table.addLine("__All Sites__", "_N/A_", "_N/A_", (combCSV.getNumRows()-1)+" records for "+uniqueEvents.size()+" events", - "["+csvFile.getName()+"]("+resourcesDir.getName()+"/"+csvFile.getName()+")"); - combCSV.writeToFile(csvFile); - - lines.addAll(table.build()); - lines.add(""); - - for (int i=0; i funcs = new ArrayList<>(); - List chars = new ArrayList<>(); - - DefaultXY_DataSet oneToOne = new DefaultXY_DataSet(); - oneToOne.set(range.getLowerBound(), range.getLowerBound()); - oneToOne.set(range.getUpperBound(), range.getUpperBound()); - - funcs.add(oneToOne); - chars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 2f, Color.GRAY)); - - funcs.add(scatter); - chars.add(new PlotCurveCharacterstics(PlotSymbol.X, 3f, Color.BLACK)); - - PlotSpec spec = new PlotSpec(funcs, chars, title, "Full Ground Motion (g)", yAxisLabel+" (g)"); - - HeadlessGraphPanel gp = PlotUtils.initHeadless(); - - gp.drawGraphPanel(spec, log, log, range, range); - - PlotUtils.writePlots(outputDir, prefix+(log ? "_log" : ""), gp, 800, false, true, false, false); - } - } - -} diff --git a/src/main/java/scratch/kevin/simulators/ruptures/subduction/CrustalAndSubductionPageGen.java b/src/main/java/scratch/kevin/simulators/ruptures/subduction/CrustalAndSubductionPageGen.java new file mode 100644 index 00000000..54bcf1b6 --- /dev/null +++ b/src/main/java/scratch/kevin/simulators/ruptures/subduction/CrustalAndSubductionPageGen.java @@ -0,0 +1,963 @@ +package scratch.kevin.simulators.ruptures.subduction; + +import java.awt.Color; +import java.awt.Font; +import java.awt.geom.Point2D; +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.EnumMap; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map; +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.function.Supplier; + +import org.apache.commons.math3.stat.StatUtils; +import org.jfree.chart.annotations.XYAnnotation; +import org.jfree.chart.annotations.XYTextAnnotation; +import org.jfree.chart.ui.TextAnchor; +import org.jfree.data.Range; +import org.opensha.commons.data.CSVFile; +import org.opensha.commons.data.Site; +import org.opensha.commons.data.function.DefaultXY_DataSet; +import org.opensha.commons.data.function.DiscretizedFunc; +import org.opensha.commons.data.function.HistogramFunction; +import org.opensha.commons.data.function.LightFixedXFunc; +import org.opensha.commons.data.function.XY_DataSet; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.Region; +import org.opensha.commons.gui.plot.GeographicMapMaker; +import org.opensha.commons.gui.plot.HeadlessGraphPanel; +import org.opensha.commons.gui.plot.PlotCurveCharacterstics; +import org.opensha.commons.gui.plot.PlotLineType; +import org.opensha.commons.gui.plot.PlotSpec; +import org.opensha.commons.gui.plot.PlotSymbol; +import org.opensha.commons.gui.plot.PlotUtils; +import org.opensha.commons.mapping.PoliticalBoundariesData; +import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; +import org.opensha.commons.param.Parameter; +import org.opensha.commons.util.ExceptionUtils; +import org.opensha.commons.util.MarkdownUtils; +import org.opensha.commons.util.DataUtils.MinMaxAveTracker; +import org.opensha.commons.util.MarkdownUtils.TableBuilder; +import org.opensha.commons.util.cpt.CPT; +import org.opensha.sha.earthquake.EqkRupture; +import org.opensha.sha.imr.AttenRelRef; +import org.opensha.sha.imr.AttenRelSupplier; +import org.opensha.sha.imr.ScalarIMR; +import org.opensha.sha.imr.attenRelImpl.nshmp.NSHMP_AttenRelSupplier; +import org.opensha.sha.imr.param.IntensityMeasureParams.SA_Param; +import org.opensha.sha.imr.param.PropagationEffectParams.DistanceRupParameter; +import org.opensha.sha.simulators.EventRecord; +import org.opensha.sha.simulators.RSQSimEvent; +import org.opensha.sha.simulators.SimulatorElement; +import org.opensha.sha.simulators.utils.RSQSimEqkRupture; + +import com.google.common.base.Preconditions; +import com.google.common.collect.HashBasedTable; +import com.google.common.collect.HashBiMap; +import com.google.common.collect.Table; +import com.google.common.collect.Table.Cell; +import com.google.common.primitives.Doubles; + +import gov.usgs.earthquake.nshmp.gmm.Gmm; +import scratch.kevin.bbp.BBP_Module.VelocityModel; +import scratch.kevin.bbp.BBP_Site; +import scratch.kevin.simCompare.MultiRupGMPE_ComparePageGen; +import scratch.kevin.simulators.RSQSimCatalog; +import scratch.kevin.simulators.RSQSimCatalog.Catalogs; +import scratch.kevin.simulators.ruptures.BBP_CatalogSimZipLoader; +import scratch.kevin.simulators.ruptures.RSQSimBBP_Config; +import scratch.kevin.simulators.ruptures.RSQSimGeographicMapMaker; + +public class CrustalAndSubductionPageGen { + + private enum Type { + FULL("Full"), + CRUSTAL("Crustal"), + SUBDUCTION("Subduction"); + + private String prefix; + + private Type(String prefix) { + this.prefix = prefix; + } + }; + + private enum Quantity { + MAX("Maximum", "Max(Subduction, Crustal)") { + + @Override + public double calc(double subductionGM, double crustalGM) { + return Double.max(subductionGM, crustalGM); + } + + }, + CRUSTAL("Crustal-Only GM", "Crustal-Only") { + + @Override + public double calc(double subductionGM, double crustalGM) { + return crustalGM; + } + + }, + SUBDUCTION("Subduction-Only GM", "Subduction-Only") { + + @Override + public double calc(double subductionGM, double crustalGM) { + return subductionGM; + } + + }, + SUM_SQ("Sum Of Squares", "sqrt(Subduction^2 + Crustal^2)") { + + @Override + public double calc(double subductionGM, double crustalGM) { + return Math.sqrt(subductionGM*subductionGM + crustalGM*crustalGM); + } + + }, + SUM("Sum", "Subduction + Crustal") { + + @Override + public double calc(double subductionGM, double crustalGM) { + return subductionGM + crustalGM; + } + + }; + + private String name; + private String axisLabel; + + private Quantity(String name, String axisLabel) { + this.name = name; + this.axisLabel = axisLabel; + } + + public abstract double calc(double subductionGM, double crustalGM); + } + + private static final String DEBUG_SITE_NAME = "WLG"; + private static final int DEBUG_EVENT_ID = 152653; + + public static void main(String[] args) throws IOException { + RSQSimCatalog fullCatalog = Catalogs.BRUCE_5566.instance(); + RSQSimCatalog crustalCatalog = Catalogs.BRUCE_5566_CRUSTAL.instance(); + RSQSimCatalog subductionCatalog = Catalogs.BRUCE_5566_SUB.instance(); + + File bbpBaseDir = new File("/data/kevin/bbp/parallel"); + + File fullBBPdir = new File(bbpBaseDir, + "2023_03_30-rundir5566-all-m6.5-skipYears5000-noHF-vmLA_BASIN_500-standardSitesNZ-griddedSitesNZ"); + File crustalBBPdir = new File(bbpBaseDir, +// "2023_06_27-rundir5566_crustal-all-m6.5-skipYears5000-noHF-vmLA_BASIN_500-standardSitesNZ-griddedSitesNZ"); + "2023_06_27-rundir5566_crustal_corupture-all-m5.0-skipYears5000-noHF-vmLA_BASIN_500-standardSitesNZ-griddedSitesNZ"); + File subductionBBPdir = new File(bbpBaseDir, +// "2023_06_27-rundir5566_subduction-all-m6.5-skipYears5000-maxDist500-noHF-vmLA_BASIN_500-standardSitesNZ-griddedSitesNZ"); + "2023_06_27-rundir5566_subduction_corupture-all-m5.0-skipYears5000-noHF-vmLA_BASIN_500-standardSitesNZ-griddedSitesNZ"); + +// AttenRelSupplier subductionGMM = AttenRelRef.ZHAO_2006; + AttenRelSupplier subductionGMM = new NSHMP_AttenRelSupplier( + Gmm.AG_20_GLOBAL_INTERFACE, "AG2020 (Global)"); + AttenRelSupplier crustalGMM = AttenRelRef.ASK_2014; + + Quantity[] mapQuantities = {Quantity.SUM_SQ}; + double mapPeriod = 10d; + int maxMaps = 20; + + HashSet highlightSiteNames = new HashSet<>(); + + highlightSiteNames.add("WLG"); + + File markdownDir = new File("/home/kevin/markdown/rsqsim-analysis/catalogs/"); + + File markdownCatalogDir = new File(markdownDir, fullCatalog.getCatalogDir().getName()); + Preconditions.checkState(markdownCatalogDir.exists() || markdownCatalogDir.mkdir()); + File outputDir = new File(markdownCatalogDir, "crustal_subduction_gms"); + Preconditions.checkState(outputDir.exists() || outputDir.mkdir()); + File resourcesDir = new File(outputDir, "resources"); + Preconditions.checkState(resourcesDir.exists() || resourcesDir.mkdir()); + + VelocityModel vm = VelocityModel.LA_BASIN_500; + double minMag = 6.5; + + File fullZip = new File(fullBBPdir, "results_rotD.zip"); + File curstalZip = new File(crustalBBPdir, "results_rotD.zip"); + File subductionZip = new File(subductionBBPdir, "results_rotD.zip"); + + List sites = BBP_Site.readFile(crustalBBPdir); + + HashBiMap sitesBBPtoGMPE = HashBiMap.create(); + List gmpeSites = new ArrayList<>(); + for (BBP_Site site : sites) { + Site gmpeSite = site.buildGMPE_Site(vm); + gmpeSite.setName(RSQSimBBP_Config.siteCleanName(site)); + gmpeSites.add(gmpeSite); + sitesBBPtoGMPE.put(site, gmpeSite); + } + + List events = fullCatalog.loader().minMag(minMag).load(); + + Map crustalMap = crustalCatalog.loader().minMag(minMag).loadMap(); + Map subductionMap = subductionCatalog.loader().minMag(minMag).loadMap(); + + Map eventsMap = new HashMap<>(); + for (RSQSimEvent event : events) + if (crustalMap.containsKey(event.getID()) && subductionMap.containsKey(event.getID())) + eventsMap.put(event.getID(), event); + + BBP_CatalogSimZipLoader fullLoader = new BBP_CatalogSimZipLoader(fullZip, sites, sitesBBPtoGMPE, eventsMap); + BBP_CatalogSimZipLoader crustalLoader = new BBP_CatalogSimZipLoader(curstalZip, sites, sitesBBPtoGMPE, crustalMap); + BBP_CatalogSimZipLoader subductionLoader = new BBP_CatalogSimZipLoader(subductionZip, sites, sitesBBPtoGMPE, subductionMap); + + double[] periods = { 2d, 3d, 5d, 10d }; +// double[] periods = { 2d, 3d, 5d }; + + ExecutorService exec = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()); + + MultiRupGMPE_ComparePageGen.checkTransSiteParams(crustalGMM, gmpeSites); + MultiRupGMPE_ComparePageGen.checkTransSiteParams(subductionGMM, gmpeSites); + System.out.println("Calculating GMMs:"); + System.out.println("Full rup, crustal GMM"); + Table fullRupCrustalGMM = calcGMM( + fullCatalog, gmpeSites, fullLoader, periods, crustalGMM, exec); + System.out.println("Full rup, subduction GMM"); + Table fullRupSubductionGMM = calcGMM( + fullCatalog, gmpeSites, fullLoader, periods, subductionGMM, exec); + System.out.println("Partial rup, crustal GMM"); + Table partialRupCrustalGMM = calcGMM( + crustalCatalog, gmpeSites, crustalLoader, periods, crustalGMM, exec); + System.out.println("Partial rup, subduction GMM"); + Table partialRupSubductionGMM = calcGMM( + subductionCatalog, gmpeSites, subductionLoader, periods, subductionGMM, exec); + System.out.println("Done calculating GMMs"); + + double[] minFracts = { 0d, 0.2, 0.5 }; + + Quantity[] quantities = Quantity.values(); + + // now build combined gmms + List> combinedGMMs = new ArrayList<>(); + List> combinedSimFuncs = new ArrayList<>(); + for (Quantity q : quantities) { + // do from GMMs now + Table table = HashBasedTable.create(); + for (int eventID : partialRupCrustalGMM.rowKeySet()) { + for (Site site : partialRupCrustalGMM.row(eventID).keySet()) { + DiscretizedFunc crustalGMs = partialRupCrustalGMM.get(eventID, site); + DiscretizedFunc subductionGMs = partialRupSubductionGMM.get(eventID, site); + if (crustalGMs != null && subductionGMs != null) { + double[] yVals = new double[periods.length]; + for (int p=0; p eventSitePeriodFracts = HashBasedTable.create(); + + for (int i=0; i combCSV = null; + + HashSet uniqueEvents = new HashSet<>(); + + List lines = new ArrayList<>(); + + lines.add("# Crustal and Subduction Co-Rupture Ground Motions"); + lines.add(""); + + int tocIndex = lines.size(); + String topLink = "*[(top)](#table-of-contents)*"; + + lines.add("## Sites and CSV Files"); + lines.add(topLink); lines.add(""); + + TableBuilder table = MarkdownUtils.tableBuilder(); + + table.addLine("Site Name", "Latitude", "Longitude", "Record Count", "CSV File"); + + List highlightSites = new ArrayList<>(); + + for (BBP_Site site : sites) { + CSVFile csv = new CSVFile<>(true); + + if (highlightSiteNames.contains(site.getName())) + highlightSites.add(site); + + List header = new ArrayList<>(); + header.add("Site Name"); + header.add("Event ID"); + + for (Type type : types) + header.add(type.prefix+" Moment"); + for (Type type : types) + header.add(type.prefix+" Magnitude"); + for (Type type : types) + header.add(type.prefix+" rRup"); + for (Type type : types) + header.add(type.prefix+" rJB"); + for (double period : periods) + for (Type type : types) + header.add(type.prefix+" "+oDF.format(period)+"s RotD50"); + + csv.addLine(header); + + if (combCSV == null) { + combCSV = new CSVFile<>(true); + combCSV.addLine(header); + } + + int matches = 0; + + Site gmmSite = sitesBBPtoGMPE.get(site); + + for (RSQSimEvent fullEvent : events) { + int eventID = fullEvent.getID(); + RSQSimEvent crustalEvent = crustalMap.get(eventID); + RSQSimEvent subductionEvent = subductionMap.get(eventID); + + if (crustalEvent != null && crustalLoader.contains(site, eventID) + && subductionEvent != null && subductionLoader.contains(site, eventID)) { + // exists for both + matches++; + uniqueEvents.add(eventID); + + Map typeEventMap = new EnumMap<>(Type.class); + Map typeEqkRupMap = new EnumMap<>(Type.class); + Map typeRDMap = new EnumMap<>(Type.class); + for (Type type : types) { + RSQSimEvent event; + RSQSimCatalog catalog; + BBP_CatalogSimZipLoader loader; + switch (type) { + case FULL: + event = fullEvent; + catalog = fullCatalog; + loader = fullLoader; + break; + case CRUSTAL: + event = crustalMap.get(eventID); + catalog = crustalCatalog; + loader = crustalLoader; + break; + case SUBDUCTION: + event = subductionMap.get(eventID); + catalog = subductionCatalog; + loader = subductionLoader; + break; + + default: + throw new IllegalStateException(); + } + typeEventMap.put(type, event); + typeEqkRupMap.put(type, catalog.getEqkRupture(event)); + typeRDMap.put(type, loader.readRotD50(site, eventID)); + + } + + List line = new ArrayList<>(header.size()); + line.add(site.getName()); + line.add(eventID+""); + for (Type type : types) + line.add((float)moment(typeEventMap.get(type))+""); + for (Type type : types) + line.add((float)typeEventMap.get(type).getMagnitude()+""); + for (Type type : types) + line.add((float)typeEqkRupMap.get(type).getRuptureSurface().getDistanceRup(site.getLoc())+""); + for (Type type : types) + line.add((float)typeEqkRupMap.get(type).getRuptureSurface().getDistanceJB(site.getLoc())+""); + for (double period : periods) + for (Type type : types) + line.add(typeRDMap.get(type).getInterpolatedY(period)+""); + + csv.addLine(line); + combCSV.addLine(line); + + double[] siteEventFracts = new double[periods.length]; + + for (int p=0; p 0, "Bad fullGM=%s for event %s, p=%s, site %s", + fullGM, eventID, periods[p], site.getName()); + Preconditions.checkState(crustalGM > 0, "Bad crustalGM=%s for event %s, p=%s, site %s", + crustalGM, eventID, periods[p], site.getName()); + Preconditions.checkState(subductionGM > 0, "Bad subductionGM=%s for event %s, p=%s, site %s", + subductionGM, eventID, periods[p], site.getName()); + + double fract = crustalGM > subductionGM ? subductionGM/crustalGM : crustalGM/subductionGM; + siteEventFracts[p] = fract; + + for (int q=0; q= minFracts[i]) + scatters[i][q][p].set(fullGM, combVal); + } + } + + eventSitePeriodFracts.put(eventID, gmmSite, siteEventFracts); + } + } + System.out.println("Found "+matches+" matches for "+site.getName()); + if (matches > 0) { + File csvFile = new File(resourcesDir, site.getName()+".csv"); + csv.writeToFile(csvFile); + table.addLine(site.getName(), (float)site.getLoc().lat, (float)site.getLoc().lon, matches+" events", + "["+csvFile.getName()+"]("+resourcesDir.getName()+"/"+csvFile.getName()+")"); + } + } + + System.out.println("Found "+(combCSV.getNumRows()-1)+" records for "+uniqueEvents.size()+" events"); + File csvFile = new File(resourcesDir, "all_sites.csv"); + table.addLine("__All Sites__", "_N/A_", "_N/A_", (combCSV.getNumRows()-1)+" records for "+uniqueEvents.size()+" events", + "["+csvFile.getName()+"]("+resourcesDir.getName()+"/"+csvFile.getName()+")"); + combCSV.writeToFile(csvFile); + + lines.addAll(table.build()); + lines.add(""); + + String crustalShortName = crustalGMM.getShortName(); + String subductionShortName = subductionGMM.getShortName(); + + for (int i=0; i combinedRupComps = combinedGMMs.get(q); + Table combinedSimRupComps = combinedSimFuncs.get(q); + + for (int p=-1; p fullRupComps; + Table partialRupComps; + BBP_CatalogSimZipLoader partialSimLoader; + String partialName; + Color color; + String gmmName; + if (useSubduction) { + lines.add("### Subduction GMM Residuals, "+fractStr); + lines.add(""); + gmPrefix = fractPrefix+"_subduction_gmm"; + partialName = "Subduction Rupture Portion"; + lines.add("Residuals between simulated ground motions and the subduction GMM, "+subductionShortName); + fullRupComps = fullRupSubductionGMM; + partialRupComps = partialRupSubductionGMM; + color = Color.BLUE.darker(); + gmmName = subductionShortName; + partialSimLoader = subductionLoader; + } else { + lines.add("### Crustal GMM Residuals, "+fractStr); + lines.add(""); + gmPrefix = fractPrefix+"_crustal_gmm"; + partialName = "Crustal Rupture Portion"; + lines.add("Residuals between simulated ground motions and the crustal GMM, "+crustalShortName); + fullRupComps = fullRupCrustalGMM; + partialRupComps = partialRupCrustalGMM; + color = Color.RED.darker(); + gmmName = crustalShortName; + partialSimLoader = crustalLoader; + } + lines.add(""); + + table = MarkdownUtils.tableBuilder(); + table.addLine("Period", partialName, "Full Rupture"); + + for (int p=-1; p 0) { + lines.add("## Event Map GMM Residual"); + lines.add(topLink); lines.add(""); + + lines.add("Debug residual maps of simulations vs each GMM, all for "+oDF.format(mapPeriod)+"s"); + lines.add(""); + + table = MarkdownUtils.tableBuilder(); + table.initNewLine().addColumns(crustalShortName+" Residuals", subductionShortName+" Residuals"); + for (Quantity q : mapQuantities) + table.addColumn("Combined GMM Residuals, "+q.name); + table.finalizeLine(); + + List mapEvents = new ArrayList<>(); + for (RSQSimEvent event : eventsMap.values()) { + if (fullRupCrustalGMM.containsRow(event.getID())) + mapEvents.add(event); + } + Collections.sort(mapEvents); + if (mapEvents.size() > maxMaps) { + System.out.println("Will plot "+maxMaps+"/"+mapEvents.size()+" events"); + mapEvents = mapEvents.subList(0, maxMaps); + } else { + System.out.println("Will plot all "+mapEvents.size()+" events"); + } + for (RSQSimEvent event : mapEvents) { + System.out.println("Plotting maps for event "+event.getID()); + table.initNewLine(); + + String prefix = "event_"+event.getID(); + plotResidualMap(resourcesDir, prefix+"_crustal", event, fullLoader, fullRupCrustalGMM.row(event.getID()), + mapPeriod, sitesBBPtoGMPE.inverse(), "Log Residuals: "+crustalShortName+" - Full Simulated"); + table.addColumn("![Crustal Map]("+resourcesDir.getName()+"/"+prefix+"_crustal.png)"); + plotResidualMap(resourcesDir, prefix+"_subduction", event, fullLoader, fullRupSubductionGMM.row(event.getID()), + mapPeriod, sitesBBPtoGMPE.inverse(), "Log Residuals: "+subductionShortName+" - Full Simulated"); + table.addColumn("![Subduction Map]("+resourcesDir.getName()+"/"+prefix+"_subduction.png)"); + for (Quantity q : mapQuantities) { + int qIndex = -1; + for (int i=0; i= 0); + String qPrefix = prefix+"_"+q.name(); + plotResidualMap(resourcesDir, qPrefix, event, fullLoader, combinedGMMs.get(qIndex).row(event.getID()), + mapPeriod, sitesBBPtoGMPE.inverse(), "Log Residuals: Combined GMM ("+q.name+") - Full Simulated"); + table.addColumn("!["+q.name+" Map]("+resourcesDir.getName()+"/"+qPrefix+".png)"); + } + + table.finalizeLine(); + } + + lines.addAll(table.build()); + lines.add(""); + } + + // add TOC + lines.addAll(tocIndex, MarkdownUtils.buildTOC(lines, 2)); + lines.add(tocIndex, "## Table Of Contents"); + + // write markdown + MarkdownUtils.writeReadmeAndHTML(lines, outputDir); + + exec.shutdown(); + } + + private static DefaultXY_DataSet union(DefaultXY_DataSet[] datas) { + DefaultXY_DataSet ret = new DefaultXY_DataSet(); + for (DefaultXY_DataSet data : datas) + for (Point2D pt : data) + ret.set(pt); + return ret; + } + + private static final DecimalFormat oDF = new DecimalFormat("0.##"); + + private static double moment(RSQSimEvent event) { + double ret = 0d; + for (EventRecord rec : event) + ret += rec.getMoment(); + return ret; + } + + private static void writeScatter(File outputDir, String prefix, double period, DefaultXY_DataSet scatter, + String yAxisLabel) throws IOException { + double min = Math.min(scatter.getMinX(), scatter.getMinY()); + double max = Math.max(scatter.getMaxX(), scatter.getMaxY()); + + String title = Double.isFinite(period) ? oDF.format(period)+"s SA" : "All Periods"; + + for (boolean log : new boolean[] {false, true}) { + Range range; + if (log) + range = new Range(Math.pow(10, Math.floor(Math.log10(min))), + Math.pow(10, Math.ceil(Math.log10(max)))); + else + range = new Range(0d, max*1.05); + + List funcs = new ArrayList<>(); + List chars = new ArrayList<>(); + + DefaultXY_DataSet oneToOne = new DefaultXY_DataSet(); + oneToOne.set(range.getLowerBound(), range.getLowerBound()); + oneToOne.set(range.getUpperBound(), range.getUpperBound()); + + funcs.add(oneToOne); + chars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 2f, Color.GRAY)); + + funcs.add(scatter); + chars.add(new PlotCurveCharacterstics(PlotSymbol.X, 3f, Color.BLACK)); + + PlotSpec spec = new PlotSpec(funcs, chars, title, "Full Ground Motion (g)", yAxisLabel+" (g)"); + + HeadlessGraphPanel gp = PlotUtils.initHeadless(); + + gp.drawGraphPanel(spec, log, log, range, range); + + PlotUtils.writePlots(outputDir, prefix+(log ? "_log" : ""), gp, 800, false, true, false, false); + } + } + + private static Table calcGMM(RSQSimCatalog catalog, + List sites, BBP_CatalogSimZipLoader simLoader, double[] periods, Supplier gmmRef, + ExecutorService exec) { + ArrayDeque gmmDeque = new ArrayDeque<>(); + + Table> futures = HashBasedTable.create(); + + for (Site site : sites) { + for (RSQSimEvent event : simLoader.getRupturesForSite(site)) { + GMM_CalcCallable call = new GMM_CalcCallable(gmmDeque, gmmRef, site, + event.getID(), catalog.getEqkRupture(event), periods); + futures.put(event.getID(), site, exec.submit(call)); + } + } + + Table ret = HashBasedTable.create(); + try { + for (Cell> cell : futures.cellSet()) { + DiscretizedFunc gms = futures.get(cell.getRowKey(), cell.getColumnKey()).get(); + ret.put(cell.getRowKey(), cell.getColumnKey(), gms); + } + } catch (InterruptedException | ExecutionException e) { + throw ExceptionUtils.asRuntimeException(e); + } + System.out.println("Calculated GMs for "+ret.rowKeySet().size()+" events, "+ret.size()+" event/site pairs"); + + return ret; + } + + private static class GMM_CalcCallable implements Callable { + + private ArrayDeque gmmDeque; + private Supplier gmmRef; + private Site site; + private int eventID; + private EqkRupture rup; + private double[] periods; + + public GMM_CalcCallable(ArrayDeque gmmDeque, Supplier gmmRef, Site site, int eventID, EqkRupture rup, double[] periods) { + this.gmmDeque = gmmDeque; + this.gmmRef = gmmRef; + this.site = site; + this.eventID = eventID; + this.rup = rup; + this.periods = periods; + } + + @Override + public DiscretizedFunc call() throws Exception { + ScalarIMR gmm = null; + synchronized (gmmDeque) { + if (!gmmDeque.isEmpty()) + gmm = gmmDeque.pop(); + } + if (gmm == null) { + gmm = gmmRef.get(); + gmm.setParamDefaults(); + gmm.setIntensityMeasure(SA_Param.NAME); + } + + Parameter im = gmm.getIntensityMeasure(); + gmm.setAll(rup, site, im); + double[] vals = new double[periods.length]; + boolean D = eventID == DEBUG_EVENT_ID && site.getName().equals(DEBUG_SITE_NAME); + for (int p=0; p eventSitePeriodFracts, + Table gmmResults, Map sitesGMPEtoBBP) throws IOException { + List residuals = new ArrayList<>(); + for (Cell cell : eventSitePeriodFracts.cellSet()) { + int eventID = cell.getRowKey(); + Site site = cell.getColumnKey(); + BBP_Site bppSite = sitesGMPEtoBBP.get(site); + double[] fracts = cell.getValue(); + Preconditions.checkState(fracts.length == allPeriods.length); + boolean D = limitPeriod == null && eventID == DEBUG_EVENT_ID && site.getName().equals(DEBUG_SITE_NAME); + if (D) System.out.println("Residual debug for "+eventID+", "+site.getName()+", "+simName+" vs "+modelName+";\tTitle: "+title); + for (int p=0; p= minFract && (limitPeriod == null || limitPeriod.doubleValue() == allPeriods[p])) { + DiscretizedFunc gmmFunc = gmmResults.get(eventID, site); + Preconditions.checkNotNull(gmmFunc); + DiscretizedFunc simFunc = simLoader.readRotD50(bppSite, eventID); + Preconditions.checkNotNull(simFunc); + double logSim = Math.log(simFunc.getY(allPeriods[p])); + double logGMM = Math.log(gmmFunc.getY(allPeriods[p])); + double residual = logGMM - logSim; + if (D) + System.out.println("\tT="+(float)allPeriods[p]+": " + +(float)logGMM+" - "+(float)logSim+" = "+(float)residual); + residuals.add(residual); + } + } + } + + double[] residualArray = Doubles.toArray(residuals); + double mean = StatUtils.mean(residualArray); + double stdDev = Math.sqrt(StatUtils.variance(residualArray)); + +// Range xRange = new Range(-3, 3); + Range xRange = new Range(-2, 2); + HistogramFunction hist = HistogramFunction.getEncompassingHistogram( + xRange.getLowerBound()+0.01, xRange.getUpperBound()-0.01, 0.1); + for (double val : residuals) + hist.add(hist.getClosestXIndex(val), 1d); + + Range yRange = new Range(0d, hist.getMaxY()*1.1); + + List funcs = new ArrayList<>(); + List chars = new ArrayList<>(); + + DefaultXY_DataSet meanLine = new DefaultXY_DataSet(); + meanLine.set(mean, yRange.getLowerBound()); + meanLine.set(mean, yRange.getUpperBound()); + + funcs.add(hist); + chars.add(new PlotCurveCharacterstics(PlotLineType.HISTOGRAM, 1f, color)); + + funcs.add(meanLine); + chars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 4f, Color.BLACK)); + + PlotSpec spec = new PlotSpec(funcs, chars, title, "Log Residual: "+modelName+" - "+simName, "Count"); + + DecimalFormat df = new DecimalFormat("0.00"); + Font font = new Font(Font.SANS_SERIF, Font.BOLD, 24); + XYTextAnnotation meanAnn = new XYTextAnnotation("Mean: "+df.format(mean), + xRange.getUpperBound()-0.1, yRange.getUpperBound()*0.975); + meanAnn.setFont(font); + meanAnn.setTextAnchor(TextAnchor.TOP_RIGHT); + spec.addPlotAnnotation(meanAnn); + XYTextAnnotation sdAnn = new XYTextAnnotation("Std. Dev: "+df.format(stdDev), + xRange.getUpperBound()-0.1, yRange.getUpperBound()*0.925); + sdAnn.setFont(font); + sdAnn.setTextAnchor(TextAnchor.TOP_RIGHT); + spec.addPlotAnnotation(sdAnn); + + HeadlessGraphPanel gp = PlotUtils.initHeadless(); + + gp.drawGraphPanel(spec, false, false, xRange, yRange); + + PlotUtils.writePlots(outputDir, prefix, gp, 800, 700, true, true, false); + } + + private static void plotResidualMap(File outputDir, String prefix, RSQSimEvent event, + BBP_CatalogSimZipLoader simLoader, Map modelGMs, + double period, Map sitesGMPEtoBBP, String label) throws IOException { + Region region = new Region(new Location(-33.5, 165), new Location(-47, 179)); + RSQSimGeographicMapMaker mapMaker = new RSQSimGeographicMapMaker(region, PoliticalBoundariesData.loadNZOutlines()); + + mapMaker.setWriteGeoJSON(false); + mapMaker.setWritePDFs(false); + + MinMaxAveTracker latTrack = new MinMaxAveTracker(); + MinMaxAveTracker lonTrack = new MinMaxAveTracker(); + MinMaxAveTracker depTrack = new MinMaxAveTracker(); + List depths = new ArrayList<>(); + for (SimulatorElement elem : event.getAllElements()) { + for (Location loc : elem.getVertices()) { + latTrack.addValue(loc.getLatitude()); + lonTrack.addValue(loc.getLongitude()); + depTrack.addValue(loc.getDepth()); + } + depths.add(elem.getAveDepth()); + } +// System.out.println("Element locations:"); +// System.out.println("\tLatitude:\t"+latTrack); +// System.out.println("\tLongitude:\t"+lonTrack); +// System.out.println("\tDepth:\t"+depTrack); + +// mapMaker.plotEvent(event, Color.LIGHT_GRAY, null, Float.NaN); + CPT depthCPT = new CPT(0d, depTrack.getMax(), Color.BLACK, Color.LIGHT_GRAY); + mapMaker.plotEventFillScalars(event, depths, depthCPT, "Depth (km)"); + mapMaker.plotEventHypocenter(Color.GREEN.darker()); + + CPT residCPT = GMT_CPT_Files.DIVERGING_VIK_UNIFORM.instance().rescale(-2d, 2d); + + List siteLocs = new ArrayList<>(); + List siteResiduals = new ArrayList<>(); + System.out.println("Have "+modelGMs.size()+" GMM sites"); + for (Site site : modelGMs.keySet()) { + BBP_Site bbpSite = sitesGMPEtoBBP.get(site); + + if (simLoader.contains(bbpSite, event.getID())) { + double simVal = Math.log(simLoader.readRotD50(bbpSite, event.getID()).getY(period)); + double gmVal = Math.log(modelGMs.get(site).getY(period)); + double residual = gmVal - simVal; + siteLocs.add(site.getLocation()); + siteResiduals.add(residual); + } + } + System.out.println("Plotting "+siteResiduals.size()+" site residuals"); + + mapMaker.plotScatterScalars(siteLocs, siteResiduals, residCPT, label); + mapMaker.setScatterSymbol(PlotSymbol.FILLED_TRIANGLE, 8f, PlotSymbol.TRIANGLE, Color.BLACK); + +// mapMaker.plotScatters(siteLocs, null); + + mapMaker.plot(outputDir, prefix, "Event "+event.getID()+", M"+oDF.format(event.getMagnitude())+", "+oDF.format(period)+"s SA"); + } + +} diff --git a/src/main/java/scratch/kevin/ucerf3/PureScratch.java b/src/main/java/scratch/kevin/ucerf3/PureScratch.java index f94ff8ab..76b04f8b 100644 --- a/src/main/java/scratch/kevin/ucerf3/PureScratch.java +++ b/src/main/java/scratch/kevin/ucerf3/PureScratch.java @@ -98,6 +98,7 @@ import org.opensha.commons.logicTree.LogicTreeLevel.RandomlySampledLevel; import org.opensha.commons.logicTree.LogicTreeNode.RandomlySampledNode; import org.opensha.commons.logicTree.LogicTreeNode; +import org.opensha.commons.mapping.PoliticalBoundariesData; import org.opensha.commons.mapping.gmt.elements.GMT_CPT_Files; import org.opensha.commons.util.DataUtils; import org.opensha.commons.util.DataUtils.MinMaxAveTracker; @@ -265,6 +266,7 @@ import scratch.kevin.simCompare.SiteHazardCurveComarePageGen; import scratch.kevin.simulators.RSQSimCatalog; import scratch.kevin.simulators.RSQSimCatalog.Catalogs; +import scratch.kevin.simulators.ruptures.RSQSimGeographicMapMaker; public class PureScratch { @@ -5143,6 +5145,43 @@ private static void test268() throws IOException { System.out.println("Converted GeoJSONFaultSection feature after reserialization:"); System.out.println(sect2.toFeature().toJSON()); System.out.println("Default surface type: "+sect2.getFaultSurface(1d).getClass().getName()); + + List plotSects = new ArrayList<>(); + FaultSection fullApproxSect = sect.clone(); + fullApproxSect.setSectionName("Approx Gridded Sect"); + fullApproxSect.setSectionId(plotSects.size()); + plotSects.add(fullApproxSect); + + Feature stirlingFeature = Feature.fromJSON(sect.toFeature().toJSON()); + stirlingFeature.properties.remove(GeoJSONFaultSection.LOWER_TRACE); + FaultSection stirlingSect = GeoJSONFaultSection.fromFeature(stirlingFeature); + stirlingSect.setSectionName("Stirling Sect"); + stirlingSect.setSectionId(plotSects.size()); + plotSects.add(stirlingSect); + + FaultSection fakeAseisApproxSect = sect.clone(); + fakeAseisApproxSect.setSectionName("Fake Aseis Approx Gridded Sect"); + fakeAseisApproxSect.setSectionId(plotSects.size()); + fakeAseisApproxSect.setAseismicSlipFactor(0.25); + plotSects.add(fakeAseisApproxSect); + + FaultSection fakeAseisStirlingSect = stirlingSect.clone(); + fakeAseisStirlingSect.setSectionName("Fake Aseis Stirling Sect"); + fakeAseisStirlingSect.setSectionId(plotSects.size()); + fakeAseisStirlingSect.setAseismicSlipFactor(0.25); + plotSects.add(fakeAseisStirlingSect); + + Region plotReg = GeographicMapMaker.buildBufferedRegion(plotSects); + GeographicMapMaker mapMaker = new GeographicMapMaker(plotReg); + mapMaker.setPlotAseisReducedSurfaces(true); + mapMaker.setWriteGeoJSON(true); + + for (int i=0; i fakeScalars = new ArrayList<>(); + for (int i=0; i treePropsList = feature.properties.getPropertiesList("mfd-tree"); + System.out.println("Loaded a list of "+treePropsList.size()+" tree properties:"); + for (FeatureProperties treeProps : treePropsList) { + String id = treeProps.getString("id"); + // need to supply a default value for primitives as null can't be returned if it's not found, thus the Double.NaN + double weight = treeProps.getDouble("weight", Double.NaN); + System.out.println("\tid="+id+" with weight="+weight); + // 'value' is a map, so we'll load it in as its own FeatureProperties + FeatureProperties value = treeProps.getProperties("value"); + String type = value.getString("type"); + double m = value.getDouble("m", Double.NaN); + System.out.println("\t\ttype: "+type+", m="+m); + } +// FeatureProperties treeProps = feature.properties.getProperties("mfd-tree", null); +// System.out.println(feature.toJSON()); + } + /** * @param args * @throws Exception */ public static void main(String[] args) throws Exception { - test268(); + test270(); } }