Complete first section of Pauli Twirling user guide

docs/source/guide/pt-1-intro.md

@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.11.1
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -26,20 +26,18 @@ import mitiq
 mitiq.SUPPORTED_PROGRAM_TYPES.keys()
 ```
 
-
 ## Problem setup
-We first define the circuit of interest. In this example, the circuit has 
-two CNOT gates and a CZ gate. We can see that when we apply Pauli Twirling,
-we will generate 
+We first define the circuit of interest, which contains Hadamard, C-NOT, and C-Z gates.
 
 ```{code-cell} ipython3
-from cirq import LineQubit, Circuit, CZ, CNOT
+from cirq import LineQubit, Circuit, CZ, CNOT, H
 
-a, b, c, d  = LineQubit.range(4)
+q0, q1, q2, q3  = LineQubit.range(4)
 circuit = Circuit(
-    CNOT.on(a, b),
-    CZ.on(b, c),
-    CNOT.on(c, d),
+    H(q0),
+    CNOT.on(q0, q1),
+    CZ.on(q1, q2),
+    CNOT.on(q2, q3),
 )
 
 print(circuit)
@@ -50,45 +48,73 @@ the circuit on a noisy simulator, and returns the probability of the ground
 state. See the [Executors](executors.md) section for more information on
 how to define more advanced executors.
 
+As per the noise model that is applied during execution by the simulator, 
+we choose a depolarizing channel applied to the output of each 2-qubit gates.
+Admittedly, this is not necessarily the most realistic noise,
+but for the sake of this tutorial, it is going to be useful for highlighting 
+the effect of Pauli Twirling.  
+
 ```{code-cell} ipython3
-import numpy as np
-from cirq import DensityMatrixSimulator, amplitude_damp
-from mitiq.interface import convert_to_mitiq
+from cirq import CZPowGate, CXPowGate, CircuitOperation, depolarize, DensityMatrixSimulator
+from cirq.devices.noise_model import GateSubstitutionNoiseModel
 
-def execute(circuit, noise_level=0.1):
-    """Returns Tr[ρ |0⟩⟨0|] where ρ is the state prepared by the circuit
-    executed with amplitude damping noise.
+def get_noise_model(noise_level: float) -> GateSubstitutionNoiseModel:
+    """Substitute each C-Z and C-NOT gate in the circuit 
+    with the gate itself followed by a depolarizing channel
     """
-    # Replace with code based on your frontend and backend.
-    mitiq_circuit, _ = convert_to_mitiq(circuit)
-    noisy_circuit = mitiq_circuit.with_noise(amplitude_damp(gamma=noise_level))
-    rho = DensityMatrixSimulator().simulate(noisy_circuit).final_density_matrix
-    return rho[0, 0].real
+    def noisy_c_gate(op):
+        if isinstance(op.gate, (CZPowGate, CXPowGate)):
+            return CircuitOperation(
+                Circuit(
+                    op.gate.on(*op.qubits), 
+                    depolarize(p=noise_level, n_qubits=2).on_each(op.qubits)
+                ).freeze())
+        return op
+
+    return GateSubstitutionNoiseModel(noisy_c_gate)
+
+def execute(circuit: Circuit, noise_level: float):
+    """Returns Tr[ρ |0⟩⟨0|] where ρ is the state prepared by the circuit"""
+    return (
+        DensityMatrixSimulator(noise=get_noise_model(noise_level=noise_level))
+        .simulate(circuit)
+        .final_density_matrix[0, 0]
+        .real
+    )
 ```
 
 The [executor](executors.md) can be used to evaluate noisy (unmitigated)
 expectation values.
 
 ```{code-cell} ipython3
+# Set the parameter of the depolarizing channel
+NOISE_LEVEL = 0.1
+
 # Compute the expectation value of the |0><0| observable.
-noisy_value = execute(circuit)
 ideal_value = execute(circuit, noise_level=0.0)
-print(f"Error without mitigation: {abs(ideal_value - noisy_value) :.3}")
+noisy_value = execute(circuit, noise_level=NOISE_LEVEL)
+
+print(f"Error without twirling: {abs(ideal_value - noisy_value) :.3}")
 ```
 
 ## Apply PT
