Update demos to use new opmath (#1067)

Remove explicit `qml.operation.Tensor` instances, as well as
`qml.ops.Hamiltonian` usage

branching against `dev` as new opmath will be default only in `v0.36`
(scheduled for 07. May)

- [x] IMPORTANT change back final cell in `tutorial_qubit_tapering` when
PennyLaneAI/pennylane#5532 is fixed

_static/demonstration_assets/barren_gadgets/barren_gadgets.py

@@ -1,6 +1,8 @@
 import pennylane as qml
 from pennylane import numpy as np
 
+def non_identity_obs(obs):
+    return [o for o in obs if not isinstance(o, qml.Identity)]
 
 class PerturbativeGadgets:
     """ Class to generate the gadget Hamiltonian corresponding to a given
@@ -28,12 +30,13 @@ def gadgetize(self, Hamiltonian, target_locality=3):
         # checking for unaccounted for situations
         self.run_checks(Hamiltonian, target_locality)
         computational_qubits, computational_locality, computational_terms = self.get_params(Hamiltonian)
+        Hamiltonian_coeffs, Hamiltonian_ops = Hamiltonian.terms()
 
         # total qubit count, updated progressively when adding ancillaries
         total_qubits = computational_qubits
         #TODO: check proper convergence guarantee
         gap = 1
-        perturbation_norm = np.sum(np.abs(Hamiltonian.coeffs)) \
+        perturbation_norm = np.sum(np.abs(Hamiltonian_coeffs)) \
                           + computational_terms * (computational_locality - 1)
         lambda_max = gap / (4 * perturbation_norm)
         l = self.perturbation_factor * lambda_max
@@ -44,7 +47,7 @@ def gadgetize(self, Hamiltonian, target_locality=3):
         obs_anc = []
         obs_pert = []
         ancillary_register_size = int(computational_locality / (target_locality - 2))
-        for str_count, string in enumerate(Hamiltonian.ops):
+        for str_count, string in enumerate(Hamiltonian_ops):
             previous_total = total_qubits
             total_qubits += ancillary_register_size
             # Generating the ancillary part
@@ -54,10 +57,10 @@ def gadgetize(self, Hamiltonian, target_locality=3):
             # Generating the perturbative part
             for anc_q in range(ancillary_register_size):
                 term = qml.PauliX(previous_total+anc_q) @ qml.PauliX(previous_total+(anc_q+1)%ancillary_register_size)
-                term = qml.operation.Tensor(term, *string.non_identity_obs[
+                term = qml.prod(term, *non_identity_obs(string.operands)[
                     (target_locality-2)*anc_q:(target_locality-2)*(anc_q+1)])
                 obs_pert.append(term)
-            coeffs_pert += [l * sign_correction * Hamiltonian.coeffs[str_count]] \
+            coeffs_pert += [l * sign_correction * Hamiltonian_coeffs[str_count]] \
                         + [l] * (ancillary_register_size - 1)
         coeffs = coeffs_anc + coeffs_pert
         obs = obs_anc + obs_pert
@@ -77,12 +80,13 @@ def get_params(self, Hamiltonian):
             computational_terms (int)     : number of terms in the sum 
                                             composing the Hamiltonian
         """
+        _, Hamiltonian_ops = Hamiltonian.terms()
         # checking how many qubits the Hamiltonian acts on
         computational_qubits = len(Hamiltonian.wires)
         # getting the number of terms in the Hamiltonian
-        computational_terms = len(Hamiltonian.ops)
+        computational_terms = len(Hamiltonian_ops)
         # getting the locality, assuming all terms have the same
-        computational_locality = max([len(Hamiltonian.ops[s].non_identity_obs) 
+        computational_locality = max([len(non_identity_obs(Hamiltonian_ops[s])) 
                                       for s in range(computational_terms)])
         return computational_qubits, computational_locality, computational_terms
 
@@ -96,6 +100,7 @@ def run_checks(self, Hamiltonian, target_locality):
         Returns:
             None
         """
+        _, Hamiltonian_ops = Hamiltonian.terms()
         computational_qubits, computational_locality, _ = self.get_params(Hamiltonian)
         computational_qubits = len(Hamiltonian.wires)
         if computational_qubits != Hamiltonian.wires[-1] + 1:
@@ -104,8 +109,8 @@ def run_checks(self, Hamiltonian, target_locality):
                             'Decomposition not implemented for this case')
         # Check for same string lengths
         localities=[]
-        for string in Hamiltonian.ops:
-            localities.append(len(string.non_identity_obs))
+        for string in Hamiltonian_ops:
+            localities.append(len(non_identity_obs(string)))
         if len(np.unique(localities)) > 1:
             raise Exception('The given Hamiltonian has terms with different locality.' +
                             ' Gadgetization not implemented for this case')

demonstrations/barren_gadgets/barren_gadgets.py

@@ -1,6 +1,8 @@
 import pennylane as qml
 from pennylane import numpy as np
 
+def non_identity_obs(obs):
+    return [o for o in obs if not isinstance(o, qml.Identity)]
 
 class PerturbativeGadgets:
     """ Class to generate the gadget Hamiltonian corresponding to a given
@@ -28,12 +30,13 @@ def gadgetize(self, Hamiltonian, target_locality=3):
         # checking for unaccounted for situations
         self.run_checks(Hamiltonian, target_locality)
         computational_qubits, computational_locality, computational_terms = self.get_params(Hamiltonian)
+        Hamiltonian_coeffs, Hamiltonian_ops = Hamiltonian.terms()
 
         # total qubit count, updated progressively when adding ancillaries
         total_qubits = computational_qubits
         #TODO: check proper convergence guarantee
         gap = 1
-        perturbation_norm = np.sum(np.abs(Hamiltonian.coeffs)) \
+        perturbation_norm = np.sum(np.abs(Hamiltonian_coeffs)) \
                           + computational_terms * (computational_locality - 1)
         lambda_max = gap / (4 * perturbation_norm)
         l = self.perturbation_factor * lambda_max
@@ -44,7 +47,7 @@ def gadgetize(self, Hamiltonian, target_locality=3):
         obs_anc = []
         obs_pert = []
         ancillary_register_size = int(computational_locality / (target_locality - 2))
-        for str_count, string in enumerate(Hamiltonian.ops):
+        for str_count, string in enumerate(Hamiltonian_ops):
             previous_total = total_qubits
             total_qubits += ancillary_register_size
             # Generating the ancillary part
@@ -54,10 +57,10 @@ def gadgetize(self, Hamiltonian, target_locality=3):
             # Generating the perturbative part
             for anc_q in range(ancillary_register_size):
                 term = qml.PauliX(previous_total+anc_q) @ qml.PauliX(previous_total+(anc_q+1)%ancillary_register_size)
-                term = qml.operation.Tensor(term, *string.non_identity_obs[
+                term = qml.prod(term, *non_identity_obs(string.operands)[
                     (target_locality-2)*anc_q:(target_locality-2)*(anc_q+1)])
                 obs_pert.append(term)
-            coeffs_pert += [l * sign_correction * Hamiltonian.coeffs[str_count]] \
+            coeffs_pert += [l * sign_correction * Hamiltonian_coeffs[str_count]] \
                         + [l] * (ancillary_register_size - 1)
         coeffs = coeffs_anc + coeffs_pert
         obs = obs_anc + obs_pert
@@ -77,12 +80,13 @@ def get_params(self, Hamiltonian):
             computational_terms (int)     : number of terms in the sum 
                                             composing the Hamiltonian
         """
+        _, Hamiltonian_ops = Hamiltonian.terms()
         # checking how many qubits the Hamiltonian acts on
         computational_qubits = len(Hamiltonian.wires)
         # getting the number of terms in the Hamiltonian
-        computational_terms = len(Hamiltonian.ops)
+        computational_terms = len(Hamiltonian_ops)
         # getting the locality, assuming all terms have the same
-        computational_locality = max([len(Hamiltonian.ops[s].non_identity_obs) 
+        computational_locality = max([len(non_identity_obs(Hamiltonian_ops[s])) 
                                       for s in range(computational_terms)])
         return computational_qubits, computational_locality, computational_terms
 
@@ -96,6 +100,7 @@ def run_checks(self, Hamiltonian, target_locality):
         Returns:
             None
         """
+        _, Hamiltonian_ops = Hamiltonian.terms()
         computational_qubits, computational_locality, _ = self.get_params(Hamiltonian)
         computational_qubits = len(Hamiltonian.wires)
         if computational_qubits != Hamiltonian.wires[-1] + 1:
@@ -104,8 +109,8 @@ def run_checks(self, Hamiltonian, target_locality):
                             'Decomposition not implemented for this case')
         # Check for same string lengths
         localities=[]
-        for string in Hamiltonian.ops:
-            localities.append(len(string.non_identity_obs))
+        for string in Hamiltonian_ops:
+            localities.append(len(non_identity_obs(string)))
         if len(np.unique(localities)) > 1:
             raise Exception('The given Hamiltonian has terms with different locality.' +
                             ' Gadgetization not implemented for this case')

demonstrations/braket-parallel-gradients.py

@@ -138,7 +138,7 @@ def circuit(params):
 
     # Measure all qubits to make sure all's good with Braket
     observables = [qml.PauliZ(n_wires - 1)] + [qml.Identity(i) for i in range(n_wires - 1)]
-    return qml.expval(qml.operation.Tensor(*observables))
+    return qml.expval(qml.prod(*observables))
 
 
 ##############################################################################

demonstrations/tutorial_barren_gadgets.py

@@ -202,7 +202,7 @@
 
 gadgetizer = PerturbativeGadgets()
 H_gadget = gadgetizer.gadgetize(H_target)
-print(H_gadget)
+H_gadget
 
 ##############################################################################
 # So, let's see what we got.

demonstrations/tutorial_classical_shadows.py

@@ -512,8 +512,8 @@ def estimate_shadow_obervable(shadow, observable, k=10):
         ), np.array([observable.wires[0]])
     else:
         target_obs, target_locs = np.array(
-            [map_name_to_int[o.name] for o in observable.obs]
-        ), np.array([o.wires[0] for o in observable.obs])
+            [map_name_to_int[o.name] for o in observable.operands]
+        ), np.array([o.wires[0] for o in observable.operands])
 
     # classical values
     b_lists, obs_lists = shadow

demonstrations/tutorial_clifford_circuit_simulations.py

@@ -332,8 +332,9 @@ def clifford_tableau(op):
     for pauli in pauli_ops:
         conjugate = qml.prod(qml.adjoint(op), pauli, op).simplify()
         decompose = qml.pauli_decompose(conjugate.matrix(), wire_order=op.wires)
-        phase = "+" if list(decompose.coeffs)[0] >= 0 else "-"
-        print(pauli, "-—>", phase, list(decompose.ops)[0])
+        decompose_coeffs, decompose_ops = decompose.terms()
+        phase = "+" if list(decompose_coeffs)[0] >= 0 else "-"
+        print(pauli, "-—>", phase, list(decompose_ops)[0])
 
 clifford_tableau(qml.X(0))
 

demonstrations/tutorial_diffable_shadows.py

@@ -376,7 +376,7 @@ def qnode_shadow():
     method="pyscf",
 )
 
-coeffs, obs = H.coeffs, H.ops
+coeffs, obs = H.terms()
 H_qwc = qml.Hamiltonian(coeffs, obs, grouping_type="qwc")
 
 groups = qml.pauli.group_observables(obs)
@@ -421,6 +421,7 @@ def circuit():
 
         # execute qwc measurements
         dev_finite = qml.device("default.qubit", wires=range(n_wires), shots=int(shots))
+
         @qml.qnode(dev_finite, interface="autograd")
         def qnode_finite(H):
             circuit()

demonstrations/tutorial_implicit_diff_susceptibility.py

@@ -449,11 +449,8 @@ def energy(z, a):
         float: The expectation value (energy).
     """
     variational_ansatz(*z, wires=range(N))
-    # here we compute the Hamiltonian coefficients and operations
-    # 'by hand' because the qml.Hamiltonian class does not support
-    # operator arithmetic with JAX device arrays.
-    coeffs = jnp.concatenate([H0.coeffs, a * A.coeffs])
-    return qml.expval(qml.Hamiltonian(coeffs, H0.ops + A.ops))
+
+    return qml.expval(H0 + a * A)
 
 
 z_init = [jnp.array(2 * np.pi * np.random.random(s)) for s in weights_shape]

demonstrations/tutorial_kernel_based_training.py

@@ -172,7 +172,6 @@
 
 import pennylane as qml
 from pennylane.templates import AngleEmbedding, StronglyEntanglingLayers
-from pennylane.operation import Tensor
 
 import matplotlib.pyplot as plt
 

demonstrations/tutorial_lcu_blockencoding.py

@@ -59,10 +59,11 @@
 )
 
 LCU = qml.pauli_decompose(A)
+LCU_coeffs, LCU_ops = LCU.terms()
 
 print(f"LCU decomposition:\n {LCU}")
-print(f"Coefficients:\n {LCU.coeffs}")
-print(f"Unitaries:\n {LCU.ops}")
+print(f"Coefficients:\n {LCU_coeffs}")
+print(f"Unitaries:\n {LCU_ops}")
 
 
 ##############################################################################
@@ -144,7 +145,7 @@
 dev1 = qml.device("default.qubit", wires=1)
 
 # normalized square roots of coefficients
-alphas = (np.sqrt(LCU.coeffs) / np.linalg.norm(np.sqrt(LCU.coeffs)))
+alphas = (np.sqrt(LCU_coeffs) / np.linalg.norm(np.sqrt(LCU_coeffs)))
 
 
 @qml.qnode(dev1)
@@ -167,7 +168,7 @@ def prep_circuit():
 dev2 = qml.device("default.qubit", wires=3)
 
 # unitaries
-ops = LCU.ops
+ops = LCU_ops
 # relabeling wires: 0 → 1, and 1 → 2
 unitaries = [qml.map_wires(op, {0: 1, 1: 2}) for op in ops]
 

demonstrations/tutorial_liealgebra.py

@@ -101,7 +101,6 @@
 import numpy as np
 import pennylane as qml
 from pennylane import X, Y, Z
-qml.operation.enable_new_opmath()
 
 su2 = [1j * X(0), 1j * Y(0), 1j * Z(0)]
 
@@ -338,8 +337,6 @@
 print(qml.commutator(SZ, SX) == (2j*SY).simplify())
 print(qml.commutator(SY, SZ) == (2j*SX).simplify())
 
-qml.operation.disable_new_opmath()
-
 ##############################################################################
 #
 # Another perspective on the inherent :math:`SU(2)` symmetry of :math:`H_\text{Heis}` is that the expectation 

demonstrations/tutorial_measurement_optimize.py

@@ -164,7 +164,7 @@ def cost_circuit(params):
 H, num_qubits = dataset.hamiltonian, len(dataset.hamiltonian.wires)
 
 print("Required number of qubits:", num_qubits)
-print("Number of Hamiltonian terms/required measurements:", len(H.ops))
+print("Number of Hamiltonian terms/required measurements:", len(H.terms()[0]))
 
 print("\n", H)
 
@@ -572,8 +572,8 @@ def circuit(weights):
 def format_pauli_word(term):
     """Convenience function that nicely formats a PennyLane
     tensor observable as a Pauli word"""
-    if isinstance(term, qml.operation.Tensor):
-        return " ".join([format_pauli_word(t) for t in term.obs])
+    if isinstance(term, qml.ops.Prod):
+        return " ".join([format_pauli_word(t) for t in term])
 
     return f"{term.name[-1]}{term.wires.tolist()[0]}"
 
@@ -738,6 +738,7 @@ def circuit(weights, group=None, **kwargs):
 # automatically optimize the measurements.
 
 H = qml.Hamiltonian(coeffs=np.ones(len(terms)), observables=terms, grouping_type="qwc")
+_, H_ops = H.terms()
 @qml.qnode(dev, interface="autograd")
 def cost_fn(weights):
     qml.StronglyEntanglingLayers(weights, wires=range(4))
@@ -754,10 +755,10 @@ def cost_fn(weights):
 
 dataset = qml.data.load('qchem', molname="H2O")[0]
 H, num_qubits = dataset.hamiltonian, len(dataset.hamiltonian.wires)
-print("Number of Hamiltonian terms/required measurements:", len(H.ops))
+print("Number of Hamiltonian terms/required measurements:", len(H_ops))
 
 # grouping
-groups = qml.pauli.group_observables(H.ops, grouping_type='qwc', method='rlf')
+groups = qml.pauli.group_observables(H_ops, grouping_type='qwc', method='rlf')
 print("Number of required measurements after optimization:", len(groups))
 
 ##############################################################################

demonstrations/tutorial_optimal_control.py

@@ -565,7 +565,6 @@ def run_adam(profit_fn, grad_fn, params, learning_rate, num_steps):
 colors = {0: "#70CEFF", 1: "#C756B2", 2: "#FDC357"}
 dashes = {"X": [10, 0], "Y": [2, 2, 10, 2], "Z": [6, 2]}
 
-
 def plot_optimal_pulses(hist, pulse_fn, ops, T, target_name):
     _, profit_hist = list(zip(*hist))
     fig, axs = plt.subplots(2, 1, figsize=(10, 9), gridspec_kw={"hspace": 0.0}, sharex=True)
@@ -574,18 +573,11 @@ def plot_optimal_pulses(hist, pulse_fn, ops, T, target_name):
     max_params, max_profit = hist[jnp.argmax(jnp.array(profit_hist))]
     plot_times = jnp.linspace(0, T, 300)
     # Iterate over pulse parameters and parametrized operators
-    for p, op in zip(max_params, ops):
+    for i, (p, op) in enumerate(zip(max_params, ops)):
         # Create label, and pick correct axis
-        label = op.name
-        ax = axs[0] if isinstance(label, str) else axs[1]
-        # Convert the label into a concise string. This differs depending on
-        # whether the operator has a single or multiple Pauli terms. Pick the line style
-        if isinstance(label, str):
-            label = f"${label[-1]}_{op.wires[0]}$"
-            dash = dashes[label[1]]
-        else:
-            label = "$" + " ".join([f"{n[-1]}_{w}" for w, n in zip(op.wires, label)]) + "$"
-            dash = [10, 0]
+        label = str(op)
+        dash = dashes[label[0]]
+        ax = axs[0] if len(op.wires) == 1 else axs[1]
 
         # Set color according to qubit the term acts on
         col = colors[op.wires[0]]

demonstrations/tutorial_pulse_programming101.py

@@ -224,9 +224,10 @@ def qnode(params):
 
 data = qml.data.load("qchem", molname="HeH+", basis="STO-3G", bondlength=1.5)[0]
 H_obj = data.tapered_hamiltonian
+H_obj_coeffs, H_obj_ops = H_obj.terms()
 
 # casting the Hamiltonian coefficients to a jax Array
-H_obj = qml.Hamiltonian(jnp.array(H_obj.coeffs), H_obj.ops)
+H_obj = qml.Hamiltonian(jnp.array(H_obj_coeffs), H_obj_ops)
 E_exact = data.fci_energy
 n_wires = len(H_obj.wires)