<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def symbolic_heisenberg_eom( self, X=None, noises=None, expand_simplify=True):
 """Compute the symbolic Heisenberg equations of motion of a system operator X. If no X is given, an OperatorSymbol is created in its place. If no noises are given, this correspnds to the ensemble-averaged Heisenberg equation of motion. Args: X (Operator):
 A system operator noises (Operator):
 A vector of noise inputs Returns: Operator: The RHS of the Heisenberg equations of motion of X. """

        L, H = self.L, self.H

        if X is None:
            X = OperatorSymbol('X', hs=(L.space | H.space))

        summands = [I * (H * X - X * H), ]
        for Lk in L.matrix.ravel():
            summands.append(adjoint(Lk) * X * Lk)
            summands.append(-(adjoint(Lk) * Lk * X + X * adjoint(Lk) * Lk) / 2)

        if noises is not None:
            if not isinstance(noises, Matrix):
                noises = Matrix(noises)
            LambdaT = (noises.adjoint().transpose() * noises.transpose()).transpose()
            assert noises.shape == L.shape
            S = self.S
            summands.append((adjoint(noises) * S.adjoint() * (X * L - L * X))
                            .expand()[0, 0])
            summand = (((L.adjoint() * X - X * L.adjoint()) * S * noises)
                       .expand()[0, 0])
            summands.append(summand)
            if len(S.space & X.space):
                comm = (S.adjoint() * X * S - X)
                summands.append((comm * LambdaT).expand().trace())

        ret = OperatorPlus.create(*summands)
        if expand_simplify:
            ret = ret.expand().simplify_scalar()
        return ret