Merge pull request #10072 from gem/median

Investigating median spectrum weights

openquake/calculators/postproc/median_spectrum.py

@@ -198,7 +198,12 @@ def main(dstore, csm):
                            site_id=N, kind=['mea', 'sig', 'wei'],
                            period=periods, poe=oq.poes)
     # sanity check on the weights
-    # np.testing.assert_allclose(tot_w, 1, rtol=.01)
+    for p, poe in enumerate(oq.poes):
+        maxw = tot_w[:, :, p].max()
+        logging.info(f'{poe=} {maxw=}')
+        if (np.abs(tot_w[:, :, p] - 1) > .01).any():
+            raise ValueError(f'The weights sum up to {maxw:.3f} != 1: perhaps the '
+                             f'hazard curve is not invertible around {poe=}')
 
     # sanity check on the rup_ids
     if N == 1 and P == 1:

openquake/commands/plot.py

@@ -83,6 +83,8 @@ def make_figure_hcurves(extractors, what):
             if (arr == 0).all():
                 logging.warning('There is a zero curve %s_%s', *ck)
             ax.loglog(imls, arr.flat, '-', label='%s_%s' % ck)
+        for poe in oq.poes:
+            ax.plot(imls, [poe]*len(imls), label=f'{poe=}')
         ax.grid(True)
         ax.legend()
     return plt

openquake/hazardlib/valid.py

@@ -929,6 +929,23 @@ def logscale(x_min, x_max, n):
     return numpy.exp(delta * numpy.arange(n) / (n - 1)) * x_min
 
 
+def linscale(x_min, x_max, n):
+    """
+    :param x_min: minumum value
+    :param x_max: maximum value
+    :param n: number of steps
+    :returns: an array of n values from x_min to x_max
+    """
+    if not (isinstance(n, int) and n > 0):
+        raise ValueError('n must be a positive integer, got %s' % n)
+    if x_min <= 0:
+        raise ValueError('x_min must be positive, got %s' % x_min)
+    if x_max <= x_min:
+        raise ValueError('x_max (%s) must be bigger than x_min (%s)' %
+                         (x_max, x_min))
+    return numpy.linspace(x_min, x_max, num=n)
+
+
 def dictionary(value):
     """
     :param value:
@@ -951,18 +968,32 @@ def dictionary(value):
     """
     if not value:
         return {}
-    value = value.replace('logscale(', '("logscale", ')  # dirty but quick
+
+    if 'logscale' in value:
+        value = value.replace('logscale(', '("logscale", ')  # dirty but quick
+    if 'linscale' in value:
+        value = value.replace('linscale(', '("linscale", ')  # dirty but quick
+
     try:
         dic = dict(ast.literal_eval(value))
     except Exception:
         raise ValueError('%r is not a valid Python dictionary' % value)
+
     for key, val in dic.items():
         try:
             has_logscale = (val[0] == 'logscale')
-        except Exception:  # no val[0]
+        except IndexError:  # no val[0]
             continue
         if has_logscale:
             dic[key] = list(logscale(*val[1:]))
+
+        try:
+            has_linscale = (val[0] == 'linscale')
+        except IndexError:  # no val[0]
+            continue
+        if has_linscale:
+            dic[key] = list(linscale(*val[1:]))
+
     return dic