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
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