---
title: Utility III
description: This lesson is about to build a GHZ circuit for 20 qubits or more so that, upon measurement, the fidelity of GHZ state over 0.5.
source: https://eu-de.quantum.cloud.ibm.com/learning/en/courses/utility-scale-quantum-computing/utility-iii
---

# Utility-scale experiment III

> **Note**
>
> Toshinari Itoko, Tamiya Onodera, Kifumi Numata (19 July 2024)
>
> [Download the pdf](https://ibm.ent.box.com/public/static/fw1538dogvyv0qbfqg8tan1k2fs27mfv.zip) of the original lecture. Note that some code snippets might become deprecated since these are static images.
>
> *Approximate QPU time to run this first experiment is 12 m 30 s. There is an additional experiment below that requires approximately 4 m.*
>
> (Note: this notebook might not evaluate in the time allowed on the Open Plan. Be sure to use quantum computing resources wisely.)

```python
import qiskit

qiskit.__version__
```

Output:

```
'2.0.2'
```

```python
import qiskit_ibm_runtime

qiskit_ibm_runtime.__version__
```

Output:

```
'0.40.1'
```

```python
import numpy as np
import rustworkx as rx

from qiskit import QuantumCircuit
from qiskit.visualization import plot_histogram
from qiskit.visualization import plot_gate_map
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.providers import BackendV2
from qiskit.quantum_info import SparsePauliOp

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Sampler, Estimator, Batch, SamplerOptions
```

## 1. Introduction

Let us briefly review GHZ states, and what sort of distribution you might expect from `Sampler` applied to one. Then we will spell out the goal of this lesson explicitly.

### 1.1 GHZ state

The GHZ state (Greenberger-Horne-Zeilinger state) for $n$ qubits is defined as

$$
\frac{1}{\sqrt 2}(|0\rangle ^ {\otimes n}+ |1\rangle^ {\otimes n})
$$

Naturally, it can be created for 6 qubits with the following quantum circuit.

```python
N = 6
qc = QuantumCircuit(N, N)

qc.h(0)
for i in range(N - 1):
    qc.cx(0, i + 1)

# qc.measure_all()
qc.barrier()
qc.measure(list(range(N)), list(range(N)))

qc.draw(output="mpl", idle_wires=False, scale=0.5)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/1bf914c5-e1c6-4664-99c2-69dea3961114-0.avif)

```python
print("Depth:", qc.depth())
```

Output:

```
Depth: 7
```

The depth is not too large, though you know from previous lessons that you can do better. Let's pick a backend and transpile this circuit.

```python
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
backend.name
# or
# backend = service.least_busy(operational=True)
# backend.name
```

Output:

```
'ibm_kingston'
```

```python
pm = generate_preset_pass_manager(3, backend=backend)
qc_transpiled = pm.run(qc)
qc_transpiled.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/2e003078-3707-44d9-86e8-bbc66b095e84-0.avif)

```python
print("Depth:", qc_transpiled.depth())
print(
    "Two-qubit Depth:",
    qc_transpiled.depth(filter_function=lambda x: x.operation.num_qubits == 2),
)
```

Output:

```
Depth: 27
Two-qubit Depth: 11
```

Again the transpiled two-qubit depth is not too large. But to work with a GHZ state on more qubits, you will clearly need to think about optimizing the circuit. Let's run this using `Sampler` and see what a real quantum computer returns.

```python
sampler = Sampler(mode=backend)
shots = 40000
job = sampler.run([qc_transpiled], shots=shots)
job_id = job.job_id()
print(job_id)
```

Output:

```
d147y20n2txg008jvv70
```

```python
job.status()
```

Output:

```
'DONE'
```

