[Bugfix] Hamiltonian is rendered in top level docs (#5640)

Due to the dispatch of `qml.Hamiltonian` that involves removing
`Hamiltonian` from the top level in `pennylane/__init__.py`, the doc was
not displayed anymore. This fix "hacks" sphinx to add the missing obect
back, at least for sphinx.

Optional separate issue: A user seeking out the `qml.Hamiltonian` docs
at the top level may be confused because this now dispatches to
`LinearCombination` but this is not mentioned there. Added a warning to
avoid that confusion.

---------

Co-authored-by: Thomas R. Bromley <49409390+trbromley@users.noreply.github.com>
Co-authored-by: Astral Cai <astral.cai@xanadu.ai>

doc/conf.py

@@ -108,6 +108,7 @@
 # built documents.
 
 import pennylane
+pennylane.Hamiltonian = pennylane.ops.Hamiltonian
 
 # The full version, including alpha/beta/rc tags.
 release = pennylane.__version__

pennylane/ops/qubit/hamiltonian.py

@@ -67,6 +67,10 @@ class Hamiltonian(Observable):
     The Hamiltonian is represented as a linear combination of other operators, e.g.,
     :math:`\sum_{k=0}^{N-1} c_k O_k`, where the :math:`c_k` are trainable parameters.
 
+    .. warning::
+
+        As of ``v0.36``, ``qml.Hamiltonian`` dispatches to :class:`~.pennylane.ops.op_math.LinearCombination`
+        by default. For further details, see :doc:`Updated Operators </news/new_opmath/>`.
 
     Args:
         coeffs (tensor_like): coefficients of the Hamiltonian expression
@@ -83,83 +87,137 @@ class Hamiltonian(Observable):
 
     **Example:**
 
-    A Hamiltonian can be created by simply passing the list of coefficients
-    as well as the list of observables:
+    .. note::
+        As of ``v0.36``, ``qml.Hamiltonian`` dispatches to :class:`~.pennylane.ops.op_math.LinearCombination`
+        by default, so the following examples assume this behaviour.
+
+    ``qml.Hamiltonian`` takes in a list of coefficients and a list of operators.
 
     >>> coeffs = [0.2, -0.543]
     >>> obs = [qml.X(0) @ qml.Z(1), qml.Z(0) @ qml.Hadamard(2)]
     >>> H = qml.Hamiltonian(coeffs, obs)
     >>> print(H)
-      (-0.543) [Z0 H2]
-    + (0.2) [X0 Z1]
+    0.2 * (X(0) @ Z(1)) + -0.543 * (Z(0) @ Hadamard(wires=[2]))
 
     The coefficients can be a trainable tensor, for example:
 
     >>> coeffs = tf.Variable([0.2, -0.543], dtype=tf.double)
     >>> obs = [qml.X(0) @ qml.Z(1), qml.Z(0) @ qml.Hadamard(2)]
     >>> H = qml.Hamiltonian(coeffs, obs)
     >>> print(H)
-      (-0.543) [Z0 H2]
-    + (0.2) [X0 Z1]
+    0.2 * (X(0) @ Z(1)) + -0.543 * (Z(0) @ Hadamard(wires=[2]))
 
-    The user can also provide custom observables:
+    A ``qml.Hamiltonian`` stores information on which commuting observables should be measured
+    together in a circuit:
 
-    >>> obs_matrix = np.array([[0.5, 1.0j, 0.0, -3j],
-                               [-1.0j, -1.1, 0.0, -0.1],
-                               [0.0, 0.0, -0.9, 12.0],
-                               [3j, -0.1, 12.0, 0.0]])
-    >>> obs = qml.Hermitian(obs_matrix, wires=[0, 1])
-    >>> H = qml.Hamiltonian((0.8, ), (obs, ))
-    >>> print(H)
-    (0.8) [Hermitian0,1]
+    >>> obs = [qml.X(0), qml.X(1), qml.Z(0)]
+    >>> coeffs = np.array([1., 2., 3.])
+    >>> H = qml.Hamiltonian(coeffs, obs, grouping_type='qwc')
+    >>> H.grouping_indices
+    ((0, 1), (2,))
 
-    Alternatively, the :func:`~.molecular_hamiltonian` function from the
-    :doc:`/introduction/chemistry` module can be used to generate a molecular
-    Hamiltonian.
+    This attribute can be used to compute groups of coefficients and observables:
 
-    In many cases, Hamiltonians can be constructed using Pythonic arithmetic operations.
-    For example:
+    >>> grouped_coeffs = [coeffs[list(indices)] for indices in H.grouping_indices]
+    >>> grouped_obs = [[H.ops[i] for i in indices] for indices in H.grouping_indices]
+    >>> grouped_coeffs
+    [array([1., 2.]), array([3.])]
+    >>> grouped_obs
+    [[X(0), X(1)], [Z(0)]]
 
-    >>> qml.Hamiltonian([1.], [qml.X(0)]) + 2 * qml.Z(0) @ qml.Z(1)
+    Devices that evaluate a ``qml.Hamiltonian`` expectation by splitting it into its local
+    observables can use this information to reduce the number of circuits evaluated.
 
-    is equivalent to the following Hamiltonian:
+    Note that one can compute the ``grouping_indices`` for an already initialized ``qml.Hamiltonian``
+    by using the :func:`compute_grouping <pennylane.ops.LinearCombination.compute_grouping>` method.
 
-    >>> qml.Hamiltonian([1, 2], [qml.X(0), qml.Z(0) @ qml.Z(1)])
+    .. details::
+        :title: Old Hamiltonian behaviour
 
-    While scalar multiplication requires native python floats or integer types,
-    addition, subtraction, and tensor multiplication of Hamiltonians with Hamiltonians or
-    other observables is possible with tensor-valued coefficients, i.e.,
+        The following code examples show the behaviour of ``qml.Hamiltonian`` using old operator
+        arithmetic. See :doc:`Updated Operators </news/new_opmath/>` for more details. The old
+        behaviour can be reactivated by calling
 
-    >>> H1 = qml.Hamiltonian(torch.tensor([1.]), [qml.X(0)])
-    >>> H2 = qml.Hamiltonian(torch.tensor([2., 3.]), [qml.Y(0), qml.X(1)])
-    >>> obs3 = [qml.X(0), qml.Y(0), qml.X(1)]
-    >>> H3 = qml.Hamiltonian(torch.tensor([1., 2., 3.]), obs3)
-    >>> H3.compare(H1 + H2)
-    True
+        >>> qml.operation.disable_new_opmath()
 
-    A Hamiltonian can store information on which commuting observables should be measured together in
-    a circuit:
+        Alternatively, ``qml.ops.Hamiltonian`` provides a permanent access point for Hamiltonian
+        behaviour before ``v0.36``.
 
-    >>> obs = [qml.X(0), qml.X(1), qml.Z(0)]
-    >>> coeffs = np.array([1., 2., 3.])
-    >>> H = qml.Hamiltonian(coeffs, obs, grouping_type='qwc')
-    >>> H.grouping_indices
-    [[0, 1], [2]]
+        >>> coeffs = [0.2, -0.543]
+        >>> obs = [qml.X(0) @ qml.Z(1), qml.Z(0) @ qml.Hadamard(2)]
+        >>> H = qml.Hamiltonian(coeffs, obs)
+        >>> print(H)
+          (-0.543) [Z0 H2]
+        + (0.2) [X0 Z1]
 
-    This attribute can be used to compute groups of coefficients and observables:
+        The coefficients can be a trainable tensor, for example:
 
-    >>> grouped_coeffs = [coeffs[indices] for indices in H.grouping_indices]
-    >>> grouped_obs = [[H.ops[i] for i in indices] for indices in H.grouping_indices]
-    >>> grouped_coeffs
-    [tensor([1., 2.], requires_grad=True), tensor([3.], requires_grad=True)]
-    >>> grouped_obs
-    [[qml.X(0), qml.X(1)], [qml.Z(0)]]
+        >>> coeffs = tf.Variable([0.2, -0.543], dtype=tf.double)
+        >>> obs = [qml.X(0) @ qml.Z(1), qml.Z(0) @ qml.Hadamard(2)]
+        >>> H = qml.Hamiltonian(coeffs, obs)
+        >>> print(H)
+          (-0.543) [Z0 H2]
+        + (0.2) [X0 Z1]
+
+        The user can also provide custom observables:
+
+        >>> obs_matrix = np.array([[0.5, 1.0j, 0.0, -3j],
+                                   [-1.0j, -1.1, 0.0, -0.1],
+                                   [0.0, 0.0, -0.9, 12.0],
+                                   [3j, -0.1, 12.0, 0.0]])
+        >>> obs = qml.Hermitian(obs_matrix, wires=[0, 1])
+        >>> H = qml.Hamiltonian((0.8, ), (obs, ))
+        >>> print(H)
+        (0.8) [Hermitian0,1]
+
+        Alternatively, the :func:`~.molecular_hamiltonian` function from the
+        :doc:`/introduction/chemistry` module can be used to generate a molecular
+        Hamiltonian.
+
+        In many cases, Hamiltonians can be constructed using Pythonic arithmetic operations.
+        For example:
+
+        >>> qml.Hamiltonian([1.], [qml.X(0)]) + 2 * qml.Z(0) @ qml.Z(1)
+
+        is equivalent to the following Hamiltonian:
+
+        >>> qml.Hamiltonian([1, 2], [qml.X(0), qml.Z(0) @ qml.Z(1)])
+
+        While scalar multiplication requires native python floats or integer types,
+        addition, subtraction, and tensor multiplication of Hamiltonians with Hamiltonians or
+        other observables is possible with tensor-valued coefficients, i.e.,
+
+        >>> H1 = qml.Hamiltonian(torch.tensor([1.]), [qml.X(0)])
+        >>> H2 = qml.Hamiltonian(torch.tensor([2., 3.]), [qml.Y(0), qml.X(1)])
+        >>> obs3 = [qml.X(0), qml.Y(0), qml.X(1)]
+        >>> H3 = qml.Hamiltonian(torch.tensor([1., 2., 3.]), obs3)
+        >>> H3.compare(H1 + H2)
+        True
+
+        A Hamiltonian can store information on which commuting observables should be measured together in
+        a circuit:
+
+        >>> obs = [qml.X(0), qml.X(1), qml.Z(0)]
+        >>> coeffs = np.array([1., 2., 3.])
+        >>> H = qml.Hamiltonian(coeffs, obs, grouping_type='qwc')
+        >>> H.grouping_indices
+        [[0, 1], [2]]
+
+        This attribute can be used to compute groups of coefficients and observables:
+
+        >>> grouped_coeffs = [coeffs[indices] for indices in H.grouping_indices]
+        >>> grouped_obs = [[H.ops[i] for i in indices] for indices in H.grouping_indices]
+        >>> grouped_coeffs
+        [tensor([1., 2.], requires_grad=True), tensor([3.], requires_grad=True)]
+        >>> grouped_obs
+        [[qml.X(0), qml.X(1)], [qml.Z(0)]]
+
+        Devices that evaluate a Hamiltonian expectation by splitting it into its local observables can
+        use this information to reduce the number of circuits evaluated.
 
-    Devices that evaluate a Hamiltonian expectation by splitting it into its local observables can
-    use this information to reduce the number of circuits evaluated.
+        Note that one can compute the ``grouping_indices`` for an already initialized Hamiltonian by
+        using the :func:`compute_grouping <pennylane.Hamiltonian.compute_grouping>` method.
 
-    Note that one can compute the ``grouping_indices`` for an already initialized Hamiltonian by
-    using the :func:`compute_grouping <pennylane.Hamiltonian.compute_grouping>` method.
     """
 
     num_wires = qml.operation.AnyWires