Ctrl+K

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 configuration file and on the 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:

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:

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):

signatures of the area source first#

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

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:

# 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:

# 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 RateMaps 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 RateMaps 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.

On this page
Show Source