```python
job = service.job(job_id)
result = job.result()
plot_histogram(result[0].data.c.get_counts(), figsize=(30, 5))
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/8ab5b37f-fb87-40a8-8e13-0c6cb3154f64-0.avif)

This is the result of the 6-qubit GHZ circuit. As you can see, the states of all $|0\rangle$'s and all $|1\rangle$'s do dominate, but the errors are substantial. Let's try to see how large a GHZ circuit you can make with an Eagle device, while still getting results where the correct states are at least more than 50% likely.

### 1.2 Your goal

Build a GHZ circuit for 20 qubits or more so that, upon measurement, **the fidelity of your GHZ state > 0.5.**
Note:

- You need to use an Eagle device (`min_num_qubits=127`) and set the shots number as 40,000.
- You should execute the GHZ circuit using the `execute_ghz_fidelity` function, and calculate the fidelity using the `check_ghz_fidelity_from_jobs` function.

This is intended as an independent exercise, in which you leverage what you have learned so far in this course.

```python
def execute_ghz_fidelity(
    ghz_circuit: QuantumCircuit,  # Quantum circuit to create GHZ state (Circuit after Routing or without Routing), Classical register name is "c"
    physical_qubits: list[int],  # Physical qubits to represent GHZ state
    backend: BackendV2,
    sampler_options: dict | SamplerOptions | None = None,
):
    N_SHOTS = 40_000
    N = len(physical_qubits)
    base_circuit = ghz_circuit.remove_final_measurements(inplace=False)
    # M_k measurement circuits
    mk_circuits = []
    for k in range(1, N + 1):
        circuit = base_circuit.copy()
        # change measurement basis
        for q in physical_qubits:
            circuit.rz(-k * np.pi / N, q)
            circuit.h(q)
        mk_circuits.append(circuit)

    obs = SparsePauliOp.from_sparse_list(
        [("Z" * N, physical_qubits, 1)], num_qubits=backend.num_qubits
    )
    job_ids = []
    pm1 = generate_preset_pass_manager(1, backend=backend)
    org_transpiled = pm1.run(ghz_circuit)
    mk_transpiled = pm1.run(mk_circuits)
    with Batch(backend=backend):
        sampler = Sampler(options=sampler_options)
        sampler.options.twirling.enable_measure = True
        job = sampler.run([org_transpiled], shots=N_SHOTS)
        job_ids.append(job.job_id())
        # print(f"Sampler job id: {job.job_id()}, shots={N_SHOTS}")
        estimator = Estimator()  # TREX is applied as default
        estimator.options.dynamical_decoupling.enable = True
        estimator.options.execution.rep_delay = 0.0005
        estimator.options.twirling.enable_measure = True
        job2 = estimator.run([(circ, obs) for circ in mk_transpiled], precision=1 / 100)
        job_ids.append(job2.job_id())
        # print("Estimator job id:", job2.job_id())
        return [job.job_id(), job2.job_id()]
```

```python
def check_ghz_fidelity_from_jobs(
    sampler_job,
    estimator_job,
    num_qubits,
    shots=40_000,
):
    N = num_qubits
    sampler_result = sampler_job.result()
    counts = sampler_result[0].data.c.get_counts()
    all_zero = counts.get("0" * N, 0) / shots
    all_one = counts.get("1" * N, 0) / shots
    top3 = sorted(counts, key=counts.get, reverse=True)[:3]
    print(
        f"N={N}: |00..0>: {counts.get('0'*N, 0)}, |11..1>: {counts.get('1'*N, 0)}, |3rd>: {counts.get(top3[2], 0)} ({top3[2]})"
    )
    print(f"P(|00..0>)={all_zero}, P(|11..1>)={all_one}")

    estimator_result = estimator_job.result()
    non_diagonal = (1 / N) * sum(
        (-1) ** k * estimator_result[k - 1].data.evs for k in range(1, N + 1)
    )
    print(f"REM: Coherence (non-diagonal): {non_diagonal:.6f}")
    fidelity = 0.5 * (all_zero + all_one + non_diagonal)
    sigma = 0.5 * np.sqrt(
        (1 - all_zero - all_one) * (all_zero + all_one) / shots
        + sum(estimator_result[k].data.stds ** 2 for k in range(N)) / (N * N)
    )
    print(f"GHZ fidelity = {fidelity:.6f} ± {sigma:.6f}")
    if fidelity - 2 * sigma > 0.5:
        print("GME (genuinely multipartite entangled) test: Passed")
    else:
        print("GME (genuinely multipartite entangled) test: Failed")
    return {
        "fidelity": fidelity,
        "sigma": sigma,
        "shots": shots,
        "job_ids": [sampler_job.job_id(), estimator_job.job_id()],
    }
```

In this notebook, we will apply three strategies to creating good GHZ states using 16 qubits and 30 qubits. These approaches build on strategies you already know from previous lessons.

## 2. Strategy 1. Noise-aware qubit selection

We first specify a backend. Because we will be working extensively with the properties of a specific backend, it is a good idea to specify one backend, as opposed to using the `least_busy` option.

```python
backend = service.backend("ibm_strasbourg")  # eagle
twoq_gate = "ecr"
print(f"Device {backend.name} Loaded with {backend.num_qubits} qubits")
print(f"Two Qubit Gate: {twoq_gate}")
```

Output:

```
Device ibm_strasbourg Loaded with 127 qubits
Two Qubit Gate: ecr
```

We are going to build a circuit involving many two-qubit gates. It makes sense for us to use the qubits that have the lowest errors when implementing those two-qubit gates. Finding the best "qubit chain" based on the reported 2q-gate errors is a nontrivial problem. But we can define a few functions to help us determine the best qubits to use.

```python
coupling_map = backend.target.build_coupling_map(twoq_gate)
G = coupling_map.graph
```

```python
def to_edges(path):  # create edges list from node paths
    edges = []
    prev_node = None
    for node in path:
        if prev_node is not None:
            if G.has_edge(prev_node, node):
                edges.append((prev_node, node))
            else:
                edges.append((node, prev_node))
        prev_node = node
    return edges


