Skip to content

Commit

Permalink
adding test for full mfd workflow to ensure with same random seed res…
Browse files Browse the repository at this point in the history
…ults are the same
  • Loading branch information
Kendra Johnson authored and Kendra Johnson committed Apr 3, 2024
1 parent e8c19b4 commit 27ca787
Show file tree
Hide file tree
Showing 9 changed files with 123 additions and 12 deletions.
43 changes: 34 additions & 9 deletions openquake/cat/completeness/mfd_eval_plots.py
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,15 @@ def _make_a_b_histos(df_all, df_best, figsdir):
ax[1].tick_params(axis='y', labelcolor=color)

for ii, idx in enumerate(set(df_all.idx)):

df_sub = df_all[df_all.idx == idx]
ax[0].hist(df_sub.agr, bins=binsA, color='gray', alpha=10/num_cats)
ax[1].hist(df_sub.bgr, bins=binsB, color='gray', alpha=10/num_cats)
if num_cats < 10:
alpha = 0.1
else:
alpha = 10/num_cats

ax[0].hist(df_sub.agr, bins=binsA, color='gray', alpha=alpha)
ax[1].hist(df_sub.bgr, bins=binsB, color='gray', alpha=alpha)
ax2a = ax[0].twinx()
ax2b = ax[1].twinx()
ax2a.hist(df_best.agr, bins=binsA, color='red', alpha=0.2)
Expand Down Expand Up @@ -159,13 +165,18 @@ def plot_best_mfds(df_best, figsdir):
linewidth=0.1, zorder=0, label='Cumulative MFD')

else:
if num <= 10:
alpha1 = 0.1
else:
alpha1 = 10/num

