def convert_to_sympy_matrix(expr, full_space=None):
    """Convert a QNET expression to an explicit ``n x n`` instance of
    `sympy.Matrix`, where ``n`` is the dimension of `full_space`. The entries
    of the matrix may contain symbols.

    Parameters:
        expr: a QNET expression
        full_space (qnet.algebra.hilbert_space_algebra.HilbertSpace): The
            Hilbert space in which `expr` is defined. If not given,
            ``expr.space`` is used. The Hilbert space must have a well-defined
            basis.

    Raises:
        qnet.algebra.hilbert_space_algebra.BasisNotSetError: if `full_space`
            does not have a defined basis
        ValueError: if `expr` is not in `full_space`, or if `expr` cannot be
            converted.
    """
    if full_space is None:
        full_space = expr.space
    if not expr.space.is_tensor_factor_of(full_space):
        raise ValueError("expr must be in full_space")
    if expr is IdentityOperator:
        return sympy.eye(full_space.dimension)
    elif expr is ZeroOperator:
        return 0
    elif isinstance(expr, LocalOperator):
        n = full_space.dimension
        if full_space != expr.space:
            all_spaces = full_space.local_factors
            own_space_index = all_spaces.index(expr.space)
            factors = [sympy.eye(s.dimension)
                       for s in all_spaces[:own_space_index]]
            factors.append(convert_to_sympy_matrix(expr, expr.space))
            factors.extend([sympy.eye(s.dimension)
                            for s in all_spaces[own_space_index + 1:]])
            return tensor(*factors)
        if isinstance(expr, (Create, Jz, Jplus)):
            return SympyCreate(n)
        elif isinstance(expr, (Destroy, Jminus)):
            return SympyCreate(n).H
        elif isinstance(expr, Phase):
            phi = expr.phase
            result = sympy.zeros(n)
            for i in range(n):
                result[i, i] = sympy.exp(sympy.I * i * phi)
            return result
        elif isinstance(expr, Displace):
            alpha = expr.operands[1]
            a = SympyCreate(n)
            return (alpha * a - alpha.conjugate() * a.H).exp()
        elif isinstance(expr, Squeeze):
            eta = expr.operands[1]
            a = SympyCreate(n)
            return ((eta/2) * a**2 - (eta.conjugate()/2) * (a.H)**2).exp()
        elif isinstance(expr, LocalSigma):
            ket = basis_state(expr.index_j, n)
            bra = basis_state(expr.index_k, n).H
            return ket * bra
        else:
            raise ValueError("Cannot convert '%s' of type %s"
                             % (str(expr), type(expr)))
    elif (isinstance(expr, Operator) and isinstance(expr, Operation)):
        if isinstance(expr, OperatorPlus):
            s = convert_to_sympy_matrix(expr.operands[0], full_space)
            for op in expr.operands[1:]:
                s += convert_to_sympy_matrix(op, full_space)
            return s
        elif isinstance(expr, OperatorTimes):
            # if any factor acts non-locally, we need to expand distributively.
            if any(len(op.space) > 1 for op in expr.operands):
                se = expr.expand()
                if se == expr:
                    raise ValueError("Cannot represent as sympy matrix: %s"
                                     % expr)
                return convert_to_sympy_matrix(se, full_space)
            all_spaces = full_space.local_factors
            by_space = []
            ck = 0
            for ls in all_spaces:
                # group factors by associated local space
                ls_ops = [convert_to_sympy_matrix(o, o.space)
                          for o in expr.operands if o.space == ls]
                if len(ls_ops):
                    # compute factor associated with local space
                    by_space.append(ls_ops[0])
                    for ls_op in ls_ops[1:]:
                        by_space[-1] *= ls_op
                    ck += len(ls_ops)
                else:
                    # if trivial action, take identity matrix
                    by_space.append(sympy.eye(ls.dimension))
            assert ck == len(expr.operands)
            # combine local factors in tensor product
            if len(by_space) == 1:
                return by_space[0]
            else:
                return tensor(*by_space)
        elif isinstance(expr, Adjoint):
            return convert_to_sympy_matrix(expr.operand, full_space).H
        elif isinstance(expr, PseudoInverse):
            raise NotImplementedError(
                'Cannot convert PseudoInverse to sympy matrix')
        elif isinstance(expr, NullSpaceProjector):
            raise NotImplementedError(
                'Cannot convert NullSpaceProjector to sympy')
        elif isinstance(expr, ScalarTimesOperator):
            return expr.coeff * convert_to_sympy_matrix(expr.term, full_space)
        else:
            raise ValueError(
                "Cannot convert '%s' of type %s" % (str(expr), type(expr)))
    else:
        raise ValueError(
            "Cannot convert '%s' of type %s" % (str(expr), type(expr)))
Convert a QNET expression to an explicit ``n x n`` instance of
    `sympy.Matrix`, where ``n`` is the dimension of `full_space`. The entries
    of the matrix may contain symbols.

    Parameters:
        expr: a QNET expression
        full_space (qnet.algebra.hilbert_space_algebra.HilbertSpace): The
            Hilbert space in which `expr` is defined. If not given,
            ``expr.space`` is used. The Hilbert space must have a well-defined
            basis.

    Raises:
        qnet.algebra.hilbert_space_algebra.BasisNotSetError: if `full_space`
            does not have a defined basis
        ValueError: if `expr` is not in `full_space`, or if `expr` cannot be
            converted.