def bss_eval(reference_sources, estimated_sources,
             window=2 * 44100, hop=1.5 * 44100,
             compute_permutation=False,
             filters_len=512,
             framewise_filters=False,
             bsseval_sources_version=False
             ):
    """BSS_EVAL version 4.

    Measurement of the separation quality for estimated source signals
    in terms of source to distortion, interference and artifacts ratios,
    (SDR, SIR, SAR) as well as the image to spatial ratio (ISR), as defined
    in [#vincent2005bssevalv3]_.

    The metrics are computed on a framewise basis, with overlap allowed between
    the windows.

    The key difference between this version 4 and BSS Eval version 3 is the
    possibility of using the same distortion filters for all windows when
    matching the sources to their estimates, instead of estimating the filters
    anew at every frame, as done in BSS Eval v3.

    This implementation is fully compatible with BSS Eval v2 and v3 written
    in MATLAB.

    Examples
    --------
    >>> # reference_sources[n] should be a 2D ndarray, with first dimension the
    >>> # samples and second dimension the channels of the n'th reference
    >>> # source estimated_sources[n] should be the same for the n'th estimated
    >>> # source
    >>> (sdr, isr, sir, sar, perm) = mir_eval.separation.bss_eval(
    >>>    reference_sources,
    >>>    estimated_sources)

    Parameters
    ----------
    reference_sources : np.ndarray, shape=(nsrc, nsampl, nchan)
        matrix containing true sources
    estimated_sources : np.ndarray, shape=(nsrc, nsampl, nchan)
        matrix containing estimated sources
    window : int, optional
        size of each window for time-varying evaluation. Picking np.inf or any
        integer greater than nsampl will compute metrics on the whole signal.
    hop : int, optional
        hop size between windows
    compute_permutation : bool, optional
        compute all permutations of estimate/source combinations to compute
        the best scores (False by default). Note that picking True will lead
        to a significant computation overhead.
    filters_len : int, optional
        maximum time lag for the computation of the distortion filters. Default
        is filters_len = 512.
    framewise_filters : bool, optional
        Compute a new distortion filter for each frame (False by default). Note
        that picking True as in BSS Eval v2 and v3 leads to a significant
        computation overhead.
    bsseval_sources_version : bool, optional
        if  ``True``, results correspond to the `bss_eval_sources` version from
        the BSS Eval v2 and v3. Note however that this is not recommended
        because this evaluation method modifies the references according to the
        estimated sources, leading to potential problems for the estimation of
        SDR. For instance, zeroing some frequencies in the estimates will lead
        those to also be zeroed in the references, and hence not evaluated,
        artificially boosting results. For this reason, SiSEC always uses
        the `bss_eval_images` version, corresponding to ``False``.

    Returns
    -------
    sdr : np.ndarray, shape=(nsrc, nwin)
        matrix of Signal to Distortion Ratios (SDR). One for each source and
        window
    isr : np.ndarray, shape=(nsrc, nwin)
        matrix of source Image to Spatial distortion Ratios (ISR)
    sir : np.ndarray, shape=(nsrc, nwin)
        matrix of Source to Interference Ratios (SIR)
    sar : np.ndarray, shape=(nsrc, nwin)
        matrix of Sources to Artifacts Ratios (SAR)
    perm : np.ndarray, shape=(nsrc, nwin)
        vector containing the best ordering of estimated sources in
        the mean SIR sense (estimated source number ``perm[j,t]`` corresponds
        to true source number ``j`` at window ``t``).
        Note: ``perm`` will be ``(0,2,...,nsrc-1)`` if ``compute_permutation``
        is ``False``.

    References
    ----------
  .. [#liutkus2018bssevalv4] Antoine Liutkus, Fabian-Robert Stöter and Nobutaka
     Ito, "The 2018 Signal Separation Evaluation Campaign," In Proceedings of
     LVA/ICA 2018.
  .. [#vincent2005bssevalv3] Emmanuel Vincent, Rémi Gribonval, and Cédric
      Févotte, "Performance measurement in blind audio source separation," IEEE
      Trans. on Audio, Speech and Language Processing, 2006."""

    # assuming input is in shape (nsampl) or (nsrc, nsampl)
    estimated_sources = np.atleast_3d(estimated_sources)
    reference_sources = np.atleast_3d(reference_sources)

    # validate input
    validate(reference_sources, estimated_sources)

    # If empty matrices were supplied, return empty lists (special case)
    if reference_sources.size == 0 or estimated_sources.size == 0:
        return (
            np.array([]),
            np.array([]),
            np.array([]),
            np.array([]),
            np.array([])
        )

    # determine shape parameters
    (nsrc, nsampl, nchan) = estimated_sources.shape

    # defines all the permutations desired by user
    if compute_permutation:
        candidate_permutations = np.array(list(
            itertools.permutations(list(range(nsrc)))))
    else:
        candidate_permutations = np.array(np.arange(nsrc))[None, :]

    # initialize variables
    framer = Framing(window, hop, nsampl)
    nwin = framer.nwin

    (SDR, ISR, SIR, SAR) = list(range(4))
    s_r = np.empty((4, nsrc, nsrc, nwin))

    # define helper functions for computing filters on windows of the signals
    def compute_GsfC(win=slice(0, nsampl)):
        # First compute the references correlations
        G, sf = _compute_reference_correlations(
            reference_sources[:, win], filters_len
        )
        # compute the interference distortion filters
        C = np.zeros((nsrc, nsrc, nchan, filters_len, nchan))
        for jtrue in range(nsrc):
            C[jtrue] = _compute_projection_filters(
                G, sf, estimated_sources[jtrue, win]
            )
        return (G, sf, C)

    def compute_Cj(win=slice(0, nsampl)):
        Cj = np.zeros((nsrc, nsrc, 1, nchan, filters_len, nchan))
        for jtrue in range(nsrc):
            for jest in candidate_permutations[:, jtrue]:
                # compute the projection filters for this combination
                Cj[jtrue, jest] = _compute_projection_filters(
                    G[jtrue, jtrue],
                    sf[jtrue],
                    estimated_sources[jest, win]
                )
        return Cj

    if not framewise_filters:
        # compute filters on whole signals if no framewise filters
        (G, sf, C) = compute_GsfC()
        Cj = compute_Cj()

    # loop over all windows
    for (t, win) in enumerate(framer):
        # if we have time-varying distortion filters
        if framewise_filters:
            (G, sf, C) = compute_GsfC(win)
            Cj = compute_Cj(win)

        # loop over all permutations
        done = np.zeros((nsrc, nsrc))

        ref_slice = reference_sources[:, win]
        est_slice = estimated_sources[:, win]
        if (
            not _any_source_silent(ref_slice) and
            not _any_source_silent(est_slice)
        ):
            for jtrue in range(nsrc):
                for (k, jest) in enumerate(candidate_permutations[:, jtrue]):
                    # if we have a silent frame set results as np.nan
                    if not done[jtrue, jest]:
                            s_true, e_spat, e_interf, e_artif = \
                                _bss_decomp_mtifilt(
                                    reference_sources[:, win],
                                    estimated_sources[jest, win],
                                    jtrue, C[jest],
                                    Cj[jtrue, jest, 0]
                                )
                            s_r[:, jtrue, jest, t] = _bss_crit(
                                s_true,
                                e_spat,
                                e_interf,
                                e_artif,
                                bsseval_sources_version
                            )
                            done[jtrue, jest] = True
        else:
            a = np.empty((4, nsrc, nsrc))
            a[:] = np.nan
            s_r[:, :, :, t] = a

    # select the best ordering
    if framewise_filters:
        # if we have framewise filters, output one permutation for each window
        mean_sir = np.empty((len(candidate_permutations), nwin))
        axis_mean = 0
    else:
        # otherwise, output one permutation for the whole signal as the best
        # average one
        mean_sir = np.empty((len(candidate_permutations), 1))
        axis_mean = None
    dum = np.arange(nsrc)
    for (i, perm) in enumerate(candidate_permutations):
        mean_sir[i] = np.mean(s_r[SIR, dum, perm, :], axis=axis_mean)
    popt = candidate_permutations[np.argmax(mean_sir, axis=0)].T

    # now prepare the output
    if not framewise_filters:
        result = s_r[:, dum, popt[:, 0], :]
    else:
        result = np.empty((4, nsrc, nwin))
        for (m, t) in itertools.product(list(range(4)), list(range(nwin))):
            result[m, :, t] = s_r[m, dum, popt[:, t], t]

    return (result[SDR], result[ISR], result[SIR], result[SAR], popt)