def path_fidelity(path, correct_by_duration: bool = True, readout_scale: float = None):
    """Compute an estimate of the total fidelity of 2-qubit gates on a path.
    If `correct_by_duration` is true, each gate fidelity is worsen by
    scale = max_duration / duration, that is, gate_fidelity^scale.
    If `readout_scale` > 0 is supplied, readout_fidelity^readout_scale
    for each qubit on the path is multiplied to the total fielity.
    The path is given in node indices form, for example, [0, 1, 2].
    An external function `to_edges` is used to obtain edge list, for example, [(0, 1), (1, 2)]."""
    path_edges = to_edges(path)
    max_duration = max(backend.target[twoq_gate][qs].duration for qs in path_edges)

    def gate_fidelity(qpair):
        duration = backend.target[twoq_gate][qpair].duration
        scale = max_duration / duration if correct_by_duration else 1.0
        # 1.25 = (d+1)/d with d = 4
        return max(0.25, 1 - (1.25 * backend.target[twoq_gate][qpair].error)) ** scale

    def readout_fidelity(qubit):
        return max(0.25, 1 - backend.target["measure"][(qubit,)].error)

    total_fidelity = np.prod(
        [gate_fidelity(qs) for qs in path_edges]
    )  # two qubits gate fidelity for each path
    if readout_scale:
        total_fidelity *= (
            np.prod([readout_fidelity(q) for q in path]) ** readout_scale
        )  # multiply readout fidelity
    return total_fidelity


def flatten(paths, cutoff=None):  # cutoff is for not making run time too large
    return [
        path
        for s, s_paths in paths.items()
        for t, st_paths in s_paths.items()
        for path in st_paths[:cutoff]
        if s < t
    ]
```

```python
N = 16  # Number of qubits to use in the GHZ circuit
num_qubits_in_chain = N
```

We will use the functions above to find all the simple paths of N qubits between all pairs of nodes in the graph (Reference: [all\_pairs\_all\_simple\_paths](https://www.rustworkx.org/apiref/rustworkx.all_pairs_all_simple_paths.html#rustworkx-all-pairs-all-simple-paths)).

Then, using the `path_fidelity` function created above, we will find the best qubit chain which has the largest path fidelity.

```python
from functools import partial

%%time
paths = rx.all_pairs_all_simple_paths(
    G.to_undirected(multigraph=False),
    min_depth=num_qubits_in_chain,
    cutoff=num_qubits_in_chain,
)
paths = flatten(paths, cutoff=25)  # If you have time, you could set a larger cutoff.
if not paths:
    raise Exception(
        f"No qubit chain with length={num_qubits_in_chain} exists in {backend.name}. Try smaller num_qubits_in_chain."
    )

print(f"Selecting the best from {len(paths)} candidate paths")

best_qubit_chain = max(
    paths, key=partial(path_fidelity, correct_by_duration=True, readout_scale=1.0)
)
assert len(best_qubit_chain) == num_qubits_in_chain
print(f"Predicted (best possible) process fidelity: {path_fidelity(best_qubit_chain)}")
```

Output:

```
Selecting the best from 6046 candidate paths
Predicted (best possible) process fidelity: 0.8929026784775056
CPU times: user 284 ms, sys: 10.9 ms, total: 295 ms
Wall time: 295 ms
```

```python
np.array(best_qubit_chain)
```

Output:

```
array([55, 49, 48, 47, 46, 45, 54, 64, 65, 66, 73, 85, 86, 87, 88, 89],
      dtype=uint64)
```

Let's plot the best qubit chain, shown in pink, in the coupling map diagram.

```python
qubit_color = []
for i in range(133):
    if i in best_qubit_chain:
        qubit_color.append("#ff00dd")  # pink
    else:
        qubit_color.append("#8c00ff")  # purple
plot_gate_map(
    backend, qubit_color=qubit_color, qubit_size=50, font_size=25, figsize=(6, 6)
)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/a1c164a4-0.avif)

### 2.1 Build a GHZ circuit on the best qubit chain

We choose a qubit in the middle of the chain to first apply the H gate to. This should reduce the depth of the circuit by about half.

```python
ghz1 = QuantumCircuit(max(best_qubit_chain) + 1, N)
ghz1.h(best_qubit_chain[N // 2])
for i in range(N // 2, 0, -1):
    ghz1.cx(best_qubit_chain[i], best_qubit_chain[i - 1])
for i in range(N // 2, N - 1, +1):
    ghz1.cx(best_qubit_chain[i], best_qubit_chain[i + 1])
ghz1.barrier()  # for visualization
ghz1.measure(best_qubit_chain, list(range(N)))
ghz1.draw(output="mpl", idle_wires=False, scale=0.5, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/1336fba0-0.avif)

```python
ghz1.depth()
```

Output:

```
10
```

```python
pm = generate_preset_pass_manager(1, backend=backend)
ghz1_transpiled = pm.run(ghz1)
ghz1_transpiled.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/71874b6d-eb9b-4ff5-b673-944191f67010-0.avif)

```python
print("Depth:", ghz1_transpiled.depth())
print(
    "Two-qubit Depth:",
    ghz1_transpiled.depth(filter_function=lambda x: x.operation.num_qubits == 2),
)
```

