docs: def get_quantile_levels(density, x, y, xp, yp, q, normalize=True):
      """Compute density levels for given quantiles by interpolation
  
      For a given 2D density, compute the density levels at which
      the resulting contours contain the fraction `1-q` of all
      data points. E.g. for a measurement of 1000 events, all
      contours at the level corresponding to a quantile of
      `q=0.95` (95th percentile) contain 50 events (5%).
  
      Parameters
      ----------
      density: 2d ndarray of shape (M, N)
          Kernel density estimate for which to compute the contours
      x: 2d ndarray of shape (M, N) or 1d ndarray of size M
          X-values corresponding to `kde`
      y: 2d ndarray of shape (M, N) or 1d ndarray of size M
          Y-values corresponding to `kde`
      xp: 1d ndarray of size D
          Event x-data from which to compute the quantile
      yp: 1d ndarray of size D
          Event y-data from which to compute the quantile
      q: array_like or float between 0 and 1
          Quantile along which to find contours in `kde` relative
          to its maximum
      normalize: bool
          Whether output levels should be normalized to the maximum
          of `density`
  
      Returns
      -------
      level: float
          Contours level corresponding to the given quantile
  
      Notes
      -----
      NaN-values events in `xp` and `yp` are ignored.
      """
      # xy coordinates
      if len(x.shape) == 2:
          assert np.all(x[:, 0] == x[:, 1])
          x = x[:, 0]
      if len(y.shape) == 2:
          assert np.all(y[0, :] == y[1, :])
          y = y[0, :]
  
      # remove bad events
      bad = get_bad_vals(xp, yp)
      xp = xp[~bad]
      yp = yp[~bad]
  
      # Normalize interpolation data such that the spacing for
      # x and y is about the same during interpolation.
      x_norm = x.max()
      x = x / x_norm
      xp = xp / x_norm
  
      y_norm = y.max()
      y = y / y_norm
      yp = yp / y_norm
  
      # Perform interpolation
      dp = spint.interpn((x, y), density,
                         (xp, yp),
                         method='linear',
                         bounds_error=False,
                         fill_value=0)
  
      if normalize:
          dp /= density.max()
  
      if not np.isscalar(q):
          q = np.array(q)
      plev = np.nanpercentile(dp, q=q*100)
      return plev