Merge pull request #10128 from gem/aristotle

Adding with_betw_ratio in Aristotle mode/2

debian/changelog

@@ -1,4 +1,5 @@
   [Michele Simionato]
+  * Added parameter `with_betw_ratio`
   * Added command `oq info peril`
   * Extended consequences to perils
   * Replaced taxonomy mapping by loss type with taxonomy mapping by peril

openquake/baselib/general.py

@@ -437,7 +437,7 @@ def gettemp(content=None, dir=None, prefix="tmp", suffix="tmp", remove=True):
     if dir is not None:
         if not os.path.exists(dir):
             os.makedirs(dir)
-    fh, path = tempfile.mkstemp(dir=dir or config.directory.custom_tmp,
+    fh, path = tempfile.mkstemp(dir=dir or config.directory.custom_tmp or None,
                                 prefix=prefix, suffix=suffix)
     if remove:
         _tmp_paths.append(path)

openquake/calculators/event_based.py

@@ -37,7 +37,7 @@
     nofilter, getdefault, get_distances, SourceFilter)
 from openquake.hazardlib.calc.gmf import GmfComputer
 from openquake.hazardlib.calc.conditioned_gmfs import ConditionedGmfComputer
-from openquake.hazardlib import valid, logictree, InvalidFile
+from openquake.hazardlib import logictree, InvalidFile
 from openquake.hazardlib.calc.stochastic import get_rup_array, rupture_dt
 from openquake.hazardlib.source.rupture import (
     RuptureProxy, EBRupture, get_ruptures)
@@ -340,7 +340,6 @@ def starmap_from_rups(func, oq, full_lt, sitecol, dstore, save_tmp=None):
         # assume scenario with a single true rupture
         rlzs_by_gsim = full_lt.get_rlzs_by_gsim(0)
         cmaker = ContextMaker(trt, rlzs_by_gsim, oq)
-        cmaker.gid = numpy.arange(len(rlzs_by_gsim))
         cmaker.scenario = True
         maxdist = oq.maximum_distance(cmaker.trt)
         srcfilter = SourceFilter(sitecol.complete, maxdist)
@@ -372,7 +371,6 @@ def starmap_from_rups(func, oq, full_lt, sitecol, dstore, save_tmp=None):
         extra = sitecol.array.dtype.names
         rlzs_by_gsim = full_lt.get_rlzs_by_gsim(trt_smr)
         cmaker = ContextMaker(trt, rlzs_by_gsim, oq, extraparams=extra)
-        cmaker.gid = numpy.arange(len(rlzs_by_gsim))
         cmaker.min_mag = getdefault(oq.minimum_magnitude, trt)
         for gsim in rlzs_by_gsim:
             toml_gsims.append(gsim._toml)
@@ -430,7 +428,6 @@ def compute_avg_gmf(gmf_df, weights, min_iml):
 
 def read_gsim_lt(oq):
     # in aristotle mode the gsim_lt is read from the exposure.hdf5 file
-    gsim_lt = readinput.get_gsim_lt(oq)
     if oq.aristotle:
         if not oq.mosaic_model:
             if oq.rupture_dict:
@@ -458,13 +455,8 @@ def read_gsim_lt(oq):
         gsim_lt = logictree.GsimLogicTree.from_hdf5(
             expo_hdf5, oq.mosaic_model,
             oq.tectonic_region_type.encode('utf8'))
-        # add with_betw_ratio when only the total stddev is defined
-        betw = {'with_betw_ratio': 1.7}
-        for gsims in gsim_lt.values.values():
-            for g, gsim in enumerate(gsims):
-                if len(gsim.DEFINED_FOR_STANDARD_DEVIATION_TYPES) == 1:
-                    gsims[g] = valid.modified_gsim(
-                        gsim, add_between_within_stds=betw)
+    else:
+        gsim_lt = readinput.get_gsim_lt(oq)
     return gsim_lt
 
 
@@ -649,7 +641,6 @@ def _read_scenario_ruptures(self):
             rup = readinput.get_rupture(oq)
             oq.mags_by_trt = {trt: ['%.2f' % rup.mag]}
             self.cmaker = ContextMaker(trt, rlzs_by_gsim, oq)
-            self.cmaker.gid = numpy.arange(len(rlzs_by_gsim))
             if self.N > oq.max_sites_disagg:  # many sites, split rupture
                 ebrs = []
                 for i in range(ngmfs):

openquake/commands/plot.py

@@ -1063,9 +1063,11 @@ def plot_wkt(wkt_string):
     poly = wkt.loads(wkt_string)
     if hasattr(poly, 'exterior'):
         coo = numpy.array(poly.exterior.coords)