-Pauli Twirling can be easily implemented with the function
-{func}`.pauli_twirl_circuit()`.
+Pauli Twirling can be applied with the function
+{func}`.pauli_twirl_circuit()` from the `mitiq.pt` module.
 
 ```{code-cell} ipython3
-from mitiq import pt
-mitigated_result = pt.pauli_twirl_circuit(
-    circuit=circuit,
-)
-```
+from functools import partial
+import numpy as np
+from mitiq.executor.executor import Executor
+from mitiq.pt import generate_pauli_twirl_variants
 
-```{code-cell} ipython3
-# print(f"Error with mitigation (PT): {abs(ideal_value - mitigated_result) :.3}")
+# Generate twirled circuits
+twirled_circuits = generate_pauli_twirl_variants(circuit)
+
+# Average results executed over twirled circuits
+pt_vals = Executor(partial(execute, noise_level=NOISE_LEVEL)).evaluate(twirled_circuits)
+mitigated_result = np.average(pt_vals)
+
+print(f"Error with mitigation (PT): {abs(ideal_value - mitigated_result) :.3}")
 ```
 
 Here we observe that the application of PT does not reduce the estimation error when compared

mitiq/pt/__init__.py

@@ -4,7 +4,7 @@
 # LICENSE file in the root directory of this source tree.
 
 from mitiq.pt.pt import (
-    pauli_twirl_circuit,
+    generate_pauli_twirl_variants,
     twirl_CNOT_gates,
     twirl_CZ_gates,
 )

mitiq/pt/pt.py

@@ -66,7 +66,7 @@
 }
 
 
-def pauli_twirl_circuit(
+def generate_pauli_twirl_variants(
     circuit: QPROGRAM,
     num_circuits: int = 10,
     noise_name: Optional[str] = None,
@@ -91,7 +91,7 @@ def pauli_twirl_circuit(
         noise on these gates.
 
     Returns:
-        The expectation value estimated with Pauli twirling.
+        A list of `num_circuits` twirled versions of `circuit`
     """
     CNOT_twirled_circuits = twirl_CNOT_gates(circuit, num_circuits)
     twirled_circuits = [

mitiq/pt/tests/test_pt.py

@@ -20,7 +20,7 @@
     PENNYLANE_NOISE_OP,
     CNOT_twirling_gates,
     CZ_twirling_gates,
-    pauli_twirl_circuit,
+    generate_pauli_twirl_variants,
     twirl_CNOT_gates,
     twirl_CZ_gates,
 )
@@ -122,7 +122,7 @@ def test_twirl_CNOT_increases_layer_count():
         assert num_gates_after == num_gates_before
 
 
-def test_pauli_twirl_circuit():
+def test_generate_pauli_twirl_variants():
     num_qubits = 3
     num_layers = 20
     circuit, _ = generate_mirror_circuit(
@@ -131,12 +131,13 @@ def test_pauli_twirl_circuit():
         connectivity_graph=nx.complete_graph(num_qubits),
     )
     num_circuits = 10
-    twirled_output = pauli_twirl_circuit(circuit, num_circuits)
+    twirled_output = generate_pauli_twirl_variants(circuit, num_circuits)
     assert len(twirled_output) == num_circuits
 
 
 @pytest.mark.parametrize(
-    "twirl_func", [pauli_twirl_circuit, twirl_CNOT_gates, twirl_CZ_gates]
+    "twirl_func",
+    [generate_pauli_twirl_variants, twirl_CNOT_gates, twirl_CZ_gates],
 )
 def test_no_CNOT_CZ_circuit(twirl_func):
     num_qubits = 2
@@ -155,7 +156,7 @@ def test_noisy_cirq(noise_name):
     p = 0.01
     a, b = cirq.LineQubit.range(2)
     circuit = cirq.Circuit(cirq.H.on(a), cirq.CNOT.on(a, b), cirq.CZ.on(a, b))
-    twirled_circuit = pauli_twirl_circuit(
+    twirled_circuit = generate_pauli_twirl_variants(
         circuit, num_circuits=1, noise_name=noise_name, p=p
     )[0]
 
@@ -175,7 +176,7 @@ def test_noisy_pennylane(noise_name):
         qml.CZ((0, 1)),
     ]
     circuit = QuantumTape(ops)
-    twirled_circuit = pauli_twirl_circuit(
+    twirled_circuit = generate_pauli_twirl_variants(
         circuit, num_circuits=1, noise_name=noise_name, p=p
     )[0]