plt.scatter(row.mags, row.rates, marker='_', color='r',
alpha=10/num)
alpha=alpha1)
plt.scatter(row.mags, row.cm_rates, marker='.', color='b',
alpha=10/num)
plt.semilogy(mfd_m, mfd_r, color='r', alpha=30/num, linewidth=0.1,
alpha=alpha1)
plt.semilogy(mfd_m, mfd_r, color='r', alpha=3*alpha1, linewidth=0.1,
zorder=0)
plt.semilogy(mfd_m, mfd_cr, color='b', alpha=50/num,
plt.semilogy(mfd_m, mfd_cr, color='b', alpha=5*alpha1,
linewidth=0.1, zorder=0)

plt.xlabel('Magnitude')
Expand Down Expand Up @@ -273,6 +284,8 @@ def plot_weighted_covariance_ellipse(df, figdir, n_std=1.0,
plt.savefig(fout, dpi=300)
plt.close()

return center_x, center_y, major_x1, major_y1, major_x2, major_y2



def plot_all_mfds(df_all, df_best, figsdir, field='rates', bins=10, bw=0.2, figname=None):
Expand Down Expand Up @@ -314,11 +327,16 @@ def plot_all_mfds(df_all, df_best, figsdir, field='rates', bins=10, bw=0.2, fign
mgrts = mfd.get_annual_occurrence_rates()
mfd_m = [m[0] for m in mgrts]
mfd_r = [m[1] for m in mgrts]
if len(df_best) <=30:
alpha = 0.1
else:
alpha=30/len(df_best)

if field == 'cm_rates':
mfd_cr = [sum(mfd_r[ii:]) for ii in range(len(mfd_r))]
plt.semilogy(mfd_m, mfd_cr, color='maroon', linewidth=0.2, zorder=10, alpha=30/len(df_best))
plt.semilogy(mfd_m, mfd_cr, color='maroon', linewidth=0.2, zorder=10, alpha=alpha)
else:
plt.semilogy(mfd_m, mfd_r, color='maroon', linewidth=0.2, zorder=10, alpha=30/len(df_best))
plt.semilogy(mfd_m, mfd_r, color='maroon', linewidth=0.2, zorder=10, alpha=alpha)

if figname==None:
figname = f'mfds_all_{field}.png'
Expand All @@ -334,6 +352,7 @@ def plot_all_mfds(df_all, df_best, figsdir, field='rates', bins=10, bw=0.2, fign


def make_all_plots(resdir_base, compdir, figsdir_base, labels):
agrs, bgrs, labs = [], [], []
for label in labels:
print(f'Running for {label}')
resdir = os.path.join(resdir_base, label)
Expand Down Expand Up @@ -368,5 +387,11 @@ def make_all_plots(resdir_base, compdir, figsdir_base, labels):
plot_all_mfds(df_all, df_thresh, figsdir, field='cm_rates', bins=60,
figname='mfds_thresh_cmrates.png')
print('plotting covariance')
plot_weighted_covariance_ellipse(df_best, figsdir)
cx, cy, mx1, my1, mx2, my2 = plot_weighted_covariance_ellipse(df_best, figsdir)
plot_weighted_covariance_ellipse(df_thresh, figsdir, figname=f'{label}-20percent.png')

labs.extend([f'{label}-center', f'{label}-low', f'{label}-high'])
agrs.extend([cx, mx1, mx2])
bgrs.extend([cy, my1, my2])

return labs, agrs, bgrs
2 changes: 1 addition & 1 deletion openquake/mbi/cat/MFDs_sample_mag_sigma.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@
import numpy as np
from openquake.baselib import sap

from openquake.mbt.mfds_sample.make_mfds import make_many_mfds
from openquake.mbt.tools.mfd_sample.make_mfds import make_many_mfds


def main(fname_config):
Expand Down
21 changes: 21 additions & 0 deletions openquake/mbt/tests/tools/mfd_sample/config/compl.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
name = "CAM test"
mmin = 4.0
bin_width = 0.2
rupture_mesh_spacing = 2.5
kernel_maximum_distance = 120.0
kernel_smoothing = [ [ 0.8, 20.0,], [ 0.15, 40.0,], [ 0.05, 80.0,],]
source_prefix = "L01"

[completeness]
mags = []
years = []
completeness_ref = [ [ 1980.0, 4.0,], [ 1964.0, 5.0,], [ 1904.0, 7.0,],]
deviation = 1.7
ref_mag = 2.0
ref_upp_mag = 10.0
min_mag_compl = 4.0
bmin = 0.7
bmax = 1.3
flexible = true
num_steps = 2
optimization_criterion = "optimize"
16 changes: 16 additions & 0 deletions openquake/mbt/tests/tools/mfd_sample/config/dec.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
[main]

catalogue = './data2/sub_cat.pkl'
tr_file = './data2/classified_test.hdf5'

create_subcatalogues = 'true'
catalogue_add_defaults = 'true'

[method1]
name = 'GardnerKnopoffType1'
params = {'time_distance_window' = 'UhrhammerWindow', 'fs_time_prop' = 0.1}
label = 'UH'

[case1]
regions = ['crustal', 'intA']
label = 'int_cru'
23 changes: 23 additions & 0 deletions openquake/mbt/tests/tools/mfd_sample/config/test.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
[main]
catalogue_filename = './data2/sub_cat.pkl'
output_directory = 'test'

[catalogues]
number = 3
sample_type = 'random'
create_catalogues = true

[decluster]
decluster_catalogues = true
decluster_settings = 'config/dec.toml'

[completeness]
generate_completeness = true
completeness_settings = 'config/compl.toml'

[mfds]
create_mfds = true
labels = ['intA','crustal']

[plot]
make_plots = true
Binary file not shown.
Binary file not shown.
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
,label,a-values,b-values
0,intA-center,3.4221566716867344,0.8477886165624219
1,intA-low,3.3341265324448983,0.8035538508320615
2,intA-high,3.5101868109285705,0.8920233822927822
3,crustal-center,3.350288725907407,0.8778145185036605
4,crustal-low,3.263790223525028,0.8485030837041351
5,crustal-high,3.436787228289786,0.907125953303186
23 changes: 21 additions & 2 deletions openquake/mbt/tests/tools/mfd_sample/mfd_sample_test.py
Original file line number Diff line number Diff line change
@@ -1,14 +1,16 @@
import os
import shutil
import unittest
import tempfile
import numpy as np
import pandas as pd

from openquake.mbt.tools.mfd_sample.make_mfds import _create_catalogue_versions
from openquake.mbt.tools.mfd_sample.make_mfds import (_create_catalogue_versions,
make_many_mfds)

BASE_PATH = os.path.dirname(__file__)

class TestWorkflow(unittest.TestCase):
class TestGenCats(unittest.TestCase):

def setUp(self):
self.catfi = os.path.join(BASE_PATH, 'data', 'catalogue.csv')
Expand All @@ -27,3 +29,20 @@ def test_generate_cats(self):
expected = pd.read_csv(mags_exp_fi)
output = pd.read_csv(mags_out_fi)
assert expected.equals(output)

shutil.rmtree(self.tmpd)

class TestWorkflow(unittest.TestCase):

def test_full_wkflow(self):
config_fi = os.path.join(BASE_PATH, 'config', 'test.toml')
make_many_mfds(config_fi)

expected_fi = os.path.join(BASE_PATH, 'expected2', 'mfd-results.csv')
output_fi = os.path.join(BASE_PATH, 'test', 'mfd-results.csv')
expected = pd.read_csv(expected_fi)
output = pd.read_csv(output_fi)
assert expected.equals(output)

shutil.rmtree('test')
os.remove('tmp-config-dcl.toml')

0 comments on commit 27ca787

Please sign in to comment.