wip: model test updates

docsrc/contents/cat.rst

@@ -194,7 +194,7 @@ This would give us a different fit to our data and a different equation to suppl
 
 Where there are not enough events to allow for a direct regression or we are unhappy with the fit for our data, there are many conversions in the literature which may be useful. This process may take some revising and iterating - it is sometimes very difficult to identify a best fit, especially where we have few datapoints or highly uncertain data. Once we are happy with the fits to our data, we can add the regression equation to the homogenisation .toml file. This process should be repeated for every magnitude we wish to convert to Mw. 
 
-The final homogenisation step itself is also controlled by a toml file, where each observed magnitude is specified individually and the regression coefficients and uncertainty are included. It is also necessary to specify a hierarchy of catalogues so that a preferred catalogue is used for the magnitude where the event has multiple entries. In the example below, we merge the ISCGEM and a local catalogue, preferring ISCGEM magnitudes where available as specified in the ranking. Because the ISCGEM already provides magnitudes in Mw, we simply retain all Mw magnitudes from ISCGEM. In this example, our local catalogue has two different magnitude types for which we have derived a regression. We specify how to convert to the standardised Mw from the local.mw and the standard deviations, which are outputs of the fitting we carried out above. 
+The final homogenisation step itself is also controlled by a toml file, where each observed magnitude is specified individually and the regression coefficients and uncertainty are included. It is also necessary to specify a hierarchy of catalogues so that a preferred catalogue is used for the magnitude where the event has multiple entries. If you have an isf format catalogue, you can also specify here the hierarchy of individual agencies (or authors in the isf format) within the catalogue. In the example below, we merge the ISCGEM and a local catalogue, preferring ISCGEM magnitudes where available as specified in the ranking. Because the ISCGEM already provides magnitudes in Mw, we simply retain all Mw magnitudes from ISCGEM. In this example, our local catalogue has two different magnitude types for which we have derived a regression. We specify how to convert to the standardised Mw from the local.mw and the standard deviations, which are outputs of the fitting we carried out above.
 
 .. code-block:: ini
 
@@ -229,6 +229,7 @@ The final homogenisation step itself is also controlled by a toml file, where ea
 	conv_eqs = ["0.1928 + 0.9757 * m"]
 	std_devs = [0.0091, 0.0016]
 
+The order of conversions in the list will determine priority for conversion, so for the local events we will first convert all events with mw magnitudes and then use mww only where the mw magnitudes are not available, and the local conversions will not be used when we have an ISCGEM Mw. In this way we can specify hierarchies for both the agencies and the magnitudes. 
 The actual homogenisation step is carried out by calling
 :code:`oqm cat homogenise $ARG1 $ARG2 $ARG3`
 as in the bash script example, where $ARG1 is the homogenisation toml file and and $ARG2 and $ARG3 are the hdf5 file outputs from the merge step, describing the origins and magnitude information for the merged catalogue respectively.

openquake/cat/completeness/analysis.py

@@ -176,7 +176,7 @@ def check_criterion(criterion, rate, previous_norm, tvars):
         norm = abs(tmp_rate - rate_ma)
 
     elif criterion == 'optimize':
-        tmp_rate = -1
+        tmp_rate = -1 
         norm = get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, 
                                  last_year, info=False)
 
@@ -190,20 +190,19 @@ def check_criterion(criterion, rate, previous_norm, tvars):
                                    mmin=ref_mag, mmax=ref_upp_mag)
     elif criterion == 'optimize_c':
         tmp_rate = -1
-        norm = get_norm_optimize_c(tcat, aval, bval, ctab, last_year)
+        norm = get_norm_optimize_c(tcat, aval, bval, ctab, last_year, ref_mag, ref_upp_mag, binw)
 
     elif criterion == 'gft':
         tmp_rate = -1
         norm = get_norm_optimize_gft(tcat, aval, bval, ctab, cmag, n_obs,
                                      t_per, last_year)
-
     elif criterion == 'weichert':
         tmp_rate = -1
         norm = get_norm_optimize_weichert(tcat, aval, bval, ctab, last_year)
 
     elif criterion == 'poisson':
