diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 507b0804700d..ea9eb835e773 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -326,7 +326,7 @@ jobs: echo "Upload OpenQuakeforAdvancedUsers (PDF)" rsync -e 'ssh -o StrictHostKeyChecking=no -o UserKnownHostsFile=/dev/null' -ax --delete build/latex/OpenQuakeforAdvancedUsers.pdf "docs@docs.openquake.org:${DOCS_BASE}${TARGET}PDF/" - ssh -v -o StrictHostKeyChecking=no -o UserKnownHostsFile=/dev/null docs@docs.openquake.org "bash -cx 'cd \"${DOCS_BASE}${PDFDOCS}\"; ln -sf \"../${TARGET}PDF/OpenQuakeforAdvancedUsers.pdf\" \"OpenQuake for Advanced Users ${PDF_VER}.pdf\" .'" + ssh -v -o StrictHostKeyChecking=no -o UserKnownHostsFile=/dev/null docs@docs.openquake.org "bash -cx 'cd \"${DOCS_BASE}${PDFDOCS}\"; ln -sf \"../${TARGET}PDF/OpenQuakeforAdvancedUsers.pdf\" \"OpenQuake for Advanced Users ${PDF_VER}.pdf\"'" echo "Make Advanced Manual (HTML)" cd ../adv-manual && make clean && make html diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index f036fcd27649..29bcf0967c7f 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -76,6 +76,7 @@ Julián Montejo Espitia May 2022 Pablo Iturrieta (@pabloitu) June 2022 Anne Hulsey (@annehulsey) July 2022 Nicole Paul (@nicolepaul) May 2023 +Claudio Schill (@claudio525) June 2023 Project management ------------------ diff --git a/debian/changelog b/debian/changelog index 0202f6d85156..851964e03b59 100644 --- a/debian/changelog +++ b/debian/changelog @@ -1,3 +1,4 @@ + [Athanasios Papadopoulos] * Adjusted the Swiss-specific implementations of the GMPEs used in the Swiss national seismic hazard (SUIhaz15) and risk (ERM-CH23) models. @@ -11,7 +12,18 @@ ch_phisss03, ch_phis2s06 site parameters to the site.py file. These parameters are required by the apply_swiss_amplification_sa method. + [Marco Pagani, Michele Simionato] + * Added support for the "direct" method in MRD calculations + + [Michele Simionato] + * Reimplemented the conditional_spectrum calculator as a post-processor + + [Paolo Tormene, Michele Simionato] + * Raise an error when the same parameter is set in different sections + of the job.ini file + [Michele Simionato] + * Fixed tau-phi inversion in LanzanoLuzi2019 * Fixed another bug with conditioned GMFs appearing as the error `array[m, n, bad] = exp(mean_covs[0, g, m, n], im) TypeError: list indices must be integers or slices, not tuple` diff --git a/debian/control b/debian/control index 6b36c7154153..a05d5f724268 100644 --- a/debian/control +++ b/debian/control @@ -10,13 +10,13 @@ Homepage: http://github.com/gem/oq-engine Package: python3-oq-engine Architecture: amd64 X-Python-Version: >= 3.6 -Depends: ${shlibs:Depends}, ${misc:Depends}, ${python3:Depends}, ${dist:Depends}, python3-oq-libs (>=3.8.0~) +Depends: ${shlibs:Depends}, ${misc:Depends}, ${python3:Depends}, ${dist:Depends}, python3-oq-libs (>=3.8.1~) Conflicts: python-noq, python-oq, python-nhlib Replaces: python-oq-risklib, python-oq-hazardlib, python-oq-engine (<< ${binary:Version}) # See https://www.debian.org/doc/debian-policy/ch-relationships.html#s-breaks Breaks: python-oq-risklib, python-oq-hazardlib, python-oq-engine (<< ${binary:Version}) Recommends: ${dist:Recommends} -Suggests: python3-oq-libs-extra (>=3.8.0~) +Suggests: python3-oq-libs-extra (>=3.8.1~) Description: computes seismic hazard and physical risk based on the hazard and risk libraries developed by the GEM foundation (http://www.globalquakemodel.org/) diff --git a/debian/control.in b/debian/control.in index 9fd9bd24347e..7425637756c7 100644 --- a/debian/control.in +++ b/debian/control.in @@ -10,13 +10,13 @@ Homepage: http://github.com/gem/oq-engine Package: python3-oq-engine Architecture: amd64 X-Python-Version: >= 3.6 -Depends: ${shlibs:Depends}, ${misc:Depends}, ${python3:Depends}, ${dist:Depends}, python3-oq-libs (>=3.8.0~) +Depends: ${shlibs:Depends}, ${misc:Depends}, ${python3:Depends}, ${dist:Depends}, python3-oq-libs (>=3.8.1~) Conflicts: python-noq, python-oq, python-nhlib Replaces: python-oq-risklib, python-oq-hazardlib, python-oq-engine (<< ${binary:Version}) # See https://www.debian.org/doc/debian-policy/ch-relationships.html#s-breaks Breaks: python-oq-risklib, python-oq-hazardlib, python-oq-engine (<< ${binary:Version}) Recommends: ${dist:Recommends} -Suggests: python3-oq-libs-extra (>=3.8.0~) +Suggests: python3-oq-libs-extra (>=3.8.1~) Description: computes seismic hazard and physical risk based on the hazard and risk libraries developed by the GEM foundation (http://www.globalquakemodel.org/) diff --git a/doc/adv-manual/classical_PSHA.rst b/doc/adv-manual/classical_PSHA.rst index 56bfe93911f8..f0928708dfac 100644 --- a/doc/adv-manual/classical_PSHA.rst +++ b/doc/adv-manual/classical_PSHA.rst @@ -851,84 +851,74 @@ are the same. 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 5-dimensional array called ``disagg_by_src`` with -dimensions (site ID, realization ID, intensity measure type, -intensity measure level, source ID). For 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 array ``disagg_by_src`` can also be read as a pandas DataFrame, -then getting something like the following:: - - >> dstore.read_df('disagg_by_src', index='src_id') - site_id rlz_id imt lvl value - ASCTRAS407 0 0 PGA 0 9.703749e-02 - IF-CFS-GRID03 0 0 PGA 0 3.720510e-02 - ASCTRAS407 0 0 PGA 1 6.735009e-02 - IF-CFS-GRID03 0 0 PGA 1 2.851081e-02 - ASCTRAS407 0 0 PGA 2 4.546237e-02 - ... ... ... ... ... ... - IF-CFS-GRID03 0 31 PGA 17 6.830692e-05 - ASCTRAS407 0 31 PGA 18 1.072884e-06 - IF-CFS-GRID03 0 31 PGA 18 1.275539e-05 - ASCTRAS407 0 31 PGA 19 1.192093e-07 - IF-CFS-GRID03 0 31 PGA 19 5.960464e-07 +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. -There is a consistency check comparing the mean hazard curves -with the value obtained by composing the probabilities in the -disagg_by_src` array, for the heighest level of each intensity measure type. - -It should be noticed that many hazard models contain thousands of -sources and as a consequence the ``disagg_by_src`` matrix can be -impossible to compute without running out of memory. Even if you have -enough memory, having a very large ``disagg_by_src`` matrix is a bad -idea, so there is a limit on the size of the matrix, hard-coded to 4 -GB. The way to circumvent the limit is to reduce the number of -sources: for instance you could convert point sources in multipoint -sources. - -In engine 3.15 we also 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 ``disagg_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 ``disagg_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 ``disagg_by_src``: you can see an example in the -test ``qa_tests_data/classical/case_79`` in the engine code base. In that -case the ``source_info`` dataset will list 6 sources -``ASCTRAS407;0``, ``ASCTRAS407;1``, ``ASCTRAS407;2``, ``ASCTRAS407;3``, -``IF-CFS-GRID03;0``, ``IF-CFS-GRID03;1`` but the matrix ``disagg_by_src`` -will see only two sources ``ASCTRAS407`` and ``IF-CFS-GRID03`` obtained -by composing together the versions of the underlying sources. - -In version 3.15 ``disagg_by_src`` was extended to work with mutually +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 @@ -975,15 +965,102 @@ NB: ``disagg_by_src`` can be set to true only if the not in the original source model, thus breaking the connection between the values of the matrix and the original sources. -The conditional spectrum calculator -======================================== +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: + +.. code-block:: + + 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 istance 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 + +.. code-block:: + + # 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: + +.. code-block:: + + >> from openquake.commonlib.datastore import read, build_dstore_log + >> from openquake.calculators.postproc import mypostproc + >> dstore, log = build_dstore_log(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 ``conditional_spectrum`` calculator is an experimental calculator -introduced in version 3.13, which is able to compute the conditional -spectrum in the sense of Baker. +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 (on top of the usual parameter of a classical calculation): +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``) @@ -1033,15 +1110,3 @@ standard deviations. Conditional spectra for individual realizations are also computed and stored for debugging purposes, but they are not exportable. - -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. -The core formula is implemented in the method -`openquake.hazardlib.contexts.get_cs_contrib`. - -The ``conditional_spectrum`` calculator, like the disaggregation calculator, -is a kind of post-calculator, i.e. you can run a regular classical calculation -and then compute the ``conditional_spectrum`` in post-processing by using -the ``--hc`` option. diff --git a/doc/adv-manual/event_based.rst b/doc/adv-manual/event_based.rst index cef853dcf339..0a8b0fba0a8b 100644 --- a/doc/adv-manual/event_based.rst +++ b/doc/adv-manual/event_based.rst @@ -710,8 +710,8 @@ Rupture sampling: how to get it wrong Rupture samplings is *much more complex than one could expect* and in many respects *surprising*. In the many years of existence of the -engine, multiple approached were tried and you can expect some -detail of the rupture sampling mechanism to be different nearly at every +engine, multiple approached were tried and you can expect the +details of the rupture sampling mechanism to be different nearly at every version of the engine. Here we will discuss some tricky points that may help you understand @@ -721,8 +721,10 @@ rupture sampling is nontrivial. We will start with the first subtlety, the *interaction between sampling and filtering*. The short version is that you should *first -sample and then filter*. Here is the long version. Consider the -following code emulating rupture sampling for poissonian ruptures: +sample and then filter*. + +Here is the long version. Consider the following code emulating +rupture sampling for poissonian ruptures: .. code-block:: @@ -756,7 +758,7 @@ feature and it is able to discard ruptures below the minimum magnitude. But how should it work? The natural approach to follow, for performance-oriented applications, would be to first discard the low magnitudes and then perform the sampling. However, that would have effects that would be surprising -for many users. Consider the following two alternative: +for most users. Consider the following two alternative: .. code-block:: @@ -791,17 +793,12 @@ are consistent with the no-filtering case: >> calc_n_occ_before_filtering(fake_ruptures, eff_time, seed, min_mag) [ 9 6 13 7 6 6 10] -Historically the engine followed the filter-early approach and it is in -the process of transiting to the filter-late approach. The transition is -extremely tricky and error prone, since the minimum magnitude feature is -spread across many modules, and it could easily go on for years. - The problem with the filtering is absolutely general and not restricted only to the magnitude filtering: it is exactly the same also for distance filtering. Suppose you have a ``maximum_distance`` of 300 km and than you decide that you want to increase it to 301 km. One would expect this -change to have a minor impact; instead, the occupation numbers will be -completely different and you sample a very different set of ruptures. +change to have a minor impact; instead, you may end up sampling a very +different set of ruptures. It is true that average quantities like the hazard curves obtained from the ground motion fields will converge for long enough effective time, @@ -809,12 +806,13 @@ however in practice you are always in situations were 1. you cannot perform the calculation for a long enough effective time since it would be computationally prohibitive -2. you are interested on quantities which are strongly sensitive to a - change in the sampling, like the Maximum Probable Loss at some return period +2. you are interested on quantities which are strongly sensitive to aany + change, like the Maximum Probable Loss at some return period In such situations changing the site collection (or changing -the maximum distance which is akin to changing the site collection) can change -the sampling of the ruptures significantly. +the maximum distance which is akin to changing the site collection) +can change the sampling of the ruptures significantly, at least for +engine versions lower than 3.17. Users wanting to compare the GMFs or the risk on different site collections should be aware of this effect; the solution is to first sample the diff --git a/doc/adv-manual/general.rst b/doc/adv-manual/general.rst index 703ed81f3913..fbac12cc1313 100644 --- a/doc/adv-manual/general.rst +++ b/doc/adv-manual/general.rst @@ -129,7 +129,10 @@ magnitudes between 5 and 6 (the maximum distance in this magnitude range will vary linearly between 0 and 100), keep ruptures up to 200 km for magnitudes between 6 and 7 (with `maximum_distance` increasing linearly from 100 to 200 km from magnitude 6 to magnitude 7), keep -ruptures up to 300 km for magnitudes over 7. +ruptures up to 300 km for magnitudes between 7 and 8 (with +`maximum_distance` increasing linearly from 200 to 300 km from +magnitude 7 to magnitude 8) and discard ruptures for magnitudes +over 8. You can have both trt-dependent and mag-dependent maximum distance:: @@ -315,7 +318,7 @@ to know git and the tools of Open Source development in general, in particular about testing. If this is not the case, you should do some study on your own and come back later. There is a huge amount of resources on the net about these topics. This -manual will focus solely on the OpenQuake engine and it assume that +manual will focus solely on the OpenQuake engine and it assumes that you already know how to use it, i.e. you have read the User Manual first. @@ -359,14 +362,14 @@ Understanding the engine ------------------------- Once you have the engine installed you can run calculations. We recommend -starting from the demos directory which contains example of hazard and +starting from the demos directory which contains examples of hazard and risk calculations. For instance you could run the area source demo with the following command:: $ oq run demos/hazard/AreaSourceClassicalPSHA/job.ini You should notice that we used here the command ``oq run`` while the engine -manual recommend the usage of ``oq engine --run``. There is no contradiction. +manual recommends the usage of ``oq engine --run``. There is no contradiction. The command ``oq engine --run`` is meant for production usage, but here we are doing development, so the recommended command is ``oq run`` which will will be easier to debug thanks to the @@ -387,7 +390,7 @@ For that one can understand the bottlenecks of a calculation and, with experience, he can understand which part of the engine he needs to optimize. Also, it is very useful to play with the parameters of the calculation (like the maximum distance, the area discretization, the magnitude binning, -etc etc) and see how the performance change. There is also a command to +etc etc) and see how the performance changes. There is also a command to plot hazard curves and a command to compare hazard curves between different calculations: it is common to be able to get big speedups simply by changing the input parameters in the `job.ini` of the model, without changing much the @@ -544,7 +547,7 @@ The ``ContextMaker`` takes care of the maximum_distance filtering, so in general the number of contexts is lower than the total number of ruptures, since some ruptures are normally discarded, being distant from the sites. -The contexts contains all the rupture, site and distance parameters. +The contexts contain all the rupture, site and distance parameters. Then you have @@ -557,7 +560,7 @@ Then you have In this example, the GMPE ``ToroEtAl2002SHARE`` does not require site parameters, so calling ``ctx.vs30`` will raise an ``AttributeError`` -but in general the contexts contains also arrays of site parameters. +but in general the contexts contain also arrays of site parameters. There is also an array of indices telling which are the sites affected by the rupture associated to the context: @@ -612,9 +615,9 @@ table has a structure like this:: 2.0,5.0,MEAN,6.43850933e-03,3.61047741e-04,4.57949482e-04,7.24558049e-04,9.44495571e-04,5.11252304e-04,2.21076069e-04,4.73435138e-05 ... -The columns starting with ``rup_`` contains rupture parameters (the +The columns starting with ``rup_`` contain rupture parameters (the magnitude in this example) while the columns starting with ``dist_`` -contains distance parameters. The column ``result_type`` is a string +contain distance parameters. The column ``result_type`` is a string in the set {"MEAN", "INTER_EVENT_STDDEV", "INTRA_EVENT_STDDEV", "TOTAL_STDDEV"}. The remaining columns are the expected results for each intensity measure type; in the the example the IMTs are PGV, PGA, @@ -659,7 +662,7 @@ code in https://github.com/gem/oq-engine/blob/master/openquake/hazardlib/tests/g Running the engine tests ---------------------------------- -If you are a hazard scientists contributing a bug fix to a GMPE (or any +If you are a hazard scientist contributing a bug fix to a GMPE (or any other kind of bug fix) you may need to run the engine tests and possibly change the expected files if there is a change in the numbers. The way to do it is to start the dbserver and then run the tests from the repository root:: @@ -684,7 +687,7 @@ Architecture of the OpenQuake engine The engine is structured as a regular scientific application: we try to perform calculations as much as possible in memory and when it is not possible intermediate results are stored in HDF5 format. -We try work as much as possible in terms of arrays which are +We try to work as much as possible in terms of arrays which are efficiently manipulated at C/Fortran speed with a stack of well established scientific libraries (numpy/scipy). @@ -815,8 +818,8 @@ same as testability: it is essential to have a suite of tests checking that the calculators are providing the expected outputs against a set of predefined inputs. Currently we have thousands of tests which are run multiple times per day in our Continuous Integration environments -r(avis, GitLab, Jenkins), split into unit tests, end-to-end tests and -long running tests. +(GitHub Actions, GitLab Pipelines), split into unit +tests, end-to-end tests and long running tests. How the parallelization works in the engine =========================================== @@ -824,8 +827,7 @@ How the parallelization works in the engine The engine builds on top of existing parallelization libraries. Specifically, on a single machine it is based on multiprocessing, which is part of the Python standard library, while on a cluster -it is based on the combination celery/rabbitmq, which are well-known -and maintained tools. +it is based on zeromq, which is a well-known and maintained library. While the parallelization used by the engine may look trivial in theory (it only addresses embarrassingly parallel problems, not true concurrency) @@ -836,18 +838,17 @@ without affecting other calculations that may be running concurrently. Because of this requirement, we abandoned `concurrent.futures`, which is also in the standard library, but is lacking the ability to kill the pool of processes, which is instead available in multiprocessing with the -`Pool.shutdown` method. For the same reason, we discarded `dask`, which -is a lot more powerful than `celery` but lacks the revoke functionality. +`Pool.shutdown` method. For the same reason, we discarded `dask`. Using a real cluster scheduling mechanism (like SLURM) would be of course better, but we do not want to impose on our users a specific cluster -architecture. celery/rabbitmq have the advantage of being simple to install +architecture. Zeromq has the advantage of being simple to install and manage. Still, the architecture of the engine parallelization library -is such that it is very simple to replace celery/rabbitmq with other +is such that it is very simple to replace zeromq with other parallelization mechanisms: people interested in doing so should just contact us. -Another tricky aspects of parallelizing large scientific calculations is +Another tricky aspect of parallelizing large scientific calculations is that the amount of data returned can exceed the 4 GB limit of Python pickles: in this case one gets ugly runtime errors. The solution we found is to make it possible to yield partial results from a task: in this way instead of @@ -855,7 +856,7 @@ returning say 40 GB from one task, one can yield 40 times partial results of 1 GB each, thus bypassing the 4 GB limit. It is up to the implementor to code the task carefully. In order to do so, it is essential to have in place some monitoring mechanism measuring how much data is returned -back from a task, as well as other essential informations like how much +back from a task, as well as other essential information like how much memory is allocated and how long it takes to run a task. To this aim the OpenQuake engine offers a ``Monitor`` class @@ -895,20 +896,14 @@ hazard curves, because the algorithm is considering one rupture at the time and it is not accumulating ruptures in memory, differently from what happens when sampling the ruptures in event based. -If you have ever coded in celery, you will see that the OpenQuake engine -concept of task is different: there is no ``@task`` decorator and while at the -end engine tasks will become celery tasks this is hidden to the developer. -The reason is that we needed a celery-independent abstraction layer to make -it possible to use different kinds of parallelization frameworks/ - From the point of view of the coder, in the engine there is no difference -between a task running on a cluster using celery and a task running locally +between a task running on a cluster using zeromq and a task running locally using ``multiprocessing.Pool``: they are coded the same, but depending -on a configuration parameter in openquake.cfg (``distribute=celery`` +on a configuration parameter in openquake.cfg (``distribute=zmq`` or ``distribute=processpool``) the engine will treat them differently. You can also set an environment variable ``OQ_DISTRIBUTE``, which takes the precedence over openquake.cfg, to specify which kind of distribution -you want to use (``celery`` or ``processpool``): this is mostly used +you want to use (``zmq`` or ``processpool``): this is mostly used when debugging, when you typically insert a breakpoint in the task and then run the calculation with @@ -946,20 +941,19 @@ a single writer). This approach bypasses all the limitations of the SWMR mode in HDF5 and did not require a large refactoring of our existing code. -Another tricky point in cluster situations is that rabbitmq is not good -at transferring gigabytes of data: it was meant to manage lots of small -messages, but here we are perverting it to manage huge messages, i.e. -the large arrays coming from a scientific calculations. - -Hence, since recent versions of the engine we are no longer returning -data from the tasks via celery/rabbitmq: instead, we use zeromq. This +Another tricky point in cluster situations is that we need +to transfer gigabytes of data, i.e. the large arrays coming from +scientific calculations, rather than lots of small +messages. Hence, since recent versions of the engine, we are returning +data from the tasks via zeromq instead of using celery/rabbitmq as we +did in the past. This is hidden from the user, but internally the engine keeps track of all tasks that were submitted and waits until they send the message that -they finished. If one tasks runs out of memory badly and never sends +they finished. If one task runs out of memory badly and never sends the message that it finished, the engine may hang, waiting for the -results of a task that does not exist anymore. You have to be +results of a task that does not exist anymore. You have to be careful. What we did in our cluster is to set some memory limit on the -celery user with the cgroups technology, so that an out of memory task +openquake user with the cgroups technology, so that an out of memory task normally fails in a clean way with a Python MemoryError, sends the message that it finished and nothing particularly bad happens. Still, in situations of very heavy load the OOM killer may enter in action @@ -978,7 +972,7 @@ Suppose you want to code a character-counting algorithm, which is a textbook exercise in parallel computing and suppose that you want to store information about the performance of the algorithm. Then you should use the OpenQuake ``Monitor`` class, as well as the utility -``openquake.baselib.commonlib.hdf5new`` that build an empty datastore +``openquake.baselib.commonlib.hdf5new`` that builds an empty datastore for you. Having done that, the ``openquake.baselib.parallel.Starmap`` class can take care of the parallelization for you as in the following example: @@ -986,7 +980,7 @@ example: .. include:: char_count_par.py :code: python -The name ``Starmap`` was chosen it looks very similar to +The name ``Starmap`` was chosen because it looks very similar to how ``multiprocessing.Pool.starmap`` works, the only apparent difference being in the additional monitor argument:: @@ -994,8 +988,8 @@ being in the additional monitor argument:: In reality the ``Starmap`` has a few other differences: -1. it does not use the multiprocessing mechanism to returns back the results, - it uses zmq instead; +1. it does not use the multiprocessing mechanism to return back the results; + it uses zmq instead 2. thanks to that, it can be extended to generator functions and can yield partial results, thus overcoming the limitations of multiprocessing 3. the ``Starmap`` has a ``.submit`` method and it is actually more similar to @@ -1028,7 +1022,7 @@ parallelizing complex situations where it is expensive to use a single starmap. However, there is limit on the number of starmaps that can be alive at the same moment. -Moreover the ``Starmap`` has a ``.shutdown`` methods that allows to +Moreover the ``Starmap`` has a ``.shutdown`` method that allows to shutdown the underlying pool. The idea is to submit the text of each file - here I am considering .rst files, @@ -1047,7 +1041,7 @@ However, sometimes you *must* be tech-savvy: for instance if you want to post-process hundreds of GB of ground motion fields produced by an event based calculation, you should *not* use the CSV output, at least if you care about efficiency. To manage this case (huge amounts of data) -there a specific solution, which is also able to manage the case of data +there is a specific solution, which is also able to manage the case of data lacking a predefined exporter: the ``Extractor`` API. There are actually two different kind of extractors: the simple @@ -1094,7 +1088,7 @@ or even all realizations: >> dic = vars(extractor.get('hcurves?kind=rlzs')) Here is an example of using the `WebExtractor` to retrieve hazard maps. -Here we assumes that there is available in a remote machine where there is +Here we assume that in a remote machine there is a WebAPI server running, a.k.a. the Engine Server. The first thing to is to set up the credentials to access the WebAPI. There are two cases: @@ -1145,7 +1139,7 @@ Here are three examples of use:: $ oq plot 'hmaps?kind=mean&imt=PGA' $ oq plot 'uhs?kind=mean&site_id=0' -The ``site_id`` is optional; if missing only the first site (``site_id=0``) +The ``site_id`` is optional; if missing, only the first site (``site_id=0``) will be plotted. If you want to plot all the realizations you can do:: $ oq plot 'hcurves?kind=rlzs&imt=PGA' @@ -1159,7 +1153,7 @@ realizations and also the mean the command to give is:: $ oq plot 'hcurves?kind=rlzs&kind=mean&imt=PGA' -If you want to plot the mean and the median the command is:: +If you want to plot the median and the mean the command is:: $ oq plot 'hcurves?kind=quantile-0.5&kind=mean&imt=PGA' @@ -1175,7 +1169,7 @@ specifying a kind that is not available you will get an error. Extracting disaggregation outputs --------------------------------- -Disaggregation outputs are particularly complex and they are stored in +Disaggregation outputs are particularly complex and they are stored in the datastore in different ways depending on the engine version. Here we will give a few examples for the Disaggregation Demo, which has the flag ``individual_rlzs`` set. If you run the demos with a recent enough version @@ -1291,7 +1285,7 @@ Example: how many events per magnitude? When analyzing an event based calculation, users are often interested in checking the magnitude-frequency distribution, i.e. to count how many -events of a given magnitude present in the stochastic event set for +events of a given magnitude are present in the stochastic event set for a fixed investigation time and a fixed ``ses_per_logic_tree_path``. You can do that with code like the following: @@ -1308,7 +1302,7 @@ You can do that with code like the following: for mag, grp in events.groupby(['mag']): print(mag, len(grp)) # number of events per group -If you want to now the number of events per realization and per stochastic +If you want to know the number of events per realization and per stochastic event set you can just refine the `groupby` clause, using the list ``['mag', 'rlz_id', 'ses_id']`` instead of simply ``['mag']``. @@ -1484,20 +1478,20 @@ to use them. You can see the full list of commands by running `oq --help`:: - $ oq --help - usage: oq [--version] - {workerpool,webui,dbserver,info,ltcsv,dump,export,celery,plot_losses,restore,plot_assets,reduce_sm,check_input,plot_ac,upgrade_nrml,shell,plot_pyro,nrml_to,postzip,show,workers,abort,engine,reaggregate,db,compare,renumber_sm,download_shakemap,importcalc,purge,tidy,zip,checksum,to_hdf5,extract,reset,run,show_attrs,prepare_site_model,sample,plot} - ... + $ oq --help + usage: oq [-h] [-v] + {shell,upgrade_nrml,reduce_smlt,show_attrs,prepare_site_model,nrml_from,shakemap2gmfs,importcalc,run,show,purge,renumber_sm,workers,postzip,plot_assets,db,dbserver,tidy,extract,sample,to_hdf5,ltcsv,reaggregate,restore,mosaic,check_input,dump,info,zip,abort,nrml_to,engine,reset,checksum,export,webui,compare,plot,reduce_sm} + ... - positional arguments: - {workerpool,webui,dbserver,info,ltcsv,dump,export,celery,plot_losses,restore,plot_assets,reduce_sm,check_input,plot_ac,upgrade_nrml,shell,plot_pyro,nrml_to,postzip,show,workers,abort,engine,reaggregate,db,compare,renumber_sm,download_shakemap,importcalc,purge,tidy,zip,checksum,to_hdf5,extract,reset,run,show_attrs,prepare_site_model,sample,plot} - available subcommands; use oq --help + positional arguments: + {shell,upgrade_nrml,reduce_smlt,show_attrs,prepare_site_model,nrml_from,shakemap2gmfs,importcalc,run,show,purge,renumber_sm,workers,postzip,plot_assets,db,dbserver,tidy,extract,sample,to_hdf5,ltcsv,reaggregate,restore,mosaic,check_input,dump,info,zip,abort,nrml_to,engine,reset,checksum,export,webui,compare,plot,reduce_sm} + available subcommands; use oq --help - optional arguments: - -h, --help show this help message and exit - -v, --version show program's version number and exit + options: + -h, --help show this help message and exit + -v, --version show program's version number and exit -This is the output that you get at the present time (engine 3.11); depending +This is the output that you get at the present time (engine 3.17); depending on your version of the engine you may get a different output. As you see, there are several commands, like `purge`, `show_attrs`, `export`, `restore`, ... You can get information about each command with `oq --help`; @@ -1528,7 +1522,7 @@ features. logic tree of the calculation. 2. When invoked with the `--report` option, it produces a `.rst` report with -several important informations about the computation. It is ESSENTIAL in the +important information about the computation. It is ESSENTIAL in the case of large calculations, since it will give you an idea of the feasibility of the computation without running it. Here is an example of usage:: @@ -1536,7 +1530,7 @@ of the computation without running it. Here is an example of usage:: Generated /tmp/report_1644.rst -You can open `/tmp/report_1644.rst` and read the informations listed there +You can open `/tmp/report_1644.rst` and read the information listed there (`1644` is the calculation ID, the number will be different each time). 3. It can be invoked without a `job.ini` file, and it that case it provides @@ -1575,7 +1569,7 @@ The list of available exports (i.e. the datastore keys and the available export formats) can be extracted with the `oq info exports` command; the number of exporters defined changes at each version:: - $ oq info exports$ oq info exports + $ oq info exports ? "aggregate_by" ['csv'] ? "disagg_traditional" ['csv'] ? "loss_curves" ['csv'] @@ -1624,8 +1618,9 @@ a specific realization (say the first one), you can use:: but currently this only works for `csv` and `xml`. The exporters are one of the most time-consuming parts on the engine, mostly because of the sheer number -of them; the are more than fifty exporters and they are always increasing. -If you need new exports, please [add an issue on GitHub](https://github.com/gem/oq-engine/issues). +of them; there are more than fifty exporters and they are always increasing. +If you need new exports, +please `add an issue on GitHub `_. oq zip ------ diff --git a/doc/adv-manual/logic_trees.rst b/doc/adv-manual/logic_trees.rst index 1d52f15fc8b6..57c2e613bbed 100644 --- a/doc/adv-manual/logic_trees.rst +++ b/doc/adv-manual/logic_trees.rst @@ -80,10 +80,10 @@ and the 6 possible paths can be extracted as follows: The empty square brackets means that the branchset should be applied to all branches in the previous branchset and correspond to the ``applyToBranches`` -tag in the XML version of the logic tree. If ``applyToBranches`` is missing +tag in the XML version of the logic tree. If ``applyToBranches`` is missing, the logic tree is multiplicative and the total number of paths can be obtained simply by multiplying the number of paths in each branchset. -When ``applyToBranches`` is used the logic tree becomes additive and the +When ``applyToBranches`` is used, the logic tree becomes additive and the total number of paths can be obtained by summing the number of paths in the different subtrees. For instance, let us extend the previous example by adding another ``extendModel`` branchset and by using ``applyToBranches``: @@ -262,7 +262,7 @@ export it, you will get a CSV file with the following structure:: 323,ACCCC~BB,3.1111853e-03 For each realization there is a ``branch_path`` string which is split in -two parts separated by a tilde. The left part describe the branches of +two parts separated by a tilde. The left part describes the branches of the source model logic tree and the right part the branches of the gmpe logic tree. In past versions of the engine the branch path was using directly the branch IDs, so it was easy to assess the correspondence @@ -322,8 +322,7 @@ and then, using the correspondence table ``abbrev->uvalue``, to:: "[ChiouYoungs2008]" "[ToroEtAl2002]" For convenience, the engine provides a simple command to display the content -of a realization, given the realization number, thus answering the -FAQ:: +of a realization, given the realization number:: $ oq show rlz:322 | uncertainty_type | uvalue | @@ -402,7 +401,7 @@ Reduction of the logic tree The simplest case of logic tree reduction is when the actual sources do not span the full range of tectonic region types in the GMPE logic tree file. This happens very often. -FOr instance, in the SHARE calculation for Europe +For instance, in the SHARE calculation for Europe the GMPE logic tree potentially contains 1280 realizations coming from 7 different tectonic region types: diff --git a/doc/installing/README.md b/doc/installing/README.md index cfe46dcd5cc7..f245a4e4b0ce 100644 --- a/doc/installing/README.md +++ b/doc/installing/README.md @@ -37,7 +37,7 @@ Check more advanced [hardware suggestions here](./hardware-suggestions.md). **On Windows** - Download OpenQuake Engine for Windows: https://downloads.openquake.org/pkgs/windows/oq-engine/OpenQuake_Engine_3.16.3-1.exe . + Download OpenQuake Engine for Windows: https://downloads.openquake.org/pkgs/windows/oq-engine/OpenQuake_Engine_3.16.4-1.exe . Then follow the wizard on screen. > **Warning**: > Administrator level access may be required. diff --git a/doc/sphinx/openquake.calculators.rst b/doc/sphinx/openquake.calculators.rst index cb184265a2ae..100b19091f29 100644 --- a/doc/sphinx/openquake.calculators.rst +++ b/doc/sphinx/openquake.calculators.rst @@ -104,14 +104,6 @@ reportwriter module :undoc-members: :show-inheritance: -conditional_spectrum module --------------------------------------------- - -.. automodule:: openquake.calculators.conditional_spectrum - :members: - :undoc-members: - :show-inheritance: - views module ---------------------------------- diff --git a/openquake/calculators/base.py b/openquake/calculators/base.py index 63a8691f3802..32e2d065b87a 100644 --- a/openquake/calculators/base.py +++ b/openquake/calculators/base.py @@ -997,7 +997,9 @@ def post_process(self): """ oq = self.oqparam if oq.postproc_func: - func = getattr(postproc, oq.postproc_func).main + modname, funcname = oq.postproc_func.rsplit('.', 1) + mod = getattr(postproc, modname) + func = getattr(mod, funcname) if 'csm' in inspect.getargspec(func).args: if hasattr(self, 'csm'): # already there csm = self.csm diff --git a/openquake/calculators/conditional_spectrum.py b/openquake/calculators/conditional_spectrum.py deleted file mode 100644 index a29d547421ee..000000000000 --- a/openquake/calculators/conditional_spectrum.py +++ /dev/null @@ -1,195 +0,0 @@ -# -*- coding: utf-8 -*- -# vim: tabstop=4 shiftwidth=4 softtabstop=4 -# -# Copyright (C) 2021 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 . - -""" -Conditional spectrum calculator, inspired by the disaggregation calculator -""" -import logging -import numpy - -from openquake.baselib import parallel, hdf5 -from openquake.baselib.python3compat import decode -from openquake.hazardlib.probability_map import compute_hazard_maps -from openquake.hazardlib.imt import from_string -from openquake.hazardlib import valid -from openquake.hazardlib.contexts import read_cmakers, read_ctx_by_grp -from openquake.hazardlib.calc.cond_spectra import get_cs_out, outdict -from openquake.calculators import base - -U16 = numpy.uint16 -U32 = numpy.uint32 -TWO24 = 2 ** 24 - -# helper function to be used when saving the spectra as an array -def to_spectra(outdic, n, p): - """ - :param outdic: dictionary rlz_id -> array MNOP - :param n: site index in the range 0..N-1 - :param p: IMLs index in the range 0..P-1 - :returns: conditional spectra as an array of shape (R, M, 2) - """ - R = len(outdic) - M = len(outdic[0]) - out = numpy.zeros((R, M, 2)) - for r in range(R): - c = outdic[r] - out[r, :, 0] = numpy.exp(c[:, n, 1, p]) # NB: seems wrong! - out[r, :, 1] = numpy.sqrt(c[:, n, 2, p]) - return out - - -def store_spectra(dstore, name, R, oq, spectra): - """ - Store the conditional spectra array - """ - N = len(spectra) - stats = list(oq.hazard_stats()) - kind = 'rlz' if 'rlzs' in name else 'stat' - dic = dict(shape_descr=['site_id', 'poe', kind, 'period']) - dic['site_id'] = numpy.arange(N) - dic[kind] = numpy.arange(R) if kind == 'rlz' else stats - dic['poe'] = list(oq.poes) - dic['period'] = [from_string(imt).period for imt in oq.imtls] - hdf5.ArrayWrapper(spectra, dic, ['mea', 'std']).save(name, dstore.hdf5) - # NB: the following would segfault - # dstore.hdf5[name] = hdf5.ArrayWrapper(spectra, dic, ['mea', 'std']) - - -@base.calculators.add('conditional_spectrum') -class ConditionalSpectrumCalculator(base.HazardCalculator): - """ - Conditional spectrum calculator, to be used for few sites only - """ - precalc = 'classical' - accept_precalc = ['classical', 'disaggregation'] - - def execute(self): - """ - Compute the conditional spectrum in a sequential way. - NB: since conditional spectra are meant to be computed only for few - sites, there is no point in parallelizing: the computation is dominated - by the time spent reading the contexts, not by the CPU. - """ - assert self.N <= self.oqparam.max_sites_disagg, ( - self.N, self.oqparam.max_sites_disagg) - logging.warning('Conditional spectrum calculations are still ' - 'experimental') - - oq = self.oqparam - self.full_lt = self.datastore['full_lt'].init() - trt_smrs = self.datastore['trt_smrs'][:] - self.trt_rlzs = self.full_lt.get_trt_rlzs(trt_smrs) - self.trts = list(self.full_lt.gsim_lt.values) - self.imts = list(oq.imtls) - imti = self.imts.index(oq.imt_ref) - self.M = len(self.imts) - dstore = (self.datastore.parent if self.datastore.parent - else self.datastore) - totrups = len(dstore['rup/mag']) - logging.info('Reading {:_d} ruptures'.format(totrups)) - rdt = [('grp_id', U32), ('idx', U32)] - rdata = numpy.zeros(totrups, rdt) - rdata['idx'] = numpy.arange(totrups) - rdata['grp_id'] = dstore['rup/grp_id'][:] - self.periods = [from_string(imt).period for imt in self.imts] - - # extract imls from the "mean" hazard map - if 'hcurves-stats' in dstore: # shape (N, S, M, L1) - curves = dstore.sel('hcurves-stats', stat='mean', imt=oq.imt_ref) - else: # there is only 1 realization - curves = dstore.sel('hcurves-rlzs', rlz_id=0, imt=oq.imt_ref) - self.imls = compute_hazard_maps( # shape (N, L) => (N, P) - curves[:, 0, 0, :], oq.imtls[oq.imt_ref], oq.poes) - N, self.P = self.imls.shape - self.cmakers = read_cmakers(self.datastore) - - # IMPORTANT!! we rely on the fact that the classical part - # of the calculation stores the ruptures in chunks of constant - # grp_id, therefore it is possible to build (start, stop) slices - - # Computing CS - toms = decode(dstore['toms'][:]) - ctx_by_grp = read_ctx_by_grp(dstore) - self.datastore.swmr_on() - smap = parallel.Starmap(get_cs_out, h5=self.datastore) - for grp_id, ctx in ctx_by_grp.items(): - tom = valid.occurrence_model(toms[grp_id]) - cmaker = self.cmakers[grp_id] - smap.submit((cmaker, ctx, imti, self.imls, tom)) - out = smap.reduce() - - # Apply weights and get two dictionaries with integer keys - # (corresponding to the rlz ID) and array values - # of shape (M, N, 3, P) where: - # M is the number of IMTs - # N is the number of sites - # 3 is the number of statistical momenta - # P is the the number of IMLs - outdic, outmean = self._apply_weights(out) - - # sanity check: the mean conditional spectrum must be close to the imls - ref_idx = list(oq.imtls).index(oq.imt_ref) - mean_cs = numpy.exp(outmean[0][ref_idx, :, 1, :]) # shape (N, P) - for n in range(self.N): - for p in range(self.P): - diff = abs(1. - self.imls[n, p]/ mean_cs[n, p]) - if diff > .03: - logging.warning('Conditional Spectrum vs mean IML\nfor ' - 'site_id=%d, poe=%s, imt=%s: %.5f vs %.5f', - n, oq.poes[p], oq.imt_ref, - mean_cs[n, p], self.imls[n, p]) - - # Computing standard deviation - smap = parallel.Starmap(get_cs_out, h5=self.datastore.hdf5) - for grp_id, ctx in ctx_by_grp.items(): - cmaker = self.cmakers[grp_id] - smap.submit((cmaker, ctx, imti, self.imls, tom, outmean[0])) - for res in smap: - for g in res: - out[g][:, :, 2] += res[g][:, :, 2] # STDDEV - return out - - def convert_and_save(self, dsetname, outdic): - """ - Save the conditional spectra - """ - spe = numpy.array([[to_spectra(outdic, n, p) for p in range(self.P)] - for n in range(self.N)]) - store_spectra(self.datastore, dsetname, self.R, self.oqparam, spe) - - def post_execute(self, acc): - # apply weights - outdic, outmean = self._apply_weights(acc) - # convert dictionaries into spectra and save them - self.convert_and_save('cs-rlzs', outdic) - self.convert_and_save('cs-stats', outmean) - - def _apply_weights(self, acc): - # build conditional spectra for each realization - outdic = outdict(self.M, self.N, self.P, 0, self.R) - for g, trs in enumerate(self.trt_rlzs): - for r in trs % TWO24: - outdic[r] += acc[g] - - # build final conditional mean and std - weights = self.datastore['weights'][:] - outmean = outdict(self.M, self.N, self.P, 0, 1) # dict with key 0 - for r, weight in enumerate(weights): - outmean[0] += outdic[r] * weight - - return outdic, outmean diff --git a/openquake/calculators/postproc/compute_mrd.py b/openquake/calculators/postproc/compute_mrd.py index c0b4229d0be8..c1193903a5dd 100644 --- a/openquake/calculators/postproc/compute_mrd.py +++ b/openquake/calculators/postproc/compute_mrd.py @@ -33,7 +33,7 @@ def __init__(self, ctx, cmaker, N): def compute_mrd(inp, crosscorr, imt1, imt2, - meabins, sigbins, monitor): + meabins, sigbins, method, monitor): """ :param inp: an Input object (contexts, contextm maker, num_sites) :param N: the total number of sites @@ -43,11 +43,12 @@ def compute_mrd(inp, crosscorr, imt1, imt2, :param str imt1: the second Intensity Measure Type :param meabins: bins for the means :param sigbins: bins for the sigmas + :param method: string 'direct' or 'indirect' :returns: 4D-matrix with shape (L1, L1, N, G) """ G = len(inp.cmaker.gsims) mrd = calc_mean_rate_dist(inp.ctx, inp.N, inp.cmaker, crosscorr, - imt1, imt2, meabins, sigbins) + imt1, imt2, meabins, sigbins, method) return {g: mrd[:, :, :, i % G] for i, g in enumerate(inp.cmaker.gid)} @@ -62,7 +63,8 @@ def combine_mrds(acc, g_weights): return out -def main(dstore, imt1, imt2, cross_correlation, seed, meabins, sigbins): +def main(dstore, imt1, imt2, cross_correlation, seed, meabins, sigbins, + method='indirect'): """ :param dstore: datastore with the classical calculation @@ -76,7 +78,7 @@ def main(dstore, imt1, imt2, cross_correlation, seed, meabins, sigbins): if L1 > 24: logging.warning('There are many intensity levels (%d), the ' 'calculation can be pretty slow', L1 + 1) - assert N <= 10, 'Too many sites: %d' % N + assert N <= oq.max_sites_disagg, 'Too many sites: %d' % N cmakers = contexts.read_cmakers(dstore) ctx_by_grp = contexts.read_ctx_by_grp(dstore) n = sum(len(ctx) for ctx in ctx_by_grp.values()) @@ -86,7 +88,7 @@ def main(dstore, imt1, imt2, cross_correlation, seed, meabins, sigbins): # NB: a trivial splitting of the contexts would case task-dependency! cmaker = cmakers[grp_id] smap.submit((Input(ctx, cmaker, N), crosscorr, imt1, imt2, - meabins, sigbins)) + meabins, sigbins, method)) acc = smap.reduce() mrd = dstore.create_dset('mrd', float, (L1, L1, N)) trt_rlzs = full_lt.get_trt_rlzs(dstore['trt_smrs'][:]) diff --git a/openquake/calculators/postproc/conditional_spectrum.py b/openquake/calculators/postproc/conditional_spectrum.py new file mode 100644 index 000000000000..91fc5f461238 --- /dev/null +++ b/openquake/calculators/postproc/conditional_spectrum.py @@ -0,0 +1,196 @@ +# -*- coding: utf-8 -*- +# vim: tabstop=4 shiftwidth=4 softtabstop=4 +# +# Copyright (C) 2023, 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 . +""" +Conditional spectrum post-processor +""" +import logging +import numpy +from openquake.baselib import parallel, hdf5, sap +from openquake.baselib.python3compat import decode +from openquake.hazardlib.probability_map import compute_hazard_maps +from openquake.hazardlib.imt import from_string +from openquake.hazardlib import valid, InvalidFile +from openquake.hazardlib.contexts import read_cmakers, read_ctx_by_grp +from openquake.hazardlib.calc.cond_spectra import get_cs_out, outdict + +U16 = numpy.uint16 +U32 = numpy.uint32 +TWO24 = 2 ** 24 + +# helper function to be used when saving the spectra as an array +def to_spectra(outdic, n, p): + """ + :param outdic: dictionary rlz_id -> array MNOP + :param n: site index in the range 0..N-1 + :param p: IMLs index in the range 0..P-1 + :returns: conditional spectra as an array of shape (R, M, 2) + """ + R = len(outdic) + M = len(outdic[0]) + out = numpy.zeros((R, M, 2)) + for r in range(R): + c = outdic[r] + out[r, :, 0] = numpy.exp(c[:, n, 1, p]) # NB: seems wrong! + out[r, :, 1] = numpy.sqrt(c[:, n, 2, p]) + return out + + +def store_spectra(dstore, name, R, oq, spectra): + """ + Store the conditional spectra array + """ + N = len(spectra) + stats = list(oq.hazard_stats()) + kind = 'rlz' if 'rlzs' in name else 'stat' + dic = dict(shape_descr=['site_id', 'poe', kind, 'period']) + dic['site_id'] = numpy.arange(N) + dic[kind] = numpy.arange(R) if kind == 'rlz' else stats + dic['poe'] = list(oq.poes) + dic['period'] = [from_string(imt).period for imt in oq.imtls] + hdf5.ArrayWrapper(spectra, dic, ['mea', 'std']).save(name, dstore.hdf5) + # NB: the following would segfault + # dstore.hdf5[name] = hdf5.ArrayWrapper(spectra, dic, ['mea', 'std']) + + +def compute_cs(dstore, oq, N, M, P): + """ + Compute the conditional spectrum in a sequential way. + NB: since conditional spectra are meant to be computed only for few + sites, there is no point in parallelizing: the computation is dominated + by the time spent reading the contexts, not by the CPU. + """ + full_lt = dstore['full_lt'].init() + trt_rlzs = full_lt.get_trt_rlzs(dstore['trt_smrs'][:]) + R = full_lt.get_num_paths() + imts = list(oq.imtls) + imti = imts.index(oq.imt_ref) + totrups = len(dstore['rup/mag']) + logging.info('Reading {:_d} contexts'.format(totrups)) + rdt = [('grp_id', U32), ('idx', U32)] + rdata = numpy.zeros(totrups, rdt) + rdata['idx'] = numpy.arange(totrups) + rdata['grp_id'] = dstore['rup/grp_id'][:] + + # extract imls from the "mean" hazard map + if 'hcurves-stats' in dstore: # shape (N, S, M, L1) + curves = dstore.sel('hcurves-stats', stat='mean', imt=oq.imt_ref) + else: # there is only 1 realization + curves = dstore.sel('hcurves-rlzs', rlz_id=0, imt=oq.imt_ref) + imls = compute_hazard_maps( # shape (N, L) => (N, P) + curves[:, 0, 0, :], oq.imtls[oq.imt_ref], oq.poes) + N, P = imls.shape + cmakers = read_cmakers(dstore) + + # Computing CS + toms = decode(dstore['toms'][:]) + ctx_by_grp = read_ctx_by_grp(dstore) + dstore.swmr_on() + smap = parallel.Starmap(get_cs_out, h5=dstore) + for grp_id, ctx in ctx_by_grp.items(): + tom = valid.occurrence_model(toms[grp_id]) + cmaker = cmakers[grp_id] + smap.submit((cmaker, ctx, imti, imls, tom)) + out = smap.reduce() + + # Apply weights and get two dictionaries with integer keys + # (corresponding to the rlz ID) and array values + # of shape (M, N, 3, P) where: + # M is the number of IMTs + # N is the number of sites + # 3 is the number of statistical momenta + # P is the the number of IMLs + outdic, outmean = _apply_weights(dstore, out, N, R, M, P, trt_rlzs) + + # sanity check: the mean conditional spectrum must be close to the imls + ref_idx = list(oq.imtls).index(oq.imt_ref) + mean_cs = numpy.exp(outmean[0][ref_idx, :, 1, :]) # shape (N, P) + for n in range(N): + for p in range(P): + diff = abs(1. - imls[n, p]/ mean_cs[n, p]) + if diff > .03: + logging.warning('Conditional Spectrum vs mean IML\nfor ' + 'site_id=%d, poe=%s, imt=%s: %.5f vs %.5f', + n, oq.poes[p], oq.imt_ref, + mean_cs[n, p], imls[n, p]) + + # Computing standard deviation + dstore.swmr_on() + smap = parallel.Starmap(get_cs_out, h5=dstore.hdf5) + for grp_id, ctx in ctx_by_grp.items(): + cmaker = cmakers[grp_id] + smap.submit((cmaker, ctx, imti, imls, tom, outmean[0])) + for res in smap: + for g in res: + out[g][:, :, 2] += res[g][:, :, 2] # STDDEV + + # apply weights + outdic, outmean = _apply_weights(dstore, out, N, R, M, P, trt_rlzs) + # convert dictionaries into spectra and save them + convert_and_save(dstore, 'cs-rlzs', outdic, oq, N, R, P) + convert_and_save(dstore, 'cs-stats', outmean, oq, N, R, P) + + +def convert_and_save(dstore, dsetname, outdic, oq, N, R, P): + """ + Save the conditional spectra + """ + spe = numpy.array([[to_spectra(outdic, n, p) for p in range(P)] + for n in range(N)]) + store_spectra(dstore, dsetname, R, oq, spe) + + +def _apply_weights(dstore, acc, N, R, M, P, trt_rlzs): + # build conditional spectra for each realization + outdic = outdict(M, N, P, 0, R) + for g, trs in enumerate(trt_rlzs): + for r in trs % TWO24: + outdic[r] += acc[g] + + # build final conditional mean and std + weights = dstore['weights'][:] + outmean = outdict(M, N, P, 0, 1) # dict with key 0 + for r, weight in enumerate(weights): + outmean[0] += outdic[r] * weight + + return outdic, outmean + + +def main(dstore): + logging.warning('Conditional spectrum calculations are still ' + 'experimental') + oq = dstore['oqparam'] + if not oq.cross_correl: + raise InvalidFile("%(job_ini)s: you must specify cross_correlation" + % oq.inputs) + if not oq.imt_ref: + raise InvalidFile("%(job_ini)s: you must specify imt_ref" % oq.inputs) + if not oq.poes: + raise InvalidFile("%(job_ini)s: you must specify the poes" % oq.inputs) + if list(oq.hazard_stats()) != ['mean']: + raise InvalidFile('%(job_ini)s: only the mean is supported' % oq.inputs) + + sitecol = dstore['sitecol'] + N = len(sitecol) + M = len(oq.imtls) + P = len(oq.poes) + assert N <= oq.max_sites_disagg, (N, oq.max_sites_disagg) + compute_cs(dstore, oq, N, M, P) + + +if __name__ == '__main__': + sap.run(main) diff --git a/openquake/calculators/tests/logictree_test.py b/openquake/calculators/tests/logictree_test.py index 91c26af22fc2..597f60002362 100644 --- a/openquake/calculators/tests/logictree_test.py +++ b/openquake/calculators/tests/logictree_test.py @@ -112,7 +112,7 @@ def test_case_04(self): calculation_mode='preclassical') hc_id = str(self.calc.datastore.calc_id) self.run_calc(case_04.__file__, 'job.ini', hazard_calculation_id=hc_id, - postproc_func='disagg_by_rel_sources') + postproc_func='disagg_by_rel_sources.main') [fname] = export(('hcurves', 'csv'), self.calc.datastore) self.assertEqualFiles('expected/curve-mean.csv', fname) diff --git a/openquake/calculators/tests/postproc_test.py b/openquake/calculators/tests/postproc_test.py index 34d0f438cba0..15e0709b099a 100644 --- a/openquake/calculators/tests/postproc_test.py +++ b/openquake/calculators/tests/postproc_test.py @@ -34,7 +34,7 @@ class PostProcTestCase(CalculatorTestCase): def test_mrd(self): # Computes the mean rate density using a simple PSHA input model - self.run_calc(case_mrd.__file__, 'job.ini', postproc_func='dummy') + self.run_calc(case_mrd.__file__, 'job.ini', postproc_func='dummy.main') fnames = export(('hcurves', 'csv'), self.calc.datastore) for fname in fnames: self.assertEqualFiles('expected/%s' % strip_calc_id(fname), fname) diff --git a/openquake/commonlib/oqvalidation.py b/openquake/commonlib/oqvalidation.py index 23f43f52a2d1..cd2f97f760f0 100644 --- a/openquake/commonlib/oqvalidation.py +++ b/openquake/commonlib/oqvalidation.py @@ -364,10 +364,10 @@ Default: no default imt_ref: - Reference intensity measure type usedto compute the conditional spectrum. + Reference intensity measure type used to compute the conditional spectrum. The imt_ref must belong to the list of IMTs of the calculation. Example: *imt_ref = SA(0.15)*. - Default: no default + Default: empty string individual_rlzs: When set, store the individual hazard curves and/or individual risk curves @@ -566,7 +566,7 @@ postproc_func: Specify a postprocessing function in calculators/postproc. - Example: *postproc_func = compute_mrd* + Example: *postproc_func = compute_mrd.main* Default: '' (no postprocessing) postproc_args: @@ -837,7 +837,6 @@ 'multi_risk', 'classical_bcr', 'preclassical', - 'conditional_spectrum', 'event_based_damage', 'scenario_damage'] @@ -966,7 +965,7 @@ class OqParam(valid.ParamSet): ignore_missing_costs = valid.Param(valid.namelist, []) ignore_covs = valid.Param(valid.boolean, False) iml_disagg = valid.Param(valid.floatdict, {}) # IMT -> IML - imt_ref = valid.Param(valid.intensity_measure_type) + imt_ref = valid.Param(valid.intensity_measure_type, '') individual_rlzs = valid.Param(valid.boolean, None) inputs = valid.Param(dict, {}) ash_wet_amplification_factor = valid.Param(valid.positivefloat, 1.0) @@ -1005,7 +1004,7 @@ class OqParam(valid.ParamSet): poes = valid.Param(valid.probabilities, []) poes_disagg = valid.Param(valid.probabilities, []) pointsource_distance = valid.Param(valid.floatdict, {'default': PSDIST}) - postproc_func = valid.Param(valid.simple_id, '') + postproc_func = valid.Param(valid.mod_func, 'dummy.main') postproc_args = valid.Param(valid.dictionary, {}) prefer_global_site_params = valid.Param(valid.boolean, None) ps_grid_spacing = valid.Param(valid.positivefloat, 0) @@ -1292,13 +1291,6 @@ def __init__(self, **names_vals): raise InvalidFile('%s: you cannot set rlzs_index and ' 'num_rlzs_disagg at the same time' % job_ini) - # checks for conditional_spectrum - if self.calculation_mode == 'conditional_spectrum': - if not self.poes: - raise InvalidFile("%s: you must specify the poes" % job_ini) - elif list(self.hazard_stats()) != ['mean']: - raise InvalidFile('%s: only the mean is supported' % job_ini) - # checks for classical_damage if self.calculation_mode == 'classical_damage': if self.conditional_loss_poes: @@ -2062,7 +2054,9 @@ def __toh5__(self): return hdf5.dumps(vars(self)), {} def __fromh5__(self, array, attrs): - if isinstance(array, numpy.ndarray): # old format <= 3.11 + if isinstance(array, numpy.ndarray): + # old format <= 3.11, tested in read_old_data, + # used to read old GMFs dd = collections.defaultdict(dict) for (name_, literal_) in array: name = python3compat.decode(name_) @@ -2073,5 +2067,11 @@ def __fromh5__(self, array, attrs): else: dd[name] = ast.literal_eval(literal) vars(self).update(dd) - else: # new format >= 3.12 + else: + # for version >= 3.12 vars(self).update(json.loads(python3compat.decode(array))) + + Idist = calc.filters.IntegrationDistance + if hasattr(self, 'maximum_distance') and not isinstance( + self.maximum_distance, Idist): + self.maximum_distance = Idist(**self.maximum_distance) diff --git a/openquake/commonlib/readinput.py b/openquake/commonlib/readinput.py index 0cc8c99f1182..5686c1d1a364 100644 --- a/openquake/commonlib/readinput.py +++ b/openquake/commonlib/readinput.py @@ -236,15 +236,15 @@ def update(params, items, base_path): params[key] = value -def _warn_about_duplicates(cp): +def check_params(cp, fname): params_sets = [ set(cp.options(section)) for section in cp.sections()] for pair in itertools.combinations(params_sets, 2): params_intersection = sorted(set.intersection(*pair)) if params_intersection: - logging.warning( - f'Parameter(s) {params_intersection} is(are) defined in' - f' multiple sections') + raise InvalidFile( + f'{fname}: parameter(s) {params_intersection} is(are) defined' + ' in multiple sections') # NB: this function must NOT log, since it is called when the logging @@ -286,7 +286,7 @@ def get_params(job_ini, kw={}): params = dict(base_path=base_path, inputs={'job_ini': job_ini}) cp = configparser.ConfigParser(interpolation=None) cp.read([job_ini], encoding='utf-8-sig') # skip BOM on Windows - _warn_about_duplicates(cp) + check_params(cp, job_ini) dic = {} for sect in cp.sections(): dic.update(cp.items(sect)) diff --git a/openquake/commonlib/tests/readinput_test.py b/openquake/commonlib/tests/readinput_test.py index 0059caf693ec..19060902c7a2 100644 --- a/openquake/commonlib/tests/readinput_test.py +++ b/openquake/commonlib/tests/readinput_test.py @@ -160,19 +160,10 @@ def test_duplicated_parameter(self): bar = foo aggregate_by = taxonomy, policy """) - with self.assertLogs() as captured: + with self.assertRaises(InvalidFile) as ctx: readinput.get_params(job_config, {}) - warning_general_bar = ( - "Parameter(s) ['aggregate_by', 'foo'] is(are) defined in" - " multiple sections") - warning_foo_bar = ("Parameter(s) ['bar'] is(are) defined in" - " multiple sections") - self.assertEqual(captured.records[0].levelname, 'WARNING') - self.assertEqual( - captured.records[0].getMessage(), warning_general_bar) - self.assertEqual(captured.records[1].levelname, 'WARNING') - self.assertEqual( - captured.records[1].getMessage(), warning_foo_bar) + msg = "['aggregate_by', 'foo'] is(are) defined in multiple sections" + self.assertIn(msg, str(ctx.exception)) def sitemodel(): diff --git a/openquake/engine/aelo.py b/openquake/engine/aelo.py index 6b42a9696689..aea9bfc0fad8 100644 --- a/openquake/engine/aelo.py +++ b/openquake/engine/aelo.py @@ -53,7 +53,7 @@ def get_params_from(inputs): params['sites'] = '%(lon)s %(lat)s' % inputs if 'vs30' in inputs: params['override_vs30'] = '%(vs30)s' % inputs - params['postproc_func'] = 'disagg_by_rel_sources' + params['postproc_func'] = 'disagg_by_rel_sources.main' # params['cachedir'] = datastore.get_datadir() return params diff --git a/openquake/hazardlib/calc/mrd.py b/openquake/hazardlib/calc/mrd.py index abea54212c7f..f5fa2a175c82 100644 --- a/openquake/hazardlib/calc/mrd.py +++ b/openquake/hazardlib/calc/mrd.py @@ -42,7 +42,7 @@ def get_uneven_bins_edges(lefts, num_bins): return numpy.array(tmp) -def update_mrd(ctxt: numpy.recarray, cm, crosscorr, mrd): +def update_mrd(ctxt: numpy.recarray, cm, crosscorr, mrd, monitor=Monitor()): """ This computes the mean rate density by means of the multivariate normal function available in scipy. @@ -85,7 +85,8 @@ def update_mrd(ctxt: numpy.recarray, cm, crosscorr, mrd): # Compute the MRD for the current rupture mvn = sts.multivariate_normal(mea[slc0], comtx) - partial = get_mrd(mvn, grids) + with monitor: # this is ultra-slow! + partial = get_mrd(mvn, grids) # Check msg = f'{numpy.max(partial):.8f}' @@ -210,10 +211,10 @@ def update_mrd_indirect(ctx, cm, corrm, be_mea, be_sig, mrd, monitor=Monitor()): mrd[:, :, gid] += arr[R] * partial -def calc_mean_rate_dist(ctxt, nsites, cmaker, crosscorr, imt1, imt2, - bins_mea, bins_sig, mon=Monitor()): +def calc_mean_rate_dist(ctx, nsites, cmaker, crosscorr, imt1, imt2, + bins_mea, bins_sig, method='indirect', mon=Monitor()): """ - :param ctxt: a sequence of parametric sources + :param ctx: a sequence of parametric sources :param num_sites: total number of sites (small) :param cmaker: a ContextMaker instance :param crosscorr: a CrossCorrelation instance @@ -221,6 +222,7 @@ def calc_mean_rate_dist(ctxt, nsites, cmaker, crosscorr, imt1, imt2, :param str imt2: second IMT to consider (must be inside cmaker.imtls) :param bins_mea: bins for the mean :param bins_sig: bins for the standard deviation + :param method: a string 'direct' or 'indirect' """ cm = cmaker.restrict([imt1, imt2]) G = len(cm.gsims) @@ -228,7 +230,11 @@ def calc_mean_rate_dist(ctxt, nsites, cmaker, crosscorr, imt1, imt2, corrm = crosscorr.get_cross_correlation_mtx(cm.imts) mrd = numpy.zeros((len1, len1, nsites, G)) for sid in range(nsites): - update_mrd_indirect( - ctxt[ctxt.sids == sid], cm, corrm, - bins_mea, bins_sig, mrd[:, :, sid], mon) + if method == 'direct': + update_mrd( + ctx[ctx.sids == sid], cm, crosscorr, mrd[:, :, sid], mon) + else: + update_mrd_indirect( + ctx[ctx.sids == sid], cm, corrm, + bins_mea, bins_sig, mrd[:, :, sid], mon) return mrd diff --git a/openquake/hazardlib/gsim/lanzano_luzi_2019.py b/openquake/hazardlib/gsim/lanzano_luzi_2019.py index 41b1881331ca..861110c9514e 100644 --- a/openquake/hazardlib/gsim/lanzano_luzi_2019.py +++ b/openquake/hazardlib/gsim/lanzano_luzi_2019.py @@ -141,7 +141,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi): _compute_distance(ctx.rhypo, C, self.kind) + _get_site_amplification(ctx, C)) - istddevs = [C['sigma'], C['phi'], C['tau']] + istddevs = [C['sigma'], C['tau'], C['phi']] # Convert units to g, but only for PGA and SA (not PGV) if imt.string.startswith(("SA", "PGA")): diff --git a/openquake/hazardlib/tests/calc/data/job.ini b/openquake/hazardlib/tests/calc/data/job.ini index 4ce21d62a6ed..75733faafadb 100644 --- a/openquake/hazardlib/tests/calc/data/job.ini +++ b/openquake/hazardlib/tests/calc/data/job.ini @@ -1,7 +1,8 @@ [general] description = Conditional spectrum -calculation_mode = conditional_spectrum +calculation_mode = classical +postproc_func = conditional_spectrum.main uniform_hazard_spectra = true [geometry] @@ -46,6 +47,7 @@ intensity_measure_types_and_levels = { truncation_level = 3 maximum_distance = 200.0 +# conditional_spectrum parameters cross_correlation = BakerJayaram2008 imt_ref = SA(0.2) poes = 0.002105 diff --git a/openquake/hazardlib/valid.py b/openquake/hazardlib/valid.py index 301c96350aae..861507eff2d2 100644 --- a/openquake/hazardlib/valid.py +++ b/openquake/hazardlib/valid.py @@ -353,6 +353,7 @@ def __call__(self, value): source_id = SimpleId(MAX_ID_LENGTH, r'^[\w\-_:]+$') nice_string = SimpleId( # nice for Windows, Linux, HDF5 and XML ASSET_ID_LENGTH, r'[a-zA-Z0-9\.`!#$%\(\)\+/,;@\[\]\^_{|}~-]+') +mod_func = SimpleId(MAX_ID_LENGTH, r'[\w_]+\.[\w_]+') def risk_id(value): diff --git a/openquake/qa_tests_data/classical_risk/case_master/job.ini b/openquake/qa_tests_data/classical_risk/case_master/job.ini index 7fa316cf889e..685dfa5a6c57 100644 --- a/openquake/qa_tests_data/classical_risk/case_master/job.ini +++ b/openquake/qa_tests_data/classical_risk/case_master/job.ini @@ -56,7 +56,6 @@ risk_investigation_time = 50 lrem_steps_per_interval = 1 [risk_outputs] -quantiles = 0.15, 0.50, 0.85 conditional_loss_poes = 0.02, 0.10 [export] diff --git a/openquake/qa_tests_data/conditional_spectrum/case_1/job.ini b/openquake/qa_tests_data/conditional_spectrum/case_1/job.ini index b96b4988565f..d7b35c434f7a 100644 --- a/openquake/qa_tests_data/conditional_spectrum/case_1/job.ini +++ b/openquake/qa_tests_data/conditional_spectrum/case_1/job.ini @@ -1,7 +1,8 @@ [general] description = Conditional spectrum -calculation_mode = conditional_spectrum +calculation_mode = classical +postproc_func = conditional_spectrum.main uniform_hazard_spectra = true [geometry] diff --git a/openquake/qa_tests_data/conditional_spectrum/case_2/job.ini b/openquake/qa_tests_data/conditional_spectrum/case_2/job.ini index 8fc0bc97d6d4..58a26c5c78b5 100644 --- a/openquake/qa_tests_data/conditional_spectrum/case_2/job.ini +++ b/openquake/qa_tests_data/conditional_spectrum/case_2/job.ini @@ -1,7 +1,8 @@ [general] description = Conditional spectrum -calculation_mode = conditional_spectrum +calculation_mode = classical +postproc_func = conditional_spectrum.main uniform_hazard_spectra = true [geometry] diff --git a/openquake/qa_tests_data/event_based/case_29/job.ini b/openquake/qa_tests_data/event_based/case_29/job.ini index 7abe46aa76e5..1fc469b83865 100644 --- a/openquake/qa_tests_data/event_based/case_29/job.ini +++ b/openquake/qa_tests_data/event_based/case_29/job.ini @@ -34,7 +34,6 @@ intensity_measure_types_and_levels = {"PGA": logscale(0.005, 3.00, 20)} truncation_level = 5 maximum_distance = {'default': 100.} minimum_magnitude = 4.0 -ses_per_logic_tree_path = 100 [output] diff --git a/openquake/qa_tests_data/event_based_risk/case_1f/job.ini b/openquake/qa_tests_data/event_based_risk/case_1f/job.ini index b2a0cf629f51..2a3cd55f380f 100644 --- a/openquake/qa_tests_data/event_based_risk/case_1f/job.ini +++ b/openquake/qa_tests_data/event_based_risk/case_1f/job.ini @@ -15,9 +15,6 @@ structural_vulnerability_file = vulnerability_model.xml [hazard] asset_hazard_distance = 20.0 -[sites] -exposure_file = exposure_model.xml - [site_params] reference_vs30_type = measured reference_vs30_value = 760.0 diff --git a/openquake/qa_tests_data/logictree/case_01/job.ini b/openquake/qa_tests_data/logictree/case_01/job.ini index ab06f67aea1e..a012a2770117 100644 --- a/openquake/qa_tests_data/logictree/case_01/job.ini +++ b/openquake/qa_tests_data/logictree/case_01/job.ini @@ -5,7 +5,7 @@ calculation_mode = classical random_seed = 1066 disagg_by_src = true use_rates = true -postproc_func = disagg_by_rel_sources +postproc_func = disagg_by_rel_sources.main [geometry] diff --git a/openquake/qa_tests_data/logictree/case_05/job.ini b/openquake/qa_tests_data/logictree/case_05/job.ini index f47fb5cf4b75..9381970ebaef 100644 --- a/openquake/qa_tests_data/logictree/case_05/job.ini +++ b/openquake/qa_tests_data/logictree/case_05/job.ini @@ -5,7 +5,7 @@ calculation_mode = classical random_seed = 23 disagg_by_src = true use_rates = true -postproc_func = disagg_by_rel_sources +postproc_func = disagg_by_rel_sources.main [geometry] diff --git a/openquake/qa_tests_data/logictree/case_05/job_bis.ini b/openquake/qa_tests_data/logictree/case_05/job_bis.ini index 9c6c859dac1b..701edffb1afd 100644 --- a/openquake/qa_tests_data/logictree/case_05/job_bis.ini +++ b/openquake/qa_tests_data/logictree/case_05/job_bis.ini @@ -3,7 +3,7 @@ description = Reducing source specific logic tree, sampling calculation_mode = classical random_seed = 23 -postproc_func = disagg_by_rel_sources +postproc_func = disagg_by_rel_sources.main [geometry] diff --git a/openquake/qa_tests_data/logictree/case_06/job.ini b/openquake/qa_tests_data/logictree/case_06/job.ini index 8afddd4dafa8..83e04e962203 100644 --- a/openquake/qa_tests_data/logictree/case_06/job.ini +++ b/openquake/qa_tests_data/logictree/case_06/job.ini @@ -5,7 +5,7 @@ calculation_mode = classical random_seed = 1066 disagg_by_src = true use_rates = true -postproc_func = disagg_by_rel_sources +postproc_func = disagg_by_rel_sources.main [geometry] diff --git a/openquake/qa_tests_data/logictree/case_07/job.ini b/openquake/qa_tests_data/logictree/case_07/job.ini index 2b57de9801fe..e876c36efd2a 100644 --- a/openquake/qa_tests_data/logictree/case_07/job.ini +++ b/openquake/qa_tests_data/logictree/case_07/job.ini @@ -5,7 +5,7 @@ calculation_mode = classical random_seed = 1066 disagg_by_src = true use_rates = true -postproc_func = disagg_by_rel_sources +postproc_func = disagg_by_rel_sources.main [geometry] diff --git a/openquake/qa_tests_data/logictree/case_07/sampling.ini b/openquake/qa_tests_data/logictree/case_07/sampling.ini index 0f905a0568af..ff2bf1aedb29 100644 --- a/openquake/qa_tests_data/logictree/case_07/sampling.ini +++ b/openquake/qa_tests_data/logictree/case_07/sampling.ini @@ -5,7 +5,7 @@ calculation_mode = classical random_seed = 1066 disagg_by_src = true use_rates = true -postproc_func = disagg_by_rel_sources +postproc_func = disagg_by_rel_sources.main [geometry] diff --git a/openquake/qa_tests_data/logictree/case_12/job.ini b/openquake/qa_tests_data/logictree/case_12/job.ini index 3e20cb51733c..ba5ece00e55e 100755 --- a/openquake/qa_tests_data/logictree/case_12/job.ini +++ b/openquake/qa_tests_data/logictree/case_12/job.ini @@ -3,7 +3,7 @@ description = 2018 PSHA model - North Africa calculation_mode = classical random_seed = 19 -postproc_func = disagg_by_rel_sources +postproc_func = disagg_by_rel_sources.main [geometry] diff --git a/openquake/qa_tests_data/logictree/case_13/job_gmv.ini b/openquake/qa_tests_data/logictree/case_13/job_gmv.ini index 0f55c2cfb159..61a9b03b51e7 100644 --- a/openquake/qa_tests_data/logictree/case_13/job_gmv.ini +++ b/openquake/qa_tests_data/logictree/case_13/job_gmv.ini @@ -34,5 +34,5 @@ intensity_measure_types_and_levels = {"PGA": [0.005, 0.007, 0.0098, 0.0137, 0.01 truncation_level = 3 maximum_distance = 100 -postproc_func = med_gmv +postproc_func = med_gmv.main postproc_args = {} diff --git a/openquake/qa_tests_data/logictree/case_19/job.ini b/openquake/qa_tests_data/logictree/case_19/job.ini index 0c8f763e088d..076792155280 100644 --- a/openquake/qa_tests_data/logictree/case_19/job.ini +++ b/openquake/qa_tests_data/logictree/case_19/job.ini @@ -6,7 +6,7 @@ random_seed = 23 disagg_by_src = true use_rates = true concurrent_tasks = 1 -postproc_func = disagg_by_rel_sources +postproc_func = disagg_by_rel_sources.main [geometry] sites = 16.2657 39.3002 diff --git a/openquake/qa_tests_data/logictree/case_28/job.ini b/openquake/qa_tests_data/logictree/case_28/job.ini index 708ff224f6be..d1b288b74198 100644 --- a/openquake/qa_tests_data/logictree/case_28/job.ini +++ b/openquake/qa_tests_data/logictree/case_28/job.ini @@ -30,10 +30,7 @@ reference_depth_to_1pt0km_per_sec = 100.0 source_model_logic_tree_file = source_model_logic_tree.xml gsim_logic_tree_file = gmpe_logic_tree.xml - investigation_time = 50.0 -poes = 0.1 0.02 - intensity_measure_types_and_levels = { "PGA": [0.005, 0.007, 0.0098, 0.0137, 0.0192, 0.0269, 0.0376, 0.0527, 0.0738, 0.103, 0.145, 0.203, 0.284, 0.397, 0.556, 0.778, 1.09, 1.52, 2.13], "SA(0.05)": [0.005, 0.007, 0.0098, 0.0137, 0.0192, 0.0269, 0.0376, 0.0527, 0.0738, 0.103, 0.145, 0.203, 0.284, 0.397, 0.556, 0.778, 1.09, 1.52, 2.13], diff --git a/openquake/qa_tests_data/postproc/case_mrd/job.ini b/openquake/qa_tests_data/postproc/case_mrd/job.ini index ce0f3c75da81..55d1965fdda6 100644 --- a/openquake/qa_tests_data/postproc/case_mrd/job.ini +++ b/openquake/qa_tests_data/postproc/case_mrd/job.ini @@ -32,14 +32,15 @@ maximum_distance = 300 [postproc] -postproc_func = compute_mrd +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]} + 'sigbins': [0.2, 0.3, 0.4, 0.5, 0.6, 0.7], + 'method': 'indirect'} [output] diff --git a/openquake/qa_tests_data/scenario_risk/case_3/job.ini b/openquake/qa_tests_data/scenario_risk/case_3/job.ini index d38e7a108815..662b3b3d9d20 100644 --- a/openquake/qa_tests_data/scenario_risk/case_3/job.ini +++ b/openquake/qa_tests_data/scenario_risk/case_3/job.ini @@ -8,8 +8,6 @@ structural_vulnerability_file = vulnerability_model.xml region = 78.0 31.5,89.5 31.5,89.5 25.5,78 25.5 -export_dir = /tmp - [geometry] exposure_file = exposure_model.xml