# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## Copyright (C) 2013-2025 GEM Foundation## OpenQuake is free software: you can redistribute it and/or modify it# under the terms of the GNU Affero General Public License as published# by the Free Software Foundation, either version 3 of the License, or# (at your option) any later version.## OpenQuake is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU Affero General Public License for more details.## You should have received a copy of the GNU Affero General Public License# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>."""Module exports :class:`NGAEastGMPE` and :class:`NGAEastGMPETotalSigma`"""importioimportosimportnumpyasnpfromcopyimportdeepcopyfromscipy.statsimportchi2fromopenquake.hazardlib.gsim.baseimportCoeffsTable,add_aliasfromopenquake.hazardlib.gsim.gmpe_tableimportGMPETable,_get_meanfromopenquake.hazardlibimportconstfromopenquake.hazardlib.imtimportPGA,SA# Common interpolation function
def_scaling(mean_tau,sd_tau2):""" Returns the chi-2 scaling factor from the mean and variance of the uncertainty model, as reported in equation 5.4 of Al Atik (2015) """returnsd_tau2**2./(2.0*mean_tau**2.)def_dof(mean_tau,sd_tau2):""" Returns the degrees of freedom for the chi-2 distribution from the mean and variance of the uncertainty model, as reported in equation 5.5 of Al Atik (2015) """return2.0*mean_tau**4./sd_tau2**2.def_at_percentile(tau,var_tau,percentile):""" Returns the value of the inverse chi-2 distribution at the given percentile from the mean and variance of the uncertainty model, as reported in equations 5.1 - 5.3 of Al Atik (2015) """assertpercentile>=0.0andpercentile<=1.0c_val=_scaling(tau,var_tau)k_val=_dof(tau,var_tau)returnnp.sqrt(c_val*chi2.ppf(percentile,df=k_val))# Mean tau values from the global model - Table 5.1GLOBAL_TAU_MEAN={"PGV":{"tau1":0.3733,"tau2":0.3639,"tau3":0.3434,"tau4":0.3236},"SA":{"tau1":0.4518,"tau2":0.4270,"tau3":0.3863,"tau4":0.3508}}# Standard deviation of tau values from the global model - Table 5.1GLOBAL_TAU_SD={"PGV":{"tau1":0.0558,"tau2":0.0554,"tau3":0.0477,"tau4":0.0449},"SA":{"tau1":0.0671,"tau2":0.0688,"tau3":0.0661,"tau4":0.0491}}# Mean tau values from the CENA constant-tau modelCENA_CONSTANT_TAU_MEAN={"PGV":{"tau":0.3441},"SA":{"tau":0.3695}}# Standard deviation of tau values from CENA constant-tau modelCENA_CONSTANT_TAU_SD={"PGV":{"tau":0.0554},"SA":{"tau":0.0688}}# Mean tau values from the CENA tau modelCENA_TAU_MEAN={"PGV":{"tau1":0.3477,"tau2":0.3281,"tau3":0.3092},"SA":{"tau1":0.3730,"tau2":0.3375,"tau3":0.3064}}# Standard deviation of tau values from the CENA tau modelCENA_TAU_SD={"PGV":{"tau1":0.0554,"tau2":0.0477,"tau3":0.0449},"SA":{"tau1":0.0688,"tau2":0.0661,"tau3":0.0491}}
[docs]defglobal_tau(imt,mag,params):""" 'Global' model of inter-event variability, as presented in equation 5.6 (p103) """ifimt.string=="PGV":C=params["PGV"]else:C=params["SA"]tau=np.full_like(mag,C["tau1"])tau[mag>6.5]=C["tau4"]idx=(mag>5.5)&(mag<=6.5)tau[idx]=ITPL(mag[idx],C["tau4"],C["tau3"],5.5,1.0)idx=(mag>5.0)&(mag<=5.5)tau[idx]=ITPL(mag[idx],C["tau3"],C["tau2"],5.0,0.5)idx=(mag>4.5)&(mag<=5.0)tau[idx]=ITPL(mag[idx],C["tau2"],C["tau1"],4.5,0.5)returntau
[docs]defcena_constant_tau(imt,mag,params):""" Returns the inter-event tau for the constant tau case """ifimt.string=="PGV":returnparams["PGV"]["tau"]else:returnparams["SA"]["tau"]
[docs]defcena_tau(imt,mag,params):""" Returns the inter-event standard deviation, tau, for the CENA case """ifimt.string=="PGV":C=params["PGV"]else:C=params["SA"]tau=np.full_like(mag,C["tau1"])tau[mag>6.5]=C["tau3"]idx=(mag>5.5)&(mag<=6.5)tau[idx]=ITPL(mag[idx],C["tau3"],C["tau2"],5.5,1.0)idx=(mag>5.0)&(mag<=5.5)tau[idx]=ITPL(mag[idx],C["tau2"],C["tau1"],5.0,0.5)returntau
[docs]defget_tau_at_quantile(mean,stddev,quantile):""" Returns the value of tau at a given quantile in the form of a dictionary organised by intensity measure """tau_model={}forimtinmean:tau_model[imt]={}forkeyinmean[imt]:ifquantileisNone:tau_model[imt][key]=mean[imt][key]else:tau_model[imt][key]=_at_percentile(mean[imt][key],stddev[imt][key],quantile)returntau_model
[docs]defget_phi_ss_at_quantile(phi_model,quantile):""" Returns the phi_ss values at the specified quantile as an instance of `class`:: openquake.hazardlib.gsim.base.CoeffsTable - applies to the magnitude-dependent cases """# Setup SA coeffs - the backward compatible Python 2.7 waycoeffs=deepcopy(phi_model.sa_coeffs)coeffs.update(phi_model.non_sa_coeffs)forimtincoeffs:ifquantileisNone:coeffs[imt]={"a":phi_model[imt]["mean_a"],"b":phi_model[imt]["mean_b"]}else:coeffs[imt]={"a":_at_percentile(phi_model[imt]["mean_a"],phi_model[imt]["var_a"],quantile),"b":_at_percentile(phi_model[imt]["mean_b"],phi_model[imt]["var_b"],quantile)}returnCoeffsTable.fromdict(coeffs)
[docs]defget_phi_s2ss_at_quantile(phi_model,quantile):""" Returns the phi_s2ss value for all periods at the specific quantile as an instance of `class`::openquake.hazardlib.gsim.base.CoeffsTable """# Setup SA coeffs - the backward compatible Python 2.7 waycoeffs=deepcopy(phi_model.sa_coeffs)coeffs.update(phi_model.non_sa_coeffs)forimtincoeffs:ifquantileisNone:coeffs[imt]={"phi_s2ss":phi_model[imt]["mean"]}else:coeffs[imt]={"phi_s2ss":_at_percentile(phi_model[imt]["mean"],phi_model[imt]["var"],quantile)}returnCoeffsTable.fromdict(coeffs)
# Gather the models to setup the phi_ss values for the given quantilePHI_SETUP={"cena":PHI_SS_CENA,"cena_constant":PHI_SS_CENA_CONSTANT,"global":PHI_SS_GLOBAL}# Phi site-to-site model for th Central & Eastern USPHI_S2SS_MODEL={"cena":PHI_S2SS_CENA}
[docs]defget_phi_ss(imt,mag,params):""" Returns the single station phi (or it's variance) for a given magnitude and intensity measure type according to equation 5.14 of Al Atik (2015) """C=params[imt]phi=C["a"]+(mag-5.0)*(C["b"]-C["a"])/1.5phi[mag<=5.0]=C["a"]phi[mag>6.5]=C["b"]returnphi
# ############################ helper functions ########################def_get_f760(C_F760,vs30,CONSTANTS,is_stddev=False):""" Returns very hard rock to hard rock (Vs30 760 m/s) adjustment factor taken as the Vs30-dependent weighted mean of two reference condition factors: for impedence and for gradient conditions. The weighting model is described by equations 5 - 7 of Stewart et al. (2019) """wimp=(CONSTANTS["wt1"]-CONSTANTS["wt2"])*\
(np.log(vs30/CONSTANTS["vw2"])/np.log(CONSTANTS["vw1"]/CONSTANTS["vw2"]))+CONSTANTS["wt2"]wimp[vs30>=CONSTANTS["vw1"]]=CONSTANTS["wt1"]wimp[vs30<CONSTANTS["vw2"]]=CONSTANTS["wt2"]wgr=1.0-wimpifis_stddev:returnwimp*C_F760["f760is"]+wgr*C_F760["f760gs"]else:returnwimp*C_F760["f760i"]+wgr*C_F760["f760g"]def_get_fv(C_LIN,vs30s,f760,CONSTANTS):""" Returns the Vs30-dependent component of the mean linear amplification model, as defined in equation 3 of Stewart et al. (2019) """const1=C_LIN["c"]*np.log(C_LIN["v1"]/CONSTANTS["vref"])const2=C_LIN["c"]*np.log(C_LIN["v2"]/CONSTANTS["vref"])f_v=C_LIN["c"]*np.log(vs30s/CONSTANTS["vref"])f_v[vs30s<=C_LIN["v1"]]=const1f_v[vs30s>C_LIN["v2"]]=const2idx=vs30s>CONSTANTS["vU"]ifnp.any(idx):const3=np.log(3000./CONSTANTS["vU"])f_v[idx]=const2-(const2+f760[idx])*\
(np.log(vs30s[idx]/CONSTANTS["vU"])/const3)idx=vs30s>=3000.ifnp.any(idx):f_v[idx]=-f760[idx]returnf_v+f760
[docs]defget_fnl(C_NL,pga_rock,vs30,period,us23=None):""" Returns the nonlinear mean amplification according to equation 2 of Hashash et al. (2019) """ifperiod<=0.4:vref=760.else:vref=3000.f_nl=np.zeros(vs30.shape)f_rk=np.log((pga_rock+C_NL["f3"])/C_NL["f3"])idx=vs30<C_NL["Vc"]ifnp.any(idx):# f2 term of the mean nonlinear amplification model# according to equation 3 of Hashash et al., (2019)c_vs=np.copy(vs30[idx])c_vs[c_vs>vref]=vrefifus23:# US 2023 NSHMPf_4=C_NL["f4"]*0.5+C_NL["f4_mod"]*0.5else:f_4=C_NL["f4"]f_2=f_4*(np.exp(C_NL["f5"]*(c_vs-360.))-np.exp(C_NL["f5"]*(vref-360.)))f_nl[idx]=f_2*f_rk[idx]returnf_nl,f_rk
[docs]defget_linear_stddev(C_LIN,vs30,CONSTANTS):""" Returns the standard deviation of the linear amplification function, as defined in equation 4 of Stewart et al., (2019) """sigma_v=C_LIN["sigma_vc"]+np.zeros(vs30.shape)idx=vs30<C_LIN["vf"]ifnp.any(idx):dsig=C_LIN["sigma_L"]-C_LIN["sigma_vc"]d_v=(vs30[idx]-CONSTANTS["vL"])/\
(C_LIN["vf"]-CONSTANTS["vL"])sigma_v[idx]=C_LIN["sigma_L"]-(2.*dsig*d_v)+\
dsig*(d_v**2.)idx=np.logical_and(vs30>C_LIN["v2"],vs30<=CONSTANTS["vU"])ifnp.any(idx):d_v=(vs30[idx]-C_LIN["v2"])/(CONSTANTS["vU"]-C_LIN["v2"])sigma_v[idx]=C_LIN["sigma_vc"]+ \
(C_LIN["sigma_U"]-C_LIN["sigma_vc"])*(d_v**2.)idx=vs30>=CONSTANTS["vU"]ifnp.any(idx):sigma_v[idx]=C_LIN["sigma_U"]*\
(1.-(np.log(vs30[idx]/CONSTANTS["vU"])/np.log(3000./CONSTANTS["vU"])))sigma_v[vs30>3000.]=0.0returnsigma_v
[docs]defget_nonlinear_stddev(C_NL,vs30):""" Returns the standard deviation of the nonlinear amplification function, as defined in equation 2.5 of Hashash et al. (2017) """sigma_f2=np.zeros(vs30.shape)sigma_f2[vs30<300.]=C_NL["sigma_c"]idx=np.logical_and(vs30>=300,vs30<1000)ifnp.any(idx):sigma_f2[idx]=(-C_NL["sigma_c"]/np.log(1000./300.))*\
np.log(vs30[idx]/300.)+C_NL["sigma_c"]returnsigma_f2
[docs]defget_hard_rock_mean(self,mag,ctx,imt):""" Returns the mean and standard deviations for the reference very hard rock condition (Vs30 = 3000 m/s) """# return Distance Tablesimls=self.mean_table['%.2f'%mag,imt.string]# Get distance vector for the given magnitudeidx=np.searchsorted(self.m_w,mag)dists=self.distances[:,0,idx-1]dst=getattr(ctx,self.distance_type)# get log(mean)returnnp.log(_get_mean(self.kind,imls,dst,dists))
[docs]defget_site_amplification(self,imt,pga_r,vs30s,us23=False):""" Returns the sum of the linear (Stewart et al., 2019) and non-linear (Hashash et al., 2019) amplification terms """# Get the coefficients for the IMTC_LIN=self.COEFFS_LINEAR[imt]C_F760=self.COEFFS_F760[imt]C_NL=self.COEFFS_NONLINEAR[imt]ifstr(imt).startswith("PGA"):period=0.01elifstr(imt).startswith("PGV"):period=0.5else:period=imt.period# Get f760f760=_get_f760(C_F760,vs30s,self.CONSTANTS)# Get the linear amplification factorf_lin=_get_fv(C_LIN,vs30s,f760,self.CONSTANTS)# Get the nonlinear amplification from Hashash et al., (2017)f_nl,f_rk=get_fnl(C_NL,pga_r,vs30s,period,us23)# Mean amplificationampl=f_lin+f_nl# If an epistemic uncertainty is required then retrieve the epistemic# sigma of both models and multiply by the input epsilonifself.site_epsilon:site_epistemic=get_site_amplification_sigma(self,vs30s,f_rk,C_LIN,C_F760,C_NL)ampl+=self.site_epsilon*site_epistemicreturnampl
[docs]defget_site_amplification_sigma(self,vs30s,f_rk,C_LIN,C_F760,C_NL):""" Returns the epistemic uncertainty on the site amplification factor """# In the case of the linear model sigma_f760 and sigma_fv are# assumed independent and the resulting sigma_flin is the root# sum of squares (SRSS)f760_stddev=_get_f760(C_F760,vs30s,self.CONSTANTS,is_stddev=True)f_lin_stddev=np.sqrt(f760_stddev**2.+get_linear_stddev(C_LIN,vs30s,self.CONSTANTS)**2)# Likewise, the epistemic uncertainty on the linear and nonlinear# model are assumed independent and the SRSS is takenf_nl_stddev=get_nonlinear_stddev(C_NL,vs30s)*f_rkreturnnp.sqrt(f_lin_stddev**2.+f_nl_stddev**2.)
[docs]defget_stddevs(self,mag,imt):""" Returns the standard deviations for either the ergodic or non-ergodic models """ifself.__class__.__name__.endswith('TotalSigma'):return[_get_total_sigma(self,imt,mag),0.,0.]tau=_get_tau(self,imt,mag)phi=_get_phi(self,imt,mag)sigma=np.sqrt(tau**2+phi**2)return[sigma,tau,phi]
def_get_tau(self,imt,mag):""" Returns the inter-event standard deviation (tau) """returnTAU_EXECUTION[self.tau_model](imt,mag,self.TAU)def_get_phi(self,imt,mag):""" Returns the within-event standard deviation (phi) """phi=get_phi_ss(imt,mag,self.PHI_SS)ifself.ergodic:C=self.PHI_S2SS[imt]phi=np.sqrt(phi**2.+C["phi_s2ss"]**2.)returnphi
[docs]defget_mean_amp(self,mag,ctx,imt,u_adj=None,cstl=None):""" Compute mean ground-motion. NOTE: If a non-zero cpa_term is computed we apply this adjustment post-application of the collapsed epistemic uncertainty site ampl model to replicate the USGS java code. Therefore the cpa_term adj is added to the collapsed site epi. uncertainty adjusted mean in the compute meth of the NGAEastUSGSGMPE class (usgs_ceus_2019.py). :param u_adj: Array containing the period-dependent bias adjustment as required for the 2023 Conterminous US NSHMP. :param cstl: Dictionary containing parameters required for the computation of the Coastal Plains site amplification model of Chapman and Guo (2021) as required for the 2023 Conterminous US NSHMP. """# Get the PGA on the reference rock conditionifPGAinself.DEFINED_FOR_INTENSITY_MEASURE_TYPES:rock_imt=PGA()else:rock_imt=SA(0.01)pga_r=get_hard_rock_mean(self,mag,ctx,rock_imt)# Apply US 2023 period-dep bias adj if requiredifisinstance(u_adj,np.ndarray):pga_r+=u_adj# Get the desired spectral acceleration on rockifimt.string!="PGA":# Calculate the ground motion at required spectral period for# the reference rockmean=get_hard_rock_mean(self,mag,ctx,imt)# Again apply US 2023 period-dep bias adj if requiredifisinstance(u_adj,np.ndarray):mean+=u_adjelse:# Avoid re-calculating PGA if that was already done!mean=np.copy(pga_r)# If applying US 2023 NSHMP bias adjustment OR Coastal Plains# site amp adjustment use the alternative f4 coeff (f4_mod)# for non-linear site term as within USGS' NGAEast java codeifisinstance(cstl,dict)orisinstance(u_adj,np.ndarray):us23=Trueelse:us23=False# Get site amplificationamp=get_site_amplification(self,imt,np.exp(pga_r),ctx.vs30,us23)# Add the site term to the meanmean+=amp# Get the Chapman and Guo (2021) Coastal Plains adjustment if req.ifisinstance(cstl,dict)andisinstance(cstl['f_cpa'],np.ndarray):# Get site amp term with vs30 ref of 1000 m/svs_cg21=np.full(len(ctx.vs30),1000.)amp_cpa=get_site_amplification(self,imt,np.exp(pga_r),vs_cg21)# Then compute adjustment per sitecpa_term=cstl['f_cpa']-amp_cpa*cstl['z_scale']else:cpa_term=np.full(len(ctx.vs30),0.)# Turn off adjustmentreturnmean,amp,pga_r,cpa_term
[docs]classNGAEastGMPE(GMPETable):""" A generalised base class for the implementation of a GMPE in which the mean values are determined from tables (input by the user) and the standard deviation model taken from Al Atik (2015). Should be common to all NGA East ground motion models. The site amplification model is developed using the model described by Stewart et al. (2019) and Hashash et al. (2019). The model contains a linear and a non-linear component of amplification. The linear model is described in Stewart et al., (2017) and Stewart et al (2019), with the latter taken as the authoritative source where differences arise: Stewart, J. P., Parker, G. A., Harmon, J. A., Atkinson, G. A., Boore, D. M., Darragh, R. B., Silva, W. J. and Hashash, Y. M. A. (2017) "Expert Panel Recommendations for Ergodic Site Amplification in Central and Eastern North America", PEER Report No. 2017/04, Pacific Earthquake Engineering Research Center, University of California, Berkeley. Stewart, J. P., Parker, G. A., Atkinson, G. M., Boore, D. M., Hashash, Y. M. A. and Silva, W. J. (2019) "Ergodic Site Amplification Model for Central and Eastern North America", Earthquake Spectra, in press The nonlinear model is described in Hashash et al. (2017) and Hashash et al. (2019), with the latter taken as the authoritarive source where differences arise: Hashash, Y. M. A., Harmon, J. A., Ilhan, O., Parker, G. and Stewart, J. P. (2017), "Recommendation for Ergonic Nonlinear Site Amplification in Central and Eastern North America", PEER Report No. 2017/05, Pacific Earthquake Engineering Research Center, University of California, Berkeley. Hashash, Y. M. A., Ilhan, O., Harmon, J. A., Parker, G. A., Stewart, J. P. Rathje, E. M., Campbell, K. W., and Silva, W. J. (2019) "Nonlinear site amplification model for ergodic seismic hazard analysis in Central and Eastern North America", Earthquake Spectra, in press Note that the uncertainty provided in this model is treated as an epistemic rather than aleatory variable. As such there is no modification of the standard deviation model used for the bedrock case. The epistemic uncertainty can be added to the model by the user input site_epsilon term, which describes the number of standard deviations by which to multiply the epistemic uncertainty model, to then be added or subtracted from the median amplification model :param str tau_model: Choice of model for the inter-event standard deviation (tau), selecting from "global" {default}, "cena" or "cena_constant" :param str phi_model: Choice of model for the single-station intra-event standard deviation (phi_ss), selecting from "global" {default}, "cena" or "cena_constant" :param str phi_s2ss_model: Choice of station-term s2ss model. Can be either "cena" or None. When None is input then the non-ergodic model is used :param TAU: Inter-event standard deviation model :param PHI_SS: Single-station standard deviation model :param PHI_S2SS: Station term for ergodic standard deviation model :param bool ergodic: True if an ergodic model is selected, False otherwise :param float site_epsilon: Number of standard deviations above or below median for the uncertainty in the site amplification model """PATH=os.path.join(os.path.dirname(__file__),"nga_east_tables")DEFINED_FOR_INTENSITY_MEASURE_COMPONENT=const.IMC.RotD50DEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL,const.StdDev.INTER_EVENT,const.StdDev.INTRA_EVENT}# Requires Vs30 only - common to all modelsREQUIRES_SITES_PARAMETERS={'vs30'}kind="nga_east"def__init__(self,gmpe_table,tau_model="global",phi_model="global",phi_s2ss_model=None,tau_quantile=None,phi_ss_quantile=None,phi_s2ss_quantile=None,site_epsilon=None):""" Instantiates the class with additional terms controlling which type of aleatory uncertainty model is preferred ('global', 'cena_constant', 'cena'), and the quantile of the epistemic uncertainty model (float in the range 0 to 1). :param float tau_quantile: Epistemic uncertainty quantile for the inter-event standard deviation models. Float in the range 0 to 1, or None (mean value used) :param float phi_ss_quantile: Epistemic uncertainty quantile for the intra-event standard deviation models. Float in the range 0 to 1, or None (mean value used) :param float phi_s2ss_quantile: Epistemic uncertainty quantile for the site-to-site standard deviation models. Float in the range 0 to 1, or None (mean value used) """self.tau_model=tau_modelself.phi_model=phi_modelself.phi_s2ss_model=phi_s2ss_modelself.TAU=Noneself.PHI_SS=Noneself.PHI_S2SS=Noneifself.phi_s2ss_model:self.ergodic=Trueelse:self.ergodic=Falseself.tau_quantile=tau_quantileself.phi_ss_quantile=phi_ss_quantileself.phi_s2ss_quantile=phi_s2ss_quantileifself.kind!='usgs':# setup tauself.TAU=get_tau_at_quantile(TAU_SETUP[self.tau_model]["MEAN"],TAU_SETUP[self.tau_model]["STD"],self.tau_quantile)# setup phiself.PHI_SS=get_phi_ss_at_quantile(PHI_SETUP[self.phi_model],self.phi_ss_quantile)# if required setup phis2ssifself.ergodic:self.PHI_S2SS=get_phi_s2ss_at_quantile(PHI_S2SS_MODEL[self.phi_s2ss_model],self.phi_s2ss_quantile)self.site_epsilon=site_epsilonifnotisinstance(gmpe_table,io.BytesIO):# real path namegmpe_table=os.path.join(self.PATH,os.path.basename(gmpe_table))assertos.path.exists(gmpe_table),gmpe_tablesuper().__init__(gmpe_table)
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" Returns the mean and standard deviations """[mag]=np.unique(np.round(ctx.mag,2))# by constructionform,imtinenumerate(imts):mean[m],_,_,_=get_mean_amp(self,mag,ctx,imt)# Get standard deviation modelsig[m],tau[m],phi[m]=get_stddevs(self,ctx.mag,imt)
# For the total standard deviation model the magnitude boundaries depend on# the model selectedMAG_LIMS_KEYS={"cena":{"mag":[5.0,5.5,6.5],"keys":["tau1","tau2","tau3"]},"cena_constant":{"mag":[np.inf],"keys":["tau"]},"global":{"mag":[4.5,5.0,5.5,6.5],"keys":["tau1","tau2","tau3","tau4"]}}def_get_sigma_at_quantile(self,sigma_quantile):""" Calculates the total standard deviation at the specified quantile """# Mean mean is found in self.TAU. Get the variance in tautau_std=TAU_SETUP[self.tau_model]["STD"]# Mean phiss is found in self.PHI_SS. Get the variance in phiphi_std=deepcopy(self.PHI_SS.sa_coeffs)phi_std.update(self.PHI_SS.non_sa_coeffs)forkeyinphi_std:phi_std[key]={"a":PHI_SETUP[self.phi_model][key]["var_a"],"b":PHI_SETUP[self.phi_model][key]["var_b"]}ifself.ergodic:# IMT list should be taken from the PHI_S2SS_MODELimt_list=list(PHI_S2SS_MODEL[self.phi_s2ss_model].non_sa_coeffs)imt_list+=list(PHI_S2SS_MODEL[self.phi_s2ss_model].sa_coeffs)else:imt_list=list(phi_std)phi_std=CoeffsTable.fromdict(phi_std)tau_bar,tau_std=_get_tau_vector(self,self.TAU,tau_std,imt_list)phi_bar,phi_std=_get_phi_vector(self,self.PHI_SS,phi_std,imt_list)sigma={}# Calculate the total standard deviationforimtinimt_list:sigma[imt]={}fori,keyinenumerate(self.tau_keys):# Calculates the expected standard deviationsigma_bar=np.sqrt(tau_bar[imt][i]**2.+phi_bar[imt][i]**2.)# Calculated the variance in the standard deviationsigma_std=np.sqrt(tau_std[imt][i]**2.+phi_std[imt][i]**2.)# The keys swap from tau to sigmanew_key=key.replace("tau","sigma")ifsigma_quantileisnotNone:sigma[imt][new_key]=_at_percentile(sigma_bar,sigma_std,sigma_quantile)else:sigma[imt][new_key]=sigma_barself.tau_keys[i]=new_keyself.SIGMA=CoeffsTable.fromdict(sigma)def_get_tau_vector(self,tau_mean,tau_std,imt_list):""" Gets the vector of mean and variance of tau values corresponding to the specific model and returns them as dictionaries """self.magnitude_limits=np.array(MAG_LIMS_KEYS[self.tau_model]["mag"])self.tau_keys=MAG_LIMS_KEYS[self.tau_model]["keys"]t_bar={}t_std={}tau_model=TAU_EXECUTION[self.tau_model]forimtinimt_list:t_bar[imt]=tau_model(imt,self.magnitude_limits,tau_mean)t_std[imt]=tau_model(imt,self.magnitude_limits,tau_std)returnt_bar,t_stddef_get_phi_vector(self,phi_mean,phi_std,imt_list):""" Gets the vector of mean and variance of phi values corresponding to the specific model and returns them as dictionaries """p_bar={}p_std={}forimtinimt_list:p_bar[imt]=ss_mean=get_phi_ss(imt,self.magnitude_limits,phi_mean)p_std[imt]=ss_std=get_phi_ss(imt,self.magnitude_limits,phi_std)ifself.ergodic:# Add on the phi_s2ss term according to Eqs. 5.15 and 5.16# of Al Atik (2015)ss_mean[:]=np.sqrt(ss_mean**2.+PHI_S2SS_MODEL[self.phi_s2ss_model][imt]["mean"]**2.)ss_std[:]=np.sqrt(ss_std**2.+PHI_S2SS_MODEL[self.phi_s2ss_model][imt]["var"]**2.)returnp_bar,p_stddef_get_total_sigma(self,imt,mag):""" Returns the estimated total standard deviation for a given intensity measure type and magnitude """[mag]=np.unique(np.round(mag,2))# by constructionC=self.SIGMA[imt]ifmag<=self.magnitude_limits[0]:# The CENA constant model is always returned herereturnC[self.tau_keys[0]]elifmag>self.magnitude_limits[-1]:returnC[self.tau_keys[-1]]else:# Needs interpolationforiinrange(len(self.tau_keys)-1):l_m=self.magnitude_limits[i]u_m=self.magnitude_limits[i+1]ifmag>l_mandmag<=u_m:returnITPL(mag,C[self.tau_keys[i+1]],C[self.tau_keys[i]],l_m,u_m-l_m)
[docs]classNGAEastGMPETotalSigma(NGAEastGMPE):""" The Al Atik (2015) standard deviation model defines mean and quantiles for the inter- and intra-event residuals. However, it also defines separately a total standard deviation model with expectation and quantiles. As the inter- and intra-event quantile cannot be recovered unambiguously from the total standard deviation quantile this form of the model is defined only for total standard deviation. Most likely it is this form that would be used for seismic hazard analysis. :param SIGMA: Total standard deviation model at quantile :param list magnitude_limits: Magnitude limits corresponding to the selected standard deviation model :param list tau_keys: Keys for the tau values corresponding to the selected standard deviation model """DEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL}def__init__(self,gmpe_table,tau_model="global",phi_model="global",phi_s2ss_model=None,tau_quantile=None,phi_ss_quantile=None,phi_s2ss_quantile=None,site_epsilon=None,sigma_quantile=None):""" Instantiates the model call the BaseNGAEastModel to return the expected TAU, PHI_SS and PHI_S2SS values then uses these to calculate the expected total standard deviation and its variance according to equations 5.15, 5.16 and/or 5.18 and 5.19 of Al Atik (2015) :param float sigma_quantile: Quantile of the epistemic uncertainty model for the total standard deviation. Should be float between 0 and 1, or None (mean value taken) """super().__init__(gmpe_table,tau_model,phi_model,phi_s2ss_model,tau_quantile,phi_ss_quantile,phi_s2ss_quantile,site_epsilon)# Upon instantiation the TAU, PHI_SS, and PHI_S2SS objects contain# the mean valuesself.SIGMA=Noneself.magnitude_limits=[]self.tau_keys=[]_get_sigma_at_quantile(self,sigma_quantile)
# populate gsim_aliases for the NGA East GMPEslines='''\Boore2015NGAEastA04 BOORE_A04_J15Boore2015NGAEastAB14 BOORE_AB14_J15Boore2015NGAEastAB95 BOORE_AB95_J15Boore2015NGAEastBCA10D BOORE_BCA10D_J15Boore2015NGAEastBS11 BOORE_BS11_J15Boore2015NGAEastSGD02 BOORE_SGD02_J15DarraghEtAl2015NGAEast1CCSP DARRAGH_1CCSPDarraghEtAl2015NGAEast1CVSP DARRAGH_1CVSPDarraghEtAl2015NGAEast2CCSP DARRAGH_2CCSPDarraghEtAl2015NGAEast2CVSP DARRAGH_2CVSPYenierAtkinson2015NGAEast YENIER_ATKINSONPezeschkEtAl2015NGAEastM1SS PEZESCHK_M1SSPezeschkEtAl2015NGAEastM2ES PEZESCHK_M2ESFrankel2015NGAEast FRANKEL_J15ShahjoueiPezeschk2015NGAEast SHAHJOUEI_PEZESCHKAlNomanCramer2015NGAEast ALNOMAN_CRAMERGraizer2015NGAEast GRAIZERHassaniAtkinson2015NGAEast HASSANI_ATKINSONHollenbackEtAl2015NGAEastGP PEER_GPHollenbackEtAl2015NGAEastEX PEER_EX'''.splitlines()forlineinlines:alias,key=line.split()add_alias(alias,NGAEastGMPE,gmpe_table=f"NGAEast_{key}.hdf5")add_alias(alias+'TotalSigma',NGAEastGMPETotalSigma,gmpe_table=f"NGAEast_{key}.hdf5")
