.. _classical-internals: How the classical calculator really works ------------------------------------------- This page is for GEM personnel and for users willing to dig into the implementation of the engine. It is *not* needed to run the engine: if you only want to compute hazard curves, read the sections on the :ref:`configuration file ` and on the :ref:`outputs ` instead. Here we explain how a classical calculation is actually organized, i.e. how the sources and the logic tree are turned into a CompositeSourceModel, how the epistemic uncertainties are managed and how the hazard rates are accumulated. Throughout the page we use the demo ``demos/hazard/LogicTreeCase2ClassicalPSHA``, which is small enough to be run in a few seconds and rich enough to contain epistemic uncertainties:: $ cd demos/hazard/LogicTreeCase2ClassicalPSHA $ oq engine --run job.ini The demo has a single site, 19 intensity measure levels, an investigation time of 50 years, ``disagg_by_src = true`` and a full enumeration of the logic tree (``number_of_logic_tree_samples = 0``). The source model contains two source groups, each with a single source: an area source with id ``first`` in the Stable Continental Crust and a simple fault source with id ``second`` in the Active Shallow Crust. The source model logic tree has a ``sourceModel`` branchset and four uncertainty branchsets, ``abGRAbsolute`` and ``maxMagGRAbsolute`` applied to the first and to the second source, each with three branches: .. code-block:: text sourceModel(1) x abGRAbsolute first(3) x maxMagGRAbsolute first(3) x abGRAbsolute second(3) x maxMagGRAbsolute second(3) = 81 source model paths The ground motion logic tree has two GMMs per tectonic region type, therefore the calculation has 81 x 4 = 324 hazard realizations. The calculation pipeline ~~~~~~~~~~~~~~~~~~~~~~~~ A classical calculation is performed in two steps, the *preclassical* and the *classical*. Before the preclassical there is the reading phase (``HazardCalculator.read_inputs``): the job parameters and the logic trees are read by the ``logictree`` module and the source models are read by ``source_reader.read_source_model`` in parallel, one task per source model file. Then the CompositeSourceModel is built (without applying the uncertainties, see below), stored in memory and saved in the datastore. The preclassical proper, i.e. ``PreClassicalCalculator.populate_csm``, then works on the sources, group by group: - it builds the ``cmakers``, i.e. one ``ContextMaker`` per source group, and saves the ``trt_smrs`` of the groups in the ``trt_smrs`` dataset and the *indices of rate attribution* in the ``core_trt_smrs`` dataset, i.e. the sets of realizations the rates are attributed to (see the section on the uncertainties) - it pre-filters the sources against the sites, using a coarse site grid when there are many sites (with ``res = 4``, i.e. 40km), since the pre-filters are only used to estimate the cost of each source - it splits the area and multi-point sources into point sources, unless they are modified by the uncertainties (see below) - for the multi-fault sources it computes the multi-surface parameters of each rupture, i.e. the area, dip, strike, width, ztor, zbot, rupture width and bounding box of the multisurface, as area-weighted averages over the sections composing it (``set_msparams``) - it filters the sources by magnitude and counts their ruptures - it estimates the weight of each source, i.e. its CPU cost, from the number of contexts generated by the source (``estimate_weight``); the sources modified by the uncertainties are not split, so their cost is simply ``num_ruptures * nsites * nsubsets`` Finally the groups are split in blocks and site tiles, stored in the ``_csm`` group of the datastore and described by the ``source_groups`` dataset. The classical reads the sources back from ``_csm`` and computes the hazard rates. Both steps send tasks to the workers; the tasks of a calculation are listed by the command:: $ oq show task_info For the demo it returns:: | operation-duration | counts | mean | stddev | min | max | slowfac | |--------------------+--------+-----------+--------+-----------+-----------+---------| | classical_bysrc | 2 | 1.7653 | 39% | 1.0849 | 2.4456 | 1.3854 | | filter_weight | 2 | 0.0024 | 04% | 0.0023 | 0.0025 | 1.0386 | | postclassical | 1 | 0.0274 | nan | 0.0274 | 0.0274 | 1.0000 | | preclassical | 2 | 1.035E-04 | 24% | 7.820E-05 | 1.287E-04 | 1.2442 | | read_source_model | 1 | 0.0012 | nan | 0.0012 | 0.0012 | 1.0000 | ``read_source_model`` is a task of the reading phase, the ``preclassical`` tasks split the sources and generate the ``filter_weight`` subtasks, then the classical sends two ``classical_bysrc`` tasks (one per source group) and finally ``postclassical`` builds the hazard curves and maps. The last column, ``slowfac``, is the ratio between the maximum and the mean duration per worker: a large value means that the work was not balanced evenly across the workers, i.e. that there were slow tasks. The CompositeSourceModel ~~~~~~~~~~~~~~~~~~~~~~~~ .. note:: The structure of the CompositeSourceModel changed completely in version 3.27: before, the uncertainties were applied when reading the source models, i.e. there was a copy of the sources for each set of uncertainties and consequently a group (and a ``ContextMaker``) per set of uncertainties; now the uncertainties are applied in the workers and the structure of the CSM depends on the source models only, as described below. The CompositeSourceModel (CSM) is the in-memory representation of the source models for a given logic tree. It is a list of ``SourceGroup`` instances, each containing a list of sources; the groups are the unit of task generation, the sources the unit of computation. Each source has a ``sampling`` attribute, i.e. an array with two columns, ``trt_smr`` and ``samples``, with one row per source model realization the source belongs to:: sampling_dt = numpy.dtype([('trt_smr', U32), ('samples', U32)]) The ``trt_smr`` is an integer encoding two pieces of information, the tectonic region type and the realization: .. code-block:: text trt_smr = trti * 2**24 + ordinal where ``trti`` is the index of the tectonic region type in the GMM logic tree (``full_lt.trti``) and ``ordinal`` is the ordinal of the *effective* source model realization, i.e. of the row of the ``sm_rlzs`` table:: $ oq show sm_rlzs | ordinal | lt_path | value | samples | weight | |---------+---------------------+-----------------------------------------------------------+---------+--------| | 0 | b11_b21_b31_b41_b51 | source_model.xml 4.6000 1.1000 3.3000 1.0000 7.0000 7.5000 | 4 | 0.0123 | | 1 | b11_b21_b31_b41_b52 | source_model.xml 4.6000 1.1000 3.3000 1.0000 7.0000 7.8000 | 4 | 0.0123 | ... In the demo there are 81 effective source model realizations, so the area source ``first`` has the ``trt_smrs`` ``16777216, ..., 16777296`` and the fault source ``second`` has the ``trt_smrs`` ``0, 1, ..., 80``: the ``trt_smrs`` of a source are offset by ``2**24`` for each tectonic region type after the first, and the ordinals are the same for all the sources. Each of them has ``samples = 4``, the four GMM paths of the GMM logic tree, so that 81 x 4 = 324 are the hazard realizations of the calculation. The CSM is built by ``source_reader.build_groups`` and it has two properties worth remembering: - there is one group per source group of the source model files (then regrouped by ``trt_smrs`` and temporal occurrence model), i.e. the structure of the CSM, the number of groups and the number of ``cmakers`` depend on the source models only and *not* on the uncertainties: for the demo there are two groups, one per tectonic region type, and two ``cmakers``, each with the two GMMs of its tectonic region type - a source appears **once**, with the ``trt_smrs`` of all the realizations it belongs to, and *not* once per realization (see the next section) The groups are stored in the ``_csm`` group of the datastore and their parameters in the ``source_groups`` dataset, one row per group:: $ oq show source_groups | grp_id | gsims | tiles | blocks | max_mb | weight | codes | trt | |--------+-------+-------+--------+-----------+--------+-------+--------------------------| | 0 | 2 | 1 | 1 | 1.450E-04 | 3_620 | S | Active Shallow Crust | | 1 | 2 | 1 | 1 | 1.450E-04 | 28_236 | A | Stable Continental Crust | ``gsims`` is the number of GMMs per group, ``tiles`` the number of site tiles, ``blocks`` the number of blocks the sources are split in, ``weight`` the estimated cost of the group and ``codes`` the source codes (``S`` simple fault, ``A`` area source). .. note:: The datastore contains two datasets with the realizations of the groups, ``trt_smrs`` and ``core_trt_smrs``, which look similar but have a different meaning. ``trt_smrs`` has one row per source group, with the ``trt_smrs`` of all the realizations of the group, and it is used to build the cmakers and the tasks; ``core_trt_smrs`` has one row per *realization set*, i.e. per subset of realizations with the same uncertainties, so it has more rows and they are shorter:: >> from openquake.commonlib import datastore >> ds = datastore.read(calc_id) >> [len(t) for t in ds['trt_smrs'][:]] [81, 81] >> [len(t) for t in ds['core_trt_smrs'][:]] [9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9] >> ds['core_trt_smrs'][:2] [array([ 0, 3, 6, 27, 30, 33, 54, 57, 60], dtype=uint32) array([ 1, 4, 7, 28, 31, 34, 55, 58, 61], dtype=uint32)] Every realization set is contained in exactly one row of ``trt_smrs``, but the converse does not hold: the trt_smrs of a group can be spread over several realization sets, since the sources of a group can be affected by different uncertainties (see the page on correlated uncertainties), so a realization can belong to more than one index. The rates cannot be attributed to a whole group, since the uncertainties change within it, hence the attribution happens at the level of the indices; see the section on the uncertainties below. Management of the uncertainties ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The engine applies the epistemic uncertainties in the classical phase, one set of realizations at a time. The mechanism is the following. 1. While reading the logic tree, for each pair (realization, source) the engine computes the *signature* of the uncertainties applying to that source, i.e. the tuple of ``(branchset id, value)`` pairs (see ``source_reader.unc_signature``). Sources not affected by a branchset have an empty signature. 2. When building the CSM, the ``trt_smrs`` of a source are grouped by signature and stored in the ``bysrc_subsets`` attribute; the sets of realizations sharing the same uncertainties are called *subsets*. They are read back with ``lt.unc_subsets(src)``. 3. In the preclassical the subsets are stored in the ``core_trt_smrs`` dataset, one row per subset: these are the *indices of rate attribution*, i.e. the indices the rates are computed and attributed with, as opposed to the ``trt_smrs`` of a group, which contain all the realizations. Concretely, a subset is turned into a set of ``gids``, i.e. the (realization, GMM) pairs of the subset, and the rates computed for it are associated to those ``gids``; see the section on the ``RateMap`` below. 4. In the workers, for each subset the engine restricts the sampling of the sources to the realizations of the subset (``lt.restrict_sampling``) and applies the uncertainties of the first realization of the subset (``lt.apply_uncertainties``, called by ``source_reader.modified_groups``). The resulting sources are then split (if needed), filtered by magnitude and used to compute the rates, which are attributed to the ``gids`` of the subset. For the demo each source has 81 ``trt_smrs`` organized in **9 subsets of 9 realizations**. In the logic tree the branchset ``bs1`` is the ``sourceModel``, while ``bs2`` and ``bs3`` are the ``abGRAbsolute`` branchsets applied to ``first`` and to ``second`` and ``bs4`` and ``bs5`` are the ``maxMagGRAbsolute`` ones. The table below lists the signatures of the area source ``first``, i.e. the branchsets and values applying to it, together with the ordinals of the realizations having that signature (the corresponding ``trt_smrs`` are the ordinals plus ``2**24``, since ``first`` is in the second tectonic region type): .. list-table:: signatures of the area source ``first`` :header-rows: 1 :widths: 30 70 * - signature - ordinals of the realizations * - ``bs2=(4.6, 1.1)``, ``bs4=7.0`` - 0, 1, 2, 9, 10, 11, 18, 19, 20 * - ``bs2=(4.6, 1.1)``, ``bs4=7.3`` - 3, 4, 5, 12, 13, 14, 21, 22, 23 * - ``bs2=(4.6, 1.1)``, ``bs4=7.6`` - 6, 7, 8, 15, 16, 17, 24, 25, 26 * - ``bs2=(4.5, 1.0)``, ``bs4=7.0`` - 27, 28, 29, 36, 37, 38, 45, 46, 47 * - ``bs2=(4.5, 1.0)``, ``bs4=7.3`` - 30, 31, 32, 39, 40, 41, 48, 49, 50 * - ``bs2=(4.5, 1.0)``, ``bs4=7.6`` - 33, 34, 35, 42, 43, 44, 51, 52, 53 * - ``bs2=(4.4, 0.9)``, ``bs4=7.0`` - 54, 55, 56, 63, 64, 65, 72, 73, 74 * - ``bs2=(4.4, 0.9)``, ``bs4=7.3`` - 57, 58, 59, 66, 67, 68, 75, 76, 77 * - ``bs2=(4.4, 0.9)``, ``bs4=7.6`` - 60, 61, 62, 69, 70, 71, 78, 79, 80 Only the branchsets ``bs2`` and ``bs4`` apply to ``first``, so its signature has 3 x 3 = 9 values, and each of them covers the 9 realizations obtained by varying the uncertainties of ``second``, which do not apply to ``first``: the ordinals of each row are the 3x3 combinations of the branches of ``bs3`` and ``bs5``. The fault source ``second`` has the analogous 9 signatures, built from ``bs3`` and ``bs5``, therefore there are 18 indices of rate attribution, 9 per source, stored in the ``core_trt_smrs`` dataset. The advantage is visible in the size of the CSM and in the number of tasks: if the uncertainties were applied when reading the source models, the demo would have 9 copies of each source and 18 groups (with 18 ``cmakers``) instead of 2, i.e. 9 times the data to transfer to the workers and 9 times the tasks, for the same hazard curves. It is also what makes the task structure independent from the logic tree: adding realizations to the logic tree does not change the number of groups. There are two consequences in the preclassical worth knowing about, both controlled by the ``bysrc_unc`` flag set on the sources modified by the uncertainties: - such sources are **not split**, since splitting an area source would destroy the geometry the uncertainties apply to (and for fault sources it would require recomputing the rupture counts); they are split in the workers after the uncertainties have been applied (``preclassical.split_modified``) - they are **not filtered by magnitude**, since the filtering depends on the occurrence rates, which the uncertainties modify; the filtering is performed in the workers (``preclassical.filter_mag``) together with the check on the maximum number of ruptures .. _ratemap: The RateMap ~~~~~~~~~~~ The classical calculator computes *rates* (annual exceedance rates) and not probabilities of exceedance, so that the contributions of the sources can be summed linearly; the conversion to PoEs happens at the end, in ``build_curves_maps``. A ``RateMap`` is a wrapper around an array of shape ``(N, L, G)`` of float32 rates, plus a dictionary mapping the ``gids`` to the columns of the array: .. code-block:: python # openquake/hazardlib/map_array.py class RateMap: def __init__(self, sids, L, gids): self.shape = len(sids), L, len(gids) self.array = numpy.zeros(self.shape, F32) self.gdic = {g: j for j, g in enumerate(gids)} Each task returns the ``RateMap`` of one subset of realizations, i.e. with ``G = 2`` columns for the demo (the two GMMs of the tectonic region type), and the master accumulates it into the ``RateMap`` of the group: .. code-block:: python # openquake/hazardlib/map_array.py def __iadd__(self, other): sidx = self.sidx[other.sids] for i, g in enumerate(other.gid): oarray = other.array[:, :, i] self.array[sidx, :, self.gdic[g]] += oarray return self The group ``RateMap``\ s are allocated once in the master, with the ``gids`` of all the subsets of the group, so that the accumulation is a ``(N, L, G)`` array addition. A ``RateMap`` is a potentially large object and the engine checks upfront that there is enough memory, logging a line like:: Requiring 12.3 GB for the RateMaps and refusing to start if the required memory is not available. The check is not a global one, since a ``RateMap`` is allocated in the master only in three cases: when there are few sites, i.e. no more than ``max_sites_disagg`` (10 by default), when ``disagg_by_src`` is set and when the sources of a group are split in blocks. The estimate in ``preclassical.get_req_gb`` is ``N * L * gsims * 4`` bytes and it is computed for the groups falling in the last two cases only, i.e. exactly the ones with many sites. With many sites, no ``disagg_by_src`` and a single block per group there is nothing to accumulate in the master, so the tasks return arrays of rates which are stored right away (``baseclassical`` with ``as_rmap = false``). In all cases the rates end up in the ``_rates`` dataset as a table of ``(sid, lid, gid, rate)`` rows, using the *global* ``gids``, and are then converted into hazard curves:: >> from openquake.commonlib import datastore >> ds = datastore.read(calc_id) >> df = ds.read_df('_rates') >> len(df) # 579 rows for 36 gids and 19 IMLs in the demo Two final remarks about the rates: the sites with zero rate are removed in the worker before the rates are returned (``remove_zeros``), which reduces the data sent back to the master and the size of the arrays stored; and the rates of a group are also used to compute the mean rates by source (``oq show mean_rates_by_src``) and, with ``disagg_by_src``, the disaggregations are attributed to the individual sources via the ``basename`` key of the ``RateMap``. From the rates to the hazard curves ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The rates are stored in the ``_rates`` dataset as a table of ``(sid, lid, gid, rate)`` rows, using the *global* gids:: >> from openquake.commonlib import datastore >> ds = datastore.read(calc_id) >> df = ds.read_df('_rates') >> len(df) # 579 rows for 36 gids and 19 IMLs in the demo A gid is just a global index, i.e. the number of a column of the rates: ``FullLogicTree.get_gids`` allocates one consecutive gid per index of rate attribution and per GMM, so there are 36 ``gids`` for the 18 indices of the demo, and a gid carries no information by itself. The realizations behind a gid are given by ``FullLogicTree.get_trt_rlzs``, which associates to each gid an array of ``rlz + 2**24 * trti``, i.e. the realizations having that GMM in that tectonic region type: in the demo each gid is shared by 18 realizations and each realization has 2 ``gids``, since the GMM logic tree of the demo has a branching level per tectonic region type. The last step of the calculation is performed by the objects returned by ``getters.map_getters(dstore, full_lt, oq)``, i.e. one ``MapGetter`` per *chunk* of sites: ``getters.get_num_chunks`` decides how many chunks to generate, one per site if there are few sites, ``concurrent_tasks / 2`` otherwise, and more if the rates are large. Each getter knows - ``trt_rlzs``, i.e. for each gid the list of realizations having it, from ``FullLogicTree.get_trt_rlzs`` - its own sites, the ones with ``sid % nchunks == chunk`` ``MapGetter.init`` reads the ``_rates`` dataframes, possibly from more than one file (the extra files are written by ``save_rates`` when the ``RateMap``\ s are too large to be kept in memory), and rebuilds a dense array of shape ``(N, L, Gt)`` of rates. Then ``get_hcurve`` reconstructs the hazard curve of a site, of shape ``(L, R)``:: # openquake/calculators/getters.py def get_hcurve(self, sid): array = self.init() # (N, L, Gt) rates r0 = numpy.zeros((self.L, self.R)) idx = self.sid2idx[sid] for g, t_rlzs in enumerate(self.trt_rlzs): rlzs = t_rlzs % TWO24 # the realizations behind the gid rates = array[idx, :, g] for rlz in rlzs: r0[:, rlz] += rates # the same rates for all of them return to_probs(r0) # rates -> PoEs i.e. the rates of a gid are added to every realization sharing it and the result is converted into probabilities of exceedance with ``to_probs``, i.e. ``1 - exp(-rate * investigation_time)``. The task ``postclassical`` loops over the sites of its chunk, calls ``get_hcurve`` and computes the statistics (mean and quantiles, weighted by the weights of the GMM logic tree) with ``getters.build_stat_curve``; with ``fastmean = true`` the weighted mean is computed directly from the rates by ``MapGetter.get_fast_mean``, skipping the per-realization curves. The resulting MapArrays are stored in the ``hcurves-rlzs`` and ``hcurves-stats`` datasets, and the maps in ``hmaps-rlzs`` and ``hmaps-stats``. .. note:: The logic tree of the demo has 324 realizations but only 36 gids: 36 is the *core size* of the logic tree, i.e. the number of distinct rate components, each one standing for a set of realizations with the same uncertainties and the same GMM, as shown by the rows of ``oq show core_trt_smrs`` (18 realization sets of 9 realizations, one per GMM, since the demo has two GMMs per tectonic region type). In general the core size is much smaller than the size of the logic tree, since with a sampled logic tree many realizations share the same set of uncertainties; in the extreme case of a single GMM and no uncertainties the two are equal. The rates of all the realizations can be reconstructed from the core ``RateMap``, since each realization is the sum of the columns of its gids; that is exactly what ``MapGetter.init`` and ``get_hcurve`` do:: >> from openquake.commonlib import datastore >> from openquake.calculators import getters >> ds = datastore.read(calc_id) >> full_lt = ds['full_lt'].init() >> oq = ds['oqparam'] >> full_lt.get_num_paths(), getters.get_rmap_gb(ds, full_lt)[1] \ ... # (324, [36 arrays of realizations]) (324, [array([0, 1, 12, ...], dtype=uint32), ...]) >> getter = getters.map_getters(ds, full_lt, oq)[0] >> rates = getter.init() # (N, L, Gt) core RateMap of rates >>> rates.shape (1, 19, 36) >> hcurves = getter.get_hcurve(getter.sids[0]) # (L, R) >>> hcurves.shape (19, 324) i.e. the core ``RateMap`` of shape ``(N, L, Gt)`` contains all the rates of the ``Gt`` components, and the ``R`` columns of the hazard curves are obtained by summing the ``G(rlz)`` columns belonging to each realization.