# Path: dppy/utils.py
# def det_ST(array, S, T=None):
#     """ Compute :math:`\\det M_{S, T} = \\det [M_{ij}]_{i\\inS, j\\in T}`
# 
#     :param array:
#         Matrix
#     :type M:
#         array_like
# 
#     :param S:
#         collection of indices
#     :type M:
#         1D list, array_like
# 
#     :param T:
#         collection of indices
#     :type M:
#         1D list, array_like, default None
# 
#     :return:
#         - if `T is None` return :math:`\\det M_{S, S}`
#         - else return :math:`\\det M_{S, T}`. If S=T=[], numpy convention = 1.0
#     :rtype:
#         float
#     """
# 
#     if T is None:  # det M_SS = det M_S
#         return det(array[np.ix_(S, S)])
# 
#     else:  # det M_ST, numpy deals with det M_[][] = 1.0
#         return det(array[np.ix_(S, T)])
# 
# 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)


import time
import numpy as np
import scipy.linalg as la
from dppy.utils import det_ST, check_random_state
        from cvxopt import matrix, spmatrix, solvers

        # Add: S1 = S0 + t
        if U < 0.5 * (1 - ratio)**2:
            S1.append(t)  # S1 = S0 + t
            # Accept_reject the move
            det_S1 = det_ST(kernel, S1)  # det K_S1
            if rng.rand() < det_S1 / det_S0 * (size_S0 + 1) / (N - size_S0):
                S0, det_S0 = S1, det_S1
                chain.append(S1)
                size_S0 += 1
            else:
                chain.append(S0)

        # Exchange: S1 = S0 - s + t
        elif (0.5 * (1 - ratio)**2 <= U) & (U < 0.5 * (1 - ratio)):
            del S1[s_ind]  # S1 = S0 - s
            S1.append(t)  # S1 = S1 + t = S0 - s + t
            # Accept_reject the move
            det_S1 = det_ST(kernel, S1)  # det K_S1
            if rng.rand() < (det_S1 / det_S0):
                S0, det_S0 = S1, det_S1
                chain.append(S1)
                # size_S0 stays the same
            else:
                chain.append(S0)

        # Delete: S1 = S0 - s
        elif (0.5 * (1 - ratio) <= U) & (U < 0.5 * (ratio**2 + (1 - ratio))):
            del S1[s_ind]  # S0 - s
            # Accept_reject the move

# Path: dppy/utils.py
# def det_ST(array, S, T=None):
#     """ Compute :math:`\\det M_{S, T} = \\det [M_{ij}]_{i\\inS, j\\in T}`
# 
#     :param array:
#         Matrix
#     :type M:
#         array_like
# 
#     :param S:
#         collection of indices
#     :type M:
#         1D list, array_like
# 
#     :param T:
#         collection of indices
#     :type M:
#         1D list, array_like, default None
# 
#     :return:
#         - if `T is None` return :math:`\\det M_{S, S}`
#         - else return :math:`\\det M_{S, T}`. If S=T=[], numpy convention = 1.0
#     :rtype:
#         float
#     """
# 
#     if T is None:  # det M_SS = det M_S
#         return det(array[np.ix_(S, S)])
# 
#     else:  # det M_ST, numpy deals with det M_[][] = 1.0
#         return det(array[np.ix_(S, T)])
# 
# 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)


# Path: dppy/mcmc_sampling.py
import time
import numpy as np
import scipy.linalg as la
from dppy.utils import det_ST, check_random_state
        from cvxopt import matrix, spmatrix, solvers


        # Add: S1 = S0 + t
        if U < 0.5 * (1 - ratio)**2:
            S1.append(t)  # S1 = S0 + t
            # Accept_reject the move
            det_S1 = det_ST(kernel, S1)  # det K_S1
            if rng.rand() < det_S1 / det_S0 * (size_S0 + 1) / (N - size_S0):
                S0, det_S0 = S1, det_S1
                chain.append(S1)
                size_S0 += 1
            else:
                chain.append(S0)

        # Exchange: S1 = S0 - s + t
        elif (0.5 * (1 - ratio)**2 <= U) & (U < 0.5 * (1 - ratio)):
            del S1[s_ind]  # S1 = S0 - s
            S1.append(t)  # S1 = S1 + t = S0 - s + t
            # Accept_reject the move
            det_S1 = det_ST(kernel, S1)  # det K_S1
            if rng.rand() < (det_S1 / det_S0):
                S0, det_S0 = S1, det_S1
                chain.append(S1)
                # size_S0 stays the same
            else:
                chain.append(S0)

        # Delete: S1 = S0 - s
        elif (0.5 * (1 - ratio) <= U) & (U < 0.5 * (ratio**2 + (1 - ratio))):
            del S1[s_ind]  # S0 - s
            # Accept_reject the move

            det_S1 = det_ST(kernel, S1)  # det K_S1