diff --git a/debian/changelog b/debian/changelog index e86d52d32..a64cafb06 100644 --- a/debian/changelog +++ b/debian/changelog @@ -1,4 +1,5 @@ [Michele Simionato] + * Avoided doing the rupture filtering twice * Optimized the case of multiple GSIMs by instantiating the ruptures/sites/distances contexts only once diff --git a/openquake/hazardlib/calc/hazard_curve.py b/openquake/hazardlib/calc/hazard_curve.py index 9b2b64f08..5d8e95e49 100644 --- a/openquake/hazardlib/calc/hazard_curve.py +++ b/openquake/hazardlib/calc/hazard_curve.py @@ -27,7 +27,7 @@ from openquake.baselib.python3compat import range, raise_ from openquake.baselib.performance import DummyMonitor from openquake.hazardlib.calc import filters -from openquake.hazardlib.gsim.base import ContextMaker +from openquake.hazardlib.gsim.base import ContextMaker, FarAwayRupture from openquake.hazardlib.imt import from_string from openquake.baselib.general import deprecated @@ -93,7 +93,8 @@ def hazard_curves( def calc_hazard_curves( sources, sites, imtls, gsim_by_trt, truncation_level=None, source_site_filter=filters.source_site_noop_filter, - rupture_site_filter=filters.rupture_site_noop_filter): + rupture_site_filter=filters.rupture_site_noop_filter, + maximum_distance=None): """ Compute hazard curves on a list of sites, given a set of seismic sources and a set of ground shaking intensity models (one per tectonic region type @@ -158,11 +159,13 @@ def calc_hazard_curves( return curves +# TODO: remove the rupture_site_filter, since its work is now done by the +# maximum_distance parameter; see what would break def hazard_curves_per_trt( sources, sites, imtls, gsims, truncation_level=None, source_site_filter=filters.source_site_noop_filter, rupture_site_filter=filters.rupture_site_noop_filter, - monitor=DummyMonitor()): + maximum_distance=None, monitor=DummyMonitor()): """ Compute the hazard curves for a set of sources belonging to the same tectonic region type for all the GSIMs associated to that TRT. @@ -175,7 +178,7 @@ def hazard_curves_per_trt( by the intensity measure types; the size of each field is given by the number of levels in ``imtls``. """ - cmaker = ContextMaker(gsims) + cmaker = ContextMaker(gsims, maximum_distance) gnames = list(map(str, gsims)) imt_dt = numpy.dtype([(imt, float, len(imtls[imt])) for imt in sorted(imtls)]) @@ -194,7 +197,11 @@ def hazard_curves_per_trt( (rupture, s_sites) for rupture in source.iter_ruptures())) for rupture, r_sites in rupture_sites: with ctx_mon: - sctx, rctx, dctx = cmaker.make_contexts(r_sites, rupture) + try: + sctx, rctx, dctx = cmaker.make_contexts( + r_sites, rupture) + except FarAwayRupture: + continue for i, gsim in enumerate(gsims): with pne_mon: for imt in imts: @@ -202,7 +209,7 @@ def hazard_curves_per_trt( sctx, rctx, dctx, imt, imts[imt], truncation_level) pno = rupture.get_probability_no_exceedance(poes) - expanded_pno = r_sites.expand(pno, placeholder=1) + expanded_pno = sctx.sites.expand(pno, 1.0) curves[i][str(imt)] *= expanded_pno except Exception as err: etype, err, tb = sys.exc_info() diff --git a/openquake/hazardlib/gsim/base.py b/openquake/hazardlib/gsim/base.py index d14f8dea9..f4a6ce677 100644 --- a/openquake/hazardlib/gsim/base.py +++ b/openquake/hazardlib/gsim/base.py @@ -108,21 +108,52 @@ def forbid_instantiation(cls): cls.instantiable = True +def get_distances(rupture, mesh, param='rjb'): + """ + :param rupture: a rupture + :param mesh: a mesh of points + :param param: the kind of distance to compute (default rjb) + :returns: an array of distances from the given mesh + """ + if param == 'rrup': + dist = rupture.surface.get_min_distance(mesh) + elif param == 'rx': + dist = rupture.surface.get_rx_distance(mesh) + elif param == 'ry0': + dist = rupture.surface.get_ry0_distance(mesh) + elif param == 'rjb': + dist = rupture.surface.get_joyner_boore_distance(mesh) + elif param == 'rhypo': + dist = rupture.hypocenter.distance_to_mesh(mesh) + elif param == 'repi': + dist = rupture.hypocenter.distance_to_mesh(mesh, with_depths=False) + elif param == 'rcdpp': + dist = rupture.get_cdppvalue(mesh) + else: + raise ValueError('Unknown distance measure %r' % param) + return dist + + +class FarAwayRupture(Exception): + """Raised if the rupture is outside the maximum distance for all sites""" + + class ContextMaker(object): """ A class to manage the creation of contexts for distances, sites, rupture. """ REQUIRES = ['DISTANCES', 'SITES_PARAMETERS', 'RUPTURE_PARAMETERS'] - def __init__(self, gsims): + def __init__(self, gsims, maximum_distance=None): self.gsims = gsims + self.maximum_distance = maximum_distance for req in self.REQUIRES: reqset = set() for gsim in gsims: reqset.update(getattr(gsim, 'REQUIRES_' + req)) setattr(self, 'REQUIRES_' + req, reqset) - def make_distances_context(self, site_collection, rupture): + def make_distances_context(self, site_collection, rupture, dist_dict=()): """ Create distances context object for given site collection and rupture. @@ -136,6 +167,9 @@ def make_distances_context(self, site_collection, rupture): :class: `~openquake.hazardlib.source.rupture.BaseProbabilisticRupture`). + :param dist_dict: + A dictionary of already computed distances, keyed by distance name + :returns: Source to site distances as instance of :class: `DistancesContext()`. Only those values that are required by GSIM @@ -146,30 +180,11 @@ def make_distances_context(self, site_collection, rupture): """ dctx = DistancesContext() for param in self.REQUIRES_DISTANCES: - if param == 'rrup': - dist = rupture.surface.get_min_distance(site_collection.mesh) - elif param == 'rx': - dist = rupture.surface.get_rx_distance(site_collection.mesh) - elif param == 'ry0': - dist = rupture.surface.get_ry0_distance(site_collection.mesh) - elif param == 'rjb': - dist = rupture.surface.get_joyner_boore_distance( - site_collection.mesh - ) - elif param == 'rhypo': - dist = rupture.hypocenter.distance_to_mesh( - site_collection.mesh - ) - elif param == 'repi': - dist = rupture.hypocenter.distance_to_mesh( - site_collection.mesh, with_depths=False - ) - elif param == 'rcdpp': - dist = rupture.get_cdppvalue(site_collection.mesh) + if param in dist_dict: # already computed distances + distances = dist_dict[param] else: - raise ValueError('%s requires unknown distance measure %r' % - (type(self).__name__, param)) - setattr(dctx, param, dist) + distances = get_distances(rupture, site_collection.mesh, param) + setattr(dctx, param, distances) return dctx def make_sites_context(self, site_collection): @@ -189,6 +204,7 @@ def make_sites_context(self, site_collection): """ sctx = SitesContext() + sctx.sites = site_collection for param in self.REQUIRES_SITES_PARAMETERS: try: value = getattr(site_collection, param) @@ -204,9 +220,8 @@ def make_rupture_context(self, rupture): :param rupture: Instance of - :class:`~openquake.hazardlib.source.rupture.Rupture` (or - subclass of - :class:`~openquake.hazardlib.source.rupture.BaseProbabilisticRupture`). + :class:`openquake.hazardlib.source.rupture.Rupture` or subclass of + :class:`openquake.hazardlib.source.rupture.BaseProbabilisticRupture` :returns: Rupture parameters as instance of :class: @@ -244,16 +259,16 @@ def make_rupture_context(self, rupture): def make_contexts(self, site_collection, rupture): """ - Create context objects for given site collection and rupture. + Filter the site collection with respect to the rupture and + create context objects. :param site_collection: Instance of :class:`openquake.hazardlib.site.SiteCollection`. :param rupture: Instance of - :class:`~openquake.hazardlib.source.rupture.Rupture` (or - subclass of - :class:`~openquake.hazardlib.source.rupture.BaseProbabilisticRupture`). + :class:`openquake.hazardlib.source.rupture.Rupture` or subclass of + :class:`openquake.hazardlib.source.rupture.BaseProbabilisticRupture` :returns: Tuple of three items: sites context, rupture context and @@ -267,9 +282,20 @@ def make_contexts(self, site_collection, rupture): If any of declared required parameters (that includes site, rupture and distance parameters) is unknown. """ - return (self.make_sites_context(site_collection), - self.make_rupture_context(rupture), - self.make_distances_context(site_collection, rupture)) + rctx = self.make_rupture_context(rupture) + distances = get_distances(rupture, site_collection.mesh, 'rjb') + sites = site_collection + if self.maximum_distance: + mask = distances <= self.maximum_distance + if mask.any(): + sites = site_collection.filter(mask) + distances = distances[mask] + else: + raise FarAwayRupture + + sctx = self.make_sites_context(sites) + dctx = self.make_distances_context(sites, rupture, {'rjb': distances}) + return (sctx, rctx, dctx) @functools.total_ordering diff --git a/openquake/hazardlib/tests/calc/hazard_curve_test.py b/openquake/hazardlib/tests/calc/hazard_curve_test.py index 088f22ff7..4c4fee116 100644 --- a/openquake/hazardlib/tests/calc/hazard_curve_test.py +++ b/openquake/hazardlib/tests/calc/hazard_curve_test.py @@ -26,6 +26,12 @@ from openquake.hazardlib.gsim.base import ContextMaker +class FakeSiteContext(object): + def __init__(self, sites): + self.sites = sites + self.mesh = sites.mesh + + class HazardCurvesTestCase(unittest.TestCase): class FakeRupture(object): def __init__(self, probability, tectonic_region_type): @@ -71,7 +77,7 @@ def __str__(self): def setUp(self): self.orig_make_contexts = ContextMaker.make_contexts ContextMaker.make_contexts = lambda self, sites, rupture: ( - sites, rupture, None) + FakeSiteContext(sites), rupture, None) self.truncation_level = 3.4 self.imts = {'PGA': [1, 2, 3], 'PGD': [2, 4]} self.time_span = 49.2 diff --git a/openquake/hazardlib/tests/gsim/base_test.py b/openquake/hazardlib/tests/gsim/base_test.py index 0991fe6bc..c28bb9b68 100644 --- a/openquake/hazardlib/tests/gsim/base_test.py +++ b/openquake/hazardlib/tests/gsim/base_test.py @@ -428,7 +428,7 @@ def make_contexts(self, site_collection, rupture): def test_unknown_site_param_error(self): self.gsim_class.REQUIRES_SITES_PARAMETERS.add('unknown!') err = "ContextMaker requires unknown site parameter 'unknown!'" - sites = SiteCollection([self.site1]) + sites = SiteCollection([self.site1, self.site2]) self._assert_value_error(self.make_contexts, err, site_collection=sites, rupture=self.rupture) @@ -441,8 +441,8 @@ def test_unknown_rupture_param_error(self): def test_unknown_distance_error(self): self.gsim_class.REQUIRES_DISTANCES.add('jump height') - err = "ContextMaker requires unknown distance measure 'jump height'" - sites = SiteCollection([self.site1]) + err = "Unknown distance measure 'jump height'" + sites = SiteCollection([self.site1, self.site2]) self._assert_value_error(self.make_contexts, err, site_collection=sites, rupture=self.rupture)