[new opmath 4] Update tests to pass under new opmath

doc/development/deprecations.rst

@@ -34,6 +34,20 @@ Pending deprecations
 
   - Added and deprecated for ``Sum`` and ``Prod`` instances in v0.35
 
+* Accessing ``qml.ops.Hamiltonian`` with new operator arithmetic enabled is deprecated. Using ``qml.Hamiltonian``
+  with new operator arithmetic enabled now returns a ``LinearCombination`` instance. Some functionality
+  may not work as expected, and use of the Hamiltonian class with the new operator arithmetic will not
+  be supported in future releases of PennyLane.
+
+  You can update your code to the new operator arithmetic by using ``qml.Hamiltonian`` instead of importing
+  the Hamiltonian class directly or via ``qml.ops.Hamiltonian``. When the new operator arithmetic is enabled, 
+  ``qml.Hamiltonian`` will access the new corresponding implementation. 
+
+  Alternatively, to continue accessing the legacy functionality, you can use 
+  ``qml.operation.disable_new_opmath()``.
+
+  - Deprecated in v0.36
+
 Completed deprecation cycles
 ----------------------------
 

doc/releases/changelog-dev.md

@@ -64,6 +64,14 @@
 * Added new function `qml.operation.convert_to_legacy_H` to convert `Sum`, `SProd`, and `Prod` to `Hamiltonian` instances.
   [(#5309)](https://github.com/PennyLaneAI/pennylane/pull/5309)
 
+<h3>Improvements 🛠</h3>
+
+* The `qml.is_commuting` function now accepts `Sum`, `SProd`, and `Prod` instances.
+  [(#5351)](https://github.com/PennyLaneAI/pennylane/pull/5351)
+
+* Operators can now be left multiplied `x * op` by numpy arrays.
+  [(#5361)](https://github.com/PennyLaneAI/pennylane/pull/5361)
+
 * Create the `qml.Reflection` operator, useful for amplitude amplification and its variants.
   [(#5159)](https://github.com/PennyLaneAI/pennylane/pull/5159)
 
@@ -151,9 +159,14 @@
 * The `molecular_hamiltonian` function calls `PySCF` directly when `method='pyscf'` is selected.
   [(#5118)](https://github.com/PennyLaneAI/pennylane/pull/5118)
 
-* All generators in the source code (except those in the `qchem` module) no longer return
-  `Hamiltonian` or `Tensor` instances. Wherever possible, these return `Sum`, `SProd`, and `Prod` instances.
+* The generators in the source code return operators consistent with the global setting for 
+  `qml.operator.active_new_opmath()` wherever possible. `Sum`, `SProd` and `Prod` instances 
+  will be returned even after disabling the new operator arithmetic in cases where they offer 
+  additional functionality not available using legacy operators.
   [(#5253)](https://github.com/PennyLaneAI/pennylane/pull/5253)
+  [(#5410)](https://github.com/PennyLaneAI/pennylane/pull/5410)
+  [(#5411)](https://github.com/PennyLaneAI/pennylane/pull/5411) 
+  [(#5421)](https://github.com/PennyLaneAI/pennylane/pull/5421)
 
 * Upgraded `null.qubit` to the new device API. Also, added support for all measurements and various modes of differentiation.
   [(#5211)](https://github.com/PennyLaneAI/pennylane/pull/5211)
@@ -164,6 +177,9 @@
 * `Hamiltonian.pauli_rep` is now defined if the hamiltonian is a linear combination of paulis.
   [(#5377)](https://github.com/PennyLaneAI/pennylane/pull/5377)
 
+* `Prod.eigvals()` is now compatible with Qudit operators.
+  [(#5400)](https://github.com/PennyLaneAI/pennylane/pull/5400)
+
 * Obtaining classical shadows using the `default.clifford` device is now compatible with
   [stim](https://github.com/quantumlib/Stim) `v1.13.0`.
   [(#5409)](https://github.com/PennyLaneAI/pennylane/pull/5409)
@@ -230,6 +246,10 @@
 * Attempting to multiply `PauliWord` and `PauliSentence` with `*` will raise an error. Instead, use `@` to conform with the PennyLane convention.
   [(#5341)](https://github.com/PennyLaneAI/pennylane/pull/5341)
 
+* When new operator arithmetic is enabled, `qml.Hamiltonian` is now an alias for `qml.ops.LinearCombination`.
+  `Hamiltonian` will still be accessible as `qml.ops.Hamiltonian`.
+  [(#5393)](https://github.com/PennyLaneAI/pennylane/pull/5393)
+
 * Since `default.mixed` does not support snapshots with measurements, attempting to do so will result in a `DeviceError` instead of getting the density matrix.
   [(#5416)](https://github.com/PennyLaneAI/pennylane/pull/5416)
 
@@ -249,6 +269,11 @@
   ...     circuit = qml.from_qasm(f.read())
   ```
 
+* Accessing `qml.ops.Hamiltonian` with new operator arithmetic is deprecated. Using `qml.Hamiltonian` with new operator arithmetic enabled now
+  returns a `LinearCombination` instance. Some functionality may not work as expected. To continue using the `Hamiltonian` class, you can use
+  `qml.operation.disable_new_opmath()` to disable the new operator arithmetic.
+  [(#5393)](https://github.com/PennyLaneAI/pennylane/pull/5393)
+
 <h3>Documentation 📝</h3>
 
 * Removed some redundant documentation for the `evolve` function.

pennylane/__init__.py

@@ -145,6 +145,18 @@ class PennyLaneDeprecationWarning(UserWarning):
     """Warning raised when a PennyLane feature is being deprecated."""
 
 
+del globals()["Hamiltonian"]
+
+
+def __getattr__(name):
+    if name == "Hamiltonian":
+        if pennylane.operation.active_new_opmath():
+            return pennylane.ops.LinearCombination
+        return pennylane.ops.Hamiltonian
+
+    raise AttributeError(f"module 'pennylane' has no attribute '{name}'")
+
+
 def _get_device_entrypoints():
     """Returns a dictionary mapping the device short name to the
     loadable entrypoint"""

pennylane/_device.py

@@ -680,7 +680,6 @@ def default_expand_fn(self, circuit, max_expansion=10):
         obs_on_same_wire &= not any(
             isinstance(o, (Hamiltonian, LinearCombination)) for o in circuit._obs_sharing_wires
         )
-
         ops_not_supported = not all(self.stopping_condition(op) for op in circuit.operations)
 
         if obs_on_same_wire:

pennylane/data/attributes/operator/operator.py

@@ -61,7 +61,13 @@ def consumes_types(cls) -> FrozenSet[Type[Operator]]:
                 qml.QubitCarry,
                 qml.QubitSum,
                 # pennylane/ops/qubit/hamiltonian.py
-                qml.Hamiltonian,
+                qml.ops.Hamiltonian,
+                # pennylane/ops/op_math/linear_combination.py
+                qml.ops.LinearCombination,
+                # pennylane/ops/op_math - prod.py, s_prod.py, sum.py
+                qml.ops.Prod,
+                qml.ops.SProd,
+                qml.ops.Sum,
                 # pennylane/ops/qubit/matrix_qml.py
                 qml.QubitUnitary,
                 qml.DiagonalQubitUnitary,
@@ -206,6 +212,8 @@ def _ops_to_hdf5(
         op_class_names = []
         for i, op in enumerate(value):
             op_key = f"op_{i}"
+            if isinstance(op, (qml.ops.Prod, qml.ops.SProd, qml.ops.Sum)):
+                op = op.simplify()
             if type(op) not in self.consumes_types():
                 raise TypeError(
                     f"Serialization of operator type '{type(op).__name__}' is not supported."
@@ -214,11 +222,19 @@ def _ops_to_hdf5(
             if isinstance(op, Tensor):
                 self._ops_to_hdf5(bind, op_key, op.obs)
                 op_wire_labels.append("null")
-            elif isinstance(op, qml.Hamiltonian):
+            elif isinstance(op, (qml.ops.Hamiltonian, qml.ops.LinearCombination)):
                 coeffs, ops = op.terms()
                 ham_grp = self._ops_to_hdf5(bind, op_key, ops)
                 ham_grp["hamiltonian_coeffs"] = coeffs
                 op_wire_labels.append("null")
+            elif isinstance(op, (qml.ops.Prod, qml.ops.Sum)):
+                self._ops_to_hdf5(bind, op_key, op.operands)
+                op_wire_labels.append("null")
+            elif isinstance(op, qml.ops.SProd):
+                coeffs, ops = op.terms()
+                sprod_grp = self._ops_to_hdf5(bind, op_key, ops)
+                sprod_grp["sprod_scalar"] = coeffs
+                op_wire_labels.append("null")
             else:
                 bind[op_key] = op.data if len(op.data) else h5py.Empty("f")
                 op_wire_labels.append(wires_to_json(op.wires))
@@ -238,21 +254,29 @@ def _hdf5_to_ops(self, bind: HDF5Group) -> List[Operator]:
         wires_bind = bind["op_wire_labels"]
         op_class_names = [] if names_bind.shape == (0,) else names_bind.asstr()
         op_wire_labels = [] if wires_bind.shape == (0,) else wires_bind.asstr()
-
         with qml.QueuingManager.stop_recording():
             for i, op_class_name in enumerate(op_class_names):
                 op_key = f"op_{i}"
 
                 op_cls = self._supported_ops_dict()[op_class_name]
                 if op_cls is Tensor:
                     ops.append(Tensor(*self._hdf5_to_ops(bind[op_key])))
-                elif op_cls is qml.Hamiltonian:
+                elif op_cls in (qml.ops.Hamiltonian, qml.ops.LinearCombination):
                     ops.append(
                         qml.Hamiltonian(
                             coeffs=list(bind[op_key]["hamiltonian_coeffs"]),
                             observables=self._hdf5_to_ops(bind[op_key]),
                         )
                     )
+                elif op_cls in (qml.ops.Prod, qml.ops.Sum):
+                    ops.append(op_cls(*self._hdf5_to_ops(bind[op_key])))
+                elif op_cls is qml.ops.SProd:
+                    ops.append(
+                        qml.ops.s_prod(
+                            scalar=bind[op_key]["sprod_scalar"][0],
+                            operator=self._hdf5_to_ops(bind[op_key])[0],
+                        )
+                    )
                 else:
                     wire_labels = json.loads(op_wire_labels[i])
                     op_data = bind[op_key]

pennylane/devices/default_clifford.py

@@ -1018,9 +1018,9 @@ def _measure_single_sample(stim_ct, meas_ops, meas_idx, meas_wire):
         """Sample a single qubit Pauli measurement from a stim circuit"""
         stim_sm = stim.TableauSimulator()
         stim_sm.do_circuit(stim_ct)
-        return stim_sm.measure_observable(
-            stim.PauliString([0] * meas_idx + meas_ops + [0] * (meas_wire - meas_idx - 1))
-        )
+        res = [0] * meas_idx + meas_ops + [0] * (meas_wire - meas_idx - 1)
+        res = [int(r) for r in res]
+        return stim_sm.measure_observable(stim.PauliString(res))
 
     def _sample_classical_shadow(self, meas, stim_circuit, shots, seed):
         """Measures classical shadows from the state of simulator device"""

pennylane/devices/default_qutrit.py

@@ -33,6 +33,7 @@
 OMEGA = qml.math.exp(2 * np.pi * 1j / 3)
 
 
+# pylint: disable=too-many-arguments
 class DefaultQutrit(QutritDevice):
     """Default qutrit device for PennyLane.
 
@@ -87,11 +88,7 @@ class DefaultQutrit(QutritDevice):
 
     # Identity is supported as an observable for qml.state() to work correctly. However, any
     # measurement types that rely on eigenvalue decomposition will not work with qml.Identity
-    observables = {
-        "THermitian",
-        "GellMann",
-        "Identity",
-    }
+    observables = {"THermitian", "GellMann", "Identity", "Prod"}
 
     # Static methods to use qml.math to allow for backprop differentiation
     _reshape = staticmethod(qml.math.reshape)

pennylane/devices/qubit/measure.py

@@ -194,7 +194,7 @@ def get_measurement_function(
                 return full_dot_products
 
             backprop_mode = math.get_interface(state, *measurementprocess.obs.data) != "numpy"
-            if isinstance(measurementprocess.obs, (LinearCombination, Hamiltonian)):
+            if isinstance(measurementprocess.obs, (Hamiltonian, LinearCombination)):
                 # need to work out thresholds for when its faster to use "backprop mode" measurements
                 return sum_of_terms_method if backprop_mode else csr_dot_products
 

pennylane/devices/qubit/sampling.py

@@ -59,7 +59,7 @@ def _group_measurements(mps: List[Union[SampleMeasurement, ClassicalShadowMP, Sh
         elif mp.obs is None:
             mp_no_obs.append(mp)
             mp_no_obs_indices.append(i)
-        elif isinstance(mp.obs, (Sum, Hamiltonian, LinearCombination, SProd, Prod)):
+        elif isinstance(mp.obs, (Hamiltonian, Sum, SProd, Prod)):
             # Sum, Hamiltonian, SProd, and Prod are treated as valid Pauli words, but
             # aren't accepted in qml.pauli.group_observables
             mp_other_obs.append([mp])

pennylane/devices/qutrit_mixed/measure.py

@@ -18,7 +18,7 @@
 from typing import Callable
 from string import ascii_letters as alphabet
 from pennylane import math
-from pennylane.ops import Sum, Hamiltonian, LinearCombination
+from pennylane.ops import Sum, Hamiltonian
 from pennylane.measurements import (
     StateMeasurement,
     MeasurementProcess,
@@ -262,7 +262,7 @@ def get_measurement_function(
         if isinstance(measurementprocess, ExpectationMP):
             # TODO add faster methods
             # TODO add support for sparce Hamiltonians
-            if isinstance(measurementprocess.obs, (Hamiltonian, Sum, LinearCombination)):
+            if isinstance(measurementprocess.obs, (Hamiltonian, Sum)):
                 return calculate_expval_sum_of_terms
             if measurementprocess.obs.has_matrix:
                 return calculate_expval

pennylane/devices/qutrit_mixed/sampling.py

@@ -18,7 +18,7 @@
 import numpy as np
 import pennylane as qml
 from pennylane import math
-from pennylane.ops import Sum, Hamiltonian
+from pennylane.ops import Sum
 from pennylane.measurements import (
     Shots,
     SampleMeasurement,
@@ -203,13 +203,15 @@ def _sum_for_single_shot(s):
                 )
             )
 
-        if isinstance(mp.obs, Hamiltonian):
+        if isinstance(mp.obs, qml.ops.Hamiltonian):
             # If Hamiltonian apply coefficients
             return sum((c * res for c, res in zip(mp.obs.terms()[0], results)))
+
         return sum(results)
 
     if shots.has_partitioned_shots:
         return tuple(_sum_for_single_shot(type(shots)(s)) for s in shots)
+
     return _sum_for_single_shot(shots)
 
 
@@ -350,7 +352,7 @@ def measure_with_samples(
         TensorLike[Any]: Sample measurement results
     """
 
-    if isinstance(mp, ExpectationMP) and isinstance(mp.obs, (Hamiltonian, Sum)):
+    if isinstance(mp, ExpectationMP) and isinstance(mp.obs, (qml.ops.Hamiltonian, Sum)):
         measure_fn = _measure_sum_with_samples
     else:
         # measure with the usual method (rotate into the measurement basis)

pennylane/gradients/gradient_transform.py

@@ -166,7 +166,7 @@ def _try_zero_grad_from_graph_or_get_grad_method(tape, param_index, use_graph=Tr
             # there is no influence of this operation on any of the observables
             return "0"
 
-    return par_info["op"].grad_method
+    return getattr(par_info["op"], "grad_method", None)
 
 
 def _find_gradient_methods(tape, trainable_param_indices, use_graph=True):

pennylane/gradients/hadamard_gradient.py

@@ -329,7 +329,8 @@ def _expval_hadamard_grad(tape, argnum, aux_wire):
                     obs_new = [qml.Z(i) for i in m.wires]
 
                 obs_new.append(qml.Y(aux_wire))
-                obs_new = qml.operation.Tensor(*obs_new)
+                obs_type = qml.prod if qml.operation.active_new_opmath() else qml.operation.Tensor
+                obs_new = obs_type(*obs_new)
 
                 if isinstance(m, qml.measurements.ExpectationMP):
                     measurements.append(qml.expval(op=obs_new))
@@ -460,9 +461,6 @@ def _get_generators(trainable_op):
     elif isinstance(trainable_op, qml.Rot):
         generators = [qml.Z(trainable_op.wires)]
         coeffs = [-0.5]
-    elif isinstance(trainable_op, (qml.RX, qml.RY, qml.RZ)):
-        generators = [trainable_op.generator().base]
-        coeffs = [trainable_op.generator().scalar]
     else:
         generators = trainable_op.generator().ops
         coeffs = trainable_op.generator().coeffs

pennylane/gradients/hamiltonian_grad.py

@@ -30,7 +30,13 @@ def hamiltonian_grad(tape, idx):
     # get position in queue
     queue_position = m_pos - len(tape.operations)
     new_measurements = list(tape.measurements)
-    new_measurements[queue_position] = qml.expval(op.ops[p_idx])
+
+    new_parameters = [0 * d for d in op.data]
+    new_parameters[p_idx] = qml.math.ones_like(op.data[p_idx])
+    new_obs = qml.ops.functions.bind_new_parameters(op, new_parameters)
+    new_obs = qml.simplify(new_obs)
+
+    new_measurements[queue_position] = qml.expval(new_obs)
 
     new_tape = qml.tape.QuantumScript(tape.operations, new_measurements, shots=tape.shots)
 
@@ -40,9 +46,7 @@ def processing_fn(results):
             res = results[0][queue_position]
             zeros = qml.math.zeros_like(res)
 
-            final = []
-            for i, _ in enumerate(tape.measurements):
-                final.append(res if i == queue_position else zeros)
+            final = [res if i == queue_position else zeros for i, _ in enumerate(tape.measurements)]
 
             return qml.math.expand_dims(qml.math.stack(final), 0)