Skip to content

Commit

Permalink
Merge pull request #364 from GEMScienceTools/cleaning
Browse files Browse the repository at this point in the history 
Cleaning smt plotting functions up
  • Loading branch information
@CB-quakemodel
CB-quakemodel authored Nov 1, 2023
2 parents 7ac6263 + 439fff7 commit 61723cb
Show file tree
Hide file tree
Showing 9 changed files with 248 additions and 485 deletions.
69 changes: 25 additions & 44 deletions openquake/smt/residuals/gmpe_residuals.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -930,7 +930,7 @@ def get_edr_values(self, bandwidth=0.01, multiplier=3.0):
edr_values[gmpe]["EDR"] = results[2]
return edr_values

def get_edr_values_wrt_spectral_period(self, bandwidth=0.01,
def get_edr_values_wrt_imt(self, bandwidth=0.01,
multiplier=3.0):
"""
Calculates the EDR values for each GMPE according to the Euclidean
Expand All @@ -947,12 +947,12 @@ def get_edr_values_wrt_spectral_period(self, bandwidth=0.01,
:param float multiplier:
"Multiplier of standard deviation (equation 8 of Kale and Akkar)
"""

self.edr_values_wrt_imt = OrderedDict([(gmpe, {}) for
gmpe in self.gmpe_list])
for gmpe in self.gmpe_list:
obs_wrt_imt, expected_wrt_imt, stddev_wrt_imt = self._get_edr_gmpe_information_wrt_spectral_period(gmpe)
results = self._get_edr_wrt_spectral_period(obs_wrt_imt,
(obs_wrt_imt, expected_wrt_imt, stddev_wrt_imt
) = self._get_edr_gmpe_information_wrt_imt(gmpe)
results = self._get_edr_wrt_imt(obs_wrt_imt,
expected_wrt_imt,
stddev_wrt_imt,
bandwidth,
Expand All @@ -965,7 +965,7 @@ def get_edr_values_wrt_spectral_period(self, bandwidth=0.01,
def _get_edr_gmpe_information(self, gmpe):
"""
Extract the observed ground motions, expected and total standard
deviation for the GMPE (aggregating over all IMS)
deviation for the GMPE (aggregating over all IMTs)
"""
obs = np.array([], dtype=float)
expected = np.array([], dtype=float)
Expand All @@ -979,35 +979,16 @@ def _get_edr_gmpe_information(self, gmpe):
context["Expected"][gmpe][imtx]["Total"]])
return obs, expected, stddev

def _get_edr_gmpe_information_wrt_spectral_period(self, gmpe):
def _get_edr_gmpe_information_wrt_imt(self, gmpe):
"""
Extract the observed ground motions, expected and total standard
deviation for the GMPE (per imt)
"""
#Remove non-acceleration imts from residuals.imts
imt_append=pd.DataFrame(self.imts,index=self.imts)
imt_append.columns=['imt']
for imt_idx in range(0,np.size(self.imts)):
if self.imts[imt_idx]=='PGV':
imt_append=imt_append.drop(imt_append.loc['PGV'])
if self.imts[imt_idx]=='PGD':
imt_append=imt_append.drop(imt_append.loc['PGD'])
if self.imts[imt_idx]=='CAV':
imt_append=imt_append.drop(imt_append.loc['CAV'])
if self.imts[imt_idx]=='Ia':
imt_append=imt_append.drop(imt_append.loc['Ia'])

imt_append_list=pd.DataFrame()
for idx in range(0,len(imt_append)):
imt_append_list[idx]=imt_append.iloc[idx]
imt_append=imt_append.reset_index()
imt_append_list.columns=imt_append.imt
self.imts_appended=list(imt_append_list)

# Get EDR values per imt
obs_wrt_imt={}
expected_wrt_imt={}
stddev_wrt_imt={}
for imtx in self.imts_appended:
for imtx in self.imts:
obs = np.array([], dtype=float)
expected = np.array([], dtype=float)
stddev = np.array([], dtype=float)
Expand All @@ -1020,7 +1001,7 @@ def _get_edr_gmpe_information_wrt_spectral_period(self, gmpe):
obs_wrt_imt[imtx]=obs
expected_wrt_imt[imtx]=expected
stddev_wrt_imt[imtx]=stddev

return obs_wrt_imt, expected_wrt_imt, stddev_wrt_imt

def _get_edr(self, obs, expected, stddev, bandwidth=0.01, multiplier=3.0):
Expand Down Expand Up @@ -1057,7 +1038,7 @@ def _get_edr(self, obs, expected, stddev, bandwidth=0.01, multiplier=3.0):
edr = np.sqrt(kappa * inv_n * np.sum(mde ** 2.))
return mde_norm, np.sqrt(kappa), edr

def _get_edr_wrt_spectral_period(self, obs_wrt_imt, expected_wrt_imt,
def _get_edr_wrt_imt(self, obs_wrt_imt, expected_wrt_imt,
stddev_wrt_imt, bandwidth=0.01,
multiplier=3.0):
"""
Expand All @@ -1068,31 +1049,31 @@ def _get_edr_wrt_spectral_period(self, obs_wrt_imt, expected_wrt_imt,
edr_wrt_imt={}
kappa_wrt_imt={}

for imt in self.imts_appended:
nvals = len(obs_wrt_imt[imt])
for imtx in self.imts:
nvals = len(obs_wrt_imt[imtx])
min_d = bandwidth / 2.
kappa_wrt_imt[imt] = self._get_edr_kappa(obs_wrt_imt[imt],
expected_wrt_imt[imt])
mu_d = obs_wrt_imt[imt] - expected_wrt_imt[imt]
d1c = np.fabs(obs_wrt_imt[imt] - (expected_wrt_imt[imt] - (
multiplier * stddev_wrt_imt[imt])))
d2c = np.fabs(obs_wrt_imt[imt] - (expected_wrt_imt[imt] + (
multiplier * stddev_wrt_imt[imt])))
kappa_wrt_imt[imtx] = self._get_edr_kappa(obs_wrt_imt[imtx],
expected_wrt_imt[imtx])
mu_d = obs_wrt_imt[imtx] - expected_wrt_imt[imtx]
d1c = np.fabs(obs_wrt_imt[imtx] - (expected_wrt_imt[imtx] - (
multiplier * stddev_wrt_imt[imtx])))
d2c = np.fabs(obs_wrt_imt[imtx] - (expected_wrt_imt[imtx] + (
multiplier * stddev_wrt_imt[imtx])))
dc_max = ceil(np.max(np.array([np.max(d1c), np.max(d2c)])))
num_d = len(np.arange(min_d, dc_max, bandwidth))
mde_wrt_imt = np.zeros(nvals)
for iloc in range(0, num_d):
d_val = (min_d + (float(iloc) * bandwidth)) * np.ones(nvals)
d_1 = d_val - min_d
d_2 = d_val + min_d
p_1 = norm.cdf((d_1 - mu_d) / stddev_wrt_imt[imt]) -\
norm.cdf((-d_1 - mu_d) / stddev_wrt_imt[imt])
p_2 = norm.cdf((d_2 - mu_d) / stddev_wrt_imt[imt]) -\
norm.cdf((-d_2 - mu_d) / stddev_wrt_imt[imt])
p_1 = norm.cdf((d_1 - mu_d) / stddev_wrt_imt[imtx]) -\
norm.cdf((-d_1 - mu_d) / stddev_wrt_imt[imtx])
p_2 = norm.cdf((d_2 - mu_d) / stddev_wrt_imt[imtx]) -\
norm.cdf((-d_2 - mu_d) / stddev_wrt_imt[imtx])
mde_wrt_imt += (p_2 - p_1) * d_val
inv_n = 1.0 / float(nvals)
mde_norm_wrt_imt[imt] = np.sqrt(inv_n * np.sum(mde_wrt_imt ** 2.))
edr_wrt_imt[imt] = np.sqrt(kappa_wrt_imt[imt] * inv_n * np.sum(
mde_norm_wrt_imt[imtx] = np.sqrt(inv_n * np.sum(mde_wrt_imt ** 2.))
edr_wrt_imt[imtx] = np.sqrt(kappa_wrt_imt[imtx] * inv_n * np.sum(
mde_wrt_imt ** 2.))

return mde_norm_wrt_imt, np.sqrt(pd.Series(kappa_wrt_imt)), edr_wrt_imt
Expand Down
Loading

0 comments on commit 61723cb

Please sign in to comment.