Output:

```
Depth: 27
Two-qubit Depth: 8
```

```python
opts = SamplerOptions()
```

```python
res = execute_ghz_fidelity(
    ghz_circuit=ghz1,
    physical_qubits=best_qubit_chain,
    backend=backend,
    sampler_options=opts,
)
```

```python
job_s = service.job(res[0])  # Use your job id showed above.
job_e = service.job(res[1])
print(job_s.status(), job_e.status())
```

Output:

```
DONE DONE
```

Be careful to execute the next cell after the above jobs statuses have become 'DONE', to show the result using the `check_ghz_fidelity_from_jobs` function.

```python
N = 16
# Check fidelity from job IDs
res = check_ghz_fidelity_from_jobs(
    sampler_job=job_s,
    estimator_job=job_e,
    num_qubits=N,
)
```

Output:

```
N=16: |00..0>: 153, |11..1>: 8681, |3rd>: 2262 (1111111111101111)
P(|00..0>)=0.003825, P(|11..1>)=0.217025
REM: Coherence (non-diagonal): 0.073809
GHZ fidelity = 0.147329 ± 0.002438
GME (genuinely multipartite entangled) test: Failed
```

```python
result = job_s.result()
plot_histogram(result[0].data.c.get_counts(), figsize=(30, 5))
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/236126cb-0.avif)

This result does not meet the criteria. Let's move to the next idea.

## 3. Strategy 2. Balanced tree of qubits

The next idea is to find a balanced tree of qubits. Using the tree rather than the chain, the circuit depth should become lower. Before that, we remove nodes with "bad" readout errors and edges with "bad" gate errors from the coupling graph.

```python
BAD_READOUT_ERROR_THRESHOLD = 0.1
BAD_ECRGATE_ERROR_THRESHOLD = 0.1
bad_readout_qubits = [
    q
    for q in range(backend.num_qubits)
    if backend.target["measure"][(q,)].error > BAD_READOUT_ERROR_THRESHOLD
]
bad_ecrgate_edges = [
    qpair
    for qpair in backend.target["ecr"]
    if backend.target["ecr"][qpair].error > BAD_ECRGATE_ERROR_THRESHOLD
]
print("Bad readout qubits:", bad_readout_qubits)
print("Bad ECR gates:", bad_ecrgate_edges)
```

Output:

```
Bad readout qubits: [19, 28, 41, 72, 91, 114, 120]
Bad ECR gates: []
```

```python
g = backend.coupling_map.graph.copy().to_undirected()
g.remove_edges_from(
    bad_ecrgate_edges
)  # remove edge first (otherwise might fail with a NoEdgeBetweenNodes error)
g.remove_nodes_from(bad_readout_qubits)
```

Let's draw the coupling map graph without the bad edges and bad qubits.

```python
qubit_color = []
for i in range(133):
    if i in bad_readout_qubits:
        qubit_color.append("#000000")  # black
    else:
        qubit_color.append("#8c00ff")  # purple
line_color = []
for e in backend.target.build_coupling_map().get_edges():
    if e in bad_ecrgate_edges:
        line_color.append("#ffffff")  # white
    else:
        line_color.append("#888888")  # gray
plot_gate_map(
    backend,
    qubit_color=qubit_color,
    line_color=line_color,
    qubit_size=50,
    font_size=25,
    figsize=(6, 6),
)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/18825a87-4679-4027-af58-31efd24a9a60-0.avif)

We try to create a 16-qubit GHZ state as before.

```python
N = 16
```

