.. _advanced-calculations:

Advanced Calculations
==========================


Large Calculations in Classical PSHA
------------------------------------

Running large hazard calculations, especially ones with large logic trees, is an art, and there are various techniques 
that can be used to reduce an impossible calculation to a feasible one.

**********************
Reducing a calculation
**********************

The first thing to do when you have a large calculation is to reduce it so that it can run in a reasonable amount of 
time. For instance you could reduce the number of sites, by considering a small geographic portion of the region 
interested, of by increasing the grid spacing. Once the calculation has been reduced, you can run it and determine what 
are the factors dominating the run time.

As we discussed in section about common mistakes, you may want to tweak the quadratic parameters (``maximum_distance``, 
``area_source_discretization``, ``rupture_mesh_spacing``, ``complex_fault_mesh_spacing``). Also, you may want to choose 
different GMPEs, since some are faster than others. You may want to play with the logic tree, to reduce the number of 
realizations: this is especially important, in particular for event based calculation were the number of generated 
ground motion fields is linear with the number of realizations.

Once you have tuned the reduced computation, you can have an idea of the time required for the full calculation. It 
will be less than linear with the number of sites, so if you reduced your sites by a factor of 100, the full computation 
will take a lot less than 100 times the time of the reduced calculation (fortunately). Still, the full calculation can 
be impossible because of the memory/data transfer requirements, especially in the case of event based calculations. 
Sometimes it is necessary to reduce your expectations. The examples below will discuss a few concrete cases. But first 
of all, we must stress an important point::

	Our experience tells us that THE PERFORMANCE BOTTLENECKS OF THE
	REDUCED CALCULATION ARE TOTALLY DIFFERENT FROM THE BOTTLENECKS OF
	THE FULL CALCULATION. Do not trust your performance intuition.

*********************************
Classical PSHA for Europe (SHARE)
*********************************

Suppose you want to run a classical PSHA calculation for the latest model for Europe and that it turns out to be too 
slow to run on your infrastructure. Let’s say it takes 4 days to run. How do you proceed to reduce the computation time?

The first thing that comes to mind is to tune the ``area_source_discretization`` parameter, since the calculation (as 
most calculations) is dominated by area sources. For instance, by doubling it (say from 10 km to 20 km) we would expect 
to reduce the calculation time from 4 days to just 1 day, a definite improvement.

But how do we check if the results are still acceptable? Also, how we check that in less than 4+1=5 days? As we said 
before we have to reduce the calculation and the engine provides several ways to do that.

If you want to reduce the number of sites, IMTs and realizations you can:

- manually change the ``sites.csv`` or ``site_model.csv`` files
- manually change the ``region_grid_spacing``
- manually change the ``intensity_measure_types_and_levels`` variable
- manually change the GMPE logic tree file by commenting out branches
- manually change the source logic tree file by commenting out branches
- use the environment variable ``OQ_SAMPLE_SITES``
- use the environment variable ``OQ_REDUCE``

Starting from engine 3.11 the simplest approach is to use the ``OQ_REDUCE`` environment variable than not only reduce 
reduces the number of sites, but also reduces the number of intensity measure types (it takes the first one only) and 
the number of realizations to just 1 (it sets ``number_of_logic_tree_samples=1``) and if you are in an event based 
calculation reduces the parameter ``ses_per_logic_tree_path`` too. For instance the command::

	$ OQ_REDUCE=.01 oq engine --run job.ini

will reduce the number of sites by 100 times by random sampling, as well a reducing to 1 the number of IMTs and 
realizations. As a result the calculation will be very fast (say 1 hour instead of 4 days) and it will possible to 
re-run it multiple times with different parameters. For instance, you can test the impact of the area source 
discretization parameter by running::

	$ OQ_REDUCE=.01 oq engine --run job.ini --param area_source_discretization=20

Then the engine provides a command oq compare to compare calculations; for instance::

	$ oq compare hmaps PGA -2 -1 --atol .01

will compare the hazard maps for PGA for the original (ID=-2, area_source_discretization=10 km) and the new calculation 
(ID=-2, area_source_discretization=20 km) on all sites, printing out the sites where the hazard values are different 
more than .01 g (``--atol`` means absolute tolerence). You can use ``oq compare --help`` to see what other options are 
available.

If the call to ``oq compare`` gives a result::

	There are no differences within the tolerances atol=0.01, rtol=0%, sids=[...]

it means that within the specified tolerance the hazard is the same on all the sites, so you can safely use the area 
discretization of 20 km. Of course, the complete calculation will contain 100 times more sites, so it could be that in 
the complete calculation some sites will have different hazard. But that’s life. If you want absolute certitude you will 
need to run the full calculation and to wait. Still, the reduced calculation is useful, because if you see that are 
already big differences there, you can immediately assess that doubling the ``area_source_discretization`` parameter is a 
no go and you can try other strategies, like for instance doubling the ``width_of_mfd_bin`` parameter.

