Source code for openquake.hmtk.seismicity.occurrence.b_maximum_likelihood
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## LICENSE## Copyright (C) 2010-2025 GEM Foundation, G. Weatherill, M. Pagani,# D. Monelli.## The Hazard Modeller's Toolkit 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.## You should have received a copy of the GNU Affero General Public License# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>## DISCLAIMER## The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein# is released as a prototype implementation on behalf of# scientists and engineers working within the GEM Foundation (Global# Earthquake Model).## It is distributed for the purpose of open collaboration and in the# hope that it will be useful to the scientific, engineering, disaster# risk and software design communities.## The software is NOT distributed as part of GEM’s OpenQuake suite# (https://www.globalquakemodel.org/tools-products) and must be considered as a# separate entity. The software provided herein is designed and implemented# by scientific staff. It is not developed to the design standards, nor# subject to same level of critical review by professional software# developers, as GEM’s OpenQuake software suite.## Feedback and contribution to the software is welcome, and can be# directed to the hazard scientific staff of the GEM Model Facility# (hazard@globalquakemodel.org).## The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed WITHOUT# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License# for more details.## The GEM Foundation, and the authors of the software, assume no# liability for use of the software.""""""importnumpyasnpfromopenquake.hmtk.seismicity.occurrence.utilsimport(input_checks,recurrence_table,)fromopenquake.hmtk.seismicity.occurrence.baseimport(SeismicityOccurrence,OCCURRENCE_METHODS,)fromopenquake.hmtk.seismicity.occurrence.aki_maximum_likelihoodimport(AkiMaxLikelihood,)
[docs]@OCCURRENCE_METHODS.add("calculate",**{"completeness":True,"reference_magnitude":0.0,"magnitude_interval":0.1,"Average Type":["Weighted","Harmonic"],},)classBMaxLikelihood(SeismicityOccurrence):"""Implements maximum likelihood calculations taking into account time variation in completeness" """
[docs]defcalculate(self,catalogue,config,completeness=None):"""Calculates recurrence parameters a_value and b_value, and their respective uncertainties :param catalogue: Earthquake Catalogue An instance of :class:`openquake.hmtk.seismicity.catalogue` :param dict config: A configuration dictionary; the only parameter that can be defined in this case if the type of average to be applied in the calculation :param list or numpy.ndarray completeness: Completeness table """# Input checkscmag,ctime,ref_mag,dmag,config=input_checks(catalogue,config,completeness)# Check the configurationifconfig["Average Type"]notin["Weighted","Harmonic"]:raiseValueError("Average type not recognised in bMaxLiklihood!")returnself._b_ml(catalogue,config,cmag,ctime,ref_mag,dmag)
def_b_ml(self,catalogue,config,cmag,ctime,ref_mag,dmag):end_year=float(catalogue.end_year)catalogue=catalogue.dataival=0mag_eq_tolerance=1e-5aki_ml=AkiMaxLikelihood()whileival<np.shape(ctime)[0]:id0=np.abs(ctime-ctime[ival])<mag_eq_tolerancem_c=np.min(cmag[id0])print("--- ctime",ctime[ival]," m_c",m_c)# Find events later than cut-off year, and with magnitude# greater than or equal to the corresponding completeness# magnitude. m_c - mag_eq_tolerance is required to correct# floating point differences.id1=np.logical_and(catalogue["year"]>=ctime[ival],catalogue["magnitude"]>=(m_c-mag_eq_tolerance),)# Get a- and b- value for the selected eventstemp_rec_table=recurrence_table(catalogue["magnitude"][id1],dmag,catalogue["year"][id1],end_year-ctime[ival]+1,)bval,sigma_b=aki_ml._aki_ml(temp_rec_table[:,0],temp_rec_table[:,1],dmag,m_c)ifival==0:gr_pars=np.array([np.hstack([bval,sigma_b])])neq=np.sum(id1)# Number of eventselse:gr_pars=np.vstack([gr_pars,np.hstack([bval,sigma_b])])neq=np.hstack([neq,np.sum(id1)])ival=ival+np.sum(id0)# Get average GR parametersbval,sigma_b=self._average_parameters(gr_pars,neq,config["Average Type"])aval=self._calculate_a_value(bval,float(np.sum(neq)),cmag,ctime,catalogue["magnitude"],end_year,dmag,)sigma_a=self._calculate_a_value(bval+sigma_b,float(np.sum(neq)),cmag,ctime,catalogue["magnitude"],end_year,dmag,)ifnotconfig["reference_magnitude"]:returnbval,sigma_b,aval,sigma_a-avalelse:rate=10.0**(aval-bval*config["reference_magnitude"])sigma_rate=(10.0**(sigma_a-bval*config["reference_magnitude"])-rate)returnbval,sigma_b,rate,sigma_ratedef_average_parameters(self,gr_params,neq,average_type="Weighted"):""" Calculates the average of a set of Gutenberg-Richter parameters depending on the average type :param numpy.ndarray gr_params: Gutenberg-Richter parameters [b, sigma_b, a, sigma_a] :param numpy.ndarray neq: """ifnp.shape(gr_params)[0]!=neq.size:raiseValueError("Number of weights does not correspond"" to number of parameters")if"Harmonic"inaverage_type:average_parameters=self._harmonic_mean(gr_params,neq)else:average_parameters=self._weighted_mean(gr_params,neq)bval=average_parameters[0]sigma_b=average_parameters[1]returnbval,sigma_bdef_calculate_a_value(self,bvalue,nvalue,cmag,cyear,magnitude,end_year,dmag):""" Calculates the a-value using the method of Weichert (1980) and McGuire (2004) """mmin=cmag[0]mmax=np.max(magnitude)ifmmax>np.max(cmag):cmag=np.hstack([cmag,mmax+dmag])target_mag=(cmag[:-1]+cmag[1:])/2.0nyear=end_year-cyear+1.0beta=bvalue*np.log(10.0)rate_mmin=(nvalue*np.sum(np.exp(-beta*target_mag))/np.sum(nyear*np.exp(-beta*target_mag)))returnnp.log10(rate_mmin)+bvalue*mmindef_weighted_mean(self,parameters,neq):"""Simple weighted mean"""weight=neq.astype(float)/np.sum(neq)ifnp.shape(parameters)[0]!=weight.size:raiseValueError("Parameter vector not same shape as weights")else:average_value=np.zeros(np.shape(parameters)[1],dtype=float)forilocinrange(0,np.shape(parameters)[1]):average_value[iloc]=np.sum(parameters[:,iloc]*weight)returnaverage_valuedef_harmonic_mean(self,parameters,neq):"""Harmonic mean"""weight=neq.astype(float)/np.sum(neq)ifnp.shape(parameters)[0]!=weight.size:raiseValueError("Parameter vector not same shape as weights")average_value=np.zeros(np.shape(parameters)[1],dtype=float)forilocinrange(0,np.shape(parameters)[1]):average_value[iloc]=1.0/np.sum((weight*(1.0/parameters[:,iloc])))returnaverage_value