code: def rdp_gaussian(logq, sigma, orders):
      """Bounds RDP from above of GNMax given an upper bound on q (Theorem 6).
  
      Args:
        logq: Natural logarithm of the probability of a non-argmax outcome.
        sigma: Standard deviation of Gaussian noise.
        orders: An array_like list of Renyi orders.
  
      Returns:
        Upper bound on RPD for all orders. A scalar if orders is a scalar.
  
      Raises:
        ValueError: If the input is malformed.
      """
      if (
          logq > 0 or sigma < 0 or np.any(orders <= 1)
      ):  # not defined for alpha=1
          raise ValueError("Inputs are malformed.")
  
      if np.isneginf(
          logq
      ):  # If the mechanism's output is fixed, it has 0-DP.
          if np.isscalar(orders):
              return 0.0
          else:
              return np.full_like(orders, 0.0, dtype=float)
  
      variance = sigma**2
  
      # Use two different higher orders: mu_hi1 and mu_hi2 computed according to
      # Proposition 10.
      mu_hi2 = math.sqrt(variance * -logq)
      mu_hi1 = mu_hi2 + 1
  
      orders_vec = np.atleast_1d(orders)
  
      ret = orders_vec / variance  # baseline: data-independent bound
  
      # Filter out entries where data-dependent bound does not apply.
      mask = np.logical_and(mu_hi1 > orders_vec, mu_hi2 > 1)
  
      rdp_hi1 = mu_hi1 / variance
      rdp_hi2 = mu_hi2 / variance
  
      log_a2 = (mu_hi2 - 1) * rdp_hi2
  
      # Make sure q is in the increasing wrt q range and A is positive.
      if (
          np.any(mask)
          and logq
          <= log_a2
          - mu_hi2
          * (
              math.log(1 + 1 / (mu_hi1 - 1))
              + math.log(1 + 1 / (mu_hi2 - 1))
          )
          and -logq > rdp_hi2
      ):
          # Use log1p(x) = log(1 + x) to avoid catastrophic cancellations when x ~ 0.
          log1q = _log1mexp(logq)  # log1q = log(1-q)
          log_a = (orders - 1) * (
              log1q - _log1mexp((logq + rdp_hi2) * (1 - 1 / mu_hi2))
          )
          log_b = (orders - 1) * (rdp_hi1 - logq / (mu_hi1 - 1))
  
          # Use logaddexp(x, y) = log(e^x + e^y) to avoid overflow for large x, y.
          log_s = np.logaddexp(log1q + log_a, logq + log_b)
          ret[mask] = np.minimum(ret, log_s / (orders - 1))[mask]
  
      assert np.all(ret >= 0)
  
      if np.isscalar(orders):
          return ret.item()
      else:
          return ret
  
imports: from __future__ import absolute_import
  from __future__ import division
  from __future__ import print_function
  import math
  from absl import app
  import numpy
  import scipy.stats
label_desc: Use image node_id 6 for calling a global function with example usage: rdp_gaussian(logq, sigma, orders) and returns: ret, ret, int, np