As of version 3.11, the ``oq compare hmaps`` command will give an output like the following, in case of differences::

	site_id calc_id 0.5     0.1     0.05    0.02    0.01    0.005
	======= ======= ======= ======= ======= ======= ======= =======
	767     -2      0.10593 0.28307 0.37808 0.51918 0.63259 0.76299
	767     -1      0.10390 0.27636 0.36955 0.50503 0.61676 0.74079
	======= ======= ======= ======= ======= ======= ======= =======
	===== =========
	poe   rms-diff
	===== =========
	0.5   1.871E-04
	0.1   4.253E-04
	0.05  5.307E-04
	0.02  7.410E-04
	0.01  8.856E-04
	0.005 0.00106
	===== =========

This is an example with 6 hazard maps, for poe = .5, .1, .05, .02, .01 and .005 respectively. Here the only site that 
shows some discrepancy if the site number 767. If that site is in Greenland where nobody lives one can decide that the 
approximation is good anyway ;-) The engine also report the RMS-differences by considering all the sites, i.e.::

	rms-diff = sqrt<(hmap1 - hmap2)^2> # mediating on all the sites

As to be expected, the differences are larger for maps with a smaller poe, i.e. a larger return period. But even in the 
worst case the RMS difference is only of 1E-3 g, which is not much. The complete calculation will have more sites, so 
the RMS difference will likely be even smaller. If you can check the few outlier sites and convince yourself that they 
are not important, you have succeeded in doubling the speed on your computation. And then you can start to work on the 
other quadratic and linear parameter and to get an ever bigger speedup!

******************************
Collapsing the GMPE logic tree
******************************

Some hazard models have GMPE logic trees which are insanely large. For instance the GMPE logic tree for the latest 
European model (ESHM20) contains 961,875 realizations. This causes two issues:

1. it is impossible to run a calculation with full enumeration, so one must use sampling
2. when one tries to increase the number of samples to study the stability of the mean hazard curves, the calculation runs out of memory

Fortunately, it is possible to compute the exact mean hazard curves by collapsing the GMPE logic tree. This is a simple 
as listing the name of the branchsets in the GMPE logic tree that one wants to collapse. For instance in the case of 
ESHM20 model there are the following 6 branchsets:

1. Shallow_Def (19 branches)
2. CratonModel (15 branches)
3. BCHydroSubIF (15 branches)
4. BCHydroSubIS (15 branches)
5. BCHydroSubVrancea (15 branches)
6. Volcanic (1 branch)

By setting in the job.ini the following parameters::

	number_of_logic_tree_samples = 0
	collapse_gsim_logic_tree = Shallow_Def CratonModel BCHydroSubIF BCHydroSubIS BCHydroSubVrancea Volcanic

it is possible to collapse completely the GMPE logic tree, i.e. going from 961,875 realizations to 1. Then the memory 
issues are solved and one can assess the correct values of the mean hazard curves. Then it is possible to compare with 
the value produce with sampling and assess how much they can be trusted.

NB: the ``collapse_gsim_logic_tree`` feature is rather old but only for engine versions >=3.13 it produces the exact 
mean curves (using the ``AvgPoeGMPE``); otherwise it will produce a different kind of collapsing (using the ``AvgGMPE``).

Disaggregation by source ``disagg_by_src``
------------------------------------------

Given a system of various sources affecting a specific site, one very common question to ask is: what are the more 
relevant sources, i.e. which sources contribute the most to the mean hazard curve? The engine is able to answer such 
question by setting the ``disagg_by_src`` flag in the job.ini file. When doing that, the engine saves in the datastore a 
4-dimensional ArrayWrapper called ``mean_rates_by_src`` with dimensions (site ID, intensity measure type, intensity measure 
level, source ID). From that it is possible to extract the contribution of each source to the mean hazard curve 
(interested people should look at the code in the function ``check_disagg_by_src``). The ArrayWrapper ``mean_rates_by_src`` 
can also be converted into a pandas DataFrame, then getting something like the following::

	>> dstore['mean_rates_by_src'].to_dframe().set_index('src_id')
	               site_id  imt  lvl         value
	ASCTRAS407           0  PGA    0  9.703749e-02
	IF-CFS-GRID03        0  PGA    0  3.720510e-02
	ASCTRAS407           0  PGA    1  6.735009e-02
	IF-CFS-GRID03        0  PGA    1  2.851081e-02
	ASCTRAS407           0  PGA    2  4.546237e-02
	...                ...  ...  ...           ...
	IF-CFS-GRID03        0  PGA   17  6.830692e-05
	ASCTRAS407           0  PGA   18  1.072884e-06
	IF-CFS-GRID03        0  PGA   18  1.275539e-05
	ASCTRAS407           0  PGA   19  1.192093e-07
	IF-CFS-GRID03        0  PGA   19  5.960464e-07

The ``value`` field here is the probability of exceedence in the hazard curve. The ``lvl`` field is an integer 
corresponding to the intensity measure level in the hazard curve.