-    else:  # LINESTRING
+    else:  # POINT or LINESTRING
         coo = numpy.array(poly.coords)
-    plt.plot(coo[:, 0], coo[:, 1], '-')
+    _fig, ax = plt.subplots()
+    ax.plot(coo[:, 0], coo[:, 1], 'o')
+    add_borders(ax, readinput.read_mosaic_df, buffer=0.)
     return plt
 
 
@@ -1108,7 +1110,7 @@ def main(what,
     if what.endswith('.csv'):
         plot_csv(what)
         return
-    if what.startswith(('POLYGON', 'LINESTRING')):
+    if what.startswith(('POINT', 'POLYGON', 'LINESTRING')):
         plt = plot_wkt(what)
         plt.show()
         return

openquake/commonlib/oqvalidation.py

@@ -868,6 +868,12 @@
   Used to specify the width of the Magnitude Frequency Distribution.
   Example: *width_of_mfd_bin = 0.2*.
   Default: None
+
+with_betw_ratio:
+  Specify the between ratio for GSIMs with only the Total StdDev.
+  This is necessary in conditioned GMFs calculations.
+  Example: *with_betw_ratio = 1.7*
+  Default: None
 """ % __version__
 
 PSDIST = float(config.performance.pointsource_distance)
@@ -1156,6 +1162,7 @@ class OqParam(valid.ParamSet):
     use_rates = valid.Param(valid.boolean, False)
     vs30_tolerance = valid.Param(int, 0)
     width_of_mfd_bin = valid.Param(valid.positivefloat, None)
+    with_betw_ratio = valid.Param(valid.positivefloat, None)
 
     @property
     def no_pointsource_distance(self):

openquake/commonlib/readinput.py

@@ -765,11 +765,8 @@ def get_gsim_lt(oqparam, trts=('*',)):
             # NB: gsim.DEFINED_FOR_TECTONIC_REGION_TYPE can be != trt,
             # but it is not an error, it is actually the most common case!
             if gmfcorr and (gsim.DEFINED_FOR_STANDARD_DEVIATION_TYPES ==
-                            {StdDev.TOTAL}):
-                modifications = list(gsim.kwargs.keys())
-                if not (type(gsim).__name__ == 'ModifiableGMPE' and
-                        'add_between_within_stds' in modifications):
-                    raise CorrelationButNoInterIntraStdDevs(gmfcorr, gsim)
+                            {StdDev.TOTAL}) and not oqparam.with_betw_ratio:
+                raise CorrelationButNoInterIntraStdDevs(gmfcorr, gsim)
     imt_dep_w = any(len(branch.weight.dic) > 1 for branch in gsim_lt.branches)
     if oqparam.number_of_logic_tree_samples and imt_dep_w:
         logging.error('IMT-dependent weights in the logic tree cannot work '

openquake/hazardlib/calc/conditioned_gmfs.py

@@ -116,7 +116,6 @@
 from openquake.hazardlib.calc.gmf import GmfComputer
 from openquake.hazardlib.const import StdDev
 from openquake.hazardlib.geo.geodetic import geodetic_distance
-from openquake.hazardlib.contexts import ContextMaker
 
 U32 = numpy.uint32
 F32 = numpy.float32
@@ -215,12 +214,11 @@ def get_mea_tau_phi(self):
         :returns: a list of arrays [mea, sig, tau, phi]
         """
         return get_mean_covs(
-            self.rupture, self.cmaker.gsims,
+            self.rupture, self.cmaker,
             self.station_sitecol, self.station_data,
             self.observed_imt_strs, self.sitecol, self.imts,
-            self.spatial_correl,
-            self.cross_correl_between, self.cross_correl_within,
-            self.cmaker.maximum_distance, sigma=False)
+            self.spatial_correl, self.cross_correl_between, self.cross_correl_within,
+            sigma=False)
 
 
 @dataclass
@@ -388,10 +386,10 @@ def compute_spatial_cross_covariance_matrix(
 
 # tested in openquake/hazardlib/tests/calc/conditioned_gmfs_test.py
 def get_mean_covs(
-        rupture, gsims, station_sitecol, station_data,
+        rupture, cmaker, station_sitecol, station_data,
         observed_imt_strs, target_sitecol, target_imts,
         spatial_correl, cross_correl_between, cross_correl_within,
-        maximum_distance, sigma=True):
+        sigma=True):
     """
     :returns: a list of arrays [mea, sig, tau, phi] or [mea, tau, phi]
     """
