From 69d5ccf096b293d935a0dd9043ad8d186c7d89b6 Mon Sep 17 00:00:00 2001 From: field Date: Sun, 1 Sep 2024 10:29:52 -0600 Subject: [PATCH] Misc TD loss modeling progress --- .../ned/FSS_Inversion2019/PlottingUtils.java | 45 + .../U3ETAS_LossSimulationAnalysis.java | 2321 +++++++++++++++-- .../U3ETAS_SimulationAnalysis.java | 7 +- 3 files changed, 2116 insertions(+), 257 deletions(-) diff --git a/src/main/java/scratch/ned/FSS_Inversion2019/PlottingUtils.java b/src/main/java/scratch/ned/FSS_Inversion2019/PlottingUtils.java index 3b45f0b3..c06996e7 100644 --- a/src/main/java/scratch/ned/FSS_Inversion2019/PlottingUtils.java +++ b/src/main/java/scratch/ned/FSS_Inversion2019/PlottingUtils.java @@ -10,10 +10,12 @@ import org.jfree.chart.axis.NumberAxis; import org.jfree.data.Range; import org.opensha.commons.data.function.XY_DataSet; +import org.opensha.commons.data.uncertainty.UncertainArbDiscFunc; import org.opensha.commons.data.xyz.EvenlyDiscrXYZ_DataSet; import org.opensha.commons.gui.plot.GraphWindow; 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.PlotPreferences; import org.opensha.commons.gui.plot.PlotSpec; import org.opensha.commons.gui.plot.jfreechart.xyzPlot.XYZPlotSpec; @@ -60,6 +62,49 @@ public static void writeAndOrPlotFuncs( } + /** + * Plotting method for a list of UncertainArbDiscFunc with default width and height + * @param funcs + * @param plotChars + * @param plotName + * @param xAxisLabel + * @param yAxisLabel + * @param xAxisRange + * @param yAxisRange + * @param logX + * @param logY + * @param fileNamePrefix - set a null if you don't want to save to files + * @param popupWindow - set as false if you don't want a pop-up windows with the plots + */ + public static void writeAndOrPlotUncertFuncs( + ArrayList uncertFuncs, + ArrayList plotChars, + String plotName, + String xAxisLabel, + String yAxisLabel, + Range xAxisRange, + Range yAxisRange, + boolean logX, + boolean logY, + String fileNamePrefix, + boolean popupWindow) { + // make the douplicates for shaded regions + ArrayList funcs = new ArrayList(); + ArrayList plotChars2 = new ArrayList(); + for(int i=0; i meanLossForNthRup, lossCOVForNthRup, magForNthRup; + HashMap fssIndexForNthRup, gridCellIndexForNthRup; + FaultSystemSolutionERF_ETAS erf=null; + + static double totalPorfolioValue = 483e9/1e3; // total portfolio value in units of $1000 static double catalogDuration=500; // hard coded static double catalogStartYear=2012; static long millisPerYr = (long)(1000*60*60*24*365.25); - + long startTimeMillis, endTimeMillis, fullDurationMillis; ArrayList catalogList; ArrayList catalogDeclusteredList; @@ -114,37 +113,926 @@ public class U3ETAS_LossSimulationAnalysis { static double deltaMag = 0.1; static int numMag = 40; - static double lossCurveLnMin = Math.log(0.001); - static double lossCurveLnMax = Math.log(10); - static int lossCurveNum = 100; //40; - static double lossCurveDelta = (lossCurveLnMax-lossCurveLnMin)/(double)(lossCurveNum-1); + // following are units of 1000 dollars + static double lossCurveLogMin = 0; // 1,000 + static double lossCurveLogMax = 9; // 1 trillion + static int lossCurveNum = 91; //40; + static double lossCurveDelta = (lossCurveLogMax-lossCurveLogMin)/(double)(lossCurveNum-1); + ArbitrarilyDiscretizedFunc curveXvals = getBlankLossCurve(0d); + LossCOV_Model covModel = LossCOV_Model.PORTER_POWER_LAW_2020_09_01_fixed; + double[] randomLossForEventID; + int totalNumEvents; - public U3ETAS_LossSimulationAnalysis() { + public U3ETAS_LossSimulationAnalysis(int seed) { + + // Note the grid cell stated for a catalog rupture does not necessarily equal that + // implied by the hypocenter (known leakage problem due to randomness added in U3ETAS) + + this.randomGen = new JDKRandomGenerator(seed); + + + File outputDir = new File(dirName); + if(!outputDir.exists()) outputDir.mkdir(); + TimeSpan tsp = new TimeSpan(TimeSpan.YEARS,TimeSpan.YEARS); + tsp.setStartTime((int)catalogStartYear); + tsp.setDuration((int)catalogDuration); + startTimeMillis = tsp.getStartTimeInMillis(); + endTimeMillis = tsp.getEndTimeCalendar().getTimeInMillis(); + fullDurationMillis = endTimeMillis-startTimeMillis; + + erf = getERF(); + + // Load catalogs & loss data + CSVFile lossDataCSV_File=null; try { -// this.catalogList = loadCatalogs(new File(fssFileName), new File(catalogsFileName), 0.0); - this.catalogList = loadCatalogs(new File(fssFileName), new File(catalogsFileName), 5.0); + catalogList = loadCatalogs(new File(catalogsFileName), 5.0, false); + lossDataCSV_File = CSVFile.readFile(new File(lossDataCSV_Filename), true); } catch (IOException | DocumentException e) { throw ExceptionUtils.asRuntimeException(e); } + if (D) { + System.out.println("num Catalogs = "+catalogList.size()); + System.out.println("Loss Data:"); + System.out.println("\tNumCols: "+lossDataCSV_File.getNumCols()); + System.out.println("\tNumRows: "+lossDataCSV_File.getNumRows()); + + } - if (D) System.out.println("num Catalogs = "+catalogList.size()); + // set HashMaps from loss data + meanLossForNthRup = new HashMap(); + lossCOVForNthRup = new HashMap(); + magForNthRup = new HashMap(); + fssIndexForNthRup = new HashMap(); + gridCellIndexForNthRup = new HashMap(); + for(int r=1;rlossCOVForNthRup.get(n)) minCOVfile=lossCOVForNthRup.get(n); + if(maxCOVfilecovCalc) minCOVcalc=covCalc; + if(maxCOVcalc mfdList = new ArrayList(); + mfdList.add(catalogMFD); + mfdList.add(catalogGKdeclMFD); + mfdList.add(catalogGKfiltMFD); + mfdList.add(catalogSpontaneustMFD); + mfdList.add(erfMFD_TD); + mfdList.add(erfMFD_TimeInd); + ArrayList plotChars = new ArrayList(); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLUE)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.RED)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLACK)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.GREEN)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.ORANGE)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.GRAY)); + PlottingUtils.writeAndOrPlotFuncs(mfdList, plotChars, "MFDs", "Mag", "Rate (/yr)", + new Range(5,8.5), new Range(1e-5,2), false, true, null, true); + + ArrayList mfdCumList = new ArrayList(); + mfdCumList.add(catalogMFD.getCumRateDistWithOffset()); + mfdCumList.add(catalogGKdeclMFD.getCumRateDistWithOffset()); + mfdCumList.add(catalogGKfiltMFD.getCumRateDistWithOffset()); + mfdCumList.add(catalogSpontaneustMFD.getCumRateDistWithOffset()); + mfdCumList.add(erfMFD_TD.getCumRateDistWithOffset()); + mfdCumList.add(erfMFD_TimeInd.getCumRateDistWithOffset()); + PlottingUtils.writeAndOrPlotFuncs(mfdCumList, plotChars, "Cumunlative MFDs", "Mag", "Cumulative Rate (/yr)", + new Range(5,8.5), new Range(1e-5,10), false, true, null, true); + + IncrementalMagFreqDist mfdGKdeclMFD_Ratio = makeMFD(); + mfdGKdeclMFD_Ratio.setName("mfdGKdeclMFD_Ratio"); + IncrementalMagFreqDist mfdGKfiltMFD_Ratio = makeMFD(); + mfdGKfiltMFD_Ratio.setName("mfdGKfiltMFD_Ratio"); + IncrementalMagFreqDist mfdSpontaneousMFD_Ratio = makeMFD(); + mfdSpontaneousMFD_Ratio.setName("mfdSpontaneousMFD_Ratio"); + for(int i=0;i ratioList = new ArrayList(); + ratioList.add(mfdGKdeclMFD_Ratio); + ratioList.add(mfdGKfiltMFD_Ratio); + ratioList.add(mfdSpontaneousMFD_Ratio); + PlottingUtils.writeAndOrPlotFuncs(ratioList, plotChars, "Incr MFD Ratios", "Mag", "Ratio", + new Range(5,8.5), new Range(0,1), false, false, null, true); } + /** + * Units: Millions of dollars + */ + public void writeAAL_ValuesBillions() { + System.out.println("\nAverage Annual losses (AAL; billion $):"); + +// double aal_Catalog = getAveAnnaulLossFromCatalogList(catalogList)/1e6d; +// System.out.println("\tAAL from catalogs: "+(float)aal_Catalog); + double[] stats = getLossStatsFromCatalogList(catalogList, true); + System.out.println("\tAAL from catalogs: "+(float)(stats[0]/1e6d)+"\t(AAL=EventRate*MeanLoss)"); + System.out.println("\t\tEventRate = : "+(float)(stats[1])); + System.out.println("\t\tMeanLoss = : "+(float)(stats[2]/1e6d)); + System.out.println("\t\tMedianLoss = : "+(float)(stats[3]/1e6d)); + System.out.println("\t\tLossCOV = : "+(float)(stats[4])); + System.out.println("\t\tMinNonZeroLoss = : "+(float)(stats[5]/1e6d)); + System.out.println("\t\tMaxLoss = : "+(float)(stats[6]/1e6d)); + System.out.println("\t\tfractZeroLossRups = : "+(float)(stats[7])); + + + + double aal_TD = getAveAnnualLossFromERF(true)/1e6d; // Time Dependent + if(D) System.out.println("\tAAL from ERF-TD = "+(float)aal_TD); + + double aal_TimeInd = getAveAnnualLossFromERF(false)/1e6d; // Time Independent + if(D) System.out.println("\tAAL from ERF-TimeInd = "+(float)aal_TimeInd); + + double aal_GKdeclCatalog = getAveAnnaulLossFromCatalogList(catalogDeclusteredList)/1e6d; + System.out.println("\tAAL from GK declustered catalogs: "+(float)aal_GKdeclCatalog+"\n"); + + double aal_GKfilteredCatalog = getAveAnnaulLossFromCatalogList(getU3_GK_FilteredCatalog(catalogList))/1e6d; + System.out.println("\tAAL from GK U3 filtered catalogs: "+(float)aal_GKfilteredCatalog+"\n"); + + + } + + + /** + * This plots the distribution of rupture losses (mean and 68% and 95% fractiles) as a function + * of magnitude (not including rup rates). + * + * @param catList + * @param mkPlot - not yet working + */ + public void plotRupLossVsMagStats(ArrayList catList, boolean randomFromCOV, boolean mkPlot) { + + IncrementalMagFreqDist mfd = makeMFD(); + ArbDiscrEmpiricalDistFunc_3D arbFunc3D = new ArbDiscrEmpiricalDistFunc_3D(mfd.getMinX(),mfd.getMaxX(),mfd.size()); + for (ObsEqkRupList catalog : catList) { + for (ObsEqkRupture obsRup : catalog) { + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double rupLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + if(randomFromCOV) { + if(rupLoss==0) { + arbFunc3D.set(mfd.getClosestXIndex(rup.getMag()), rupLoss, 1.0); + } + else { + double randLoss = randomLossForEventID[rup.getID()]; // covModel.getDistribution(rupLoss).sample(); // get random sample + arbFunc3D.set(mfd.getClosestXIndex(rup.getMag()), randLoss, 1.0); + } + } + else + arbFunc3D.set(mfd.getClosestXIndex(rup.getMag()), rupLoss, 1.0); + } + } + +// ArbDiscrEmpiricalDistFunc[] array = arbFunc3D.getArbDiscrEmpDistFuncArray(); +// for(int i=0;i funcList = new ArrayList(); + funcList.add(arbFunc3D.getMeanCurve()); + funcList.add(arbFunc3D.getInterpolatedFractileCurve(0.5)); + funcList.add(arbFunc3D.getInterpolatedFractileCurve(0.025)); + funcList.add(arbFunc3D.getInterpolatedFractileCurve(0.975)); + funcList.add(arbFunc3D.getInterpolatedFractileCurve(0.16)); + funcList.add(arbFunc3D.getInterpolatedFractileCurve(0.84)); + funcList.get(0).setName("Mean"); + funcList.get(1).setName("Median"); + funcList.get(2).setName("2.5% fractile"); + funcList.get(3).setName("97.5% fractile"); + funcList.get(4).setName("16% fractile"); + funcList.get(5).setName("84% fractile"); + ArrayList plotChars = new ArrayList(); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLUE)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.DOTTED_AND_DASHED, 2f, Color.BLUE)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.DOTTED, 2f, Color.BLUE)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.DOTTED, 2f, Color.BLUE)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 2f, Color.BLUE)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 2f, Color.BLUE)); + + PlottingUtils.writeAndOrPlotFuncs(funcList, plotChars, "", "Magnitude", "Mean Loss ($1000)", + new Range(5d,8.5), new Range(1e-2,1e8), false, true, 6d, 5d, null, true); + } + + /** + * This computes various stats from catalogs assuming all have duration = this.catalogDuration + * Output: + * stats[0] = AAL (average annual loss); equals EventRate*MeanLoss + * stats[1] = Event Rate (per year) + * stats[2] = Mean Loss + * stats[3] = Median Loss + * stats[4] = Loss COV + * stats[5] = minLoss; // minimum non-zero loss + * stats[6] = maxLoss; + * stats[7] = fractZeroLosses; // fraction of ruptures that are zero-loss + * @param catList + * @return stats + */ + public double[] getLossStatsFromCatalogList(ArrayList catList, boolean mkPlot) { + HistogramFunction hist = getBlankLossCurveLogX(); + double[] stats = new double[8]; + ArbDiscrEmpiricalDistFunc distFunc = new ArbDiscrEmpiricalDistFunc(); + double aal =0; + int n=0; + int numZeroLosses = 0; + double minLoss = Double.MAX_VALUE, maxLoss=0; + for (ObsEqkRupList catalog : catList) { + for (ObsEqkRupture obsRup : catalog) { + n+=1; + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double rupLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + + if(rupLoss==0.0) + numZeroLosses+=1; + else + if(minLoss>rupLoss) minLoss=rupLoss; + + if(maxLoss 1e-6) + throw new RuntimeException("AAL != EventRate*MeanLoss"); + + if(mkPlot) { + hist.scale(1.0/(catList.size()*catalogDuration)); +// // this should about equal the aal; it does +// double aalTest =0; +// for(int i=1;i catList) { + +// // commented out stuff was looking at empirical distributions +// ArbDiscrEmpiricalDistFunc distFunc = new ArbDiscrEmpiricalDistFunc(); +// +// int numPts = 10000000; +// double maxX = Math.pow(10,lossCurveLogMax); +// double delta = maxX/(double)numPts; +// double minX = 0.5*delta; +// HistogramFunction histFunc = new HistogramFunction(minX, numPts, delta); + + double aal =0; + int n=0; + for (ObsEqkRupList catalog : catList) { + for (ObsEqkRupture obsRup : catalog) { + n+=1; + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double rupLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + aal += rupLoss; +// distFunc.set(rupLoss, 1.0); +// histFunc.add(rupLoss, 1.0); + } + } + aal /= (double)catList.size()*catalogDuration; +// double meanRateEvents = n/((double)catList.size()*catalogDuration); +// System.out.println("HERE: mean = "+distFunc.getMean()+"; median = "+ +// distFunc.getMedian()+"; mode = "+distFunc.getMostCentralMode()+"; COV = "+distFunc.getCOV() +// +"; AAL = "+distFunc.getMean()*meanRateEvents); +// histFunc.normalizeBySumOfY_Vals(); +//// double median = histFunc.getCumulativeDistFunctionWithHalfBinOffset().getFirstInterpolatedX(0.5); +//// double median = histFunc.getCumulativeDistFunctionWithHalfBinOffset().getFirstInterpolatedX_inLogYDomain(0.5); +// double median = histFunc.getCumulativeDistFunctionWithHalfBinOffset().getFirstInterpolatedX_inLogXLogYDomain(0.5); +// System.out.println("HERE2: mean = "+histFunc.computeMean()+"; median = "+median+ +// "; mode = "+histFunc.getMode()+"; AAL = "+histFunc.computeMean()*meanRateEvents); +// System.out.println("HERE2: first X bin: "+(histFunc.getX(0)-histFunc.getDelta()/2d)+" to "+(histFunc.getX(0)+histFunc.getDelta()/2d)); + return aal; + } + + + + public double getAveAnnualLossFromERF(boolean timeDep) { + double aal =0; + if(erf==null) + erf = getERF(); + + if(timeDep) + erf.getParameter(ProbabilityModelParam.NAME).setValue(ProbabilityModelOptions.U3_BPT); + else + erf.getParameter(ProbabilityModelParam.NAME).setValue(ProbabilityModelOptions.POISSON); + erf.updateForecast(); + +// if(D)System.out.println("Prob Model: "+erf.getParameter(ProbabilityModelParam.NAME).getValue()); + + double duration = erf.getTimeSpan().getDuration(); + double maxMagRejected = 0; + for(int n=0;n funcList = new ArrayList(); + funcList.add(lossExceedProbCurve); + funcList.add(lossExceedProbCurve2); + funcList.add(lossExceedProbCurve3); + return funcList; + } + + + + + public ArbitrarilyDiscretizedFunc getLossRateCurveFromERF() { + ArbitrarilyDiscretizedFunc lossRateCurve = getBlankLossCurve(0d); + if(erf==null) + erf = getERF(); + erf.getParameter(ProbabilityModelParam.NAME).setValue(ProbabilityModelOptions.POISSON); + erf.updateForecast(); + double duration = erf.getTimeSpan().getDuration(); + for(int n=0;n catList) { + if(D) System.out.println("Working on getLossRateCurveFromCatalogs(); num catalogs = "+catList.size()); + ArbitrarilyDiscretizedFunc lossRateCurve = getBlankLossCurve(0d); + for (ObsEqkRupList catalog : catList) { + for (ObsEqkRupture obsRup : catalog) { + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + if(aveLoss == 0.0) + continue; + DiscretizedFunc condCurve = covModel.calcLossExceedanceProbs(lossRateCurve, aveLoss); + for(int i=0;i catList) { + ArbitrarilyDiscretizedFunc lossRateCurve = getBlankLossCurve(0d); + for (ObsEqkRupList catalog : catList) { + for (ObsEqkRupture obsRup : catalog) { + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + if(aveLoss == 0.0) + continue; + for(int i=0;i catList, double duration) { + + if(D) System.out.println("Working on getLossExceedProbCurveFromCatalogs(); num catalogs = "+catList.size()); + + ArbitrarilyDiscretizedFunc lossProbCurve = getBlankLossCurve(0d); + ArrayList subCatList = getSubcatalogList(catList, duration); + + for (ObsEqkRupList catalog : subCatList) { + ArbitrarilyDiscretizedFunc lossProbCurveForSubCat = getBlankLossCurve(1d); + for (ObsEqkRupture obsRup : catalog) { + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + if(aveLoss == 0.0) + continue; + DiscretizedFunc condCurve = covModel.calcLossExceedanceProbs(lossProbCurveForSubCat, aveLoss); + for(int i=0;ilossProbCurve.getY(i)) { // numerical problem + if(D) System.out.println("ISSUE: Low bound above mean:\nmean="+lossProbCurve.getY(i)+"\nlowBound="+confArray[0]+"\nN="+subCatList.size()); + confArray[0] = 0; + } + low95confCurve.set(xVal,confArray[0]); + high95confCurve.set(xVal,confArray[1]); + } + + UncertainArbDiscFunc uncertFunc = new UncertainArbDiscFunc(lossProbCurve,low95confCurve,high95confCurve); + return uncertFunc; + } + + + + + + + public ArbitrarilyDiscretizedFunc getLossExceedProbCurveFromCatalogsZeroCOV(ArrayList catList, double duration) { + + ArbitrarilyDiscretizedFunc lossProbCurve = getBlankLossCurve(0d); + ArrayList subCatList = getSubcatalogList(catList, duration); + + for (ObsEqkRupList catalog : subCatList) { + ArbitrarilyDiscretizedFunc lossProbCurveForSubCat = getBlankLossCurve(0d); + double maxLoss = 0; + for (ObsEqkRupture obsRup : catalog) { + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + if(maxLoss catList, double duration) { + + ArbitrarilyDiscretizedFunc lossProbCurve = getBlankLossCurve(0d); + ArrayList subCatList = getSubcatalogList(catList, duration); + + // TEMP + double lossExceedThresh = 4.0e4; + int[] numExceedForCatArray = new int[subCatList.size()]; + + // TEMP + int c=-1; + for (ObsEqkRupList catalog : subCatList) { + c+=1; + + ArbitrarilyDiscretizedFunc lossProbCurveForSubCat = getBlankLossCurve(0d); + double maxLoss = 0; + for (ObsEqkRupture obsRup : catalog) { + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + if(aveLoss==0) + continue; + double randLoss = randomLossForEventID[rup.getID()]; // covModel.getDistribution(aveLoss).sample(); // get randome sample + if(maxLoss=lossExceedThresh) + numExceedForCatArray[c] += 1; +// if(randLoss>=Math.pow(10, 0.95) && randLoss250) { +// for (ObsEqkRupture obsRup : subCatList.get(c-1)) { // previous catalog +// ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; +// double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); +// if(aveLoss==0) +// continue; +// System.out.println(rup.getMag()+"\t"+rup.getGeneration()+"\tPrevious catalog"); // +"\t"+rup.getOriginTimeCal()); +// } +// for (ObsEqkRupture obsRup : catalog) { // this catalog +// ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; +// double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); +// if(aveLoss==0) +// continue; +// System.out.println(rup.getMag()+"\t"+rup.getGeneration()); // +"\t"+rup.getOriginTimeCal()); +// } +// System.exit(0); +// } + + + } + lossProbCurve.scale(1.0/((double)subCatList.size())); + ArbitrarilyDiscretizedFunc low95confCurve = getBlankLossCurve(0d); + ArbitrarilyDiscretizedFunc high95confCurve = getBlankLossCurve(0d); + for(int i=0;ilossProbCurve.getY(i)) { // numerical problem + if(D) System.out.println("ISSUE: Low bound above mean:\nmean="+lossProbCurve.getY(i)+"\nlowBound="+confArray[0]+"\nN="+subCatList.size()); + confArray[0] = 0; + } + low95confCurve.set(xVal,confArray[0]); + high95confCurve.set(xVal,confArray[1]); + } + + UncertainArbDiscFunc uncertFunc = new UncertainArbDiscFunc(lossProbCurve,low95confCurve,high95confCurve); + + // TEMP + int maxNum = 0; + for(int num:numExceedForCatArray) + if(num>maxNum) maxNum=num; +// System.out.println("maxNum: "+maxNum); + HistogramFunction hist = new HistogramFunction(0d,(double)maxNum,maxNum+1); + for(int num:numExceedForCatArray) + hist.add(num, 1d); + // System.out.println("\nhist.calcSumOfY_Vals():\n"+hist.calcSumOfY_Vals()); + + hist.normalizeBySumOfY_Vals(); + System.out.println("\nnumExceedPDF:\n"+hist); + + return uncertFunc; + } + + + + public UncertainArbDiscFunc getLossIncrProbDistFromCatalogsRandomFromCOV(ArrayList catList, double duration) { + + ArbitrarilyDiscretizedFunc lossProbCurve = getBlankLossCurve(0d); + HistogramFunction lossHistLogX = getBlankLossCurveLogX(); + ArrayList subCatList = getSubcatalogList(catList, duration); + + for (ObsEqkRupList catalog : subCatList) { + + ArbitrarilyDiscretizedFunc lossProbCurveForSubCat = getBlankLossCurve(0d); + double maxLoss = 0; + for (ObsEqkRupture obsRup : catalog) { + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + if(aveLoss==0) + continue; + double randLoss = randomLossForEventID[rup.getID()]; // covModel.getDistribution(aveLoss).sample(); // get randome sample +// if(maxLosslossProbCurve.getY(i)) { // numerical problem + if(D) System.out.println("ISSUE: Low bound above mean:\nmean="+lossProbCurve.getY(i)+"\nlowBound="+confArray[0]+"\nN="+subCatList.size()); + confArray[0] = 0; + } + low95confCurve.set(xVal,confArray[0]); + high95confCurve.set(xVal,confArray[1]); + } + + UncertainArbDiscFunc uncertFunc = new UncertainArbDiscFunc(lossProbCurve,low95confCurve,high95confCurve); + + return uncertFunc; + } + + + + + + public UncertainArbDiscFunc getAggrLossExceedProbCurveFromCatalogsRandomFromCOV(ArrayList catList, double duration) { + + ArbitrarilyDiscretizedFunc lossProbCurve = getBlankLossCurve(0d); + ArrayList subCatList = getSubcatalogList(catList, duration); + + double aveAggrLoss = 0; + for (ObsEqkRupList catalog : subCatList) { + ArbitrarilyDiscretizedFunc lossProbCurveForSubCat = getBlankLossCurve(0d); + double totLoss = 0; + for (ObsEqkRupture obsRup : catalog) { + ETAS_EqkRupture rup = (ETAS_EqkRupture)obsRup; + double aveLoss = meanLossForNthRup.get(rup.getNthERF_Index()); + if(aveLoss==0) + continue; + double randLoss = randomLossForEventID[rup.getID()]; // covModel.getDistribution(aveLoss).sample(); // get random sample + totLoss+=randLoss; +// totLoss+=aveLoss; + } + aveAggrLoss+=totLoss; + for(int i=0;i declusteredList = new ArrayList(); for(ObsEqkRupList rupList: catalogList) { ObsEqkRupList declusteredCatalog = GardnerKnopoffDeclustering.getDeclusteredCatalog(rupList); declusteredList.add(declusteredCatalog); } + if(D) System.out.println("Done Declustering"); return declusteredList; } @@ -175,7 +1063,7 @@ public static ArrayList getSpontaneousEventsCatalog(ArrayList getU3_GK_FilteredCatalog(ArrayList catalogList) { + public ArrayList getU3_GK_FilteredCatalog(ArrayList catalogList) { int numGriddedRups=0; int numFSS_Rups=0; GardnerKnopoffAftershockFilter gkFilter = new GardnerKnopoffAftershockFilter(0.05, 9.95, 100); @@ -183,15 +1071,15 @@ public static ArrayList getU3_GK_FilteredCatalog(ArrayList getU3_GK_FilteredCatalog(ArrayList getRandomizedCatalogs(Random random) { - if(catalogRandomizedList==null) { - if(D) System.out.println("Making randomized catalogs"); - if(random==null) random = new Random(); - catalogRandomizedList = getRandomizedCatalogs(this.getCatalogs(), random); - } - return catalogRandomizedList; - } - - - /** - * This not only randomizes event times within each catalog, but also among catalogs - * so there isn't a high variance of total rates among catalogs. This clones the - * ruptures and modifies the origin times. + * This randomizes the original/full catalogs (creating catalogRandomizedList). This not only randomizes + * event times within each catalog, but also among catalogs so there isn't a high variance of total rates + * among catalogs. This clones the ruptures and modifies the origin times. + * @param random - supply if reproducibility is desired (null otherwise) + * @return */ - public static ArrayList getRandomizedCatalogs(ArrayList catalogList, Random random) { - if(random==null) - random = new Random(); - ArrayList catalogRandomizedList = new ArrayList(); - for(int i=0;icatalogStartYear+catalogDuration) - throw new RuntimeException("ERROR: "+randomTimeYr); -// rup.setOriginTime((long)((double)millisPerYr*(randomTimeYr-1970))); - ObsEqkRupture newRup = (ObsEqkRupture) rup.clone(); - newRup.setOriginTime((long)((double)millisPerYr*(randomTimeYr-1970))); -// catalogRandomizedList.get(ithCatalog).add(rup); - catalogRandomizedList.get(ithCatalog).add(newRup); + public ArrayList getRandomizedCatalogs() { + if(catalogRandomizedList==null) { + if(D) System.out.println("Making randomized catalogs"); + catalogRandomizedList = new ArrayList(); + for(int i=0;i getRandomizedCatalogs(ArrayList makeMFDs(boolean mkPopUpPlots) { + public ArrayList old_makeCatalogMFDs(boolean mkPopUpPlots) { - IncrementalMagFreqDist mfd = makeMFD(catalogList); + IncrementalMagFreqDist mfd = makeCatalogMFD(catalogList); double bVal = Math.log10(mfd.getY(5.05)/mfd.getY(6.95))/(6.95-5.05); mfd.setName("Full TD Catalog"); mfd.setInfo("bVal = "+(float)bVal); - IncrementalMagFreqDist mfdDeclustered = makeMFD(catalogDeclusteredList); + IncrementalMagFreqDist mfdDeclustered = makeCatalogMFD(catalogDeclusteredList); bVal = Math.log10(mfdDeclustered.getY(5.05)/mfdDeclustered.getY(6.95))/(6.95-5.05); mfdDeclustered.setInfo("GK Declustered"); mfdDeclustered.setInfo("bVal = "+(float)bVal); @@ -273,7 +1150,7 @@ public ArrayList makeMFDs(boolean mkPopUpPlots) { // mfdDeclustered_U3_GKfilter.setName("Declustered with U3 GK Filter"); // mfdDeclustered_U3_GKfilter.setInfo("bVal = "+(float)bVal); - IncrementalMagFreqDist mfdSpontaneousEventsOnly = makeMFD(getSpontaneousEventsCatalog(catalogList)); + IncrementalMagFreqDist mfdSpontaneousEventsOnly = makeCatalogMFD(getSpontaneousEventsCatalog(catalogList)); bVal = Math.log10(mfdSpontaneousEventsOnly.getY(5.05)/mfdSpontaneousEventsOnly.getY(6.95))/(6.95-5.05); mfdSpontaneousEventsOnly.setName("Spontaneous Events Only"); mfdSpontaneousEventsOnly.setInfo("bVal = "+(float)bVal); @@ -362,120 +1239,354 @@ public ArrayList makeMFDs(boolean mkPopUpPlots) { return returnList; } - - public static ArrayList loadCatalogs(File fssFile, File catalogsFile, double minMag) throws IOException, DocumentException { - + + 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)); + } catch (IOException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } + + FaultSystemSolutionERF_ETAS erf = ETAS_Launcher.buildERF(sol,false, 1d, 2012); + erf.updateForecast(); + return erf; + } + + + /** + * The check against the ERF does not work because indexing between ERF and catalogs somehow got screwed up, + * but this shouldn't matter here. The FSS solution ruptures do not have the finite surfaces populated unless + * checkAgainstERF is set as true. + * @param fssFile + * @param catalogsFile + * @param minMag + * @param checkAgainstERF - checks things against ERF and adds the finite rupture surfaces + * @return + * @throws IOException + * @throws DocumentException + */ + public ArrayList loadCatalogs(File catalogsFile, double minMag, boolean checkAgainstERF) throws IOException, DocumentException { + + + List catalogs = ETAS_CatalogIO.loadCatalogsBinary(catalogsFile, minMag); + + System.out.println("catalogs.size()="+catalogs.size()); + + // This checks things against ERF and adds the finite rupture surfaces + if(checkAgainstERF & D) { + + if(erf==null) + erf = getERF(); + + int numFSS_Rups = erf.getTotNumRupsFromFaultSystem(); + + double minERF_Mag = 10000; + for(int r=0;rerf.getNthRupture(numFSS_Rups).getMag()) + minERF_Mag=erf.getNthRupture(numFSS_Rups).getMag(); + + System.out.println("ERF totNumRups: "+erf.getTotNumRups()); + System.out.println("ERF getTotNumRupsFromFaultSystem: "+numFSS_Rups); + System.out.println("ERF getNumFaultSystemSources: "+erf.getNumFaultSystemSources()); + System.out.println("ERF getNumSources: "+erf.getNumSources()); + System.out.println("ERF minMag: "+minERF_Mag); + System.out.println("ERF numGridCells: "+erf.getGridSourceProvider().getGriddedRegion().getNumLocations()); + + double[] numProbForGridIndex = new double[erf.getGridSourceProvider().getGriddedRegion().getNumLocations()]; + + // set the finite rupture surfaces + int numProbMags =0; + int numProbNodeIndices =0; + double minCatalogMag = 10000; + int maxCatalogNthIndex = 0; + int totNumCatRups=0; + int totFSS_catRups=0; + int totGridSeisCatRups=0; + HashMap numProbForGeneration = new HashMap(); + HistogramFunction histFunc = new HistogramFunction(0.5, 10, 1.0); + for (ETAS_Catalog catalog : catalogs) { + for (ETAS_EqkRupture rup : catalog) { + int nth = rup.getNthERF_Index(); + int rupFSS_Index = rup.getFSSIndex(); + + totNumCatRups+=1; + if(rupFSS_Index != -1) + totFSS_catRups += 1; + else + totGridSeisCatRups += 1; + + + if(minCatalogMag>rup.getMag()) + minCatalogMag=rup.getMag(); + if(maxCatalogNthIndex0) { +// System.out.println("\t"+i+"\t"+numProbForGridIndex[i]); + n+=1; + } + System.out.println(histFunc); + System.out.println("\t% of grid nodes with leakage: "+(double)n/(double)numProbForGridIndex.length); + + System.out.println("\tThis shows the problem is for triggered events and not spontaneous ones"); + for(int gen:numProbForGeneration.keySet()) + System.out.println("\tgen "+gen+":\t"+numProbForGeneration.get(gen)); + + System.out.println("\ntotNumCatRups="+totNumCatRups); + System.out.println("totFSS_catRups="+totFSS_catRups); + System.out.println("totGridSeisCatRups="+totGridSeisCatRups); + System.out.println("num problematic mags="+numProbMags); + System.out.println("minCatalogMag="+minCatalogMag); + System.out.println("minERF_Mag="+minERF_Mag); + System.out.println("maxCatalogNthIndex="+maxCatalogNthIndex); + + } + + // convert to a list of ObsEqkRupList objects (and filter M≤5 events out) + ArrayList obsEqkRupListList = new ArrayList(); + + int id = -1; // reset the IDs so I can store random loss samples by this. + for (ETAS_Catalog catalog : catalogs) { + ObsEqkRupList obsEqkRupList = new ObsEqkRupList(); + for (ETAS_EqkRupture rup : catalog) { + if(rup.getMag()>=minMag) { + id+=1; + rup.setID(id); + obsEqkRupList.add(rup); + } + } + obsEqkRupListList.add(obsEqkRupList); + } + totalNumEvents = id+1; + + return obsEqkRupListList; + } + + + private static HistogramFunction getBlankLossCurveLogX() { + return new HistogramFunction(lossCurveLogMin, lossCurveLogMax, lossCurveNum); + } + + private static ArbitrarilyDiscretizedFunc getBlankLossCurve(double initYvalues) { + ArbitrarilyDiscretizedFunc func = new ArbitrarilyDiscretizedFunc(); + HistogramFunction histFunc = getBlankLossCurveLogX(); + for (int i=0; i catalogs = ETAS_CatalogIO.loadCatalogsBinary(catalogsFile, minMag); + hist.setName("Histogram from "+curve.getName()); + double mean =0; + for(int i=0;i0 && y2>0) + diff2 = Math.abs(1.0-y1/y2); + else + diff2 = Math.abs(y1-y2); + if(diff1>0.00001 || diff2>0.00001) + throw new RuntimeException("Problem with convertExceedCurveToLog10_PDF_Hist()\t"+x1+"\t"+x2+"\t"+y1+"\t"+y2+"\t"+diff1+"\t"+diff2); + } +// System.exit(0); - FaultSystemSolutionERF_ETAS erf = ETAS_Launcher.buildERF(sol,false, 1d, 2012); - erf.updateForecast(); - int numFSS_Rups = erf.getTotNumRupsFromFaultSystem(); - - double minERF_Mag = 10000; - for(int r=0;rerf.getNthRupture(numFSS_Rups).getMag()) - minERF_Mag=erf.getNthRupture(numFSS_Rups).getMag(); - - System.out.println("Sol has RupMFDs? "+sol.hasModule(RupMFDsModule.class)); - System.out.println("ERF totNumRups: "+erf.getTotNumRups()); - System.out.println("ERF getTotNumRupsFromFaultSystem: "+numFSS_Rups); - System.out.println("ERF getNumFaultSystemSources: "+erf.getNumFaultSystemSources()); - System.out.println("ERF getNumSources: "+erf.getNumSources()); - System.out.println("ERF minMag: "+minERF_Mag); - - // set the finite rupture surfaces - String errorMessage = null; - int numProbMags =0; - int numProbNodeIndices =0; - double minCatalogMag = 10000; - int maxCatalogNthIndex = 0; - for (ETAS_Catalog catalog : catalogs) { - for (ETAS_EqkRupture rup : catalog) { - int nth = rup.getNthERF_Index(); - int rupFSS_Index = rup.getFSSIndex(); - - if(minCatalogMag>rup.getMag()) - minCatalogMag=rup.getMag(); - if(maxCatalogNthIndex mfdList = new ArrayList(); + mfdList.add(lossRate_cat); +// mfdList.add(lossRate_erf); + mfdList.add(lossRate_catDecl); + mfdList.add(lossRate_spontaneous); + ArrayList plotChars = new ArrayList(); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLUE)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLACK)); + plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.GREEN)); + // new Range(1e3,2e4) + PlottingUtils.writeAndOrPlotFuncs(mfdList, plotChars, "Loss Rate vs Mag", "Mag", "Loss Rate (/yr)", + null, new Range(1e3,2e4), false, true, null, true); + return mfdList; + } + - private static IncrementalMagFreqDist makeMagNumDist(ObsEqkRupList rupList) { + private static IncrementalMagFreqDist old_makeMagNumDist(ObsEqkRupList rupList) { IncrementalMagFreqDist mfd = makeMFD(); for(ObsEqkRupture rup:rupList) mfd.add(rup.getMag(), 1.0); @@ -514,22 +1686,14 @@ private static void addToMFD(ObsEqkRupList rupList, IncrementalMagFreqDist mfd) mfd.add(rup.getMag(), 1.0); } } - - - - - - - - - public void computeNumEventDistribution(ArrayList catalogList, double duration, double[] magArray) { + public void old_computeNumEventDistribution(ArrayList catalogList, double duration, double[] magArray) { // get subcatalogs ArrayList eqkRupListList = getSubcatalogList(catalogList, duration); - EvenlyDiscretizedFunc tempFunc = makeMagNumDist(eqkRupListList.get(0)).getCumRateDistWithOffset(); + EvenlyDiscretizedFunc tempFunc = old_makeMagNumDist(eqkRupListList.get(0)).getCumRateDistWithOffset(); ArbDiscrEmpiricalDistFunc_3D cumMFD_FromAllCatalogsFunc_3D = new ArbDiscrEmpiricalDistFunc_3D(tempFunc.getMinX(), tempFunc.size(), tempFunc.getDelta()); int showProgressAt = eqkRupListList.size()/10; @@ -538,7 +1702,7 @@ public void computeNumEventDistribution(ArrayList catalogList, do if(D) System.out.println(i); showProgressAt+=eqkRupListList.size()/10; } - IncrementalMagFreqDist func = makeMagNumDist(eqkRupListList.get(i)); + IncrementalMagFreqDist func = old_makeMagNumDist(eqkRupListList.get(i)); cumMFD_FromAllCatalogsFunc_3D.set(func.getCumRateDistWithOffset(), 1.0); } @@ -568,12 +1732,12 @@ public void computeNumEventDistribution(ArrayList catalogList, do ArrayList plottingFuncsArray = new ArrayList(); plottingFuncsArray.add(numDist); plottingFuncsArray.add(poissNumDist); - this.quickPlot(plottingFuncsArray,"Num","Probability","M="+mag+"; Duration="+duration); + this.quickPlot(plottingFuncsArray,"Num","Probability","M="+mag+"; Duration="+duration, false, false); } } - + /** * This divides the catalogs up into subcatalogs of the given duration @@ -581,30 +1745,94 @@ public void computeNumEventDistribution(ArrayList catalogList, do * @param duration - yrs * @return */ - public static ArrayList getSubcatalogList(ArrayList catalogList, double duration) { + public static ArrayList old_getSubcatalogList(ArrayList catalogList, double duration) { int numSubCatPerCat = (int)Math.floor((catalogDuration+0.003)/duration); // add a day (0.03) to make sure we get the last window if(D) System.out.println("numSubCatPerCat = "+numSubCatPerCat); + int numRupsNotUsed = 0; + for(ObsEqkRupList catalog : catalogList) + numRupsNotUsed += catalog.size(); + ArrayList eqkRupListList = new ArrayList(); for(ObsEqkRupList catalog : catalogList) { long endEpoch = (long)((catalogStartYear-1970)*millisPerYr) + (long)(duration*millisPerYr); ObsEqkRupList currEqkList = new ObsEqkRupList(); eqkRupListList.add(currEqkList); for(ObsEqkRupture rup: catalog) { - if(rup.getOriginTime() funcList1 = new ArrayList(); + funcList1.add(normEvenPDF); + funcList1.add(normArbPDF); + quickPlot(funcList1, "X", "Prob", "Normal Dist", false, false); + + ArrayList funcList2 = new ArrayList(); + funcList2.add(logNormEvenPDF); + funcList2.add(logNormArbPDF); + quickPlot(funcList2, "X", "Prob", "LogNormal Dist", false, false); + + } + + + + + /** + * + * @param funcList + * @param includeUncert - if true, columns with normalized standard deviation (std/loss) are added + * @return + */ + public static String getTableStringOfCEA_Values(ArrayList funcList, boolean includeUncert) { + String tableString = ""; + for(int row = 0; row funcList, boolean includeUncert) { + String tableString = ""; + for(int row = 0; row funcList = new ArrayList(); + + String randCOVincluded = ""; + if(randCOV) + randCOVincluded = " (cond loss PDF for each rup randomly sampled)"; + + if(randCOV) { + UncertainArbDiscFunc lossExceedProbFromTD_CatalogsRandFromCOV = getLossExceedProbCurveFromCatalogsRandomFromCOV(getCatalogs(), 1d); + lossExceedProbFromTD_CatalogsRandFromCOV.setName("lossExceedProbFromTD_CatalogsRandFromCOV"); + UncertainArbDiscFunc lossExceedProbFromPoisCatalogsRandFromCOV = getLossExceedProbCurveFromCatalogsRandomFromCOV(getRandomizedCatalogs(), 1d); + lossExceedProbFromPoisCatalogsRandFromCOV.setName("lossExceedProbFromPoisCatalogsRandFromCOV"); + UncertainArbDiscFunc lossExceedProbFromDeclusteredCatalogsRandFromCOV = getLossExceedProbCurveFromCatalogsRandomFromCOV(getCatalogsDeclustered(), 1d); + lossExceedProbFromDeclusteredCatalogsRandFromCOV.setName("lossExceedProbFromDeclusteredCatalogsRandFromCOV"); + UncertainArbDiscFunc lossExceedProbFromSpontaneousCatalogsRandFromCOV = getLossExceedProbCurveFromCatalogsRandomFromCOV(getSpontaneousEventsCatalog(getCatalogs()), 1d); + lossExceedProbFromSpontaneousCatalogsRandFromCOV.setName("lossExceedProbFromSpontaneousCatalogsRandFromCOV"); + funcList.add(lossExceedProbFromTD_CatalogsRandFromCOV); + funcList.add(lossExceedProbFromPoisCatalogsRandFromCOV); + funcList.add(lossExceedProbFromDeclusteredCatalogsRandFromCOV); + funcList.add(lossExceedProbFromSpontaneousCatalogsRandFromCOV); + } + else { + UncertainArbDiscFunc lossExceedProbFromTD_Catalogs = getLossExceedProbCurveFromCatalogs(getCatalogs(), 1d); + lossExceedProbFromTD_Catalogs.setName("lossExceedProbFromTD_Catalogs"); + UncertainArbDiscFunc lossExceedProbFromPoisCatalogs = getLossExceedProbCurveFromCatalogs(getRandomizedCatalogs(), 1d); + lossExceedProbFromPoisCatalogs.setName("lossExceedProbFromPoisCatalogs"); + UncertainArbDiscFunc lossExceedProbFromDeclusteredCatalogs = getLossExceedProbCurveFromCatalogs(getCatalogsDeclustered(), 1d); + lossExceedProbFromDeclusteredCatalogs.setName("lossExceedProbFromDeclusteredCatalogs"); + UncertainArbDiscFunc lossExceedProbFromSpontaneousCatalogs = getLossExceedProbCurveFromCatalogs(getSpontaneousEventsCatalog(getCatalogs()), 1d); + lossExceedProbFromSpontaneousCatalogs.setName("lossExceedProbFromSpontaneousCatalogs"); + funcList.add(lossExceedProbFromTD_Catalogs); + funcList.add(lossExceedProbFromPoisCatalogs); + funcList.add(lossExceedProbFromDeclusteredCatalogs); + funcList.add(lossExceedProbFromSpontaneousCatalogs); + } + System.out.println("\n1yr Loss for CEA probability levels"+randCOVincluded+":\n"); + System.out.println("Probability\tTD_loss\tTD_uncert\tTI_loss\tTI_uncert\tGK_loss\tGK_uncert\tSP_loss\tSP_uncert"); + System.out.println(getTableStringOfCEA_Values(funcList,true)); + System.out.println("\n1yr Loss ratios, relative to TD, for CEA probability levels"+randCOVincluded+":\n"); + System.out.println("Probability\tTIvsTD_loss\tTIvsTD_uncert\tGKvsTD_loss\tGKvsTD_uncert\tSPvsTD_loss\tSPvsTD_uncert"); + System.out.println(getTableStringOfCEA_Ratios(funcList, true)); + System.out.println("TI = Time Dependent\nTI = Time Independent (Poisson)\nGK = Gardner Knopoff declustered\nSP = Spontaneous ETAS events\nuncert = 1-sigma uncertainty normalized by the mean\n"); + ArrayList plotChars2 = new ArrayList(); + plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLUE)); + plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.RED)); + plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLACK)); + plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.GREEN)); + String fileNamePrefix = dirName+"/FigureA_Losses"; + if(randCOV) fileNamePrefix += "_randCOV"; + PlottingUtils.writeAndOrPlotUncertFuncs(funcList, plotChars2, "", "Loss ($1000)", "Probability", new Range(1e3,1e9), new Range(1e-6,1.0), true, true, fileNamePrefix, true); + + // Make ratio plots + UncertainArbDiscFunc denomCurve = funcList.get(0); + ArrayList funcListRatios = new ArrayList(); + + double xThresh = denomCurve.getClosestXtoY(1e-3); + for(int i=0;ixThresh) + break; + double ratio = funcList.get(i).getY(j)/denomCurve.getY(j); + if(ratio<1e-9 || ratio>1e9) ratio = 1e-9; + ratioCurve.set(denomCurve.getX(j), ratio); + } + funcListRatios.add(ratioCurve); + } + fileNamePrefix = dirName+"/FigureA_LossRatios"; + if(randCOV) fileNamePrefix += "_randCOV"; + PlottingUtils.writeAndOrPlotFuncs(funcListRatios, plotChars2, "test", "Loss ($1000)", "Loss Ratio", new Range(1e3,xThresh), new Range(0.0,2.0), true, false, fileNamePrefix, true); + + ArrayList pdfList = new ArrayList(); + for(int i=0;i<2;i++) { + pdfList.add(convertExceedCurveToLog10_PDF_Hist(funcList.get(i))); + } + PlottingUtils.writeAndOrPlotFuncs(pdfList, plotChars2, "", "Loss ($1000)", "Probability", null, null, false, false, null, true); + + } + + + + + public void makeFigAndTableSetB_1yr_Exceedances() { + + ArrayList funcList = new ArrayList(); + + boolean randCOV = true; + String randCOVincluded = ""; + if(randCOV) + randCOVincluded = " (cond loss PDF for each rup randomly sampled)"; + + UncertainArbDiscFunc aggrLossExceedProbFromTD_CatalogsRandFromCOV = getAggrLossExceedProbCurveFromCatalogsRandomFromCOV(getCatalogs(), 1d); + aggrLossExceedProbFromTD_CatalogsRandFromCOV.setName("aggrLossExceedProbFromTD_CatalogsRandFromCOV"); + UncertainArbDiscFunc aggrLossExceedProbFromPoisCatalogsRandFromCOV = getAggrLossExceedProbCurveFromCatalogsRandomFromCOV(getRandomizedCatalogs(), 1d); + aggrLossExceedProbFromPoisCatalogsRandFromCOV.setName("aggrLossExceedProbFromPoisCatalogsRandFromCOV"); + UncertainArbDiscFunc aggrLossExceedProbFromDeclusteredCatalogsRandFromCOV = getAggrLossExceedProbCurveFromCatalogsRandomFromCOV(getCatalogsDeclustered(), 1d); + aggrLossExceedProbFromDeclusteredCatalogsRandFromCOV.setName("aggrLossExceedProbFromDeclusteredCatalogsRandFromCOV"); + UncertainArbDiscFunc aggrLossExceedProbFromSpontaneousCatalogsRandFromCOV = getAggrLossExceedProbCurveFromCatalogsRandomFromCOV(getSpontaneousEventsCatalog(getCatalogs()), 1d); + aggrLossExceedProbFromSpontaneousCatalogsRandFromCOV.setName("aggrLossExceedProbFromSpontaneousCatalogsRandFromCOV"); + funcList.add(aggrLossExceedProbFromTD_CatalogsRandFromCOV); + funcList.add(aggrLossExceedProbFromPoisCatalogsRandFromCOV); + funcList.add(aggrLossExceedProbFromDeclusteredCatalogsRandFromCOV); + funcList.add(aggrLossExceedProbFromSpontaneousCatalogsRandFromCOV); + + System.out.println("\n1yr Loss for CEA probability levels"+randCOVincluded+":\n"); + System.out.println("Probability\tTD_loss\tTD_uncert\tTI_loss\tTI_uncert\tGK_loss\tGK_uncert\tSP_loss\tSP_uncert"); + System.out.println(getTableStringOfCEA_Values(funcList,true)); + System.out.println("\n1yr Loss ratios, relative to TD, for CEA probability levels"+randCOVincluded+":\n"); + System.out.println("Probability\tTIvsTD_loss\tTIvsTD_uncert\tGKvsTD_loss\tGKvsTD_uncert\tSPvsTD_loss\tSPvsTD_uncert"); + System.out.println(getTableStringOfCEA_Ratios(funcList, true)); + System.out.println("TI = Time Dependent\nTI = Time Independent (Poisson)\nGK = Gardner Knopoff declustered\nSP = Spontaneous ETAS events\nuncert = 1-sigma uncertainty normalized by the mean\n"); + ArrayList plotChars2 = new ArrayList(); + plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLUE)); + plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.RED)); + plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLACK)); + plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.GREEN)); + String fileNamePrefix = dirName+"/FigureB_AggrLosses"; + if(randCOV) fileNamePrefix += "_randCOV"; + PlottingUtils.writeAndOrPlotUncertFuncs(funcList, plotChars2, "", "Aggregate Loss ($1000)", "Probability", new Range(1e3,1e9), new Range(1e-6,1.0), true, true, fileNamePrefix, true); + + // Make ratio plots + UncertainArbDiscFunc denomCurve = funcList.get(0); + ArrayList funcListRatios = new ArrayList(); + + double xThresh = denomCurve.getClosestXtoY(1e-3); + for(int i=0;ixThresh) + break; + double ratio = funcList.get(i).getY(j)/denomCurve.getY(j); + if(ratio<1e-9 || ratio>1e9) ratio = 1e-9; + ratioCurve.set(denomCurve.getX(j), ratio); + } + funcListRatios.add(ratioCurve); + } + fileNamePrefix = dirName+"/FigureB_AggrLossRatios"; + if(randCOV) fileNamePrefix += "_randCOV"; + PlottingUtils.writeAndOrPlotFuncs(funcListRatios, plotChars2, "", "AggregateLoss ($1000)", "Loss Ratio", new Range(1e3,xThresh), new Range(0.0,2.0), true, false, fileNamePrefix, true); + + ArrayList pdfList = new ArrayList(); + for(int i=0;i<2;i++) { + pdfList.add(convertExceedCurveToLog10_PDF_Hist(funcList.get(i))); + } + PlottingUtils.writeAndOrPlotFuncs(pdfList, plotChars2, "", "Aggregate Loss ($1000)", "Probability", null, null, false, false, null, true); + } + + + /** + * This stores a random loss sample for each event (so we can have the same sampled value across different methods) + */ + public void setRandomLossForEvents() { + randomLossForEventID = new double[totalNumEvents]; + for (ObsEqkRupList catalog : catalogList) { + for (ObsEqkRupture rup : catalog) { + ETAS_EqkRupture etasRup = (ETAS_EqkRupture)rup; + double aveLoss = meanLossForNthRup.get(etasRup.getNthERF_Index()); + if(aveLoss==0) + randomLossForEventID[etasRup.getID()] = 0; + else { + double randLoss = covModel.getDistribution(aveLoss, randomGen).sample(); // get randome sample + randomLossForEventID[etasRup.getID()] = randLoss; + } + } + } + } + + + public static void main(String[] args) throws IOException, DocumentException { + + + /** + * NOTES: + * + * Modal and Median loss will always be zero if we go down to very small earthquakes (or at least it's + * always strongly magnitude dependent). + * + * AAL is an aggregate metric (summing all losses) + * + * + * + */ + + + + // make results reproducible + int seed = 123470; + U3ETAS_LossSimulationAnalysis analysis = new U3ETAS_LossSimulationAnalysis(seed); + +// analysis.writeAAL_ValuesBillions(); + +// UncertainArbDiscFunc lossExceedProbFromTD_CatalogsRandFromCOV = analysis.getLossExceedProbCurveFromCatalogsRandomFromCOV(analysis.getCatalogs(), 1d); +// UncertainArbDiscFunc lossExceedProbFromPoisCatalogsRandFromCOV = analysis.getLossExceedProbCurveFromCatalogsRandomFromCOV(analysis.getRandomizedCatalogs(), 1d); +// +// // trying to figure out what's wrong with the incremental distributions +// UncertainArbDiscFunc incrTestTD = analysis.getLossIncrProbDistFromCatalogsRandomFromCOV(analysis.getCatalogs(), 1d); +// UncertainArbDiscFunc incrTestPoiss = analysis.getLossIncrProbDistFromCatalogsRandomFromCOV(analysis.getRandomizedCatalogs(), 1d); +// ArrayList funcList = new ArrayList(); +// funcList.add(incrTestTD); +// funcList.add(incrTestPoiss); +// ArrayList plotChars2 = new ArrayList(); +// plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLUE)); +// plotChars2.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.RED)); +// plotChars2.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 2f, Color.BLUE)); +// plotChars2.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 2f, Color.RED)); +// String fileNamePrefix = null; +// PlottingUtils.writeAndOrPlotUncertFuncs(funcList, plotChars2, "Incr Test", "Loss ($1000)", "Probability", new Range(1e3,1e9), new Range(1e-6,1.0), true, true, fileNamePrefix, true); +// +// ArbitrarilyDiscretizedFunc exceedTestTD = getBlankLossCurve(0d); +// ArbitrarilyDiscretizedFunc exceedTestPois = getBlankLossCurve(0d); +// for(int i=0;i funcList2 = new ArrayList(); // AbstractDiscretizedFunc +// funcList2.add(exceedTestTD); +// funcList2.add(exceedTestPois); +// funcList2.add(lossExceedProbFromTD_CatalogsRandFromCOV); +// funcList2.add(lossExceedProbFromPoisCatalogsRandFromCOV); +// +// PlottingUtils.writeAndOrPlotFuncs(funcList2, plotChars2, "Exceed Test", "Loss ($1000)", "Probability", new Range(1e3,1e9), new Range(1e-6,1.0), true, true, fileNamePrefix, true); + + + + + // this plots the distribution of rupture losses as a function of magnitude (not including rup rates) +// analysis.plotRupLossVsMagStats(analysis.getCatalogs(), false, false); +// analysis.plotRupLossVsMagStats(analysis.getCatalogs(), true, false); + + // these are long term MFDs so COV not needed + // if I want these for the paper see scratch.ned.GK_Declustering.MakeFigures.makeFigure2_Parts(*) +// analysis.plotMFDs(); -// File lossDataCSV_Filename = new File("/Users/field/Library/CloudStorage/OneDrive-DOI/Field_Other/CEA_WGCEP/UCERF3/DeclusteringOnLossAnalysisPaper2024/Data/average_nth_rup_losses.csv"); -// CSVFile lossDataCSV_File = CSVFile.readFile(lossDataCSV_Filename, true); -// System.out.println("Loss Data:"); -// System.out.println("\tNumCols: "+lossDataCSV_File.getNumCols()); -// System.out.println("\tNumRows: "+lossDataCSV_File.getNumRows()); + // this assumes COV = 0 +// analysis.plotLossRateVsMag(); + + analysis.makeFigAndTableSetA_1yr_Exceedances(true); +// analysis.makeFigAndTableSetB_1yr_Exceedances(); + + +// // Plotting +// UncertainArbDiscFunc lossExceedProbFromTD_Catalogs = analysis.getLossExceedProbCurveFromCatalogs(analysis.getCatalogs(), 1d); +// lossExceedProbFromTD_Catalogs.setName("lossExceedProbFromTD_Catalogs"); +// UncertainArbDiscFunc lossExceedProbFromTD_CatalogsRandFromCOV = analysis.getLossExceedProbCurveFromCatalogsRandomFromCOV(analysis.getCatalogs(), 1d); +// lossExceedProbFromTD_CatalogsRandFromCOV.setName("lossExceedProbFromTD_CatalogsRandFromCOV"); +// UncertainArbDiscFunc aggrLossExceedProbFromTD_CatalogsRandFromCOV = analysis.getAggrLossExceedProbCurveFromCatalogsRandomFromCOV(analysis.getCatalogs(), 1d); +// aggrLossExceedProbFromTD_CatalogsRandFromCOV.setName("aggrLossExceedProbFromTD_CatalogsRandFromCOV"); +// UncertainArbDiscFunc aggrLossExceedProbFromTI_CatalogsRandFromCOV = analysis.getAggrLossExceedProbCurveFromCatalogsRandomFromCOV(analysis.getRandomizedCatalogs(random), 1d); +// aggrLossExceedProbFromTI_CatalogsRandFromCOV.setName("aggrLossExceedProbFromTI_CatalogsRandFromCOV"); +// ArrayList testFuncList = new ArrayList(); +// testFuncList.add(lossExceedProbFromTD_Catalogs); +// testFuncList.add(lossExceedProbFromTD_Catalogs); // for conf bounds +// testFuncList.add(lossExceedProbFromTD_CatalogsRandFromCOV); +// testFuncList.add(lossExceedProbFromTD_CatalogsRandFromCOV); // for conf bounds +// testFuncList.add(aggrLossExceedProbFromTD_CatalogsRandFromCOV); +// testFuncList.add(aggrLossExceedProbFromTD_CatalogsRandFromCOV); // for conf bounds +// testFuncList.add(aggrLossExceedProbFromTI_CatalogsRandFromCOV); +// testFuncList.add(aggrLossExceedProbFromTI_CatalogsRandFromCOV); // for conf bounds +// ArrayList plotChars = new ArrayList(); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLUE)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SHADED_UNCERTAIN_TRANS, 1f, Color.BLUE)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.RED)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SHADED_UNCERTAIN_TRANS, 1f, Color.RED)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLACK)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SHADED_UNCERTAIN_TRANS, 1f, Color.BLACK)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.ORANGE)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SHADED_UNCERTAIN_TRANS, 1f, Color.ORANGE)); +// PlottingUtils.writeAndOrPlotFuncs(testFuncList, plotChars, "test", "loss", "prob", null, null, false, true, null, true); + + + + + + // this demonstrates to be careful when log-plotting PDF x-axes (the distribution is distorted) +// testTakingLogOfPDF_Xaxis(); +// System.exit(0); + + // This was a test showing that LossCOV_Model was not giving the correct distribution. The two cases produced + // below are now equivalent because the bug was fixed. +// double mean = 1e5; +// LossCOV_Model covModel = LossCOV_Model.PORTER_POWER_LAW_2020_09_01; +// double cov = covModel.getCOV(mean); +// double samples[] = covModel.getDistribution(mean).sample(1000000); +// DescriptiveStatistics stats = new DescriptiveStatistics(samples); +// double meanFromDist = stats.getMean(); +// double covFromDist = stats.getStandardDeviation()/meanFromDist; +// System.out.println("mean="+mean+";\tmeanFromDist="+meanFromDist+";\tcov="+cov+";\tcovFromDist="+covFromDist+ +// "; meanRatio="+(float)(meanFromDist/mean)+"; covRatio="+(float)(covFromDist/cov)); +// double sigma = Math.sqrt(Math.log(cov*cov+1)); +// double mu = Math.log(mean)-(sigma*sigma/2); +// LogNormalDistribution dist = new LogNormalDistribution(mu, sigma); +// double samples2[] = dist.sample(1000000); +// DescriptiveStatistics stats2 = new DescriptiveStatistics(samples2); +// double meanFromDist2 = stats2.getMean(); +// double covFromDist2 = stats2.getStandardDeviation()/meanFromDist2; +// System.out.println("mean="+mean+";\tmeanFromDist2="+meanFromDist2+";\tcov="+cov+";\tcovFromDist2="+covFromDist2+ +// "; meanRatio2="+(float)(meanFromDist2/mean)+"; covRatio2="+(float)(covFromDist2/cov)); +// System.exit(0); + + + + + +// // This tests the covModel fix in Aug 2024; blue and black should be identical +// ArrayList testFuncArray = analysis.testCOV_ModelChangesInAug2024(); +// ArrayList plotChars = new ArrayList(); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLUE)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.RED)); +// plotChars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, 2f, Color.BLACK)); +// PlottingUtils.writeAndOrPlotFuncs(testFuncArray, plotChars, "test covModel changes", "loss", "prob", null, null, false, true, null, true); + + + +// System.exit(0); +// +// +// + +// // plot PDF in log10 space & compare to pure gaussian (this ignores spike at zero loss) +// ArbitrarilyDiscretizedFunc testCurve = analysis.getLossExceedProbCurveFromERF(); +// HistogramFunction testHist = convertExceedCurveToLog10_PDF_Hist(testCurve); +// testHist.setName("testHist"); +//// testHist.set(0,0.0); +// double mean = testHist.computeMean(); +// double stdDev = testHist.computeStdDev(); +// testHist.setInfo("Mean = "+mean+"\nStdDev = "+stdDev); +// org.opensha.commons.calc.GaussianDistCalc gaussDistCalc = new org.opensha.commons.calc.GaussianDistCalc(); +// EvenlyDiscretizedFunc testHistFit = testHist.deepClone(); +// for(int i=1;i testFuncList2 = new ArrayList(); +// testFuncList2.add(testHist); +// testFuncList2.add(testHistFit); +// analysis.quickPlot(testFuncList2, "Loss (thousand $)", "Prob", "PDF", false, false); + + + + + + +// // plot long-term loss rate curves +// ArbitrarilyDiscretizedFunc lossRateCurveFromERF = analysis.getLossRateCurveFromERF(); +// lossRateCurveFromERF.setName("lossRateCurveFromERF"); +// ArbitrarilyDiscretizedFunc lossRateCurveFromERF_ZeroCOV = analysis.getLossRateCurveFromERF_zeroCOV(); +// lossRateCurveFromERF_ZeroCOV.setName("lossRateCurveFromERF_ZeroCOV"); +// ArbitrarilyDiscretizedFunc lossRateCurveFromCatalogs = analysis.getLossRateCurveFromCatalogs(analysis.getCatalogs()); +// lossRateCurveFromCatalogs.setName("lossRateCurveFromCatalogs"); +// ArrayList funcs = new ArrayList(); +// funcs.add(lossRateCurveFromERF); +// funcs.add(lossRateCurveFromERF_ZeroCOV); +// funcs.add(lossRateCurveFromCatalogs); +// analysis.quickPlot(funcs, "Loss (thousand $)", "Rate (/yr)", "Loss Rate from ERF", true, true); + + +// // plot 1-yr loss exceed from ERF +// ArbitrarilyDiscretizedFunc lossExceedProbFromERF = analysis.getLossExceedProbCurveFromERF(); +// lossExceedProbFromERF.setName("lossExceedProbFromERF"); +// ArbitrarilyDiscretizedFunc testFunc = analysis.getBlankLossCurve(0.0); +// ArbitrarilyDiscretizedFunc lossExceedProbFromERF_ZeroCOV = analysis.getBlankLossCurve(0.0); +// ArbitrarilyDiscretizedFunc lossExceedProbFromCatalogs = analysis.getBlankLossCurve(0.0); +// +// for(int i=0;i funcProb = new ArrayList(); +// funcProb.add(lossExceedProbFromERF); +// funcProb.add(lossExceedProbFromERF_ZeroCOV); +// funcProb.add(lossExceedProbFromCatalogs); +// funcProb.add(testFunc); +// analysis.quickPlot(funcProb, "Loss (thousand $)", "Prob (/yr)", "Loss Exceed Curves", true, true); +// +// +// UncertainArbDiscFunc lossExceedCurveFromRandCats = analysis.getLossExceedProbCurveFromCatalogs(analysis.getRandomizedCatalogs(random), 1d); +// lossExceedCurveFromRandCats.setName("lossExceedCurveFromRandCats"); +// ArrayList funcProb2 = new ArrayList(); +// funcProb2.add(lossExceedProbFromCatalogs); +// funcProb2.add(lossExceedCurveFromRandCats); +// analysis.quickPlot(funcProb2, "Loss (thousand $)", "Prob (/yr)", "Test", true, true); + + + + + + +// analysis.old_writeAndOrPlotRateVersusTime(analysis.getCatalogs(), null, true, "Orig Catalogs"); +// analysis.old_writeAndOrPlotRateVersusTime(analysis.getRandomizedCatalogs(null), null, true, "Randomized Catalogs"); +// +// ArbitrarilyDiscretizedFunc lossRateCurveFromCatalogs = analysis.getLossRateCurveFromCatalogs(analysis.getCatalogs()); +// ArbitrarilyDiscretizedFunc lossExceedProbFromRateCatalogs = analysis.getBlankLossCurve(0.0); +// lossExceedProbFromRateCatalogs.setName("lossExceedProbFromRateCatalogs"); +// for(int i=0;i funcProb3 = new ArrayList(); +// funcProb3.add(lossExceedProbFromRateCatalogs); +// funcProb3.add(lossExceedProbFromRandCatalogs); +// funcProb3.add(lossExceedProbFromCatalogs); +// funcProb3.add(lossExceedProbFromTD_Catalogs); +// funcProb3.add(lossExceedProbFromGKdeclusteredCatalogs); +// analysis.quickPlot(funcProb3, "Loss (thousand $)", "Prob (/yr)", "Test", true, true); + + + + + +// // VERIFYING EQUIVALENCE OF VARIOUS TIME-IND APPROACHES FOR ZERO COV CASE (DURING DEBUGGIN) +// +// int numEvents = 0; +// for(ObsEqkRupList list:analysis.getCatalogs()) +// numEvents+=list.size(); +// System.out.println("Catalogs numEvents: "+numEvents); +// +// numEvents = 0; +// ArrayList randCatalogs = analysis.getRandomizedCatalogs(null); +// for(ObsEqkRupList list:randCatalogs) +// numEvents+=list.size(); +// System.out.println("randCatalogs numEvents: "+numEvents); +// ArrayList randSubCatalogs = analysis.getSubcatalogList(randCatalogs, 1.0); +// numEvents = 0; +// for(ObsEqkRupList list:randSubCatalogs) +// numEvents+=list.size(); +// System.out.println("randSubCatalogs numEvents: "+numEvents); +// +// +// // Any one of the following three work +// ArbitrarilyDiscretizedFunc lossRateCurveFromCatalogsZeroCOV = analysis.getLossRateCurveFromCatalogsZeroCOV(analysis.getCatalogs()); +//// ArbitrarilyDiscretizedFunc lossRateCurveFromCatalogsZeroCOV = analysis.getLossRateCurveFromCatalogsZeroCOV(randCatalogs); +//// ArbitrarilyDiscretizedFunc lossRateCurveFromCatalogsZeroCOV = analysis.getLossRateCurveFromCatalogsZeroCOV(randSubCatalogs); +//// lossRateCurveFromCatalogsZeroCOV.scale(500d); +// ArbitrarilyDiscretizedFunc lossExceedProbFromRateCatalogsZeroCOV = analysis.getBlankLossCurve(0.0); +// lossExceedProbFromRateCatalogsZeroCOV.setName("lossExceedProbFromRateCatalogsZeroCOV"); +// for(int i=0;i funcProb3 = new ArrayList(); +// funcProb3.add(lossExceedProbFromRateCatalogsZeroCOV); +// funcProb3.add(lossExceedProbFromCatalogsZeroCOV_RandCat); +//// funcProb3.add(lossExceedProbFromCatalogsZeroCOV); +// analysis.quickPlot(funcProb3, "Loss (thousand $)", "Prob (/yr)", "Test", true, true); - U3ETAS_LossSimulationAnalysis analysis = new U3ETAS_LossSimulationAnalysis(); + + + // REALLY OLD STUFF BELOW // Random random = new Random(102553864); // for reproducibility; change argument to get different results // analysis.getRandomizedCatalogs(random); // make and save the randomized catalogs 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 fc4a4171..f7992d0b 100644 --- a/src/main/java/scratch/ned/GK_Declustering/U3ETAS_SimulationAnalysis.java +++ b/src/main/java/scratch/ned/GK_Declustering/U3ETAS_SimulationAnalysis.java @@ -362,13 +362,12 @@ 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); - - // temporary hack - AbstractGridSourceProvider.SOURCE_MIN_MAG_CUTOFF = 2.55; - + FaultSystemSolutionERF_ETAS erf = ETAS_Launcher.buildERF(sol,false, 1d, 2012); erf.updateForecast();