In engine 3.15 we introduced the so-called “colon convention” on source IDs: if you have many sources that for some 
reason should be collected together - for instance because they all account for seismicity in the same tectonic region, 
or because they are components of a same source but are split into separate sources by magnitude - you can tell the 
engine to collect them into one source in the ``mean_rates_by_src`` matrix. The trick is to use IDs with the same 
prefix, a colon, and then a numeric index. For instance, if you had 3 sources with IDs ``src_mag_6.65``, ``src_mag_6.75``, 
``src_mag_6.85``, fragments of the same source with different magnitudes, you could change their IDs to something like 
``src:0``, ``src:1``, ``src:2`` and that would reduce the size of the matrix mean_rates_by_src by 3 times by collecting 
together the contributions of each source. There is no restriction on the numeric indices to start from 0, so using the 
names ``src:665``, ``src:675``, ``src:685`` would work too and would be clearer: the IDs should be unique, however.

If the IDs are not unique and the engine determines that the underlying sources are different, then an extension 
“semicolon + incremental index” is automatically added. This is useful when the hazard modeler wants to define a model 
where the more than one version of the same source appears in one source model, having changed some of the parameters, 
or when varied versions of a source appear in each branch of a logic tree. In that case, the modeler should use always 
the exact same ID (i.e. without the colon and numeric index): the engine will automatically distinguish the sources 
during the calculation of the hazard curves and consider them the same when saving the array ``mean_rates_by_src``: you 
can see an example in the test ``qa_tests_data/classical/case_20/job_bis.ini`` in the engine code base. In that case 
the ``source_info`` dataset will list 7 sources ``CHAR1;0 CHAR1;1 CHAR1;2 COMFLT1;0 COMFLT1;1 SFLT1;0 SFLT1;1`` but the 
matrix ``mean_rates_by_src`` will see only three sources ``CHAR1 COMFLT1 SFLT1`` obtained by composing together the 
versions of the underlying sources.

In version 3.15 ``mean_rates_by_src`` was extended to work with mutually exclusive sources, i.e. for the Japan model. 
You can see an example in the test ``qa_tests_data/classical/case_27``. However, the case of mutually exclusive ruptures 
- an example is the New Madrid cluster in the USA model - is not supported yet.

In some cases it is tricky to discern whether use of the colon convention or identical source IDs is appropriate. The 
following list indicates several possible cases that a user may encounter, and the appropriate approach to assigning 
source IDs. Note that this list includes the cases that have been tested so far, and is not a comprehensive list of all 
cases that may arise.

1. Sources in the same source group/source model are scaled
   alternatives of each other. For example, this occurs when for a
   given source, epistemic uncertainties such as occurrence rates or
   geometries are considered, but the modeller has pre-scaled the
   rates rather than including the alternative hypothesis in separate
   logic tree branches.

   **Naming approach**: identical IDs.

2. Sources in different files are alternatives of each other,
   e.g. each is used in a different branch of the source model logic
   tree.

   **Naming approach**: identical IDs.

3. A source is defined in OQ by numerous sources, either in the same
   file or different ones. For example, one could have a set of
   non-parametric sources, each with many ruptures, that are grouped
   together into single files by magnitude. Or, one could have many
   point sources that together represent the seismicity from one
   source.

   **Naming approach**: colon convention

4. One source consists of many mutually exclusive sources, as in qa_tests_data/classical/case_27.

   **Naming approach**: colon convention

Cases 1 and 2 could include include more than one source typology, as in ``qa_tests_data/classical/case_79``.

NB: ``disagg_by_src`` can be set to true only if the ``ps_grid_spacing`` approximation is disabled. The reason is that 
the ``ps_grid_spacing`` approximation builds effective sources which are not in the original source model, thus breaking 
the connection between the values of the matrix and the original sources.

The post-processing framework and Vector-valued PSHA calculations
-----------------------------------------------------------------

Since version 3.17 the OpenQuake engine has special support for custom postprocessors. A postprocessor is a Python 
module located in the directory ``openquake/calculators/postproc`` and containing a ``main`` function with signature::

	def main(dstore, [csm], ...):
	    ...

Post-processors are called after a classical or preclassical calculation: the ``dstore`` parameter is a DataStore 
instance corresponding to the calculation, while the ``csm`` parameter is a CompositeSourceModel instance (it can be 
omitted if not needed).

The ``main`` function is called when the user sets in the job.ini file the parameters ``postproc_func`` and ``postproc_args``. 
``postproc_func`` is the dotted name of the postprocessing function (in the form ``modulename.funcname`` where ``funcname`` 
is normally ``main``) and ``postproc_args`` is a dictionary of literal arguments that get passed to the function; if not 
specified the empty dictionary is passed. This happens for instance for the conditional spectrum post-processor since it 
does not require additional arguments with respect to the ones in ``dstore['oqparam']``.

The post-processing framework was put in place in order to run VPSHA calculations. The user can find an example in 
``qa_tests_data/postproc/case_mrd``. In the job.ini file there are the lines::

	postproc_func = compute_mrd.main
	postproc_args = {
	  'imt1': 'PGA',
	  'imt2': 'SA(0.05)',
	  'cross_correlation': 'BakerJayaram2008',
	  'seed': 42,
	  'meabins': [0.1, 0.2, 0.3, 0.4, 0.5, 0.6],
	  'sigbins': [0.2, 0.3, 0.4, 0.5, 0.6, 0.7],
	  'method': 'indirect'}

