# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2019, GEM Foundation
#
# OpenQuake 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.
#
# OpenQuake 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 OpenQuake. If not, see <http://www.gnu.org/licenses/>.
import time
import zlib
import copy
import json
import os.path
import pickle
import operator
import logging
import numpy
from openquake.baselib import parallel, performance, general, hdf5
from openquake.hazardlib import (
geo, nrml, source, sourceconverter, InvalidFile, calc)
from openquake.hazardlib.source_group import (
CompositeSourceModel, SourceGroup, get_unique)
from openquake.hazardlib.source.multi_fault import save_and_split
from openquake.hazardlib.lt import (
apply_uncertainties, check_correlated, get_bset_value,
restrict_sampling, sampling_dt, unc_subsets)
from openquake.hazardlib.contexts import get_unique_inverse
from openquake.hazardlib.valid import basename
TWO24 = 2**24
TWO32 = 2**32
# the calculations building the rates with RateMap\s, i.e. the ones
# attributing the rates to the core trt_smrs (core_trt_smrs)
CLASSICAL_MODES = ('classical', 'classical_risk', 'classical_damage',
'classical_bcr', 'disaggregation', 'preclassical')
U16 = numpy.uint16
U32 = numpy.uint32
F32 = numpy.float32
bybranch = operator.attrgetter('branch')
checksum = operator.attrgetter('checksum')
source_info_dt = numpy.dtype([
('source_id', hdf5.vstr), # 0
('grp_id', U16), # 1
('code', (numpy.bytes_, 1)), # 2
('calc_time', F32), # 3
('num_ctxs', numpy.uint64), # 4
('est_ctxs', numpy.uint64), # 5
('num_ruptures', U32), # 6
('weight', F32), # 7
('mul', U16), # 8
])
[docs]def sampling(samples, trt_smr):
"""
:returns: a structured array (trt_smr, samples) of length 1
"""
return numpy.array([(trt_smr, samples)], sampling_dt)
# NB: blocksize is chosen so that event_based/case_35 works
[docs]def splitMF(sources, disagg_by_src, blocksize=1000):
"""
Split the MultiPointSources to avoid oceanic sources with 160M ruptures
hanging during rupture sampling
"""
splits = []
for src in sources:
if src.code == b'M' and len(src) > blocksize:
for i, slc in enumerate(general.gen_slices(0, len(src), blocksize)):
segment = source.MultiPointSource(
source_id=f'{src.source_id}:{i}',
name=src.name,
tectonic_region_type=src.tectonic_region_type,
mfd=src.mfd[slc],
magnitude_scaling_relationship=(
src.magnitude_scaling_relationship),
rupture_aspect_ratio=src.rupture_aspect_ratio,
upper_seismogenic_depth=src.upper_seismogenic_depth,
lower_seismogenic_depth=src.lower_seismogenic_depth,
nodal_plane_distribution=src.nodal_plane_distribution,
hypocenter_distribution=src.hypocenter_distribution,
mesh=geo.Mesh(src.mesh.lons[slc], src.mesh.lats[slc]),
temporal_occurrence_model=src.temporal_occurrence_model)
segment.sampling = src.sampling
splits.append(segment)
elif src.code == b'F' and not disagg_by_src:
# use the colon convention only in absence of kendra-splitting
for segment in src:
segment.sampling = src.sampling
splits.append(segment)
else:
splits.append(src)
sources[:] = splits
[docs]def check_unique(ids, msg='', strict=True):
"""
Raise a DuplicatedID exception if there are duplicated IDs
"""
if isinstance(ids, dict): # ids by key
all_ids = sum(ids.values(), [])
unique, counts = numpy.unique(all_ids, return_counts=True)
for dupl in unique[counts > 1]:
keys = [k for k in ids if dupl in ids[k]]
if keys:
errmsg = '%r appears in %s %s' % (dupl, keys, msg)
if strict:
raise nrml.DuplicatedID(errmsg)
else:
logging.info('*' * 60 + ' DuplicatedID:\n' + errmsg)
return
unique, counts = numpy.unique(ids, return_counts=True)
for u, c in zip(unique, counts):
if c > 1:
errmsg = '%r appears %d times %s' % (u, c, msg)
if strict:
raise nrml.DuplicatedID(errmsg)
else:
logging.info('*' * 60 + ' DuplicatedID:\n' + errmsg)
[docs]def create_source_info(csm, h5):
"""
Creates source_info, trt_smrs, toms
"""
csm.source_info = csm.get_source_info()
num_srcs = len(csm.source_info)
# avoid hdf5 damned bug by creating source_info in advance
h5.create_dataset('source_info', (num_srcs,), source_info_dt)
[docs]def trt_smrs(src):
return tuple(src.trt_smrs)
[docs]def read_source_model(fname, branch, converter, applied, sample, monitor):
"""
:param fname: path to a source model XML file
:param branch: source model logic tree branch ID
:param converter: SourceConverter
:param applied: list of source IDs within applyToSources
:param sample: a string with the sampling factor (if any)
:param monitor: a Monitor instance
:returns: a SourceModel instance
"""
t0 = time.time()
[sm] = nrml.read_source_models([fname], converter)
sm.branch = branch
for sg in sm.src_groups:
if sample and not sg.atomic:
srcs = []
for src in sg:
if src.source_id in applied:
srcs.append(src)
else:
srcs.extend(calc.filters.split_source(src))
if srcs:
kept = [src for src in srcs if src.source_id in applied]
rand = general.random_filter(srcs, float(sample))
sg.sources = (kept + rand) or [srcs[0]]
else:
sg.sources = []
sm.rtime = time.time() - t0 # save the read time
return {fname: sm}
# NB:this is called after reduce_sources, so ";" is not added
# if the same source appears multiple times, i.e. len(srcs) == 1
[docs]def add_semicolons(src_groups):
"""
Add semicolons to differentiate different sources with the same source_id
and then sort the sources in each group by extended source_id
"""
sources = general.AccumDict(accum=[])
for sg in src_groups:
for src in sg:
sources[src.source_id].append(src)
for src_id, srcs in sources.items():
srcs = get_unique(srcs)
if len(srcs) > 1:
# happens in logictree/case_01/rup.ini
for i, src in enumerate(srcs):
src.source_id = '%s;%d' % (src.source_id, i)
for sg in src_groups:
sg.sources.sort(key=operator.attrgetter('source_id'))
# tested in logictree/case_05 with 9 variations of
# AreaSource 1 and SimpleFaultSource 2
[docs]def check_branchID(branchID):
"""
Forbids invalid characters .:; used in fragmentno
"""
if '.' in branchID:
raise InvalidFile('branchID %r contains an invalid "."' % branchID)
elif ':' in branchID:
raise InvalidFile('branchID %r contains an invalid ":"' % branchID)
elif ';' in branchID:
raise InvalidFile('branchID %r contains an invalid ";"' % branchID)
[docs]def check_duplicates(smdict, strict):
# check_duplicates in the same file
for sm in smdict.values():
srcids = []
for sg in sm.src_groups:
srcids.extend(src.source_id for src in sg)
if sg.src_interdep == 'mutex':
# mutex sources in the same group must have all the same
# basename, i.e. the colon convention must be used
basenames = set(map(basename, sg))
assert len(basenames) == 1, basenames
check_unique(srcids, 'in ' + sm.fname, strict)
# check duplicates in different files but in the same branch
# the problem was discovered in the DOM model
for branch, sms in general.groupby(smdict.values(), bybranch).items():
srcids = general.AccumDict(accum=[])
fnames = []
for sm in sms:
if isinstance(sm, nrml.GeometryModel):
# the section IDs are not checked since they not count
# as real sources
continue
for sg in sm.src_groups:
srcids[sm.fname].extend(src.source_id for src in sg)
fnames.append(sm.fname)
check_unique(srcids, 'in branch %s' % branch, strict=strict)
[docs]def save_read_times(dstore, source_models):
"""
Store how many seconds it took to read each source model file
in a table (fname, rtime)
"""
dt = [('fname', hdf5.vstr), ('rtime', float)]
arr = numpy.array([(sm.fname, sm.rtime) for sm in source_models], dt)
dstore.create_dset('source_model_read_times', arr)
[docs]def get_csm(oq, full_lt, dstore=None):
"""
Build a CompositeSourceModel without applying the uncertainties,
that are applied in the workers, see modified_groups.
"""
converter = sourceconverter.SourceConverter(
oq.investigation_time, oq.rupture_mesh_spacing,
oq.complex_fault_mesh_spacing, oq.width_of_mfd_bin,
oq.area_source_discretization, oq.minimum_magnitude,
oq.source_id,
discard_trts=[s.strip() for s in oq.discard_trts.split(',')],
floating_x_step=oq.floating_x_step,
floating_y_step=oq.floating_y_step,
source_nodes=oq.source_nodes,
infer_occur_rates=oq.infer_occur_rates,
filter_sourcecodes=oq.filter_sourcecodes)
full_lt.ses_seed = oq.ses_seed
logging.info('Reading the source model(s) in parallel')
# NB: the source models file must be in the shared directory
# NB: dstore is None in logictree_test.py
allargs = []
sdata = full_lt.source_model_lt.source_data
allpaths = set(full_lt.source_model_lt.info.smpaths)
dic = general.group_array(sdata, 'fname')
smpaths = []
ss = os.environ.get('OQ_SAMPLE_SOURCES')
applied = set()
for srcs in full_lt.source_model_lt.info.applytosources.values():
applied.update(srcs)
for fname, rows in dic.items():
path = os.path.abspath(
os.path.join(full_lt.source_model_lt.basepath, fname))
smpaths.append(path)
allargs.append((path, rows[0]['branch'], converter, applied, ss))
for path in allpaths - set(smpaths): # geometry models
allargs.append((path, '', converter, applied, ss))
smdict = parallel.Starmap(read_source_model, allargs,
h5=dstore if dstore else None).reduce()
parallel.Starmap.shutdown() # save memory
smdict = {k: smdict[k] for k in sorted(smdict)}
if dstore:
save_read_times(dstore, smdict.values())
check_duplicates(smdict, strict=oq.disagg_by_src)
found = add_bangs(smdict)
if found:
logging.info('Found different sources with same ID %s',
general.shortlist(found))
# checking ps_grid_spacing
pointlike_sources = 0
for sm in smdict.values():
for sg in sm.src_groups:
for src in sg:
if src.code in b'PAM':
pointlike_sources += 1
break
if (oq.strict and oq.mosaic_model and pointlike_sources and
'classical' in oq.calculation_mode and oq.ps_grid_spacing == 0
and not oq.sites and not oq.disagg_by_src):
raise InvalidFile(f'{oq.inputs["job_ini"]}: '
'missing ps_grid_spacing')
return build_csm(oq, full_lt, smdict, dstore)
[docs]def get_bset_values(full_lt, sources):
"""
:param full_lt: a FullLogicTree instance
:param sources: a SourceGroup or a list of sources of the same group
:returns: the dictionary of uncertainties to apply expected by
modified_groups, i.e. one entry for each set of realizations
with the same uncertainties; the uncertainties of a set are the
ones of its first realization
NB: only the uncertainties relevant for the given sources are
returned, since the logic tree can be huge and the dictionary is
sent to the workers
"""
ordinals = {trt_smrs[0] % TWO24 for src in sources
for trt_smrs in unc_subsets(src)}
return {ordinal: full_lt.get_bset_values(ordinal)
for ordinal in sorted(ordinals)}
[docs]def modified_groups(sources, bset_values):
"""
Apply the uncertainties to a group of sources built *without* them,
as needed by the workers computing the rates or the ruptures: this
is done one set of realizations at a time, i.e. one set with the
same uncertainties at a time (see build_groups).
:param sources: a SourceGroup or a list of sources of the same group
:param bset_values:
the uncertainties to apply, as returned by get_bset_values; it
can be empty, if there are no uncertainties at all
:returns:
a generator of (trt_smrs, group) pairs, one for each set of
realizations with the same uncertainties, with the uncertainties
applied and the sampling restricted to the set of realizations
"""
for trt_smrs, srcs in _subsets_by_unc(sources).items():
grp = _restricted_group(sources, srcs, trt_smrs)
# NB: the trt_smrs are trti * TWO24 + ordinal, see gen_groups
bvals = bset_values[trt_smrs[0] % TWO24] if bset_values else []
# NB: check=False since the group is a fragment of the original
# one (split by weight in the preclassical), so the correlated
# branchsets were already checked at build time, see build_groups
grp = apply_uncertainties(bvals, grp, check=False)
for src in grp:
# the sources are modified after the preclassical, so the
# cached geometry must be discarded; it depends on the
# occurrence rates (see PointSource.
# _get_max_rupture_projection_radius)
if hasattr(src, 'radius'):
del src.radius
yield trt_smrs, grp
def _subsets_by_unc(sources):
"""
:returns: a dictionary trt_smrs -> sources, i.e. the sources grouped
by set of realizations with the same uncertainties
"""
if getattr(sources, 'atomic', False):
# the sources of an atomic group are mutually exclusive (or belong
# to a cluster), so they must be kept together, see sample_cluster
# and cmakers_groups; the sets of realizations are the same for
# all of them, since they belong to the same source model
return {trt_smrs: list(sources)
for trt_smrs in unc_subsets(sources[0])}
subsets = {}
for src in sources:
for trt_smrs in unc_subsets(src):
subsets.setdefault(trt_smrs, []).append(src)
return subsets
def _restricted_group(sources, srcs, trt_smrs):
"""
:returns: a group with the given sources, i.e. copies of the sources
of `sources` with the sampling restricted to trt_smrs
"""
restricted = [restrict_sampling(src, trt_smrs) for src in srcs]
if hasattr(sources, 'sources'): # keep the attributes of the group
grp = copy.copy(sources)
grp.sources = restricted
else: # a plain list of sources, e.g. a block of sources
grp = SourceGroup(sources[0].tectonic_region_type, restricted)
return grp
[docs]def unc_signature(bset_values, src):
"""
:returns: a tuple identifying the uncertainties applied to src
"""
sig = []
for bset, value in bset_values:
ok, val = get_bset_value(bset, value, src)
if ok:
sig.append((bset.id, str(val)))
return tuple(sig)
[docs]def collect_sources(full_lt, rlz_groups):
"""
:param full_lt: a FullLogicTree instance
:param rlz_groups: an iterator of (sm_rlz, source group) pairs
:returns: a dictionary id(grp) -> [grp, {id(src): (src, pairs)}], where
pairs is a list of (trt_smr, samples, signature) tuples, one per
realization of the source
"""
dic = {}
for rlz, grp in rlz_groups:
trti = full_lt.trti.get(grp.trt, 0)
trt_smr = trti * TWO24 + rlz.ordinal
bset_values = full_lt.get_bset_values(rlz.ordinal)
# NB: the uncertainties are applied later, in the workers, on
# groups split by weight, so the correlated branchsets are checked
# here, where the groups are still whole
check_correlated(bset_values, grp)
# NB: the sources are collected by id(grp), since the group objects
# are shared by all the realizations selecting the same source model
# file (see gen_groups); the groups built from them are then merged
# by trt_smrs below. The dicts preserve the order of first
# appearance, making the groups and the sources inside them
# reproducible.
srcs = dic.setdefault(id(grp), [grp, {}])[1]
for src in grp:
sig = unc_signature(bset_values, src)
pairs = srcs.setdefault(id(src), (src, []))[1]
pairs.append((trt_smr, rlz.samples, sig))
return dic
[docs]def source_with_subsets(src, pairs, sigrows):
"""
:param src: a source
:param pairs: a list of (trt_smr, samples, signature) tuples, one per
realization of the source
:param sigrows: a list of rows describing the uncertainty signatures,
extended with the signatures of the source
:returns: a copy of the source with the attributes sampling,
bysrc_subsets and bysrc_unc set, i.e. the sets of realizations
with the same uncertainties
"""
arrays, sigdict = [], {}
for trt_smr, samples, sig in pairs:
arrays.append((trt_smr, samples))
sigdict.setdefault(sig, []).append(trt_smr)
new_src = copy.copy(src)
new_src.sampling = numpy.array(
sorted(arrays), sampling_dt) # sorted by trt_smr
# NB: the subsets are stored only if the uncertainties are not the
# same in all the realizations; a source with no uncertainties (or with
# the same uncertainties everywhere) keeps its sampling as it is
new_src.bysrc_subsets = [
numpy.array(sorted(t), U32) for t in sigdict.values()
] if len(sigdict) > 1 else []
# flag the sources which will be modified in the workers: they must
# not be split in the preclassical, since the splitting destroys the
# geometry (and the MFD of the fault sources)
new_src.bysrc_unc = any(sigdict) # NB: () means no uncertainty
# the uncertainty signatures of the source, i.e. the sets of
# realizations with the same uncertainties, which are the indices of
# rate attribution for its rates (see unc_subsets)
for sig, trt_smrs in sigdict.items():
sigrows.append(dict(source_id=src.source_id, realizations=len(arrays),
signature=dict(sig), count=len(trt_smrs)))
return new_src
[docs]def build_groups(full_lt, rlz_groups, oq, dstore=None):
"""
Build the source groups without applying the uncertainties, as needed
by the workers: there is one group per source group in the source
model files, and the sources keep the trt_smrs of all the realizations
they belong to, not just the ones with a given set of uncertainties, so
that the number of groups depends on the source models only.
:param dstore: a DataStore instance or None; if given, the indices of
rate attribution are saved in it (the event based calculations do
not use them)
The uncertainties to be applied in each realization are not known
until the workers, so the realizations with different uncertainties are
stored in the bysrc_subsets attribute of each source (see
unc_subsets), and the rates are computed and attributed one set at a
time.
NB: the same source_id can be used by different sources, i.e. in
different source models, so the sources are keyed by id(src) and not
by source_id; the ids are disambiguated at the end, by adding a
semicolon, see add_semicolons
"""
# NB: the uncertainty signatures are stored in the datastore, so that
# `oq check_input` can print them, see store_unc_signatures
sigrows = [] # dicts with keys source_id, realizations, signature, count
dic = collect_sources(full_lt, rlz_groups)
out, atomic, acc = [], [], general.AccumDict(accum=[])
for grp, srcs in dic.values():
new_srcs = [source_with_subsets(src, pairs, sigrows)
for src, pairs in srcs.values()]
if grp.atomic:
# the atomic groups are never merged with the other groups,
# since their sources must be computed together
new = copy.copy(grp)
new.sources = new_srcs
atomic.append(new)
else:
acc[grp.trt].extend(new_srcs)
if atomic:
logging.info('Found %d atomic groups', len(atomic))
# NB: the sources are grouped by trt_smrs and TOM, so that there is
# one cmaker for each set of realizations, see get_cmakers
red_sources = 0
for trt, sources in acc.items():
grps, red = _group_sources(trt, sources)
out.extend(grps)
red_sources += red
if red_sources:
logging.info('reduce_sources was called %d times', red_sources)
out.extend(atomic)
for grp in out:
splitMF(grp.sources, oq.disagg_by_src)
add_semicolons(out) # else sources with the same id are lost
store_unc_signatures(out, sigrows, full_lt, oq, dstore)
return out
[docs]def store_unc_signatures(groups, sigrows, full_lt, oq, dstore=None):
"""
Store the uncertainty signatures of the sources in the datastore, so
that they can be printed by `oq check_input` (see the check_input
command) and read back with `oq show unc_signatures`. NB: the
signatures are stored as JSON objects, so that the unc_signatures view
can display a column for each branchset.
:param groups: a list of SourceGroups built without the uncertainties
:param sigrows: a list of dicts with keys source_id, realizations,
signature (a dictionary branchset_id -> value) and count
:param full_lt: a FullLogicTree instance
:param oq: an OqParam instance
:param dstore: a DataStore instance or None
"""
if dstore is not None and any(row['signature'] for row in sigrows):
# NB: the signatures are stored only if there are uncertainties,
# i.e. only if the rates must be split in indices of rate
# attribution
dt = [('source_id', hdf5.vstr), ('realizations', int),
('signature', hdf5.vstr), ('count', int)]
data = [(row['source_id'], row['realizations'], json.dumps(
row['signature']), row['count']) for row in sigrows]
dstore.create_df('unc_signatures', numpy.array(data, dt))
if oq.calculation_mode not in CLASSICAL_MODES:
# the event based calculators attribute the rates to the
# realizations of the group, not to the subsets
return
core_trt_smrs = get_core_trt_smrs(groups)
if dstore is not None:
dstore.hdf5.save_vlen('core_trt_smrs', core_trt_smrs)
Gt = get_core_size(core_trt_smrs, full_lt)
if Gt >= TWO32:
# NB: the gids are stored as uint32 in the _rates dataset, so
# there cannot be more than 2**32 columns in the RateMap
raise ValueError(
'The core size Gt=%d is too large (the maximum is %d), '
'you must reduce the logic tree' % (Gt, TWO32))
logging.warning('Core size Gt=%d out of R=%d realizations',
Gt, full_lt.get_num_paths())
if dstore is not None and 'sitecol' in dstore:
# NB: the sites are read after the sources, so the size in bytes
# is logged only if they are already in the datastore
N = len(dstore['sitecol'])
imtls = oq.imtls
L = imtls.size if hasattr(imtls, 'size') else len(imtls)
logging.warning('Global RateMap of %s for %d sites and %d levels',
general.humansize(4 * N * L * Gt), N, L)
[docs]def get_core_trt_smrs(groups):
"""
:param groups: a list of SourceGroups built without applying the
uncertainties (i.e. the src_groups of a CompositeSourceModel)
:returns: a sorted list of trt_smrs, the units of rate attribution
(to be stored as an hdf5.vuint32 array)
The uncertainties are applied in the workers, so the realizations
with different uncertainties are not separated at build time (see
build_groups). The rates are nevertheless computed separately for each
set of uncertainties and must be attributed to the right
realizations, hence this extra list; a rate is attributed to the unit
containing its trt_smrs. NB: the subsets are deduplicated, i.e. groups
with the same trt_smrs (for instance the same trt with different TOMs)
share the same indices of rate attribution.
"""
all_trt_smrs = [trt_smrs for sg in groups for src in sg
for trt_smrs in unc_subsets(src)]
unique, _ = get_unique_inverse(all_trt_smrs)
return [numpy.array(trt_smrs, numpy.uint32) for trt_smrs in unique]
[docs]def get_core_size(core_trt_smrs, full_lt):
"""
:param core_trt_smrs: the sets of realizations with the same
uncertainties, i.e. the units of rate attribution, as returned
by get_core_trt_smrs (already deduplicated)
:param full_lt: a FullLogicTree instance
:returns: the core size Gt = Σ_i G(trt_i) of the logic tree, i.e. the
number of distinct rate components, i.e. of gids (see
FullLogicTree.get_gids), which is >= len(core_trt_smrs) since
each unit contributes G(trt_i) >= 1 gids
"""
return sum(len(full_lt.gsim_lt.values[full_lt.trts[t[0] // TWO24]])
for t in core_trt_smrs)
[docs]def read_core_trt_smrs(dstore):
"""
:param dstore: a DataStore instance, possibly closed
:returns: a list of trt_smrs tuples, one per realization set, i.e.
len(dstore['core_trt_smrs']) tuples, less than the core size Gt
"""
with dstore: # NB: the datastore is closed when passed to a task
return [tuple(t) for t in dstore['core_trt_smrs'][:]]
[docs]def build_csm(oq, full_lt, smdict, dstore):
"""
:param oq: OqParam instance
:param full_lt: FullLogicTree instance
:param smdict: dictionary source_model_path -> SourceModel instance
:param dstore: DataStore instance
:returns: a CompositeSourceModel instance
"""
mon = performance.Monitor('_build_groups', measuremem=True)
with mon:
rlz_groups = []
for rlz in full_lt.sm_rlzs:
rlz_groups.extend((rlz, grp) for grp in
gen_groups(full_lt, smdict, rlz))
logging.info(mon)
logging.info('Building CompositeSourceModel')
groups = build_groups(full_lt, rlz_groups, oq, dstore)
csm = CompositeSourceModel(oq, full_lt, groups)
store_data(oq, smdict, csm, dstore)
return csm
[docs]def store_data(oq, smdict, csm, dstore):
"""
Create src_mutex, grp_probability in calc_XXX.hdf5 and sources
and mf_sections in calc_XXX_tmp.hdf5
"""
out = []
probs = []
for sg in csm.src_groups:
if sg.src_interdep == 'mutex' and 'src_mutex' not in dstore:
segments = []
for src in sg:
segments.append(src.source_id.split(':')[1])
t = (src.source_id, src.grp_id,
src.num_ruptures, src.mutex_weight,
sg.rup_interdep == 'mutex')
out.append(t)
probs.append((src.grp_id, sg.grp_probability))
assert len(segments) == len(set(segments)), segments
if out:
dtlist = [('src_id', hdf5.vstr), ('grp_id', int),
('num_ruptures', int), ('mutex_weight', float),
('rup_mutex', bool)]
dstore.create_dset('src_mutex', numpy.array(out, dtlist))
lst = [('grp_id', int), ('probability', float)]
dstore.create_dset('grp_probability', numpy.array(probs, lst))
# add rupids_by_tag to multifault sources if there is a single site
try:
sitecol = dstore['sitecol']
except (KeyError, TypeError): # 'NoneType' object is not subscriptable
sitecol = None
else:
# NB: in AELO we can have multiple vs30 on the same location
lonlats = set(zip(sitecol.lons, sitecol.lats))
if len(lonlats) > 1:
sitecol = None
# must be called *after* add_semicolons
t0 = time.time()
secparams = fix_geometry_sections(
smdict, csm.src_groups, dstore.tempname if dstore else '',
sitecol if oq.disagg_by_src and oq.use_rates else None,
split=not oq.calculation_mode.startswith('event_based'))
if secparams is not None and len(secparams):
logging.info('Spent %.1f seconds in fix_geometry_sections',
time.time()-t0)
# called by reduce_sources
[docs]def add_checksums(srcs):
"""
Build and attach a checksum to each source
"""
for src in srcs:
dic = {k: v for k, v in vars(src).items()
if k not in 'source_id sampling branch'}
src.checksum = zlib.adler32(pickle.dumps(dic, protocol=4))
# called before add_semicolons
[docs]def add_bangs(smdict):
"""
Discriminate different sources with same ID (false duplicates)
and put an exclamation mark in their source ID
"""
acc = general.AccumDict(accum=[])
atomic = set()
for smodel in smdict.values():
for sgroup in smodel.src_groups:
for src in sgroup:
src.branch = smodel.branch
srcid = (src.source_id if sgroup.atomic
else basename(src))
acc[srcid].append(src)
if sgroup.atomic:
atomic.add(src.source_id)
found = []
for srcid, srcs in acc.items():
if len(srcs) > 1: # duplicated ID
if any(src.source_id in atomic for src in srcs):
raise RuntimeError('Sources in atomic groups cannot be '
'duplicated: %s', srcid)
if any(getattr(src, 'mutex_weight', 0) for src in srcs):
raise RuntimeError('Mutually exclusive sources cannot be '
'duplicated: %s', srcid)
add_checksums(srcs)
gb = general.AccumDict(accum=[])
for src in srcs:
gb[checksum(src)].append(src)
if len(gb) > 1:
for same_checksum in gb.values():
for src in same_checksum:
check_branchID(src.branch)
src.source_id += '!%s' % src.branch
found.append(srcid)
return found
[docs]def fix_geometry_sections(smdict, src_groups, hdf5path='', site1=None,
split=True):
"""
If there are MultiFaultSources, fix the sections according to the
GeometryModels (if any).
"""
gmodels = []
gfiles = []
for fname, mod in smdict.items():
if isinstance(mod, nrml.GeometryModel):
gmodels.append(mod)
gfiles.append(fname)
# merge and reorder the sections
sec_ids = []
sections = {}
for gmod in gmodels:
sec_ids.extend(gmod.sections)
sections.update(gmod.sections)
check_unique(sec_ids, 'section ID in files ' + ' '.join(gfiles))
if sections:
# save in the temporary file sources and sections
assert hdf5path, ('You forgot to pass the dstore to '
'get_composite_source_model')
mfsources = []
for sg in src_groups:
for src in sg:
if src.code == b'F':
mfsources.append(src)
if mfsources:
split_dic, secparams = save_and_split(
mfsources, sections, hdf5path, site1, split=split)
for sg in src_groups:
new = []
for src in sg.sources:
tag = src.source_id
if tag in split_dic:
new.extend(split_dic[tag])
else:
new.append(src)
sg.sources[:] = new
return secparams
return None
def _groups_ids(smlt_dir, smdict, fnames):
# extract the source groups and ids from a sequence of source files
groups = []
for fname in fnames:
fullname = os.path.abspath(os.path.join(smlt_dir, fname))
groups.extend(smdict[fullname].src_groups)
return groups, set(src.source_id for grp in groups for src in grp)
def _add_sampling(src, rlz, trti):
# associate the source to the sampling parameters of the realization;
# the same source can appear in multiple realizations, hence the list.
# NB: the multiplicity of the source is len(src.sampling) and it enters
# the classical calculations too, via SourceGroup.fix_src_offset and
# the source_info rows, so the sampling must be always set
sampl = sampling(rlz.samples, trti * TWO24 + rlz.ordinal)
if src.sampling is None:
# the first time
src.sampling = [sampl]
else:
# if the same source belongs to multiple realizations
src.sampling.append(sampl)
[docs]def gen_groups(full_lt, smdict, rlz):
# yield all the possible source groups from the given rlz
smlt_file = full_lt.source_model_lt.filename
smlt_dir = os.path.dirname(smlt_file)
src_groups, source_ids = _groups_ids(
smlt_dir, smdict, rlz.value[0].split())
bset_values = full_lt.source_model_lt.bset_values(rlz.lt_path)
if rlz.ordinal % 100 == 0:
logging.info('Building source groups for rlz'
f'#{rlz.ordinal}: {"_".join(rlz.lt_path)}')
while (bset_values and
bset_values[0][0].uncertainty_type == 'extendModel'):
(_bset, value), *bset_values = bset_values
extra, extra_ids = _groups_ids(smlt_dir, smdict, value.split())
common = source_ids & extra_ids
if common:
raise InvalidFile(
'%s contains source(s) %s already present in %s' %
(value, common, rlz.value))
src_groups.extend(extra)
# NB: the uncertainties are not applied here, but in the workers,
# one set of realizations at a time, see modified_groups; the
# sampling info is set anyway, since it determines the multiplicity
for src_group in src_groups:
trti = full_lt.trti.get(src_group.trt, 0)
for src in src_group: # tested in case_83_eb
_add_sampling(src, rlz, trti)
yield src_group
# check applyToSources
sm_branch = rlz.lt_path[0]
src_id = full_lt.source_model_lt.info.applytosources[sm_branch]
for srcid in src_id:
if srcid not in source_ids:
if full_lt.source_model_lt.branchID:
continue
raise ValueError(
"The source %s is not in the source model,"
" please fix applyToSources in %s or the "
"source model(s) %s" % (srcid, smlt_file,
rlz.value[0].split()))
[docs]def reduce_sources(sources_with_same_id):
"""
:param sources_with_same_id: a list of sources with the same source_id
:returns: a list of truly unique sources
"""
# first reduce identical sources having the same id(src)
# tested in LogictreeTestCase.test_case_08, where <PoinstSource 2>
# appears 3 times
unique = get_unique(sources_with_same_id)
out = []
add_checksums(unique)
# in LogicTreeCase2ClassicalPSHA there 81 unique sources
# grouped in 9 groups of 9 sources each with the same checksum
for srcs in general.groupby(unique, checksum).values():
# NB: the simplest test featuring the same source in two
# different source models is logictree/case_01
src = srcs[0]
if len(srcs) > 1 and len(src.sampling) == 1:
src.sampling = numpy.concatenate([s.sampling for s in srcs])
out.append(src)
return out
[docs]def split_by_tom(sources):
"""
Groups together sources with the same TOM and collect multifault sources
"""
def key(src):
tom = getattr(src, 'temporal_occurrence_model', None)
return (tom.__class__.__name__, src.code == b'F')
return general.groupby(sources, key).values()
def _group_sources(trt, sources):
"""
Reduce identical sources, regroup by trt_smrs and TOM,
then return (source_groups, reduction_count).
"""
key = operator.attrgetter('source_id', 'code')
lst = []
red = 0
for srcs in general.groupby(sources, key).values():
if len(srcs) > 1:
srcs = reduce_sources(srcs)
red += 1
lst.extend(srcs)
src_groups = []
for sources in general.groupby(lst, trt_smrs).values():
for grp in split_by_tom(sources):
src_groups.append(sourceconverter.SourceGroup(trt, grp))
return src_groups, red
