def nquad(func, ranges, args=None, opts=None):
depth = len(ranges)<EOL>ranges = [rng if callable(rng) else _RangeFunc(rng) for rng in ranges]<EOL>if args is None:<EOL><INDENT>args = ()<EOL><DEDENT>if opts is None:<EOL><INDENT>opts = [dict([])] * depth<EOL><DEDENT>if isinstance(opts, dict):<EOL><INDENT>opts = [opts] * depth<EOL><DEDENT>else:<EOL><INDENT>opts = [opt if callable(opt) else _OptFunc(opt) for opt in opts]<EOL><DEDENT>return _NQuad(func, ranges, opts).integrate(*args)<EOL>
Integration over multiple variables.

Wraps `quad` to enable integration over multiple variables.
Various options allow improved integration of discontinuous functions, as
well as the use of weighted integration, and generally finer control of the
integration process.

Parameters
----------
func : callable
    The function to be integrated. Has arguments of ``x0, ... xn``,
    ``t0, tm``, where integration is carried out over ``x0, ... xn``, which
    must be floats.  Function signature should be
    ``func(x0, x1, ..., xn, t0, t1, ..., tm)``.  Integration is carried out
    in order.  That is, integration over ``x0`` is the innermost integral,
    and ``xn`` is the outermost.
    If performance is a concern, this function may be a ctypes function of
    the form::

        f(int n, double args[n])

    where ``n`` is the number of extra parameters and args is an array
    of doubles of the additional parameters.  This function may then
    be compiled to a dynamic/shared library then imported through
    ``ctypes``, setting the function's argtypes to ``(c_int, c_double)``,
    and the function's restype to ``(c_double)``.  Its pointer may then be
    passed into `nquad` normally.
    This allows the underlying Fortran library to evaluate the function in
    the innermost integration calls without callbacks to Python, and also 
    speeds up the evaluation of the function itself.
ranges : iterable object
    Each element of ranges may be either a sequence  of 2 numbers, or else
    a callable that returns such a sequence.  ``ranges[0]`` corresponds to
    integration over x0, and so on.  If an element of ranges is a callable,
    then it will be called with all of the integration arguments available.
    e.g. if ``func = f(x0, x1, x2)``, then ``ranges[0]`` may be defined as
    either ``(a, b)`` or else as ``(a, b) = range0(x1, x2)``.
args : iterable object, optional
    Additional arguments ``t0, ..., tn``, required by `func`.
opts : iterable object or dict, optional
    Options to be passed to `quad`.  May be empty, a dict, or
    a sequence of dicts or functions that return a dict.  If empty, the
    default options from scipy.integrate.quadare used.  If a dict, the same
    options are used for all levels of integraion.  If a sequence, then each
    element of the sequence corresponds to a particular integration. e.g.
    opts[0] corresponds to integration over x0, and so on. The available
    options together with their default values are:

      - epsabs = 1.49e-08
      - epsrel = 1.49e-08
      - limit  = 50
      - points = None
      - weight = None
      - wvar   = None
      - wopts  = None

    The ``full_output`` option from `quad` is unavailable, due to the
    complexity of handling the large amount of data such an option would
    return for this kind of nested integration.  For more information on
    these options, see `quad` and `quad_explain`.

Returns
-------
result : float
    The result of the integration.
abserr : float
    The maximum of the estimates of the absolute error in the various
    integration results.

See Also
--------
quad : 1-dimensional numerical integration
dblquad, tplquad : double and triple integrals
fixed_quad : fixed-order Gaussian quadrature
quadrature : adaptive Gaussian quadrature

Examples
--------
>>> from scipy import integrate
>>> func = lambda x0,x1,x2,x3 : x0**2 + x1*x2 - x3**3 + np.sin(x0) + (
...                                 1 if (x0-.2*x3-.5-.25*x1>0) else 0)
>>> points = [[lambda (x1,x2,x3) : 0.2*x3 + 0.5 + 0.25*x1], [], [], []]
>>> def opts0(*args, **kwargs):
...     return {'points':[0.2*args[2] + 0.5 + 0.25*args[0]]}
>>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]],
...                 opts=[opts0,{},{},{}])
(1.5267454070738633, 2.9437360001402324e-14)

>>> scale = .1
>>> def func2(x0, x1, x2, x3, t0, t1):
...     return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0+t1*x1-t0>0 else 0)
>>> def lim0(x1, x2, x3, t0, t1):
...     return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
...             scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]
>>> def lim1(x2, x3, t0, t1):
...     return [scale * (t0*x2 + t1*x3) - 1,
...             scale * (t0*x2 + t1*x3) + 1]
>>> def lim2(x3, t0, t1):
...     return [scale * (x3 + t0**2*t1**3) - 1,
...             scale * (x3 + t0**2*t1**3) + 1]
>>> def lim3(t0, t1):
...     return [scale * (t0+t1) - 1, scale * (t0+t1) + 1]
>>> def opts0(x1, x2, x3, t0, t1):
...     return {'points' : [t0 - t1*x1]}
>>> def opts1(x2, x3, t0, t1):
...     return {}
>>> def opts2(x3, t0, t1):
...     return {}
>>> def opts3(t0, t1):
...     return {}
>>> integrate.nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0),
                    opts=[opts0, opts1, opts2, opts3])
(25.066666666666666, 2.7829590483937256e-13)