while the postprocessor module ``openquake.calculators.postproc.compute_mrd`` contains the function::

	# inside openquake.calculators.postproc.compute_mrd
	def main(dstore, imt1, imt2, cross_correlation, seed, meabins, sigbins,
	         method='indirect'):
	    ...

Inside ``main`` there is code to create the dataset ``mrd`` which contains the Mean Rate Distribution as an array of 
shape L1 x L1 x N where L1 is the number of levels per IMT minus 1 and N the number of sites (normally 1).

While the postprocessing part for VPSHA calculations is computationally intensive, it is much more common to have a 
light postprocessing, i.e. faster than the classical calculation it depends on. In such situations the postprocessing 
framework really shines, since it is possible to reuse the original calculation via the standard ``--hc`` switch, i.e. you 
can avoid repeating multiple times the same classical calculation if you are interested in running the postprocessor 
with different parameters. In that situation the ``main`` function will get a DataStore instance with an attribute ``parent`` 
corresponding to the DataStore of the original calculation.

The postprocessing framework also integrates very well with interactive development (think of Jupyter notebooks). The 
following lines are all you need to create a child datastore where the postprocessing function can store its results 
after reading the data from the calculation datastore::

	>> from openquake.commonlib.datastore import read, create_job_dstore
	>> from openquake.calculators.postproc import mypostproc
	>> log, dstore = create_job_dstore(parent=read(calc_id))
	>> with log:
	..     mypostproc.main(dstore)

***************************************
The conditional spectrum post-processor
***************************************

Since version 3.17 the engine includes an experimental post-processor which is able to compute the conditional spectrum.

The implementation was adapted from the paper *Conditional Spectrum Computation Incorporating Multiple Causal 
Earthquakes and Ground-Motion Prediction Models* by Ting Lin, Stephen C. Harmsen, Jack W. Baker, and Nicolas Luco 
(http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.845.163&rep=rep1&type=pdf) and it is rather sophisticated.

In order to perform a conditional spectrum calculation you need to specify, in addition to the usual parameter of a 
classical calculation:

1. a reference intensity measure type (i.e. imt_ref = SA(0.2))
2. a cross correlation model (i.e. cross_correlation = BakerJayaram2008)
3. a set of poes (i.e. poes = 0.01 0.1)

The engine will compute a mean conditional spectrum for each ``poe`` and site, as well as the usual mean uniform hazard 
spectra. The following restrictions are enforced:

1. the IMTs can only be of type ``SA`` and ``PGA``
2. the source model cannot contain mutually exclusive sources (i.e. you cannot compute the conditional spectrum for the Japan model)

An example can be found in the engine repository, in the directory openquake/qa_tests_data/conditional_spectrum/case_1. 
If you run it, you will get something like the following::

	$ oq engine --run job.ini
	...
	 id | name
	261 | Full Report
	262 | Hazard Curves
	260 | Mean Conditional Spectra
	263 | Realizations
	264 | Uniform Hazard Spectra

Exporting the output 260 will produce two files ``conditional-spectrum-0.csv`` and ``conditional-spectrum-1.csv``; the 
first will refer to the first poe, the second to the second poe. Each file will have a structure like the following::

	 #,,,,"generated_by='OpenQuake engine 3.13.0-gitd78d717e66', start_date='2021-10-13T06:15:20', checksum=3067457643, imls=[0.99999, 0.61470], site_id=0, lon=0.0, lat=0.0"
	sa_period,val0,std0,val1,std1
	0.00000E+00,1.02252E+00,2.73570E-01,7.53388E-01,2.71038E-01
	1.00000E-01,1.99455E+00,3.94498E-01,1.50339E+00,3.91337E-01
	2.00000E-01,2.71828E+00,9.37914E-09,1.84910E+00,9.28588E-09
	3.00000E-01,1.76504E+00,3.31646E-01,1.21929E+00,3.28540E-01
	1.00000E+00,3.08985E-01,5.89767E-01,2.36533E-01,5.86448E-01

The number of columns will depend from the number of sites. The conditional spectrum calculator, like the disaggregation 
calculator, is meant to be run on a very small number of sites, normally one. In this example there are two sites 0 and 1 
and the columns ``val0`` and ``val`` give the value of the conditional spectrum on such sites respectively, while the 
columns ``std0`` and ``std1`` give the corresponding standard deviations.

Conditional spectra for individual realizations are also computed and stored for debugging purposes, but they are not 
exportable.

Event Based and Scenarios
-------------------------

Scenario risk calculations usually do not pose a performance problem, since they involve a single rupture and limited 
geographical region for analysis. Some event-based risk calculations, however, may involve millions of ruptures and 
exposures spanning entire countries or even continents. This section offers some practical tips for running large 
event-based risk calculations, especially ones involving large logic trees, and proposes techniques that might be used 
to make an otherwise intractable calculation tractable.

************************
Understanding the hazard
************************

Event-based calculations are typically dominated by the hazard component (unless there are lots of assets aggregated on 
a few hazard sites) and therefore the first thing to do is to estimate the size of the hazard, i.e. the number of GMFs 
that will be produced. Since we are talking about a large calculation, first of all, we need to reduce it to a size 
that is guaranteed to run quickly. The simplest way to do that is to reduce the parameters directly affecting the 
number of ruptures generated, i.e.

