Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Adding tests + improving apriori conditions #287

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Adding tests + improving apriori conditions #287

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b8b662f
Adding tests + improving apriori conditions
mmpagani Feb 1, 2023
c66a6e1
More work on the apriori constraints
mmpagani Feb 2, 2023
2e662c3
Fixing the use of the prefix for multipoint sources
mmpagani Feb 2, 2023
c06a2e5
Slight changes to the declustering
mmpagani Feb 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
121 changes: 95 additions & 26 deletions openquake/cat/completeness/generate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,21 @@
import numpy as np
import multiprocessing
from openquake.wkf.utils import create_folder
from openquake.cat.completeness.analysis import clean_completeness


def get_compl(perms, mags, years, idxperm, flexible):
"""
example: get_compl(perms, mags, years, 0, True)
"""
if flexible:
idx = np.where(perms[idxperm, :] == len(mags) + 1, -1, perms)
else:
idx = perms[idxperm, :]
idx = idx.astype(int)
tmp = np.array([(y, m) for y, m in zip(years, mags[idx])])
ctab = clean_completeness(tmp)
return ctab


def mm(a):
Expand All @@ -52,24 +67,71 @@ def get_completenesses(fname_config, folder_out):
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', {})
# This is left for back compatibility
apriori_conditions_out = config[key].get('apriori_conditions', {})
apriori_conditions_out = config[key].get('apriori_conditions_out',
apriori_conditions_out)
apriori_conditions_in = config[key].get('apriori_conditions_in', {})

print(apriori_conditions_out)
print(apriori_conditions_in)

step = config[key].get('step', 8)
flexible = config[key].get('flexible', False)

_get_completenesses(mags, years, folder_out, num_steps, min_mag_compl,
apriori_conditions, step, flexible)
apriori_conditions_out, apriori_conditions_in, step,
flexible)


def _get_completenesses(mags, years, folder_out=None, num_steps=0,
min_mag_compl=None, apriori_conditions={},
step=6, flexible=False):
min_mag_compl=None, apriori_conditions_out={},
apriori_conditions_in={}, step=6, flexible=False):
"""
This function generates a set of completeness windows starting from a
set of magnitude and time thresholds included in the `mags` and `years`
vectors.

The output is a triple

:param mags:
A list or numpy array in increasing order
:param years:
A list or numpy array in increasing order
:param folder_out:
A string with the output folder. When `None` the code does not save
the produced output.
:param num_steps:
The minimum number of steps that each completeness window must have.
Zero corresponds to a flat completeness window.
:param min_mag_compl:
The minimum value of magnitude used to define the completeness. Note
that this is also the minimum magnitude that will be used to compute
the magnitude-frequency distribution.
:param apriori_conditions:
This is a dictionary where keys are years and values are magnitudes.
Each year-magnitude combination defines points (i.e. earthquakes)
that must be left out by the selected completeness windows.
:param step:
The step used to parallelize the calculation [default is 6]
:param flexible:
A boolean. When False the completeness window is always starting on
the right with the smallest magnitude value
:returns:
A triple of :class:`numpy.ndarray` instances. The first one, `perms`,
is an array of size <number of perms> x <number of time intervals>.
The array contains integer indexes indicating the magnitude value
of completeness for the corresponding time period. For example,
if `years` is [1900, 2000], `mags` is [3, 4, 5], the row of
perms corresponding to [0, 2], states that the catalogue is complete
above 5 since 1900 and above 3 since 2000.
"""

if isinstance(mags, list):
mags = np.array(mags)
if isinstance(years, list):
years = np.array(years)

