docs: def compute_dkl(f, x, samples, prior_samples, **kwargs):
      r"""
      Compute the Kullback-Leibler divergence at each value of `x` for the prior
      and posterior defined by `prior_samples` and `samples`.
  
      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, prior_samples: 2D array-like
          :math:`\theta` samples (or list of :math:`\theta` samples) from
          posterior and prior 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, prior_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, prior_cache: str, optional
          File roots for saving previous calculations for re-use
  
      parallel, tqdm_args:
          see docstring for :func:`fgivenx.parallel.parallel_apply`
  
      kwargs: further keyword arguments
          Any further keyword arguments are plotting keywords that are passed to
          :func:`fgivenx.plot.plot`.
  
      Returns
      -------
      1D numpy array:
          dkl values at each value of `x`.
      """
  
      logZ = kwargs.pop('logZ', None)
      weights = kwargs.pop('weights', None)
      prior_weights = kwargs.pop('prior_weights', None)
      ntrim = kwargs.pop('ntrim', None)
      cache = kwargs.pop('cache', '')
      prior_cache = kwargs.pop('prior_cache', '')
      parallel = kwargs.pop('parallel', False)
      tqdm_kwargs = kwargs.pop('tqdm_kwargs', {})
  
      if kwargs:
          raise TypeError('Unexpected **kwargs: %r' % kwargs)
  
      if logZ is None:
          logZ = [0]
          f = [f]
          samples = [samples]
          prior_samples = [prior_samples]
          weights = [weights]
          prior_weights = [prior_weights]
          cache = [cache]
          prior_cache = [prior_cache]
  
      DKLs = []
  
      for fi, c, pc, s, w, ps, pw in zip(f, cache, prior_cache, samples, weights,
                                         prior_samples, prior_weights):
  
          fsamps = compute_samples(fi, x, s, weights=w, ntrim=ntrim,
                                   parallel=parallel, cache=c,
                                   tqdm_kwargs=tqdm_kwargs)
  
          fsamps_prior = compute_samples(fi, x, ps, weights=pw, ntrim=ntrim,
                                         parallel=parallel, cache=pc,
                                         tqdm_kwargs=tqdm_kwargs)
  
          dkls = fgivenx.dkl.compute_dkl(fsamps, fsamps_prior,
                                         parallel=parallel, cache=c,
                                         tqdm_kwargs=tqdm_kwargs)
          DKLs.append(dkls)
  
      logZ = numpy.array(logZ)
      DKLs = numpy.array(DKLs)
  
      Zs = numpy.exp(logZ-logZ.max())
      Zs /= Zs.sum()
      return numpy.sum(Zs * DKLs.transpose(), axis=1)