- investigation_time
- ses_per_logic_tree_path
- number_of_logic_tree_samples

For instance, if you have ``ses_per_logic_tree_path = 10,000`` reduce it to 10, run the calculation and you will see in 
the log something like this::

	[2018-12-16 09:09:57,689 #35263 INFO] Received
	{'gmfdata': '752.18 MB', 'hcurves': '224 B', 'indices': '29.42 MB'}

The amount of GMFs generated for the reduced calculation is 752.18 MB; and since the calculation has been reduced by a 
factor of 1,000, the full computation is likely to generate around 750 GB of GMFs. Even if you have sufficient disk 
space to store this large quantity of GMFs, most likely you will run out of memory. Even if the hazard part of the 
calculation manages to run to completion, the risk part of the calculation is very likely to fail — managing 750 GB of 
GMFs is beyond the current capabilities of the engine. Thus, you will have to find ways to reduce the size of the 
computation.

A good start would be to carefully set the parameters ``minimum_magnitude`` and ``minimum_intensity``:

- ``minimum_magnitude`` is a scalar or a dictionary keyed by tectonic region; the engine will discard ruptures with magnitudes below the given thresholds
- ``minimum_intensity`` is a scalar or a dictionary keyed by the intensity measure type; the engine will discard GMFs below the given intensity thresholds

Choosing reasonable cutoff thresholds with these parameters can significantly reduce the size of your computation when 
there are a large number of small magnitude ruptures or low intensity GMFs being generated, which may have a negligible 
impact on the damage or losses, and thus could be safely discarded.

For instance, at the time of this writing (summer 2024) the hazard model for
South America (SAM) contains 539,831 sites and a full realistic simulation
for 100,000 years without setting ``minimum_magnitude`` and ``minimum_intensity``
would be impossible, generating an estimated 3.5 TB of ground motion fields.
In order to give this estimate, the trick is to reduce the effective investigation
time to something manageable (say 1000 years), run the calculation and measure
the size of the generated GMFs. The numbers are as follows::

 eff_time = 1000
 num_events = 159_559
 num_ruptures = 151_119
 mean mag = 4.50
 gmf_data = 35.04 GB
 time(event_based) = 3_725

It is clear that the calculation is dominated by the small magnitude ruptures which
however are expected to have a minimum impact on the damage. Usually we consider
a ``minimum_magnitude`` of 5; with this setting the situation improves drastically::

 eff_time = 1000
 minimum_magnitude = 5
 num_events = 20_514
 num_ruptures = 20_381
 mean mag = 5.68
 gmf_data = 4.39 GB
 time(event_based) = 467

We produce 8x less events and 8x less GMFs in 8x less time. Most sites, however, are affected
by very small shaking values, so setting ``minimum_intensity`` will reduce a lot
more the size of the GMFs (6.5x)::

 eff_time = 1000
 minimum_magnitude = 5
 minimum_intensity = .05
 num_events = 20_514
 num_ruptures = 20_381
 mean mag = 5.68
 gmf_data = 0.67 GB
 time(event_based) = 439

In this example with 1000 years setting ``minimum_magnitude`` and
``minimum_intensity`` reduced the size of the generated GMFs by 8 *
6.5 = 52 times (!) In the realistic case of 100,000 years the saving
will be even larger, i.e. you can easily gain two orders of magnitude
in the size of the generated GMFs. It is not only that: setting those
parameters can make the difference between being able to run the
calculation and running out of memory.
This is why, starting from engine v1.21, such parameters are mandatory except
in toy calculations.

*******************
region_grid_spacing
*******************

In our experience, the most common error made by our users is to compute the hazard at the sites of the exposure. The 
issue is that it is possible to have exposures with millions of assets on millions of distinct hazard sites. Computing 
the GMFs for millions of sites is hard or even impossible (there is a limit of 4 billion rows on the size of the GMF 
table in the datastore). Even in the cases when computing the hazard is possible, then computing the risk starting from 
an extremely large amount of GMFs will likely be impossible, due to memory/runtime constraints.

The second most common error is using an extremely fine grid for the site model. Remember that if you have a resolution 
of 250 meters, a square of 250 km x 250 km will contain one million sites, which is definitely too much. The engine was 
designed when the site models had resolutions around 5-10 km, i.e. of the same order of the hazard grid, while nowadays 
the vs30 fields have a much larger resolution.

Both problems can be solved in a simple way by specifying the ``region_grid_spacing`` parameter. Make it large enough 
that the resulting number of sites becomes reasonable and you are done. You will lose some precision, but that is 
preferable to not being able to run the calculation. You will need to run a sensitivity analysis with different values 
of region_grid_spacing parameter to make sure that you get consistent results, but that’s it.

Once a ``region_grid_spacing`` is specified, the engine computes the convex hull of the exposure sites and builds a 
grid of hazard sites, associating the site parameters from the closest site in the site model and discarding sites in 
the region where there are no assets (i.e. more distant than ``region_grid_spacing * sqrt(2)``). The precise logic is 
encoded in the function ``openquake.commonlib.readinput.get_sitecol_assetcol``, if you want to know the specific 
implementation details.

Our recommendation is to use the command ``oq prepare_site_model`` to apply such logic before starting a calculation and 
thus producing a custom site model file tailored to your exposure (see the section :ref:`prepare_site_model <prepare-site-model>`).

**********************
Collapsing of branches
**********************

When one is not interested in the uncertainty around the loss estimates and cares more about the mean estimates, all of 
the source model branches can be “collapsed” into one branch. Using the collapsed source model should yield the same 
mean hazard or loss estimates as using the full source model logic tree and then computing the weighted mean of the 
individual branch results.

Similarly, the GMPE logic tree for each tectonic region can also be “collapsed” into a single branch. Using a single 
collapsed GMPE for each TRT should also yield the same mean hazard estimates as using the full GMPE logic tree and then 
computing the weighted mean of the individual branch results. This has become possible through the introduction of 
`AvgGMPE <https://github.com/gem/oq-engine/blob/engine-3.9/openquake/qa_tests_data/classical/case_19/gmpe_logic_tree.xml#L26-L40>`_ feature in version 3.9.

*****************************************
Splitting the calculation into subregions
*****************************************

If one is interested in propagating the full uncertainty in the source models or ground motion models to the hazard or 
loss estimates, collapsing the logic trees into a single branch to reduce computational expense is not an option. But 
before going through the effort of trimming the logic trees, there is an interim step that must be explored, at least 
for large regions, like the entire continental United States. This step is to geographically divide the large region 
into logical smaller subregions, such that the contribution to the hazard or losses in one subregion from the other 
subregions is negligibly small or even zero. The effective realizations in each of the subregions will then be much 
fewer than when trying to cover the entire large region in a single calculation.

*******************************************************
Trimming of the logic-trees or sampling of the branches
*******************************************************

Trimming or sampling may be necessary if the following two conditions hold:

1. You are interested in propagating the full uncertainty to the hazard and loss estimates; only the mean or quantile results are not sufficient for your analysis requirements, AND
2. The region of interest cannot be logically divided further as described above; the logic-tree for your chosen region of interest still leads to a very large number of effective realizations.

Sampling is the easier of the two options now. You only need to ensure that you sample a sufficient number of branches 
to capture the underlying distribution of the hazard or loss results you are interested in. The drawback of random 
sampling is that you may still need to sample hundreds of branches to capture well the underlying distribution of the 
results.

Trimming can be much more efficient than sampling, because you pick a few branches such that the distribution of the 
hazard or loss results obtained from a full-enumeration of these branches is nearly the same as the distribution of the 
hazard or loss results obtained from a full-enumeration of the entire logic-tree.

Extra tips specific to event based calculations
-----------------------------------------------

Event based calculations differ from classical calculations because they produce visible ruptures, which can be 
exported and made accessible to the user. In classical calculations, instead, the underlying ruptures only live in 
memory and are normally not saved in the datastore, nor are exportable. The limitation is fundamentally a technical one: 
in the case of an event based calculation only a small fraction of the ruptures contained in a source are actually 
generated, so it is possible to store them. In a classical calculation all ruptures are generated and there are so many 
millions of them that it is impractical to save them, unless there are very few sites. For this reason they live in 
memory, they are used to produce the hazard curves and immediately discarded right after. The exception if for the case 
of few sites, i.e. if the number of sites is less than the parameter ``max_sites_disagg`` which by default is 10.

***************************************************
Convergency of the GMFs for non-trivial logic trees
***************************************************

In theory, the hazard curves produced by an event based calculation should converge to the curves produced by an 
equivalent classical calculation. In practice, if the parameters ``number_of_logic_tree_samples`` and 
``ses_per_logic_tree_path`` (the product of them is the relevant one) are not large enough they may be different. The 
engine is able to compare the mean hazard curves and to see how well they converge. This is done automatically if the 
option ``mean_hazard_curves = true`` is set. Here is an example of how to generate and plot the curves for one of our 
QA tests (a case with bad convergence was chosen on purpose)::

	$ oq engine --run event_based/case_7/job.ini
	<snip>
	WARNING:root:Relative difference with the classical mean curves for IMT=SA(0.1): 51%
	WARNING:root:Relative difference with the classical mean curves for IMT=PGA: 49%
	<snip>
	$ oq plot /tmp/cl/hazard.pik /tmp/hazard.pik --sites=0,1,2

.. figure:: _images/ebcl-convergency.png

The relative difference between the classical and event based curves is computed by computing the relative difference 
between each point of the curves for each curve, and by taking the maximum, at least for probabilities of exceedence 
larger than 1% (for low values of the probability the convergency may be bad). For the details I suggest you to look at 
the code.

***************************************************
Risk profiles
***************************************************

The OpenQuake engine can produce risk profiles, i.e. estimates of average losses
and maximum probable losses for all countries in the world. Even if you
are interested in a single country, you can still use this feature
to compute risk profiles for each province in your country.

However, the calculation of the risk profiles is tricky and there are
actually several different ways to do it.

1. The least-recommended way is to run independent calculations, one
   for each country. The issue with this approach is that even if the
   hazard model is the same for all the countries (say you are
   interested in the 13 countries of South America), due to the nature
   of event based calculations, different ruptures will be sampled in
   different countries. In practice, when comparing Chile with Peru you will
   see differences due to the fact that the random sampling picked
   different ruptures in the two contries and not real differences. In
   theory, the effect should disappear if the calculations have sufficiently
   long investigation times, when all possible ruptures are sampled, 
   but in practice, for finite investigation times there will always be 
   different ruptures.

2. To avoid such issues, the country-specific calculations should
   ideally all start from the same set of precomputed ruptures. You can
   compute the whole stochastic event set by running an event based
   calculation without specifying the sites and with the parameter
   ``ground_motion_fields`` set to false. Currently, one must specify
   a few global site parameters in the precalculation to make the
   engine checker happy, but they will not be used since the
   ground motion fields will not be generated in the
   precalculation. The ground motion fields will be generated on-the-fly  
   in the subsequent individual country calculations, but not stored 
   in the file system. This approach is fine if you do not have a lot of disk
   space at your disposal, but it is still inefficient since it is
   quite prone to the slow tasks issue.

3. If you have plenty of disk space it is better to also generate the
   ground motion fields in the precalculation, and then run the
   contry-specific risk calculations starting from there. This is
   particularly convenient if you foresee the need to run the risk
   part of the calculations multiple times, while the hazard part remains
   unchanged. Using a precomputed set of GMFs removes the need to rerun
   the hazard part of the calculations each time. This workflow has
   been particularly optimized since version 3.16 of the engine and it is
   now quite efficient.

4. If you have a really powerful machine, the simplest is to
   run a single calculation considering all countries in a single
   job.ini file. The risk profiles can be obtained by using the
   ``aggregate_by`` and ``reaggregate_by`` parameters. This approach
   can be much faster than the previous ones. However, approaches #2
   and #3 are cloud-friendly and can be preferred if you have access
   to cloud-computing resources, since then you can spawn a different
   machine for each country and parallelize horizontally.

Here are some tips on how to prepare the required job.ini files:

When using approach #1 you will have 13 different files (in the example
of South America) with a format like the following::

 $ cat job_Argentina.ini
 calculation_mode = event_based_risk
 source_model_logic_tree_file = ssmLT.xml
 gsim_logic_tree_file = gmmLTrisk.xml
 site_model_file = Site_model_Argentina.csv
 exposure_file = Exposure_Argentina.xml
 ...
 $ cat job_Bolivia.ini
 calculation_mode = event_based_risk
 source_model_logic_tree_file = ssmLT.xml
 gsim_logic_tree_file = gmmLTrisk.xml
 site_model_file = Site_model_Bolivia.csv
 exposure_file = Exposure_Bolivia.xml
 ...

