Kijko

openquake/hmtk/seismicity/occurrence/kijko_smit.py

@@ -101,15 +101,16 @@ def calculate(self, catalogue, config, completeness=None):
                 nyr[ival] = ctime[ival - 1] - ctime[ival]
             neq[ival] = np.sum(id1)
             # Get a- and b- value for the selected events
-            temp_rec_table = recurrence_table(catalogue.data['magnitude'][id1],
-                                              dmag,
-                                              catalogue.data['year'][id1])
+            if len(id1) > 0:
+                temp_rec_table = recurrence_table(
+                    catalogue.data['magnitude'][id1], dmag,
+                    catalogue.data['year'][id1])
 
-            aki_ml = AkiMaxLikelihood()
-            b_est[ival] = aki_ml._aki_ml(temp_rec_table[:, 0],
-                                         temp_rec_table[:, 1],
-                                         dmag, m_c)[0]
-            ival += 1
+                aki_ml = AkiMaxLikelihood()
+                b_est[ival] = aki_ml._aki_ml(temp_rec_table[:, 0],
+                                             temp_rec_table[:, 1],
+                                             dmag, m_c)[0]
+            #ival += 1
         total_neq = np.float(np.sum(neq))
         bval = self._harmonic_mean(b_est, neq)
         sigma_b = bval / np.sqrt(total_neq)
@@ -141,7 +142,24 @@ def _calculate_a_value(self, bval, nvalue, nyr, cmag, ref_mag):
         """
         Calculates the rate of events >= ref_mag using the b-value estimator
         and Eq. 10 of Kijko & Smit
-        """
 
-        denominator = np.sum(nyr * np.exp(-bval * (cmag - ref_mag)))
-        return nvalue / denominator
+        :param bval:
+            b-value
+        :param nvalue:
+            Number of earthquakes within the completeness window
+        :param nyr:
+            A vector with the duration [yr] of each completeness interval
+        :param cmag:
+            A vector with the magnitude lower limit of each completeness
+            interval
+        :param
+        """
+        # Computing the rate for eqs above the min magnitude included in the
+        # completeness window
+        mmin = min(cmag)
+        denominator = np.sum(nyr * np.exp(-bval * (cmag - mmin)))
+        rate_above_mmin = nvalue / denominator
+        # Computing aGR
+        agr = np.log10(rate_above_mmin) + bval * mmin
+        # Returning rate above the reference magnitude provided
+        return 10**(agr-bval*ref_mag)

openquake/hmtk/seismicity/occurrence/utils.py

@@ -126,14 +126,14 @@ def input_checks(catalogue, config, completeness):
         config = {'reference_magnitude': None,
                   'magnitude_interval': 0.1}
     else:
-        if (not 'reference_magnitude' in config.keys()) or\
+        if ('reference_magnitude' not in config.keys()) or\
                 (config['reference_magnitude'] is None):
             ref_mag = 0.
             config['reference_magnitude'] = None
         else:
             ref_mag = config['reference_magnitude']
 
-        if (not 'magnitude_interval' in config.keys()) or \
+        if ('magnitude_interval' not in config.keys()) or \
                 not config['magnitude_interval']:
             dmag = 0.1
         else:

openquake/hmtk/tests/seismicity/occurrence/kijko_smit_test.py

@@ -105,13 +105,15 @@ def test_kijko_smit_maximum_likelihood(self):
         """
         bval, sigma_b, aval, sigma_a = self.ks_ml.calculate(
             self.catalogue, self.config, self.compl)
-        print(bval, sigma_b)
         self.assertAlmostEqual(self.bval, bval, 1)
+        self.assertAlmostEqual(5.81, aval, 2)
 
     def test_kijko_smit_set_reference_magnitude(self):
         completeness_table = np.array([[1900, 1.0]])
         catalogue = Catalogue.make_from_dict(
             {'magnitude': np.array([5.0, 6.0]),
              'year': np.array([2000, 2000])})
         config = {'reference_magnitude': 0.0}
-        self.ks_ml.calculate(catalogue, config, completeness_table)
+        bval, sigma_b, aval, sigma_a = self.ks_ml.calculate(
+            catalogue, config, completeness_table)
+        self.assertAlmostEqual(-0.9136, aval, 2)