Skip to content

Commit

Permalink
Merge pull request #373 from GEMScienceTools/tidy_and_fix_dict_issue
Browse files Browse the repository at this point in the history 
tidy + fix a dict storage issue
  • Loading branch information
@CB-quakemodel
CB-quakemodel authored Nov 17, 2023
2 parents 241a269 + 1e661d8 commit d8286e9
Show file tree
Hide file tree
Showing 8 changed files with 333 additions and 164 deletions.
6 changes: 3 additions & 3 deletions docsrc/contents/smt.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ We can specify the inputs to perform a residual analysis with as follows:
with_betw_ratio = 1.7 # add between-event and within-event sigma using ratio of 1.7 to partition total sigma
[models.3-CampbellBozorgnia2014]
set_between_epsilon = 1.5 # set between-event epsilon (i.e. tau epsilon)
set_between_epsilon = 0.5 # Shift the mean with formula mean --> mean + epsilon_tau * between event
[models.1-AbrahamsonEtAl2014]
median_scaling_scalar = 1.4 # scale median by factor of 1.4 over all imts
Expand Down Expand Up @@ -153,7 +153,7 @@ We can specify the inputs to perform a residual analysis with as follows:
[models.KothaEtAl2020ESHM20] # ESHM20 model
sigma_mu_epsilon = 2.85697
c3_epsilon = 1.72
region = 4 # Note that within residuals specify region here, whereas in comparison module toml (below) we specify using the eshm20_region param
region = 4 # Note that within the residuals toml we specify the region here, whereas in the comparison module toml (below) we specify the region for all ESHM20 GMMs uniformly using the eshm20_region param
[imts]
imt_list = ['PGA', 'SA(0.2)', 'SA(0.5)', 'SA(1.0']
Expand Down Expand Up @@ -466,7 +466,6 @@ Comparing GMPEs
maxR = 300 # max dist. used in trellis, Sammon's, clusters and matrix plots
dist_type = 'rrup' # or rjb, repi or rhypo (dist type used in trellis plots)
dist_list = [10, 100, 250] # distance intervals for use in spectra plots
region = 0 # for NGAWest2 GMPE regionalisation
eshm20_region = 2 # for KothaEtAl2020 ESHM20 GMPE regionalisation
Nstd = 1 # num. of std. dev. to sample sigma for in median prediction (0, 1, 2 or 3)
Expand All @@ -476,6 +475,7 @@ Comparing GMPEs
Z1 = -999
Z25 = -999
up_or_down_dip = 1 # 1 = up-dip, 0 = down-dip
region = 'Global' # get region specific z1pt0 and zpt50 ('Global' or 'Japan')
# Characterise earthquake for the region of interest as finite rupture
[source_properties]
Expand Down
281 changes: 200 additions & 81 deletions openquake/smt/comparison/compare_gmpes.py
View file

Large diffs are not rendered by default.