Notice that the ``source_model_logic_tree_file`` and ``gsim_logic_tree_file``
will be the same for all countries since the hazard model is the same;
the same sources will be read 13 times and the ruptures will be sampled
and filtered 13 times. This is inefficient. Also, hazard parameters like

::

 truncation_level = 3
 investigation_time = 1
 number_of_logic_tree_samples = 1000
 ses_per_logic_tree_path = 100
 maximum_distance = 300

must be the same in all 13 files to ensure the consistency of the
calculation. Ensuring this consistency can be prone to human error.

When using approach #2 you will have 14 different files: 13 files for
the individual countries and a special file for precomputing the ruptures::

 $ cat job_rup.ini 
 calculation_mode = event_based
 source_model_logic_tree_file = ssmLT.xml
 gsim_logic_tree_file = gmmLTrisk.xml
 reference_vs30_value = 760
 reference_depth_to_1pt0km_per_sec = 440
 ground_motion_fields = false
 ...

The files for the individual countries will be as before, except for
the parameter ``source_model_logic_tree_file`` which should be
removed. That will avoid reading 13 times the same source model files,
which are useless anyway, since the calculation now starts from
precomputed ruptures. There are still a lot of repetitions in the
files and the potential for making mistakes.

Approach #3 is very similar to approach #2: the only differences will be
in the initial file, the one used to precompute the GMFs. Obviously it
will require setting ``ground_motion_fields = true``; moreover, it will
require specifying the full site model as follows::

  site_model_file =
    Site_model_Argentina.csv
    Site_model_Bolivia.csv
    ...

