Source code for openquake.hmtk.seismicity.occurrence.weichert
# -*- 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.importwarningsimportnumpyasnpfromopenquake.hmtk.seismicity.occurrence.baseimport(SeismicityOccurrence,OCCURRENCE_METHODS,)fromopenquake.hmtk.seismicity.occurrence.utilsimport(input_checks,get_completeness_counts,)
[docs]@OCCURRENCE_METHODS.add("calculate",completeness=True,reference_magnitude=0.0,magnitude_interval=0.1,bvalue=1.0,itstab=1e-5,maxiter=1000,)classWeichert(SeismicityOccurrence):"""Class to Implement Weichert Algorithm"""
[docs]defcalculate(self,catalogue,config,completeness=None):"""Calculates b value and rate for mag ref"""bval,sigma_b,rate,sigma_rate,_aval,_sigma_a=self._calculate(catalogue,config,completeness)returnbval,sigma_b,rate,sigma_rate
[docs]defcalc(self,catalogue,config,completeness=None):"""Calculates GR params"""bval,sigma_b,_rate,_sigma_rate,aval,sigma_a=self._calculate(catalogue,config,completeness)returnbval,sigma_b,aval,sigma_a
def_calculate(self,catalogue,config,completeness=None):"""Calculates a, b values + rate for mag ref"""# Input checkscmag,ctime,ref_mag,_,config=input_checks(catalogue,config,completeness)if"dtime"notincatalogue.data.keys()ornotlen(catalogue.data["dtime"]):catalogue.data["dtime"]=catalogue.get_decimal_time()ifnotcatalogue.end_year:catalogue.update_end_year()ifcompletenessisNone:# start_year = float(np.min(catalogue.data["year"]))completeness=np.column_stack([ctime,cmag])# Apply Weichert preparationcent_mag,t_per,n_obs=get_completeness_counts(catalogue,completeness,config["magnitude_interval"])# A few more Weichert checkskey_list=config.keys()if("bvalue"notinkey_list)or(notconfig["bvalue"]):config["bvalue"]=1.0if("itstab"notinkey_list)or(notconfig["itstab"]):config["itstab"]=1e-5if("maxiter"notinkey_list)or(notconfig["maxiter"]):config["maxiter"]=1000bval,sigma_b,rate,sigma_rate,fn0,stdfn0=self.weichert_algorithm(t_per,cent_mag,n_obs,ref_mag,config["bvalue"],config["itstab"],config["maxiter"],)# if not config['reference_magnitude']:agr=np.log10(fn0)agr_sigma=np.log10(fn0+stdfn0)-np.log10(fn0)returnbval,sigma_b,rate,sigma_rate,agr,agr_sigma
[docs]defweichert_algorithm(self,tper,fmag,nobs,mrate=0.0,bval=1.0,itstab=1e-5,maxiter=1000):""" Weichert algorithm :param tper: length of observation period corresponding to magnitude :type tper: numpy.ndarray (float) :param fmag: central magnitude :type fmag: numpy.ndarray (float) :param nobs: number of events in magnitude increment :type nobs: numpy.ndarray (int) :keyword mrate: reference magnitude :type mrate: float :keyword bval: initial value for b-value :type beta: float :keyword itstab: stabilisation tolerance :type itstab: float :keyword maxiter: Maximum number of iterations :type maxiter: Int :returns: b-value, sigma_b, a-value, sigma_a :rtype: float """beta=bval*np.log(10.0)d_m=fmag[1]-fmag[0]itbreak=0snm=np.sum(nobs*fmag)nkount=np.sum(nobs)iteration=1whileitbreak!=1:beta_exp=np.exp(-beta*fmag)tjexp=tper*beta_exptmexp=tjexp*fmagsumexp=np.sum(beta_exp)stmex=np.sum(tmexp)sumtex=np.sum(tjexp)stm2x=np.sum(fmag*tmexp)dldb=stmex/sumtexifnp.isnan(stmex)ornp.isnan(sumtex):warnings.warn("NaN occurs in Weichert iteration")returnnp.nan,np.nan,np.nan,np.nan,np.nan,np.nan# raise ValueError('NaN occers in Weichert iteration')d2ldb2=nkount*((dldb**2.0)-(stm2x/sumtex))dldb=(dldb*nkount)-snmbetl=np.copy(beta)beta=beta-(dldb/d2ldb2)sigbeta=np.sqrt(-1.0/d2ldb2)ifnp.abs(beta-betl)<=itstab:# Iteration has reached convergencebval=beta/np.log(10.0)sigb=sigbeta/np.log(10.0)fngtm0=nkount*(sumexp/sumtex)fn0=fngtm0*np.exp((beta)*(fmag[0]-(d_m/2.0)))stdfn0=fn0/np.sqrt(nkount)a_m=fngtm0*np.exp((-beta)*(mrate-(fmag[0]-(d_m/2.0))))siga_m=a_m/np.sqrt(nkount)itbreak=1else:iteration+=1ifiteration>maxiter:warnings.warn("Maximum Number of Iterations reached")returnnp.nan,np.nan,np.nan,np.nan,np.nan,np.nanreturnbval,sigb,a_m,siga_m,fn0,stdfn0