-        tmp_rate = -1
-        norm = get_norm_optimize_c(tcat, aval, bval, ctab, last_year, ref_mag)
+        tmp_rate = -1        
+        norm = get_norm_optimize_c(tcat, aval, bval, ctab, last_year, ref_mag, ref_upp_mag, binw)
 
     if norm is None or np.isnan(norm):
         return False, -1, previous_norm
@@ -233,7 +232,7 @@ def _make_ctab(prm, years, mags):
         return 'skip'
 
 
-def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag,
+def _completeness_analysis(fname, years, mags, binw, ref_mag, mag_low, ref_upp_mag,
                            bgrlim, criterion, compl_tables, src_id=None,
                            folder_out_figs=None, rewrite=False,
                            folder_out=None):
@@ -247,6 +246,9 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag,
     :param ref_mag:
         The reference magnitude used to compute the rate and select the
         completeness table
+    :param mag_low:
+        The lowest magnitude to include in rate calculations.
+        Will default to ref_mag if not provided.
     :param ref_upp_mag:
         The reference upper magnitude limit used to compute the rate and
         select the completeness table
@@ -264,7 +266,7 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag,
 
     # Checking input
     if criterion not in ['match_rate', 'largest_rate', 'optimize', 'weichert',
-                         'poisson', 'optimize_a', 'optimize_b', 'optimize_d']:
+                         'poisson', 'optimize_a', 'optimize_b', 'optimize_c','optimize_d', 'gft']:
         raise ValueError('Unknown optimization criterion')
 
     tcat = _load_catalogue(fname)
@@ -273,7 +275,7 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag,
 
     # Info
     # Should have option to specify a mag_low != ref_mag
-    mag_low = ref_mag
+    #mag_low = lower_mag
     idx = tcat.data["magnitude"] >= mag_low
     fmt = 'Catalogue contains {:d} events equal or above {:.1f}'
     print('\nSOURCE:', src_id)
@@ -310,7 +312,7 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag,
         print(f'Iteration: {iper:05d} norm: {norm:12.6e}', end="\r")
 
         ctab = _make_ctab(prm, years, mags)
-        #print(ctab)
+
         if isinstance(ctab, str):
             continue
 
@@ -327,8 +329,9 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag,
         assert np.all(np.diff(ctab[:, 1]) >= 0)
 
         # Compute occurrence
-
+        #print('binwidth for completeness counts: ', binw)
         cent_mag, t_per, n_obs = get_completeness_counts(tcat, ctab, binw)
+        #print("cmag = ", cent_mag)
         if len(cent_mag) == 0:
             continue
         wei_conf['reference_magnitude'] = min(ctab[:, 1])
@@ -466,15 +469,19 @@ def read_compl_params(config):
     key = 'completeness'
     ms = np.array(config[key]['mags'], dtype=float)
     yrs = np.array(config[key]['years'])
-    bw = config.get('bin_width', 0.1)
+    try: 
+        bw = np.array(config[key]['bin_width'], dtype =float)
+    except: 
+    	bw = config.get('bin_width', 0.1)
     r_m = config[key].get('ref_mag', 5.0)
+    m_low = config[key].get('mag_low', r_m)
     r_up_m = config[key].get('ref_upp_mag', None)
     bmin = config[key].get('bmin', 0.8)
     bmax = config[key].get('bmax', 1.2)
     # Options: 'largest_rate', 'match_rate', 'optimize'
     crit = config[key].get('optimization_criterion', 'optimize')
 
-    return ms, yrs, bw, r_m, r_up_m, bmin, bmax, crit
+    return ms, yrs, bw, r_m, m_low, r_up_m, bmin, bmax, crit
 
 
 def read_compl_data(folder_in):
