Merge pull request #402 from GEMScienceTools/fix_samplign

[WIP] Fix sampling of profiles

openquake/cat/completeness/analysis.py

@@ -177,7 +177,8 @@ def check_criterion(criterion, rate, previous_norm, tvars):
 
     elif criterion == 'optimize':
         tmp_rate = -1
-        norm = get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, info=False)
+        norm = get_norm_optimize(tcat, aval, bval, ctab, cmag, n_obs, t_per, 
+                                 last_year, info=False)
 
     elif criterion == 'optimize_a':
         tmp_rate = -1
@@ -309,20 +310,10 @@ 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
 
-        #tmp = []
-        #for yea, j in zip(years, prm):
-        #    if j >= -1e-10:
-        #        tmp.append([yea, mags[int(j)]])
-        #tmp = np.array(tmp)
-
-        #if len(tmp) > 0:
-        #    ctab = clean_completeness(tmp)
-        #else:
-        #    continue
-
         # Check compatibility between catalogue and completeness table. This
         # function finds in each magnitude interval defined in the completeness
         # table the earliest year since the occurrence of a number X of
@@ -467,8 +458,42 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag,
 
     return save
 
+def read_compl_params(config):
+    """
+    """
+
+    # Read parameters for completeness analysis
+    key = 'completeness'
+    ms = np.array(config[key]['mags'], dtype=float)
+    yrs = np.array(config[key]['years'])
+    bw = config.get('bin_width', 0.1)
+    r_m = config[key].get('ref_mag', 5.0)
+    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
+
 
-def completeness_analysis(fname_input_pattern, fname_config, folder_out_figs,
+def read_compl_data(folder_in):
+    """
+    """
+    # Reading completeness data
+    print(f'Reading completeness data from: {folder_in:s}')
+    fname_disp = os.path.join(folder_in, 'dispositions.npy')
+    perms = np.load(fname_disp)
+    mags_chk = np.load(os.path.join(folder_in, 'mags.npy'))
+    years_chk = np.load(os.path.join(folder_in, 'years.npy'))
+    compl_tables = {'perms': perms, 'mags_chk': mags_chk,
+                    'years_chk': years_chk}
+
+    return compl_tables, mags_chk, years_chk
+
+
+
+def completeness_analysis(fname_input_pattern, f_config, folder_out_figs,
                           folder_in, folder_out, skip=''):
     """
     :param fname_input_pattern:
@@ -486,36 +511,18 @@ def completeness_analysis(fname_input_pattern, fname_config, folder_out_figs,
     """
 
     # Loading configuration
-    config = toml.load(fname_config)
+    config = toml.load(f_config)
 
-    # Read parameters for completeness analysis
-    key = 'completeness'
-    mags = np.array(config[key]['mags'])
-    years = np.array(config[key]['years'])
-    binw = config.get('bin_width', 0.1)
-    ref_mag = config[key].get('ref_mag', 5.0)
-    ref_upp_mag = 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'
-    criterion = config[key].get('optimization_criterion', 'optimize')
-    print(criterion)
+    ms, yrs, bw, r_m, r_up_m, bmin, bmax, crit = read_compl_params(config)
 
-    # Reading completeness data
-    print(f'Reading completeness data from: {folder_in:s}')
-    fname_disp = os.path.join(folder_in, 'dispositions.npy')
-    perms = np.load(fname_disp)
-    mags_chk = np.load(os.path.join(folder_in, 'mags.npy'))
-    years_chk = np.load(os.path.join(folder_in, 'years.npy'))
-    compl_tables = {'perms': perms, 'mags_chk': mags_chk,
-                    'years_chk': years_chk}
+    compl_tables, mags_chk, years_chk = read_compl_data(folder_in)
 
     # Fixing sorting of years
-    if np.all(np.diff(years)) >= 0:
-        years = np.flipud(years)
+    if np.all(np.diff(yrs)) >= 0:
+        yrs = np.flipud(yrs)
 
-    np.testing.assert_array_equal(mags, mags_chk)
-    np.testing.assert_array_equal(years, years_chk)
+    np.testing.assert_array_almost_equal(ms, mags_chk)
+    np.testing.assert_array_almost_equal(yrs, years_chk)
 
     # Info
     if len(skip) > 0:
@@ -537,8 +544,8 @@ def completeness_analysis(fname_input_pattern, fname_config, folder_out_figs,
         else:
             var = {}
 
-        res = _completeness_analysis(fname, years, mags, binw, ref_mag,
-                                     ref_upp_mag, [bmin, bmax], criterion,
+        res = _completeness_analysis(fname, yrs, ms, bw, r_m,
+                                     r_up_m, [bmin, bmax], crit,
                                      compl_tables, src_id,
                                      folder_out_figs=folder_out_figs,
                                      folder_out=folder_out,
@@ -563,6 +570,6 @@ def completeness_analysis(fname_input_pattern, fname_config, folder_out_figs,
         # Updating configuration
         config['sources'][src_id] = var
 
-    with open(fname_config, 'w', encoding='utf-8') as fou:
+    with open(f_config, 'w', encoding='utf-8') as fou:
         fou.write(toml.dumps(config))
-        print(f'Updated {fname_config:s}')
+        print(f'Updated {f_config:s}')

openquake/cat/completeness/generate.py

@@ -50,21 +50,22 @@ def get_completenesses(fname_config, folder_out):
     create_folder(folder_out)
     config = toml.load(fname_config)
     key = 'completeness'
-    mags = np.array(config[key]['mags'])
+    mags = np.array(config[key]['mags'], dtype=np.float32)
     years = np.array(config[key]['years'])
     num_steps = config[key].get('num_steps', 0)
     min_mag_compl = config[key].get('min_mag_compl', None)
     apriori_conditions = config[key].get('apriori_conditions', {})
     cref = config[key].get('completeness_ref', None)
     step = config[key].get('step', 8)
+    mrange = config['completeness'].get('deviation', 1.0) 
 
     _get_completenesses(mags, years, folder_out, num_steps, min_mag_compl,
-                        apriori_conditions, cref, step)
+                        apriori_conditions, cref, step, mrange)
 
 
 def _get_completenesses(mags, years, folder_out=None, num_steps=0,
                         min_mag_compl=None, apriori_conditions={},
-                        completeness_ref=None, step=6):
+                        completeness_ref=None, step=6, mrange=1.0):
     """
     :param mags:
         A list or numpy array in increasing order
@@ -194,21 +195,15 @@ def _get_completenesses(mags, years, folder_out=None, num_steps=0,
         mags_ref = [c[1] for c in completeness_ref]
         rem = []
         for iper, prm in enumerate(perms):
-            #tmp = []
-            #for yea, j in zip(years, prm):
-            #    if j >= -1e-10:
-            #        tmp.append([yea, mags[int(j)]])
-#
-#            tmp = np.array(tmp)
-#            ctab = clean_completeness(tmp)
+
             ctab = _make_ctab(prm, years, mags)
             for yr, mg in ctab:
                 if not yr in years_ref:
                     rem.append(iper)
                     continue
                 index = years_ref.index(yr)
                 mdiff = abs(mags_ref[index] - mg)
-                if mdiff > 1:
+                if mdiff > mrange:
                     rem.append(iper)
                     continue