@@ -409,20 +407,15 @@ def get_mean_covs(
 
     # Generate the contexts and calculate the means and
     # standard deviations at the *station* sites ("_D")
-    cmaker_D = ContextMaker(
-        rupture.tectonic_region_type, gsims,
-        dict(imtls=observed_imtls, maximum_distance=maximum_distance))
+    cmaker_D = cmaker.copy(imtls=observed_imtls)
 
     [ctx_D] = cmaker_D.get_ctx_iter([rupture], station_sitecol)
     mean_stds_D = cmaker_D.get_mean_stds([ctx_D])
     # shape (4, G, M, N) where 4 means (mean, sig, tau, phi)
 
     # Generate the contexts and calculate the means and 
     # standard deviations at the *target* sites ("_Y")
-    cmaker_Y = ContextMaker(
-        rupture.tectonic_region_type, gsims, dict(
-            imtls={target_imts[0].string: [0]},
-            maximum_distance=maximum_distance))
+    cmaker_Y = cmaker.copy(imtls={target_imts[0].string: [0]})
 
     [ctx_Y] = cmaker_Y.get_ctx_iter([rupture], target_sitecol)
     mean_stds_Y = cmaker_D.get_mean_stds([ctx_Y])
@@ -436,13 +429,13 @@ def get_mean_covs(
     compute_cov = partial(compute_spatial_cross_covariance_matrix,
                           spatial_correl, cross_correl_within)
 
-    G = len(gsims)
+    G = len(cmaker.gsims)
     M = len(target_imts)
     N = len(ctx_Y)
     me = numpy.zeros((G, M, N, 1))
     ta = numpy.zeros((G, M, N, N))
     ph = numpy.zeros((G, M, N, N))
-    for g, gsim in enumerate(gsims):
+    for g, gsim in enumerate(cmaker.gsims):
         if gsim.DEFINED_FOR_STANDARD_DEVIATION_TYPES == {StdDev.TOTAL}:
             if not (type(gsim).__name__ == "ModifiableGMPE"
                     and "add_between_within_stds" in gsim.kwargs):

openquake/hazardlib/calc/disagg.py

@@ -25,7 +25,6 @@
 import operator
 import collections
 import itertools
-from unittest.mock import Mock
 from functools import lru_cache
 import numpy
 import scipy.stats
@@ -43,7 +42,7 @@
 from openquake.hazardlib.tom import get_pnes
 from openquake.hazardlib.site import Site, SiteCollection
 from openquake.hazardlib.gsim.base import to_distribution_values
-from openquake.hazardlib.contexts import ContextMaker, FarAwayRupture
+from openquake.hazardlib.contexts import ContextMaker, Oq, FarAwayRupture
 from openquake.hazardlib.calc.mean_rates import (
     calc_rmap, calc_mean_rates, to_rates, to_probs)
 
@@ -634,19 +633,19 @@ def disaggregation(
     mags_by_trt = AccumDict(accum=set())
     dists = []
     tom = sources[0].temporal_occurrence_model
-    oq = Mock(imtls={str(imt): [iml]},
-              poes=[None],
-              rlz_index=[0],
-              epsstar=epsstar,
-              truncation_level=truncation_level,
-              investigation_time=tom.time_span,
-              maximum_distance=source_filter.integration_distance,
-              mags_by_trt=mags_by_trt,
-              num_epsilon_bins=n_epsilons,
-              mag_bin_width=mag_bin_width,
-              distance_bin_width=dist_bin_width,
-              coordinate_bin_width=coord_bin_width,
-              disagg_bin_edges=bin_edges)
+    oq = Oq(imtls={str(imt): [iml]},
+            poes=[None],
+            rlz_index=[0],
+            epsilon_star=epsstar,
+            truncation_level=truncation_level,
+            investigation_time=tom.time_span,
+            maximum_distance=source_filter.integration_distance,
+            mags_by_trt=mags_by_trt,
+            num_epsilon_bins=n_epsilons,
+            mag_bin_width=mag_bin_width,
+            distance_bin_width=dist_bin_width,
+            coordinate_bin_width=coord_bin_width,
+            disagg_bin_edges=bin_edges)
     for trt, srcs in by_trt.items():
         cmaker[trt] = cm = ContextMaker(trt, rlzs_by_gsim, oq)
         ctxs[trt].extend(cm.from_srcs(srcs, sitecol))

openquake/hazardlib/contexts.py

@@ -196,10 +196,15 @@ class Oq(object):
     """
     A mock for OqParam
     """
+    af = None
+    aristotle = False
+    cross_correl = None
     mea_tau_phi = False
     split_sources = True
     use_rates = False
-    af = None
+    with_betw_ratio = None
+    infer_occur_rates = False
+    inputs = ()
 
     def __init__(self, **hparams):
         vars(self).update(hparams)
@@ -512,6 +517,20 @@ def _set_poes(mean_std, loglevels, phi_b, out):
 
 # ############################ ContextMaker ############################### #
 
+
+def _fix(gsimdict, betw):
+    if betw:
+        out = {}
+        for gsim, uints in gsimdict.items():
+            if len(gsim.DEFINED_FOR_STANDARD_DEVIATION_TYPES) == 1:
+                out[valid.modified_gsim(gsim, add_between_within_stds=betw)] \
+                    = uints
+            else:
+                out[gsim] = uints
+        return out
+    return gsimdict
+
+
 class ContextMaker(object):
     """
     A class to manage the creation of contexts and to compute mean/stddevs
@@ -536,13 +555,8 @@ class ContextMaker(object):
 
     def __init__(self, trt, gsims, oq, monitor=Monitor(), extraparams=()):
         self.trt = trt
-        if isinstance(gsims, dict):
-            self.gsims = gsims
-        else:
-            self.gsims = {gsim: U32([i]) for i, gsim in enumerate(gsims)}
-        # NB: the gid array can be overridden later on
-        self.gid = numpy.arange(len(gsims), dtype=numpy.uint16)
         if isinstance(oq, dict):
+            # this happens when instantiating RuptureData in extract.py
             param = oq
             oq = Oq(**param)
             self.mags = param.get('mags', ())  # list of strings %.2f
@@ -560,6 +574,19 @@ def __init__(self, trt, gsims, oq, monitor=Monitor(), extraparams=()):
             except KeyError:  # missing TRT but there is only one
                 [(_, self.mags)] = oq.mags_by_trt.items()
 
+        if oq.with_betw_ratio:
+            betw_ratio = {'with_betw_ratio': oq.with_betw_ratio}
+        elif oq.aristotle:
+            betw_ratio = {'with_betw_ratio': 1.7}  # same as in GEESE
+        else:
+            betw_ratio = {}
+        if isinstance(gsims, dict):
+            self.gsims = _fix(gsims, betw_ratio)
+        else:
+            self.gsims = _fix({gsim: U32([i]) for i, gsim in enumerate(gsims)},
+                              betw_ratio)
+        # NB: the gid array can be overridden later on
+        self.gid = numpy.arange(len(gsims), dtype=numpy.uint16)
         self.oq = oq
         self.monitor = monitor
         self._init1(param)
@@ -576,8 +603,7 @@ def _init1(self, param):
                     raise TypeError('Expected string, got %s' % type(imt))
             self.imtls = param['imtls']
         elif 'hazard_imtls' in param:
-            self.imtls = DictArray(
-                imt_module.sort_by_imt(param['hazard_imtls']))
+            self.imtls = imt_module.dictarray(param['hazard_imtls'])
         elif not hasattr(self, 'imtls'):
             raise KeyError('Missing imtls in ContextMaker!')
         self.cache_distances = param.get('cache_distances', False)
@@ -671,6 +697,17 @@ def init_monitoring(self, monitor):
         self.out_no = getattr(monitor, 'out_no', self.task_no)
         self.cfactor = numpy.zeros(2)
 
+    def copy(self, **kw):
+        """
+        :returns: a copy of the ContextMaker with modified attributes
+        """
+        new = copy.copy(self)
+        for k, v in kw.items():
+            setattr(new, k, v)
+        if 'imtls' in kw:
+            new.set_imts_conv()
+        return new
+
     def restrict(self, imts):
         """
         :param imts: a list of IMT strings subset of the full list
@@ -684,7 +721,7 @@ def restrict(self, imts):
     def set_imts_conv(self):
         """
         Set the .imts list and .conv dictionary for the horizontal component
-        conversion (if any).
+        conversion (if any). Also set the .loglevels.
         """
         self.loglevels = DictArray(self.imtls) if self.imtls else {}
         with warnings.catch_warnings():

openquake/hazardlib/imt.py

@@ -23,6 +23,7 @@
 import re
 import collections
 import numpy
+from openquake.baselib.general import DictArray
 
 FREQUENCY_PATTERN = '^(EAS|FAS|DRVT|AvgSA)\\((\\d+\\.*\\d*)\\)'
 
@@ -97,6 +98,13 @@ def sort_by_imt(imtls):
     return {imt: imtls[imt] for imt in imts}
 
 
+def dictarray(imtls):
+    """
+    :returns: a DictArray sorted by IMT
+    """
+    return DictArray(sort_by_imt(imtls))
+
+
 def repr(self):
     if self.period and self.damping != 5.0:
         if self.string.startswith('SDi'):