We call the `betweenness_centrality` function to find a qubit for the root node. The node with the highest value of betweenness centrality is at the center of the graph. Reference: [https://www.rustworkx.org/tutorial/betweenness\_centrality.html](https://www.rustworkx.org/tutorial/betweenness_centrality.html)

Or you can select it manually.

```python
# central = 65 #Select the center node manually
c_degree = dict(rx.betweenness_centrality(g))
central = max(c_degree, key=c_degree.get)
central
```

Output:

```
66
```

Starting from the root node, generate a tree by breadth first search (BFS). Reference: [https://qiskit.org/ecosystem/rustworkx/apiref/rustworkx.bfs\_search.html#rustworkx-bfs-search](https://qiskit.org/ecosystem/rustworkx/apiref/rustworkx.bfs_search.html#rustworkx-bfs-search)

```python
class TreeEdgesRecorder(rx.visit.BFSVisitor):
    def __init__(self, N):
        self.edges = []
        self.N = N

    def tree_edge(self, edge):
        self.edges.append(edge)
        if len(self.edges) >= self.N - 1:
            raise rx.visit.StopSearch()


vis = TreeEdgesRecorder(N)
rx.bfs_search(g, [central], vis)
best_qubits = sorted(list(set(q for e in vis.edges for q in (e[0], e[1]))))
# print('Tree edges:', vis.edges)
```

```python
print("Qubits selected:", best_qubits)
```

Output:

```
Qubits selected: [54, 55, 63, 64, 65, 66, 67, 68, 69, 70, 73, 83, 84, 85, 86, 87]
```

Let us plot the selected qubits, shown in pink, in the coupling map diagram.

```python
qubit_color = []
for i in range(133):
    if i in bad_readout_qubits:
        qubit_color.append("#000000")  # black
    elif i in best_qubits:
        qubit_color.append("#ff00dd")  # pink
    else:
        qubit_color.append("#8c00ff")  # purple
plot_gate_map(
    backend,
    qubit_color=qubit_color,
    line_color=line_color,
    qubit_size=50,
    font_size=25,
    figsize=(6, 6),
)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/ef18e037-a3ec-4370-bee7-bd71dbc12795-0.avif)

Let us show the tree structure of qubits.

```python
from rustworkx.visualization import graphviz_draw

tree = rx.PyDiGraph()
tree.extend_from_weighted_edge_list(vis.edges)
tree.remove_nodes_from([n for n in range(max(best_qubits) + 1) if n not in best_qubits])

graphviz_draw(tree, method="dot")
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/8758d3a6-0.avif)

```python
ghz2 = QuantumCircuit(max(best_qubits) + 1, N)

ghz2.h(tree.edge_list()[0][0])  # apply H-gate to the root node
# Apply CNOT from the root node to the each edge.
for u, v in tree.edge_list():
    ghz2.cx(u, v)
ghz2.barrier()  # for visualization
ghz2.measure(best_qubits, list(range(N)))
ghz2.draw(output="mpl", idle_wires=False, scale=0.5)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/3bc4b071-0.avif)

```python
ghz2.depth()
```

Output:

```
8
```

```python
pm = generate_preset_pass_manager(1, backend=backend)
ghz2_transpiled = pm.run(ghz2)
ghz2_transpiled.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/c490e37a-309a-43e3-9b18-3be7c7d8a957-0.avif)

```python
print("Depth:", ghz2_transpiled.depth())
print(
    "Two-qubit Depth:",
    ghz2_transpiled.depth(filter_function=lambda x: x.operation.num_qubits == 2),
)
```

Output:

```
Depth: 22
Two-qubit Depth: 6
```

The depth of the circuit has now become much lower than that in the chain structure.

```python
res = execute_ghz_fidelity(
    ghz_circuit=ghz2,
    physical_qubits=best_qubits,
    backend=backend,
    sampler_options=opts,
)
```

```python
job_s = service.job(res[0])  # Use your job id showed above.
job_e = service.job(res[1])
print(job_s.status(), job_e.status())
```

Output:

```
DONE DONE
```

```python
N = 16
# Check fidelity from job IDs
res = check_ghz_fidelity_from_jobs(
    sampler_job=job_s,
    estimator_job=job_e,
    num_qubits=N,
)
```

Output:

```
N=16: |00..0>: 9509, |11..1>: 10978, |3rd>: 1795 (1111110111111111)
P(|00..0>)=0.237725, P(|11..1>)=0.27445
REM: Coherence (non-diagonal): 0.606515
GHZ fidelity = 0.559345 ± 0.003188
GME (genuinely multipartite entangled) test: Passed
```

We have successfully passed the criteria with the balanced tree structure!

```python
result = job_s.result()
plot_histogram(result[0].data.c.get_counts(), figsize=(30, 5))
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/cf8b4dac-0.avif)

Now, let us try to create a larger GHZ state: a 30-qubit GHZ state.

### 3.1 N = 30

We will follow the [Qiskit patterns](/docs/guides/intro-to-patterns) framework.

- Step 1: Map problem to quantum circuits and operators
- Step 2: Optimize for target hardware
- Step 3: Execute on target hardware
- Step 4: Post-process results

#### Step 1: Map problem to quantum circuits and operators and Step 2: Optimize for target hardware

Here we select the root node manually.

```python
central = 62  # Select the center node manually
# c_degree = dict(rx.betweenness_centrality(g))
# central = max(c_degree, key=c_degree.get)
# central
```

```python
N = 30

vis = TreeEdgesRecorder(N)
rx.bfs_search(g, [central], vis)
best_qubits = sorted(list(set(q for e in vis.edges for q in (e[0], e[1]))))
print("Qubits selected:", best_qubits)
```

Output:

```
Qubits selected: [34, 35, 42, 43, 44, 45, 46, 47, 48, 53, 54, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 71, 73, 77, 84, 85, 86]
```

```python
qubit_color = []
for i in range(133):
    if i in bad_readout_qubits:
        qubit_color.append("#000000")
    elif i in best_qubits:
        qubit_color.append("#ff00dd")
    else:
        qubit_color.append("#8c00ff")
line_color = []
for e in backend.target.build_coupling_map().get_edges():
    if e in bad_ecrgate_edges:
        line_color.append("#ffffff")
    else:
        line_color.append("#888888")