The engine will automatically concatenate the site model files for all
13 countries and produce a single site collection. The site parameters
will be extracted from such files, so the dummy global parameters
``reference_vs30_value``, ``reference_depth_to_1pt0km_per_sec``, etc
can be removed.

It is FUNDAMENTAL FOR PERFORMANCE to have reasonable site model files,
i.e. you should not compute the hazard at the location of every single asset,
but rather you should use a variable-size grid fitting the exposure.

The engine provides a command ``oq prepare_site_model``
which is meant to generate sensible site model files starting from
the country exposures and the global USGS vs30 grid.
It works by using a hazard grid so that the number of sites
can be reduced to a manageable number. Please refer to the manual in
the section about the oq commands to see how to use it, or try
``oq prepare_site_model --help``.

For reference, we were able to compute the hazard for all of South
America on a grid of half million sites and 1 million years of effective time
in a few hours in a machine with 120 cores, generating half terabyte of GMFs.
The difficult part is avoiding running out memory when running the risk
calculation; huge progress in this direction was made in version 3.16 of
the engine.

Approach #4 is the best, when applicable, since there is only a single file,
thus avoiding entirely the possibily of having inconsistent parameters
in different files. It is also the faster approach, not to mention the
most convenient one, since you have to manage a single calculation and
not 13. That makes the task of managing any kind of post-processing a lot
simpler. Unfortunately, it is also the option that requires more
memory and it can be infeasable if the model is too large and you do not
have enough computing resources. In that case your best bet might be to
go back to options #2 or #3. If you have access to multiple small machines,
approaches #2 and #3 can be more attractive than #4, since then you 
can scale horizontally. If you decide to use approach #4, 
in the single file you must specify the ``site_model_file`` as done in
approach #3, and also the ``exposure_file`` as follows::

 exposure_file =
   Exposure_Argentina.xml
   Exposure_Bolivia.xml
   ...

