def observable(operator, rho, unfolding, complex=False):
    r"""Return an observable ammount.

    INPUT:

    -  ``operator`` - An square matrix representing a hermitian operator \
    in thesame basis as the density matrix.
    -  ``rho`` - A density matrix in unfolded format, or a list of such \
    density matrices.
    -  ``unfolding`` - A mapping from matrix element indices to unfolded \
    indices.

    >>> Ne = 2
    >>> unfolding = Unfolding(Ne, True, True, True)
    >>> rho = np.array([[0.6, 1+2j], [1-2j, 0.4]])
    >>> rho = unfolding(rho)
    >>> sx = np.array([[0, 1], [1, 0]])
    >>> print(observable(sx, rho, unfolding))
    2.0

    """
    if len(rho.shape) == 2:
        return np.array([observable(operator, i, unfolding) for i in rho])

    Ne = unfolding.Ne
    Mu = unfolding.Mu
    obs = 0

    if unfolding.normalized:
        rho11 = 1 - sum([rho[Mu(1, i, i)] for i in range(1, Ne)])

    for i in range(Ne):
        for k in range(Ne):
            if unfolding.real:
                if k == 0 and i == 0:
                    obs += operator[i, k]*rho11
                else:
                    if k < i:
                        u, v = (i, k)
                    else:
                        u, v = (k, i)
                    obs += operator[i, k]*rho[Mu(1, u, v)]
                    if k != i:
                        if k < i:
                            obs += 1j*operator[i, k]*rho[Mu(-1, u, v)]
                        else:
                            obs += -1j*operator[i, k]*rho[Mu(-1, u, v)]

            else:
                if k == 0 and i == 0:
                    obs += operator[i, k]*rho11
                else:
                    obs += operator[i, k]*rho[Mu(0, k, i)]

    if not complex:
        obs = np.real(obs)
    return obs
def observable ( operator , rho , unfolding , complex = False ) : if len ( rho . shape ) == 2 : return np . array ( [ observable ( operator , i , unfolding ) for i in rho ] ) Ne = unfolding . Ne Mu = unfolding . Mu obs = 0 if unfolding . normalized : rho11 = 1 - sum ( [ rho [ Mu ( 1 , i , i ) ] for i in range ( 1 , Ne ) ] ) for i in range ( Ne ) : for k in range ( Ne ) : if unfolding . real : if k == 0 and i == 0 : obs += operator [ i , k ] * rho11 else : if k < i : u , v = ( i , k ) else : u , v = ( k , i ) obs += operator [ i , k ] * rho [ Mu ( 1 , u , v ) ] if k != i : if k < i : obs += 1j * operator [ i , k ] * rho [ Mu ( - 1 , u , v ) ] else : obs += - 1j * operator [ i , k ] * rho [ Mu ( - 1 , u , v ) ] else : if k == 0 and i == 0 : obs += operator [ i , k ] * rho11 else : obs += operator [ i , k ] * rho [ Mu ( 0 , k , i ) ] if not complex : obs = np . real ( obs ) return obs
r"""Return an observable ammount.

    INPUT:

    -  ``operator`` - An square matrix representing a hermitian operator \
    in thesame basis as the density matrix.
    -  ``rho`` - A density matrix in unfolded format, or a list of such \
    density matrices.
    -  ``unfolding`` - A mapping from matrix element indices to unfolded \
    indices.

    >>> Ne = 2
    >>> unfolding = Unfolding(Ne, True, True, True)
    >>> rho = np.array([[0.6, 1+2j], [1-2j, 0.4]])
    >>> rho = unfolding(rho)
    >>> sx = np.array([[0, 1], [1, 0]])
    >>> print(observable(sx, rho, unfolding))
    2.0