From 8d8634c0b36818b1d35325763e047f98f2b5d474 Mon Sep 17 00:00:00 2001 From: Kevin Milner Date: Wed, 3 Jul 2024 11:12:04 -0700 Subject: [PATCH] better point source experiments --- ..._LogicTreeInversionRunnerScriptWriter.java | 754 +++++++++--------- .../HazardEquivalentRJBCorrPointSurface.java | 33 + .../HazardTargetedDistCorrCalc.java | 260 ++++++ .../InvCDF_RJBCorrPointSurface.java | 192 +++++ .../pointSources/InvCDF_RJBTableWriter.java | 234 ++++++ .../PointSourceHazardComparison.java | 525 ++++++++++-- .../pointSources/PointSurfaceBuilder.java | 357 ++++++++- .../TableBackedDistCorrPointSurface.java | 117 +++ .../scratch/kevin/ucerf3/PureScratch.java | 19 +- 9 files changed, 2031 insertions(+), 460 deletions(-) create mode 100644 src/main/java/scratch/kevin/pointSources/HazardEquivalentRJBCorrPointSurface.java create mode 100644 src/main/java/scratch/kevin/pointSources/HazardTargetedDistCorrCalc.java create mode 100644 src/main/java/scratch/kevin/pointSources/InvCDF_RJBCorrPointSurface.java create mode 100644 src/main/java/scratch/kevin/pointSources/InvCDF_RJBTableWriter.java create mode 100644 src/main/java/scratch/kevin/pointSources/TableBackedDistCorrPointSurface.java diff --git a/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java b/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java index 1ad9a8b6..2c8c1fc4 100644 --- a/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java +++ b/src/main/java/scratch/kevin/nshm23/MPJ_LogicTreeInversionRunnerScriptWriter.java @@ -246,388 +246,202 @@ public static void main(String[] args) throws IOException { * NSHM23 logic tree * TODO (this is a just a marker to find this part quickly, not an actual todo) */ -// List> levels = NSHM23_U3_HybridLogicTreeBranch.levels; -// dirName += "-nshm23_u3_hybrid_branches"; -// double avgNumRups = 325000; - - List> levels = NSHM23_LogicTreeBranch.levelsOnFault; - dirName += "-nshm23_branches"; - double avgNumRups = 600000; - -// List> levels = NSHM18_LogicTreeBranch.levels; -// dirName += "-nshm18_branches-wc_94"; -// double avgNumRups = 500000; - -// List> levels = NSHM18_LogicTreeBranch.levelsNewScale; -// dirName += "-nshm18_branches-new_scale"; -// double avgNumRups = 500000; - -// levels = new ArrayList<>(levels); -// for (int i=levels.size(); --i>=0;) -// if (levels.get(i).getType().isAssignableFrom(ShawSegmentationModels.class) -// || levels.get(i).getType().isAssignableFrom(NSHM23_SegmentationModels.class) -// || levels.get(i).getType().isAssignableFrom(SegmentationMFD_Adjustment.class) -// || levels.get(i).getType().isAssignableFrom(DistDependSegShift.class)) -// levels.remove(i); -// dirName += "-no_seg"; -//// levels.add(NSHM23_LogicTreeBranch.RUPS_THROUGH_CREEPING); -//// dirName += "-creep_branches"; -//// levels.add(NSHM23_LogicTreeBranch.MAX_DIST); -//// dirName += "-strict_cutoff_seg"; strictSeg = true; - - -// dirName += "-pre_zero_slip_parent_fix"; -// dirName += "-reweight_seg_2_3_4"; - -// levels = new ArrayList<>(levels); -// int origSize = levels.size(); -// for (int i=levels.size(); --i>=0;) -// if (levels.get(i).getType().isAssignableFrom(ScalingRelationships.class)) -// levels.remove(i); -// Preconditions.checkState(levels.size() < origSize); -// levels.add(NSHM23_LogicTreeBranch.SCALE); -// dirName += "-new_scale_rels"; -// dirName += "-full_set"; - - Class factoryClass = NSHM23_InvConfigFactory.class; - -// Class factoryClass = NSHM23_InvConfigFactory.MFDUncert0p1.class; -// dirName += "-mfd_uncert_0p1"; - -// Class factoryClass = NSHM23_InvConfigFactory.ConstantSlipRateStdDev0p1.class; -// dirName += "-const_slip_sd_0p1"; - -// Class factoryClass = NSHM23_InvConfigFactory.ConstantSlipRateStdDev0p2.class; -// dirName += "-const_slip_sd_0p2"; - -// Class factoryClass = NSHM23_InvConfigFactory.FullSysInv.class; -// dirName += "-full_sys_inv"; - -// Class factoryClass = NSHM23_InvConfigFactory.ClusterSpecific.class; -// dirName += "-cluster_specific_inversion"; - -// Class factoryClass = NSHM23_InvConfigFactory.SegWeight100.class; -// dirName += "-seg_weight_100"; - -// Class factoryClass = NSHM23_InvConfigFactory.SegWeight1000.class; -// dirName += "-seg_weight_1000"; - -// Class factoryClass = NSHM23_InvConfigFactory.SegWeight10000.class; -// dirName += "-seg_weight_10000"; - -// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevWeightAdjust.class; -// dirName += "-no_reweight_use_prev"; - -// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevWeightAdjustFullSys.class; -// dirName += "-full_sys_inv-no_reweight_use_prev"; - -// Class factoryClass = NSHM23_InvConfigFactory.HardcodedOrigWeights.class; -// dirName += "-no_reweight_use_orig"; - -// Class factoryClass = NSHM23_InvConfigFactory.HardcodedOrigWeightsFullSys.class; -// dirName += "-full_sys_inv-no_reweight_use_orig"; - -// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevAvgWeights.class; -// dirName += "-no_reweight_use_prev_avg"; - -// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevAvgWeightsFullSys.class; -// dirName += "-full_sys_inv-no_reweight_use_prev_avg"; - -// Class factoryClass = NSHM23_InvConfigFactory.NoPaleoParkfield.class; -// dirName += "-no_paleo_parkfield"; - -// Class factoryClass = NSHM23_InvConfigFactory.NoMFDScaleAdjust.class; -// dirName += "-no_scale_adj_mfds"; - -// Class factoryClass = NSHM23_InvConfigFactory.NoIncompatibleDataAdjust.class; -// dirName += "-no_mfd_sigma_data_adj"; - -// Class factoryClass = NSHM23_InvConfigFactory.ScaleLowerDepth1p3.class; -// dirName += "-scaleLowerDepth1.3"; - -// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevAsInitial.class; -// dirName += "-prev_as_initial"; - -// Class factoryClass = NSHM23_InvConfigFactory.NoAvg.class; -// dirName += "-no_avg"; - -// Class factoryClass = NSHM23_InvConfigFactory.ForceNewPaleo.class; -// dirName += "-new_paleo"; - -// Class factoryClass = NSHM23_InvConfigFactory.NewScaleUseOrigWidths.class; -// dirName += "-use_orig_widths"; - - // also set nonzero weights! -// Class factoryClass = NSHM23_InvConfigFactory.ForceWideSegBranches.class; -// dirName += "-wide_seg_branches"; - -// Class factoryClass = NSHM23_InvConfigFactory.ForceNoGhostTransient.class; -// dirName += "-no_ghost_trans"; - -// Class factoryClass = NSHM23_InvConfigFactory.ScaleSurfSlipUseActualWidths.class; -// dirName += "-surf_slip_use_actual_w"; - -// Class factoryClass = NSHM23_InvConfigFactory.RemoveIsolatedFaults.class; -// dirName += "-remove_isolated_faults"; - -// Class factoryClass = NSHM23_InvConfigFactory.RemoveProxyFaults.class; -// dirName += "-remove_proxy_faults"; - -// Class factoryClass = NSHM23_InvConfigFactory.NoPaleoSlip.class; -// dirName += "-no_paleo_slip"; - -// Class factoryClass = NSHM23_InvConfigFactory.PaleoSlipInequality.class; -// dirName += "-paleo_slip_ineq"; - -// Class factoryClass = NSHM23_InvConfigFactory.TenThousandItersPerRup.class; -// dirName += "-10000ip"; - -// Class factoryClass = NSHM23_InvConfigFactory.DM_OriginalWeights.class; -// dirName += "-dm_orig_weights"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; - -// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierlMinimizationWeights.class; -// dirName += "-dm_outlier_minimize_weights"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = false; - -// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierReplacementYc2p0.class; -// dirName += "-dm_outlier_sub_yc_2"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; - -// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierReplacementYc3p5.class; -// dirName += "-dm_outlier_sub_yc_3p5"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; - -// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierReplacementYc5p0.class; -// dirName += "-dm_outlier_sub_yc_5"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; - -// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierLogReplacementYc2p0.class; -// dirName += "-dm_outlier_log_sub_yc_2"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; - -// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierLogReplacementYc3p5.class; -// dirName += "-dm_outlier_log_sub_yc_3p5"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; - -// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierLogReplacementYc5p0.class; -// dirName += "-dm_outlier_log_sub_yc_5"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; - -// Class factoryClass = NSHM23_InvConfigFactory.SegModelLimitMaxLen.class; -// dirName += "-seg_limit_max_length"; - -// Class factoryClass = NSHM23_InvConfigFactory.SlipRateStdDevCeil0p1.class; -// dirName += "-slip_rate_sd_ceil_0p1"; - -// Class factoryClass = NSHM23_InvConfigFactory.SegModelMaxLen600.class; -// dirName += "-seg_limit_max_length_600"; - -// Class factoryClass = NSHM23_InvConfigFactory.SparseGRDontSpreadSingleToMulti.class; -// dirName += "-sparse_gr_dont_spread_single_multi"; - -// Class factoryClass = NSHM23_InvConfigFactory.ModDepthGV08.class; -// dirName += "-gv_08_mod_depth"; - -// Class factoryClass = NSHM23_InvConfigFactory.OrigDraftScaling.class; -// dirName += "-orig_draft_scaling"; - -// Class factoryClass = NSHM23_InvConfigFactory.ModScalingAdd4p3.class; -// dirName += "-mod_scaling"; - -// Class factoryClass = NSHM23_InvConfigFactory.NSHM18_UseU3Paleo.class; -// dirName += "-u3_paleo"; - -// Class factoryClass = NSHM23_InvConfigFactory.ModPitasPointDDW.class; -// dirName += "-mod_pitas_ddw"; - -// Class factoryClass = DefModSamplingEnabledInvConfig.ConnDistB0p5MidSegCorr.class; -// dirName += "-dm_sampling"; -// individualRandomLevels.add(new RandomDefModSampleLevel()); - -// Class factoryClass = DefModSamplingEnabledInvConfig.ConnDistB0p5MidSegCorrCapSigma.class; -// dirName += "-dm_sampling_cap_sigma"; -// individualRandomLevels.add(new RandomDefModSampleLevel()); - - if (!factoryClass.equals(NSHM23_InvConfigFactory.class)) { - // try instantiate it to make sure we get any static modifiers that might change branch weights - try { - System.out.println("Instantiating factory class: "+factoryClass.getName()); - factoryClass.getDeclaredConstructor().newInstance(); - } catch (Exception e) { - e.printStackTrace(); - } - } - -// levels = new ArrayList<>(levels); -// boolean randB = true; -// boolean randSeg = true; -// int origSize = levels.size(); -// for (int i=levels.size(); --i>=0;) { -// if (randB && SupraSeisBValues.class.isAssignableFrom(levels.get(i).getType())) -// levels.remove(i); -// if (randSeg && SegmentationModelBranchNode.class.isAssignableFrom(levels.get(i).getType())) -// levels.remove(i); -// } -// Preconditions.checkState(levels.size() < origSize); -// if (randB) { -// samplingBranchCountMultiplier *= 5; // there were originally 5 each -// dirName += "-randB"; -// individualRandomLevels.add(new RandomBValSampler.Level()); -// } -// if (randSeg) { -// samplingBranchCountMultiplier *= 5; // there were originally 5 each -// dirName += "-randSeg"; -// individualRandomLevels.add(new RandomSegModelSampler.Level()); -// } - -// dirName += "-mini_one_fifth"; -// samplingBranchCountMultiplier /= 5; - -// dirName += "-u3_perturb"; -// extraArgs.add("--perturb "+GenerationFunctionType.UNIFORM_0p001.name()); -// dirName += "-exp_perturb"; -// extraArgs.add("--perturb "+GenerationFunctionType.EXPONENTIAL_SCALE.name()); -// dirName += "-limit_zeros"; -// extraArgs.add("--non-negativity "+NonnegativityConstraintType.LIMIT_ZERO_RATES.name()); -// dirName += "-classic_sa"; -// extraArgs.add("--cooling-schedule "+CoolingScheduleType.CLASSICAL_SA.name()); - -// levels = new ArrayList<>(levels); -// levels.add(NSHM23_LogicTreeBranch.SINGLE_STATES); -// dirName += "-single_state"; - -// dirName += "-mod_west_valley_ddw"; - -// dirName += "-mod_dm_weights"; - - forceHazardGridSpacing = 0.1; - - forceRequiredNonzeroWeight = true; - LogicTreeNode[] required = { - // FAULT MODELS -// FaultModels.FM3_1, -// FaultModels.FM3_2, -// NSHM18_FaultModels.NSHM18_WUS_NoCA, -// NSHM18_FaultModels.NSHM18_WUS_PlusU3_FM_3p1, -// NSHM23_FaultModels.FM_v1p4, -// NSHM23_FaultModels.FM_v2, - NSHM23_FaultModels.WUS_FM_v3, -// PRVI25_FaultModels.PRVI_FM_INITIAL, - -// // SINGLE STATE -// NSHM23_SingleStates.NM, -// NSHM23_SingleStates.UT, - - // RUPTURE SETS -// RupturePlausibilityModels.COULOMB, // default -// RupturePlausibilityModels.COULOMB_5km, -// RupturePlausibilityModels.AZIMUTHAL, -// RupturePlausibilityModels.SEGMENTED, -// RupturePlausibilityModels.UCERF3, -// RupturePlausibilityModels.UCERF3_REDUCED, - - // DEFORMATION MODELS -// U3_UncertAddDeformationModels.U3_ZENG, -// U3_UncertAddDeformationModels.U3_MEAN, -// NSHM18_DeformationModels.BRANCH_AVERAGED, - NSHM23_DeformationModels.AVERAGE, -// NSHM23_DeformationModels.GEOLOGIC, -// NSHM23_DeformationModels.EVANS, -// NSHM23_DeformationModels.MEDIAN, - - // SCALING RELATIONSHIPS -// ScalingRelationships.SHAW_2009_MOD, -// ScalingRelationships.MEAN_UCERF3, -// NSHM23_ScalingRelationships.AVERAGE, -// NSHM23_ScalingRelationships.LOGA_C4p2_SQRT_LEN, -// NSHM23_ScalingRelationships.WIDTH_LIMITED_CSD, - - // SLIP ALONG RUPTURE -// NSHM23_SlipAlongRuptureModels.UNIFORM, -// NSHM23_SlipAlongRuptureModels.TAPERED, -// SlipAlongRuptureModels.UNIFORM, -// SlipAlongRuptureModels.TAPERED, - - // SUB-SECT CONSTRAINT -// SubSectConstraintModels.TOT_NUCL_RATE, // default -// SubSectConstraintModels.NUCL_MFD, - - // SUB-SEIS MO REDUCTION -// SubSeisMoRateReductions.SUB_B_1, -// SubSeisMoRateReductions.NONE, // default -// SubSeisMoRateReductions.SYSTEM_AVG, -// SubSeisMoRateReductions.SYSTEM_AVG_SUB_B_1, - - // SUPRA-SEIS-B -// SupraSeisBValues.B_0p5, - SupraSeisBValues.AVERAGE, - - // PALEO UNCERT -// NSHM23_PaleoUncertainties.EVEN_FIT, - - // SEGMENTATION -// SegmentationModels.SHAW_R0_3, - NSHM23_SegmentationModels.AVERAGE, -// NSHM23_SegmentationModels.MID, -// NSHM23_SegmentationModels.CLASSIC, -// NSHM23_SegmentationModels.CLASSIC_FULL, - - // SEG-SHIFT -// DistDependSegShift.NONE, -// DistDependSegShift.ONE_KM, -// DistDependSegShift.TWO_KM, -// DistDependSegShift.THREE_KM, - - // SEG ADJUSTMENT -// SegmentationMFD_Adjustment.NONE, -// SegmentationMFD_Adjustment.JUMP_PROB_THRESHOLD_AVG, -// SegmentationMFD_Adjustment.REL_GR_THRESHOLD_AVG_SINGLE_ITER, -// SegmentationMFD_Adjustment.REL_GR_THRESHOLD_AVG, // default -// SegmentationMFD_Adjustment.CAPPED_REDIST, -// SegmentationMFD_Adjustment.CAPPED_REDIST_SELF_CONTAINED, -// SegmentationMFD_Adjustment.GREEDY, -// SegmentationMFD_Adjustment.GREEDY_SELF_CONTAINED, -// SegmentationMFD_Adjustment.JUMP_PROB_THRESHOLD_AVG_MATCH_STRICT, - - // CREEPING SECTION -// RupsThroughCreepingSect.INCLUDE, -// RupsThroughCreepingSect.EXCLUDE, - }; -// LogicTreeNode[] required = { FaultModels.FM3_1, SubSeisMoRateReductionNode.SYSTEM_AVG }; -// LogicTreeNode[] required = { FaultModels.FM3_1, SubSeisMoRateReductionNode.FAULT_SPECIFIC }; -// Class sortBy = SubSectConstraintModels.class; - Class sortBy = NSHM23_SegmentationModels.class; - /* - * END NSHM23 logic tree - */ - - /* - * PRVI25 logic tree - * TODO (this is a just a marker to find this part quickly, not an actual todo) - */ -//// List> levels = PRVI25_LogicTreeBranch.levelsOnFault; -//// dirName += "-prvi25_crustal_branches"; -//// double avgNumRups = 50000; +//// List> levels = NSHM23_U3_HybridLogicTreeBranch.levels; +//// dirName += "-nshm23_u3_hybrid_branches"; +//// double avgNumRups = 325000; +// +// List> levels = NSHM23_LogicTreeBranch.levelsOnFault; +// dirName += "-nshm23_branches"; +// double avgNumRups = 600000; +// +//// List> levels = NSHM18_LogicTreeBranch.levels; +//// dirName += "-nshm18_branches-wc_94"; +//// double avgNumRups = 500000; +// +//// List> levels = NSHM18_LogicTreeBranch.levelsNewScale; +//// dirName += "-nshm18_branches-new_scale"; +//// double avgNumRups = 500000; // -// // random DM sampling //// levels = new ArrayList<>(levels); -//// int origNumLevels = levels.size(); //// for (int i=levels.size(); --i>=0;) -//// if (levels.get(i).getNodes().get(0) instanceof PRVI25_CrustalDeformationModels) +//// if (levels.get(i).getType().isAssignableFrom(ShawSegmentationModels.class) +//// || levels.get(i).getType().isAssignableFrom(NSHM23_SegmentationModels.class) +//// || levels.get(i).getType().isAssignableFrom(SegmentationMFD_Adjustment.class) +//// || levels.get(i).getType().isAssignableFrom(DistDependSegShift.class)) //// levels.remove(i); -//// Preconditions.checkState(levels.size() == origNumLevels -1); -//// individualRandomLevels.add(new PRVI25_CrustalRandomlySampledDeformationModelLevel()); -//// samplingBranchCountMultiplier = 5; // 5 for each branch -//// dirName += "-dmSample"; -//// if (samplingBranchCountMultiplier > 1) -//// dirName += samplingBranchCountMultiplier+"x"; +//// dirName += "-no_seg"; +////// levels.add(NSHM23_LogicTreeBranch.RUPS_THROUGH_CREEPING); +////// dirName += "-creep_branches"; +////// levels.add(NSHM23_LogicTreeBranch.MAX_DIST); +////// dirName += "-strict_cutoff_seg"; strictSeg = true; // -// List> levels = PRVI25_LogicTreeBranch.levelsSubduction; -// dirName += "-prvi25_subduction_branches"; -// double avgNumRups = 10000; -// gmpe = AttenRelRef.AG_2020_GLOBAL; +// +//// dirName += "-pre_zero_slip_parent_fix"; +//// dirName += "-reweight_seg_2_3_4"; // //// levels = new ArrayList<>(levels); -//// levels.add(NSHM23_LogicTreeBranch.SUB_SECT_CONSTR); +//// int origSize = levels.size(); +//// for (int i=levels.size(); --i>=0;) +//// if (levels.get(i).getType().isAssignableFrom(ScalingRelationships.class)) +//// levels.remove(i); +//// Preconditions.checkState(levels.size() < origSize); +//// levels.add(NSHM23_LogicTreeBranch.SCALE); +//// dirName += "-new_scale_rels"; +//// dirName += "-full_set"; +// +// Class factoryClass = NSHM23_InvConfigFactory.class; +// +//// Class factoryClass = NSHM23_InvConfigFactory.MFDUncert0p1.class; +//// dirName += "-mfd_uncert_0p1"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.ConstantSlipRateStdDev0p1.class; +//// dirName += "-const_slip_sd_0p1"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.ConstantSlipRateStdDev0p2.class; +//// dirName += "-const_slip_sd_0p2"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.FullSysInv.class; +//// dirName += "-full_sys_inv"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.ClusterSpecific.class; +//// dirName += "-cluster_specific_inversion"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.SegWeight100.class; +//// dirName += "-seg_weight_100"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.SegWeight1000.class; +//// dirName += "-seg_weight_1000"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.SegWeight10000.class; +//// dirName += "-seg_weight_10000"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevWeightAdjust.class; +//// dirName += "-no_reweight_use_prev"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevWeightAdjustFullSys.class; +//// dirName += "-full_sys_inv-no_reweight_use_prev"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.HardcodedOrigWeights.class; +//// dirName += "-no_reweight_use_orig"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.HardcodedOrigWeightsFullSys.class; +//// dirName += "-full_sys_inv-no_reweight_use_orig"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevAvgWeights.class; +//// dirName += "-no_reweight_use_prev_avg"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevAvgWeightsFullSys.class; +//// dirName += "-full_sys_inv-no_reweight_use_prev_avg"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.NoPaleoParkfield.class; +//// dirName += "-no_paleo_parkfield"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.NoMFDScaleAdjust.class; +//// dirName += "-no_scale_adj_mfds"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.NoIncompatibleDataAdjust.class; +//// dirName += "-no_mfd_sigma_data_adj"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.ScaleLowerDepth1p3.class; +//// dirName += "-scaleLowerDepth1.3"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.HardcodedPrevAsInitial.class; +//// dirName += "-prev_as_initial"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.NoAvg.class; +//// dirName += "-no_avg"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.ForceNewPaleo.class; +//// dirName += "-new_paleo"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.NewScaleUseOrigWidths.class; +//// dirName += "-use_orig_widths"; +// +// // also set nonzero weights! +//// Class factoryClass = NSHM23_InvConfigFactory.ForceWideSegBranches.class; +//// dirName += "-wide_seg_branches"; // -//// dirName += "-proxyGriddedTests"; +//// Class factoryClass = NSHM23_InvConfigFactory.ForceNoGhostTransient.class; +//// dirName += "-no_ghost_trans"; // -// Class factoryClass = PRVI25_InvConfigFactory.class; +//// Class factoryClass = NSHM23_InvConfigFactory.ScaleSurfSlipUseActualWidths.class; +//// dirName += "-surf_slip_use_actual_w"; // -// if (!factoryClass.equals(PRVI25_InvConfigFactory.class)) { +//// Class factoryClass = NSHM23_InvConfigFactory.RemoveIsolatedFaults.class; +//// dirName += "-remove_isolated_faults"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.RemoveProxyFaults.class; +//// dirName += "-remove_proxy_faults"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.NoPaleoSlip.class; +//// dirName += "-no_paleo_slip"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.PaleoSlipInequality.class; +//// dirName += "-paleo_slip_ineq"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.TenThousandItersPerRup.class; +//// dirName += "-10000ip"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.DM_OriginalWeights.class; +//// dirName += "-dm_orig_weights"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; +// +//// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierlMinimizationWeights.class; +//// dirName += "-dm_outlier_minimize_weights"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = false; +// +//// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierReplacementYc2p0.class; +//// dirName += "-dm_outlier_sub_yc_2"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; +// +//// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierReplacementYc3p5.class; +//// dirName += "-dm_outlier_sub_yc_3p5"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; +// +//// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierReplacementYc5p0.class; +//// dirName += "-dm_outlier_sub_yc_5"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; +// +//// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierLogReplacementYc2p0.class; +//// dirName += "-dm_outlier_log_sub_yc_2"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; +// +//// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierLogReplacementYc3p5.class; +//// dirName += "-dm_outlier_log_sub_yc_3p5"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; +// +//// Class factoryClass = NSHM23_InvConfigFactory.DM_OutlierLogReplacementYc5p0.class; +//// dirName += "-dm_outlier_log_sub_yc_5"; NSHM23_DeformationModels.ORIGINAL_WEIGHTS = true; +// +//// Class factoryClass = NSHM23_InvConfigFactory.SegModelLimitMaxLen.class; +//// dirName += "-seg_limit_max_length"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.SlipRateStdDevCeil0p1.class; +//// dirName += "-slip_rate_sd_ceil_0p1"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.SegModelMaxLen600.class; +//// dirName += "-seg_limit_max_length_600"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.SparseGRDontSpreadSingleToMulti.class; +//// dirName += "-sparse_gr_dont_spread_single_multi"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.ModDepthGV08.class; +//// dirName += "-gv_08_mod_depth"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.OrigDraftScaling.class; +//// dirName += "-orig_draft_scaling"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.ModScalingAdd4p3.class; +//// dirName += "-mod_scaling"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.NSHM18_UseU3Paleo.class; +//// dirName += "-u3_paleo"; +// +//// Class factoryClass = NSHM23_InvConfigFactory.ModPitasPointDDW.class; +//// dirName += "-mod_pitas_ddw"; +// +//// Class factoryClass = DefModSamplingEnabledInvConfig.ConnDistB0p5MidSegCorr.class; +//// dirName += "-dm_sampling"; +//// individualRandomLevels.add(new RandomDefModSampleLevel()); +// +//// Class factoryClass = DefModSamplingEnabledInvConfig.ConnDistB0p5MidSegCorrCapSigma.class; +//// dirName += "-dm_sampling_cap_sigma"; +//// individualRandomLevels.add(new RandomDefModSampleLevel()); +// +// if (!factoryClass.equals(NSHM23_InvConfigFactory.class)) { // // try instantiate it to make sure we get any static modifiers that might change branch weights // try { // System.out.println("Instantiating factory class: "+factoryClass.getName()); @@ -637,40 +451,226 @@ public static void main(String[] args) throws IOException { // } // } // +//// levels = new ArrayList<>(levels); +//// boolean randB = true; +//// boolean randSeg = true; +//// int origSize = levels.size(); +//// for (int i=levels.size(); --i>=0;) { +//// if (randB && SupraSeisBValues.class.isAssignableFrom(levels.get(i).getType())) +//// levels.remove(i); +//// if (randSeg && SegmentationModelBranchNode.class.isAssignableFrom(levels.get(i).getType())) +//// levels.remove(i); +//// } +//// Preconditions.checkState(levels.size() < origSize); +//// if (randB) { +//// samplingBranchCountMultiplier *= 5; // there were originally 5 each +//// dirName += "-randB"; +//// individualRandomLevels.add(new RandomBValSampler.Level()); +//// } +//// if (randSeg) { +//// samplingBranchCountMultiplier *= 5; // there were originally 5 each +//// dirName += "-randSeg"; +//// individualRandomLevels.add(new RandomSegModelSampler.Level()); +//// } +// +//// dirName += "-mini_one_fifth"; +//// samplingBranchCountMultiplier /= 5; +// +//// dirName += "-u3_perturb"; +//// extraArgs.add("--perturb "+GenerationFunctionType.UNIFORM_0p001.name()); +//// dirName += "-exp_perturb"; +//// extraArgs.add("--perturb "+GenerationFunctionType.EXPONENTIAL_SCALE.name()); +//// dirName += "-limit_zeros"; +//// extraArgs.add("--non-negativity "+NonnegativityConstraintType.LIMIT_ZERO_RATES.name()); +//// dirName += "-classic_sa"; +//// extraArgs.add("--cooling-schedule "+CoolingScheduleType.CLASSICAL_SA.name()); +// +//// levels = new ArrayList<>(levels); +//// levels.add(NSHM23_LogicTreeBranch.SINGLE_STATES); +//// dirName += "-single_state"; +// +//// dirName += "-mod_west_valley_ddw"; +// +//// dirName += "-mod_dm_weights"; +// // forceHazardGridSpacing = 0.1; // // forceRequiredNonzeroWeight = true; // LogicTreeNode[] required = { // // FAULT MODELS -//// PRVI25_CrustalFaultModels.PRVI_FM_INITIAL, -//// PRVI25_SubductionFaultModels.PRVI_SUB_FM_INITIAL, +//// FaultModels.FM3_1, +//// FaultModels.FM3_2, +//// NSHM18_FaultModels.NSHM18_WUS_NoCA, +//// NSHM18_FaultModels.NSHM18_WUS_PlusU3_FM_3p1, +//// NSHM23_FaultModels.FM_v1p4, +//// NSHM23_FaultModels.FM_v2, +// NSHM23_FaultModels.WUS_FM_v3, +//// PRVI25_FaultModels.PRVI_FM_INITIAL, +// +//// // SINGLE STATE +//// NSHM23_SingleStates.NM, +//// NSHM23_SingleStates.UT, // // // RUPTURE SETS //// RupturePlausibilityModels.COULOMB, // default //// RupturePlausibilityModels.COULOMB_5km, +//// RupturePlausibilityModels.AZIMUTHAL, +//// RupturePlausibilityModels.SEGMENTED, +//// RupturePlausibilityModels.UCERF3, +//// RupturePlausibilityModels.UCERF3_REDUCED, // // // DEFORMATION MODELS -//// PRVI25_CrustalDeformationModels.GEOLOGIC, +//// U3_UncertAddDeformationModels.U3_ZENG, +//// U3_UncertAddDeformationModels.U3_MEAN, +//// NSHM18_DeformationModels.BRANCH_AVERAGED, +// NSHM23_DeformationModels.AVERAGE, +//// NSHM23_DeformationModels.GEOLOGIC, +//// NSHM23_DeformationModels.EVANS, +//// NSHM23_DeformationModels.MEDIAN, // // // SCALING RELATIONSHIPS +//// ScalingRelationships.SHAW_2009_MOD, +//// ScalingRelationships.MEAN_UCERF3, +//// NSHM23_ScalingRelationships.AVERAGE, +//// NSHM23_ScalingRelationships.LOGA_C4p2_SQRT_LEN, +//// NSHM23_ScalingRelationships.WIDTH_LIMITED_CSD, +// +// // SLIP ALONG RUPTURE +//// NSHM23_SlipAlongRuptureModels.UNIFORM, +//// NSHM23_SlipAlongRuptureModels.TAPERED, +//// SlipAlongRuptureModels.UNIFORM, +//// SlipAlongRuptureModels.TAPERED, // // // SUB-SECT CONSTRAINT //// SubSectConstraintModels.TOT_NUCL_RATE, // default //// SubSectConstraintModels.NUCL_MFD, // +// // SUB-SEIS MO REDUCTION +//// SubSeisMoRateReductions.SUB_B_1, +//// SubSeisMoRateReductions.NONE, // default +//// SubSeisMoRateReductions.SYSTEM_AVG, +//// SubSeisMoRateReductions.SYSTEM_AVG_SUB_B_1, +// // // SUPRA-SEIS-B //// SupraSeisBValues.B_0p5, -//// SupraSeisBValues.AVERAGE, +// SupraSeisBValues.AVERAGE, +// +// // PALEO UNCERT +//// NSHM23_PaleoUncertainties.EVEN_FIT, // // // SEGMENTATION -//// NSHM23_SegmentationModels.AVERAGE, +//// SegmentationModels.SHAW_R0_3, +// NSHM23_SegmentationModels.AVERAGE, //// NSHM23_SegmentationModels.MID, //// NSHM23_SegmentationModels.CLASSIC, +//// NSHM23_SegmentationModels.CLASSIC_FULL, +// +// // SEG-SHIFT +//// DistDependSegShift.NONE, +//// DistDependSegShift.ONE_KM, +//// DistDependSegShift.TWO_KM, +//// DistDependSegShift.THREE_KM, +// +// // SEG ADJUSTMENT +//// SegmentationMFD_Adjustment.NONE, +//// SegmentationMFD_Adjustment.JUMP_PROB_THRESHOLD_AVG, +//// SegmentationMFD_Adjustment.REL_GR_THRESHOLD_AVG_SINGLE_ITER, +//// SegmentationMFD_Adjustment.REL_GR_THRESHOLD_AVG, // default +//// SegmentationMFD_Adjustment.CAPPED_REDIST, +//// SegmentationMFD_Adjustment.CAPPED_REDIST_SELF_CONTAINED, +//// SegmentationMFD_Adjustment.GREEDY, +//// SegmentationMFD_Adjustment.GREEDY_SELF_CONTAINED, +//// SegmentationMFD_Adjustment.JUMP_PROB_THRESHOLD_AVG_MATCH_STRICT, +// +// // CREEPING SECTION +//// RupsThroughCreepingSect.INCLUDE, +//// RupsThroughCreepingSect.EXCLUDE, // }; //// LogicTreeNode[] required = { FaultModels.FM3_1, SubSeisMoRateReductionNode.SYSTEM_AVG }; //// LogicTreeNode[] required = { FaultModels.FM3_1, SubSeisMoRateReductionNode.FAULT_SPECIFIC }; //// Class sortBy = SubSectConstraintModels.class; // Class sortBy = NSHM23_SegmentationModels.class; + /* + * END NSHM23 logic tree + */ + + /* + * PRVI25 logic tree + * TODO (this is a just a marker to find this part quickly, not an actual todo) + */ + List> levels = PRVI25_LogicTreeBranch.levelsOnFault; + dirName += "-prvi25_crustal_branches"; + double avgNumRups = 50000; + + // random DM sampling + levels = new ArrayList<>(levels); + int origNumLevels = levels.size(); + for (int i=levels.size(); --i>=0;) + if (levels.get(i).getNodes().get(0) instanceof PRVI25_CrustalDeformationModels) + levels.remove(i); + Preconditions.checkState(levels.size() == origNumLevels -1); + individualRandomLevels.add(new PRVI25_CrustalRandomlySampledDeformationModelLevel()); + samplingBranchCountMultiplier = 5; // 5 for each branch + dirName += "-dmSample"; + if (samplingBranchCountMultiplier > 1) + dirName += samplingBranchCountMultiplier+"x"; + +// List> levels = PRVI25_LogicTreeBranch.levelsSubduction; +// dirName += "-prvi25_subduction_branches"; +// double avgNumRups = 10000; +// gmpe = AttenRelRef.AG_2020_GLOBAL; + +// levels = new ArrayList<>(levels); +// levels.add(NSHM23_LogicTreeBranch.SUB_SECT_CONSTR); + +// dirName += "-proxyGriddedTests"; + + Class factoryClass = PRVI25_InvConfigFactory.class; + + if (!factoryClass.equals(PRVI25_InvConfigFactory.class)) { + // try instantiate it to make sure we get any static modifiers that might change branch weights + try { + System.out.println("Instantiating factory class: "+factoryClass.getName()); + factoryClass.getDeclaredConstructor().newInstance(); + } catch (Exception e) { + e.printStackTrace(); + } + } + + forceHazardGridSpacing = 0.1; + + forceRequiredNonzeroWeight = true; + LogicTreeNode[] required = { + // FAULT MODELS +// PRVI25_CrustalFaultModels.PRVI_FM_INITIAL, +// PRVI25_SubductionFaultModels.PRVI_SUB_FM_INITIAL, + + // RUPTURE SETS +// RupturePlausibilityModels.COULOMB, // default +// RupturePlausibilityModels.COULOMB_5km, + + // DEFORMATION MODELS +// PRVI25_CrustalDeformationModels.GEOLOGIC, + + // SCALING RELATIONSHIPS + + // SUB-SECT CONSTRAINT +// SubSectConstraintModels.TOT_NUCL_RATE, // default +// SubSectConstraintModels.NUCL_MFD, + + // SUPRA-SEIS-B +// SupraSeisBValues.B_0p5, +// SupraSeisBValues.AVERAGE, + + // SEGMENTATION +// NSHM23_SegmentationModels.AVERAGE, +// NSHM23_SegmentationModels.MID, +// NSHM23_SegmentationModels.CLASSIC, + }; +// LogicTreeNode[] required = { FaultModels.FM3_1, SubSeisMoRateReductionNode.SYSTEM_AVG }; +// LogicTreeNode[] required = { FaultModels.FM3_1, SubSeisMoRateReductionNode.FAULT_SPECIFIC }; +// Class sortBy = SubSectConstraintModels.class; + Class sortBy = NSHM23_SegmentationModels.class; /* * END PRVI25 logic tree */ diff --git a/src/main/java/scratch/kevin/pointSources/HazardEquivalentRJBCorrPointSurface.java b/src/main/java/scratch/kevin/pointSources/HazardEquivalentRJBCorrPointSurface.java new file mode 100644 index 00000000..ec99bae3 --- /dev/null +++ b/src/main/java/scratch/kevin/pointSources/HazardEquivalentRJBCorrPointSurface.java @@ -0,0 +1,33 @@ +package scratch.kevin.pointSources; + +import java.io.File; +import java.io.IOException; +import java.util.Random; + +import org.opensha.commons.data.CSVFile; +import org.opensha.commons.data.xyz.EvenlyDiscrXYZ_DataSet; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.LocationUtils; + +import com.google.common.base.Preconditions; + +public class HazardEquivalentRJBCorrPointSurface extends FiniteApproxPointSurface { + + private EvenlyDiscrXYZ_DataSet rJBtable; + private double mag; + private Location loc; + + public HazardEquivalentRJBCorrPointSurface(EvenlyDiscrXYZ_DataSet rJBtable, Location loc, double dip, double zTop, + double zBot, double length, double mag, boolean footwall) { + super(loc, dip, zTop, zBot, footwall, length); + this.rJBtable = rJBtable; + this.loc = loc; + this.mag = mag; + } + + @Override + public double getDistanceJB(Location siteLoc) { + return TableBackedDistCorrPointSurface.getDist(rJBtable, LocationUtils.horzDistanceFast(loc, siteLoc), mag); + } + +} diff --git a/src/main/java/scratch/kevin/pointSources/HazardTargetedDistCorrCalc.java b/src/main/java/scratch/kevin/pointSources/HazardTargetedDistCorrCalc.java new file mode 100644 index 00000000..cb2e659d --- /dev/null +++ b/src/main/java/scratch/kevin/pointSources/HazardTargetedDistCorrCalc.java @@ -0,0 +1,260 @@ +package scratch.kevin.pointSources; + +import java.io.File; +import java.io.IOException; +import java.text.DecimalFormat; +import java.util.ArrayDeque; +import java.util.ArrayList; +import java.util.Deque; +import java.util.List; +import java.util.Random; +import java.util.stream.IntStream; + +import org.opensha.commons.calc.GaussianDistCalc; +import org.opensha.commons.data.CSVFile; +import org.opensha.commons.data.Site; +import org.opensha.commons.data.function.DiscretizedFunc; +import org.opensha.commons.data.function.EvenlyDiscretizedFunc; +import org.opensha.commons.data.xyz.EvenlyDiscrXYZ_DataSet; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.LocationUtils; +import org.opensha.sha.earthquake.ProbEqkRupture; +import org.opensha.sha.earthquake.ProbEqkSource; +import org.opensha.sha.earthquake.faultSysSolution.util.SolHazardMapCalc.ReturnPeriods; +import org.opensha.sha.faultSurface.RuptureSurface; +import org.opensha.sha.imr.AttenRelRef; +import org.opensha.sha.imr.ScalarIMR; +import org.opensha.sha.imr.param.IntensityMeasureParams.PGA_Param; +import org.opensha.sha.imr.param.IntensityMeasureParams.SA_Param; +import org.opensha.sha.magdist.GutenbergRichterMagFreqDist; +import org.opensha.sha.magdist.IncrementalMagFreqDist; +import org.opensha.sha.util.FocalMech; + +import com.google.common.base.Preconditions; + +import scratch.kevin.pointSources.PointSourceHazardComparison.PointSourceCalcERF; +import scratch.kevin.pointSources.PointSourceHazardComparison.PointSourceType; + +public class HazardTargetedDistCorrCalc { + + public static void main(String[] args) throws IOException { + // calculate the rJB for each mag/horz dist that gives you equivalent ground motion exceedance probs using + // PointSourceNshm as you would doing it right + + // we'll do this for a single period + double period = 0d; + // and for this conditional probability of exceedance + double condExceedProb = 0.5; +// double condExceedProb = GaussianDistCalc.getExceedProb(1d); // +1 sigma + // and for this GMM + AttenRelRef gmmRef = AttenRelRef.NGAWest_2014_AVG_NOIDRISS; + + Location centerLoc = new Location(0d, 0d); + + GutenbergRichterMagFreqDist mfd = new GutenbergRichterMagFreqDist(5.05, 30, 0.1); + mfd.setAllButTotMoRate(mfd.getMinX(), mfd.getMaxX(), 1d, 1d); + + EvenlyDiscretizedFunc distFunc = new EvenlyDiscretizedFunc(0d, 301, 1d); + int numEachDist = 10; + List> distLocs = new ArrayList<>(); + for (int i=0; i locs = new ArrayList<>(numEachDist); + double deltaEach = 2d*Math.PI/numEachDist; + for (int j=0; j gmmDeque = new ArrayDeque<>(); + + for (FocalMech mech : mechs) { + System.out.println("Building ERFs for "+mech); + PointSourceCalcERF refERF = new PointSourceCalcERF(refType, centerLoc, mfd, + mech.rake(), mech.dip(), numPerRefType, rand); + PointSourceCalcERF targetERF = new PointSourceCalcERF(targetType, centerLoc, mfd, + mech.rake(), mech.dip(), numPerTargetType, rand); + + System.out.println("Calculating rJB as a function of center location for "+targetType); + EvenlyDiscrXYZ_DataSet targetRJBs = calcRJBs(targetERF, distFunc, distLocs, mfd); +// printTable(targetRJBs); + + System.out.println("Calculating IMs at condExceedProb="+(float)condExceedProb+" for "+refType); + EvenlyDiscrXYZ_DataSet refIMs = calcIMsAtExceedProb(refERF, distFunc, distLocs, mfd, gmmRef, period, gmmDeque, condExceedProb); +// printTable(refIMs); + System.out.println("Calculating IMs at condExceedProb="+(float)condExceedProb+" for "+targetType); + EvenlyDiscrXYZ_DataSet targetIMs = calcIMsAtExceedProb(targetERF, distFunc, distLocs, mfd, gmmRef, period, gmmDeque, condExceedProb); +// printTable(targetIMs); + + EvenlyDiscrXYZ_DataSet corrRJBs = new EvenlyDiscrXYZ_DataSet(distFunc.size(), mfd.size(), + distFunc.getMinX(), mfd.getMinX(), distFunc.getDelta(), mfd.getDelta()); + for (int m=0; m minIM, "maxIM=%s should be greater than minIM=%s; mag=%s", maxIM, minIM, mag); + for (int d=0; d maxIM) { + // this value is greater than any of our IMs, use our smallest rJB + corrRJB = minRJB; + } else if (refIM < minIM) { + // this value is less than any of our IMs, use our largest rJB + corrRJB = maxRJB; + } else { + // interpolate + double horzDistForIM = targetIMsFunc.getFirstInterpolatedX(refIM); + corrRJB = targetRJBs.bilinearInterpolation(horzDistForIM, mag); + } + corrRJBs.set(d, m, corrRJB); + } + } + + printTable(corrRJBs); + + File csvFile = new File(outputDir, mech.name()+"_equiv_rJBs.csv"); + System.out.println("Writing "+csvFile.getAbsolutePath()); + CSVFile csv = new CSVFile<>(true); + List header = new ArrayList<>(mfd.size()+1); + header.add("Horizontal Distance"); + for (int m=0; m line = new ArrayList<>(header.size()); + line.add((float)distFunc.getX(d)+""); + for (int m=0; m> distLocs, IncrementalMagFreqDist mfd) { + double[][] rJBs = new double[mfd.size()][distFunc.size()]; // y is outer index + + IntStream.range(0, distFunc.size()).parallel().forEach(d -> { + int[] numEntries = new int[mfd.size()]; + + for (Location loc : distLocs.get(d)) { + for (ProbEqkSource source : erf) { + for (ProbEqkRupture rup : source) { + double mag = rup.getMag(); + int m = mfd.getClosestXIndex(mag); + RuptureSurface surf = rup.getRuptureSurface(); + rJBs[m][d] += surf.getDistanceJB(loc); + numEntries[m]++; + } + } + } + + for (int m=0; m> distLocs, IncrementalMagFreqDist mfd, AttenRelRef gmmRef, double period, + Deque gmmDeque, double condExceedProb) { + double[][] ims = new double[mfd.size()][distFunc.size()]; // y is outer index + + List doneIndexes = new ArrayList<>(distFunc.size()); + IntStream.range(0, distFunc.size()).parallel().forEach(d -> { + int[] numEntries = new int[mfd.size()]; + + ScalarIMR gmm = null; + synchronized (gmmDeque) { + if (!gmmDeque.isEmpty()) + gmm = gmmDeque.pop(); + } + if (gmm == null) + gmm = gmmRef.get(); + if (period == 0d) { + gmm.setIntensityMeasure(PGA_Param.NAME); + } else { + Preconditions.checkState(period > 0d); + gmm.setIntensityMeasure(SA_Param.NAME); + SA_Param.setPeriodInSA_Param(gmm.getIntensityMeasure(), period); + } + + for (Location loc : distLocs.get(d)) { + Site site = new Site(loc); + site.addParameterList(gmm.getSiteParams()); + gmm.setSite(site); + for (ProbEqkSource source : erf) { + for (ProbEqkRupture rup : source) { + gmm.setEqkRupture(rup); + double mag = rup.getMag(); + int m = mfd.getClosestXIndex(mag); + ims[m][d] += gmm.getIML_AtExceedProb(condExceedProb); + numEntries[m]++; + } + } + } + + for (int m=0; m csv) { + // first figure out distance and magnitude spacing + double minDist = csv.getDouble(1, 0); + double maxDist = csv.getDouble(csv.getNumRows()-1, 0); + double distDelta = 0d; + for (int row=2; row (float)minDist) { + distDelta = dist - minDist; + break; + } + } + int numDist = (int)((maxDist - minDist)/distDelta + 0.5)+1; + distFunc = new EvenlyDiscretizedFunc(minDist, maxDist, numDist); + System.out.println("Detected dist range: ["+minDist+", "+maxDist+"]; size="+numDist); + + double minMag = csv.getDouble(1, 1); + double magDelta = csv.getDouble(2, 1) - minMag; + double maxMag = csv.getDouble(csv.getNumRows()-1, 1); + int numMag = (int)((maxMag - minMag)/magDelta + 0.5)+1; + magFunc = new EvenlyDiscretizedFunc(minMag, maxMag, numMag); + System.out.println("Detected mag range: ["+minMag+", "+maxMag+"]; size="+numMag); + + int expectedRows = numMag*numDist + 1; + Preconditions.checkState(expectedRows == csv.getNumRows(), + "Bad row count; expected %s, have %s", expectedRows, csv.getNumRows()); + + int numProbs = csv.getNumCols()-2; + EvenlyDiscretizedFunc probVals = new EvenlyDiscretizedFunc(0d, 1d, numProbs); + for (int i=0; i distFunc.getMaxX()) + dist = distFunc.getMaxX(); + if (mag < magFunc.getMinX()) + mag = magFunc.getMinX(); + else if (mag > magFunc.getMaxX()) + mag = magFunc.getMaxX(); + + int d = distFunc.getClosestXIndex(dist); + int m = magFunc.getClosestXIndex(mag); + double ret = invCmlPDFs[d][m].getInterpolatedY(prob); +// if (d == 10 && (float)mag == 7.15f) { +// synchronized (InvCDF_RJBCorrPointSurface.class) { +// System.out.println("Debug for "+(float)dist+" km, M"+(float)mag+", p="+(float)prob); +// System.out.println("\trJB["+(float)invCmlPDFs[d][m].getMinX()+"]: "+(float)invCmlPDFs[d][m].getY(0)); +// System.out.println("\trJB["+(float)prob+"]: "+(float)ret); +// System.out.println("\trJB["+(float)invCmlPDFs[d][m].getMaxX()+"]: "+(float)invCmlPDFs[d][m].getY(invCmlPDFs[d][m].size()-1)); +// } +// } + return ret; +// // this seems to be buggy, and would be slow anyway +// int distIndex1 = getXIndexBefore(distFunc, dist); +// Preconditions.checkState(distIndex1 >= 0); +// int distIndex2; +// if (distIndex1 == distFunc.size()-1 || (float)dist == (float)distFunc.getX(distIndex1)) +// distIndex2 = distIndex1; +// else +// distIndex2 = distIndex1+1; +// +// int magIndex1 = getXIndexBefore(magFunc, mag); +// Preconditions.checkState(magIndex1 >= 0); +// int magIndex2; +// if (magIndex1 == magFunc.size()-1 || (float)mag == (float)magFunc.getX(magIndex1)) +// magIndex2 = magIndex1; +// else +// magIndex2 = magIndex1+1; +// +// double val11 = invCmlPDFs[distIndex1][magIndex1].getInterpolatedY(prob); +// double val21; +// if (distIndex2 == distIndex1) +// val21 = val11; +// else +// val21 = invCmlPDFs[distIndex2][magIndex1].getInterpolatedY(prob); +// double val12, val22; +// if (magIndex2 == magIndex1) { +// val12 = val11; +// val22 = val21; +// } else { +// val12 = invCmlPDFs[distIndex1][magIndex2].getInterpolatedY(prob); +// if (distIndex2 == distIndex1) +// val22 = val12; +// else +// val22 = invCmlPDFs[distIndex2][magIndex2].getInterpolatedY(prob); +// } +// +// double distFrac = (dist - distFunc.getX(distIndex1))/distFunc.getDelta(); +// double magFrac = (mag - magFunc.getX(magIndex1))/magFunc.getDelta(); +// +// return (1 - magFrac) * ((1 - distFrac)*val11 + distFrac*val12) + +// magFrac * ((1 - distFrac)*val21 + distFrac*val22); + } + + protected int getXIndexBefore(EvenlyDiscretizedFunc func, double x) { + return (int)Math.floor((x-func.getMinX())/func.getDelta()); + } + } + private RJBCorrInvPDFs invCmlPDFs; + private Location loc; + private double mag; + private double rJB_prob; + + public InvCDF_RJBCorrPointSurface(RJBCorrInvPDFs invCmlPDFs, Location loc, double dip, double zTop, double zBot, + boolean footwall, double length, double mag, double rJB_prob) { + super(loc, dip, zTop, zBot, footwall, length); + this.invCmlPDFs = invCmlPDFs; + this.loc = loc; + this.mag = mag; + this.rJB_prob = rJB_prob; + } + + @Override + public double getDistanceJB(Location siteLoc) { + double dist = LocationUtils.horzDistanceFast(loc, siteLoc); + return invCmlPDFs.calcRJB(dist, mag, rJB_prob); + } + + public static DiscretizedFunc getEvenlySpacedProbs(int num) { + EvenlyDiscretizedFunc edgeFunc = new EvenlyDiscretizedFunc(0d, 1d, num+1); + EvenlyDiscretizedFunc centerFunc = new EvenlyDiscretizedFunc(0.5*edgeFunc.getDelta(), 1d-0.5*edgeFunc.getDelta(), num); + for (int i=0; i 0) + edgeFunc.set(1d-cdf, 0d); + } + edgeFunc.set(1d, 0d); + ArbitrarilyDiscretizedFunc func = new ArbitrarilyDiscretizedFunc(); + // now calculate weights + for (int i=1; i> distLocs = new ArrayList<>(); + for (int i=0; i locs = new ArrayList<>(numEachDist); + double deltaEach = 2d*Math.PI/numEachDist; + for (int j=0; j csv = new CSVFile<>(true); + List header = new ArrayList<>(mfd.size()+1); + header.add("Horizontal Distance"); + header.add("Magnitude"); + for (int i=0; i line = new ArrayList<>(header.size()); + line.add((float)distFunc.getX(d)+""); + line.add((float)mfd.getX(m)+""); + for (int i=0; i> distLocs, IncrementalMagFreqDist mfd, int numDistPts) { + EvenlyDiscretizedFunc[][] ret = new EvenlyDiscretizedFunc[distFunc.size()][mfd.size()]; + + System.out.println("Calculating for "+distLocs.size()+" distances"); + List doneIndexes = new ArrayList<>(); +// IntStream.range(0, distFunc.size()).forEach(d -> { + IntStream.range(0, distFunc.size()).parallel().forEach(d -> { + List> magRJBLists = new ArrayList<>(mfd.size()); + for (int m=0; m()); + + for (Location loc : distLocs.get(d)) { + for (ProbEqkSource source : erf) { + for (ProbEqkRupture rup : source) { + double mag = rup.getMag(); + int m = mfd.getClosestXIndex(mag); + RuptureSurface surf = rup.getRuptureSurface(); + magRJBLists.get(m).add(surf.getDistanceJB(loc)); + } + } + } + + for (int m=0; m rJBList = magRJBLists.get(m); + Preconditions.checkState(rJBList.size() > 1); + + LightFixedXFunc normCDF = ArbDiscrEmpiricalDistFunc.calcQuickNormCDF(rJBList, null); +// System.out.println("NormCDF for horzDist="+(float)distFunc.getX(d)+", M"+(float)mfd.getX(m) +// +": "+rJBList.size()+" rJBs, normCDF.size()="+normCDF.size()); +// System.out.println("\tFirst point: "+normCDF.get(0)); +// System.out.println("\tLast point: "+normCDF.get(normCDF.size()-1)); + + double minDist = normCDF.getX(0); + double minProb = normCDF.getY(0); + double maxDist = normCDF.getX(normCDF.size()-1); + double maxProb = normCDF.getY(normCDF.size()-1); + + EvenlyDiscretizedFunc invCDF = new EvenlyDiscretizedFunc(0d, 1d, numDistPts); + for (int i=0; i= maxProb) + invCDF.set(i, maxDist); + else + invCDF.set(i, normCDF.getFirstInterpolatedX(prob)); + } + ret[d][m] = invCDF; + } + synchronized (doneIndexes) { + System.out.print("."); + doneIndexes.add(d); + if (doneIndexes.size() % 100 == 0) + System.out.println(" "+doneIndexes.size()); + } + }); + + if (doneIndexes.size() % 100 != 0) + System.out.println(" "+doneIndexes.size()); + + return ret; + } + +} diff --git a/src/main/java/scratch/kevin/pointSources/PointSourceHazardComparison.java b/src/main/java/scratch/kevin/pointSources/PointSourceHazardComparison.java index 678fff38..1c04b1fa 100644 --- a/src/main/java/scratch/kevin/pointSources/PointSourceHazardComparison.java +++ b/src/main/java/scratch/kevin/pointSources/PointSourceHazardComparison.java @@ -29,12 +29,14 @@ import org.jfree.chart.ui.TextAnchor; import org.jfree.data.Range; import org.opensha.commons.calc.magScalingRelations.magScalingRelImpl.WC1994_MagLengthRelationship; +import org.opensha.commons.data.CSVFile; import org.opensha.commons.data.Site; import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc; import org.opensha.commons.data.function.DefaultXY_DataSet; import org.opensha.commons.data.function.DiscretizedFunc; import org.opensha.commons.data.function.EvenlyDiscretizedFunc; import org.opensha.commons.data.function.XY_DataSet; +import org.opensha.commons.data.xyz.EvenlyDiscrXYZ_DataSet; import org.opensha.commons.data.xyz.GriddedGeoDataSet; import org.opensha.commons.geo.GriddedRegion; import org.opensha.commons.geo.Location; @@ -86,6 +88,11 @@ import com.google.common.base.Preconditions; import com.google.common.base.Stopwatch; +import com.google.common.collect.HashBasedTable; +import com.google.common.collect.Table; + +import scratch.kevin.pointSources.InvCDF_RJBCorrPointSurface.RJBCorrInvPDFs; +import scratch.kevin.pointSources.TableBackedDistCorrPointSurface.DistanceCorrTables; public class PointSourceHazardComparison { @@ -165,67 +172,134 @@ public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, return new RupListSource(rups, source0); } }, + SIMPLE_APPROX_FINITE_POINT_NO_CORR("Simple Approx. Finite, No Dist. Corr.", true, Type.NONE) { + @Override + public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, + double aveDip, boolean isSupersample, Random r) { + return buildPointSource(centerLoc, null, mfd, aveRake, aveDip, r, getDistCorrType()); + } + }, + SIMPLE_APPROX_FINITE_POINT_NSHMP_CORR("Simple Approx. Finite, NSHMP08 Dist. Corr.", true, Type.NSHMP08) { + @Override + public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, + double aveDip, boolean isSupersample, Random r) { + return buildPointSource(centerLoc, null, mfd, aveRake, aveDip, r, getDistCorrType()); + } + }, SINGLE_QUAD_TRACE_CENTERED("Single Quad Surface, Trace Centered", true, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, null, mfd, aveRake, aveDip, r, true, 1); + PointSurfaceBuilder builder = new PointSurfaceBuilder(centerLoc); + builder.dip(aveDip); + builder.fractionalHypocentralDepth(0d); + builder.random(r); + return buildQuadSource(builder, mfd, aveRake, aveDip, 1); } }, SINGLE_QUAD_SURF_CENTER("Single Quad Surface", true, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, null, mfd, aveRake, aveDip, r, false, 1); + PointSurfaceBuilder builder = new PointSurfaceBuilder(centerLoc); + builder.dip(aveDip); + builder.fractionalHypocentralDepth(0.5); + builder.random(r); + return buildQuadSource(builder, mfd, aveRake, aveDip, 1); } }, CROSSHAIR_QUAD("Crosshair Quad Surface", true, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, null, mfd, aveRake, aveDip, r, false, 2); + PointSurfaceBuilder builder = new PointSurfaceBuilder(centerLoc); + builder.dip(aveDip); + builder.fractionalHypocentralDepth(0.5d); + builder.random(r); + return buildQuadSource(builder, mfd, aveRake, aveDip, 2); } }, - TRIPLE_QUAD("3x Quad Surface", true, Type.NONE) { + OCT_QUAD("8x Quad Surface", true, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, null, mfd, aveRake, aveDip, r, false, 3); + PointSurfaceBuilder builder = new PointSurfaceBuilder(centerLoc); + builder.dip(aveDip); + builder.fractionalHypocentralDepth(0.5d); + builder.random(r); + return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 8); } }, - QUAD_QUAD("4x Quad Surface", true, Type.NONE) { + OCT_QUAD_RAND_DAS_DD("8x Quad Surface, Random DAS & DD", true, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, null, mfd, aveRake, aveDip, r, false, 4); + PointSurfaceBuilder builder = new PointSurfaceBuilder(centerLoc); + builder.dip(aveDip); + builder.sampleDASs().sampleHypocentralDepths(); + builder.random(r); + return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 8); } }, - QUAD_QUAD_RAND_CELL("4x Quad Surface, Random Cell Locations", true, Type.NONE) { + OCT_QUAD_RAND_CELL("8x Quad Surface, Random Cell Locations", true, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, isSupersample ? null : cellRegion(centerLoc), mfd, aveRake, aveDip, r, false, 4); + PointSurfaceBuilder builder = new PointSurfaceBuilder(centerLoc); + builder.dip(aveDip); + builder.fractionalHypocentralDepth(0.5d); + if (!isSupersample) + builder.sampleFromCell(cellRegion(centerLoc)); + builder.random(r); + return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 8); } }, - SEXT_QUAD_RAND_CELL("6x Quad Surface, Random Cell Locations", true, Type.NONE) { + OCT_QUAD_RAND_DAS_DD_CELL("8x Quad Surface, Random DAS, DD, & Cell Locations", true, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, isSupersample ? null : cellRegion(centerLoc), mfd, aveRake, aveDip, r, false, 6); + PointSurfaceBuilder builder = new PointSurfaceBuilder(centerLoc); + builder.dip(aveDip); + builder.sampleDASs().sampleHypocentralDepths(); + if (!isSupersample) + builder.sampleFromCell(cellRegion(centerLoc)); + builder.random(r); + return buildQuadSource(builder, mfd, aveRake, aveDip, isSupersample ? 2 : 8); } }, - OCT_QUAD_RAND_CELL("8x Quad Surface, Random Cell Locations", true, Type.NONE) { + TABLE_FINITE_APPROX_POINT_SOURCE("Table-Based Finite Approx", false, Type.NONE) { + @Override + public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, + double aveDip, boolean isSupersample, Random r) { + return buildTableBackedSource(centerLoc, mfd, aveRake, aveDip, false, r); + } + }, + TABLE_FINITE_APPROX_POINT_SOURCE_SUPERSAMPLE("Table-Based Finite Supersample Approx", false, Type.NONE) { + @Override + public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, + double aveDip, boolean isSupersample, Random r) { + return buildTableBackedSource(centerLoc, mfd, aveRake, aveDip, isSupersample ? false : true, r); + } + }, + HAZ_EQUIV_APPROX_POINT_SOURCE("Hazard Equivalent Finite Approx", false, Type.NONE) { + @Override + public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, + double aveDip, boolean isSupersample, Random r) { + return buildHazEquivSource(centerLoc, mfd, aveRake, aveDip, this); + } + }, + CDF_BASED_RJB_CORR_APPROX_POINT_SOURCE("RJB Distribution Approx Finite Point Source", false, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, isSupersample ? null : cellRegion(centerLoc), mfd, aveRake, aveDip, r, false, 8); + return buildRJBCorrSource(centerLoc, mfd, aveRake, aveDip, OCT_QUAD, 5); } }, - DODEC_QUAD_RAND_CELL("12x Quad Surface, Random Cell Locations", true, Type.NONE) { + CDF_BASED_RJB_CORR_APPROX_SUPERSAMPLED_POINT_SOURCE("RJB Distribution Approx Supersampled Finite Point Source", false, Type.NONE) { @Override public ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r) { - return buildQuadSource(centerLoc, isSupersample ? null : cellRegion(centerLoc), mfd, aveRake, aveDip, r, false, 12); + return buildRJBCorrSource(centerLoc, mfd, aveRake, aveDip, isSupersample ? OCT_QUAD : OCT_QUAD_RAND_CELL, 5); } }; @@ -240,8 +314,17 @@ private PointSourceType(String label, boolean stochastic, Type distCorrType) { this.distCorrType = distCorrType; } + public Type getDistCorrType() { + return distCorrType; + } + public abstract ProbEqkSource buildSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, double aveDip, boolean isSupersample, Random r); + + @Override + public String toString() { + return label; + } } private static Map mechWtMapForRake(double rake) { @@ -277,15 +360,39 @@ private static Region cellRegion(Location loc) { return new Region(new Location(loc.lat-half, loc.lon-half), new Location(loc.lat+half, loc.lon+half)); } - private static ProbEqkSource buildQuadSource(Location centerLoc, Region cell, IncrementalMagFreqDist mfd, - double aveRake, double aveDip, Random r, boolean traceCentered, int numEach) { + private static ProbEqkSource buildQuadSource(PointSurfaceBuilder builder, IncrementalMagFreqDist mfd, + double aveRake, double aveDip, int numEach) { + double dipRad = Math.toRadians(aveDip); + List rups = new ArrayList<>(mfd.size()*numEach); + for (int i=0; i rups = new ArrayList<>(mfd.size()*numEach); + boolean[] footwalls = aveDip == 90 ? new boolean[] {true} : new boolean[] {false, true}; + List rups = new ArrayList<>(mfd.size()*footwalls.length); for (int i=0; i rups = new ArrayList<>(); + ProbEqkSource source0 = null; + for (int i=0; i refHazEquivRJBs = HashBasedTable.create(); + private static Map refHazEquivDirs = new EnumMap<>(PointSourceType.class); + static { + refHazEquivDirs.put(PointSourceType.HAZ_EQUIV_APPROX_POINT_SOURCE, + new File(REF_MAIN_DIR, "SIMPLE_APPROX_FINITE_POINT_NO_CORR_vs_OCT_QUAD/hazard_equiv_rJBs")); + } + + public static ProbEqkSource buildHazEquivSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, + double aveDip, PointSourceType type) { + EvenlyDiscrXYZ_DataSet distCorr = getHazEquivRJBTable(type, aveDip); + + double dipRad = Math.toRadians(aveDip); + List rups = new ArrayList<>(); + ProbEqkSource source0 = null; + for (int i=0; i refRJBCorrInvPDFs = HashBasedTable.create(); + + public static ProbEqkSource buildRJBCorrSource(Location centerLoc, IncrementalMagFreqDist mfd, double aveRake, + double aveDip, PointSourceType type, int numProbs) { + RJBCorrInvPDFs corrTable = getRJBCorrTable(type, aveDip); + + DiscretizedFunc cdfProbs = InvCDF_RJBCorrPointSurface.getEvenlySpacedProbs(numProbs); +// DiscretizedFunc cdfProbs = InvCDF_RJBCorrPointSurface.getEvenlySpacedProbs(501); +// DiscretizedFunc cdfProbs = InvCDF_RJBCorrPointSurface.getSigmaSpacedProbs(5); + + boolean[] footwalls = aveDip == 90d ? new boolean[] {true} : new boolean[] {false,true}; + + double dipRad = Math.toRadians(aveDip); + List rups = new ArrayList<>(); + ProbEqkSource source0 = null; + for (int i=0; i csv = CSVFile.readFile(csvFile, true); + ret = new RJBCorrInvPDFs(csv); + } catch (IOException e) { + throw ExceptionUtils.asRuntimeException(e); + } + refRJBCorrInvPDFs.put(type, aveDip, ret); + return ret; + } + private static class RupListSource extends ProbEqkSource { private ProbEqkSource source0; @@ -392,7 +688,7 @@ public ProbEqkRupture getRupture(int nRupture) { } - private static class PointSourceCalcERF extends AbstractERF { + static class PointSourceCalcERF extends AbstractERF { private Type distCorrType; private List sources; @@ -628,17 +924,27 @@ private RupSurfMapProps(String label) { private static final double SUPERSAMPLE_SPACING = 0.01; + // TODO: this is just a marker to make it easy to find the top of main public static void main(String[] args) throws IOException { // PointSourceType mainType = PointSourceType.POINT_SOURCE_13b_NSHMP_CORR; PointSourceType mainType = PointSourceType.POINT_SOURCE_NSHM; // PointSourceType mainType = PointSourceType.QUAD_QUAD; // PointSourceType mainType = PointSourceType.QUAD_QUAD_RAND_CELL; // PointSourceType mainType = PointSourceType.OCT_QUAD_RAND_CELL; -// PointSourceType mainType = PointSourceType.DODEC_QUAD_RAND_CELL; +// PointSourceType mainType = PointSourceType.OCT_QUAD; // PointSourceType mainType = PointSourceType.APROX_SUPERSAMPLE_POINT_SOURCE_NSHM; // PointSourceType mainType = PointSourceType.POINT_TO_FINITE; // PointSourceType compType = null; - PointSourceType compType = PointSourceType.OCT_QUAD_RAND_CELL; +// PointSourceType compType = PointSourceType.TABLE_FINITE_APPROX_POINT_SOURCE; +// PointSourceType compType = PointSourceType.HAZ_EQUIV_APPROX_POINT_SOURCE; +// PointSourceType compType = PointSourceType.CDF_BASED_RJB_CORR_APPROX_POINT_SOURCE; +// PointSourceType compType = PointSourceType.CDF_BASED_RJB_CORR_APPROX_SUPERSAMPLED_POINT_SOURCE; +// PointSourceType compType = PointSourceType.SIMPLE_APPROX_FINITE_POINT_NSHMP_CORR; +// PointSourceType compType = PointSourceType.SIMPLE_APPROX_FINITE_POINT_NO_CORR; +// PointSourceType compType = PointSourceType.TABLE_FINITE_APPROX_POINT_SOURCE_SUPERSAMPLE; +// PointSourceType compType = PointSourceType.OCT_QUAD; +// PointSourceType compType = PointSourceType.OCT_QUAD_RAND_CELL; + PointSourceType compType = PointSourceType.OCT_QUAD_RAND_DAS_DD; // PointSourceType compType = PointSourceType.POINT_SOURCE_NSHM; // PointSourceType compType = PointSourceType.POINT_SOURCE_13b_NSHMP_CORR; // PointSourceType compType = PointSourceType.POINT_SOURCE_13b_NO_CORR; @@ -646,11 +952,18 @@ public static void main(String[] args) throws IOException { // PointSourceType compType = PointSourceType.POINT_TO_SIMPLE_QUAD; // PointSourceType compType = PointSourceType.QUAD_QUAD_RAND_CELL; +// int sleepMins = 60*6; +// System.out.println("Sleeping for "+sleepMins+" mins"); +// try { +// Thread.sleep(sleepMins*60l*1000l); +// } catch (InterruptedException e) {} + boolean doSingleCellHazard = true; boolean doNSHMModelHazard = true; boolean doHighRes = false; boolean doSupersample = true; - boolean forceWriteIntermediate = true; + boolean forceWriteIntermediate = false; + boolean writeTables = false; AttenRelRef gmmRef = AttenRelRef.NGAWest_2014_AVG_NOIDRISS; double[] periods = {0d, 1d}; @@ -688,21 +1001,31 @@ public static void main(String[] args) throws IOException { EnumMap siteGrids = new EnumMap<>(GridType.class); siteGrids.put(GridType.COLOCATED, colocatedGridded); - EvenlyDiscretizedFunc distFunc = new EvenlyDiscretizedFunc(0d, 100d, 100); - List> distFuncLocs = new ArrayList<>(); + EvenlyDiscretizedFunc tableDistFunc = writeTables ? + new EvenlyDiscretizedFunc(0d, 300d, 301) : new EvenlyDiscretizedFunc(0d, 100d, 101); + EvenlyDiscretizedFunc distFunc = new EvenlyDiscretizedFunc(0d, 100d, 101); + Range distPlotRange = new Range(0d, 100d); + List> tableDistFuncLocs = new ArrayList<>(); int numEachDist = 10; - for (int i=0; i locs = new ArrayList<>(numEachDist); double deltaEach = 2d*Math.PI/numEachDist; for (int j=0; j> distFuncLocs = tableDistFuncLocs.subList(0, distFunc.size()); + System.out.println("Maximum relocated distance was "+(float)maxActualDist+", max intended was "+(float)tableDistFunc.getMaxX()); GriddedRegion offsetCalcGridded = new GriddedRegion(calcRegion, spacing, new Location(halfSpacing, halfSpacing)); siteGrids.put(GridType.OFFSET, offsetCalcGridded); @@ -919,7 +1242,6 @@ public static void main(String[] args) throws IOException { table.addLine("DistanceJB", "DistanceRup"); for (double mag : distPlotMags) { - Range distRange = new Range(0d, distFunc.getMaxX()); Range ratioRange = new Range(0, 2); Range diffRange = new Range(-30d, 10d); @@ -1010,7 +1332,7 @@ public static void main(String[] args) throws IOException { HeadlessGraphPanel gp = PlotUtils.initHeadless(); - gp.drawGraphPanel(List.of(ratioSpec, diffSpec), false, false, List.of(distRange), List.of(ratioRange, diffRange)); + gp.drawGraphPanel(List.of(ratioSpec, diffSpec), false, false, List.of(distPlotRange), List.of(ratioRange, diffRange)); gp.setRenderingOrder(DatasetRenderingOrder.REVERSE); @@ -1157,6 +1479,13 @@ public static void main(String[] args) throws IOException { if (writeIntermediate) writeMarkdown(outputDir, lines, tocIndex); + if (writeTables) { + System.out.println("Calculating mag/dist tables"); + writeMagDistFuncs(centerERF, tableDistFunc, tableDistFuncLocs, mfd, resourcesDir, mechPrefix+"_centered"); + if (doSupersample) + writeMagDistFuncs(supersampledERF, tableDistFunc, tableDistFuncLocs, mfd, resourcesDir, mechPrefix+"_supersampled"); + } + lines.add("### "+mechLabel+" Surface Properties"); lines.add(topLink); lines.add(""); @@ -1303,7 +1632,7 @@ public static void main(String[] args) throws IOException { if (compType != null) table.addColumn("![Map]("+resourcesDir.getName()+"/"+gridPrefix+"_supersampled_comp.png)"); table.finalizeLine(); - } else if (compType != null) { + } else if (compType == null) { table.finalizeLine(); } } @@ -1488,7 +1817,7 @@ public static void main(String[] args) throws IOException { "Distance (km)", rps[r].label+", "+perLabels[p]+" ("+perUnits[p]+")"); spec.setLegendInset(RectangleAnchor.TOP_RIGHT); - gp.drawGraphPanel(spec, false, true, new Range(0d, xy.getMaxX()), logBufferedYRange(funcs)); + gp.drawGraphPanel(spec, false, true, distPlotRange, logBufferedYRange(funcs)); gp.setRenderingOrder(DatasetRenderingOrder.REVERSE); @@ -1989,8 +2318,7 @@ private static List calcCurves(List siteLocs, PointSo private static List calcAverageCurves(List> siteLocs, PointSourceCalcERF erf, AttenRelRef gmmRef, Deque gmmDeque, Deque calcDeque, - double period, DiscretizedFunc xVals, - DiscretizedFunc logXVals, ExecutorService exec) { + double period, DiscretizedFunc xVals, DiscretizedFunc logXVals, ExecutorService exec) { List>> futures = new ArrayList<>(siteLocs.size()); for (List siteLocBundle : siteLocs) { @@ -2198,17 +2526,20 @@ private static void calcDistFuncsAtMag(PointSourceCalcERF erf, List rupsAtMag = getRupsAtMag(erf, mag); IntStream.range(0, distLocs.size()).parallel().forEach(i -> { - List rJBs = new ArrayList<>(rupsAtMag.size()); - List rRups = new ArrayList<>(rupsAtMag.size()); + double sumWeight = 0d; + double sumWeightedRJB = 0d; + double sumWeightedRRup = 0d; for (Location loc : distLocs.get(i)) { for (ProbEqkRupture rup : rupsAtMag) { + double weight = rup.getMeanAnnualRate(1); RuptureSurface surf = rup.getRuptureSurface(); - rJBs.add(surf.getDistanceJB(loc)); - rRups.add(surf.getDistanceRup(loc)); + sumWeight += weight; + sumWeightedRJB += weight*surf.getDistanceJB(loc); + sumWeightedRRup += weight*surf.getDistanceRup(loc); } } - rJB.set(i, rJBs.stream().mapToDouble(D->D).average().getAsDouble()); - rRup.set(i, rRups.stream().mapToDouble(D->D).average().getAsDouble()); + rJB.set(i, sumWeightedRJB/sumWeight); + rRup.set(i, sumWeightedRRup/sumWeight); }); } @@ -2216,18 +2547,20 @@ private static void calcMagFuncsAtDist(PointSourceCalcERF erf, List di EvenlyDiscretizedFunc rJB, EvenlyDiscretizedFunc rRup, double dist) { IntStream.range(0, rJB.size()).parallel().forEach(i -> { List rupsAtMag = getRupsAtMag(erf, rJB.getX(i)); - - List rJBs = new ArrayList<>(rupsAtMag.size()); - List rRups = new ArrayList<>(rupsAtMag.size()); - for (Location distLoc : distLocs) { + double sumWeight = 0d; + double sumWeightedRJB = 0d; + double sumWeightedRRup = 0d; + for (Location loc : distLocs) { for (ProbEqkRupture rup : rupsAtMag) { + double weight = rup.getMeanAnnualRate(1); RuptureSurface surf = rup.getRuptureSurface(); - rJBs.add(surf.getDistanceJB(distLoc)); - rRups.add(surf.getDistanceRup(distLoc)); + sumWeight += weight; + sumWeightedRJB += weight*surf.getDistanceJB(loc); + sumWeightedRRup += weight*surf.getDistanceRup(loc); } } - rJB.set(i, rJBs.stream().mapToDouble(D->D).average().getAsDouble()); - rRup.set(i, rRups.stream().mapToDouble(D->D).average().getAsDouble()); + rJB.set(i, sumWeightedRJB/sumWeight); + rRup.set(i, sumWeightedRRup/sumWeight); }); } @@ -2236,10 +2569,14 @@ private static EvenlyDiscretizedFunc calcMagProps(PointSourceCalcERF erf, Increm EvenlyDiscretizedFunc ret = new EvenlyDiscretizedFunc(mfd.getMinX(), mfd.size(), mfd.getDelta()); IntStream.range(0, mfd.size()).parallel().forEach(i -> { List rupsAtMag = getRupsAtMag(erf, mfd.getX(i)); - List props = new ArrayList<>(rupsAtMag.size()); - for (ProbEqkRupture rup : rupsAtMag) - props.add(prop.calcForRupture(rup)); - ret.set(i, props.stream().mapToDouble(D->D).average().getAsDouble()); + double sumWeight = 0d; + double sumWeightedProps = 0d; + for (ProbEqkRupture rup : rupsAtMag) { + double weight = rup.getMeanAnnualRate(1); + sumWeight += weight; + sumWeightedProps += weight*prop.calcForRupture(rup); + } + ret.set(i, sumWeightedProps/sumWeight); }); return ret; } @@ -2249,11 +2586,16 @@ private static GriddedGeoDataSet calcMapProps(PointSourceCalcERF erf, GriddedReg GriddedGeoDataSet xyz = new GriddedGeoDataSet(gridReg); IntStream.range(0, gridReg.getNodeCount()).parallel().forEach(i -> { Location loc = gridReg.getLocation(i); - List props = new ArrayList<>(); - for (ProbEqkSource source : erf) - for (ProbEqkRupture rup : source) - props.add(prop.calcForRuptureAndLocation(rup, loc)); - xyz.set(i, props.stream().mapToDouble(D->D).average().getAsDouble()); + double sumWeight = 0d; + double sumWeightedProps = 0d; + for (ProbEqkSource source : erf) { + for (ProbEqkRupture rup : source) { + double weight = rup.getMeanAnnualRate(1); + sumWeight += weight; + sumWeightedProps += weight*prop.calcForRuptureAndLocation(rup, loc); + } + } + xyz.set(i, sumWeightedProps/sumWeight); }); return xyz; } @@ -2275,5 +2617,68 @@ private static EvenlyDiscretizedFunc asDistDiff(EvenlyDiscretizedFunc distFunc) return ret; } + + private static void writeMagDistFuncs(PointSourceCalcERF erf, EvenlyDiscretizedFunc distFunc, + List> distLocs, IncrementalMagFreqDist mfd, File outputDir, String prefix) throws IOException { + System.out.println("Calculating mag/dist funcs, will write with prefix "+prefix); + Preconditions.checkState(distLocs.size() == distFunc.size()); + double[][] rJBs = new double[distFunc.size()][mfd.size()]; + double[][] rRups = new double[distFunc.size()][mfd.size()]; + double[][] rXs = new double[distFunc.size()][mfd.size()]; + double[][] rXFractPos = new double[distFunc.size()][mfd.size()]; + + IntStream.range(0, distFunc.size()).parallel().forEach(d -> { + int[] numEntries = new int[mfd.size()]; + + for (Location loc : distLocs.get(d)) { + for (ProbEqkSource source : erf) { + for (ProbEqkRupture rup : source) { + double mag = rup.getMag(); + int m = mfd.getClosestXIndex(mag); + RuptureSurface surf = rup.getRuptureSurface(); + rJBs[d][m] += surf.getDistanceJB(loc); + rRups[d][m] += surf.getDistanceRup(loc); + double rX = surf.getDistanceX(loc); + rXs[d][m] += Math.abs(rX); + if (rX >= 0d) + rXFractPos[d][m] += 1d; + numEntries[m]++; + } + } + } + + for (int m=0; m csv = new CSVFile<>(true); + List header = new ArrayList<>(mfd.size()+1); + header.add("Horizontal Distance"); + for (int m=0; m line = new ArrayList<>(header.size()); + line.add((float)distFunc.getX(d)+""); + for (int m=0; m seeds = new ArrayList<>(); + seeds.add(Double.doubleToLongBits(loc.lat)); + seeds.add(Double.doubleToLongBits(loc.lon)); + seeds.add(Double.doubleToLongBits(loc.depth)); + seeds.add(Double.doubleToLongBits(zTop)); + seeds.add(Double.doubleToLongBits(zBot)); + if (sampleFromCell != null) + seeds.add((long)sampleFromCell.hashCode()); + if (Double.isFinite(mag)) + seeds.add(Double.doubleToLongBits(mag)); + if (Double.isFinite(strike)) + seeds.add(Double.doubleToLongBits(strike)); + seeds.add(Double.doubleToLongBits(dip)); + if (Double.isFinite(length)) + seeds.add(Double.doubleToLongBits(length)); + if (footwall) + seeds.add(1l); + if (Double.isFinite(zHyp)) + seeds.add(Double.doubleToLongBits(zHyp)); + if (Double.isFinite(zHypFract)) + seeds.add(Double.doubleToLongBits(zHypFract)); + if (scale != WC94 && scale != null && scale.getName() != null) + seeds.add((long)scale.getName().hashCode()); + if (Double.isFinite(gridSpacing)) + seeds.add(Double.doubleToLongBits(gridSpacing)); + rand = new Random(uniqueSeedCombination(seeds)); + } return rand; } + /** + * Generates a 64-bit random seed that is a repeatable and psuedo-unique combination of the input seeds. + * + * This is based on {@link Arrays#hashCode(int[])}, but modified for longs. + * + * @param seeds + * @return + */ + private static long uniqueSeedCombination(List seeds) { + Preconditions.checkState(!seeds.isEmpty()); + + long result = 1; + for (long element : seeds) + result = 31l * result + element; + return result; + } + private Location getLoc() { if (sampleFromCell != null) { double minLat = sampleFromCell.getMinLat(); @@ -140,6 +215,254 @@ public PointSurfaceBuilder lowerDepth(double zBot) { return this; } + /** + * Sets the fractional hypocentral depth to this value, such that {@code zHyp = zTop + zHypFract*(zBot-zTop)}. Values of + * zero are trace-centered (the trace is coincident with the cell location), and values of 0.5 are surface-centered. + * + * Setting this will clear any already-set fixed hypocentral depth value (see {@link #hypocentralDepth(double)} or + * distribution. + * @param zHyp + * @return + */ + public PointSurfaceBuilder fractionalHypocentralDepth(double zHypFract) { + Preconditions.checkArgument(zHypFract >= 0d && zHypFract <= 1d, "zHypFract=%s is not within the range [0,1]", zHypFract); + this.zHypFract = zHypFract; + this.zHyp = Double.NaN; + this.zHypSample = false; + this.zHypCDF = null; + this.zHypFractCDF = null; + return this; + } + + /** + * Sets the hypocentral depth to this value in km. If this value is outside the range of the upper and lower depths, + * an exception will be thrown when finite ruptures are built. + * + * Setting this will clear any already-set fractional hypocentral depth value (see {@link #fractionalHypocentralDepth(double)} + * or distribution. + * @param zHyp + * @return + */ + public PointSurfaceBuilder hypocentralDepth(double zHyp) { + FaultUtils.assertValidDepth(zHyp); + this.zHyp = zHyp; + this.zHypFract = Double.NaN; + this.zHypSample = false; + this.zHypCDF = null; + this.zHypFractCDF = null; + return this; + } + + /** + * If set, hypocentral depths will be sampled from a uniform distribution between zTop and zBot. + * + * Setting this will clear any already-set fractional or absolute hypocentral depth values or distributions. + * @return + */ + public PointSurfaceBuilder sampleHypocentralDepths() { + this.zHypFract = Double.NaN; + this.zHyp = Double.NaN; + this.zHypSample = true; + this.zHypCDF = null; + this.zHypFractCDF = null; + return this; + } + + /** + * If set, hypocentral depths will be sampled from this distribution of fractional values. + * + * Setting this will clear any already-set fractional or absolute hypocentral depth values or distributions. + * @param zHypFractDistribution distribution of fractional zHyp values; y values must sum to 1, and x values must + * not exceed the range [0,1] + * @return + */ + public PointSurfaceBuilder sampleFractionalHypocentralDepths(EvenlyDiscretizedFunc zHypFractDistribution) { + Preconditions.checkNotNull(zHypFractDistribution, "Cannot set the distribution to null. Clear it by setting " + + "zHyp another way."); + Preconditions.checkState((float)zHypFractDistribution.getMinX() >= 0f && (float)zHypFractDistribution.getMaxX() <= 1f, + "Distribution sum of x values must be in the range [0,1]: [%s,%s]", + (float)zHypFractDistribution.getMinX(), (float)zHypFractDistribution.getMaxX()); + this.zHypFract = Double.NaN; + this.zHyp = Double.NaN; + this.zHypSample = true; + this.zHypCDF = null; + this.zHypFractCDF = toCDF(zHypFractDistribution); + return this; + } + + /** + * If set, hypocentral depths will be sampled from this distribution. If this distribution contains values is + * outside the range of the upper and lower depths, an exception will be thrown when finite ruptures are built. + * + * Setting this will clear any already-set fractional or absolute hypocentral depth values or distributions. + * @param zHypFractDistribution distribution of absolute zHyp values; y values must sum to 1. + * @return + */ + public PointSurfaceBuilder sampleHypocentralDepths(EvenlyDiscretizedFunc zHypDistribution) { + Preconditions.checkNotNull(zHypDistribution, "Cannot set the distribution to null. Clear it by setting " + + "zHyp another way."); + this.zHypFract = Double.NaN; + this.zHyp = Double.NaN; + this.zHypSample = true; + this.zHypCDF = toCDF(zHypDistribution); + this.zHypFractCDF = null; + return this; + } + + /** + * Sets the distance along strike to this value, such that {@code das = dasFract*length}. Values of 0.5 are surface-centered. + * + * Setting this will clear any already-set fixed DAS value (see {@link #hypocentralDepth(double)} or + * distribution. + * @param zHyp + * @return + */ + public PointSurfaceBuilder fractionalDAS(double dasFract) { + Preconditions.checkArgument(dasFract >= 0d && dasFract <= 1d, "dasFract=%s is not within the range [0,1]", dasFract); + this.das = Double.NaN; + this.dasFract = dasFract; + this.dasSample = false; + this.dasCDF = null; + this.dasFractCDF = null; + return this; + } + + /** + * Sets the distance along strike to this value in km. If this value is outside the range of [0, length], + * an exception will be thrown when finite ruptures are built. + * + * Setting this will clear any already-set fractional DAS value (see {@link #fractionalHypocentralDepth(double)} + * or distribution. + * @param zHyp + * @return + */ + public PointSurfaceBuilder das(double das) { + Preconditions.checkState(das >= 0d); + this.das = das; + this.dasFract = Double.NaN; + this.dasSample = false; + this.dasCDF = null; + this.dasFractCDF = null; + return this; + } + + /** + * If set, distances along strike will be sampled from a uniform distribution between 0 and length. + * + * Setting this will clear any already-set fractional or absolute DAS values or distributions. + * @return + */ + public PointSurfaceBuilder sampleDASs() { + this.dasFract = Double.NaN; + this.das = Double.NaN; + this.dasSample = true; + this.dasCDF = null; + this.dasFractCDF = null; + return this; + } + + /** + * If set, distances along strike will be sampled from this distribution of fractional values. + * + * Setting this will clear any already-set fractional or absolute DAS values or distributions. + * @param dasFractDistribution distribution of fractional DAS values; y values must sum to 1, and x values must + * not exceed the range [0,1] + * @return + */ + public PointSurfaceBuilder sampleFractionalDASs(EvenlyDiscretizedFunc dasFractDistribution) { + Preconditions.checkNotNull(dasFractDistribution, "Cannot set the distribution to null. Clear it by setting " + + "DAS another way."); + Preconditions.checkState((float)dasFractDistribution.getMinX() >= 0f && (float)dasFractDistribution.getMaxX() <= 1f, + "Distribution sum of x values must be in the range [0,1]: [%s,%s]", + (float)dasFractDistribution.getMinX(), (float)dasFractDistribution.getMaxX()); + this.zHypFract = Double.NaN; + this.zHyp = Double.NaN; + this.zHypSample = true; + this.zHypCDF = null; + this.zHypFractCDF = toCDF(dasFractDistribution); + return this; + } + + /** + * If set, distances along strike will be sampled from this distribution. If this distribution contains values is + * outside the range of the [0, length], an exception will be thrown when finite ruptures are built. + * + * Setting this will clear any already-set fractional or absolute DAS values or distributions. + * @param dasDistribution distribution of fractional DAS values; y values must sum to 1. + * @return + */ + public PointSurfaceBuilder sampleDASs(EvenlyDiscretizedFunc dasDistribution) { + Preconditions.checkNotNull(dasDistribution, "Cannot set the distribution to null. Clear it by setting " + + "DAS another way."); + this.dasFract = Double.NaN; + this.das = Double.NaN; + this.dasSample = true; + this.dasCDF = toCDF(dasDistribution); + this.dasFractCDF = null; + return this; + } + + private static EvenlyDiscretizedFunc toCDF(EvenlyDiscretizedFunc dist) { + Preconditions.checkState((float)dist.calcSumOfY_Vals() == 1f, "Distribution sum of y values must " + + "sum to 1: %s", (float)dist.calcSumOfY_Vals()); + EvenlyDiscretizedFunc cdf = new EvenlyDiscretizedFunc(dist.getMinX(), dist.getMaxX(), dist.size()); + double sumY = 0d; + for (int i=0; i= (float)lowerBound && (float)absValue <= (float)upperBound, + "Supplied absolute value %s is not in the allowed range: [%s, %s]", + (float)absValue, (float)lowerBound, (float)upperBound); + return (absValue-lowerBound)/(upperBound-lowerBound); + } + Preconditions.checkState(sample); + double randDouble = getRand().nextDouble(); + if (fractDist != null) + return sampleFromCDF(fractDist, randDouble); + if (absDist != null) { + absValue = sampleFromCDF(absDist, randDouble); + Preconditions.checkState((float)absValue >= (float)lowerBound && (float)absValue <= (float)upperBound, + "Supplied absolute value %s is not in the allowed range: [%s, %s]", + (float)absValue, (float)lowerBound, (float)upperBound); + return (absValue-lowerBound)/(upperBound-lowerBound); + } + // uniform distribution + return randDouble; + } + /** * Sets the magnitude, used to infer the length if the length is not explicitly set * @param mag @@ -202,17 +525,6 @@ public PointSurfaceBuilder gridSpacing(double gridSpacing) { return this; } - /** - * Sets if the trace should be centered on the grid node (true), or if the middle of the surface should be (false). - * Only affects dipping ruptures. - * @param traceCentered - * @return - */ - public PointSurfaceBuilder traceCentered(boolean traceCentered) { - this.traceCentered = traceCentered; - return this; - } - private double getCalcLength() { if (Double.isFinite(length)) return length; @@ -265,18 +577,19 @@ private FaultTrace buildTrace(double strike) { double length = getCalcLength(); Preconditions.checkState(length > 0, "Can't build finite surface because length=%s; " + "set magnitude to infer length from scaling relationship", length); - double halfLen = 0.5*length; + double dasFract = getFractionalValue(0d, length, this.dasFract, das, dasSample, dasFractCDF, dasCDF); double strikeRad = Math.toRadians(strike); Location loc = getLoc(); - Location l0 = LocationUtils.location(loc, strikeRad-Math.PI, halfLen); - Location l1 = LocationUtils.location(loc, strikeRad, halfLen); - if (!traceCentered && zBot > zTop && dip < 90) { - // translate it so that the surface is centered rather than the trace + Location l0 = LocationUtils.location(loc, strikeRad-Math.PI, length*dasFract); + Location l1 = LocationUtils.location(loc, strikeRad, length*(1d-dasFract)); + if (zBot > zTop && dip < 90) { + // translate it for the given zHyp + double horzFract = getFractionalValue(zTop, zBot, zHypFract, zHyp, zHypSample, zHypFractCDF, zHypCDF); double horzWidth = (zBot-zTop)/Math.tan(Math.toRadians(dip)); - // move to the left (so that it dips to the right) + // move to the left (so that it follows the RHR and dips to the right) double transAz = strikeRad - 0.5*Math.PI; - l0 = LocationUtils.location(l0, transAz, 0.5*horzWidth); - l1 = LocationUtils.location(l1, transAz, 0.5*horzWidth); + l0 = LocationUtils.location(l0, transAz, horzFract*horzWidth); + l1 = LocationUtils.location(l1, transAz, horzFract*horzWidth); } l0 = new Location(l0.lat, l0.lon, zTop); l1 = new Location(l1.lat, l1.lon, zTop); diff --git a/src/main/java/scratch/kevin/pointSources/TableBackedDistCorrPointSurface.java b/src/main/java/scratch/kevin/pointSources/TableBackedDistCorrPointSurface.java new file mode 100644 index 00000000..f45ba434 --- /dev/null +++ b/src/main/java/scratch/kevin/pointSources/TableBackedDistCorrPointSurface.java @@ -0,0 +1,117 @@ +package scratch.kevin.pointSources; + +import java.io.File; +import java.io.IOException; +import java.util.Random; + +import org.opensha.commons.data.CSVFile; +import org.opensha.commons.data.xyz.EvenlyDiscrXYZ_DataSet; +import org.opensha.commons.geo.Location; +import org.opensha.commons.geo.LocationUtils; + +import com.google.common.base.Preconditions; + +public class TableBackedDistCorrPointSurface extends FiniteApproxPointSurface { + + public static class DistanceCorrTables { + public final EvenlyDiscrXYZ_DataSet rJB; + public final EvenlyDiscrXYZ_DataSet rRup; + public final EvenlyDiscrXYZ_DataSet rX; + public final EvenlyDiscrXYZ_DataSet rXPositive; + private DistanceCorrTables(EvenlyDiscrXYZ_DataSet rJB, EvenlyDiscrXYZ_DataSet rRup, EvenlyDiscrXYZ_DataSet rX, + EvenlyDiscrXYZ_DataSet rXPositive) { + super(); + this.rJB = rJB; + this.rRup = rRup; + this.rX = rX; + this.rXPositive = rXPositive; + } + } + + static DistanceCorrTables loadCorrTables(File dir, String prefix) throws IOException { + EvenlyDiscrXYZ_DataSet rJB = loadXYZ_CSV(new File(dir, prefix+"_rJB.csv")); + EvenlyDiscrXYZ_DataSet rRup = loadXYZ_CSV(new File(dir, prefix+"_rRup.csv")); + EvenlyDiscrXYZ_DataSet rX = loadXYZ_CSV(new File(dir, prefix+"_rX.csv")); + EvenlyDiscrXYZ_DataSet rXPositive = loadXYZ_CSV(new File(dir, prefix+"_rX_positive.csv")); + return new DistanceCorrTables(rJB, rRup, rX, rXPositive); + } + + static EvenlyDiscrXYZ_DataSet loadXYZ_CSV(File csvFile) throws IOException { + CSVFile csv = CSVFile.readFile(csvFile, true); + double minMag = csv.getDouble(0, 1); + double magSpacing = csv.getDouble(0, 2) - minMag; + int numMag = csv.getNumCols()-1; + + double minDist = csv.getDouble(1, 0); + double distSpacing = csv.getDouble(2, 0) - minDist; + int numDist = csv.getNumRows()-1; + + EvenlyDiscrXYZ_DataSet xyz = new EvenlyDiscrXYZ_DataSet(numDist, numMag, minDist, minMag, distSpacing, magSpacing); + for (int d=0; d 0 || (float)xyz.getY(m) == (float)csv.getDouble(0, m+1)); + xyz.set(d, m, csv.getDouble(d+1, m+1)); + } + } + return xyz; + } + + private DistanceCorrTables corrTables; + private double mag; + private Location loc; + + private Random rand; + + public TableBackedDistCorrPointSurface(DistanceCorrTables corrTables, Location loc, double dip, double zTop, + double zBot, double length, double mag, Random rand) { + super(loc, dip, zTop, zBot, true, length); + this.corrTables = corrTables; + this.loc = loc; + this.mag = mag; + this.rand = rand; + } + + static double getDist(EvenlyDiscrXYZ_DataSet distMagFunc, double dist, double mag) { + // buffer this check by 10km as we can have supersampled locations that are a little further + Preconditions.checkState((float)dist >= (float)distMagFunc.getMinX() + && (float)dist <= (float)(distMagFunc.getMaxX()+10), "Bad dist=%s with range=[%s, %s]", + (float)dist, (float)distMagFunc.getMinX(), (float)distMagFunc.getMaxX()); + if (dist < distMagFunc.getMinX()) + dist = distMagFunc.getMinX(); + else if (dist > distMagFunc.getMaxX()) + dist = distMagFunc.getMaxX(); + Preconditions.checkState((float)mag >= (float)distMagFunc.getMinY() + && (float)mag <= (float)distMagFunc.getMaxY(), "Bad mag=%s with range=[%s, %s]", + (float)mag, (float)distMagFunc.getMinY(), (float)distMagFunc.getMaxY()); + if (mag < distMagFunc.getMinY()) + mag = distMagFunc.getMinY(); + else if (mag > distMagFunc.getMaxY()) + mag = distMagFunc.getMaxY(); + return distMagFunc.bilinearInterpolation(dist, mag); + } + + @Override + public double getDistanceX(Location siteLoc) { + double dist = LocationUtils.horzDistanceFast(this.loc, siteLoc); + double rX = getDist(corrTables.rX, dist, mag); + boolean positive = rand.nextDouble() < getDist(corrTables.rXPositive, dist, mag); + return positive ? rX : -rX; + } + + @Override + public double getDistanceJB(Location siteLoc) { + return getDist(corrTables.rJB, LocationUtils.horzDistanceFast(loc, siteLoc), mag); + } + + @Override + public double getDistanceRup(Location siteLoc) { + return getDist(corrTables.rRup, LocationUtils.horzDistanceFast(loc, siteLoc), mag); + } + + @Override + public double getDistanceSeis(Location loc) { + throw new IllegalStateException(); + } + +} diff --git a/src/main/java/scratch/kevin/ucerf3/PureScratch.java b/src/main/java/scratch/kevin/ucerf3/PureScratch.java index 73bc71d3..52901bbe 100644 --- a/src/main/java/scratch/kevin/ucerf3/PureScratch.java +++ b/src/main/java/scratch/kevin/ucerf3/PureScratch.java @@ -41,6 +41,7 @@ import org.jfree.chart.ui.RectangleEdge; import org.jfree.data.Range; import org.opensha.commons.calc.FaultMomentCalc; +import org.opensha.commons.calc.GaussianDistCalc; import org.opensha.commons.calc.WeightedSampler; import org.opensha.commons.calc.magScalingRelations.magScalingRelImpl.Ellsworth_B_WG02_MagAreaRel; import org.opensha.commons.calc.magScalingRelations.magScalingRelImpl.Somerville_2006_MagAreaRel; @@ -298,6 +299,7 @@ import scratch.UCERF3.utils.UCERF2_A_FaultMapper; import scratch.UCERF3.utils.paleoRateConstraints.UCERF3_PaleoRateConstraintFetcher; import scratch.kevin.nshm23.dmCovarianceTests.RandomDefModSampleLevel; +import scratch.kevin.pointSources.InvCDF_RJBCorrPointSurface; import scratch.kevin.simCompare.SiteHazardCurveComarePageGen; import scratch.kevin.simulators.RSQSimCatalog; import scratch.kevin.simulators.RSQSimCatalog.Catalogs; @@ -6130,13 +6132,28 @@ private static void test299() throws IOException { "; meanRatio="+(float)(meanFromDist/mean)+"; covRatio="+(float)(covFromDist/cov)); } + private static void test300() throws IOException { + NormalDistribution normDist = new NormalDistribution(); + for (int i=0; i<5; i++) { + System.out.println("Sigma="+i); + System.out.println("\tCDF1: "+(float)GaussianDistCalc.getCDF(i)); + System.out.println("\tCDF2: "+(float)normDist.cumulativeProbability(i)); + } + for (int i=1; i<6; i++) { + System.out.println("Sigmas="+i); + DiscretizedFunc func = InvCDF_RJBCorrPointSurface.getSigmaSpacedProbs(i); + for (Point2D pt : func) + System.out.println("\t"+(float)pt.getX()+": "+(float)pt.getY()); + } + } + /** * @param args * @throws Exception */ public static void main(String[] args) throws Exception { try { - test299(); + test300(); } catch (Throwable t) { t.printStackTrace(); System.exit(1);