The engine will automatically build a single asset collection for the
entire continent of South America. In order to use this approach, you
need to collect all the vulnerability functions in a single file and
the taxonomy mapping file must cover the entire exposure for all
countries.  Moreover, the exposure must contain the associations
between asset<->country; in GEM's exposure models, this is typically
encoded in a field called ``ID_0``. Then the aggregation by country
can be done with the option

::

   aggregate_by = ID_0

Sometimes, one is interested in finer aggregations, for instance by country
and also by occupancy (Residential, Industrial or Commercial); then you have
to set

::

 aggregate_by = ID_0, OCCUPANCY
 reaggregate_by = ID_0

``reaggregate_by`` is a new feature of engine 3.13 which allows to go
from a finer aggregation (i.e. one with more tags, in this example 2)
to a coarser aggregation (i.e. one with fewer tags, in this example 1).
Actually the command ``oq reaggregate`` has been there for more than one
year; the new feature is that it is automatically called at the end of
a calculation, by spawning a subcalculation to compute the reaggregation.
Without ``reaggregate_by`` the aggregation by country would be lost,
since only the result of the finer aggregation would be stored.

Single-line commands
--------------------

When using approach #1 your can run all of the required calculations
with the command::

 $ oq engine --run job_Argentina.csv job_Bolivia.csv ...

When using approach #2 your can run all of the required calculations
with the command::

 $ oq engine --run job_rup.ini job_Argentina.csv job_Bolivia.csv ...

When using approach #3 your can run all of the required calculations
with the command::

 $ oq engine --run job_gmf.ini job_Argentina.csv job_Bolivia.csv ...

When using approach #4 your can run all of the required calculations
with the command::

 $ oq engine --run job_all.ini

Here ``job_XXX.ini`` are the country specific configuration files,
``job_rup.ini`` is the file generating the ruptures, ``job_rup.ini``
is the file generating the ruptures, ``job_gmf.ini`` is the file
generating the ground motion files and ``job_all.ini`` is the
file encompassing all countries.

Finally, if you have a file ``job_haz.ini`` generating the full GMFs,
a file ``job_weak.ini`` generating the losses with a weak building code
and a file ``job_strong.ini`` generating the losses with a strong building
code, you can run the entire an analysis with a single command as follows::

 $ oq engine --run job_haz.ini job_weak.ini job_strong.ini

This will generate three calculations and the GMFs will be reused.
This is as efficient as possible for this kind of problem.

Caveat: GMFs are split-dependent
--------------------------------

You should not expect the results of approach #4 to match exactly the
results of approaches #3 or #2, since splitting a calculation by
countries is a tricky operation. In general, if you have a set of
sites and you split it in disjoint subsets, and then you compute the
ground motion fields for each subset, you will get different results
than if you do not split.

To be concrete, if you run a calculation for Chile and then one for
Argentina, you will get different results than running a single
calculation for Chile+Argentina, *even if you have precomputed the
ruptures for both countries, even if the random seeds are the same and
even if there is no spatial correlation*. Many users are surprised but
this fact, but it is obvious if you know how the GMFs are
computed. Suppose you are considering 3 sites in Chile and 2 sites in
Argentina, and that the value of the random seed in 123456: if you
split, assuming there is a single event, you will produce the
following 3+2 normally distributed random numbers:

>>> np.random.default_rng(123456).normal(size=3)  # for Chile
array([ 0.1928212 , -0.06550702,  0.43550665])
>>> np.random.default_rng(123456).normal(size=2)  # for Argentina
array([ 0.1928212 , -0.06550702])

If you do not split, you will generate the following 5 random numbers
instead:

>>> np.random.default_rng(123456).normal(size=5)
array([ 0.1928212 , -0.06550702,  0.43550665,  0.88235875,  0.37132785])

They are unavoidably different. You may argue that not splitting is
the correct way of proceeding, since the splitting causes some
random numbers to be repeated (the numbers 0.1928212 and -0.0655070
in this example) and actually breaks the normal distribution.

In practice, if there is a sufficiently large event-set and if you are
interested in statistical quantities, things work out and you should
see similar results with and without splitting. But you will
*never produce identical results*. Only the classical calculator does
not depend on the splitting of the sites, for event based and scenario
calculations there is no way out.