20 changes: 9 additions & 11 deletions openquake/smt/comparison/utils_compare_gmpes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,13 +212,11 @@ def plot_spectra_util(trt, ztor, rake, strike, dip, depth, Z1, Z25, Vs30,
tau_store.append(tau_dist)

if Nstd != 0:
rs_plus_sigma_dist = np.exp(f(i)+(Nstd*sigma_dist))
rs_minus_sigma_dist = np.exp(f(i)-(Nstd*sigma_dist))

if Nstd != 0:
rs_plus_sigma_dist = np.exp(f(i)+(Nstd*sigma_dist))
rs_minus_sigma_dist = np.exp(f(i)-(Nstd*sigma_dist))
rs_plus_sigma.append(rs_plus_sigma_dist)
rs_minus_sigma.append(rs_minus_sigma_dist)

# Plot individual GMPEs
if 'plot_lt_only' not in str(gmpe):
ax1.plot(period, rs_50p, color = col, linewidth=2,
Expand Down Expand Up @@ -806,9 +804,9 @@ def spectra_data(gmpe, Nstd, rs_50p, rs_plus_sigma, rs_minus_sigma,

# And if Nstd != 0 store these weighted branches too
if Nstd != 0:
lt_vals_plus_sig_gmc1[gmpe] = {
lt_vals_plus_sig_gmc1[gmpe,'p_sig'] = {
'plus_sigma': rs_plus_sigma_weighted_gmc1}
lt_vals_minus_sig_gmc1[gmpe] = {
lt_vals_minus_sig_gmc1[gmpe,'m_sig'] = {
'minus_sigma': rs_minus_sigma_weighted_gmc1}

# Logic tree #2
Expand All @@ -831,9 +829,9 @@ def spectra_data(gmpe, Nstd, rs_50p, rs_plus_sigma, rs_minus_sigma,

# And if Nstd != 0 store these weighted branches too
if Nstd != 0:
lt_vals_plus_sig_gmc2[gmpe] = {
lt_vals_plus_sig_gmc2[gmpe,'p_sig'] = {
'plus_sigma': rs_plus_sigma_weighted_gmc2}
lt_vals_minus_sig_gmc2[gmpe] = {
lt_vals_minus_sig_gmc2[gmpe,'m_sig'] = {
'minus_sigma': rs_minus_sigma_weighted_gmc2}

return (lt_vals_gmc1, lt_vals_plus_sig_gmc1, lt_vals_minus_sig_gmc1,
Expand Down Expand Up @@ -872,10 +870,10 @@ def lt_spectra(ax1, gmpe, gmpe_list, Nstd, period, gmc1_or_gmc2,
lt_df_gmc[gmpe].loc['median']))
if Nstd != 0:
wt_plus_sigma_per_gmpe_gmc[gmpe] = np.array(
pd.Series(lt_df_plus_sigma_gmc[gmpe].loc[
pd.Series(lt_df_plus_sigma_gmc[gmpe,'p_sig'].loc[
'plus_sigma']))
wt_minus_sigma_per_gmpe_gmc[gmpe] = np.array(
pd.Series(lt_df_minus_sigma_gmc[gmpe].loc[
pd.Series(lt_df_minus_sigma_gmc[gmpe,'m_sig'].loc[
'minus_sigma']))

lt_df_gmc = pd.DataFrame(wt_per_gmpe_gmc, index=period)
Expand Down
27 changes: 9 additions & 18 deletions openquake/smt/comparison/utils_gmpes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,32 +216,23 @@ def att_curves(gmpe, orig_gmpe, depth, mag, aratio, strike, dip, rake, Vs30,

def _get_z1(Vs30, region):
"""
:param region:
Choose among: region= 0 for global; 1 for California; 2 for Japan;
3 for China; 4 for Italy; 5 for Turkey (locally = 3); 6 for Taiwan (
locally = 0)
Get z1pt0 using Chiou and Youngs (2014) relationship.
"""
if region == 2: # in California and non-Japan region
Z1 = (np.exp(-5.23 / 2 * np.log((Vs30**2 + 412.39**2) /
(1360**2 + 412.39**2))))
else:
Z1 = (np.exp(-7.15 / 4 * np.log((Vs30**4 + 570.94**4) /
(1360**4 + 570.94**4))))

if region == 'Global': # California and non-Japan regions
Z1 = np.exp(-7.15/4*np.log((Vs30**4 + 570.94**4)/(1360**4 + 570.94**4)))
else: # Japan
Z1 = np.exp(-5.23/2*np.log((Vs30**2 + 412.39**2)/(1360**2 + 412.39**2)))
return Z1


def _get_z25(Vs30, region):
"""
:param region:
Choose among: region= 0 for global; 1 for California; 2 for Japan;
3 for China; 4 for Italy; 5 for Turkey (locally = 3); 6 for Taiwan (
locally = 0)
Get z2pt5 using Campbell and Bozorgnia (2014) relationship.
"""
if region == 2: # in Japan
Z25 = np.exp(5.359 - 1.102 * np.log(Vs30))
else:
if region == 'Global': # California and non-Japan regions
Z25 = np.exp(7.089 - 1.144 * np.log(Vs30))
else: # Japan
Z25 = np.exp(5.359 - 1.102 * np.log(Vs30))

return Z25

Expand Down
157 changes: 109 additions & 48 deletions openquake/smt/tests/comparison/comparison_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@

# Defines the target values for each run in the inputted .toml file
TARGET_VS30 = 800
TARGET_REGION = 0
TARGET_REGION = 'Global'
TARGET_TRELLIS_DEPTHS = [20,25,30]
TARGET_RMAX = 300
TARGET_NSTD = 2
Expand Down Expand Up @@ -118,15 +118,24 @@ def test_mtxs_median_calculation(self):
# Check each parameter matches target
config = comp.Configurations(self.input_file)

mtxs_medians = compute_matrix_gmpes(config.trt, config.ztor,
config.imt_list, config.mag_list,
config.gmpes_list, config.rake,
config.strike, config.dip,
mtxs_medians = compute_matrix_gmpes(config.trt,
config.ztor,
config.imt_list,
config.mag_list,
config.gmpes_list,
config.rake,
config.strike,
config.dip,
config.depth_for_non_trel_or_rs_fun,
config.Z1, config.Z25, config.Vs30,
config.region, config.maxR,
config.aratio, config.eshm20_region,
config.dist_type, mtxs_type='median',
config.Z1,
config.Z25,
config.Vs30,
config.region,
config.maxR,
config.aratio,
config.eshm20_region,
config.dist_type,
mtxs_type='median',
up_or_down_dip=config.up_or_down_dip)

# Check correct number of imts
Expand All @@ -147,15 +156,24 @@ def test_sammons_and_euclidean_distance_matrix_functions(self):
config = comp.Configurations(self.input_file)

# Get medians
mtxs_medians = compute_matrix_gmpes(config.trt, config.ztor,
mtxs_medians = compute_matrix_gmpes(config.trt,
config.ztor,
config.imt_list,
config.mag_list, config.gmpes_list,
config.rake, config.strike, config.dip,
config.mag_list,
config.gmpes_list,
config.rake,
config.strike,
config.dip,
config.depth_for_non_trel_or_rs_fun,
config.Z1, config.Z25, config.Vs30,
config.region, config.maxR,
config.aratio, config.eshm20_region,
config.dist_type, mtxs_type='median',
config.Z1,
config.Z25,
config.Vs30,
config.region,
config.maxR,
config.aratio,
config.eshm20_region,
config.dist_type,
mtxs_type='median',
up_or_down_dip=config.up_or_down_dip)

# Sammons checks
Expand Down Expand Up @@ -196,14 +214,21 @@ def test_clustering_median(self):
config = comp.Configurations(self.input_file)

# Median clustering checks
mtxs_medians = compute_matrix_gmpes(config.trt, config.ztor,
config.imt_list, config.mag_list,
config.gmpes_list, config.rake,
config.strike, config.dip,
mtxs_medians = compute_matrix_gmpes(config.trt,
config.ztor,
config.imt_list,
config.mag_list,
config.gmpes_list,
config.rake,
config.strike,
config.dip,
config.depth_for_non_trel_or_rs_fun,
config.Z1, config.Z25,
config.Vs30, config.region,
config.maxR, config.aratio,
config.Z1,
config.Z25,
config.Vs30,
config.region,
config.maxR,
config.aratio,
config.eshm20_region,
config.dist_type,
mtxs_type='median',
Expand Down Expand Up @@ -232,14 +257,21 @@ def test_clustering_84th_perc(self):
config = comp.Configurations(self.input_file)

# Median clustering checks
mtxs_medians = compute_matrix_gmpes(config.trt, config.ztor,
config.imt_list, config.mag_list,
config.gmpes_list, config.rake,
config.strike, config.dip,
mtxs_medians = compute_matrix_gmpes(config.trt,
config.ztor,
config.imt_list,
config.mag_list,
config.gmpes_list,
config.rake,
config.strike,
config.dip,
config.depth_for_non_trel_or_rs_fun,
config.Z1, config.Z25,
config.Vs30, config.region,
config.maxR, config.aratio,
config.Z1,
config.Z25,
config.Vs30,
config.region,
config.maxR,
config.aratio,
config.eshm20_region,
config.dist_type,
mtxs_type='84th_perc',
Expand Down Expand Up @@ -268,27 +300,56 @@ def test_trellis_and_spectra_functions(self):
config = comp.Configurations(self.input_file)

# Trellis plots
plot_trellis_util(config.trt, config.ztor, config.rake, config.strike,
config.dip, config.trellis_and_rs_depth, config.Z1,
config.Z25, config.Vs30, config.region, config.imt_list,
config.trellis_and_rs_mag_list, config.maxR, config.gmpes_list,
config.aratio, config.Nstd, self.output_directory,
config.custom_color_flag, config.custom_color_list,
config.eshm20_region, config.dist_type,
config.lt_weights_gmc1, config.lt_weights_gmc2,
plot_trellis_util(config.trt,
config.ztor,
config.rake, config.strike,
config.dip,
config.trellis_and_rs_depth,
config.Z1,
config.Z25,
config.Vs30,
config.region,
config.imt_list,
config.trellis_and_rs_mag_list,
config.maxR,
config.gmpes_list,
config.aratio,
config.Nstd,
self.output_directory,
config.custom_color_flag,
config.custom_color_list,
config.eshm20_region,
config.dist_type,
config.lt_weights_gmc1,
config.lt_weights_gmc2,
up_or_down_dip=config.up_or_down_dip)

# Spectra plots
plot_spectra_util(config.trt, config.ztor, config.rake, config.strike,
config.dip, config.trellis_and_rs_depth, config.Z1,
config.Z25, config.Vs30, config.region,
config.max_period, config.trellis_and_rs_mag_list,
config.dist_list, config.gmpes_list, config.aratio,
config.Nstd, self.output_directory,
config.custom_color_flag, config.custom_color_list,
config.eshm20_region, config.dist_type,
config.lt_weights_gmc1, config.lt_weights_gmc2,
obs_spectra=None, up_or_down_dip=config.up_or_down_dip)
plot_spectra_util(config.trt,
config.ztor,
config.rake,
config.strike,
config.dip,
config.trellis_and_rs_depth,
config.Z1,
config.Z25,
config.Vs30,
config.region,
config.max_period,
config.trellis_and_rs_mag_list,
config.dist_list,
config.gmpes_list,
config.aratio,
config.Nstd,
self.output_directory,
config.custom_color_flag,
config.custom_color_list,
config.eshm20_region,
config.dist_type,
config.lt_weights_gmc1,
config.lt_weights_gmc2,
obs_spectra=None,
up_or_down_dip=config.up_or_down_dip)

# Specify target files
target_file_trellis = (os.path.join(
Expand Down
2 changes: 1 addition & 1 deletion openquake/smt/tests/comparison/data/compare_gmpe_inputs.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ max_period = 2 # max period for response spectra (can't exceed maximum period in
maxR = 300 # max distance to consider
dist_type = 'rrup' # specify distance metric used in trellis plots
dist_list = [10, 50, 200] # distances for use in response spectra + sigma wrt period plots
region = 0 # region for NGAWest2 GMMs
eshm20_region = 0 # region for ESHM20 GMMs
Nstd = 2 # number of standard deviations to sample sigma for

Expand All @@ -15,6 +14,7 @@ vs30 = 800
Z1 = -999
Z25 = -999
up_or_down_dip = 1 # 1 = up-dip, 0 = down-dip
region = 'Global' # region specific basin effects

[source_properties] # Characterise EQ as finite rupture
trt = 'ASCR'
Expand Down
2 changes: 1 addition & 1 deletion openquake/smt/tests/comparison/data/mgmpe_inputs.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ max_period = 2 # max period for response spectra (can't exceed maximum period in
maxR = 300 # max distance to consider
dist_type = 'rrup' # specify distance metric used in trellis plots
dist_list = [10, 50, 200] # distances for use in response spectra + sigma wrt period plots
region = 0 # region for NGAWest2 GMMs
eshm20_region = 0 # region for ESHM20 GMMs
Nstd = 1 # number of standard deviations to sample sigma for

Expand All @@ -15,6 +14,7 @@ vs30 = 800
Z1 = -999
Z25 = -999
up_or_down_dip = 1 # 1 = up-dip, 0 = down-dip
region = 'Global' # region specific basin effects

[source_properties] # Characterise EQ as finite rupture
trt = 'ASCR'
Expand Down
2 changes: 1 addition & 1 deletion openquake/smt/tests/file_samples/example_comparison_module_inputs.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ max_period = 2 # max period for spectra plots
maxR = 300 # max dist. used in trellis, Sammon's, clusters and matrix plots
dist_type = 'rrup' # or rjb, repi or rhypo (dist type used in trellis plots)
dist_list = [10, 100, 250] # distance intervals for use in spectra plots
region = 0 # for NGAWest2 GMPE regionalisation
eshm20_region = 2 # for KothaEtAl2020 ESHM20 GMPE regionalisation
Nstd = 1 # num. of std. dev. to sample sigma for in median prediction (0, 1, 2 or 3)

Expand All @@ -16,6 +15,7 @@ vs30 = 800
Z1 = -999
Z25 = -999
up_or_down_dip = 1 # 1 = up-dip, 0 = down-dip
region = 'Global' # region specific basin effects

# Characterise earthquake for the region of interest as finite rupture
[source_properties]
Expand Down

0 comments on commit d8286e9

Please sign in to comment.