From 1c2ca51c8e4b38fbd392a24713774f4b90d16522 Mon Sep 17 00:00:00 2001
From: Marco Pagani <marco.pagani@globalquakemodel.org>
Date: Wed, 18 Aug 2021 12:45:24 +0200
Subject: [PATCH 1/2] Fixing bug in calculation of a

---
 .../hmtk/seismicity/occurrence/kijko_smit.py  | 23 ++++++++++++++++---
 .../seismicity/occurrence/kijko_smit_test.py  |  6 +++--
 2 files changed, 24 insertions(+), 5 deletions(-)

diff --git a/openquake/hmtk/seismicity/occurrence/kijko_smit.py b/openquake/hmtk/seismicity/occurrence/kijko_smit.py
index 666cda15ca74..7a8771c4a069 100644
--- a/openquake/hmtk/seismicity/occurrence/kijko_smit.py
+++ b/openquake/hmtk/seismicity/occurrence/kijko_smit.py
@@ -141,7 +141,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)
diff --git a/openquake/hmtk/tests/seismicity/occurrence/kijko_smit_test.py b/openquake/hmtk/tests/seismicity/occurrence/kijko_smit_test.py
index 60a9f8054472..a9281dacfc97 100644
--- a/openquake/hmtk/tests/seismicity/occurrence/kijko_smit_test.py
+++ b/openquake/hmtk/tests/seismicity/occurrence/kijko_smit_test.py
@@ -105,8 +105,8 @@ 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]])
@@ -114,4 +114,6 @@ def test_kijko_smit_set_reference_magnitude(self):
             {'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)

From 3f4c15dc0d2e96664f5311777a5f7ae450ec2329 Mon Sep 17 00:00:00 2001
From: Marco Pagani <marco.pagani@globalquakemodel.org>
Date: Fri, 20 Aug 2021 19:01:27 +0200
Subject: [PATCH 2/2] minor changes

---
 .../hmtk/seismicity/occurrence/kijko_smit.py    | 17 +++++++++--------
 openquake/hmtk/seismicity/occurrence/utils.py   |  4 ++--
 2 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/openquake/hmtk/seismicity/occurrence/kijko_smit.py b/openquake/hmtk/seismicity/occurrence/kijko_smit.py
index 7a8771c4a069..252071446db6 100644
--- a/openquake/hmtk/seismicity/occurrence/kijko_smit.py
+++ b/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)
diff --git a/openquake/hmtk/seismicity/occurrence/utils.py b/openquake/hmtk/seismicity/occurrence/utils.py
index 6a6174af910f..2e6c76997354 100644
--- a/openquake/hmtk/seismicity/occurrence/utils.py
+++ b/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: