def dpp_sampler_generic_kernel(K, random_state=None):
    """ Sample from generic :math:`\\operatorname{DPP}(\\mathbf{K})` with potentially non hermitian correlation kernel :math:`\\operatorname{DPP}(\\mathbf{K})` based on :math:`LU` factorization procedure.

    :param K:
        Correlation kernel (potentially non hermitian).
    :type K:
        array_like

    :return:
        A sample :math:`\\mathcal{X}` from :math:`\\operatorname{DPP}(K)` and
        the in-place :math:`LU factorization of :math:`K − I_{\\mathcal{X}^{c}}` where :math:`I_{\\mathcal{X}^{c}}` is the diagonal indicator matrix for the entries not in the sample :math:`\\mathcal{X}`.
    :rtype:
        list and array_like

    .. seealso::

        - :cite:`Pou19` Algorithm 1
    """

    rng = check_random_state(random_state)

    A = K.copy()
    sample = []

    for j in range(len(A)):

        if rng.rand() < A[j, j]:
            sample.append(j)
        else:
            A[j, j] -= 1

        A[j + 1:, j] /= A[j, j]
        A[j + 1:, j + 1:] -= np.outer(A[j + 1:, j], A[j, j + 1:])
        # A[j+1:, j+1:] -=  np.einsum('i,j', A[j+1:, j], A[j, j+1:])

    return sample, A def proj_dpp_sampler_kernel(kernel, mode='GS', size=None, random_state=None):
    """
    .. seealso::

        - :func:`proj_dpp_sampler_kernel_GS <proj_dpp_sampler_kernel_GS>`
        - :func:`proj_dpp_sampler_kernel_Schur <proj_dpp_sampler_kernel_Schur>`
        - :func:`proj_dpp_sampler_kernel_Chol <proj_dpp_sampler_kernel_Chol>`
    """

    rng = check_random_state(random_state)

    if size:
        rank = np.rint(np.trace(kernel)).astype(int)
        if size > rank:
            raise ValueError('size k={} > rank={}'. format(size, rank))

    # Sample from orthogonal projection kernel K = K^2 = K.H K
    if mode == 'GS':  # Gram-Schmidt equiv Cholesky
        sampl = proj_dpp_sampler_kernel_GS(kernel, size, rng)

    elif mode == 'Chol':  # Cholesky updates of Pou19
        sampl = proj_dpp_sampler_kernel_Chol(kernel, size, rng)[0]

    elif mode == 'Schur':  # Schur complement
        sampl = proj_dpp_sampler_kernel_Schur(kernel, size, rng)

    else:
        str_list = ['Invalid sampling mode, choose among:',
                    '- "GS (default)',
                    '- "Chol"',
                    '- "Schur"',
                    'Given "{}"'.format(mode)]
        raise ValueError('\n'.join(str_list))

    return sampl def proj_dpp_sampler_eig(eig_vecs, mode='GS', size=None,
                         random_state=None):
    """ Sample from projection :math:`\\operatorname{DPP}(K)` using the eigendecomposition of the projection kernel :math:`K=VV^{\\top}` where :math:`V^{\\top}V = I_r` and :math:`r=\\operatorname{rank}(\\mathbf{K})`.

    .. seealso::

        Phase 1:

        - :func:`dpp_eig_vecs_selector <dpp_eig_vecs_selector>`

        Phase 2:

        - :func:`proj_dpp_sampler_eig_GS <proj_dpp_sampler_eig_GS>`
        - :func:`proj_dpp_sampler_eig_GS_bis <proj_dpp_sampler_eig_GS_bis>`
        - :func:`proj_dpp_sampler_eig_KuTa12 <proj_dpp_sampler_eig_KuTa12>`
    """

    rng = check_random_state(random_state)

    if eig_vecs.shape[1]:
        # Phase 2: Sample from projection kernel VV.T
        # Chain rule, conditionals are updated using:
        if mode == 'GS':  # Gram-Schmidt
            sampl = proj_dpp_sampler_eig_GS(eig_vecs, size, rng)

        elif mode == 'GS_bis':  # Slight modif of 'GS'
            sampl = proj_dpp_sampler_eig_GS_bis(eig_vecs, size, rng)

        elif mode == 'KuTa12':  # cf Kulesza-Taskar
            sampl = proj_dpp_sampler_eig_KuTa12(eig_vecs, size, rng)

        else:
            str_list = ['Invalid sampling mode, choose among:',
                        '- "GS" (default)',
                        '- "GS_bis"',
                        '- "KuTa12"',
                        'Given "{}"'.format(mode)]
            raise ValueError('\n'.join(str_list))
    else:
        sampl = []

    return sampl def dpp_vfx_sampler(intermediate_sample_info,
                    X_data,
                    eval_L,
                    random_state=None,
                    **params):
    """ First pre-compute quantities necessary for the vfx rejection sampling loop, such as the inner Nystrom approximation, and the RLS of all elements in :math:`\\mathbf{L}`.
    Then, given the pre-computed information,run a rejection sampling loop to generate DPP samples.

    :param intermediate_sample_info:
        If available, the pre-computed information necessary for the vfx rejection sampling loop.
        If ``None``, this function will compute and return an ``_IntermediateSampleInfo`` with fields

        - ``.alpha_star``: appropriate rescaling such that the expected sample size of :math:`\\operatorname{DPP}(\\alpha^* \\mathbf{L})` is equal to a user-indicated constant ``params['desired_expected_size']``, or 1.0 if no such constant was specified by the user.
        - ``.logdet_I_A``: :math:`\\log \\det` of the Nystrom approximation of :math:`\\mathbf{L} + I`
        - ``.q``: placeholder q constant used for vfx sampling, to be replaced by the user before the sampling loop
        - ``.s`` and ``.z``: approximations of the expected sample size of :math:`\\operatorname{DPP}(\\alpha^* \\mathbf{L})` to be used in the sampling loop. For more details see :cite:`DeCaVa19`
        - ``.rls_estimate``: approximations of the RLS of all elements in X (i.e. in :math:`\\mathbf{L}`)

    :type intermediate_sample_info:
        ``_IntermediateSampleInfo`` or ``None``, default ``None``

    :param array_like X_data:
        dataset such that :math:`\\mathbf{L}=` ``eval_L(X_data)``, out of which we are sampling objects according to a DPP

    :param callable eval_L:
        Likelihood function.
        Given two sets of n points X and m points Y, ``eval_L(X, Y)`` should compute the :math:`n x m` matrix containing the likelihood between points.
        The function should also accept a single argument X and return ``eval_L(X) = eval_L(X, X)``.
        As an example, see the implementation of any of the kernels provided by scikit-learn (e.g. `PairwiseKernel <https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.PairwiseKernel.html>`_).

    :param random_state:
        random source used for sampling, if None a RandomState is automatically generated

    :type random_state:
        RandomState or None, default None

    :param dict params:
        Dictionary including optional parameters:

        - ``'desired_expected_size'`` (float or None, default None)

            Desired expected sample size for the DPP.
            If None, use the natural DPP expected sample size.
            The vfx sampling algorithm can approximately adjust the expected sample size of the DPP by rescaling the :math:`\\mathbf{L}` matrix with a scalar :math:`\\alpha^*\\leq 1` .
            Adjusting the expected sample size can be useful to control downstream complexity, and it is necessary to improve the probability of drawing a sample with exactly :math:`k` elements when using vfx for k-DPP sampling.
            Currently only reducing the sample size is supported, and the sampler will return an exception if the DPP sample has already a natural expected size smaller than ``params['desired_expected_size'``.]

        - ``'rls_oversample_dppvfx'`` (float, default 4.0)

            Oversampling parameter used to construct dppvfx's internal Nystrom approximation.
            The ``rls_oversample_dppvfx``:math:`\\geq 1` parameter is used to increase the rank of the approximation by a ``rls_oversample_dppvfx`` factor.
            This makes each rejection round slower and more memory intensive, but reduces variance and the number of rounds of rejections, so the actual runtime might increase or decrease.
            Empirically, a small factor ``rls_oversample_dppvfx``:math:`\\in [2,10]` seems to work.
            It is suggested to start with a small number and increase if the algorithm fails to terminate.

        - ``'rls_oversample_bless'`` (float, default 4.0)
            Oversampling parameter used during bless's internal Nystrom approximation.
            Note that this is a different Nystrom approximation than the one related to :func:`rls_oversample_dppvfx`, and can be tuned separately.
            The ``rls_oversample_bless``:math:`\\geq 1` parameter is used to increase the rank of the approximation by a ``rls_oversample_bless`` factor.
            This makes the one-time pre-processing slower and more memory intensive, but reduces variance and the number of rounds of rejections, so the actual runtime might increase or decrease.
            Empirically, a small factor ``rls_oversample_bless``:math:`\\in [2,10]` seems to work.
            It is suggested to start with a small number and increase if the algorithm fails to terminate or is not accurate.

        - ``'q_func'`` (function, default x: x*x)
            Mapping from estimate expected size of the DPP to Poisson intensity used to choose size of the intermediate sample.
            Larger intermediate sampler cause less efficient iterations but higher acceptance probability.

        - ``'nb_iter_bless'`` (int or None, default None)

            Iterations for inner BLESS execution, if None it is set to log(n)

        - ``'verbose'`` (bool, default True)

            Controls verbosity of debug output, including progress bars.
            If intermediate_sample_info is not provided, the first progress bar reports the inner execution of
            the bless algorithm, showing:

                - lam: lambda value of the current iteration
                - m: current size of the dictionary (number of centers contained)
                - m_expected: expected size of the dictionary before sampling
                - probs_dist: (mean, max, min) of the approximate rlss at the current iteration

            Subsequent progress bars show the execution of each rejection sampling loops (i.e. once per sample generated)
                - acc_thresh: latest computed probability of acceptance
                - rej_iter: iteration of the rejection sampling loop (i.e. rejections so far)

        - ``'max_iter'`` (int, default 1000)

            Maximum number of intermediate sample rejections before giving up.

    :return:
        Sample from a DPP (as a list) and updated intermediate_sample_info

    :rtype:
        tuple(list, _IntermediateSampleInfo)
    """
    rng = check_random_state(random_state)

    if intermediate_sample_info is None:
        intermediate_sample_info = vfx_sampling_precompute_constants(
                                    X_data=X_data,
                                    eval_L=eval_L,
                                    rng=rng,
                                    **params)

        q_func = params.get('q_func', lambda s: s * s)
        intermediate_sample_info = intermediate_sample_info._replace(q=q_func(intermediate_sample_info.s))

    sampl, rej_count = vfx_sampling_do_sampling_loop(X_data, eval_L, intermediate_sample_info, rng, **params)

    return sampl, intermediate_sample_info def alpha_dpp_sampler(intermediate_sample_info,
                      X_data,
                      eval_L,
                      random_state=None,
                      **params):
    """ First pre-compute quantities necessary for the alpha-dpp rejection sampling loop, such as the inner Nystrom
    approximation, and the and the initial rescaling alpha_hat for the binary search.
    Then, given the pre-computed information,run a rejection sampling loop to generate samples from DPP(alpha * L).

    :param intermediate_sample_info:
        If available, the pre-computed information necessary for the alpha-dpp rejection sampling loop.
        If ``None``, this function will compute and return an ``_IntermediateSampleInfoAlphaRescale`` (see :func:`alpha_dpp_sampling_precompute_constants`)

    :type intermediate_sample_info:
        ``_IntermediateSampleInfoAlphaRescale`` or ``None``, default ``None``

    :param array_like X_data:
        dataset such that :math:`\\mathbf{L}=` ``eval_L(X_data)``, out of which we are sampling objects according to a DPP

    :param callable eval_L:
        Likelihood function.
        Given two sets of n points X and m points Y, ``eval_L(X, Y)`` should compute the :math:`n x m` matrix containing the likelihood between points.
        The function should also accept a single argument X and return ``eval_L(X) = eval_L(X, X)``.
        As an example, see the implementation of any of the kernels provided by scikit-learn (e.g. `PairwiseKernel <https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.PairwiseKernel.html>`_).

    :param random_state:
        random source used for sampling, if None a RandomState is automatically generated

    :type random_state:
        RandomState or None, default None

    :param dict params:
        Dictionary including optional parameters:

        - ``'desired_expected_size'`` (float or None, default None)

            Desired expected sample size for the rescaled DPP.
            If None, use the natural DPP expected sample size.
            The alpha-dpp sampling algorithm can approximately adjust the expected sample size of the DPP by rescaling the :math:`\\mathbf{L}` matrix with a scalar :math:`\\alpha^*\\leq 1` .
            Adjusting the expected sample size can be useful to control downstream complexity, and it is necessary to improve the probability of drawing a sample with exactly :math:`k` elements when using alpha-dpp for k-DPP sampling.
            Currently only reducing the sample size is supported, and the sampler will return an exception if the DPP sample has already a natural expected size smaller than ``params['desired_expected_size'``.]

        - ``'rls_oversample_alphadpp'`` (float, default 4.0)

            Oversampling parameter used to construct alphadpp's internal Nystrom approximation.
            The ``rls_oversample_alphadpp``:math:`\\geq 1` parameter is used to increase the rank of the approximation by a ``rls_oversample_alphadpp`` factor.
            This makes each rejection round slower and more memory intensive, but reduces variance and the number of rounds of rejections, so the actual runtime might increase or decrease.
            Empirically, a small factor ``rls_oversample_alphadpp``:math:`\\in [2,10]` seems to work.
            It is suggested to start with a small number and increase if the algorithm fails to terminate.

        - ``'rls_oversample_bless'`` (float, default 4.0)
            Oversampling parameter used during bless's internal Nystrom approximation.
            Note that this is a different Nystrom approximation than the one related to :func:`rls_oversample_alphadpp`, and can be tuned separately.
            The ``rls_oversample_bless``:math:`\\geq 1` parameter is used to increase the rank of the approximation by a ``rls_oversample_bless`` factor.
            This makes the one-time pre-processing slower and more memory intensive, but reduces variance and the number of rounds of rejections, so the actual runtime might increase or decrease.
            Empirically, a small factor ``rls_oversample_bless``:math:`\\in [2,10]` seems to work.
            It is suggested to start with a small number and increase if the algorithm fails to terminate or is not accurate.

        - ``'r_func'`` (function, default x: x)
            Mapping from estimate expected size of the rescaled alpha-DPP to Poisson intensity used to choose size
            of the intermediate sample. Larger intermediate sampler cause less efficient iterations but higher
            acceptance probability.

        - ``'nb_iter_bless'`` (int or None, default None)

            Iterations for inner BLESS execution, if None it is set to log(n)

        - ``'verbose'`` (bool, default True)

            Controls verbosity of debug output, including progress bars.
            If intermediate_sample_info is not provided, the first progress bar reports the inner execution of
            the bless algorithm, showing:

                - lam: lambda value of the current iteration
                - m: current size of the dictionary (number of centers contained)
                - m_expected: expected size of the dictionary before sampling
                - probs_dist: (mean, max, min) of the approximate rlss at the current iteration

            Subsequent progress bars show the execution of each rejection sampling loops (i.e. once per sample generated)
                - acc_thresh: latest computed probability of acceptance
                - rej_iter: iteration of the rejection sampling loop (i.e. rejections so far)

        - ``'max_iter'`` (int, default 1000)

            Maximum number of intermediate sample rejections before giving up.

    :return:
        Sample from a DPP (as a list) and updated intermediate_sample_info

    :rtype:
        tuple(list, _IntermediateSampleInfoAlphaRescale)
    """
    rng = check_random_state(random_state)

    if intermediate_sample_info is None:
        intermediate_sample_info = alpha_dpp_sampling_precompute_constants(
                                    X_data=X_data,
                                    eval_L=eval_L,
                                    rng=rng,
                                    **params)

        r_func = params.get('r_func', lambda r: r)
        intermediate_sample_info = intermediate_sample_info._replace(r=r_func(intermediate_sample_info.deff_alpha_L_hat))

    sampl, rej_count, intermediate_sample_info = alpha_dpp_sampling_do_sampling_loop(X_data,
                                            eval_L,
                                            intermediate_sample_info,
                                            rng,
                                            **params)

    return sampl, intermediate_sample_info def dpp_eig_vecs_selector(ber_params, eig_vecs,
                          random_state=None):
    """ Phase 1 of exact sampling procedure. Subsample eigenvectors :math:`V` of the initial kernel (correlation :math:`K`, resp. likelihood :math:`L`) to build a projection DPP with kernel :math:`U U^{\\top}` from which sampling is easy.
    The selection is made based on a realization of Bernoulli variables with parameters to the eigenvalues of :math:`K`.

    :param ber_params:
        Parameters of Bernoulli variables
        :math:`\\lambda^K=\\lambda^L/(1+\\lambda^L)
    :type ber_params:
        list, array_like

    :param eig_vecs:
        Collection of eigenvectors of the kernel :math:`K`, resp. :math:`L`
    :type eig_vecs:
        array_like

    :return:
        selected eigenvectors
    :rtype:
        array_like

    .. seealso::

        - :func:`dpp_sampler_eig <dpp_sampler_eig>`
    """
    rng = check_random_state(random_state)

    # Realisation of Bernoulli random variables with params ber_params
    ind_sel = rng.rand(ber_params.size) < ber_params

    return eig_vecs[:, ind_sel] def k_dpp_vfx_sampler(size,
                      intermediate_sample_info,
                      X_data,
                      eval_L,
                      random_state=None,
                      **params):
    """ First pre-compute quantities necessary for the vfx rejection sampling loop, such as the inner Nystrom approximation, and the RLS of all elements in :math:`\\mathbf{L}`.
    Then, given the pre-computed information,run a rejection sampling loop to generate DPP samples.
    To guarantee that the returned sample has size ``size``, we internally set desired_expected_size=size and
    then repeatedly invoke dpp_vfx_sampler until a sample of the correct size is returned,
    or exit with an error after a chosen number of rejections is reached.

    :param int size: The size of the sample (i.e. the k of k-DPPs)

    :param intermediate_sample_info:
        If available, the pre-computed information necessary for the vfx rejection sampling loop.
        If ``None``, this function will compute and return an ``_IntermediateSampleInfo`` with fields

        - ``.alpha_star``: appropriate rescaling such that the expected sample size of :math:`\\operatorname{DPP}(\\alpha^* \\mathbf{L})` is equal to a user-indicated constant ``params['desired_expected_size']``, or 1.0 if no such constant was specified by the user.
        - ``.logdet_I_A``: :math:`\\log \\det` of the Nystrom approximation of :math:`\\mathbf{L} + I`
        - ``.q``: placeholder q constant used for vfx sampling, to be replaced by the user before the sampling loop
        - ``.s`` and ``.z``: approximations of the expected sample size of :math:`\\operatorname{DPP}(\\alpha^* \\mathbf{L})` to be used in the sampling loop. For more details see :cite:`DeCaVa19`
        - ``.rls_estimate``: approximations of the RLS of all elements in X (i.e. in :math:`\\mathbf{L}`)

    :type intermediate_sample_info:
        ``_IntermediateSampleInfo`` or ``None``, default ``None``

    :param array_like X_data:
        dataset such that :math:`\\mathbf{L}=` ``eval_L(X_data)``, out of which we are sampling objects according to a DPP

    :param callable eval_L:
        Likelihood function.
        Given two sets of n points X and m points Y, ``eval_L(X, Y)`` should compute the :math:`n x m` matrix containing the likelihood between points.
        The function should also accept a single argument X and return ``eval_L(X) = eval_L(X, X)``.
        As an example, see the implementation of any of the kernels provided by scikit-learn (e.g. `PairwiseKernel <https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.PairwiseKernel.html>`_).

    :param random_state:
        random source used for sampling, if None a RandomState is automatically generated

    :type random_state:
        RandomState or None, default None

    :param dict params:
        Dictionary including optional parameters:

        - ``'rls_oversample_dppvfx'`` (float, default 4.0)

            Oversampling parameter used to construct dppvfx's internal Nystrom approximation.
            The ``rls_oversample_dppvfx``:math:`\\geq 1` parameter is used to increase the rank of the approximation by a ``rls_oversample_dppvfx`` factor.
            This makes each rejection round slower and more memory intensive, but reduces variance and the number of rounds of rejections, so the actual runtime might increase or decrease.
            Empirically, a small factor ``rls_oversample_dppvfx``:math:`\\in [2,10]` seems to work.
            It is suggested to start with a small number and increase if the algorithm fails to terminate.

        - ``'rls_oversample_bless'`` (float, default 4.0)
            Oversampling parameter used during bless's internal Nystrom approximation.
            Note that this is a different Nystrom approximation than the one related to :func:`rls_oversample_dppvfx`, and can be tuned separately.
            The ``rls_oversample_bless``:math:`\\geq 1` parameter is used to increase the rank of the approximation by a ``rls_oversample_bless`` factor.
            This makes the one-time pre-processing slower and more memory intensive, but reduces variance and the number of rounds of rejections, so the actual runtime might increase or decrease.
            Empirically, a small factor ``rls_oversample_bless``:math:`\\in [2,10]` seems to work.
            It is suggested to start with a small number and increase if the algorithm fails to terminate or is not accurate.

        - ``'q_func'`` (function, default x: x*x)
            Mapping from estimate expected size of the DPP to Poisson intensity used to choose size of the intermediate sample.
            Larger intermediate sampler cause less efficient iterations but higher acceptance probability.

        - ``'nb_iter_bless'`` (int or None, default None)

            Iterations for inner BLESS execution, if None it is set to log(n)

        - ``'verbose'`` (bool, default True)

            Controls verbosity of debug output, including progress bars.
            If intermediate_sample_info is not provided, the first progress bar reports the inner execution of
            the bless algorithm, showing:

                - lam: lambda value of the current iteration
                - m: current size of the dictionary (number of centers contained)
                - m_expected: expected size of the dictionary before sampling
                - probs_dist: (mean, max, min) of the approximate rlss at the current iteration

            Subsequent progress bars show the execution of each rejection sampling loops (i.e. once per sample generated)
                - acc_thresh: latest computed probability of acceptance
                - rej_iter: iteration of the rejection sampling loop (i.e. rejections so far)

        - ``'max_iter'`` (int, default 1000)

            Maximum number of intermediate sample rejections before giving up.

        - ``'max_iter_size_rejection'`` (int, default 100)

            Maximum number of size-based rejections before giving up.


    :return:
        Sample from a DPP (as a list) and updated intermediate_sample_info

    :rtype:
        tuple(list, _IntermediateSampleInfo)
    """
    rng = check_random_state(random_state)

    if (intermediate_sample_info is None
        or not np.isclose(intermediate_sample_info.s, size).item()):
        intermediate_sample_info = vfx_sampling_precompute_constants(
                                    X_data=X_data,
                                    eval_L=eval_L,
                                    desired_expected_size=size,
                                    rng=rng,
                                    **params)

        q_func = params.get('q_func', lambda s: s * s)

        intermediate_sample_info = intermediate_sample_info._replace(q=q_func(intermediate_sample_info.s))

    max_iter_size_rejection = params.get('max_iter_size_rejection', 100)

    for size_rejection_iter in range(max_iter_size_rejection):
        sampl, rej_count = vfx_sampling_do_sampling_loop(
                                X_data,
                                eval_L,
                                intermediate_sample_info,
                                rng,
                                **params)

        intermediate_sample_info = intermediate_sample_info._replace(rej_to_first_sample=intermediate_sample_info.rej_to_first_sample + rej_count)
        if len(sampl) == size:
            break
    else:
        raise ValueError('The vfx sampler reached the maximum number of rejections allowed '
                         'for the k-DPP size rejection ({}), try to increase the q factor '
                         '(see q_func parameter) or the Nystrom approximation accuracy '
                         'see rls_oversample_* parameters).'.format(max_iter_size_rejection))

    return sampl, intermediate_sample_info def alpha_k_dpp_sampler(size,
                        intermediate_sample_info,
                        X_data,
                        eval_L,
                        random_state=None,
                        **params):
    """ First pre-compute quantities necessary for the alpha-dpp rejection sampling loop, such as the inner Nystrom
    approximation, the and the initial rescaling alpha_hat for the binary search.
    Then, given the pre-computed information,run a rejection sampling loop to generate k-DPP samples.
    To guarantee that the returned sample has size ``size``, we internally set desired_expected_size=size and
    then repeatedly invoke alpha_dpp_sampler until a sample of the correct size is returned,
    or exit with an error after a chosen number of rejections is reached.

    :param int size: The size of the sample (i.e. the k of k-DPPs)

    :param intermediate_sample_info:
        If available, the pre-computed information necessary for the alpha-dpp rejection sampling loop.
        If ``None``, this function will compute and return an ``_IntermediateSampleInfoAlphaRescale`` (see :func:`alpha_dpp_sampling_precompute_constants`)

    :type intermediate_sample_info:
        ``_IntermediateSampleInfoAlphaRescale`` or ``None``, default ``None``

    :param array_like X_data:
        dataset such that :math:`\\mathbf{L}=` ``eval_L(X_data)``, out of which we are sampling objects according to a DPP

    :param callable eval_L:
        Likelihood function.
        Given two sets of n points X and m points Y, ``eval_L(X, Y)`` should compute the :math:`n x m` matrix containing the likelihood between points.
        The function should also accept a single argument X and return ``eval_L(X) = eval_L(X, X)``.
        As an example, see the implementation of any of the kernels provided by scikit-learn (e.g. `PairwiseKernel <https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.PairwiseKernel.html>`_).

    :param random_state:
        random source used for sampling, if None a RandomState is automatically generated

    :type random_state:
        RandomState or None, default None

    :param dict params:
        Dictionary including optional parameters:

        - ``'rls_oversample_alphadpp'`` (float, default 4.0)

            Oversampling parameter used to construct alphadpp's internal Nystrom approximation.
            The ``rls_oversample_alphadpp``:math:`\\geq 1` parameter is used to increase the rank of the approximation by a ``rls_oversample_alphadpp`` factor.
            This makes each rejection round slower and more memory intensive, but reduces variance and the number of rounds of rejections, so the actual runtime might increase or decrease.
            Empirically, a small factor ``rls_oversample_alphadpp``:math:`\\in [2,10]` seems to work.
            It is suggested to start with a small number and increase if the algorithm fails to terminate.

        - ``'rls_oversample_bless'`` (float, default 4.0)
            Oversampling parameter used during bless's internal Nystrom approximation.
            Note that this is a different Nystrom approximation than the one related to :func:`rls_oversample_alphadpp`, and can be tuned separately.
            The ``rls_oversample_bless``:math:`\\geq 1` parameter is used to increase the rank of the approximation by a ``rls_oversample_bless`` factor.
            This makes the one-time pre-processing slower and more memory intensive, but reduces variance and the number of rounds of rejections, so the actual runtime might increase or decrease.
            Empirically, a small factor ``rls_oversample_bless``:math:`\\in [2,10]` seems to work.
            It is suggested to start with a small number and increase if the algorithm fails to terminate or is not accurate.

        - ``'r_func'`` (function, default x: x)
            Mapping from estimate expected size of the rescaled alpha-DPP to Poisson intensity used to choose size
            of the intermediate sample. Larger intermediate sampler cause less efficient iterations but higher
            acceptance probability.

        - ``'nb_iter_bless'`` (int or None, default None)

            Iterations for inner BLESS execution, if None it is set to log(n)

        - ``'verbose'`` (bool, default True)

            Controls verbosity of debug output, including progress bars.
            If intermediate_sample_info is not provided, the first progress bar reports the inner execution of
            the bless algorithm, showing:

                - lam: lambda value of the current iteration
                - m: current size of the dictionary (number of centers contained)
                - m_expected: expected size of the dictionary before sampling
                - probs_dist: (mean, max, min) of the approximate rlss at the current iteration

            Subsequent progress bars show the execution of each rejection sampling loops (i.e. once per sample generated)
                - acc_thresh: latest computed probability of acceptance
                - rej_iter: iteration of the rejection sampling loop (i.e. rejections so far)

        - ``'early_stop'`` (bool, default False)

            Wheter to return as soon as a first sample is accepted. If True, the sampling loop is interrupted
            as soon as a k-DPP sample is generated. If False, the algorithm continues the binary search until
            of a sufficiently good rescaling alpha is found. While this makes subsequent sampling faster, it is wasteful
            in the case where a single k-DPP sample is desired.

        - ``'max_iter_size_rejection'`` (int, default 100)

            Maximum number of size-based rejections before giving up.

        - ``'max_iter_size_rejection'`` (int, default 100)

            Maximum number of size-based rejections before giving up.


    :return:
        Sample from a DPP (as a list) and updated intermediate_sample_info

    :rtype:
        tuple(list, _IntermediateSampleInfoAlphaRescale)
    """
    rng = check_random_state(random_state)

    if intermediate_sample_info is None or intermediate_sample_info.k != size:
        intermediate_sample_info = alpha_dpp_sampling_precompute_constants(
                                        X_data=X_data,
                                        eval_L=eval_L,
                                        desired_expected_size=size,
                                        rng=rng,
                                        **params)

        r_func = params.get('r_func', lambda r: r)

        intermediate_sample_info = intermediate_sample_info._replace(r=r_func(intermediate_sample_info.deff_alpha_L_hat))

    max_iter_size_rejection = params.get('max_iter_size_rejection', 100)
    number_trial_search = np.ceil(np.sqrt(size)).astype('int')
    stopping_ratio = (1 + 1 / (size + 3) ** 2)

    sample_count = 0
    trial_count = 0
    under_k_count = 0
    over_k_count = 0

    ratio_alpha = intermediate_sample_info.alpha_max / intermediate_sample_info.alpha_min
    found_good_alpha = ratio_alpha <= stopping_ratio

    prog_bar = get_progress_bar(disable=not params.get('verbose', False))
    verbose_outer = None
    if 'verbose' in params:
        verbose_outer = params.pop('verbose')
    params['verbose'] = False

    early_stop = params.get('early_stop', False)

    trial_count_overall = 0
    for size_rejection_iter in range(max_iter_size_rejection):
        sampl, rej_count, intermediate_sample_info = alpha_dpp_sampling_do_sampling_loop(X_data,
                                                eval_L,
                                                intermediate_sample_info,
                                                rng,
                                                **params)

        trial_count += 1
        trial_count_overall += 1

        prog_bar.set_postfix(trial_count=trial_count,
                             alpha="{:.4}".format(intermediate_sample_info.alpha_hat),
                             alpha_switch=intermediate_sample_info.alpha_switches,
                             k=size,
                             k_emp=len(sampl),
                             rej_count=rej_count)
        prog_bar.update()

        if len(sampl) == size:
            sampl_out = sampl
            if intermediate_sample_info.trial_to_first_sample == 0:
                intermediate_sample_info = intermediate_sample_info._replace(trial_to_first_sample=trial_count_overall)
            sample_count += 1
            if early_stop:
                break
        if len(sampl) < size:
            under_k_count += 1
        if len(sampl) > size:
            over_k_count += 1

        if intermediate_sample_info.trial_to_first_sample == 0:
            intermediate_sample_info = intermediate_sample_info._replace(rej_to_first_sample=intermediate_sample_info.rej_to_first_sample + rej_count)

        if sample_count == 2:
            found_good_alpha = True
            break

        if trial_count == number_trial_search:
            if under_k_count > over_k_count:
                intermediate_sample_info = intermediate_sample_info._replace(alpha_min=intermediate_sample_info.alpha_hat)
            else:
                intermediate_sample_info = intermediate_sample_info._replace(alpha_max=intermediate_sample_info.alpha_hat)

            geom_mean_alpha = np.sqrt(intermediate_sample_info.alpha_min * intermediate_sample_info.alpha_max)
            diag_L = intermediate_sample_info.diag_L
            intermediate_sample_info = intermediate_sample_info._replace(alpha_hat=geom_mean_alpha)
            intermediate_sample_info = intermediate_sample_info._replace(rls_upper_bound=geom_mean_alpha * diag_L)
            intermediate_sample_info = intermediate_sample_info._replace(rls_upper_bound_valid=np.full((diag_L.shape[0],), False))

            ratio_alpha = intermediate_sample_info.alpha_max/intermediate_sample_info.alpha_min
            if ratio_alpha <= stopping_ratio and sample_count > 0:
                found_good_alpha = True
                break

            intermediate_sample_info = intermediate_sample_info._replace(alpha_switches=intermediate_sample_info.alpha_switches + 1)
            trial_count = 0
            under_k_count = 0
            over_k_count = 0
    else:
        raise ValueError('The alpha sampler reached the maximum number of rejections allowed '
                         'for the k-DPP size rejection ({}), try to increase the r factor '
                         '(see r_func parameter) or the Nystrom approximation accuracy '
                         'see rls_oversample_* parameters).'.format(max_iter_size_rejection))
    if found_good_alpha:
        intermediate_sample_info = intermediate_sample_info._replace(alpha_min=intermediate_sample_info.alpha_hat)
        intermediate_sample_info = intermediate_sample_info._replace(alpha_max=intermediate_sample_info.alpha_hat)
        intermediate_sample_info = intermediate_sample_info._replace(alpha_switches=intermediate_sample_info.alpha_switches + 1)

    if verbose_outer:
        params['verbose'] = verbose_outer
    else:
        params.pop('verbose')

    return sampl_out, intermediate_sample_info def k_dpp_eig_vecs_selector(eig_vals, eig_vecs, size,
                            E_poly=None, random_state=None):
    """ Subsample eigenvectors V of the 'L' kernel to build a projection DPP with kernel V V.T from which sampling is easy. The selection is made based a realization of Bernoulli variables with parameters the eigenvalues of 'L' and evalutations of the elementary symmetric polynomials.

    :param eig_vals:
        Collection of eigen values of 'L' (likelihood) kernel.
    :type eig_vals:
        list, array_like

    :param eig_vecs:
        Collection of eigenvectors of 'L' kernel.
    :type eig_vecs:
        array_like

    :param size:
        Size :math:`k` of :math:`k`-DPP
    :type size:
        int

    :param E_poly:
        Evaluation of symmetric polynomials in the eigenvalues
    :type E_poly:
        array_like

    :return:
        Selected eigenvectors
    :rtype:
        array_like

    .. seealso::

        - :cite:`KuTa12` Algorithm 8
        - :func:`elementary_symmetric_polynomials <elementary_symmetric_polynomials>`
    """

    rng = check_random_state(random_state)

    # Size of: ground set / sample
    N, k = eig_vecs.shape[0], size

    # as in np.linalg.matrix_rank
    tol = np.max(eig_vals) * N * np.finfo(np.float).eps
    rank = np.count_nonzero(eig_vals > tol)
    if k > rank:
        raise ValueError('size k={} > rank={}'.format(k, rank))

    if E_poly is None:
        E_poly = elementary_symmetric_polynomials(eig_vals, k)

    ind_selected = np.zeros(k, dtype=int)
    for n in range(eig_vals.size, 0, -1):

        if rng.rand() < eig_vals[n - 1] * E_poly[k - 1, n - 1] / E_poly[k, n]:
            k -= 1
            ind_selected[k] = n - 1
            if k == 0:
                break

    return eig_vecs[:, ind_selected] def elementary_symmetric_polynomials(eig_vals, size):
    """ Evaluate the elementary symmetric polynomials :math:`e_k` in the eigenvalues :math:`(\\lambda_1, \\cdots, \\lambda_N)`.

    :param eig_vals:
        Collection of eigenvalues :math:`(\\lambda_1, \\cdots, \\lambda_N)` of the similarity kernel :math:`L`.
    :type eig_vals:
        list

    :param size:
        Maximum degree of elementary symmetric polynomial.
    :type size:
        int

    :return:
        :math:`[E_{kn}]_{k=0, n=0}^{\text{size}, N}`
        :math:`E_{kn} = e_k(\\lambda_1, \\cdots, \\lambda_n)`
    :rtype:
        array_like

    .. seealso::

        - :cite:`KuTa12` Algorithm 7
        - `Wikipedia <https://en.wikipedia.org/wiki/Elementary_symmetric_polynomial>`_
    """

    # Initialize output array
    N = eig_vals.size
    E = np.zeros((size + 1, N + 1))
    E[0, :] = 1.0

    # Recursive evaluation
    for k in range(1, size + 1):
        for n in range(1, N + 1):
            E[k, n] = E[k, n - 1] + eig_vals[n - 1] * E[k - 1, n - 1]

    return E def dpp_sampler_mcmc(kernel, mode='AED', **params):
    """ Interface function with initializations and samplers for MCMC schemes.

    .. seealso::

        - :ref:`finite_dpps_mcmc_sampling_add_exchange_delete`
        - :func:`add_exchange_delete_sampler <add_exchange_delete_sampler>`
        - :func:`initialize_AED_sampler <initialize_AED_sampler>`
        - :func:`add_delete_sampler <add_delete_sampler>`
        - :func:`basis_exchange_sampler <basis_exchange_sampler>`
        - :func:`initialize_AD_and_E_sampler <initialize_AD_and_E_sampler>`
    """

    rng = check_random_state(params.get('random_state', None))

    s_init = params.get('s_init', None)
    nb_iter = params.get('nb_iter', 10)
    T_max = params.get('T_max', None)
    size = params.get('size', None)  # = Tr(K) for projection correlation K

    if mode == 'AED':  # Add-Exchange-Delete S'=S+t, S-t+u, S-t
        if s_init is None:
            s_init = initialize_AED_sampler(kernel, random_state=rng)
        sampl = add_exchange_delete_sampler(kernel, s_init, nb_iter, T_max,
                                            random_state=rng)

    elif mode == 'AD':  # Add-Delete S'=S+t, S-t
        if s_init is None:
            s_init = initialize_AD_and_E_sampler(kernel, random_state=rng)
        sampl = add_delete_sampler(kernel, s_init, nb_iter, T_max,
                                   random_state=rng)

    elif mode == 'E':  # Exchange S'=S-t+u
        if s_init is None:
            s_init = initialize_AD_and_E_sampler(kernel, size,
                                                 random_state=rng)
        sampl = basis_exchange_sampler(kernel, s_init, nb_iter, T_max,
                                       random_state=rng)

    return sampl def zonotope_sampler(A_zono, **params):
    """ MCMC based sampler for projection DPPs.
    The similarity matrix is the orthogonal projection matrix onto
    the row span of the feature vector matrix.
    Samples are of size equal to the ransampl_size of the projection matrix
    also equal to the rank of the feature matrix (assumed to be full row rank).

    :param A_zono:
        Feature vector matrix, feature vectors are stacked columnwise.
        It is assumed to be full row rank.
    :type A_zono:
        array_like

    :param params: Dictionary containing the parameters

        - ``'lin_obj'`` (list): Linear objective (:math:`c`) of the linear program used to identify the tile in which a point lies. Default is a random Gaussian vector.
        - ``'x_0'` (list): Initial point.
        - ``'nb_iter'`` (int): Number of iterations of the MCMC chain. Default is 10.
        - ``'T_max'`` (float): Maximum running time of the algorithm (in seconds).
        Default is None.
        - ``'random_state`` (default None)
    :type params: dict

    :return:
        MCMC chain of approximate samples (stacked row_wise i.e. nb_iter rows).
    :rtype:
        array_like

    .. seealso::

        Algorithm 5 in :cite:`GaBaVa17`

        - :func:`extract_basis <extract_basis>`
        - :func:`basis_exchange_sampler <basis_exchange_sampler>`
    """
    # For zonotope sampler
    try:
        from cvxopt import matrix, spmatrix, solvers
    except ImportError:
        raise ValueError('The cvxopt package is required to use the zonotype sampler (see setup.py).')

    solvers.options['show_progress'] = params.get('show_progress', False)
    solvers.options['glpk'] = {'msg_lev': params.get('show_progress', 'GLP_MSG_OFF')}

    rng = check_random_state(params.get('random_state', None))

    r, N = A_zono.shape  # Sizes of r=samples=rank(A_zono), N=ground set
    # Linear objective
    c = matrix(params.get('lin_obj', rng.randn(N)))
    # Initial point x0 = A*u, u~U[0,1]^n
    x0 = matrix(params.get('x_0', A_zono.dot(rng.rand(N))))

    nb_iter = params.get('nb_iter', 10)
    T_max = params.get('T_max', None)

    ###################
    # Linear problems #
    ###################
    # Canonical form
    # min       c.T*x         min     c.T*x
    # s.t.  G*x <= h    <=>   s.t.    G*x + s = h
    #        A*x = b                      A*x = b
    #                                      s >= 0
    # CVXOPT
    # =====> solvers.lp(c, G, h, A, b, solver='glpk')
    #################################################

    # To access the tile Z(B_x)
    # Solve P_x(A,c)
    ######################################################
    # y^* =
    # argmin  c.T*y               argmin  c.T*y
    # s.t.  A*y = x         <=>   s.t.  A  *y  = x
    #       0 <= y <= 1             [ I_n] *y <= [1^n]
    #                               [-I_n]       [0^n]
    ######################################################
    # Then B_x = \{ i ; y_i^* \in ]0,1[ \}

    A = spmatrix(0.0, [], [], (r, N))
    A[:, :] = A_zono

    G = spmatrix(0.0, [], [], (2 * N, N))
    G[:N, :] = spmatrix(1.0, range(N), range(N))
    G[N:, :] = spmatrix(-1.0, range(N), range(N))

    # Endpoints of segment
    # D_x \cap Z(A) = [x+alpha_m*d, x-alpha_M*d]
    ###########################################################################
    # alpha_m/_M = argmin  +/-alpha      argmin [+/-1 0^N].T * [alpha,lambda]
    # s.t.    x + alpha d = A lambda <=> s.t.  [-d A] *[alpha, lambda] = x
    #         0 <= lambda <= 1             [0^N I_N] *[alpha, lambda] <= [1^N]
    #                                      [0^N -I_N]                    [0^N]
    ##########################################################################

    c_mM = matrix(0.0, (N + 1, 1))
    c_mM[0] = 1.0

    A_mM = spmatrix(0.0, [], [], (r, N + 1))
    A_mM[:, 1:] = A

    G_mM = spmatrix(0.0, [], [], (2 * N, N + 1))
    G_mM[:, 1:] = G

    # Common h to both kind of LP
    # cf. 0 <= y <= 1 and 0 <= lambda <= 1
    h = matrix(0.0, (2 * N, 1))
    h[:N, :] = 1.0

    ##################
    # Initialization #
    ##################
    B_x0 = []
    while len(B_x0) != r:
        # Initial tile B_x0
        # Solve P_x0(A,c)
        y_star = solvers.lp(c, G, h, A, x0, solver='glpk')['x']
        # Get the tile
        B_x0 = extract_basis(np.asarray(y_star))

    # Initialize sequence of sample
    chain = np.zeros((nb_iter, r), dtype=int)
    chain[0] = B_x0

    # Compute the det of the tile (Vol(B)=abs(det(B)))
    det_B_x0 = la.det(A_zono[:, B_x0])

    t_start = time.time() if T_max else 0

    for it in range(1, nb_iter):

        # Take uniform direction d defining D_x0
        d = matrix(rng.randn(r, 1))

        # Define D_x0 \cap Z(A) = [x0 + alpha_m*d, x0 - alpha_M*d]
        # Update the constraint [-d A] * [alpha,lambda] = x
        A_mM[:, 0] = -d
        # Find alpha_m/M
        alpha_m = solvers.lp(c_mM, G_mM, h, A_mM, x0, solver='glpk')['x'][0]
        alpha_M = solvers.lp(-c_mM, G_mM, h, A_mM, x0, solver='glpk')['x'][0]

        # Propose x1 ~ U_{[x0+alpha_m*d, x0-alpha_M*d]}
        x1 = x0 + (alpha_m + (alpha_M - alpha_m) * rng.rand()) * d
        # Proposed tile B_x1
        # Solve P_x1(A,c)
        y_star = solvers.lp(c, G, h, A, x1, solver='glpk')['x']
        # Get the tile
        B_x1 = extract_basis(np.asarray(y_star))

        # Accept/Reject the move with proba Vol(B1)/Vol(B0)
        if len(B_x1) != r:  # if extract_basis returned smtg ill conditioned
            chain[it] = B_x0
        else:
            det_B_x1 = la.det(A_zono[:, B_x1])
            if rng.rand() < abs(det_B_x1 / det_B_x0):
                x0, B_x0, det_B_x0 = x1, B_x1, det_B_x1
                chain[it] = B_x1
            else:
                chain[it] = B_x0

        if T_max:
            if time.time() - t_start < T_max:
                break

    return chain.tolist() def check_random_state(seed: object) -> np.random.RandomState:
    """Turn seed into a np.random.RandomState instance
    Parameters
    ----------
    seed : None | int | instance of RandomState
        If seed is None, return the RandomState singleton used by np.random.
        If seed is an int, return a new RandomState instance seeded with seed.
        If seed is already a RandomState instance, return it.
        Otherwise raise ValueError.

    Returns
    -------
    np.random.RandomState

    .. seealso::

        `Scikit learn source code <https://github.com/scikit-learn/scikit-learn/blob/7813f7efb/sklearn/utils/validation.py#L763>`_
    """
    if seed is None or seed is np.random:
        return np.random.mtrand._rand
    if isinstance(seed, (int, np.integer)):
        return np.random.RandomState(seed)
    if isinstance(seed, np.random.RandomState):
        return seed
    raise ValueError('%r cannot be used to seed a numpy.random.RandomState'
                     ' instance' % seed) def is_symmetric(array):
    # Cheap test to check symmetry M^T = M

    if array is None:
        return None

    array = is_square(array)

    idx = np.arange(min(20, array.shape[0]))
    M = array[np.ix_(idx, idx)]
    if np.allclose(M.T, M):
        return array
    else:
        raise ValueError('array not symmetric: M.T != M') def is_projection(array, col_idx=None):
    # Cheap test to check reproducing property: M^2 = M

    if array is None:
        return None

    array = is_square(array)

    if col_idx is None:
        col_idx = np.arange(min(20, array.shape[0]))

    M_j = array[:, col_idx]
    Mjj = array[col_idx, col_idx]

    if np.allclose(inner1d(M_j), Mjj):
        return array
    else:
        raise ValueError('array not seem to be a projection: M^2 != M') def is_orthonormal_columns(array, col_idx=None):
    # Cheap test for checking orthonormality of columns of array: M.T M = I

    if array is None:
        return None

    if col_idx is None:
        col_idx = np.arange(np.min([5, array.shape[1]]))

    U = array[:, col_idx]

    if np.allclose(U.T.dot(U), np.eye(len(col_idx))):
        return array
    else:
        raise ValueError('array does not seem orthonormal: M.T M != I') def is_full_row_rank(array):
    # Check rank(M) = #rows

    if array is None:
        return None

    d, N = array.shape
    err_print = 'array (size = dxN) is not full row rank'

    if d > N:
        raise ValueError(err_print + 'd(={}) > N(={})'.format(d, N))
    else:
        rank = matrix_rank(array)
        if rank == d:
            return array
        else:
            raise ValueError(err_print + 'd(={}) != rank(={})'.format(d, rank)) def is_in_01(array, tol=1e-8):
    """Check if entries are **all** in :math:`[0, 1]`, for a given tolerance"""

    if array is None:
        return None
    elif np.all((-tol <= array) & (array <= 1.0 + tol)):
        return array
    else:
        raise ValueError('array with entries not all in [0,1]') def is_geq_0(array, tol=1e-8):
    """Check if entries are **all** :math:`\\geq0`, for a given tolerance"""

    if array is None:
        return None
    elif np.all(array >= -tol):
        return array
    else:
        raise ValueError('array with entries not all >= 0') def is_equal_to_O_or_1(array, tol=1e-8):
    """Check if entries are **all** in :math:`\\{0, 1\\}`, for a given tolerance"""

    if array is None:
        return None

    equal_0 = np.abs(array) <= tol
    equal_1 = np.abs(1 - array) <= tol
    equal_0_or_1 = equal_0 ^ equal_1  # ^ = xor

    if np.all(equal_0_or_1):
        return array
    else:
        raise ValueError('array with entries not all in {0,1}')
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
from warnings import warn
from dppy.exact_sampling import (dpp_sampler_generic_kernel,
                                 proj_dpp_sampler_kernel,
                                 proj_dpp_sampler_eig,
                                 dpp_vfx_sampler,
                                 alpha_dpp_sampler,
                                 dpp_eig_vecs_selector,
                                 k_dpp_vfx_sampler,
                                 alpha_k_dpp_sampler,
                                 k_dpp_eig_vecs_selector,
                                 elementary_symmetric_polynomials)
from dppy.mcmc_sampling import (dpp_sampler_mcmc,
                                zonotope_sampler)
from dppy.utils import (check_random_state,
                        is_symmetric,
                        is_projection,
                        is_orthonormal_columns,
                        is_full_row_rank,
                        is_in_01,
                        is_geq_0,
                        is_equal_to_O_or_1)
                                        E_poly=self.E_poly,
                                        random_state=rng)
            # Phase 2
            self.size_k_dpp = size
            sampl = proj_dpp_sampler_eig(V, self.sampling_mode,
                                         random_state=rng)

        elif self.K_eig_vals is not None:
            np.seterr(divide='raise')
            self.L_eig_vals = self.K_eig_vals / (1.0 - self.K_eig_vals)
            return self.sample_exact_k_dpp(size, self.sampling_mode,
                                           random_state=rng)

        # Otherwise eigendecomposition is necessary
        elif self.L_dual is not None:
            # L_dual = Phi Phi.T = W Theta W.T
            # L = Phi.T Phi = V Gamma V.T
            # implies Gamma = Theta and V = Phi.T W Theta^{-1/2}
            self.L_eig_vals, L_dual_eig_vecs = la.eigh(self.L_dual)
            self.L_eig_vals = is_geq_0(self.L_eig_vals)
            self.eig_vecs =self.L_gram_factor.T.dot(L_dual_eig_vecs
                                                    / np.sqrt(self.L_eig_vals))
            return self.sample_exact_k_dpp(size, mode=self.sampling_mode,
                                           random_state=rng)

        elif self.L is not None:
            self.L_eig_vals, self.eig_vecs = la.eigh(self.L)
            self.L_eig_vals = is_geq_0(self.L_eig_vals)
            return self.sample_exact_k_dpp(size, self.sampling_mode,
                                           random_state=rng)

        elif self.K is not None:
            self.K_eig_vals, self.eig_vecs = la.eigh(self.K)
            self.K_eig_vals = is_in_01(self.K_eig_vals)
            return self.sample_exact_k_dpp(size, self.sampling_mode,
                                           random_state=rng)

        elif self.eval_L is not None and self.X_data is not None:
            # In case mode!='vfx'
            self.compute_L()
            return self.sample_exact_k_dpp(size, self.sampling_mode,
                                           random_state=rng)

        else:
            raise ValueError('None of the available samplers could be used based on the current DPP representation. This should never happen, please consider rasing an issue on github at https://github.com/guilgautier/DPPy/issues')

        self.list_of_samples.append(sampl)
        return sampl

    # Approximate sampling
    def sample_mcmc(self, mode, **params):
        """ Run a MCMC with stationary distribution the corresponding :class:`FiniteDPP <FiniteDPP>` object.

        :param string mode:

            - ``'AED'`` Add-Exchange-Delete
            - ``'AD'`` Add-Delete
            - ``'E'`` Exchange
            - ``'zonotope'`` Zonotope sampling

        :param dict params:
            Dictionary containing the parameters for MCMC samplers with keys

            ``'random_state'`` (default None)

            - If ``mode='AED','AD','E'``

                + ``'s_init'`` (default None) Starting state of the Markov chain
                + ``'nb_iter'`` (default 10) Number of iterations of the chain
                + ``'T_max'`` (default None) Time horizon
                + ``'size'`` (default None) Size of the initial sample for ``mode='AD'/'E'``

                    * :math:`\\operatorname{rank}(\\mathbf{K})=\\operatorname{trace}(\\mathbf{K})` for projection :math:`\\mathbf{K}` (correlation) kernel and ``mode='E'``

            - If ``mode='zonotope'``:

                + ``'lin_obj'`` linear objective in main optimization problem (default np.random.randn(N))
                + ``'x_0'`` initial point in zonotope (default A*u, u~U[0,1]^n)
                + ``'nb_iter'`` (default 10) Number of iterations of the chain
                + ``'T_max'`` (default None) Time horizon

        :return:
            The last sample of the trajectory of Markov chain.

            In any case, the full trajectory of the Markov chain, made of ``params['nb_iter']`` samples, is appended to the :py:attr:`~FiniteDPP.list_of_samples` attribute as a list of lists.

        :rtype:
            list

        .. note::

            Each time you call this method, the full trajectory of the Markov chain, made of ``params['nb_iter']`` samples, is appended to the :py:attr:`~FiniteDPP.list_of_samples` attribute as a list of lists.

            The :py:attr:`~FiniteDPP.list_of_samples` attribute can be emptied using :py:meth:`~FiniteDPP.flush_samples`

        .. seealso::

            - :ref:`finite_dpps_mcmc_sampling`
            - :py:meth:`~FiniteDPP.sample_exact`
            - :py:meth:`~FiniteDPP.flush_samples`
        """

        self.sampling_mode = mode

        if self.sampling_mode == 'zonotope':
            if self.A_zono is not None:
                chain = zonotope_sampler(self.A_zono, **params)
            else:
                err_print = ['Invalid `mode=zonotope` parameter',
                             'DPP must be defined via `A_zono`',
                             'Given: {}'.format(self.params_keys)]
                raise ValueError(' '.join(err_print))

        elif self.sampling_mode == 'E':
            if self.kernel_type == 'correlation' and self.projection:
                self.compute_K()
                size = params.get('size', None)
                rank = np.rint(np.trace(self.K)).astype(int)
                # |sample| = Tr(K) = rank(K) a.s. for projection DPP(K)
                if size == rank:

                    chain = dpp_sampler_mcmc(self.K,