Issue computing NGAEast GMMs with vectorized input

[smithUSGS]

Hi GEM:

If I vectorize the GMM inputs and run, say BSSA14, everything works as expected (and is much faster compared to previous versions). I've attached a minimally reproducible example:

import numpy as np

from openquake.hazardlib.imt import PGA
from openquake.hazardlib.gsim.boore_2014 import BooreEtAl2014

n = 10
ctx = np.recarray(
    n,
    dtype=np.dtype(
        [
            ("mag", float),
            ("vs30", float),
            ("rrup", float),
            ("rjb", float),
            ("rake", float),
        ]
    ),
)

ctx["mag"] = np.linspace(5, 7, n)
ctx["vs30"] = np.linspace(200, 300, n)
ctx["rjb"] = np.linspace(1, 100, n)
ctx["rrup"] = np.linspace(1, 100, n)
ctx["rake"] = np.linspace(0, 90, n)

mean = np.zeros([1, n])
sigma = np.zeros_like(mean)
tau = np.zeros_like(mean)
phi = np.zeros_like(mean)

gsim = BooreEtAl2014()
gsim.compute(ctx, [PGA()], mean, sigma, tau, phi)

However, if I attempt to vectorize the input to NGAEast GMMs I get an error and I don't know if it is an issue with the NGAEast models requiring the magnitude to be inputted as a scalar or if the code needs to be patched.

from gmm_weighting.scenarios.nga_east import NGAEast

nsc = len(sc.mag)
gmm_mods, gmm_names = nga_east.get_gmms()
nmod = len(gmm_names)

gmms = np.full((nmod, nsc), np.nan)
ctx = np.recarray(nsc, dtype=np.dtype([("mag", float), ("rrup", float), ("vs30", float),]),)
ctx["mag"] = [4., 4.5, 5.]
ctx["rrup"] = [100., 150., 200.]
ctx["vs30"] = 3000.0

mean = np.zeros([1, nsc])
sigma = np.zeros_like(mean)
tau = np.zeros_like(mean)
phi = np.zeros_like(mean)

gmm_mods[0].compute(ctx, [SA(1.0)], mean, sigma, tau, phi)

The error I get is from the line highlighted below:

    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
        """
        Returns the mean and standard deviations
        """
### >       [mag] = np.unique(np.round(ctx.mag, 2))
       ValueError: too many values to unpack (expected 1)

If I look at other GMM implementations, they pass the recarray (ctx), which makes sense. But here that is not done and I don't know if it's related to these GMMs specifically or some other reason.

Please advise and thank you for your work on this python package.

[micheles]

Using a GSIM is non-obvious. First of all, it should not be instantiated directly, but through the valid.gsim function. Secondly, contexts objects should be instantiated by using a ContextMaker instance, which has to be instantiated carefully. Thirdly, the .compute method should not be called directly. In short, here is how to do what you want:

import numpy as np
from openquake.hazardlib import valid
from openquake.hazardlib.contexts import simple_cmaker

gsim = valid.gsim('''[NGAEastGMPE]
gmpe_table="NGAEast_PEER_EX.hdf5"''')
n = 10
mags = np.linspace(5, 7, n)
cmaker = simple_cmaker([gsim], ['PGA'], mags=['%.2f' % mag for mag in mags])
ctx = cmaker.new_ctx(n)
ctx["mag"] = mags
ctx["vs30"] = np.linspace(200, 300, n)
ctx["rrup"] = np.linspace(1, 100, n)

mean, sigma, tau, phi = cmaker.get_mean_stds([ctx])

[smithUSGS]

Thank you, I adapted your example and it passed the test. Note that I didn't need to instantiate the GMM using valid.gsim (although I can and it is what you recommend). If you don't mind explaining why using valid.gsim is a good practice, I would appreciate that. Either way, I think this issue can be closed now.

[micheles]

Using valid.gsim is a good practice because it is the official way the engine convert strings into GMPEs, so it is guaranteed to work forever, while manual instantiation changed in the past and could potentially change again in the future.

[emthompson-usgs]

Thanks @micheles, we weren't aware of this. Good to know.