NGAEastUSGSGMPE in 3.15

[emthompson-usgs]

Hi, the way we (including @cbworden) have been calling NGA East GMMs no longer seems to work in 3.15 (the last version that we were using this code and it was successful is 3.13). Here is a minimal example that reflects how we have been calling the module:

import os
import numpy as np

from openquake.hazardlib.gsim.usgs_ceus_2019 import NGAEastUSGSGMPE
from openquake.hazardlib.contexts import RuptureContext
from openquake.hazardlib.imt import PGA
from openquake.hazardlib import const

TABLE_PATHS = os.listdir(NGAEastUSGSGMPE.PATH)
SIGMA_MODS = ["EPRI", "PANEL"]

imt = PGA()
std = const.StdDev.TOTAL

ctx = RuptureContext()
ctx.rrup = np.array([10.0])
ctx.vs30 = np.array([340.0])
ctx.mag = 6.0
ctx.sids = np.array([1], dtype=int)
ctx.occurrence_rate = 1e-5

gmpe = NGAEastUSGSGMPE(gmpe_table=TABLE_PATHS[0], sigma_model=SIGMA_MODS[0])
gmpe.get_mean_and_stddevs(ctx, ctx, ctx, imt, std)

Here is the output with OQ 3.13:

Out[8]: (array([-0.0437124]), [])

Here is the error we get with 3.15:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [2], in <cell line: 2>()
      1 gmpe = NGAEastUSGSGMPE(gmpe_table=TABLE_PATHS[0], sigma_model=SIGMA_MODS[0])
----> 2 gmpe.get_mean_and_stddevs(ctx, ctx, ctx, imt, std)

File ~/miniconda3/envs/shakemap/lib/python3.9/site-packages/openquake/hazardlib/gsim/base.py:360, in GroundShakingIntensityModel.get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types)
    358     ctx = rup  # rup is already a good object
    359 if self.compute.__annotations__.get("ctx") is numpy.recarray:
--> 360     cmaker = ContextMaker('*', [self], {'imtls': {imt.string: [0]}})
    361     if not isinstance(ctx, numpy.ndarray):
    362         ctx = cmaker.recarray([ctx])

File ~/miniconda3/envs/shakemap/lib/python3.9/site-packages/openquake/hazardlib/contexts.py:403, in ContextMaker.__init__(self, trt, gsims, oq, monitor, extraparams)
    401     if hasattr(gsim, 'set_tables'):
    402         if not self.mags and not is_modifiable(gsim):
--> 403             raise ValueError(
    404                 'You must supply a list of magnitudes as 2-digit '
    405                 'strings, like mags=["6.00", "6.10", "6.20"]')
    406         gsim.set_tables(self.mags, self.imtls)
    407 self.maximum_distance = _interp(param, 'maximum_distance', trt)

ValueError: You must supply a list of magnitudes as 2-digit strings, like mags=["6.00", "6.10", "6.20"]

Is there a different/new way to call the NGAEastUSGSGMPE function without having to deal with the set_tables issue?

Any guidance on how we need to adjust our code to make this work will be greatly appreciated.

[micheles]

Yeah, I remember that. There were a few GMPEs for Canada which were table-based and the engine was interpolating the rupture magnitude versus the magnitudes in the HDF5 table for each rupture. Therefore with 1 million ruptures it was doing 1 million interpolations and it was slow. The change was to pass to the ContextMaker the list of unique magnitudes in the set of ruptures (say 20) so that the interpolation could be done in advance, in the master core, even before starting the calculation, and only 20 times instead of 1 million times.

You are using the obsolete interface get_mean_and_stddevs that does not have the ability to pass the list of magnitudes to the underlying ContextMaker. So I would suggest to use the ContextMaker directly instead, something like:

from openquake.hazardlib import valid
from openquake.hazardlib.contexts import ContextMaker

gsim = valid.gsim('NGAEastUSGSSeedSP15')
param = dict(maximum_distance=200, imtls={'PGA': [.1, .2, .3]},
             truncation_level=3., investigation_time=1.,
             mags=["6.00", "6.10", "6.20"])
cmaker = ContextMaker('TRT', [gsim], param)

ctx = cmaker.new_ctx(size=1)  # recarray of size 1
ctx.mag = 6.
ctx.vs30 = 700.
ctx.rrup = 100.
print(ctx.dtype.names)
print(ctx)
mean, sig, tau, phi = cmaker.get_mean_stds([ctx])
print(mean, sig)

NB: the method new_ctx is in the engine-3.15 branch but not in the released code, it will be once we make the first patch to v3.15. If you are using git, just do a git pull.

[micheles]

I am not hearing back, so I am assuming all is good and closing the issue.

[emthompson-usgs]

Yes, sorry I hadn't responded back yet. Thank you for the help on this.