@@ -513,7 +520,7 @@ def completeness_analysis(fname_input_pattern, f_config, folder_out_figs,
     # Loading configuration
     config = toml.load(f_config)
 
-    ms, yrs, bw, r_m, r_up_m, bmin, bmax, crit = read_compl_params(config)
+    ms, yrs, bw, r_m, m_low, r_up_m, bmin, bmax, crit = read_compl_params(config)
 
     compl_tables, mags_chk, years_chk = read_compl_data(folder_in)
 
@@ -544,13 +551,13 @@ def completeness_analysis(fname_input_pattern, f_config, folder_out_figs,
         else:
             var = {}
 
-        res = _completeness_analysis(fname, yrs, ms, bw, r_m,
+        res = _completeness_analysis(fname, yrs, ms, bw, r_m, m_low,
                                      r_up_m, [bmin, bmax], crit,
                                      compl_tables, src_id,
                                      folder_out_figs=folder_out_figs,
                                      folder_out=folder_out,
                                      rewrite=False)
-        #print(len(res))
+
         if len(res) == 0:
             continue
 

openquake/cat/completeness/generate.py

@@ -108,7 +108,8 @@ def _get_completenesses(mags, years, folder_out=None, num_steps=0,
         min_mag = min(mags)
     else:
         min_mag = min_mag_compl
-
+
+    print("minimum magnitude is ", min_mag)
     if len(np.where(min_mag <= mags)) < 1:
         msg = 'None of the magnitude intervals above the min_mag_compl'
         raise ValueError(msg)

openquake/cat/completeness/norms.py

@@ -98,8 +98,7 @@ def get_completeness_matrix(tcat, ctab, mbinw, ybinw):
 
     return oin, out, cmags, cyeas
 
-def get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year,
-                      info=False):
+def get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, last_year, info=False):
     """
     :param aval: 
         GR a-value
@@ -341,7 +340,8 @@ def get_norm_optimize_c(cat, agr, bgr, compl, last_year, ref_mag, mmax=None, bin
     mvals = np.arange(ref_mag, mmax+binw/10, binw)
     rates = list(10**(agr-bgr * mvals[:-1]) - 10**(agr - bgr * mvals[1:]))
 
-    pocc = rates / sum(rates)
+
+    #pocc = rates / sum(rates)
 
     #prob = 1
     # If using log (and not multiplicative) set initial prob to 0
@@ -429,7 +429,7 @@ def get_norm_optimize_poisson(cat, agr, bgr, compl, last_year, mmax=None, binw=0
 
     rates = list(10**(agr-bgr * mvals[:-1]) - 10**(agr - bgr * mvals[1:]))
 
-    pocc = rates / sum(rates)
+    #pocc = rates / sum(rates)
 
     # Using log (and not multiplicative) set initial prob to 0
     prob = 0
@@ -495,7 +495,7 @@ def get_norm_optimize_d(cat, agr, bgr, compl, last_year, mmax=None, binw=0.1):
 
     rates = list(10**(agr-bgr * mvals[:-1]) - 10**(agr - bgr * mvals[1:]))
     #print(rates)
-    pocc = rates / sum(rates)
+    #pocc = rates / sum(rates)
 
     # Using log (and not multiplicative) set initial prob to 0
     prob = 0

openquake/cat/isf_catalogue.py

@@ -1205,13 +1205,14 @@ def get_origin_mag_tables(self):
                 else:
                     print('Location:', orig.location.depthSolution)
                     raise ValueError("Unsupported case")
-
+                    
                 if (orig.location.depth == 'None' or
                         orig.location.depth == '' or
                         orig.location.depth is None):
                     depth = np.NaN
-                elif orig.location.depth < 0.01:
+                elif orig.location.depth <= 0.0:
                     depth = orig.location.depth
+                    print('Depth:', orig.location.depth)
                     fmt = "Warning, depth <= 0.0 (id:{:s})"
                     warnings.warn(fmt.format(eq.id))
                 elif orig.location.depth:

openquake/cat/tests/completeness/analysis_rates_test.py

@@ -169,3 +169,124 @@ def test_filter_completeness(self):
         assert min(set(magsA)) >= cref[0][1] - 1.0
         assert max(set(magsB)) <= cref[1][1] + 1.0
         assert min(set(magsB)) >= cref[1][1] - 1.0
+
+class ComputeGRParametersTest_sim2(unittest.TestCase):
+    """ 
+    Tests for completeness and FMD parameters for a synthetic catalogue
+    using the optimize criteria
+    Simulated with completeness [[1900, 7.0], [1960, 5.0]]
+    with b = 1 and a = 7
+    """
+
+    def setUp(self):
+
+        # Temp folder
+        tmp_folder = tempfile.mkdtemp()
+
+        # Folder with the catalogue
+        self.fname_input_pattern = os.path.join(DATA_PATH, 'examp_1960M5.csv')
+        ref_config = os.path.join(DATA_PATH, 'examp_config.toml')
+
+        # Load the config template
+        conf_txt = toml.load(ref_config)
+
+        # Create the config file for the test
+        self.fname_config = os.path.join(tmp_folder, 'config.toml')
+        with open(self.fname_config, 'w', encoding='utf-8') as tmpf:
+            toml.dump(conf_txt, tmpf)
+
+        # Output folder
+        self.folder_out = tmp_folder
+
+        # Create completeness files
+        get_completenesses(self.fname_config, self.folder_out)
+
+    def test_compute_gr_param(self):
+        """ Testing the calculation """
+
+        conf = toml.load(self.fname_config)
+        print(conf)
+        completeness_analysis(self.fname_input_pattern,
+                              self.fname_config,
+                              self.folder_out,
+                              self.folder_out,
+                              self.folder_out)
+
+        # Load updated configuration file
+        conf = toml.load(self.fname_config)
+        print(conf)
+
+        # Tests
+
+        expected = 7.051
+        computed = conf['sources']['1960M5']['agr_weichert']
+        self.assertAlmostEqual(computed, expected, msg='aGR', delta = 0.2)
+
+        expected = 1.0118
+        computed = conf['sources']['1960M5']['bgr_weichert']
+        self.assertAlmostEqual(computed, expected, msg='bGR', delta = 0.1)
+
+        expected = 5.0
+        computed = conf['sources']['1960M5']['rmag']
+        self.assertAlmostEqual(computed, expected, msg='rmag', places=5)
+
+
+class ComputeGRParametersTest_Poisson_v1(unittest.TestCase):
+    """ 
+    Tests for completeness and FMD parameters for a synthetic catalogue
+    using the poisson criteria
+    Simulated with completeness [[1900, 7.0], [1960, 5.0]]
+    with b = 1 and a = 7
+    """
+
+    def setUp(self):
+
+        # Temp folder
+        tmp_folder = tempfile.mkdtemp()
+
+        # Folder with the catalogue
+        self.fname_input_pattern = os.path.join(DATA_PATH, 'examp_1960M5.csv')
+        ref_config = os.path.join(DATA_PATH, 'examp_config_poisson.toml')
+
+        # Load the config template
+        conf_txt = toml.load(ref_config)
+
+        # Create the config file for the test
+        self.fname_config = os.path.join(tmp_folder, 'config.toml')
+        with open(self.fname_config, 'w', encoding='utf-8') as tmpf:
+            toml.dump(conf_txt, tmpf)
+
+        # Output folder
+        self.folder_out = tmp_folder
+
+        # Create completeness files
+        get_completenesses(self.fname_config, self.folder_out)
+
+    def test_compute_gr_param(self):
+        """ Testing the calculation """
+
+        conf = toml.load(self.fname_config)
+
+        completeness_analysis(self.fname_input_pattern,
+                              self.fname_config,
+                              self.folder_out,
+                              self.folder_out,
+                              self.folder_out)
+
+        # Load updated configuration file
+        conf = toml.load(self.fname_config)
+        print(conf)
+
+        # Tests
+
+        expected = 7.0
+        computed = conf['sources']['1960M5']['agr_weichert']
+        self.assertAlmostEqual(computed, expected, msg='aGR', delta = 0.2)
+
+        expected = 1.0000
+        computed = conf['sources']['1960M5']['bgr_weichert']
+        self.assertAlmostEqual(computed, expected, msg='bGR', delta = 0.1)
+
+        expected = 5.0
+        computed = conf['sources']['1960M5']['rmag']
+        self.assertAlmostEqual(computed, expected, msg='rmag', places=5)