plot_gate_map(
    backend,
    qubit_color=qubit_color,
    line_color=line_color,
    qubit_size=50,
    font_size=25,
    figsize=(6, 6),
)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/12d3ac00-7dfc-4734-9881-4e6385f2d880-0.avif)

```python
from rustworkx.visualization import graphviz_draw

tree = rx.PyDiGraph()
tree.extend_from_weighted_edge_list(vis.edges)
tree.remove_nodes_from([n for n in range(max(best_qubits) + 1) if n not in best_qubits])

graphviz_draw(tree, method="dot")
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/6afd6c54-b427-4f34-8c4c-285c0122aebd-0.avif)

The depth of this tree is 5.

```python
ghz3 = QuantumCircuit(max(best_qubits) + 1, N)

ghz3.h(tree.edge_list()[0][0])  # apply H-gate to the root node
# Apply CNOT from the root node to the each edge.
for u, v in tree.edge_list():
    ghz3.cx(u, v)
ghz3.barrier()  # for visualization
ghz3.measure(best_qubits, list(range(N)))
ghz3.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/d36cdb58-59fe-474e-9a31-a2d9e806ab59-0.avif)

```python
ghz3.depth()
```

Output:

```
11
```

```python
pm = generate_preset_pass_manager(1, backend=backend)
ghz3_transpiled = pm.run(ghz3)
ghz3_transpiled.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/46df385d-dd28-4fab-8344-6ef2fd53906d-0.avif)

```python
print("Depth:", ghz3_transpiled.depth())
print(
    "Two-qubit Depth:",
    ghz3_transpiled.depth(filter_function=lambda x: x.operation.num_qubits == 2),
)
```

Output:

```
Depth: 31
Two-qubit Depth: 9
```

### 3.2 Select a different root node manually

```python
central = 54

vis = TreeEdgesRecorder(N)
rx.bfs_search(g, [central], vis)
best_qubits = sorted(list(set(q for e in vis.edges for q in (e[0], e[1]))))
print("Qubits selected:", best_qubits)
```

Output:

```
Qubits selected: [23, 24, 25, 34, 35, 42, 43, 44, 45, 46, 47, 48, 49, 50, 54, 55, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 73, 84, 85, 86]
```

```python
from rustworkx.visualization import graphviz_draw

tree = rx.PyDiGraph()
tree.extend_from_weighted_edge_list(vis.edges)
tree.remove_nodes_from([n for n in range(max(best_qubits) + 1) if n not in best_qubits])

graphviz_draw(tree, method="dot")
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/5e2e2af5-34b5-4bb5-b393-18e87d56db08-0.avif)

The depth of this tree is 6.

```python
ghz3 = QuantumCircuit(max(best_qubits) + 1, N)

ghz3.h(tree.edge_list()[0][0])  # apply H-gate to the root node
# Apply CNOT from the root node to the each edge.
for u, v in tree.edge_list():
    ghz3.cx(u, v)
ghz3.barrier()  # for visualization
ghz3.measure(best_qubits, list(range(N)))
ghz3.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/5d81d890-a939-455f-81a6-c7e671671ba5-0.avif)

```python
ghz3.depth()
```

Output:

```
11
```

```python
pm = generate_preset_pass_manager(1, backend=backend)
ghz3_transpiled = pm.run(ghz3)
ghz3_transpiled.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/b9823b6d-a658-40ef-a509-f03e33534af6-0.avif)

```python
print("Depth:", ghz3_transpiled.depth())
print(
    "Two-qubit Depth:",
    ghz3_transpiled.depth(filter_function=lambda x: x.operation.num_qubits == 2),
)
```

Output:

```
Depth: 30
Two-qubit Depth: 9
```

Surprisingly, while the tree depth increased from 5 to 6, the two-qubit depth decreased from 9 to 8! So let us use the latter circuit.

#### Step 3: Execute on target hardware

```python
res = execute_ghz_fidelity(
    ghz_circuit=ghz3,
    physical_qubits=best_qubits,
    backend=backend,
    sampler_options=opts,
)
```

```python
job_s = service.job(res[0])  # Use your job id showed above.
job_e = service.job(res[1])
print(job_s.status(), job_e.status())
```

Output:

```
DONE DONE
```

#### Step 4: Post-process results

```python
N = 30
# Check fidelity from job IDs
res = check_ghz_fidelity_from_jobs(
    sampler_job=job_s,
    estimator_job=job_e,
    num_qubits=N,
)
```

Output:

```
N=30: |00..0>: 4, |11..1>: 218, |3rd>: 265 (111111111111111011111111111111)
P(|00..0>)=0.0001, P(|11..1>)=0.00545
REM: Coherence (non-diagonal): 0.187073
GHZ fidelity = 0.096312 ± 0.003254
GME (genuinely multipartite entangled) test: Failed
```

As you can see, this result has not met the criteria.

