Attempt to automatically do adiabatic elimination on an SLH object

    This will project the `Y` operator onto a truncated basis with dimension
    specified by `fock_trunc`.  `sub_P0` controls whether an attempt is made to
    replace the kernel projector P0 by an :class:`.IdentityOperator`.

def try_adiabatic_elimination(slh, k=None, fock_trunc=6, sub_P0=True):
    """Attempt to automatically do adiabatic elimination on an SLH object

    This will project the `Y` operator onto a truncated basis with dimension
    specified by `fock_trunc`.  `sub_P0` controls whether an attempt is made to
    replace the kernel projector P0 by an :class:`.IdentityOperator`.
    """
    ops = prepare_adiabatic_limit(slh, k)
    Y = ops[0]
    if isinstance(Y.space, LocalSpace):
        try:
            b = Y.space.basis_labels
            if len(b) > fock_trunc:
                b = b[:fock_trunc]
        except BasisNotSetError:
            b = range(fock_trunc)
        projectors = set(LocalProjector(ll, hs=Y.space) for ll in b)
        Id_trunc = sum(projectors, ZeroOperator)
        Yprojection = (
            ((Id_trunc * Y).expand() * Id_trunc)
            .expand().simplify_scalar())
        termcoeffs = get_coeffs(Yprojection)
        terms = set(termcoeffs.keys())

        for term in terms - projectors:
            cannot_eliminate = (
                not isinstance(term, LocalSigma) or
                not term.operands[1] == term.operands[2])
            if cannot_eliminate:
                raise CannotEliminateAutomatically(
                    "Proj. Y operator has off-diagonal term: ~{}".format(term))
        P0 = sum(projectors - terms, ZeroOperator)
        if P0 == ZeroOperator:
            raise CannotEliminateAutomatically("Empty null-space of Y!")

        Yinv = sum(t/termcoeffs[t] for t in terms & projectors)
        assert (
            (Yprojection*Yinv).expand().simplify_scalar() ==
            (Id_trunc - P0).expand())
        slhlim = eval_adiabatic_limit(ops, Yinv, P0)

        if sub_P0:
            # TODO for non-unit rank P0, this will not work
            slhlim = slhlim.substitute(
                {P0: IdentityOperator}).expand().simplify_scalar()
        return slhlim

    else:
        raise CannotEliminateAutomatically(
            "Currently only single degree of freedom Y-operators supported")