# Source code for qutip.entropy

    __all__ = ['entropy_vn', 'entropy_linear', 'entropy_mutual', 'negativity',
               'concurrence', 'entropy_conditional', 'entangling_power',
               'entropy_relative']

    from .core.numpy_backend import np
    from .partial_transpose import partial_transpose
    from . import (ptrace, tensor, sigmay, ket2dm,
                   expand_operator)
    from .core import data as _data
    from .core.dimensions import Dimensions

    def entropy_vn(rho, base=np.e, sparse=False):
        """
        Von-Neumann entropy of density matrix

        Parameters
        ----------
        rho : qobj
            Density matrix.
        base : {e, 2}, default: e
            Base of logarithm.
        sparse : bool, default: False
            Use sparse eigensolver.

        Returns
        -------
        entropy : float
            Von-Neumann entropy of `rho`.

        Examples
        --------
        >>> rho=0.5*fock_dm(2,0)+0.5*fock_dm(2,1)
        >>> entropy_vn(rho,2)
        1.0

        """
        if rho.type == 'ket' or rho.type == 'bra':
            rho = ket2dm(rho)
        vals = rho.eigenenergies(sparse=sparse)
        threshold = 1e-17
        nzvals = np.where(vals < threshold, threshold, vals)
        if base == 2:
            logvals = np.log2(nzvals)
        elif base == np.e:
            logvals = np.log(nzvals)
        else:
            raise ValueError("Base must be 2 or e.")
        return np.real(-sum(nzvals * logvals))

    def entropy_linear(rho):
        """
        Linear entropy of a density matrix.

        Parameters
        ----------
        rho : qobj
            sensity matrix or ket/bra vector.

        Returns
        -------
        entropy : float
            Linear entropy of rho.

        Examples
        --------
        >>> rho=0.5*fock_dm(2,0)+0.5*fock_dm(2,1)
        >>> entropy_linear(rho)
        0.5

        """
        if rho.type == 'ket' or rho.type == 'bra':
            rho = ket2dm(rho)
        return np.real(1.0 - (rho ** 2).tr())

    def concurrence(rho):
        """
        Calculate the concurrence entanglement measure for a two-qubit state.

        Parameters
        ----------
        state : qobj
            Ket, bra, or density matrix for a two-qubit state.

        Returns
        -------
        concur : float
            Concurrence

        References
        ----------

        .. [1] `https://en.wikipedia.org/wiki/Concurrence_(quantum_computing)`

        """
        if rho.isket and rho._dims != Dimensions([[2, 2], [1, 1]]):
            raise Exception("Ket must be tensor product of two qubits.")

        elif rho.isbra and rho._dims != Dimensions([[1, 1], [2, 2]]):
            raise Exception("Bra must be tensor product of two qubits.")

        elif rho.isoper and rho._dims != Dimensions([[2, 2], [2, 2]]):
            raise Exception("Density matrix must be tensor product of two qubits.")

        if rho.isket or rho.isbra:
            rho = ket2dm(rho)

        sysy = tensor(sigmay(), sigmay())

        rho_tilde = (rho * sysy) * (rho.conj() * sysy)

        evals = rho_tilde.eigenenergies()

        # abs to avoid problems with sqrt for very small negative numbers
        evals = abs(np.sort(np.real(evals)))

        sqrt_evals = np.sqrt(evals)
        lsum = sqrt_evals[3] - sqrt_evals[2] - sqrt_evals[1] - sqrt_evals[0]
        return np.maximum(0, lsum)

    def negativity(rho, subsys, method='tracenorm', logarithmic=False):
        """
        Compute the negativity for a multipartite quantum system described
        by the density matrix rho. The subsys argument is an index that
        indicates which system to compute the negativity for.

        .. note::

            Experimental.
        """
        rho._dims._require_pure_dims("negativity")
        if rho.isket or rho.isbra:
            rho = ket2dm(rho)
        mask = [idx == subsys for idx, n in enumerate(rho.dims[0])]
        rho_pt = partial_transpose(rho, mask)

        if method == 'tracenorm':
            N = ((rho_pt.dag() * rho_pt).sqrtm().tr().real - 1)/2.0
        elif method == 'eigenvalues':
            l = rho_pt.eigenenergies()
            N = ((abs(l)-l)/2).sum()
        else:
            raise ValueError(f"Unknown method {method}")

    # Return the negativity value (or its logarithm if specified)
        if logarithmic:
            return np.log2(2 * N + 1)
        else:
            return N

    def entropy_mutual(rho, selA, selB, base=np.e, sparse=False):
        """
        Calculates the mutual information S(A:B) between selection
        components of a system density matrix.

        Parameters
        ----------
        rho : qobj
            Density matrix for composite quantum systems
        selA : int/list
            `int` or `list` of first selected density matrix components.
        selB : int/list
            `int` or `list` of second selected density matrix components.
        base : {e, 2}, default: e
            Base of logarithm.
        sparse : bool, default: False
            Use sparse eigensolver.

        Returns
        -------
        ent_mut : float
           Mutual information between selected components.

        """
        rho._dims._require_pure_dims("mutual information")
        if isinstance(selA, int):
            selA = [selA]
        if isinstance(selB, int):
            selB = [selB]
        if rho.type != 'oper':
            raise TypeError("Input must be a density matrix.")
        if (len(selA) + len(selB)) != len(rho.dims[0]):
            raise TypeError("Number of selected components must match " +
                            "total number.")

        rhoA = ptrace(rho, selA)
        rhoB = ptrace(rho, selB)
        out = (entropy_vn(rhoA, base, sparse=sparse) +
               entropy_vn(rhoB, base, sparse=sparse) -
               entropy_vn(rho, base, sparse=sparse))
        return out

    def entropy_relative(rho, sigma, base=np.e, sparse=False, tol=1e-12):
        """
        Calculates the relative entropy S(rho||sigma) between two density
        matrices.

        Parameters
        ----------
        rho : :class:`.Qobj`
            First density matrix (or ket which will be converted to a density
            matrix).
        sigma : :class:`.Qobj`
            Second density matrix (or ket which will be converted to a density
            matrix).
        base : {e, 2}, default: e
            Base of logarithm. Defaults to e.
        sparse : bool, default: False
            Flag to use sparse solver when determining the eigenvectors
            of the density matrices. Defaults to False.
        tol : float, default: 1e-12
            Tolerance to use to detect 0 eigenvalues or dot producted between
            eigenvectors. Defaults to 1e-12.

        Returns
        -------
        rel_ent : float
            Value of relative entropy. Guaranteed to be greater than zero
            and should equal zero only when rho and sigma are identical.

        Examples
        --------

        First we define two density matrices:

        >>> rho = qutip.ket2dm(qutip.ket("00"))
        >>> sigma = rho + qutip.ket2dm(qutip.ket("01"))
        >>> sigma = sigma.unit()

        Then we calculate their relative entropy using base 2 (i.e. ``log2``)
        and base e (i.e. ``log``).

        >>> qutip.entropy_relative(rho, sigma, base=2)
        1.0
        >>> qutip.entropy_relative(rho, sigma)
        0.6931471805599453

        References
        ----------

        See Nielsen & Chuang, "Quantum Computation and Quantum Information",
        Section 11.3.1, pg. 511 for a detailed explanation of quantum relative
        entropy.
        """
        if rho.isket:
            rho = ket2dm(rho)
        if sigma.isket:
            sigma = ket2dm(sigma)
        if not rho.isoper or not sigma.isoper:
            raise TypeError("Inputs must be density matrices.")
        if rho._dims != sigma._dims:
            raise ValueError("Inputs must have the same shape and dims.")
        if base == 2:
            log_base = np.log2
        elif base == np.e:
            log_base = np.log
        else:
            raise ValueError("Base must be 2 or e.")
        # S(rho || sigma) = sum_i(p_i log p_i) - sum_ij(p_i P_ij log q_i)
        #
        # S is +inf if the kernel of sigma (i.e. svecs[svals == 0]) has non-trivial
        # intersection with the support of rho (i.e. rvecs[rvals != 0]).
        rvals, rvecs = _data.eigs(rho.data, rho.isherm, True)
        rvecs = rvecs.to_array().T
        if any(abs(np.imag(rvals)) >= tol):
            raise ValueError("Input rho has non-real eigenvalues.")
        rvals = np.real(rvals)
        svals, svecs = _data.eigs(sigma.data, sigma.isherm, True)
        svecs = svecs.to_array().T
        if any(abs(np.imag(svals)) >= tol):
            raise ValueError("Input sigma has non-real eigenvalues.")
        svals = np.real(svals)
        # Calculate inner products of eigenvectors and return +inf if kernel
        # of sigma overlaps with support of rho.
        P = abs(np.inner(rvecs, np.conj(svecs))) ** 2
        if (rvals >= tol) @ (P >= tol) @ (svals < tol):
            return np.inf
        # Avoid -inf from log(0) -- these terms will be multiplied by zero later
        # anyway
        svals[abs(svals) < tol] = 1
        nzrvals = rvals[abs(rvals) >= tol]
        # Calculate S
        S = nzrvals @ log_base(nzrvals) - rvals @ P @ log_base(svals)
        # the relative entropy is guaranteed to be >= 0, so we clamp the
        # calculated value to 0 to avoid small violations of the lower bound.
        return np.maximum(0, S)

    def entropy_conditional(rho, selB, base=np.e, sparse=False):
        """
        Calculates the conditional entropy :math:`S(A|B)=S(A,B)-S(B)`
        of a selected density matrix component.

        Parameters
        ----------
        rho : qobj
            Density matrix of composite object
        selB : int/list
            Selected components for density matrix B
        base : {e, 2}, default: e
            Base of logarithm.
        sparse : bool, default: False
            Use sparse eigensolver.

        Returns
        -------
        ent_cond : float
            Value of conditional entropy

        """
        rho._dims._require_pure_dims("conditional entropy")

        if rho.type != 'oper':
            raise TypeError("Input must be density matrix.")
        if isinstance(selB, int):
            selB = [selB]
        B = ptrace(rho, selB)
        out = (entropy_vn(rho, base, sparse=sparse) -
               entropy_vn(B, base, sparse=sparse))
        return out

    def participation_ratio(rho):
        """
        Returns the effective number of states for a density matrix.

        The participation is unity for pure states, and maximally N,
        where N is the Hilbert space dimensionality, for completely
        mixed states.

        Parameters
        ----------
        rho : qobj
            Density matrix

        Returns
        -------
        pr : float
            Effective number of states in the density matrix

        """
        if rho.type == 'ket' or rho.type == 'bra':
            return 1.0
        else:
            return 1.0 / (rho ** 2).tr()

    def entangling_power(U):
        """
        Calculate the entangling power of a two-qubit gate U, which
        is zero of nonentangling gates and 2/9 for maximally entangling gates.

        Parameters
        ----------
        U : qobj
            Qobj instance representing a two-qubit gate.

        Returns
        -------
        ep : float
            The entanglement power of U (real number between 0 and 2/9)

        References:

            Explorations in Quantum Computing, Colin P. Williams (Springer, 2011)
        """

        if not U.isoper:
            raise Exception("U must be an operator.")

        if U._dims != Dimensions([[2, 2], [2, 2]]):
            raise Exception("U must be a two-qubit gate.")

        from qutip.core.gates import swap
        swap13 =  expand_operator(swap(dtype=U.dtype), [2, 2, 2, 2], [1, 3])
        a = tensor(U, U).dag() * swap13 * tensor(U, U) * swap13
        Uswap = swap(dtype=U.dtype) * U
        b = tensor(Uswap, Uswap).dag() * swap13 * tensor(Uswap, Uswap) * swap13

        return 5.0/9 - 1.0/36 * (a.tr() + b.tr()).real