```python
# It will take some time
result = job_s.result()
plot_histogram(result[0].data.c.get_counts(), figsize=(30, 5))
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/66c87870-0.avif)

## 4. Strategy 3. Run with the error suppression options

You can set the error suppression options in Sampler V2. Refer to the [Sampler noise management](/docs/guides/sampler-noise-management) guide and the [`ExecutionOptionsV2`](/docs/api/qiskit-ibm-runtime/options-execution-options-v2) API reference for more information.

```python
opts = SamplerOptions()
opts.dynamical_decoupling.enable = True
opts.execution.rep_delay = 0.0005
opts.twirling.enable_gates = True
```

```python
res = execute_ghz_fidelity(
    ghz_circuit=ghz3,
    physical_qubits=best_qubits,
    backend=backend,
    sampler_options=opts,
)
```

```python
job_s = service.job(res[0])  # Use your job id showed above.
job_e = service.job(res[1])
print(job_s.status(), job_e.status())
```

Output:

```
DONE DONE
```

```python
N = 30
```

```python
# Check fidelity from job IDs
res = check_ghz_fidelity_from_jobs(
    sampler_job=job_s,
    estimator_job=job_e,
    num_qubits=N,
)
```

Output:

```
N=30: |00..0>: 1459, |11..1>: 1543, |3rd>: 359 (111111111111111111111111111110)
P(|00..0>)=0.036475, P(|11..1>)=0.038575
REM: Coherence (non-diagonal): 0.165532
GHZ fidelity = 0.120291 ± 0.003369
GME (genuinely multipartite entangled) test: Failed
```

```python
# It will take some time
result = job_s.result()
plot_histogram(result[0].data.c.get_counts(), figsize=(30, 5))
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/166856be-ad3a-40d6-837c-3eedda374891-0.avif)

The result has improved but still has not met the criteria.

We have seen three ideas so far. You can combine and expand these ideas or you come up with your own ideas to create a better GHZ circuit. Now let's review the goal again.

## 5. Your goal (recap)

Build a GHZ circuit for 20 qubits or more so that the measurement result meets the criteria: The fidelity of your GHZ state > 0.5.

- You need to use an Eagle device (such as `ibm_brisbane`) and set the shots number as 40,000.
- You should execute the GHZ circuit using the `execute_ghz_fidelity` function, and calculate the fidelity using the `check_ghz_fidelity_from_jobs` function.
  You need to find the biggest qubits - GHZ circuit which meet the criteria. Write your code below, show the result with the function `check_ghz_fidelity_from_jobs` .

Now we implement the same GHZ workflow as in the preceding material, but on a Heron device. This gives you some experience with the layout and features of the Heron processors. No new strategies are introduced.

*Approximate QPU time to run this next experiment is 4 m 40 s.*

```python
service = QiskitRuntimeService()
backend = service.backend("ibm_kingston")
# backend = service.backend("ibm_fez")

twoq_gate = "cz"
print(f"Device {backend.name} Loaded with {backend.num_qubits} qubits")
print(f"Two Qubit Gate: {twoq_gate}")
```

Output:

```
Device ibm_kingston Loaded with 156 qubits
Two Qubit Gate: cz
```

```python
BAD_READOUT_ERROR_THRESHOLD = 0.1
BAD_CZGATE_ERROR_THRESHOLD = 0.1
bad_readout_qubits = [
    q
    for q in range(backend.num_qubits)
    if backend.target["measure"][(q,)].error > BAD_READOUT_ERROR_THRESHOLD
]
bad_czgate_edges = [
    qpair
    for qpair in backend.target["cz"]
    if backend.target["cz"][qpair].error > BAD_CZGATE_ERROR_THRESHOLD
]
print("Bad readout qubits:", bad_readout_qubits)
print("Bad CZ gates:", bad_czgate_edges)
```

Output:

```
Bad readout qubits: [112, 113, 120, 131, 146]
Bad CZ gates: [(111, 112), (112, 111), (112, 113), (113, 112), (120, 121), (121, 120), (130, 131), (131, 130), (145, 146), (146, 145), (146, 147), (147, 146)]
```

```python
g = backend.coupling_map.graph.copy().to_undirected()
g.remove_edges_from(
    bad_czgate_edges
)  # remove edge first (otherwise might fail with a NoEdgeBetweenNodes error)
g.remove_nodes_from(bad_readout_qubits)
```

```python
qubit_color = []
for i in range(backend.num_qubits):
    if i in bad_readout_qubits:
        qubit_color.append("#000000")  # black
    else:
        qubit_color.append("#8c00ff")  # purple
line_color = []
for e in backend.target.build_coupling_map().get_edges():
    if e in bad_czgate_edges:
        line_color.append("#ffffff")  # white
    else:
        line_color.append("#888888")  # gray
plot_gate_map(
    backend,
    qubit_color=qubit_color,
    line_color=line_color,
    qubit_size=60,
    font_size=30,
    figsize=(10, 10),
)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/e4336cf4-95e6-41a2-b24c-74d1cec53c36-0.avif)

```python
N = 40
central = 100  # Select the center node manually
# c_degree = dict(rx.betweenness_centrality(g))
# central = max(c_degree, key=c_degree.get)
# central
```

```python
class TreeEdgesRecorder(rx.visit.BFSVisitor):
    def __init__(self, N):
        self.edges = []
        self.N = N

    def tree_edge(self, edge):
        self.edges.append(edge)
        if len(self.edges) >= self.N - 1:
            raise rx.visit.StopSearch()


