Return the fidelity F(rho0, rho1) between two mixed quantum states.

    Note: Fidelity cannot be calculated entirely within the tensor backend.

def fidelity(rho0: Density, rho1: Density) -> float:
    """Return the fidelity F(rho0, rho1) between two mixed quantum states.

    Note: Fidelity cannot be calculated entirely within the tensor backend.
    """
    assert rho0.qubit_nb == rho1.qubit_nb   # FIXME

    rho1 = rho1.permute(rho0.qubits)
    op0 = asarray(rho0.asoperator())
    op1 = asarray(rho1.asoperator())

    fid = np.real((np.trace(sqrtm(sqrtm(op0) @ op1 @ sqrtm(op0)))) ** 2)
    fid = min(fid, 1.0)
    fid = max(fid, 0.0)     # Correct for rounding errors

    return fid