msg = 'Years must be in ascending order'
assert np.all(np.diff(years) > 0), msg
msg = 'Mags must be in ascending order'
Expand All @@ -87,16 +149,18 @@ def _get_completenesses(mags, years, folder_out=None, num_steps=0,
if min_mag_compl is None:
min_mag_compl = min(mags)

# Check if the catalogue contains magnitudes above the minimum threshold
if len(np.where(min_mag_compl <= mags)) < 1:
msg = 'None of the magnitude intervals above the min_mag_compl'
raise ValueError(msg)
max_first_idx = np.min(np.where(min_mag_compl <= mags))

# Info
print('Total number of combinations : {:,d}'.format(len(mags)**len(years)))
print(f'Index of first magnitude : {max_first_idx:,d}')
tmp = len(mags)**len(years)
print('Total number of combinations : {:,d}'.format(tmp))
print(f'Index of the min magnitude of completeness : {max_first_idx:,d}')

# Creating the possible completenesses
# Creating the all the possible completenes windows
perms = []
for y in [years[i:min(i+step, len(years))] for i in range(0, len(years),
step)]:
Expand All @@ -108,11 +172,7 @@ def _get_completenesses(mags, years, folder_out=None, num_steps=0,
# with time
p = p[np.diff(p, axis=1).min(axis=1) >= 0, :]

# Selecting combinations within min_mag_compl
# if max(years) in y:
# p = p[p[:, 0] <= max_first_idx, :]

# Updating
# Update the list with the permutations
if len(perms):
new = []
for x in perms:
Expand All @@ -122,7 +182,7 @@ def _get_completenesses(mags, years, folder_out=None, num_steps=0,
else:
perms = p

# Full set of possible completeness windows. Each perm in perms is a
# Full set of possible completeness windows. Each perm in `perms` is a
# list of indexes for the magnitudes in the mags vector (in increasing
# order), starting from the most recent time interval. So 0 in the
# first position means the lowest value of magnitude for the most
Expand All @@ -142,24 +202,33 @@ def _get_completenesses(mags, years, folder_out=None, num_steps=0,
# lower or equal than a threshold
perms = perms[perms[:, 0] <= max_first_idx, :]

# Applying a-priori conditions
for yea_str in apriori_conditions.keys():
yea = float(yea_str)
mag = float(apriori_conditions[yea_str])
idx_yea = np.minimum(np.max(np.where(years >= yea)), len(years)-1)
idx_mag = np.max(np.where(mag >= mags))
# idxs = np.minimum(perms[:, idx_yea], len(mags)-1-dlt)
# perms = perms[idxs <= idx_mag, :]
perms = perms[perms[:, idx_yea] <= idx_mag, :]

print(f'Total number selected : {len(perms):,d}')

# Replacing the index of the 'buffer;' magnitude
if flexible:
perms = np.where(perms == len(mags), -1, perms)

# Applying a-priori conditions
print(f'Total number selected completeness windows : {len(perms):,d}')
for yea_str in apriori_conditions_out.keys():
yea = float(yea_str)
mag = float(apriori_conditions_out[yea_str])
idx_yea = np.min(np.where(yea >= years))
idx_mag = np.min(np.where(mag < mags))
perms = perms[perms[:, idx_yea] >= idx_mag, :]
print(f'Total number selected completeness windows : {len(perms):,d}')

# Applying a-priori conditions IN
for yea_str in apriori_conditions_in.keys():
yea = float(yea_str)
mag = float(apriori_conditions_in[yea_str])
idx_yea = np.min(np.where(yea >= years))
#idx_mag = np.max(np.where(mags > mag))
idx_mag = np.min(np.where(mags > mag))
perms = perms[perms[:, idx_yea] < idx_mag, :]

print(f'Total number selected completeness windows : {len(perms):,d}')

if folder_out is not None:
print(f'Saving completeness tables in: {folder_out:s}')
print(f'Saving completeness tables in : {folder_out:s}')
np.save(os.path.join(folder_out, 'dispositions.npy'), perms)
np.save(os.path.join(folder_out, 'mags.npy'), mags)
np.save(os.path.join(folder_out, 'years.npy'), years)
Expand Down
44 changes: 36 additions & 8 deletions openquake/cat/completeness/plot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,23 +59,36 @@ def get_xy(ctab, ymax=2015, rndx=0.0, rndy=0.0):
return np.array([(a, b) for a, b in zip(x, y)])


def main(folder, ymin=1900, ymax=2020, mmin=4.0, mmax=8.0):

perms = np.load(os.path.join(folder, 'dispositions.npy'))
mags = np.load(os.path.join(folder, 'mags.npy'))
years = np.load(os.path.join(folder, 'years.npy'))
def plot_completeness(
perms, mags, years, ymin=1900, ymax=2020, mmin=4.0, mmax=8.0, fname='',
apriori={}, apriori_in={}):
"""
Plots the set of completeness windows defined in by the `perms`, `mags`,
`years` triple.
"""

fig, ax = plt.subplots()
lines = []
colors = np.random.rand(len(perms), 3)

labls = []
for i, prm in enumerate(perms):
idx = prm.astype(int)
tmp = np.array([(y, m) for y, m in zip(years, mags[idx])])
ctab = clean_completeness(tmp)
coo = get_xy(ctab, rndx=0.0, rndy=0.1)
# ax.plot(coo[:, 0], coo[:, 1], color=colors[i, :])
lines.append(tuple([(x, y) for x, y in zip(coo[:, 0], coo[:, 1])]))
labls.append(f'{i}')

for key in apriori:
mag = float(apriori[key])
yea = float(key)
plt.plot(yea, mag, 'o', color='red')

for key in apriori_in:
mag = float(apriori_in[key])
yea = float(key)
print(mag, yea)
plt.plot(yea, mag, 'o', markersize=10., color='green')

coll = matplotlib.collections.LineCollection(lines, colors=colors)
ax.add_collection(coll)
Expand All @@ -86,10 +99,25 @@ def main(folder, ymin=1900, ymax=2020, mmin=4.0, mmax=8.0):
ax.set_xlim(ymin, ymax)
ax.set_ylim(mmin, mmax)

plt.savefig('completeness.png')
if len(fname):
plt.savefig(fname)
plt.show()


def main(folder, ymin=1900, ymax=2020, mmin=4.0, mmax=8.0):
"""
Plots the set of completeness windows defined by the `perms`, `mags`,
`years` triple files in a given folder.
"""

perms = np.load(os.path.join(folder, 'dispositions.npy'))
mags = np.load(os.path.join(folder, 'mags.npy'))
years = np.load(os.path.join(folder, 'years.npy'))

fname = 'completeness.png'
plot_completeness(perms, mags, years, ymin, ymax, mmin, mmax, fname)


main.folder = 'Folder containing the completeness files'
main.ymin = 'Minimum year in the plot'
main.ymax = 'Maximum year in the plot'
Expand Down
69 changes: 69 additions & 0 deletions openquake/cat/tests/completeness/completeness_matrix_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
# ------------------- The OpenQuake Model Building Toolkit --------------------
# Copyright (C) 2022 GEM Foundation
# _______ _______ __ __ _______ _______ ___ _
# | || | | |_| || _ || || | | |
# | _ || _ | ____ | || |_| ||_ _|| |_| |
# | | | || | | ||____|| || | | | | _|
# | |_| || |_| | | || _ | | | | |_
# | || | | ||_|| || |_| | | | | _ |
# |_______||____||_| |_| |_||_______| |___| |___| |_|
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Affero General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more
# details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# -----------------------------------------------------------------------------
# vim: tabstop=4 shiftwidth=4 softtabstop=4
# coding: utf-8

import unittest
import numpy as np
from openquake.hmtk.seismicity.catalogue import Catalogue
from openquake.hmtk.seismicity.occurrence.utils import get_completeness_counts
from openquake.mbt.tools.model_building.dclustering import _add_defaults
from openquake.cat.completeness.norms import get_completeness_matrix


class CompletenessMatrixTest(unittest.TestCase):

def setUp(self):
dat = [[1900, 6.0],
[1980, 6.0],
[1970, 5.0],
[1980, 5.0],
[1990, 5.0]]
dat = np.array(dat)
cat = Catalogue()
cat.load_from_array(['year', 'magnitude'], dat)
cat = _add_defaults(cat)
cat.data["dtime"] = cat.get_decimal_time()
self.cat = cat
self.compl = np.array([[1980, 5.0], [1950, 5.9]])

def test_case01(self):
binw = 0.5
oin, out, cmags, cyeas = get_completeness_matrix(self.cat, self.compl,
0.5, 10.0)
oin_expected = np.array([[-1., -1., -1., -1., -1., -1., -1., -1., 1., 1.],
[-1., -1., -1., -1., -1., 0., 0., 0., 0., 0.],
[-1., -1., -1., -1., -1., 0., 0., 0., 1., 0.]])
out_expected = np.array([[ 0., 0., 0., 0., 0., 0., 0., 1., -1., -1.],
[ 0., 0., 0., 0., 0., -1., -1., -1., -1., -1.],
[ 1., 0., 0., 0., 0., -1., -1., -1., -1., -1.]])
np.testing.assert_array_equal(oin, oin_expected)
np.testing.assert_array_equal(out, out_expected)

# Check the consistency of results with the ones provided by the
# completeness count
cmag, t_per, n_obs = get_completeness_counts(self.cat, self.compl,
binw)
oin[oin < 0] = 0.0
np.testing.assert_array_equal(np.sum(oin, axis=1), n_obs)
Loading