r"""
    Apply the function(s) :math:`f(x;\theta)` to the arrays defined in `x` and
    `samples`.  Has options for weighting, trimming, cacheing & parallelising.

    Additionally, if a list of log-evidences are passed, along with list of
    functions, samples and optional weights it marginalises over the models
    according to the evidences.

    Parameters
    ----------
    f: function
        function :math:`f(x;\theta)` (or list of functions for each model) with
        dependent variable :math:`x`, parameterised by :math:`\theta`.

    x: 1D array-like
        `x` values to evaluate :math:`f(x;\theta)` at.

    samples: 2D array-like
        :math:`\theta` samples (or list of :math:`\theta` samples) to evaluate
        :math:`f(x;\theta)` at.
        `shape = (nsamples, npars)`

    logZ: 1D array-like, optional
        log-evidences of each model if multiple models are passed.
        Should be same length as the list `f`, and need not be normalised.
        Default: `numpy.ones_like(f)`

    weights: 1D array-like, optional
        sample weights (or list of weights), if desired. Should have length
        same as `samples.shape[0]`.
        Default: `numpy.ones_like(samples)`

    ntrim: int, optional
        Approximate number of samples to trim down to, if desired. Useful if
        the posterior is dramatically oversampled.
        Default: None

    cache: str, optional
        File root for saving previous calculations for re-use.
        Default: None

    parallel, tqdm_args:
        see docstring for :func:`fgivenx.parallel.parallel_apply`

    Returns
    -------
    2D numpy.array
        Evaluate the function `f` at each x value and each theta.
        Equivalent to `[[f(x_i,theta) for theta in samples] for x_i in x]`

def compute_samples(f, x, samples, **kwargs):
    r"""
    Apply the function(s) :math:`f(x;\theta)` to the arrays defined in `x` and
    `samples`.  Has options for weighting, trimming, cacheing & parallelising.

    Additionally, if a list of log-evidences are passed, along with list of
    functions, samples and optional weights it marginalises over the models
    according to the evidences.

    Parameters
    ----------
    f: function
        function :math:`f(x;\theta)` (or list of functions for each model) with
        dependent variable :math:`x`, parameterised by :math:`\theta`.

    x: 1D array-like
        `x` values to evaluate :math:`f(x;\theta)` at.

    samples: 2D array-like
        :math:`\theta` samples (or list of :math:`\theta` samples) to evaluate
        :math:`f(x;\theta)` at.
        `shape = (nsamples, npars)`

    logZ: 1D array-like, optional
        log-evidences of each model if multiple models are passed.
        Should be same length as the list `f`, and need not be normalised.
        Default: `numpy.ones_like(f)`

    weights: 1D array-like, optional
        sample weights (or list of weights), if desired. Should have length
        same as `samples.shape[0]`.
        Default: `numpy.ones_like(samples)`

    ntrim: int, optional
        Approximate number of samples to trim down to, if desired. Useful if
        the posterior is dramatically oversampled.
        Default: None

    cache: str, optional
        File root for saving previous calculations for re-use.
        Default: None

    parallel, tqdm_args:
        see docstring for :func:`fgivenx.parallel.parallel_apply`

    Returns
    -------
    2D numpy.array
        Evaluate the function `f` at each x value and each theta.
        Equivalent to `[[f(x_i,theta) for theta in samples] for x_i in x]`
    """
    logZ = kwargs.pop('logZ', None)
    weights = kwargs.pop('weights', None)
    ntrim = kwargs.pop('ntrim', None)
    cache = kwargs.pop('cache', '')
    parallel = kwargs.pop('parallel', False)
    tqdm_kwargs = kwargs.pop('tqdm_kwargs', {})
    if kwargs:
        raise TypeError('Unexpected **kwargs: %r' % kwargs)

    logZ, f, x, samples, weights = _check_args(logZ, f, x, samples, weights)

    logZ, weights = _normalise_weights(logZ, weights, ntrim)

    for i, (s, w) in enumerate(zip(samples, weights)):
        samples[i] = _equally_weight_samples(s, w)

    return fgivenx.samples.compute_samples(f, x, samples,
                                           parallel=parallel, cache=cache,
                                           tqdm_kwargs=tqdm_kwargs)