def unc_wrapper_args(*covariance_keys):
    """
    Wrap function, calculate its Jacobian and calculate the covariance of the
    outputs given the covariance of the specified inputs.

    :param covariance_keys: indices and names of arguments corresponding to
        covariance
    :return: wrapped function bound to specified covariance keys

    This is the outer uncertainty wrapper that allows you to specify the
    arguments in the original function that correspond to the covariance. The
    inner wrapper takes the original function to be wrapped. ::

        def f(a, b, c, d, kw1='foo', *args, **kwargs):
            pass

        # arguments a, c, d and kw1 correspond to the covariance matrix
        f_wrapped = unc_wrapper_args(0, 2, 3, 'kw1')(f)

        cov = np.array([[0.0001, 0., 0., 0.], [0., 0.0001, 0., 0.],
                        [0., 0., 0.0001, 0.], [0., 0., 0., 0.0001])
        y, cov, jac = f_wrapped(a, b, c, d, kw1='bar', __covariance__=cov)

    The covariance keys can be indices of positional arguments or the names of
    keywords argument used in calling the function. If no covariance keys are
    specified then the arguments that correspond to the covariance shoud be
    grouped into a sequence. If ``None`` is anywhere in ``covariance_keys`` then
    all of the arguments will be used to calculate the Jacobian.

    The covariance matrix must be a symmetrical matrix with positive numbers on
    the diagonal that correspond to the square of the standard deviation, second
    moment around the mean or root-mean-square(RMS) of the function with respect
    to the arguments specified as covariance keys. The other elements are the
    covariances corresponding to the arguments intersecting at that element.
    Pass the covariance matrix with the keyword ``__covariance__`` and it will
    be popped from the dictionary of keyword arguments provided to the wrapped
    function.

    The wrapped function will return the evaluation of the original function,
    its Jacobian, which is the sensitivity of the return output to each
    argument specified as a covariance key and the covariance propagated using
    the first order terms of a Taylor series expansion around the arguments.

    An optional keyword argument ``__method__`` can also be passed to the
    wrapped function (not the wrapper) that specifies the method used to
    calculate the dot product. The default method is ``'loop'``. The other
    methods are ``'dense'``, ``'sparse'`` and ``'pool'``.

    If the arguments specified as covariance keys are arrays, they should all be
    the same size. These dimensions will be considered as separate observations.
    Another argument, not in the covariance keys, may also create observations.
    The resulting Jacobian will have dimensions of number of observations (nobs)
    by number of return output (nf) by number of covariance keys (nargs). The
    resulting covariance will be nobs x nf x nf.
    """
    def wrapper(f):
        @wraps(f)
        def wrapped_function(*args, **kwargs):
            cov = kwargs.pop('__covariance__', None)  # pop covariance
            method = kwargs.pop('__method__', 'loop')  # pop covariance
            # covariance keys cannot be defaults, they must be in args or kwargs
            cov_keys = covariance_keys
            # convert args to kwargs by index
            kwargs.update({n: v for n, v in enumerate(args)})
            args = ()  # empty args
            if None in cov_keys:
                # use all keys
                cov_keys = kwargs.keys()
            # group covariance keys
            if len(cov_keys) > 0:
                # uses specified keys
                x = [np.atleast_1d(kwargs.pop(k)) for k in cov_keys]
            else:
                # arguments already grouped
                x = kwargs.pop(0)  # use first argument
            # remaining args
            args_dict = {}

            def args_from_kwargs(kwargs_):
                """unpack positional arguments from keyword arguments"""
                # create mapping of positional arguments by index
                args_ = [(n, v) for n, v in kwargs_.iteritems()
                         if not isinstance(n, basestring)]
                # sort positional arguments by index
                idx, args_ = zip(*sorted(args_, key=lambda m: m[0]))
                # remove args_ and their indices from kwargs_
                args_dict_ = {n: kwargs_.pop(n) for n in idx}
                return args_, args_dict_

            if kwargs:
                args, args_dict = args_from_kwargs(kwargs)

            def f_(x_, *args_, **kwargs_):
                """call original function with independent variables grouped"""
                args_dict_ = args_dict
                if cov_keys:
                    kwargs_.update(zip(cov_keys, x_), **args_dict_)
                if kwargs_:
                    args_, _ = args_from_kwargs(kwargs_)
                    return np.array(f(*args_, **kwargs_))
                # assumes independent variables already grouped
                return f(x_, *args_, **kwargs_)

            # evaluate function and Jacobian
            avg = f_(x, *args, **kwargs)
            # number of returns and observations
            if avg.ndim > 1:
                nf, nobs = avg.shape
            else:
                nf, nobs = avg.size, 1
            jac = jacobian(f_, x, nf, nobs, *args, **kwargs)
            # calculate covariance
            if cov is not None:
                # covariance must account for all observations
                # scale covariances by x squared in each direction
                if cov.ndim == 3:
                    x = np.array([np.repeat(y, nobs) if len(y)==1
                                  else y for y in x])
                    LOGGER.debug('x:\n%r', x)
                    cov = np.array([c * y * np.row_stack(y)
                                    for c, y in zip(cov, x.T)])
                else: # x are all only one dimension
                    x = np.asarray(x)
                    cov = cov * x * x.T
                    assert jac.size / nf / nobs == cov.size / len(x)
                    cov = np.tile(cov, (nobs, 1, 1))
                # propagate uncertainty using different methods
                if method.lower() == 'dense':
                    j, c = jflatten(jac), jflatten(cov)
                    cov = prop_unc((j, c))
                # sparse
                elif method.lower() == 'sparse':
                    j, c = jtosparse(jac), jtosparse(cov)
                    cov = j.dot(c).dot(j.transpose())
                    cov = cov.todense()
                # pool
                elif method.lower() == 'pool':
                    try:
                        p = Pool()
                        cov = np.array(p.map(prop_unc, zip(jac, cov)))
                    finally:
                        p.terminate()
                # loop is the default
                else:
                    cov = np.array([prop_unc((jac[o], cov[o]))
                                    for o in xrange(nobs)])
                # dense and spares are flattened, unravel them into 3-D list of
                # observations
                if method.lower() in ['dense', 'sparse']:
                    cov = np.array([
                        cov[(nf * o):(nf * (o + 1)), (nf * o):(nf * (o + 1))]
                        for o in xrange(nobs)
                    ])
            # unpack returns for original function with ungrouped arguments
            if None in cov_keys or len(cov_keys) > 0:
                return tuple(avg.tolist() + [cov, jac])
            # independent variables were already grouped
            return avg, cov, jac
        return wrapped_function
    return wrapper