vis = TreeEdgesRecorder(N)
rx.bfs_search(g, [central], vis)
best_qubits = sorted(list(set(q for e in vis.edges for q in (e[0], e[1]))))
print("Qubits selected:", best_qubits)
```

Output:

```
Qubits selected: [61, 65, 76, 77, 80, 81, 82, 83, 84, 85, 86, 87, 96, 97, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 116, 117, 121, 122, 123, 124, 125, 126, 127, 136, 140, 141, 142, 143, 144, 145]
```

```python
qubit_color = []
for i in range(backend.num_qubits):
    if i in bad_readout_qubits:
        qubit_color.append("#000000")
    elif i in best_qubits:
        qubit_color.append("#ff00dd")
    else:
        qubit_color.append("#8c00ff")
line_color = []
for e in backend.target.build_coupling_map().get_edges():
    if e in bad_czgate_edges:
        line_color.append("#ffffff")
    else:
        line_color.append("#888888")
plot_gate_map(
    backend,
    qubit_color=qubit_color,
    line_color=line_color,
    qubit_size=60,
    font_size=30,
    figsize=(10, 10),
)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/e684af6c-aa99-460d-b0f5-9d5042fb80b0-0.avif)

```python
from rustworkx.visualization import graphviz_draw

tree = rx.PyDiGraph()
tree.extend_from_weighted_edge_list(vis.edges)
tree.remove_nodes_from([n for n in range(max(best_qubits) + 1) if n not in best_qubits])

graphviz_draw(tree, method="dot")
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/f3ff2dd9-a304-409c-a0e4-3b23e85f282e-0.avif)

```python
ghz_h = QuantumCircuit(max(best_qubits) + 1, N)

ghz_h.h(tree.edge_list()[0][0])  # apply H-gate to the root node
# Apply CNOT from the root node to the each edge.
for u, v in tree.edge_list():
    ghz_h.cx(u, v)
ghz_h.barrier()  # for visualization
ghz_h.measure(best_qubits, list(range(N)))
ghz_h.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/60b6526c-55ed-4621-b91b-4ceeef34de62-0.avif)

```python
ghz_h.depth()
```

Output:

```
15
```

```python
pm = generate_preset_pass_manager(1, backend=backend)
ghz_h_transpiled = pm.run(ghz_h)
ghz_h_transpiled.draw(output="mpl", idle_wires=False, fold=-1)
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/b96e6342-e08b-4728-8985-7a015ee16797-0.avif)

```python
print("Depth:", ghz_h_transpiled.depth())
print(
    "Two-qubit Depth:",
    ghz_h_transpiled.depth(filter_function=lambda x: x.operation.num_qubits == 2),
)
```

Output:

```
Depth: 45
Two-qubit Depth: 13
```

```python
opts = SamplerOptions()
opts.dynamical_decoupling.enable = True
opts.execution.rep_delay = 0.0005
opts.twirling.enable_gates = True
```

```python
res = execute_ghz_fidelity(
    ghz_circuit=ghz_h,
    physical_qubits=best_qubits,
    backend=backend,
    sampler_options=opts,
)
```

```python
job_s = service.job(res[0])  # Use your job id showed above.
job_e = service.job(res[1])
print(job_s.status(), job_e.status())
```

Output:

```
RUNNING RUNNING
```

```python
# Check fidelity from job IDs
N = 40
res = check_ghz_fidelity_from_jobs(
    sampler_job=job_s,
    estimator_job=job_e,
    num_qubits=N,
)
```

Output:

```
N=40: |00..0>: 3186, |11..1>: 3277, |3rd>: 620 (1111111011111111111111111111111111111111)
P(|00..0>)=0.07965, P(|11..1>)=0.081925
REM: Coherence (non-diagonal): 0.029987
GHZ fidelity = 0.095781 ± 0.002619
GME (genuinely multipartite entangled) test: Failed
```

```python
# It will take some time
result = job_s.result()
plot_histogram(result[0].data.c.get_counts(), figsize=(30, 5))
```

Output:

![Output of the previous code cell](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/utility-iii/extracted-outputs/671f57ea-7117-4d17-998f-73baa6f4463e-0.avif)

Congratulations! You have completed your introduction to utility-scale quantum computing! You are now poised to make meaningful contributions in the quest for quantum advantage! Thank you for making IBM Quantum® part of your personal quantum journey.

## Post-course survey

Congratulations on completing this course! Please take a moment to help us improve our course by filling out the following [quick survey](https://your.feedback.ibm.com/jfe/form/SV_5vvmdT5yFFetkgK). Your feedback will be used to enhance our content offering and user experience. Thank you!