code: def qmc_quad(
      func, a, b, *, n_estimates=8, n_points=1024, qrng=None, log=False
  ):
      """
      Compute an integral in N-dimensions using Quasi-Monte Carlo quadrature.
  
      Parameters
      ----------
      func : callable
          The integrand. Must accept a single argument ``x``, an array which
          specifies the point(s) at which to evaluate the scalar-valued
          integrand, and return the value(s) of the integrand.
          For efficiency, the function should be vectorized to accept an array of
          shape ``(d, n_points)``, where ``d`` is the number of variables (i.e.
          the dimensionality of the function domain) and `n_points` is the number
          of quadrature points, and return an array of shape ``(n_points,)``,
          the integrand at each quadrature point.
      a, b : array-like
          One-dimensional arrays specifying the lower and upper integration
          limits, respectively, of each of the ``d`` variables.
      n_estimates, n_points : int, optional
          `n_estimates` (default: 8) statistically independent QMC samples, each
          of `n_points` (default: 1024) points, will be generated by `qrng`.
          The total number of points at which the integrand `func` will be
          evaluated is ``n_points * n_estimates``. See Notes for details.
      qrng : `~scipy.stats.qmc.QMCEngine`, optional
          An instance of the QMCEngine from which to sample QMC points.
          The QMCEngine must be initialized to a number of dimensions ``d``
          corresponding with the number of variables ``x1, ..., xd`` passed to
          `func`.
          The provided QMCEngine is used to produce the first integral estimate.
          If `n_estimates` is greater than one, additional QMCEngines are
          spawned from the first (with scrambling enabled, if it is an option.)
          If a QMCEngine is not provided, the default `scipy.stats.qmc.Halton`
          will be initialized with the number of dimensions determine from
          the length of `a`.
      log : boolean, default: False
          When set to True, `func` returns the log of the integrand, and
          the result object contains the log of the integral.
  
      Returns
      -------
      result : object
          A result object with attributes:
  
          integral : float
              The estimate of the integral.
          standard_error :
              The error estimate. See Notes for interpretation.
  
      Notes
      -----
      Values of the integrand at each of the `n_points` points of a QMC sample
      are used to produce an estimate of the integral. This estimate is drawn
      from a population of possible estimates of the integral, the value of
      which we obtain depends on the particular points at which the integral
      was evaluated. We perform this process `n_estimates` times, each time
      evaluating the integrand at different scrambled QMC points, effectively
      drawing i.i.d. random samples from the population of integral estimates.
      The sample mean :math:`m` of these integral estimates is an
      unbiased estimator of the true value of the integral, and the standard
      error of the mean :math:`s` of these estimates may be used to generate
      confidence intervals using the t distribution with ``n_estimates - 1``
      degrees of freedom. Perhaps counter-intuitively, increasing `n_points`
      while keeping the total number of function evaluation points
      ``n_points * n_estimates`` fixed tends to reduce the actual error, whereas
      increasing `n_estimates` tends to decrease the error estimate.
  
      Examples
      --------
      QMC quadrature is particularly useful for computing integrals in higher
      dimensions. An example integrand is the probability density function
      of a multivariate normal distribution.
  
      >>> import numpy as np
      >>> from scipy import stats
      >>> dim = 8
      >>> mean = np.zeros(dim)
      >>> cov = np.eye(dim)
      >>> def func(x):
      ...     # `multivariate_normal` expects the _last_ axis to correspond with
      ...     # the dimensionality of the space, so `x` must be transposed
      ...     return stats.multivariate_normal.pdf(x.T, mean, cov)
  
      To compute the integral over the unit hypercube:
  
      >>> from scipy.integrate import qmc_quad
      >>> a = np.zeros(dim)
      >>> b = np.ones(dim)
      >>> rng = np.random.default_rng()
      >>> qrng = stats.qmc.Halton(d=dim, seed=rng)
      >>> n_estimates = 8
      >>> res = qmc_quad(func, a, b, n_estimates=n_estimates, qrng=qrng)
      >>> res.integral, res.standard_error
      (0.00018429555666024108, 1.0389431116001344e-07)
  
      A two-sided, 99% confidence interval for the integral may be estimated
      as:
  
      >>> t = stats.t(df=n_estimates-1, loc=res.integral,
      ...             scale=res.standard_error)
      >>> t.interval(0.99)
      (0.0001839319802536469, 0.00018465913306683527)
  
      Indeed, the value reported by `scipy.stats.multivariate_normal` is
      within this range.
  
      >>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
      0.00018430867675187443
  
      """
      args = _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log)
      func, a, b, n_points, n_estimates, qrng, rng, log, stats = args
  
      def sum_product(integrands, dA, log=False):
          if log:
              return logsumexp(integrands) + np.log(dA)
          else:
              return np.sum(integrands * dA)
  
      def mean(estimates, log=False):
          if log:
              return logsumexp(estimates) - np.log(n_estimates)
          else:
              return np.mean(estimates)
  
      def std(estimates, m=None, ddof=0, log=False):
          m = m or mean(estimates, log)
          if log:
              estimates, m = np.broadcast_arrays(estimates, m)
              temp = np.vstack((estimates, m + np.pi * 1j))
              diff = logsumexp(temp, axis=0)
              return np.real(
                  0.5
                  * (logsumexp(2 * diff) - np.log(n_estimates - ddof))
              )
          else:
              return np.std(estimates, ddof=ddof)
  
      def sem(estimates, m=None, s=None, log=False):
          m = m or mean(estimates, log)
          s = s or std(estimates, m, ddof=1, log=log)
          if log:
              return s - 0.5 * np.log(n_estimates)
          else:
              return s / np.sqrt(n_estimates)
  
      # The sign of the integral depends on the order of the limits. Fix this by
      # ensuring that lower bounds are indeed lower and setting sign of resulting
      # integral manually
      if np.any(a == b):
          message = (
              "A lower limit was equal to an upper limit, so the value "
              "of the integral is zero by definition."
          )
          warnings.warn(message, stacklevel=2)
          return QMCQuadResult(-np.inf if log else 0, 0)
  
      i_swap = b < a
      sign = (-1) ** (i_swap.sum(axis=-1))  # odd # of swaps -> negative
      a[i_swap], b[i_swap] = b[i_swap], a[i_swap]
  
      A = np.prod(b - a)
      dA = A / n_points
  
      estimates = np.zeros(n_estimates)
      rngs = _rng_spawn(qrng.rng, n_estimates)
      for i in range(n_estimates):
          # Generate integral estimate
          sample = qrng.random(n_points)
          # The rationale for transposing is that this allows users to easily
          # unpack `x` into separate variables, if desired. This is consistent
          # with the `xx` array passed into the `scipy.integrate.nquad` `func`.
          x = stats.qmc.scale(sample, a, b).T  # (n_dim, n_points)
          integrands = func(x)
          estimates[i] = sum_product(integrands, dA, log)
  
          # Get a new, independently-scrambled QRNG for next time
          qrng = type(qrng)(seed=rngs[i], **qrng._init_quad)
  
      integral = mean(estimates, log)
      standard_error = sem(estimates, m=integral, log=log)
      integral = (
          integral + np.pi * 1j
          if (log and sign < 0)
          else integral * sign
      )
      return QMCQuadResult(integral, standard_error)
  
imports: from __future__ import annotations
  from typing import TYPE_CHECKING, Callable, Any, cast
  import numpy
  import numpy.typing
  import math
  import warnings
  from collections import namedtuple
  from scipy.special import roots_legendre
  from scipy.special import gammaln, logsumexp
  from scipy._lib._util import _rng_spawn
  from scipy._lib.deprecation import _NoValue, _deprecate_positional_args, _deprecated
label_desc: Use image node_id 26 for calling a global function with example usage: qmc_quad(func, a, b) and returns: QMCQuadResult, QMCQuadResult, unknown, np, unknown, np, np, np, unknown, unknown