# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2018-2026 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 os
import copy
import operator
import collections
import numpy
from openquake.baselib import general, hdf5
from openquake.hazardlib.map_array import MapArray
from openquake.hazardlib.rmap_maker import get_unique_inverse
from openquake.hazardlib.calc.disagg import to_rates, to_probs
from openquake.hazardlib.source.rupture import BaseRupture
from openquake.hazardlib.source.rupture import get_ebrupture # noqa: F401
from openquake.commonlib.calc import get_proxies
U16 = numpy.uint16
U32 = numpy.uint32
I64 = numpy.int64
F32 = numpy.float32
F64 = numpy.float64
TWO24 = 2 ** 24
RATES_BLOCK = 2 ** 22 # rows of _rates read at once, to bound the memory
by_taxonomy = operator.attrgetter('taxonomy')
code2cls = BaseRupture.init()
weight = operator.itemgetter('n_occ')
slice_dt = numpy.dtype([('idx', U32), ('start', int), ('stop', int)])
[docs]class NotFound(Exception):
pass
[docs]def build_stat_curve(hcurve, imtls, stat, wget, use_rates=False):
"""
Build statistics by taking into account IMT-dependent weights
"""
weights = wget.weights
poes = hcurve.T # shape R, L
assert len(poes) == len(weights), (len(poes), len(weights))
L = imtls.size
array = numpy.zeros((L, 1))
if weights.shape[1] > 1: # IMT-dependent weights
# this is slower since the arrays are shorter
for imt in imtls:
slc = imtls(imt)
ws = wget(None, imt)
if not ws.sum(): # expect no data for this IMT
continue
if use_rates:
array[slc, 0] = to_probs(stat(to_rates(poes[:, slc]), ws))
else:
array[slc, 0] = stat(poes[:, slc], ws)
else:
if use_rates:
array[:, 0] = to_probs(stat(to_rates(poes), weights[:, -1]))
else:
array[:, 0] = stat(poes, weights[:, -1])
return array
[docs]def sig_eps_dt(imts):
"""
:returns: a composite data type for the sig_eps output
"""
lst = [('eid', U32), ('rlz_id', U16)]
for imt in imts:
lst.append(('sig_inter_' + imt, F32))
for imt in imts:
lst.append(('eps_inter_' + imt, F32))
return numpy.dtype(lst)
[docs]class HcurvesGetter(object):
"""
Read the contribution to the hazard curves coming from each source
in a calculation with a source specific logic tree
"""
def __init__(self, dstore):
self.dstore = dstore
self.imtls = dstore['oqparam'].imtls
self.full_lt = dstore['full_lt'].init()
self.sslt = self.full_lt.source_model_lt.decompose()
self.source_info = dstore['source_info'][:]
[docs] def get_hcurve(self, src_id, imt=None, site_id=0, gsim_idx=None):
"""
Return the curve associated to the given src_id, imt and gsim_idx
as an array of length L
"""
assert ';' in src_id, src_id # must be a realization specific src_id
imt_slc = self.imtls(imt) if imt else slice(None)
start, gsims, weights = self.bysrc[src_id]
dset = self.dstore['_rates']
if gsim_idx is None:
curves = dset[start:start + len(gsims), site_id, imt_slc]
return weights @ curves
return to_probs(dset[start + gsim_idx, site_id, imt_slc])
# NB: not used right now
[docs] def get_hcurves(self, src, imt=None, site_id=0, gsim_idx=None):
"""
Return the curves associated to the given src, imt and gsim_idx
as an array of shape (R, L)
"""
assert ';' not in src, src # not a rlz specific source ID
curves = []
for i in range(self.sslt[src].num_paths):
src_id = '%s;%d' % (src, i)
curves.append(self.get_hcurve(src_id, imt, site_id, gsim_idx))
return numpy.array(curves)
[docs] def get_mean_hcurve(self, src=None, imt=None, site_id=0, gsim_idx=None):
"""
Return the mean curve associated to the given src, imt and gsim_idx
as an array of shape L
"""
if src is None:
hcurves = [self.get_mean_hcurve(src) for src in self.sslt]
return general.agg_probs(*hcurves)
weights = [rlz.weight for rlz in self.sslt[src]]
curves = self.get_hcurves(src, imt, site_id, gsim_idx)
return weights @ curves
# NB: using 32 bit ratemaps
[docs]def get_rmap_gb(dstore, full_lt=None):
"""
:returns: (size_of_the_global_RateMap_in_GB, trt_rlzs, trt_smrs)
"""
N = len(dstore['sitecol/sids'])
L = dstore['oqparam'].imtls.size
full_lt = full_lt or dstore['full_lt'].init()
if 'core_trt_smrs' in dstore:
# for classical calculations the rates are attributed to the sets of
# realizations with the same uncertainties, not to the trt_smrs of
# the groups, since the groups are not split by the uncertainties,
# which are applied in the workers
key = 'core_trt_smrs'
else:
key = 'trt_smrs'
if key not in dstore: # starting from hazard_curves.csv
trt_smrs = [[0]]
else:
trt_smrs, _ = get_unique_inverse(dstore[key][:])
trt_rlzs = full_lt.get_trt_rlzs(trt_smrs)
max_gb = len(trt_rlzs) * N * L * 4 / 1024**3
return max_gb, trt_rlzs, trt_smrs
[docs]def get_num_chunks(dstore, full_lt=None):
"""
:returns: number of chunks to generate (determine postclassical tasks)
For performance, it is important to generate few chunks.
There are three regimes:
- classical_risk from hazard_curves.csv, num_chunks=1
- few sites, num_chunks=N
- regular, num_chunks=concurrent_tasks/2
- lots of data, num_chunks=req_gb
"""
oq = dstore['oqparam']
if 'hazard_curves' in oq.inputs:
return 1
N = len(dstore['sitecol/sids'])
ct2 = oq.concurrent_tasks // 2 or 1
if N < ct2 or oq.calculation_mode == 'disaggregation':
return N # one chunk per site
req_gb, _, _ = get_rmap_gb(dstore, full_lt)
ntiles = int(numpy.ceil(req_gb))
return ntiles if ntiles > ct2 else ct2
# for EUR on cole concurrent_tasks=256
# req_gb=202, N=260,000 => 202
[docs]def get_imt_gweights(full_lt, trt_smrs, imtls):
"""
:param full_lt: an initialized FullLogicTree
:param trt_smrs: a list of trt_smrs
:param imtls: a DictArray of IMTs
:returns: an array of shape (Gt, M) with the weights of each IMT,
the columns being in the same order as the IMTs in imtls
NB: the weights of the logic tree are stored in the order of the
IMTs appearing in the branches, which is not necessarily the order
of the IMTs in the calculation, so the mapping must be done with
the IMTWeigher (see IMTWeigher.__call__)
"""
gweights = full_lt.g_weights(trt_smrs) # (Gt, M+1)
wget = full_lt.gsim_lt.wget
return numpy.array([wget(gweights, imt) for imt in imtls]).T
[docs]def map_getters(dstore, full_lt=None, oq=None, disagg=False):
"""
:returns: a list of pairs (MapGetter, weights)
"""
oq = oq or dstore['oqparam']
n = get_num_chunks(dstore, full_lt)
# full_lt is None in classical_risk, classical_damage
full_lt = full_lt or dstore['full_lt']
if getattr(full_lt, 'weights', None) is None: # not initialized
full_lt.init()
R = full_lt.get_num_paths()
_req_gb, trt_rlzs, trt_smrs = get_rmap_gb(dstore, full_lt)
attrs = vars(full_lt)
if oq.fastmean:
gweights = [get_imt_gweights(full_lt, trt_smrs, oq.imtls)]
else:
wgets = [full_lt.gsim_lt.wget]
for label in oq.site_labels:
flt = copy.copy(full_lt)
flt.__dict__.update(attrs)
flt.gsim_lt = dstore['gsim_lt' + label]
flt.init()
if oq.fastmean:
gweights.append(get_imt_gweights(flt, trt_smrs, oq.imtls))
else:
wgets.append(flt.gsim_lt.wget)
fnames = [dstore.filename]
calc_dir = dstore.filename[:-5]
if os.path.exists(calc_dir):
for f in os.listdir(calc_dir):
if f.endswith('.hdf5'):
fnames.append(os.path.join(calc_dir, f))
out = []
sids = dstore['sitecol/sids'][:]
for chunk in range(n):
tile = sids[sids % n == chunk]
getter = MapGetter(fnames, chunk, trt_rlzs, tile, R, oq)
if oq.fastmean:
getter.gweights = gweights
else:
getter.wgets = wgets
if oq.site_labels:
getter.ilabels = dstore['sitecol'].ilabel
out.append(getter)
return out
[docs]class ZeroGetter(object):
"""
Return an array of zeros of shape (L, R)
"""
def __init__(self, L, R):
self.L = L
self.R = R
[docs] def get_hazard(self):
return numpy.zeros((self.L, self.R))
[docs]class CurveGetter(object):
"""
Hazard curve builder used in classical_risk/classical_damage.
:param sid: site index
:param rates: array of shape (L, G) for the given site
"""
[docs] @classmethod
def build(cls, dstore, full_lt):
"""
:returns: a dictionary sid -> CurveGetter
"""
rates = {}
for mgetter in map_getters(dstore, full_lt):
array = mgetter.init()
for sid, idx in mgetter.sid2idx.items():
rates[sid] = array[idx] # shape (L, G)
dic = collections.defaultdict(lambda: ZeroGetter(mgetter.L, mgetter.R))
for sid in rates:
dic[sid] = cls(sid, rates[sid], mgetter.trt_rlzs, mgetter.R)
return dic
def __init__(self, sid, rates, trt_rlzs, R):
self.sid = sid
self.rates = rates
self.trt_rlzs = trt_rlzs
self.R = R
[docs] def get_hazard(self):
r0 = numpy.zeros((len(self.rates), self.R))
for g, t_rlzs in enumerate(self.trt_rlzs):
rlzs = t_rlzs % TWO24
rates = self.rates[:, g]
for rlz in rlzs:
r0[:, rlz] += rates
return to_probs(r0)
[docs]class DeltaRatesGetter(object):
"""
Read the delta rates from an aftershock datastore
"""
def __init__(self, dstore):
self.dstore = dstore
def __call__(self, src_id):
with self.dstore: # opened and closed here, if not already open
return self.dstore['delta_rates'][src_id]
[docs]class MapGetter(object):
"""
Read hazard curves from the datastore for all realizations or for a
specific realization.
"""
def __init__(self, filenames, chunk, trt_rlzs, sids, R, oq):
self.filenames = filenames
self.chunk = chunk
self.trt_rlzs = trt_rlzs
self.sids = sids
self.R = R
self.imtls = oq.imtls
self.poes = oq.poes
self.use_rates = oq.use_rates
self.eids = None
self.ilabels = () # overridden in case of ilabels
self.array = None
@property
def imts(self):
return list(self.imtls)
@property
def Gt(self):
return len(self.trt_rlzs)
@property
def L(self):
return self.imtls.size
@property
def N(self):
return len(self.sids)
@property
def M(self):
return len(self.imtls)
[docs] def init(self):
"""
Build the array from the underlying dataframes
"""
if self.array is not None:
return self.array
sid2idx = {sid: idx for idx, sid in enumerate(self.sids)}
self.array = numpy.zeros((self.N, self.L, self.Gt)) # move to 32 bit
lookup = numpy.full(max(sid2idx) + 1, -1, int)
for sid, idx in sid2idx.items():
lookup[sid] = idx
for fname in self.filenames:
self.read_chunk(fname, lookup)
self.sid2idx = sid2idx
return self.array
[docs] def read_chunk(self, fname, lookup):
"""
Add the rates of the chunk stored in the file fname to self.array
:param fname: the path to a datastore with a _rates group
:param lookup: an array mapping site IDs to indices of self.array
"""
with hdf5.File(fname) as dstore:
slices = dstore['_rates/slice_by_idx'][:]
slices = slices[slices['idx'] == self.chunk]
if len(slices) == 0:
return
# read each column block by block (reading slice by slice
# is slow) and select the rows of the slices of this chunk
n = len(dstore['_rates/rate'])
for a in range(0, n, RATES_BLOCK):
b = min(a + RATES_BLOCK, n)
# clip the slices to the block and mark their rows
start = numpy.clip(slices['start'], a, b) - a
stop = numpy.clip(slices['stop'], a, b) - a
marks = numpy.zeros(b - a + 1, int)
numpy.add.at(marks, start, 1)
numpy.add.at(marks, stop, -1)
mask = numpy.cumsum(marks[:-1]) > 0
if not mask.any():
continue
sids = dstore['_rates/sid'][a:b][mask]
lid = dstore['_rates/lid'][a:b][mask]
gid = dstore['_rates/gid'][a:b][mask]
rate = dstore['_rates/rate'][a:b][mask]
idxs = lookup[sids]
if (idxs < 0).any():
raise KeyError('unknown site IDs in _rates')
# NB: numpy.add.at is needed since the same
# (sid, lid, gid) can appear several times, and
# fancy-index += would drop duplicates
numpy.add.at(self.array, (idxs, lid, gid), rate)
[docs] def get_hcurve(self, sid): # used in classical
"""
:param sid: a site ID
:returns: an array of shape (L, R) for the given site ID
"""
array = self.init()
r0 = numpy.zeros((self.L, self.R))
idx = self.sid2idx[sid]
for g, t_rlzs in enumerate(self.trt_rlzs):
rlzs = t_rlzs % TWO24
rates = array[idx, :, g]
for rlz in rlzs:
r0[:, rlz] += rates
return to_probs(r0)
[docs] def get_fast_mean(self):
"""
:returns: a MapArray of shape (N, M, L1) with the mean hcurves
"""
M = self.M
L1 = self.L // M
means = MapArray(U32(self.sids), M, L1).fill(0)
for sid in self.sids:
gweights = self.gweights[0] # (Gt, M), one column per IMT
if len(self.ilabels):
gweights = self.gweights[self.ilabels[sid]]
rates = self.array[self.sid2idx[sid]] # shape (L, G)
sidx = means.sidx[sid]
for m in range(M):
means.array[sidx, m] = rates[m*L1: m*L1+L1] @ gweights[:, m]
# NB: the conversion must be done in float64, otherwise the
# cancellation in 1 - exp(-rate) is catastrophic for small rates
means.array[:] = to_probs(means.array.astype(F64))
return means
[docs]def get_ebruptures(dstore):
"""
Extract EBRuptures from the datastore
"""
ebrs = []
trts = list(dstore['full_lt/gsim_lt'].values)
for proxy in get_proxies(dstore.filename):
ebrs.append(proxy.to_ebr(trts[0]))
return ebrs
[docs]def line(points):
return '(%s)' % ', '.join('%.5f %.5f %.5f' % tuple(p) for p in points)
[docs]def multiline(array3RC):
"""
:param array3RC: array of shape (3, R, C)
:returns: a MULTILINESTRING
"""
D, R, _C = array3RC.shape
assert D == 3, D
lines = 'MULTILINESTRING(%s)' % ', '.join(
line(array3RC[:, r, :].T) for r in range(R))
return lines
